# HG changeset patch
# User Markus Mützel <markus.muetzel@gmx.de>
# Date 1638189809 -3600
#      Mon Nov 29 13:43:29 2021 +0100
# Node ID e40a599d68cf3061d04a9dac30e36751ed20acb2
# Parent  71f2c8fde0c59f4942dd22da89bba162e8ae6c97
Remove usage of `error_state` (bug #).

* src/pcregexp.cc (many files): Remove usage of `error_state`. It was
unconditionally set to 0 since about 6 years ago and will finally be removed in
Octave 8.

diff -r 71f2c8fde0c5 -r e40a599d68cf src/__boxcount__.cc
--- a/src/__boxcount__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__boxcount__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -194,89 +194,85 @@
 
       octave_idx_type length=LENGTH-(maxembed-1)*DELAY;
 
-      if ( ! error_state)
+      // Calculate output
+      double heps              = EPSMAX*EPSFAKTOR;
+      octave_idx_type epsi_old = 0;
+      for (octave_idx_type k=0;k<EPSCOUNT;k++)
         {
-          // Calculate output
-          double heps              = EPSMAX*EPSFAKTOR;
-          octave_idx_type epsi_old = 0;
-          for (octave_idx_type k=0;k<EPSCOUNT;k++)
-            {
+
+          octave_idx_type epsi_test;
+          do {
+            heps /= EPSFAKTOR;
+            epsi_test=(octave_idx_type)(1./heps);
+          } while (epsi_test <= epsi_old);
+
+          octave_idx_type epsi = epsi_test;
+          epsi_old             = epsi;
+          deps[k]              = heps;
+
+          std::vector <double> histo (maxembed * dimension, 0.0);
+          start_box(series, which_dims, histo, maxembed, dimension, DELAY,
+                    length, epsi, Q);
 
-              octave_idx_type epsi_test;
-              do {
-                heps /= EPSFAKTOR;
-                epsi_test=(octave_idx_type)(1./heps);
-              } while (epsi_test <= epsi_old);
+          if (Q != 1.0)
+            for (std::vector<double>::iterator it = histo.begin ();
+                 it != histo.end (); it++)
+              *it = log(*it)/(1.0-Q);
+
+          histo_list.push_back (histo);
+        }
+
+      // Create and assign output
+      dim_vector dv (maxembed, dimension);
+      string_vector keys;
+      keys.append (std::string("dim"));
+      keys.append (std::string("entropy"));
+      octave_map output (dv, keys);
 
-              octave_idx_type epsi = epsi_test;
-              epsi_old             = epsi;
-              deps[k]              = heps;
+      for (octave_idx_type i=0;i<maxembed*dimension;i++)
+        {
+          octave_scalar_map tmp (keys);
+          // old fprintf(fHq,"#component = %d embedding = %d\n",
+          //             which_dims[i][0]+1, which_dims[i][1]+1);
+
+          tmp.setfield ("dim", which_dims[i][1]+1);
+
+          // Create entropy output 
+          Matrix entropy_out (EPSCOUNT, 3);
 
-              std::vector <double> histo (maxembed * dimension, 0.0);
-              start_box(series, which_dims, histo, maxembed, dimension, DELAY,
-                        length, epsi, Q);
-
-              if (Q != 1.0)
-                for (std::vector<double>::iterator it = histo.begin ();
-                     it != histo.end (); it++)
-                  *it = log(*it)/(1.0-Q);
-
-              histo_list.push_back (histo);
+          std::list<std::vector<double>>::const_iterator it_hist;
+          it_hist = histo_list.cbegin ();
+          for (octave_idx_type j=0;j<EPSCOUNT;j++)
+            {
+              if (i == 0)
+                {
+                  // old fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
+                  //             histo_el->hist[i],histo_el->hist[i]);
+                  entropy_out(j,0) = deps[j]*maxinterval;
+                  entropy_out(j,1) = (*it_hist)[i];
+                  entropy_out(j,2) = (*it_hist)[i];
+                }
+              else
+                {
+                  // old fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
+                  //             histo_el->hist[i],
+                  //             histo_el->hist[i]-histo_el->hist[i-1]);
+                  entropy_out(j,0) = deps[j]*maxinterval;
+                  entropy_out(j,1) = (*it_hist)[i];
+                  entropy_out(j,2) = (*it_hist)[i]
+                                     - (*it_hist)[i-1];
+                }
+              it_hist++;
             }
 
-          // Create and assign output
-          dim_vector dv (maxembed, dimension);
-          string_vector keys;
-          keys.append (std::string("dim"));
-          keys.append (std::string("entropy"));
-          octave_map output (dv, keys);
-
-          for (octave_idx_type i=0;i<maxembed*dimension;i++)
-            {
-              octave_scalar_map tmp (keys);
-              // old fprintf(fHq,"#component = %d embedding = %d\n",
-              //             which_dims[i][0]+1, which_dims[i][1]+1);
-
-              tmp.setfield ("dim", which_dims[i][1]+1);
-
-              // Create entropy output 
-              Matrix entropy_out (EPSCOUNT, 3);
+          tmp.setfield ("entropy",entropy_out);
 
-              std::list<std::vector<double>>::const_iterator it_hist;
-              it_hist = histo_list.cbegin ();
-              for (octave_idx_type j=0;j<EPSCOUNT;j++)
-                {
-                  if (i == 0)
-                    {
-                      // old fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
-                      //             histo_el->hist[i],histo_el->hist[i]);
-                      entropy_out(j,0) = deps[j]*maxinterval;
-                      entropy_out(j,1) = (*it_hist)[i];
-                      entropy_out(j,2) = (*it_hist)[i];
-                    }
-                  else
-                    {
-                      // old fprintf(fHq,"%e %e %e\n",deps[j]*maxinterval,
-                      //             histo_el->hist[i],
-                      //             histo_el->hist[i]-histo_el->hist[i-1]);
-                      entropy_out(j,0) = deps[j]*maxinterval;
-                      entropy_out(j,1) = (*it_hist)[i];
-                      entropy_out(j,2) = (*it_hist)[i]
-                                         - (*it_hist)[i-1];
-                    }
-                  it_hist++;
-                }
+          output.assign (idx_vector(which_dims[i][1]),
+                         idx_vector(which_dims[i][0]),
+                         tmp);
+        }
 
-              tmp.setfield ("entropy",entropy_out);
-
-              output.assign (idx_vector(which_dims[i][1]),
-                             idx_vector(which_dims[i][0]),
-                             tmp);
-            }
-
-
-          retval(0) = output;
-        }
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__c1__.cc
--- a/src/__c1__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__c1__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -85,65 +85,61 @@
       F77_INT iverb          = TO_F77_INT (verbose);
 
 
-      if (! error_state)
-        {
-
-          F77_INT lines_read   = TO_F77_INT (input.rows ()); //nmax in d1()
-          F77_INT columns_read = TO_F77_INT (input.columns ());
+      F77_INT lines_read   = TO_F77_INT (input.rows ()); //nmax in d1()
+      F77_INT columns_read = TO_F77_INT (input.columns ());
 
 
-          dim_vector dv (maxdim - mindim + 1, 1);
-          string_vector keys;
-          keys.append (std::string("dim"));
-          keys.append (std::string("c1"));
-          octave_map output (dv, keys);
+      dim_vector dv (maxdim - mindim + 1, 1);
+      string_vector keys;
+      keys.append (std::string("dim"));
+      keys.append (std::string("c1"));
+      octave_map output (dv, keys);
+
+      // Seed the rand() function for d1()
+      F77_XFCN (rand, RAND, (sqrt(seed)));
 
-          // Seed the rand() function for d1()
-          F77_XFCN (rand, RAND, (sqrt(seed)));
+      for (F77_INT m = mindim; m <= maxdim; m++)
+        {
+          octave_scalar_map tmp (keys);
+          tmp.setfield ("dim", m);
 
-          for (F77_INT m = mindim; m <= maxdim; m++)
+          // Creat c1 output
+          Matrix c1_out ((octave_idx_type) ((0 - log (1./(lines_read 
+                                                     -(m-1) * delay)) +
+                                            (log (2.) /resolution)) 
+                                                 / (log (2.) /resolution))
+                         , 2);
+
+          double pr = 0.0;
+          octave_idx_type current_row = 0;
+          for (double pl = log (1./(lines_read - (m-1)*delay));
+               pl <= 0.0; pl += log (2.) / resolution)
             {
-              octave_scalar_map tmp (keys);
-              tmp.setfield ("dim", m);
+              double pln = pl;
+              double rln;
 
-              // Creat c1 output
-              Matrix c1_out ((octave_idx_type) ((0 - log (1./(lines_read 
-                                                         -(m-1) * delay)) +
-                                                (log (2.) /resolution)) 
-                                                     / (log (2.) /resolution))
-                             , 2);
+              F77_XFCN (d1, D1,
+                        (lines_read, columns_read, lines_read,
+                         input.fortran_vec (), delay, m, cmin,
+                         pr, pln, rln, tmin, kmax, iverb));
 
-              double pr = 0.0;
-              octave_idx_type current_row = 0;
-              for (double pl = log (1./(lines_read - (m-1)*delay));
-                   pl <= 0.0; pl += log (2.) / resolution)
+              if (pln != pr)
                 {
-                  double pln = pl;
-                  double rln;
-
-                  F77_XFCN (d1, D1,
-                            (lines_read, columns_read, lines_read,
-                             input.fortran_vec (), delay, m, cmin,
-                             pr, pln, rln, tmin, kmax, iverb));
+                  pr = pln;
+                  c1_out(current_row,0) = exp (rln);
+                  c1_out(current_row,1) = exp (pln);
+                  current_row += 1;
+                }
 
-                  if (pln != pr)
-                    {
-                      pr = pln;
-                      c1_out(current_row,0) = exp (rln);
-                      c1_out(current_row,1) = exp (pln);
-                      current_row += 1;
-                    }
+            }
+          // Resize output
+          c1_out.resize (current_row, 2);
+          tmp.setfield ("c1", c1_out);
 
-                }
-              // Resize output
-              c1_out.resize (current_row, 2);
-              tmp.setfield ("c1", c1_out);
+          output.assign (idx_vector(m-mindim), tmp);
+        }
 
-              output.assign (idx_vector(m-mindim), tmp);
-            }
-
-          retval(0) = output;
-        }
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__d2__.cc
--- a/src/__d2__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__d2__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -365,251 +365,247 @@
       time_t lasttime;
       time(&lasttime);
 
-      if (! error_state)
+      bool imin_too_large = false;
+      bool pause_calc     = false;
+      // Calculate the outputs
+      for (octave_idx_type n = 1 + counter; n < nmax && !imin_too_large
+                                            && !pause_calc; n++)
         {
-
-          bool imin_too_large = false;
-          bool pause_calc     = false;
-          // Calculate the outputs
-          for (octave_idx_type n = 1 + counter; n < nmax && !imin_too_large
-                                                && !pause_calc; n++)
+          counter           += 1;
+          bool smaller       = 0;
+          octave_idx_type sn = scr[n-1];
+          double epsinv      = 1.0 / EPSMAX;
+          octave_idx_type x,y;
+          if (DIM > 1)
             {
-              counter           += 1;
-              bool smaller       = 0;
-              octave_idx_type sn = scr[n-1];
-              double epsinv      = 1.0 / EPSMAX;
-              octave_idx_type x,y;
-              if (DIM > 1)
-                {
-                  x=(octave_idx_type)(series[0][sn]*epsinv)&(NMAX-1);
-                  y=(octave_idx_type)(series[1][sn]*epsinv)&(NMAX-1);
-                }
-              else
-                {
-                  x=(octave_idx_type)(series[0][sn]*epsinv)&(NMAX-1);
-                  y=(octave_idx_type)(series[0][sn+DELAY]*epsinv)&(NMAX-1);
-                }
-              list[sn]=box[x][y];
-              box[x][y]=sn;
-              listc1[sn]=boxc1[x];
-              boxc1[x]=sn;
+              x=(octave_idx_type)(series[0][sn]*epsinv)&(NMAX-1);
+              y=(octave_idx_type)(series[1][sn]*epsinv)&(NMAX-1);
+            }
+          else
+            {
+              x=(octave_idx_type)(series[0][sn]*epsinv)&(NMAX-1);
+              y=(octave_idx_type)(series[0][sn+DELAY]*epsinv)&(NMAX-1);
+            }
+          list[sn]=box[x][y];
+          box[x][y]=sn;
+          listc1[sn]=boxc1[x];
+          boxc1[x]=sn;
 
-              octave_idx_type i=imin;
-              while (found[maxembed][i] >= MAXFOUND)
-                {
-                  smaller=1;
-                  if (++i > (HOWOFTEN-1))
-                    break;
-                }
-              if (smaller)
-                {
-                  imin=i;
-                  if (imin <= (HOWOFTEN-1))
-                    {
-                      EPSMAX        = epsm[imin];
-                      double epsinv = 1.0/EPSMAX;
-                      for (octave_idx_type i1=0;i1<NMAX;i1++)
-                        {
-                          boxc1[i1]= -1;
-                          for (octave_idx_type j1=0;j1<NMAX;j1++)
-                            box[i1][j1]= -1;
-                        }
-                      for (octave_idx_type i1=0;i1<n;i1++)
-                        {
-                          sn=scr[i1];
-                          octave_idx_type x,y;
-                          if (DIM > 1)
-                            {
-                              x=(octave_idx_type)(series[0][sn]*epsinv)
-                                &(NMAX-1);
-                              y=(octave_idx_type)(series[1][sn]*epsinv)
-                                &(NMAX-1);
-                            }
-                          else
-                            {
-                              x=(octave_idx_type)(series[0][sn]*epsinv)
-                                &(NMAX-1);
-                              y=(octave_idx_type)(series[0][sn+DELAY]*epsinv)
-                                &(NMAX-1);
-                            }
-                          list[sn]=box[x][y];
-                          box[x][y]=sn;
-                          listc1[sn]=boxc1[x];
-                          boxc1[x]=sn;
-                        }
-                    }
-                }
-
+          octave_idx_type i=imin;
+          while (found[maxembed][i] >= MAXFOUND)
+            {
+              smaller=1;
+              if (++i > (HOWOFTEN-1))
+                break;
+            }
+          if (smaller)
+            {
+              imin=i;
               if (imin <= (HOWOFTEN-1))
                 {
-                  octave_idx_type lnorm=n;
-                  if (MINDIST > 0)
+                  EPSMAX        = epsm[imin];
+                  double epsinv = 1.0/EPSMAX;
+                  for (octave_idx_type i1=0;i1<NMAX;i1++)
                     {
-                      octave_idx_type sn=scr[n];
-                      octave_idx_type n1=(sn-MINDIST>=0)?sn-MINDIST:0;
-                      octave_idx_type n2=(sn+MINDIST<length-(EMBED-1)*DELAY)
-                                         ?sn+MINDIST:length-(EMBED-1)*DELAY-1;
-                      for (octave_idx_type i1=n1;i1<=n2;i1++)
-                        if ((oscr[i1] < n))
-                          lnorm--;
+                      boxc1[i1]= -1;
+                      for (octave_idx_type j1=0;j1<NMAX;j1++)
+                        box[i1][j1]= -1;
                     }
-
-                  if (EMBED*DIM > 1)
-                    make_c2_dim(series, found, list, box, DIM, imin, EPSMAX1,
-                                EPSMAX, EPSMIN, EMBED, DELAY, MINDIST,
-                                HOWOFTEN, scr[n]);
-                  make_c2_1(series, found, listc1, boxc1, imin, EPSMAX1,
-                            EPSMAX, EPSMIN, MINDIST, HOWOFTEN, scr[n]);
-                  for (octave_idx_type i=imin;i<HOWOFTEN;i++)
-                    norm[i] += (double)(lnorm);
-                }
-
-
-              // If any of the below occurs: pause or end.
-              if (((time(NULL)-lasttime) > WHEN) || (n == (nmax-1)) || 
-                  (imin > (HOWOFTEN-1)) || (counter % it_pause == 0))
-                {
-                  time(&lasttime);
-
-                  if (imin > (HOWOFTEN-1))
+                  for (octave_idx_type i1=0;i1<n;i1++)
                     {
-                      // old exit(0);
-                      imin_too_large = true;
+                      sn=scr[i1];
+                      octave_idx_type x,y;
+                      if (DIM > 1)
+                        {
+                          x=(octave_idx_type)(series[0][sn]*epsinv)
+                            &(NMAX-1);
+                          y=(octave_idx_type)(series[1][sn]*epsinv)
+                            &(NMAX-1);
+                        }
+                      else
+                        {
+                          x=(octave_idx_type)(series[0][sn]*epsinv)
+                            &(NMAX-1);
+                          y=(octave_idx_type)(series[0][sn+DELAY]*epsinv)
+                            &(NMAX-1);
+                        }
+                      list[sn]=box[x][y];
+                      box[x][y]=sn;
+                      listc1[sn]=boxc1[x];
+                      boxc1[x]=sn;
                     }
-                  pause_calc = true;
                 }
             }
 
-          // Create vars output
-          octave_scalar_map vars;
-
-
-          // Create vars output
-          // old fprintf(fstat,"Center points treated so far= %ld\n",n);
-          vars.setfield ("treated", counter);
-          // old fprintf(fstat,"Maximum epsilon in the moment= %e\n",
-          //             epsm[imin]);
-          vars.setfield ("eps", epsm[imin]);
-
-          if (counter < nmax - 1 && imin_too_large == false)
+          if (imin <= (HOWOFTEN-1))
             {
-              vars.setfield ("counter", counter);
-              vars.setfield ("found", found_matrix);
-              vars.setfield ("norm", norm_matrix);
-              vars.setfield ("boxc1", boxc1_matrix);
-              vars.setfield ("box", box_matrix);
-              vars.setfield ("list", list_matrix);
-              vars.setfield ("listc1", listc1_matrix);
-              vars.setfield ("imin", imin);
-              vars.setfield ("EPSMAX1",EPSMAX1);
-              vars.setfield ("EPSMAX", EPSMAX);
-              vars.setfield ("EPSMIN", EPSMIN);
-            }
-
-          // Create values output
-          dim_vector dv (DIM * EMBED, 1);
-          string_vector keys;
-          keys.append (std::string("dim"));
-          keys.append (std::string("c2"));
-          keys.append (std::string("d2"));
-          keys.append (std::string("h2"));
-          octave_map values (dv, keys);
-
-          for (octave_idx_type i=0;i<DIM*EMBED;i++)
-            {
-
-              octave_scalar_map tmp (keys);
-
-              // old fprintf(fout,"#dim= %ld\n",i+1);
-              tmp.setfield ("dim", i+1);
-
-              // Allocate d2 output
-              Matrix d2_out (HOWOFTEN - 1, 2);
-              octave_idx_type d2_row = 0;
-
-              // Allocate h2 output
-              Matrix h2_out (HOWOFTEN, 2);
-              octave_idx_type h2_row = 0;
-
-              // Allocate c2 output
-              Matrix c2_out (HOWOFTEN, 2);
-              octave_idx_type c2_row = 0;
-
-              double eps = EPSMAX1 * epsfactor;
-
-              for (octave_idx_type j=0;j<HOWOFTEN;j++)
+              octave_idx_type lnorm=n;
+              if (MINDIST > 0)
                 {
-                  eps /= epsfactor;
+                  octave_idx_type sn=scr[n];
+                  octave_idx_type n1=(sn-MINDIST>=0)?sn-MINDIST:0;
+                  octave_idx_type n2=(sn+MINDIST<length-(EMBED-1)*DELAY)
+                                     ?sn+MINDIST:length-(EMBED-1)*DELAY-1;
+                  for (octave_idx_type i1=n1;i1<=n2;i1++)
+                    if ((oscr[i1] < n))
+                      lnorm--;
+                }
 
-                  // Calculate d2 output
-                  if ((j > 0) && (found[i][j] > 0.0)
-                      && (found[i][j-1] > 0.0))
-                    {
-                      // old fprintf(fout,"%e %e\n",eps,
-                      //      log(found[i][j-1]/found[i][j]/norm[j-1]
-                      //          *norm[j])
-                      //      /log(epsfactor));
-                      d2_out(d2_row,0) = eps;
-                      d2_out(d2_row,1) = log(found[i][j-1]/found[i][j]
-                                             /norm[j-1]*norm[j])
-                                         /log(epsfactor);
-                      d2_row          += 1;
-                    }
-
-                  // Calculate h2 output
-                  if (i < 1)
-                    {
-                      if (found[0][j] > 0.0)
-                        {
-                          // old fprintf(fout,"%e %e\n",eps,
-                          //             -log(found[0][j]/norm[j]));
-                          h2_out(h2_row,0) = eps;
-                          h2_out(h2_row,1) = -log(found[0][j]/norm[j]);
-                          h2_row          += 1;
-                        }
-                    }
-                  else
-                    {
-                      if ((found[i-1][j] > 0.0) && (found[i][j] > 0.0))
-                        {
-                        // old fprintf(fout,"%e %e\n",eps,
-                        //             log(found[i-1][j]/found[i][j]));
-                          h2_out(h2_row,0) = eps;
-                          h2_out(h2_row,1) = log(found[i-1][j]
-                                                 /found[i][j]);
-                          h2_row          += 1;
-                        }
-                    }
-
-                  // Calculate c2 output
-                  if (norm[j] > 0.0)
-                    {
-                      // old fprintf(fout,"%e %e\n",eps,
-                      //             found[i][j]/norm[j]);
-                      c2_out(c2_row,0) = eps;
-                      c2_out(c2_row,1) = found[i][j]/norm[j];
-                      c2_row          += 1;
-                    }
-                }
-              // Prepare d2 output
-              d2_out.resize (d2_row, 2);
-              tmp.setfield ("d2", d2_out);
-              // Prepare h2 output
-              h2_out.resize (h2_row, 2);
-              tmp.setfield ("h2", h2_out);
-              // Prepare c2 output
-              c2_out.resize (c2_row, 2);
-              tmp.setfield ("c2", c2_out);
-
-              values.assign (idx_vector(i), tmp);
+              if (EMBED*DIM > 1)
+                make_c2_dim(series, found, list, box, DIM, imin, EPSMAX1,
+                            EPSMAX, EPSMIN, EMBED, DELAY, MINDIST,
+                            HOWOFTEN, scr[n]);
+              make_c2_1(series, found, listc1, boxc1, imin, EPSMAX1,
+                        EPSMAX, EPSMIN, MINDIST, HOWOFTEN, scr[n]);
+              for (octave_idx_type i=imin;i<HOWOFTEN;i++)
+                norm[i] += (double)(lnorm);
             }
 
 
-          // Assign outputs
-          retval(0) = values;
-          retval(1) = vars;
+          // If any of the below occurs: pause or end.
+          if (((time(NULL)-lasttime) > WHEN) || (n == (nmax-1)) || 
+              (imin > (HOWOFTEN-1)) || (counter % it_pause == 0))
+            {
+              time(&lasttime);
+
+              if (imin > (HOWOFTEN-1))
+                {
+                  // old exit(0);
+                  imin_too_large = true;
+                }
+              pause_calc = true;
+            }
+        }
+
+      // Create vars output
+      octave_scalar_map vars;
+
+
+      // Create vars output
+      // old fprintf(fstat,"Center points treated so far= %ld\n",n);
+      vars.setfield ("treated", counter);
+      // old fprintf(fstat,"Maximum epsilon in the moment= %e\n",
+      //             epsm[imin]);
+      vars.setfield ("eps", epsm[imin]);
+
+      if (counter < nmax - 1 && imin_too_large == false)
+        {
+          vars.setfield ("counter", counter);
+          vars.setfield ("found", found_matrix);
+          vars.setfield ("norm", norm_matrix);
+          vars.setfield ("boxc1", boxc1_matrix);
+          vars.setfield ("box", box_matrix);
+          vars.setfield ("list", list_matrix);
+          vars.setfield ("listc1", listc1_matrix);
+          vars.setfield ("imin", imin);
+          vars.setfield ("EPSMAX1",EPSMAX1);
+          vars.setfield ("EPSMAX", EPSMAX);
+          vars.setfield ("EPSMIN", EPSMIN);
         }
 
+      // Create values output
+      dim_vector dv (DIM * EMBED, 1);
+      string_vector keys;
+      keys.append (std::string("dim"));
+      keys.append (std::string("c2"));
+      keys.append (std::string("d2"));
+      keys.append (std::string("h2"));
+      octave_map values (dv, keys);
+
+      for (octave_idx_type i=0;i<DIM*EMBED;i++)
+        {
+
+          octave_scalar_map tmp (keys);
+
+          // old fprintf(fout,"#dim= %ld\n",i+1);
+          tmp.setfield ("dim", i+1);
+
+          // Allocate d2 output
+          Matrix d2_out (HOWOFTEN - 1, 2);
+          octave_idx_type d2_row = 0;
+
+          // Allocate h2 output
+          Matrix h2_out (HOWOFTEN, 2);
+          octave_idx_type h2_row = 0;
+
+          // Allocate c2 output
+          Matrix c2_out (HOWOFTEN, 2);
+          octave_idx_type c2_row = 0;
+
+          double eps = EPSMAX1 * epsfactor;
+
+          for (octave_idx_type j=0;j<HOWOFTEN;j++)
+            {
+              eps /= epsfactor;
+
+              // Calculate d2 output
+              if ((j > 0) && (found[i][j] > 0.0)
+                  && (found[i][j-1] > 0.0))
+                {
+                  // old fprintf(fout,"%e %e\n",eps,
+                  //      log(found[i][j-1]/found[i][j]/norm[j-1]
+                  //          *norm[j])
+                  //      /log(epsfactor));
+                  d2_out(d2_row,0) = eps;
+                  d2_out(d2_row,1) = log(found[i][j-1]/found[i][j]
+                                         /norm[j-1]*norm[j])
+                                     /log(epsfactor);
+                  d2_row          += 1;
+                }
+
+              // Calculate h2 output
+              if (i < 1)
+                {
+                  if (found[0][j] > 0.0)
+                    {
+                      // old fprintf(fout,"%e %e\n",eps,
+                      //             -log(found[0][j]/norm[j]));
+                      h2_out(h2_row,0) = eps;
+                      h2_out(h2_row,1) = -log(found[0][j]/norm[j]);
+                      h2_row          += 1;
+                    }
+                }
+              else
+                {
+                  if ((found[i-1][j] > 0.0) && (found[i][j] > 0.0))
+                    {
+                    // old fprintf(fout,"%e %e\n",eps,
+                    //             log(found[i-1][j]/found[i][j]));
+                      h2_out(h2_row,0) = eps;
+                      h2_out(h2_row,1) = log(found[i-1][j]
+                                             /found[i][j]);
+                      h2_row          += 1;
+                    }
+                }
+
+              // Calculate c2 output
+              if (norm[j] > 0.0)
+                {
+                  // old fprintf(fout,"%e %e\n",eps,
+                  //             found[i][j]/norm[j]);
+                  c2_out(c2_row,0) = eps;
+                  c2_out(c2_row,1) = found[i][j]/norm[j];
+                  c2_row          += 1;
+                }
+            }
+          // Prepare d2 output
+          d2_out.resize (d2_row, 2);
+          tmp.setfield ("d2", d2_out);
+          // Prepare h2 output
+          h2_out.resize (h2_row, 2);
+          tmp.setfield ("h2", h2_out);
+          // Prepare c2 output
+          c2_out.resize (c2_row, 2);
+          tmp.setfield ("c2", c2_out);
+
+          values.assign (idx_vector(i), tmp);
+        }
+
+
+      // Assign outputs
+      retval(0) = values;
+      retval(1) = vars;
     }
+
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__delay__.cc
--- a/src/__delay__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__delay__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -60,56 +60,52 @@
 
       OCTAVE_LOCAL_BUFFER (octave_idx_type, inddelay, alldim);
 
-      if (! error_state)
-        {
-          octave_idx_type rundel=0;
-          octave_idx_type runmdel=0;
-
-          unsigned int delsum;
-          for (octave_idx_type i = 0; i < indim; i++)
-            {
-              delsum             = 0;
-              inddelay[rundel++] = delsum;
-
-              for (octave_idx_type j = 1; j < formatdelay(i); j++)
-                {
-                  delsum            += delaylist(runmdel++);
-                  inddelay[rundel++] = delsum;
-                }
-            }
-
-          octave_idx_type maxemb = 0;
-          for (octave_idx_type i = 0; i < alldim; i++)
-            maxemb = (maxemb < inddelay[i])? inddelay[i] : maxemb;
+      octave_idx_type rundel=0;
+      octave_idx_type runmdel=0;
 
-          octave_idx_type outdim = 0;
-          for (octave_idx_type i = 0; i < indim; i++)
-             outdim += formatdelay(i);
-
-          octave_idx_type out_rows = (length > maxemb) ? length - maxemb : 0;
-
-          Matrix series (out_rows, outdim);
-          unsigned int embsum;
-          for (octave_idx_type i = maxemb; i < length; i++)
-            {
-              rundel = 0;
-              embsum = 0;
+      unsigned int delsum;
+      for (octave_idx_type i = 0; i < indim; i++)
+        {
+          delsum             = 0;
+          inddelay[rundel++] = delsum;
 
-              for (octave_idx_type j = 0; j < indim; j++)
-                {
-                  octave_idx_type emb = formatdelay(j);
-
-                  for (octave_idx_type k = 0; k < emb; k++)
-                    series(i-maxemb, embsum+k) = data(i-inddelay[rundel++], j);
-
-                    // previously fprintf(stdout,"%e ",series[j][i-inddelay[rundel++]]);
-                  embsum += emb;
-                }
-              // previously fprintf(stdout,"\n");
+          for (octave_idx_type j = 1; j < formatdelay(i); j++)
+            {
+              delsum            += delaylist(runmdel++);
+              inddelay[rundel++] = delsum;
             }
-          retval(0) = series;
         }
 
+      octave_idx_type maxemb = 0;
+      for (octave_idx_type i = 0; i < alldim; i++)
+        maxemb = (maxemb < inddelay[i])? inddelay[i] : maxemb;
+
+      octave_idx_type outdim = 0;
+      for (octave_idx_type i = 0; i < indim; i++)
+         outdim += formatdelay(i);
+
+      octave_idx_type out_rows = (length > maxemb) ? length - maxemb : 0;
+
+      Matrix series (out_rows, outdim);
+      unsigned int embsum;
+      for (octave_idx_type i = maxemb; i < length; i++)
+        {
+          rundel = 0;
+          embsum = 0;
+
+          for (octave_idx_type j = 0; j < indim; j++)
+            {
+              octave_idx_type emb = formatdelay(j);
+
+              for (octave_idx_type k = 0; k < emb; k++)
+                series(i-maxemb, embsum+k) = data(i-inddelay[rundel++], j);
+
+              // previously fprintf(stdout,"%e ",series[j][i-inddelay[rundel++]]);
+              embsum += emb;
+            }
+          // previously fprintf(stdout,"\n");
+        }
+      retval(0) = series;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__false_nearest__.cc
--- a/src/__false_nearest__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__false_nearest__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -179,78 +179,75 @@
       check_alloc(vcomp=(unsigned int*)malloc(sizeof(int)*(maxdim)));
       check_alloc(vemb=(unsigned int*)malloc(sizeof(int)*(maxdim)));
 
-      if ( ! error_state)
-        {
-          for (i=0;i<maxdim;i++) {
-            if (comp == 1) {
-              vcomp[i]=0;
-              vemb[i]=i*delay;
-            }
-            else {
-              vcomp[i]=i%comp;
-              vemb[i]=(i/comp)*delay;
-            }
-          }
+      for (i=0;i<maxdim;i++) {
+        if (comp == 1) {
+          vcomp[i]=0;
+          vemb[i]=i*delay;
+        }
+        else {
+          vcomp[i]=i%comp;
+          vemb[i]=(i/comp)*delay;
+        }
+      }
 
-          // Create output matrix
-          Matrix output (maxemb-minemb+1, 4);
+      // Create output matrix
+      Matrix output (maxemb-minemb+1, 4);
 
-          // Compute output
-          for (emb=minemb;emb<=maxemb;emb++) 
-            {
-              dim=emb*comp-1;
-              epsilon=eps0;
-              toolarge=0;
-              alldone=0;
-              donesofar=0;
-              aveps=0.0;
-              vareps=0.0;
-              for (i=0;i<length;i++)
-                nearest[i]=0;
-              if (verbosity)
-                octave_stdout << "Start for dimension=" << dim+1 << "\n";
-              while (!alldone && (epsilon < 2.*varianz/rt)) {
-                alldone=1;
-                mmb(input, vcomp[dim],vemb[dim],epsilon);
-                for (i=0;i<length-maxemb*delay;i++)
-                if (!nearest[i]) {
-                  nearest[i]=find_nearest(input, i,dim,epsilon);
-                  alldone &= nearest[i];
-                  donesofar += (unsigned long)nearest[i];
-                }
-                if (verbosity)
-                octave_stdout << "Found " << donesofar << " up to epsilon=" << \
-                              epsilon*inter << "\n";
-                epsilon*=sqrt(2.0);
-                if (!donesofar)
-                eps0=epsilon;
-              }
-              if (donesofar == 0) {
-                error_with_id ("Octave:tisean", "Not enough points found");
-              }
-              aveps *= (1./(double)donesofar);
-              vareps *= (1./(double)donesofar);
+      // Compute output
+      for (emb=minemb;emb<=maxemb;emb++) 
+        {
+          dim=emb*comp-1;
+          epsilon=eps0;
+          toolarge=0;
+          alldone=0;
+          donesofar=0;
+          aveps=0.0;
+          vareps=0.0;
+          for (i=0;i<length;i++)
+            nearest[i]=0;
+          if (verbosity)
+            octave_stdout << "Start for dimension=" << dim+1 << "\n";
+          while (!alldone && (epsilon < 2.*varianz/rt)) {
+            alldone=1;
+            mmb(input, vcomp[dim],vemb[dim],epsilon);
+            for (i=0;i<length-maxemb*delay;i++)
+            if (!nearest[i]) {
+              nearest[i]=find_nearest(input, i,dim,epsilon);
+              alldone &= nearest[i];
+              donesofar += (unsigned long)nearest[i];
+            }
+            if (verbosity)
+            octave_stdout << "Found " << donesofar << " up to epsilon=" << \
+                          epsilon*inter << "\n";
+            epsilon*=sqrt(2.0);
+            if (!donesofar)
+            eps0=epsilon;
+          }
+          if (donesofar == 0) {
+            error_with_id ("Octave:tisean", "Not enough points found");
+          }
+          aveps *= (1./(double)donesofar);
+          vareps *= (1./(double)donesofar);
 
 
-          //  old fprintf(file,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
-          //              aveps*inter,sqrt(vareps)*inter);
-              int id = emb-minemb;
-              output(id,0) = dim + 1;
-              output(id,1) = (double)toolarge/(double)donesofar;
-              output(id,2) = aveps*inter;
-              output(id,3) = sqrt(vareps)*inter;
-            }
+      //  old fprintf(file,"%u %e %e %e\n",dim+1,(double)toolarge/(double)donesofar,
+      //              aveps*inter,sqrt(vareps)*inter);
+          int id = emb-minemb;
+          output(id,0) = dim + 1;
+          output(id,1) = (double)toolarge/(double)donesofar;
+          output(id,2) = aveps*inter;
+          output(id,3) = sqrt(vareps)*inter;
+        }
 
-          delete[] series;
-          delete[] list;
-          delete[] nearest;
-          for (i=0;i<BOX;i++)
-            delete[] box[i];
-          delete[] box;
+      delete[] series;
+      delete[] list;
+      delete[] nearest;
+      for (i=0;i<BOX;i++)
+        delete[] box[i];
+      delete[] box;
 
-          for (i = 0; i < 4; i++)
-            retval(i) = output.column(i);
-        }
+      for (i = 0; i < 4; i++)
+        retval(i) = output.column(i);
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__ghkss__.cc
--- a/src/__ghkss__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__ghkss__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -203,42 +203,39 @@
 
   OCTAVE_LOCAL_BUFFER (double, hav, comp);
   OCTAVE_LOCAL_BUFFER (double, hsigma, comp);
-  if ( ! error_state)
-    {
-      for (j=0;j<comp;j++)
-        hav[j]=hsigma[j]=0.0;
+  for (j=0;j<comp;j++)
+    hav[j]=hsigma[j]=0.0;
 
-      for (i=0;i<length;i++)
-        for (j=0;j<comp;j++) {
-          hav[j] += (help=delta[j][i]);
-          hsigma[j] += help*help;
-        }
+  for (i=0;i<length;i++)
+    for (j=0;j<comp;j++) {
+      hav[j] += (help=delta[j][i]);
+      hsigma[j] += help*help;
+    }
 
-      for (j=0;j<comp;j++) {
-        hav[j] /= length;
-        hsigma[j]=sqrt(fabs(hsigma[j]/length-hav[j]*hav[j]));
-      }
-      if (verbosity) {
-        for (i=0;i<comp;i++) {
-          octave_stdout << "Average shift of component " << i+1 << " = "
-                        << hav[i]*d_max[i] << "\n";
-          octave_stdout << "Average rms correction of comp. " << i+1 << " = "
-                        << hsigma[i]*d_max[i] << "\n\n";
-        }
-      }
-      for (i=0;i<length;i++)
-        for (j=0;j<comp;j++)
-          series(i,j) -= delta[j][i];
+  for (j=0;j<comp;j++) {
+    hav[j] /= length;
+    hsigma[j]=sqrt(fabs(hsigma[j]/length-hav[j]*hav[j]));
+  }
+  if (verbosity) {
+    for (i=0;i<comp;i++) {
+      octave_stdout << "Average shift of component " << i+1 << " = "
+                    << hav[i]*d_max[i] << "\n";
+      octave_stdout << "Average rms correction of comp. " << i+1 << " = "
+                    << hsigma[i]*d_max[i] << "\n\n";
+    }
+  }
+  for (i=0;i<length;i++)
+    for (j=0;j<comp;j++)
+      series(i,j) -= delta[j][i];
 
-      if (resize_eps) {
-        mineps /= epsfac;
-        if (verbosity)
-          octave_stdout << "Reset minimal neighbourhood size to "
-                        << mineps*d_max_max << "\n";
-      }
+  if (resize_eps) {
+    mineps /= epsfac;
+    if (verbosity)
+      octave_stdout << "Reset minimal neighbourhood size to "
+                    << mineps*d_max_max << "\n";
+  }
 
-      resize_eps=0;
-    }
+  resize_eps=0;
 }
 
 DEFUN_DLD (__ghkss__, args, , HELPTEXT)
@@ -338,92 +335,90 @@
         mat[i]=(double*)(matarray+dim*i);
       check_alloc(hser=(double**)malloc(sizeof(double*)*comp));
 
-      if ( ! error_state)
-        {
-          // Create output matrix
-          Matrix output (length, comp);
+      // Create output matrix
+      Matrix output (length, comp);
 
-          for (i=0;i<dim;i++) {
-            index_comp[i]=i%comp;
-            index_embed[i]=(i/comp)*delay;
-          }
+      for (i=0;i<dim;i++) {
+        index_comp[i]=i%comp;
+        index_embed[i]=(i/comp)*delay;
+      }
 
-          // Calculate the noise reduction
-          resize_eps=0;
-          for (iter=1;iter<=iterations;iter++) 
-            {
-              for (i=0;i<length;i++) {
-                ok[i]=0;
-                for (j=0;j<dim;j++)
-                corr[i][j]=0.0;
-                for (j=0;j<comp;j++)
-                delta[j][i]=0.0;
+      // Calculate the noise reduction
+      resize_eps=0;
+      for (iter=1;iter<=iterations;iter++) 
+        {
+          for (i=0;i<length;i++) {
+            ok[i]=0;
+            for (j=0;j<dim;j++)
+            corr[i][j]=0.0;
+            for (j=0;j<comp;j++)
+            delta[j][i]=0.0;
+          }
+          epsilon=mineps;
+          all_done=0;
+          epscount=1;
+          allfound=0;
+          if (verbosity)
+            octave_stdout << "Starting iteration " << iter << "\n";
+          while(!all_done) {
+            mmb(input, epsilon);
+            all_done=1;
+            for (n=emb_offset;n<length;n++)
+            if (!ok[n]) {
+              nfound=fmn(input,n,epsilon);
+              if (nfound >= minn) {
+                make_correction(input,n,nfound);
+                ok[n]=epscount;
+                if (epscount == 1)
+                  resize_eps=1;
+                allfound++;
               }
-              epsilon=mineps;
-              all_done=0;
-              epscount=1;
-              allfound=0;
-              if (verbosity)
-                octave_stdout << "Starting iteration " << iter << "\n";
-              while(!all_done) {
-                mmb(input, epsilon);
-                all_done=1;
-                for (n=emb_offset;n<length;n++)
-                if (!ok[n]) {
-                  nfound=fmn(input,n,epsilon);
-                  if (nfound >= minn) {
-                    make_correction(input,n,nfound);
-                    ok[n]=epscount;
-                    if (epscount == 1)
-                      resize_eps=1;
-                    allfound++;
+              else
+                all_done=0;
+            }
+            if (verbosity)
+              octave_stdout << "Corrected " << allfound << " points with epsilon= "
+                            << epsilon*d_max_max << "\n";
+            if (std::isinf(epsilon*d_max_max))
+              {
+                error_with_id ("Octave:tisean", "cannot reduce noise on input data");
+                return retval;
+              }
+            epsilon *= epsfac;
+            epscount++;
+          }
+          if (verbosity)
+            octave_stdout << "Start evaluating the trend\n";
+
+          epsilon=mineps;
+          allfound=0;
+          for (i=1;i<epscount;i++) {
+            mmb(input,epsilon);
+            for (n=emb_offset;n<length;n++)
+            if (ok[n] == i) {
+              nfound=fmn(input,n,epsilon);
+              handle_trend(n,nfound);
+              allfound++;
+            }
+            if (verbosity)
+            octave_stdout << "Trend subtracted for " << allfound << " points with epsilon= " 
+                          << epsilon*d_max_max << "\n";
+            epsilon *= epsfac;
+          }
+          set_correction(input);
+
+          if (iter == iterations)
+            for (i=0;i<length;i++) 
+              {
+                for (j=0;j<comp;j++) 
+                  {
+                  // old  fprintf(file,"%e ",series[j][i]*d_max[j]+d_min[j]);
+                    output(i,j) = input(i,j)*d_max[j]+d_min[j];
                   }
-                  else
-                    all_done=0;
-                }
-                if (verbosity)
-                  octave_stdout << "Corrected " << allfound << " points with epsilon= "
-                                << epsilon*d_max_max << "\n";
-                if (std::isinf(epsilon*d_max_max))
-                  {
-                    error_with_id ("Octave:tisean", "cannot reduce noise on input data");
-                    return retval;
-                  }
-                epsilon *= epsfac;
-                epscount++;
               }
-              if (verbosity)
-                octave_stdout << "Start evaluating the trend\n";
+        }
+      retval(0) = output;
 
-              epsilon=mineps;
-              allfound=0;
-              for (i=1;i<epscount;i++) {
-                mmb(input,epsilon);
-                for (n=emb_offset;n<length;n++)
-                if (ok[n] == i) {
-                  nfound=fmn(input,n,epsilon);
-                  handle_trend(n,nfound);
-                  allfound++;
-                }
-                if (verbosity)
-                octave_stdout << "Trend subtracted for " << allfound << " points with epsilon= " 
-                              << epsilon*d_max_max << "\n";
-                epsilon *= epsfac;
-              }
-              set_correction(input);
-
-              if (iter == iterations)
-                for (i=0;i<length;i++) 
-                  {
-                    for (j=0;j<comp;j++) 
-                      {
-                      // old  fprintf(file,"%e ",series[j][i]*d_max[j]+d_min[j]);
-                        output(i,j) = input(i,j)*d_max[j]+d_min[j];
-                      }
-                  }
-            }
-          retval(0) = output;
-        }
     // Deallocate of all the memory
       delete[] d_min;
       delete[] d_max;
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lfo_ar__.cc
--- a/src/__lfo_ar__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lfo_ar__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -168,12 +168,6 @@
     Matrix solved_vec = mat.solve(vec);
     double *solved_vec_arr = solved_vec.fortran_vec ();
 
-    // If errors were raised, there is no sense in countinueing
-    if (error_state)
-      {
-        return ;
-      }
-
     double cast=foreav[i];
     for (octave_idx_type j=0;j<dim;j++) {
 
@@ -262,110 +256,101 @@
 
       octave_idx_type clength=(CLENGTH <= LENGTH) ? CLENGTH-STEP : LENGTH-STEP;
 
-      if ( ! error_state)
-        {
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
 
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
-
-          // Estimate maximum possible output size
-          octave_idx_type output_rows = (octave_idx_type)
-                                        ((log(EPS1) - log(EPS0)) / log (EPSF));
-          output_rows += 2;
+      // Estimate maximum possible output size
+      octave_idx_type output_rows = (octave_idx_type)
+                                    ((log(EPS1) - log(EPS0)) / log (EPSF));
+      output_rows += 2;
 
-          // Create output
-          Matrix output (output_rows,dim+4);
-          octave_idx_type count = 0;
+      // Create output
+      Matrix output (output_rows,dim+4);
+      octave_idx_type count = 0;
 
-          if (verbose)
-            printf("\nStarting new dataset\n\n");
+      if (verbose)
+        printf("\nStarting new dataset\n\n");
 
-          for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) 
+      for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) 
+        {
+          if (verbose)
             {
-              if (verbose)
-                {
-                  printf ("For epsilon = %e, the count = %lu\n",
-                          epsilon, (unsigned long)count);
-                  fflush (stdout);
-                }
+              printf ("For epsilon = %e, the count = %lu\n",
+                      epsilon, (unsigned long)count);
+              fflush (stdout);
+            }
 
-              long pfound=0;
-              for (octave_idx_type i=0;i<dim;i++)
-                error_array[i]=hrms[i]=hav[i]=0.0;
-              double avfound=0.0;
-              make_multi_box(series,box,list,LENGTH-STEP,NMAX,dim,
-                             embed,delay,epsilon);
-              for (octave_idx_type i=(embed-1)*delay;i<clength;i++) 
-                {
-                  for (octave_idx_type j=0;j<dim;j++)
-                    hser[j] = series[j] + i;
+          long pfound=0;
+          for (octave_idx_type i=0;i<dim;i++)
+            error_array[i]=hrms[i]=hav[i]=0.0;
+          double avfound=0.0;
+          make_multi_box(series,box,list,LENGTH-STEP,NMAX,dim,
+                         embed,delay,epsilon);
+          for (octave_idx_type i=(embed-1)*delay;i<clength;i++) 
+            {
+              for (octave_idx_type j=0;j<dim;j++)
+                hser[j] = series[j] + i;
 
-                  octave_idx_type actfound;
-                  actfound=find_multi_neighbors(series,box,list,hser,NMAX,
-                                                dim,embed,delay, epsilon,
-                                                hfound);
-                  actfound=exclude_interval(actfound,i-causal+1,
-                                            i+causal+(embed-1) *delay-1,
-                                            hfound,found);
-                  if (actfound > 2*(dim*embed+1))
+              octave_idx_type actfound;
+              actfound=find_multi_neighbors(series,box,list,hser,NMAX,
+                                            dim,embed,delay, epsilon,
+                                            hfound);
+              actfound=exclude_interval(actfound,i-causal+1,
+                                        i+causal+(embed-1) *delay-1,
+                                        hfound,found);
+              if (actfound > 2*(dim*embed+1))
+                {
+                  make_fit (series, found, error_array,
+                            i,dim, embed, delay, STEP, actfound);
+                  // Checking if the fit was correct
+                  pfound++;
+                  avfound += (double)(actfound-1);
+                  for (octave_idx_type j=0;j<dim;j++) 
                     {
-                      make_fit (series, found, error_array,
-                                i,dim, embed, delay, STEP, actfound);
-                      // Checking if the fit was correct
-                      // If any errors were raised: end function
-                      if (error_state)
-                        {
-                          return retval;
-                        }
-                      pfound++;
-                      avfound += (double)(actfound-1);
-                      for (octave_idx_type j=0;j<dim;j++) 
-                        {
-                          hrms[j] += series[j][i+STEP] * series[j][i+STEP];
-                          hav[j] += series[j][i+STEP];
-                        }
+                      hrms[j] += series[j][i+STEP] * series[j][i+STEP];
+                      hav[j] += series[j][i+STEP];
                     }
                 }
-              if (pfound > 1) 
+            }
+          if (pfound > 1) 
+            {
+              double sumerror=0.0;
+              for (octave_idx_type j=0;j<dim;j++)
+                {
+                  hav[j] /= pfound;
+                  hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]*pfound
+                                    /(pfound-1)));
+                  error_array[j]=sqrt(error_array[j]/pfound)/hrms[j];
+                  sumerror += error_array[j];
+                }
+         
+            
+          // old fprintf(stdout,"%e %e ",epsilon*interval,
+          //                          sumerror/(double)dim);
+              output(count, 0) = epsilon*interval;
+              output(count, 1) = sumerror/(double)dim;
+              for (octave_idx_type j=0;j<dim;j++)
                 {
-                  double sumerror=0.0;
-                  for (octave_idx_type j=0;j<dim;j++)
-                    {
-                      hav[j] /= pfound;
-                      hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]*pfound
-                                        /(pfound-1)));
-                      error_array[j]=sqrt(error_array[j]/pfound)/hrms[j];
-                      sumerror += error_array[j];
-                    }
-             
-                
-              // old fprintf(stdout,"%e %e ",epsilon*interval,
-              //                          sumerror/(double)dim);
-                  output(count, 0) = epsilon*interval;
-                  output(count, 1) = sumerror/(double)dim;
-                  for (octave_idx_type j=0;j<dim;j++)
-                    {
-                    //old fprintf(stdout,"%e ",error_array[j]);
-                      output(count, 2 + j) = error_array[j];
-                    }
-              // old  fprintf(stdout,"%e %e\n",(double)pfound
-              //          /(clength-(embed-1)*delay),
-              //          avfound/pfound);
-                  output(count, 2 + dim) = (double)pfound
-                                           /(clength-(embed-1)*delay);
-                  output(count, 2 + dim + 1) = avfound/pfound;
+                //old fprintf(stdout,"%e ",error_array[j]);
+                  output(count, 2 + j) = error_array[j];
+                }
+          // old  fprintf(stdout,"%e %e\n",(double)pfound
+          //          /(clength-(embed-1)*delay),
+          //          avfound/pfound);
+              output(count, 2 + dim) = (double)pfound
+                                       /(clength-(embed-1)*delay);
+              output(count, 2 + dim + 1) = avfound/pfound;
 
-                  count += 1;
-                }
+              count += 1;
             }
-          // Resize output to fit actual results instead of
-          // an educated guess
-          // if count == 0 then the output will be an 0x4+dim matrix
-          output.resize (count, dim + 4);
+        }
+      // Resize output to fit actual results instead of
+      // an educated guess
+      // if count == 0 then the output will be an 0x4+dim matrix
+      output.resize (count, dim + 4);
 
-          retval(0) = output;
-        }
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lfo_run__.cc
--- a/src/__lfo_run__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lfo_run__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -226,13 +226,6 @@
       Matrix solved_vec      = mat.solve (vec);
       double *solved_vec_arr = solved_vec.fortran_vec ();
 
-    // If errors were raised (a singular matrix was encountered), 
-    // there is no sense in countinuing
-    if (error_state)
-      {
-        return ;
-      }
-
       newcast[i]=foreav[i];
       for (octave_idx_type j=0;j<dim;j++) {
         for (octave_idx_type j1=0;j1<embed;j1++) {
@@ -330,87 +323,76 @@
         for (octave_idx_type i=0;i<hdim;i++)
           cast[i][j]=series[j][LENGTH-hdim+i];
 
-      if ( ! error_state)
-        {
-
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
 
-          // Calculate the maximum epsilon that makes sense
-          // On the basis of 'i' and 'j' from put_in_boxes ()
-          NDArray input_max = input.max ();
-          double maximum_epsilon = (input_max(0) > input_max(dim-1))
-                                   ? input_max(0) : input_max(dim-1);
-          maximum_epsilon *= EPSF;
+      // Calculate the maximum epsilon that makes sense
+      // On the basis of 'i' and 'j' from put_in_boxes ()
+      NDArray input_max = input.max ();
+      double maximum_epsilon = (input_max(0) > input_max(dim-1))
+                               ? input_max(0) : input_max(dim-1);
+      maximum_epsilon *= EPSF;
 
-          // Create output
-          Matrix output (FLENGTH, dim);
-          for (octave_idx_type i=0;i<FLENGTH;i++)
+      // Create output
+      Matrix output (FLENGTH, dim);
+      for (octave_idx_type i=0;i<FLENGTH;i++)
+        {
+          bool done=0;
+          double epsilon=EPS0/EPSF;
+          while (!done) 
             {
-              bool done=0;
-              double epsilon=EPS0/EPSF;
-              while (!done) 
+              // If epsilon became too large 
+              // there is no sense in continuing
+              if (epsilon > maximum_epsilon)
                 {
-                  // If epsilon became too large 
-                  // there is no sense in continuing
-                  if (epsilon > maximum_epsilon)
+                  error_with_id ("Octave:tisean", "The neighbourhood size"
+                                 " became too large during search,"
+                                 " no sense continuing");
+                  return retval;
+                }
+
+              epsilon*=EPSF;
+              put_in_boxes(series, LENGTH, list, box, epsilon, dim, embed,
+                           DELAY);
+              octave_idx_type actfound;
+              actfound=hfind_neighbors (series, cast, found, list, box,
+                                        epsilon, dim, embed, DELAY);
+              if (actfound >= MINN)
+                {
+                  if (!do_zeroth)
+                    make_fit(series, cast, found, dim, embed, DELAY,
+                             actfound, newcast);
+                  else
+                    make_zeroth(series, found, dim, actfound,newcast);
+
+                  for (octave_idx_type j=0;j<dim;j++)
                     {
-                      error_with_id ("Octave:tisean", "The neighbourhood size"
-                                     " became too large during search,"
-                                     " no sense continuing");
-                      return retval;
+                  // old printf("%e ",newcast[j]*interval[j]+min_array[j]);
+                      output(i,j) = newcast[j]*interval[j]+min_array[j];
                     }
 
-                  epsilon*=EPSF;
-                  put_in_boxes(series, LENGTH, list, box, epsilon, dim, embed,
-                               DELAY);
-                  octave_idx_type actfound;
-                  actfound=hfind_neighbors (series, cast, found, list, box,
-                                            epsilon, dim, embed, DELAY);
-                  if (actfound >= MINN)
+                  done=1;
+                  for (octave_idx_type j=0;j<dim;j++) 
                     {
-                      if (!do_zeroth)
-                        make_fit(series, cast, found, dim, embed, DELAY,
-                                 actfound, newcast);
-                      else
-                        make_zeroth(series, found, dim, actfound,newcast);
-
-                      // Checking if the fit was correct
-                      // If any errors were raised: end function
-                      if (error_state)
+                      // If this occurs there is no sense to continue
+                      if ((newcast[j] > 2.0) || (newcast[j] < -1.0)) 
                         {
+                          error_with_id("Octave:tisean","forecast failed, "
+                                        "escaping data region");
                           return retval;
                         }
-
-                      for (octave_idx_type j=0;j<dim;j++)
-                        {
-                      // old printf("%e ",newcast[j]*interval[j]+min_array[j]);
-                          output(i,j) = newcast[j]*interval[j]+min_array[j];
-                        }
-
-                      done=1;
-                      for (octave_idx_type j=0;j<dim;j++) 
-                        {
-                          // If this occurs there is no sense to continue
-                          if ((newcast[j] > 2.0) || (newcast[j] < -1.0)) 
-                            {
-                              error_with_id("Octave:tisean","forecast failed, "
-                                            "escaping data region");
-                              return retval;
-                            }
-                        }
-                      double *swap=cast[0];
-                      for (octave_idx_type j=0;j<hdim-1;j++)
-                        cast[j]=cast[j+1];
-                      cast[hdim-1]=swap;
-                      for (octave_idx_type j=0;j<dim;j++)
-                        cast[hdim-1][j]=newcast[j];
                     }
+                  double *swap=cast[0];
+                  for (octave_idx_type j=0;j<hdim-1;j++)
+                    cast[j]=cast[j+1];
+                  cast[hdim-1]=swap;
+                  for (octave_idx_type j=0;j<dim;j++)
+                    cast[hdim-1][j]=newcast[j];
                 }
             }
-          retval(0) = output;
         }
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lfo_test__.cc
--- a/src/__lfo_test__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lfo_test__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -212,13 +212,6 @@
       Matrix solved_vec      = mat.solve (vec);
       double *solved_vec_arr = solved_vec.fortran_vec ();
 
-      // If errors were raised (a singular matrix was encountered), 
-      // there is no sense in countinueing
-      if (error_state)
-      {
-        return ;
-      }
-
       newcast[i]=foreav[i];
       for (octave_idx_type j=0;j<DIM;j++) {
         octave_idx_type hcj=indexes[0][j];
@@ -332,90 +325,81 @@
       for (octave_idx_type i=0;i<COMP;i++)
         error_array[i]=0.0;
 
-      if ( ! error_state)
-        {
-
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
 
-          // Calculate the maximum epsilon that makes sense
-          // On the basis of 'i' and 'j' from put_in_boxes ()
-          NDArray input_max = input.max ();
-          double maximum_epsilon = (input_max(0) > input_max(COMP-1))
-                                   ? input_max(0) : input_max(COMP-1);
-          maximum_epsilon *= EPSF;
+      // Calculate the maximum epsilon that makes sense
+      // On the basis of 'i' and 'j' from put_in_boxes ()
+      NDArray input_max = input.max ();
+      double maximum_epsilon = (input_max(0) > input_max(COMP-1))
+                               ? input_max(0) : input_max(COMP-1);
+      maximum_epsilon *= EPSF;
 
-          // Calculate output
-          bool alldone=0;
-          while (!alldone) {
-            alldone=1;
+      // Calculate output
+      bool alldone=0;
+      while (!alldone) {
+        alldone=1;
 
-            // If epsilon became too large there is no sense in continuing
-            if (epsilon > maximum_epsilon)
-              {
-                error_with_id ("Octave:tisean", "The neighbourhood size became"
-                               " too large during search, no sense"
-                               " continuing");
-                return retval;
-              }
-            epsilon*=EPSF;
-            put_in_boxes(series, LENGTH, STEP, COMP, hdim, epsilon, list, box);
-            for (octave_idx_type i=(EMBED-1)*DELAY;i<clength;i++)
-              if (!done[i]) 
-                {
+        // If epsilon became too large there is no sense in continuing
+        if (epsilon > maximum_epsilon)
+          {
+            error_with_id ("Octave:tisean", "The neighbourhood size became"
+                           " too large during search, no sense"
+                           " continuing");
+            return retval;
+          }
+        epsilon*=EPSF;
+        put_in_boxes(series, LENGTH, STEP, COMP, hdim, epsilon, list, box);
+        for (octave_idx_type i=(EMBED-1)*DELAY;i<clength;i++)
+          if (!done[i]) 
+            {
 
-                  octave_idx_type actfound;
-                  actfound=hfind_neighbors(series, indexes, COMP, hdim, DIM,
-                                           epsilon, hfound, list, box, i);
-                  actfound=exclude_interval(actfound,i-causal+1,
-                                      i+causal+(EMBED-1)*DELAY-1,hfound,found);
-                  if (actfound > MINN) 
-                    {
-                      make_fit(series, indexes, found, STEP, DIM, COMP,
-                               actfound,i,newcast);
-                      // Checking if the fit was correct
-                      // If any errors were raised: end function
-                      if (error_state)
-                        {
-                          return retval;
-                        }
-                      for (octave_idx_type j=0;j<COMP;j++)
-                        error_array[j] += sqr(newcast[j]-series[j][i+STEP]);
+              octave_idx_type actfound;
+              actfound=hfind_neighbors(series, indexes, COMP, hdim, DIM,
+                                       epsilon, hfound, list, box, i);
+              actfound=exclude_interval(actfound,i-causal+1,
+                                  i+causal+(EMBED-1)*DELAY-1,hfound,found);
+              if (actfound > MINN) 
+                {
+                  make_fit(series, indexes, found, STEP, DIM, COMP,
+                           actfound,i,newcast);
+                  // Checking if the fit was correct
+                  for (octave_idx_type j=0;j<COMP;j++)
+                    error_array[j] += sqr(newcast[j]-series[j][i+STEP]);
 
-                      for (octave_idx_type j=0;j<COMP;j++)
-                        individual[j][i]=(newcast[j]-series[j][i+STEP])
-                                          *interval[j];
-                      done[i]=1;
-                    }
-                  alldone &= done[i];
+                  for (octave_idx_type j=0;j<COMP;j++)
+                    individual[j][i]=(newcast[j]-series[j][i+STEP])
+                                      *interval[j];
+                  done[i]=1;
                 }
-          }
+              alldone &= done[i];
+            }
+      }
 
-          double norm=((double)clength-(double)((EMBED-1)*DELAY));
+      double norm=((double)clength-(double)((EMBED-1)*DELAY));
 
-          // Create relative forecast error output
-          Matrix rel (COMP, 1);
-          for (octave_idx_type i=0;i<COMP;i++) 
-            {
-             // old fprintf(stdout,"# %e\n",sqrt(error_array[i]/norm)/rms[i]);
-              rel(i,0) = sqrt(error_array[i]/norm)/rms[i];
-            }
+      // Create relative forecast error output
+      Matrix rel (COMP, 1);
+      for (octave_idx_type i=0;i<COMP;i++) 
+        {
+         // old fprintf(stdout,"# %e\n",sqrt(error_array[i]/norm)/rms[i]);
+          rel(i,0) = sqrt(error_array[i]/norm)/rms[i];
+        }
 
-          // Create individual forecast error output
-          Matrix ind (clength - (EMBED-1)*DELAY,COMP);
-          for (octave_idx_type i=(EMBED-1)*DELAY;i<clength;i++)
+      // Create individual forecast error output
+      Matrix ind (clength - (EMBED-1)*DELAY,COMP);
+      for (octave_idx_type i=(EMBED-1)*DELAY;i<clength;i++)
+        {
+          for (octave_idx_type j=0;j<COMP;j++)
             {
-              for (octave_idx_type j=0;j<COMP;j++)
-                {
-                 // old fprintf(stdout,"%e ",individual[j][i]);
-                  ind(i-(EMBED-1)*DELAY, j) = individual[j][i];
-                }
+             // old fprintf(stdout,"%e ",individual[j][i]);
+              ind(i-(EMBED-1)*DELAY, j) = individual[j][i];
             }
+        }
 
-          retval(0) = rel;
-          retval(1) = ind;
-        }
+      retval(0) = rel;
+      retval(1) = ind;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lyap_k__.cc
--- a/src/__lyap_k__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lyap_k__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -211,70 +211,67 @@
         eps_fak=pow(epsmax/epsmin,1.0/(double)(epscount-1));
 
 
-      if (! error_state)
+      // Calculate exponents
+      dim_vector dv (epscount ,((int)maxdim-(int)mindim + 1));
+      string_vector keys;
+      keys.append (std::string("eps"));
+      keys.append (std::string("dim"));
+      keys.append (std::string("exp"));
+      octave_map output (dv,keys);
+
+      for (octave_idx_type l=0;l<epscount;l++)
         {
-          // Calculate exponents
-          dim_vector dv (epscount ,((int)maxdim-(int)mindim + 1));
-          string_vector keys;
-          keys.append (std::string("eps"));
-          keys.append (std::string("dim"));
-          keys.append (std::string("exp"));
-          octave_map output (dv,keys);
-
-          for (octave_idx_type l=0;l<epscount;l++)
+          double epsilon=epsmin*pow(eps_fak,(double)l);
+          for (octave_idx_type i=0;i<maxdim-1;i++)
+            for (octave_idx_type j=0;j<=maxiter;j++) {
+               count[i][j]=0;
+               lyap[i][j]=0.0;
+            }
+          put_in_boxes(series, liste, box, length, maxdim, delay, maxiter, 
+                       epsilon);
+          for (octave_idx_type i=0;i<reference;i++)
             {
-              double epsilon=epsmin*pow(eps_fak,(double)l);
-              for (octave_idx_type i=0;i<maxdim-1;i++)
-                for (octave_idx_type j=0;j<=maxiter;j++) {
-                   count[i][j]=0;
-                   lyap[i][j]=0.0;
-                }
-              put_in_boxes(series, liste, box, length, maxdim, delay, maxiter, 
-                           epsilon);
-              for (octave_idx_type i=0;i<reference;i++)
-                {
-                  lfind_neighbors(series, lfound, found, liste, box, window,
-                                  maxdim, delay, i,epsilon);
-                  iterate_points(series, lyap, count, lfound, found,
-                                 maxdim, mindim, delay, maxiter, i);
-                }
+              lfind_neighbors(series, lfound, found, liste, box, window,
+                              maxdim, delay, i,epsilon);
+              iterate_points(series, lyap, count, lfound, found,
+                             maxdim, mindim, delay, maxiter, i);
+            }
 
-              if (verbose)
-                printf("epsilon= %e\n",epsilon*max_val);
-              // Assign output
-              for (octave_idx_type i=mindim-2;i<maxdim-1;i++)
-                {
+          if (verbose)
+            printf("epsilon= %e\n",epsilon*max_val);
+          // Assign output
+          for (octave_idx_type i=mindim-2;i<maxdim-1;i++)
+            {
 
-                  // old fprintf(fout,"#epsilon= %e  dim= %d\n", 
-                  //             epsilon*max_val,i+2);
-                  octave_scalar_map tmp (keys);
-                  tmp.setfield ("eps", epsilon*max_val);
-                  tmp.setfield ("dim", i+2);
+              // old fprintf(fout,"#epsilon= %e  dim= %d\n", 
+              //             epsilon*max_val,i+2);
+              octave_scalar_map tmp (keys);
+              tmp.setfield ("eps", epsilon*max_val);
+              tmp.setfield ("dim", i+2);
 
-                  // Create matrix for the exponent data
-                  Matrix lyap_exp (maxiter + 1,3);
-                  octave_idx_type counter = 0;
-                  for (octave_idx_type j=0;j<=maxiter;j++)
-                    if (count[i][j])
-                      {
-                        // old fprintf(fout,"%d %e %ld\n",j,
-                        //             lyap[i][j]/count[i][j],count[i][j]);
+              // Create matrix for the exponent data
+              Matrix lyap_exp (maxiter + 1,3);
+              octave_idx_type counter = 0;
+              for (octave_idx_type j=0;j<=maxiter;j++)
+                if (count[i][j])
+                  {
+                    // old fprintf(fout,"%d %e %ld\n",j,
+                    //             lyap[i][j]/count[i][j],count[i][j]);
 
-                          lyap_exp(counter, 0) = j;
-                          lyap_exp(counter, 1) = lyap[i][j]/count[i][j];
-                          lyap_exp(counter, 2) = count[i][j];
-                          counter             += 1;
-                       }
-                  // Resize output to fit actual number of found exponents
-                  lyap_exp.resize(counter, 3);
-                  tmp.setfield ("exp", lyap_exp);
-                  output.assign(idx_vector(l),idx_vector(i-(int)(mindim-2)),
-                                tmp);
-                }
+                      lyap_exp(counter, 0) = j;
+                      lyap_exp(counter, 1) = lyap[i][j]/count[i][j];
+                      lyap_exp(counter, 2) = count[i][j];
+                      counter             += 1;
+                   }
+              // Resize output to fit actual number of found exponents
+              lyap_exp.resize(counter, 3);
+              tmp.setfield ("exp", lyap_exp);
+              output.assign(idx_vector(l),idx_vector(i-(int)(mindim-2)),
+                            tmp);
             }
-          // Assign output
-          retval(0) = output;
         }
+      // Assign output
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lyap_r__.cc
--- a/src/__lyap_r__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lyap_r__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -167,64 +167,60 @@
       octave_idx_type maxlength=length-delay*(dim-1)-steps-1-mindist;
       bool alldone=0;
 
-      if (! error_state)
+      // Calculate the maximum epsilon that makes sense
+      // On the basis of 'i' and 'j' from put_in_boxes ()
+      NDArray input_max = input.max ();
+      double maximum_epsilon = (input_max(0) > input_max(dim-1))
+                               ? input_max(0) : input_max(dim-1);
+      maximum_epsilon *= 1.1;
+
+      // Calculate lyapunov exponents
+      for (double eps=eps0;!alldone;eps*=1.1)
         {
 
-          // Calculate the maximum epsilon that makes sense
-          // On the basis of 'i' and 'j' from put_in_boxes ()
-          NDArray input_max = input.max ();
-          double maximum_epsilon = (input_max(0) > input_max(dim-1))
-                                   ? input_max(0) : input_max(dim-1);
-          maximum_epsilon *= 1.1;
-
-          // Calculate lyapunov exponents
-          for (double eps=eps0;!alldone;eps*=1.1)
+          // If epsilon became too large 
+          // there is no sense in continuing
+          if (eps > maximum_epsilon)
             {
-
-              // If epsilon became too large 
-              // there is no sense in continuing
-              if (eps > maximum_epsilon)
-                {
-                  error_with_id ("Octave:tisean", "The neighbourhood size"
-                                 " became too large during search,"
-                                 " no sense continuing");
-                  return retval;
-                }
-
-              put_in_boxes(series, box, list, eps, length, dim, delay, steps);
-              alldone=1;
-              for (octave_idx_type n=0;n<=maxlength;n++)
-                {
-                  if (!done[n])
-                     done[n]=make_iterate(series, box, list, found, lyap, eps,
-                                          length, dim, delay, steps, mindist,
-                                          n);
-                  alldone &= done[n];
-                }
-              if (verbose)
-                printf("epsilon: %e already found: %ld\n",eps*max_val,
-                       found[0]);
+              error_with_id ("Octave:tisean", "The neighbourhood size"
+                             " became too large during search,"
+                             " no sense continuing");
+              return retval;
             }
 
-          // Create output
-          Matrix output (steps+1, 2);
-          octave_idx_type count = 0;
-          for (octave_idx_type i=0;i<=steps;i++)
-            if (found[i])
-              {
-                // old fprintf(file,"%d %e\n",i,lyap[i]/found[i]/2.0);
-                output(count,0) = i;
-                output(count,1) = lyap[i]/found[i]/2.0;
-                count          += 1;
-              }
+          put_in_boxes(series, box, list, eps, length, dim, delay, steps);
+          alldone=1;
+          for (octave_idx_type n=0;n<=maxlength;n++)
+            {
+              if (!done[n])
+                 done[n]=make_iterate(series, box, list, found, lyap, eps,
+                                      length, dim, delay, steps, mindist,
+                                      n);
+              alldone &= done[n];
+            }
+          if (verbose)
+            printf("epsilon: %e already found: %ld\n",eps*max_val,
+                   found[0]);
+        }
 
-          // Resize output to match number of found points
-          if (count < output.numel ())
-            output.resize (count, 2);
+      // Create output
+      Matrix output (steps+1, 2);
+      octave_idx_type count = 0;
+      for (octave_idx_type i=0;i<=steps;i++)
+        if (found[i])
+          {
+            // old fprintf(file,"%d %e\n",i,lyap[i]/found[i]/2.0);
+            output(count,0) = i;
+            output(count,1) = lyap[i]/found[i]/2.0;
+            count          += 1;
+          }
 
-          // Assign output
-          retval(0) = output;
-        }
+      // Resize output to match number of found points
+      if (count < output.numel ())
+        output.resize (count, 2);
+
+      // Assign output
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lyap_spec__.cc
--- a/src/__lyap_spec__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lyap_spec__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -220,13 +220,6 @@
       Matrix solved      = mat.solve (vec);
       double *solved_arr = solved.fortran_vec ();
 
-    // If errors were raised (a singular matrix was encountered), 
-    // there is no sense in countinuing
-    if (error_state)
-      {
-        return ;
-      }
-
       double new_vec = solved_arr[0];
       for (octave_idx_type i=1;i<=alldim;i++)
         dynamics[d][i-1] = solved_arr[i];
@@ -419,149 +412,138 @@
         }
      // end old indexes = make_multi_index();
 
-      if (!error_state)
-        {
-
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
 
-          // Prepare data for running first time
-          if (count == 0)
-            {
-              averr_matrix.fill (0.0);
+      // Prepare data for running first time
+      if (count == 0)
+        {
+          averr_matrix.fill (0.0);
 
-              if (epsset)
-                epsmin /= maxinterval;
+          if (epsset)
+            epsmin /= maxinterval;
 
-              // old rnd_init(0x098342L);
-              TISEAN_rand generator (0x098342L);
-              for (octave_idx_type i=0;i<10000;i++)
-                generator.rnd_long();
-              for (octave_idx_type i=0;i<alldim;i++) {
-                factor[i]=0.0;
-                for (octave_idx_type j=0;j<alldim;j++)
-                  delta[i][j] = (double)generator.rnd_long()
-                                / (double)std::numeric_limits<octave_idx_type>
-                                          ::max ();
-              }
-              gram_schmidt(alldim, delta,lfactor);
+          // old rnd_init(0x098342L);
+          TISEAN_rand generator (0x098342L);
+          for (octave_idx_type i=0;i<10000;i++)
+            generator.rnd_long();
+          for (octave_idx_type i=0;i<alldim;i++) {
+            factor[i]=0.0;
+            for (octave_idx_type j=0;j<alldim;j++)
+              delta[i][j] = (double)generator.rnd_long()
+                            / (double)std::numeric_limits<octave_idx_type>
+                                      ::max ();
+          }
+          gram_schmidt(alldim, delta,lfactor);
 
-              avneig = 0.0;
-              aveps  = 0.0;
-            }
+          avneig = 0.0;
+          aveps  = 0.0;
+        }
 
-          // Create output
-          Matrix lyap_exp (1, 1 + alldim);
-          time_t lasttime;
-          time(&lasttime);
-          bool pause_calc = false;
-          for (octave_idx_type i = count + (EMBED-1) * DELAY;
-               i < ITERATIONS && !pause_calc; i++)
+      // Create output
+      Matrix lyap_exp (1, 1 + alldim);
+      time_t lasttime;
+      time(&lasttime);
+      bool pause_calc = false;
+      for (octave_idx_type i = count + (EMBED-1) * DELAY;
+           i < ITERATIONS && !pause_calc; i++)
+        {
+          count++;
+          make_dynamics(series, box, indexes, epsmin, epsset, EPSSTEP,
+                        EMBED, MINNEIGHBORS, LENGTH, DIMENSION, count,
+                        avneig, aveps, dynamics, averr, i);
+          make_iteration(DIMENSION, alldim, dynamics, delta);
+          gram_schmidt(alldim, delta,lfactor);
+          for (octave_idx_type j=0;j<alldim;j++) {
+            factor[j] += log(lfactor[j])/(double)DELAY;
+          }
+
+          if (((time(NULL)-lasttime) > OUT) || (i == (ITERATIONS-1))
+              || (count % it_pause == 0) )
             {
-              count++;
-              make_dynamics(series, box, indexes, epsmin, epsset, EPSSTEP,
-                            EMBED, MINNEIGHBORS, LENGTH, DIMENSION, count,
-                            avneig, aveps, dynamics, averr, i);
-              // If there was an error
-              // (matrix singularity or not enough neighbors)
-              // No sense continuing
-              if (error_state)
+              time(&lasttime);
+
+              // Create spectrum output
+              // old fprintf(stdout,"%ld ",count);
+              lyap_exp(0,0) = count;
+              for (octave_idx_type j=0;j<alldim;j++)
                 {
-                  return retval;
+                  // old fprintf(stdout,"%e ",factor[j]/count);
+                  lyap_exp(0, 1+j) = factor[j]/count;
                 }
-              make_iteration(DIMENSION, alldim, dynamics, delta);
-              gram_schmidt(alldim, delta,lfactor);
-              for (octave_idx_type j=0;j<alldim;j++) {
-                factor[j] += log(lfactor[j])/(double)DELAY;
-              }
 
-              if (((time(NULL)-lasttime) > OUT) || (i == (ITERATIONS-1))
-                  || (count % it_pause == 0) )
-                {
-                  time(&lasttime);
-
-                  // Create spectrum output
-                  // old fprintf(stdout,"%ld ",count);
-                  lyap_exp(0,0) = count;
-                  for (octave_idx_type j=0;j<alldim;j++)
-                    {
-                      // old fprintf(stdout,"%e ",factor[j]/count);
-                      lyap_exp(0, 1+j) = factor[j]/count;
-                    }
+              pause_calc = true;
+            }
+        }
 
-                  pause_calc = true;
-                }
-            }
-
-          // Create pause output
-          if (count < (ITERATIONS - (EMBED-1)*DELAY))
-            {
-              octave_scalar_map pause_vars;
-              pause_vars.setfield ("averr", averr_matrix);
-              pause_vars.setfield ("delta", delta_matrix);
-              pause_vars.setfield ("count", count);
-              pause_vars.setfield ("avneig", avneig);
-              pause_vars.setfield ("aveps", aveps);
-              pause_vars.setfield ("epsmin", epsmin);
-              retval(0) = lyap_exp; 
-              retval(1) = pause_vars;
-            }
+      // Create pause output
+      if (count < (ITERATIONS - (EMBED-1)*DELAY))
+        {
+          octave_scalar_map pause_vars;
+          pause_vars.setfield ("averr", averr_matrix);
+          pause_vars.setfield ("delta", delta_matrix);
+          pause_vars.setfield ("count", count);
+          pause_vars.setfield ("avneig", avneig);
+          pause_vars.setfield ("aveps", aveps);
+          pause_vars.setfield ("epsmin", epsmin);
+          retval(0) = lyap_exp; 
+          retval(1) = pause_vars;
+        }
 
-          // Create final output
-          if (count == (ITERATIONS - (EMBED-1)*DELAY))
-            {
-              double dim=0.0;
-              octave_idx_type i;
-              for (i=0;i<alldim;i++) {
-                dim += factor[i];
-                if (dim < 0.0)
-                  break;
-              }
-              if (i < alldim)
-                dim=i+(dim-factor[i])/fabs(factor[i]);
-              else
-                dim=alldim;
+      // Create final output
+      if (count == (ITERATIONS - (EMBED-1)*DELAY))
+        {
+          double dim=0.0;
+          octave_idx_type i;
+          for (i=0;i<alldim;i++) {
+            dim += factor[i];
+            if (dim < 0.0)
+              break;
+          }
+          if (i < alldim)
+            dim=i+(dim-factor[i])/fabs(factor[i]);
+          else
+            dim=alldim;
 
-              // Create output pars
-              octave_scalar_map pars;
+          // Create output pars
+          octave_scalar_map pars;
 
-              // Create rel_err
-              Matrix rel_err (1, DIMENSION);
-              // old fprintf(stdout,"#Average relative forecast errors:= ");
-              for (octave_idx_type i=0;i<DIMENSION;i++)
-                {
-                  // old fprintf(stdout,"%e ",sqrt(averr[i]/count)/var[i]);
-                  rel_err(0,i) = sqrt(averr[i]/count)/var[i];
-                }
-              pars.setfield ("rel_err", rel_err);
+          // Create rel_err
+          Matrix rel_err (1, DIMENSION);
+          // old fprintf(stdout,"#Average relative forecast errors:= ");
+          for (octave_idx_type i=0;i<DIMENSION;i++)
+            {
+              // old fprintf(stdout,"%e ",sqrt(averr[i]/count)/var[i]);
+              rel_err(0,i) = sqrt(averr[i]/count)/var[i];
+            }
+          pars.setfield ("rel_err", rel_err);
 
-              // Create abs_err
-              Matrix abs_err (1, DIMENSION);
-              // old fprintf(stdout,"#Average absolute forecast errors:= ");
-              for (octave_idx_type i=0;i<DIMENSION;i++)
-                {
-                  // old fprintf(stdout,"%e ",sqrt(averr[i]/count)
-                  //                          *interval[i]);
-                  abs_err(0,i) = sqrt(averr[i]/count)*interval[i];
-                }
-              pars.setfield ("abs_err", abs_err);
+          // Create abs_err
+          Matrix abs_err (1, DIMENSION);
+          // old fprintf(stdout,"#Average absolute forecast errors:= ");
+          for (octave_idx_type i=0;i<DIMENSION;i++)
+            {
+              // old fprintf(stdout,"%e ",sqrt(averr[i]/count)
+              //                          *interval[i]);
+              abs_err(0,i) = sqrt(averr[i]/count)*interval[i];
+            }
+          pars.setfield ("abs_err", abs_err);
 
-              // old fprintf(stdout,"#Average Neighborhood Size= %e\n",
-              //             aveps*maxinterval/count);
-              pars.setfield ("nsize", aveps*maxinterval/count);
-
-              // old fprintf(stdout,"#Average num. of neighbors= %e\n",
-              //             avneig/count);
-              pars.setfield ("nno", avneig/count);
+          // old fprintf(stdout,"#Average Neighborhood Size= %e\n",
+          //             aveps*maxinterval/count);
+          pars.setfield ("nsize", aveps*maxinterval/count);
 
-              // old fprintf(stdout,"#estimated KY-Dimension= %f\n",dim);
-              pars.setfield ("ky_dim", dim);
+          // old fprintf(stdout,"#Average num. of neighbors= %e\n",
+          //             avneig/count);
+          pars.setfield ("nno", avneig/count);
 
-              // Assign output
-              retval(0) = lyap_exp;
-              retval(1) = pars;
-            }
+          // old fprintf(stdout,"#estimated KY-Dimension= %f\n",dim);
+          pars.setfield ("ky_dim", dim);
+
+          // Assign output
+          retval(0) = lyap_exp;
+          retval(1) = pars;
         }
     }
   return retval;
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lzo_gm__.cc
--- a/src/__lzo_gm__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lzo_gm__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -110,92 +110,89 @@
 
       MArray<octave_idx_type> box (dim_vector(NMAX,NMAX));
 
-      if ( ! error_state)
-        {
-          // Estimate maximum possible output size
-          octave_idx_type output_rows = (octave_idx_type)
-                                        ((log(EPS1) - log(EPS0)) / log (EPSF));
-          output_rows += 2;
+      // Estimate maximum possible output size
+      octave_idx_type output_rows = (octave_idx_type)
+                                    ((log(EPS1) - log(EPS0)) / log (EPSF));
+      output_rows += 2;
 
-          // Create output
-          Matrix output (output_rows,dim+4);
-          octave_idx_type count = 0;
+      // Create output
+      Matrix output (output_rows,dim+4);
+      octave_idx_type count = 0;
 
-          for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) 
+      for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) 
+        {
+          long pfound=0;
+          for (octave_idx_type i=0;i<dim;i++)
+            error_array[i]=hrms[i]=hav[i]=0.0;
+          double avfound=0.0;
+
+          make_multi_box(input,box,list,LENGTH-STEP,NMAX,dim,
+                         embed,delay,epsilon);
+          for (octave_idx_type i=(embed-1)*delay;i<clength;i++) 
             {
-              long pfound=0;
-              for (octave_idx_type i=0;i<dim;i++)
-                error_array[i]=hrms[i]=hav[i]=0.0;
-              double avfound=0.0;
-
-              make_multi_box(input,box,list,LENGTH-STEP,NMAX,dim,
-                             embed,delay,epsilon);
-              for (octave_idx_type i=(embed-1)*delay;i<clength;i++) 
+              for (octave_idx_type j=0;j<dim;j++)
                 {
-                  for (octave_idx_type j=0;j<dim;j++)
-                    {
-                      // old hser[j]=series[j]+i;
-                      hser[j] = input.fortran_vec() + j * LENGTH + i;
-                    }
-                  octave_idx_type actfound;
-                  actfound = find_multi_neighbors(input,box,list,hser,
-                                                  NMAX,dim,embed,delay,
-                                                  epsilon,hfound);
-                  actfound = exclude_interval(actfound,i-causal+1,
-                                             i+causal+(embed-1)*delay-1,
-                                             hfound,found);
-                  if (actfound > 2*(dim*embed+1)) 
-                    {
-                      make_fit (input, dim, i, actfound, STEP, found,
-                               error_array);
-                      pfound++;
-                      avfound += (double)(actfound-1);
-                      for (octave_idx_type j=0;j<dim;j++) {
-                        hrms[j] += input(i+STEP,j) * input(i+STEP,j);
-                        hav[j] += input(i+STEP,j);
-                      }
-                    }
+                  // old hser[j]=series[j]+i;
+                  hser[j] = input.fortran_vec() + j * LENGTH + i;
                 }
-              if (pfound > 1) 
+              octave_idx_type actfound;
+              actfound = find_multi_neighbors(input,box,list,hser,
+                                              NMAX,dim,embed,delay,
+                                              epsilon,hfound);
+              actfound = exclude_interval(actfound,i-causal+1,
+                                         i+causal+(embed-1)*delay-1,
+                                         hfound,found);
+              if (actfound > 2*(dim*embed+1)) 
                 {
-                  double sumerror=0.0;
-                  for (octave_idx_type j=0;j<dim;j++) 
-                    {
-                      hav[j] /= pfound;
-                      hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]
-                                   * pfound/(pfound-1)));
-                      error_array[j]=sqrt(error_array[j]/pfound)/hrms[j];
-                      sumerror += error_array[j];
-                    }
-
-                  // Write output
-        // old fprintf(stdout,"%e %e ",epsilon*interval,sumerror/(double)dim);
-                  output(count, 0) = epsilon * interval;
-                  output(count, 1) = sumerror / (double) dim;
-                for (octave_idx_type j=0;j<dim;j++)
-                  {
-                    // old fprintf(stdout,"%e ",error_array[j]);
-                    output(count, 2 + j) = error_array[j];
+                  make_fit (input, dim, i, actfound, STEP, found,
+                           error_array);
+                  pfound++;
+                  avfound += (double)(actfound-1);
+                  for (octave_idx_type j=0;j<dim;j++) {
+                    hrms[j] += input(i+STEP,j) * input(i+STEP,j);
+                    hav[j] += input(i+STEP,j);
                   }
-            
-  //    old fprintf(stdout,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
-  //                                avfound/pfound);
-                output(count,2 + dim) = (double)pfound / 
-                                        (clength-(embed-1)*delay);
-                output(count,2 + dim + 1) = avfound/pfound;
-
-                count += 1;
                 }
             }
+          if (pfound > 1) 
+            {
+              double sumerror=0.0;
+              for (octave_idx_type j=0;j<dim;j++) 
+                {
+                  hav[j] /= pfound;
+                  hrms[j]=sqrt(fabs(hrms[j]/(pfound-1)-hav[j]*hav[j]
+                               * pfound/(pfound-1)));
+                  error_array[j]=sqrt(error_array[j]/pfound)/hrms[j];
+                  sumerror += error_array[j];
+                }
 
-          // Resize output to fit actual results instead of
-          // an educated guess
-          // if count == 0 then the output will be an 0x4+dim matrix
-          output.resize (count, dim + 4);
+              // Write output
+    // old fprintf(stdout,"%e %e ",epsilon*interval,sumerror/(double)dim);
+              output(count, 0) = epsilon * interval;
+              output(count, 1) = sumerror / (double) dim;
+            for (octave_idx_type j=0;j<dim;j++)
+              {
+                // old fprintf(stdout,"%e ",error_array[j]);
+                output(count, 2 + j) = error_array[j];
+              }
+        
+// old fprintf(stdout,"%e %e\n",(double)pfound/(clength-(embed-1)*delay),
+//                              avfound/pfound);
+            output(count,2 + dim) = (double)pfound / 
+                                    (clength-(embed-1)*delay);
+            output(count,2 + dim + 1) = avfound/pfound;
 
-          retval(0) = output;
+            count += 1;
+            }
+        }
 
-        }
+      // Resize output to fit actual results instead of
+      // an educated guess
+      // if count == 0 then the output will be an 0x4+dim matrix
+      output.resize (count, dim + 4);
+
+      retval(0) = output;
+
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lzo_run__.cc
--- a/src/__lzo_run__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lzo_run__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -209,90 +209,87 @@
 
       MArray<octave_idx_type> box (dim_vector(NMAX,NMAX));
 
-      if ( ! error_state)
-        {
-          for (octave_idx_type j=0;j<dim;j++)
-            for (octave_idx_type i=0;i<hdim;i++)
-              cast[i][j]=input(LENGTH-hdim+i,j);
+      for (octave_idx_type j=0;j<dim;j++)
+        for (octave_idx_type i=0;i<hdim;i++)
+          cast[i][j]=input(LENGTH-hdim+i,j);
 
-        // old  indexes=make_multi_index(dim,embed,DELAY);
+    // old  indexes=make_multi_index(dim,embed,DELAY);
 
-          octave_idx_type alldim=dim * embed;
+      octave_idx_type alldim=dim * embed;
 
-          MArray<octave_idx_type> indexes (dim_vector (alldim, 2));
-          for (octave_idx_type i=0;i<alldim;i++) 
-            {
-              indexes(i,0)=i%dim;
-              indexes(i,1)=(i/dim)*DELAY;
-            }
+      MArray<octave_idx_type> indexes (dim_vector (alldim, 2));
+      for (octave_idx_type i=0;i<alldim;i++) 
+        {
+          indexes(i,0)=i%dim;
+          indexes(i,1)=(i/dim)*DELAY;
+        }
 
-         // end old index = make_multi_index();
+     // end old index = make_multi_index();
 
-      //   old rnd_init(seed);
-          TISEAN_rand generator (seed);
+  //   old rnd_init(seed);
+      TISEAN_rand generator (seed);
 
-          double epsilon0=EPS0/EPSF;
+      double epsilon0=EPS0/EPSF;
 
-          if (setnoise) 
-            Q /= 100.0;
+      if (setnoise) 
+        Q /= 100.0;
 
-          Matrix output (FLENGTH, dim);
-          octave_idx_type row_count = 0;
-          octave_idx_type count = 1;
-          for (octave_idx_type i=0;i<FLENGTH;i++)
-            {
-              bool done=0;
-              double epsilon;
-              if (setsort)
-                epsilon= epsilon0/((double)count*EPSF);
-              else
-                epsilon=epsilon0;
-              while (!done) {
-                epsilon*=EPSF;
-                put_in_boxes(input, box, list, epsilon,(embed-1)*DELAY);
-                octave_idx_type actfound=hfind_neighbors(input, indexes, box,
-                                                         list, cast, found,
-                                                         epsilon,
-                                                         (embed-1) * DELAY);
-                if (actfound >= MINN) {
-                  if (setsort) {
-                    epsilon0 += epsilon;
-                    count++;
-                    sort(input, found, abstand, cast, actfound, (embed-1)*DELAY);
-                    actfound=MINN;
-                  }
-                  make_zeroth(input, generator, setnoise, var, Q, found,
-                              actfound, newcast);
+      Matrix output (FLENGTH, dim);
+      octave_idx_type row_count = 0;
+      octave_idx_type count = 1;
+      for (octave_idx_type i=0;i<FLENGTH;i++)
+        {
+          bool done=0;
+          double epsilon;
+          if (setsort)
+            epsilon= epsilon0/((double)count*EPSF);
+          else
+            epsilon=epsilon0;
+          while (!done) {
+            epsilon*=EPSF;
+            put_in_boxes(input, box, list, epsilon,(embed-1)*DELAY);
+            octave_idx_type actfound=hfind_neighbors(input, indexes, box,
+                                                     list, cast, found,
+                                                     epsilon,
+                                                     (embed-1) * DELAY);
+            if (actfound >= MINN) {
+              if (setsort) {
+                epsilon0 += epsilon;
+                count++;
+                sort(input, found, abstand, cast, actfound, (embed-1)*DELAY);
+                actfound=MINN;
+              }
+              make_zeroth(input, generator, setnoise, var, Q, found,
+                          actfound, newcast);
 
-                  for (octave_idx_type j=0;j<dim-1;j++)
-                  {
-                    // old printf("%e ",newcast[j]*interval[j]+min[j]);
-                    output(row_count, j) = newcast[j]*interval[j]+min[j];
-                  }
-                  // old printf("%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
-                  output(row_count, dim-1) = newcast[dim-1]*interval[dim-1]
-                                             + min[dim-1];
-                  row_count += 1;
+              for (octave_idx_type j=0;j<dim-1;j++)
+                {
+                  // old printf("%e ",newcast[j]*interval[j]+min[j]);
+                  output(row_count, j) = newcast[j]*interval[j]+min[j];
+                }
+              // old printf("%e\n",newcast[dim-1]*interval[dim-1]+min[dim-1]);
+              output(row_count, dim-1) = newcast[dim-1]*interval[dim-1]
+                                         + min[dim-1];
+              row_count += 1;
 
-                  done=1;
-                  double *swap=cast[0];
-                  for (octave_idx_type j=0;j<hdim-1;j++)
-                    cast[j]=cast[j+1];
-                  cast[hdim-1]=swap;
-                  for (octave_idx_type j=0;j<dim;j++)
-                    cast[hdim-1][j]=newcast[j];
-                }
-              }
+              done=1;
+              double *swap=cast[0];
+              for (octave_idx_type j=0;j<hdim-1;j++)
+                cast[j]=cast[j+1];
+              cast[hdim-1]=swap;
+              for (octave_idx_type j=0;j<dim;j++)
+                cast[hdim-1][j]=newcast[j];
             }
+          }
+        }
 
-            if (row_count != 0)
-            {
-              output.resize (row_count, dim);
-              retval(0) = output;
-            }
-            else 
-              retval(0) = Matrix (0,0);
-        }
+        if (row_count != 0)
+          {
+            output.resize (row_count, dim);
+            retval(0) = output;
+          }
+        else 
+          retval(0) = Matrix (0,0);
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__lzo_test__.cc
--- a/src/__lzo_test__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__lzo_test__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -150,93 +150,90 @@
       MArray<octave_idx_type> box (dim_vector(NMAX, NMAX));
 
       // Compute forecast error
-      if ( ! error_state)
-        {
-          for (octave_idx_type i=0;i<LENGTH;i++)
-            done[i]=0;
+      for (octave_idx_type i=0;i<LENGTH;i++)
+        done[i]=0;
+
+      bool alldone = false;
+      if (epsset)
+        EPS0 /= interval;
 
-          bool alldone = false;
-          if (epsset)
-            EPS0 /= interval;
+      double epsilon = EPS0 / EPSF;
 
-          double epsilon = EPS0 / EPSF;
+      if (!clengthset)
+        CLENGTH=LENGTH;
+      octave_idx_type clength = ((CLENGTH*refstep+STEP) <= LENGTH) 
+                                ? CLENGTH : (LENGTH-STEP)/refstep;
 
-          if (!clengthset)
-            CLENGTH=LENGTH;
-          octave_idx_type clength = ((CLENGTH*refstep+STEP) <= LENGTH) 
-                                    ? CLENGTH : (LENGTH-STEP)/refstep;
+      // Compute estimates
+      octave_idx_type actfound;
+      octave_idx_type hi;
+      while (!alldone) {
+        alldone=1;
+        epsilon*=EPSF;
+        make_multi_box(input,box,list,LENGTH-(long)STEP,NMAX,dim,
+                       embed,DELAY,epsilon);
+        for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++)
+          if (!done[i]) {
+          hi=i*refstep;
+          for (octave_idx_type j=0;j<dim;j++)
+            {
+//          old  hser[j]=series[j]+hi;
+              hser[j] = input.fortran_vec() + j * LENGTH + hi;
+            }
+          actfound=find_multi_neighbors(input,box,list,hser,NMAX,
+                                        dim,embed,DELAY,epsilon,hfound);
+          actfound=exclude_interval(actfound,hi-(long)causal+1,
+                              hi+causal+(embed-1)*DELAY-1,hfound,found);
+          if (actfound >= MINN)
+            {
+              if (setsort)
+                {
+                  sort(input, found, abstand, embed, DELAY, MINN,
+                       actfound, hser);
+                  actfound=MINN;
+                }
+              for (octave_idx_type j=1;j<=STEP;j++) {
+                make_fit(input,dim,hi,actfound,j,found,error_array,diffs);
+              }
+              done[i]=1;
+            }
+          alldone &= done[i];
+          }
+      }
 
-          // Compute estimates
-          octave_idx_type actfound;
-          octave_idx_type hi;
-          while (!alldone) {
-            alldone=1;
-            epsilon*=EPSF;
-            make_multi_box(input,box,list,LENGTH-(long)STEP,NMAX,dim,
-                           embed,DELAY,epsilon);
-            for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++)
-              if (!done[i]) {
+      // Create relative forecast error output
+      Matrix rel_forecast_err (STEP, dim + 1);
+      for (octave_idx_type i=0;i<STEP;i++) 
+        {
+          rel_forecast_err(i,0) = i + 1;
+          for (octave_idx_type j=0;j<dim;j++) 
+            {
+         // old  fprintf(stdout,"%e ",
+         //         sqrt(error[j][i]/(clength-(embed-1)*DELAY))/rms[j]);
+              rel_forecast_err(i,j+1) = sqrt(error_array(i,j)
+                                      /(clength-(embed-1)*DELAY))/rms[j];
+            }
+        }
+
+      // Create individual forecast error output
+      Matrix ind_forecast_err (1,1);
+      if (nargout > 1)
+        {
+          ind_forecast_err.resize(clength - (embed-1)*DELAY, dim);
+          for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++) 
+            {
               hi=i*refstep;
               for (octave_idx_type j=0;j<dim;j++)
                 {
-//              old  hser[j]=series[j]+hi;
-                  hser[j] = input.fortran_vec() + j * LENGTH + hi;
-                }
-              actfound=find_multi_neighbors(input,box,list,hser,NMAX,
-                                            dim,embed,DELAY,epsilon,hfound);
-              actfound=exclude_interval(actfound,hi-(long)causal+1,
-                                  hi+causal+(embed-1)*DELAY-1,hfound,found);
-              if (actfound >= MINN)
-                {
-                  if (setsort)
-                    {
-                      sort(input, found, abstand, embed, DELAY, MINN,
-                           actfound, hser);
-                      actfound=MINN;
-                    }
-                  for (octave_idx_type j=1;j<=STEP;j++) {
-                    make_fit(input,dim,hi,actfound,j,found,error_array,diffs);
-                  }
-                  done[i]=1;
-                }
-              alldone &= done[i];
-              }
-          }
-
-          // Create relative forecast error output
-          Matrix rel_forecast_err (STEP, dim + 1);
-          for (octave_idx_type i=0;i<STEP;i++) 
-            {
-              rel_forecast_err(i,0) = i + 1;
-              for (octave_idx_type j=0;j<dim;j++) 
-                {
-//            old  fprintf(stdout,"%e ",
-//                    sqrt(error[j][i]/(clength-(embed-1)*DELAY))/rms[j]);
-                  rel_forecast_err(i,j+1) = sqrt(error_array(i,j)
-                                          /(clength-(embed-1)*DELAY))/rms[j];
+                  // old fprintf(stdout,"%e ",diffs[j][hi]*hinter[j]);
+                  ind_forecast_err(i-(embed-1)*DELAY,j) = \
+                  diffs(hi,j)*hinter[j];
                 }
             }
+        }
 
-          // Create individual forecast error output
-          Matrix ind_forecast_err (1,1);
-          if (nargout > 1)
-            {
-              ind_forecast_err.resize(clength - (embed-1)*DELAY, dim);
-              for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++) 
-                {
-                  hi=i*refstep;
-                  for (octave_idx_type j=0;j<dim;j++)
-                    {
-//                   old fprintf(stdout,"%e ",diffs[j][hi]*hinter[j]);
-                      ind_forecast_err(i-(embed-1)*DELAY,j) = \
-                      diffs(hi,j)*hinter[j];
-                    }
-                }
-            }
-
-          retval(0) = rel_forecast_err;
-          retval(1) = ind_forecast_err;
-        }
+      retval(0) = rel_forecast_err;
+      retval(1) = ind_forecast_err;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__poincare__.cc
--- a/src/__poincare__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__poincare__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -129,19 +129,16 @@
 
 
 
-        if ( ! error_state)
-          {
-            Matrix output (out_size, embdim);
+      Matrix output (out_size, embdim);
 
-            octave_idx_type count = poincare (input.fortran_vec(), 
-                                              input.numel(), embdim,
-                                              delay, comp, where, direction,
-                                              output);
-            // Resize output to fit sections found
-            output.resize (count, embdim);
+      octave_idx_type count = poincare (input.fortran_vec(), 
+                                        input.numel(), embdim,
+                                        delay, comp, where, direction,
+                                        output);
+      // Resize output to fit sections found
+      output.resize (count, embdim);
 
-            retval(0) = output;
-          }
+      retval(0) = output;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__polynom__.cc
--- a/src/__polynom__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__polynom__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -222,88 +222,79 @@
       make_coding(coding_vec,N,N,DIM-1,1,0);
       octave_idx_type *coding = coding_vec.data();
 
-      if (! error_state)
-        {
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
- 
-          // Make the fit
-          make_fit (series, coding, results, INSAMPLE, N, DIM, DELAY,
-                    maxencode, pars);
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
+
+      // Make the fit
+      make_fit (series, coding, results, INSAMPLE, N, DIM, DELAY,
+                maxencode, pars);
+
+      // Create outputs
 
-          // If error encountered during the fit there is no sense to continue
-          if (error_state)
-            {
-              return retval;
-            }
+      // Create output that contains the number of free parameters
+      // old fprintf(file,"#number of free parameters= %d\n\n",pars);
+      octave_idx_type free_par = pars;
 
-          // Create outputs
+      // Create output that contains the norm used for the fit
+      // old fprintf(file,"#used norm for the fit= %e\n",std_dev);
+      double fit_norm = std_dev;
 
-          // Create output that contains the number of free parameters
-          // old fprintf(file,"#number of free parameters= %d\n\n",pars);
-          octave_idx_type free_par = pars;
-
-          // Create output that contains the norm used for the fit
-          // old fprintf(file,"#used norm for the fit= %e\n",std_dev);
-          double fit_norm = std_dev;
-
-          // Create coefficients output
-          Matrix coeffs (pars, DIM + 1);
-          for (octave_idx_type j=0;j<pars;j++)
+      // Create coefficients output
+      Matrix coeffs (pars, DIM + 1);
+      for (octave_idx_type j=0;j<pars;j++)
+        {
+          decode(N,opar,DIM-1,coding[j],maxencode);
+          octave_idx_type sumpar=0;
+          for (octave_idx_type k=0;k<DIM;k++)
             {
-              decode(N,opar,DIM-1,coding[j],maxencode);
-              octave_idx_type sumpar=0;
-              for (octave_idx_type k=0;k<DIM;k++)
-                {
-                  sumpar += opar[k];
-                //  old fprintf(file,"%d ",opar[k]);
-                  coeffs(j, k) = opar[k];
-                }
-          //  old fprintf(file,"%e\n",results[j]
-          //                          /pow(std_dev,(double)(sumpar-1)));
-              coeffs(j,DIM) = results[j]/pow(std_dev,(double)(sumpar-1));
-
+              sumpar += opar[k];
+            //  old fprintf(file,"%d ",opar[k]);
+              coeffs(j, k) = opar[k];
             }
+      //  old fprintf(file,"%e\n",results[j]
+      //                          /pow(std_dev,(double)(sumpar-1)));
+          coeffs(j,DIM) = results[j]/pow(std_dev,(double)(sumpar-1));
 
-          // Create sample error
-          // 1st element of sample_error is the insample error
-          // 2nd element of sample_error is the out of sample error (if exists)
-          NDArray sample_err (dim_vector (1,1));
+        }
+
+      // Create sample error
+      // 1st element of sample_error is the insample error
+      // 2nd element of sample_error is the out of sample error (if exists)
+      NDArray sample_err (dim_vector (1,1));
 
-    // old in_error = make_error((unsigned long)0,INSAMPLE)
-    //     fprintf(file,"#average insample error= %e\n",
-    //             sqrt(in_error)*std_dev);
-          sample_err(0) = sqrt(make_error(series, coding, results, N, DIM,
-                                          DELAY, maxencode, pars,
-                                          0,INSAMPLE))
+      // old in_error = make_error((unsigned long)0,INSAMPLE)
+      //     fprintf(file,"#average insample error= %e\n",
+      //             sqrt(in_error)*std_dev);
+      sample_err(0) = sqrt(make_error(series, coding, results, N, DIM,
+                                      DELAY, maxencode, pars,
+                                      0,INSAMPLE))
+                      * std_dev;
+
+      if (INSAMPLE < LENGTH) 
+        {
+        // old out_error=make_error(INSAMPLE,LENGTH);
+        //     fprintf(file,"#average out of sample error= %e\n",
+        //             sqrt(out_error)*std_dev);
+          sample_err.resize (dim_vector (2,1));
+          sample_err(1) = sqrt (make_error (series, coding, results, N,
+                                            DIM, DELAY, maxencode, pars,
+                                            INSAMPLE, LENGTH)) 
                           * std_dev;
-
-          if (INSAMPLE < LENGTH) 
-            {
-            // old out_error=make_error(INSAMPLE,LENGTH);
-            //     fprintf(file,"#average out of sample error= %e\n",
-            //             sqrt(out_error)*std_dev);
-              sample_err.resize (dim_vector (2,1));
-              sample_err(1) = sqrt (make_error (series, coding, results, N,
-                                                DIM, DELAY, maxencode, pars,
-                                                INSAMPLE, LENGTH)) 
-                              * std_dev;
-            }
+        }
 
-          // Create forecast
-          NDArray forecast (dim_vector (0,0));
-          if (CLENGTH > 0)
-            make_cast(forecast, series, coding, results, LENGTH, CLENGTH, N,
-                      DIM, DELAY, maxencode, pars, std_dev);
+      // Create forecast
+      NDArray forecast (dim_vector (0,0));
+      if (CLENGTH > 0)
+        make_cast(forecast, series, coding, results, LENGTH, CLENGTH, N,
+                  DIM, DELAY, maxencode, pars, std_dev);
 
-          // Assign outputs
-          retval(0) = free_par;
-          retval(1) = fit_norm;
-          retval(2) = coeffs;
-          retval(3) = sample_err;
-          retval(4) = forecast;
-        }
+      // Assign outputs
+      retval(0) = free_par;
+      retval(1) = fit_norm;
+      retval(2) = coeffs;
+      retval(3) = sample_err;
+      retval(4) = forecast;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__rbf__.cc
--- a/src/__rbf__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__rbf__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -152,169 +152,161 @@
         for (octave_idx_type j=0;j<DIM;j++)
           center[i][j]=series[(DIM-1)*DELAY-j*DELAY+(i*cstep)/(CENTER-1)];
 
-      if (! error_state)
-        {
+      // Promote warnings connected with singular matrixes to errors
+      set_warning_state ("Octave:nearly-singular-matrix","error");
+      set_warning_state ("Octave:singular-matrix","error");
 
-          // Promote warnings connected with singular matrixes to errors
-          set_warning_state ("Octave:nearly-singular-matrix","error");
-          set_warning_state ("Octave:singular-matrix","error");
-
-          // Calculate coefficients
-          if (setdrift)
-            drift(CENTER, DIM, center);
-          varianz=avdistance(CENTER, DIM, center);
+      // Calculate coefficients
+      if (setdrift)
+        drift(CENTER, DIM, center);
+      varianz=avdistance(CENTER, DIM, center);
 
-          // old make_fit();
-          Matrix mat (CENTER + 1, CENTER + 1);
-          OCTAVE_LOCAL_BUFFER (double *, mat_arr, CENTER + 1);
-          for (octave_idx_type i=0;i <CENTER + 1;i++)
-            mat_arr[i]=mat.fortran_vec () + (CENTER + 1) * i;
+      // old make_fit();
+      Matrix mat (CENTER + 1, CENTER + 1);
+      OCTAVE_LOCAL_BUFFER (double *, mat_arr, CENTER + 1);
+      for (octave_idx_type i=0;i <CENTER + 1;i++)
+        mat_arr[i]=mat.fortran_vec () + (CENTER + 1) * i;
 
-          for (octave_idx_type i=0;i<=CENTER;i++) {
-            coefs_arr[i]=0.0;
-            for (octave_idx_type j=0;j<=CENTER;j++)
-              mat_arr[i][j]=0.0;
-          }
+      for (octave_idx_type i=0;i<=CENTER;i++) {
+        coefs_arr[i]=0.0;
+        for (octave_idx_type j=0;j<=CENTER;j++)
+          mat_arr[i][j]=0.0;
+      }
 
-          OCTAVE_LOCAL_BUFFER (double, hcen, CENTER);
+      OCTAVE_LOCAL_BUFFER (double, hcen, CENTER);
 
-          for (octave_idx_type n=(DIM-1)*DELAY;n<INSAMPLE-STEP;n++) {
-            octave_idx_type nst=n+STEP;
-            for (octave_idx_type i=0;i<CENTER;i++)
-              hcen[i]=rbf(DELAY, DIM, varianz, &series[n],center[i]);
-            coefs_arr[0] += series[nst];
-            mat_arr[0][0] += 1.0;
-            for (octave_idx_type i=1;i<=CENTER;i++)
-              mat_arr[i][0] += hcen[i-1];
-            for (octave_idx_type i=1;i<=CENTER;i++) {
-              double h = hcen[i-1];
-              coefs_arr[i] += series[nst] * h;
-              for (octave_idx_type j=1;j<=i;j++)
-                mat_arr[i][j] += h*hcen[j-1];
-            }
-          }
-          
-          double h=(double)(INSAMPLE-STEP-(DIM-1)*DELAY);
-          for (octave_idx_type i=0;i<=CENTER;i++) {
-            coefs_arr[i] /= h;
-            for (octave_idx_type j=0;j<=i;j++) {
-              mat_arr[i][j] /= h;
-              mat_arr[j][i]=mat_arr[i][j];
-            }
-          }
+      for (octave_idx_type n=(DIM-1)*DELAY;n<INSAMPLE-STEP;n++) {
+        octave_idx_type nst=n+STEP;
+        for (octave_idx_type i=0;i<CENTER;i++)
+          hcen[i]=rbf(DELAY, DIM, varianz, &series[n],center[i]);
+        coefs_arr[0] += series[nst];
+        mat_arr[0][0] += 1.0;
+        for (octave_idx_type i=1;i<=CENTER;i++)
+          mat_arr[i][0] += hcen[i-1];
+        for (octave_idx_type i=1;i<=CENTER;i++) {
+          double h = hcen[i-1];
+          coefs_arr[i] += series[nst] * h;
+          for (octave_idx_type j=1;j<=i;j++)
+            mat_arr[i][j] += h*hcen[j-1];
+        }
+      }
+      
+      double h=(double)(INSAMPLE-STEP-(DIM-1)*DELAY);
+      for (octave_idx_type i=0;i<=CENTER;i++) {
+        coefs_arr[i] /= h;
+        for (octave_idx_type j=0;j<=i;j++) {
+          mat_arr[i][j] /= h;
+          mat_arr[j][i]=mat_arr[i][j];
+        }
+      }
 
-          // old solvele(mat_arr,coefs_arr, CENTER+1);
-          coefs = mat.solve(coefs);
-          coefs_arr = coefs.fortran_vec ();// coefs takes up new memory space
+      // old solvele(mat_arr,coefs_arr, CENTER+1);
+      coefs = mat.solve(coefs);
+      coefs_arr = coefs.fortran_vec ();// coefs takes up new memory space
 
-          // If solving the matrix generated errors do not continue
-          if (error_state)
-            return retval;
-
-          // end make_fit()
+      // end make_fit()
 
 
-          // Create outputs
-
-          // Create centers
-          Matrix centers (CENTER, DIM); 
-          for (octave_idx_type i=0;i<CENTER;i++) 
-              for (octave_idx_type j=0;j<DIM;j++)
-                {
-                  // old fprintf(stdout," %e",center[i][j]*interval+min_val);
-                  centers(i,j) = center[i][j]*interval+min_val;
-                }
+      // Create outputs
 
-          // Create variance
-          // old fprintf(stdout,"#variance= %e\n",varianz*interval);
-          double variance_val = varianz*interval;
-
-          // Create coefficients
-          NDArray coeff (dim_vector(CENTER +1, 1));
-          //  old fprintf(stdout,"#%e\n",coefs[0]*interval+min_val);
-          coeff(0) = coefs_arr[0]*interval+min_val;
-          for (octave_idx_type i=1;i<=CENTER;i++)
+      // Create centers
+      Matrix centers (CENTER, DIM); 
+      for (octave_idx_type i=0;i<CENTER;i++) 
+          for (octave_idx_type j=0;j<DIM;j++)
             {
-            // old fprintf(stdout,"#%e\n",coefs[i]*interval);
-              coeff(i) = coefs_arr[i]*interval;
+              // old fprintf(stdout," %e",center[i][j]*interval+min_val);
+              centers(i,j) = center[i][j]*interval+min_val;
             }
 
-          // Calculate insample error
-          double sigma = 0.0;
-          av = 0.0;
-          for (octave_idx_type i=0;i<INSAMPLE;i++) {
+      // Create variance
+      // old fprintf(stdout,"#variance= %e\n",varianz*interval);
+      double variance_val = varianz*interval;
+
+      // Create coefficients
+      NDArray coeff (dim_vector(CENTER +1, 1));
+      //  old fprintf(stdout,"#%e\n",coefs[0]*interval+min_val);
+      coeff(0) = coefs_arr[0]*interval+min_val;
+      for (octave_idx_type i=1;i<=CENTER;i++)
+        {
+        // old fprintf(stdout,"#%e\n",coefs[i]*interval);
+          coeff(i) = coefs_arr[i]*interval;
+        }
+
+      // Calculate insample error
+      double sigma = 0.0;
+      av = 0.0;
+      for (octave_idx_type i=0;i<INSAMPLE;i++) {
+        av += series[i];
+        sigma += series[i]*series[i];
+      }
+      av /= INSAMPLE;
+      sigma=sqrt(fabs(sigma/INSAMPLE-av*av));
+
+      // Create sample error
+      // 1st element of sample_error is the insample error
+      // 2nd element of sample_error is the out of sample error (if exists)
+      NDArray sample_error (dim_vector(1,1));
+      // old oldfprintf(stdout,"#insample error= %e\n",
+      //                forecast_error(0LU,INSAMPLE)/sigma);
+      sample_error(0) = forecast_error(series, center, coefs_arr, CENTER, 
+                                       DELAY, DIM, STEP, varianz,
+                                       0,INSAMPLE)
+                        / sigma;
+
+      if (INSAMPLE < LENGTH)
+        {
+          // Calculate out of sample error
+          av=sigma=0.0;
+          for (octave_idx_type i=INSAMPLE;i<LENGTH;i++) {
             av += series[i];
             sigma += series[i]*series[i];
           }
-          av /= INSAMPLE;
-          sigma=sqrt(fabs(sigma/INSAMPLE-av*av));
-
-          // Create sample error
-          // 1st element of sample_error is the insample error
-          // 2nd element of sample_error is the out of sample error (if exists)
-          NDArray sample_error (dim_vector(1,1));
-          // old oldfprintf(stdout,"#insample error= %e\n",
-          //                forecast_error(0LU,INSAMPLE)/sigma);
-          sample_error(0) = forecast_error(series, center, coefs_arr, CENTER, 
-                                           DELAY, DIM, STEP, varianz,
-                                           0,INSAMPLE)
-                            / sigma;
+          av /= (LENGTH-INSAMPLE);
+          sigma=sqrt(fabs(sigma/(LENGTH-INSAMPLE)-av*av));
 
-          if (INSAMPLE < LENGTH)
-            {
-              // Calculate out of sample error
-              av=sigma=0.0;
-              for (octave_idx_type i=INSAMPLE;i<LENGTH;i++) {
-                av += series[i];
-                sigma += series[i]*series[i];
-              }
-              av /= (LENGTH-INSAMPLE);
-              sigma=sqrt(fabs(sigma/(LENGTH-INSAMPLE)-av*av));
+          // old fprintf(stdout,"#out of sample error= %e\n",
+          //             forecast_error(INSAMPLE,LENGTH)/sigma);
+          sample_error.resize (dim_vector(2,1));
+          sample_error(1) = forecast_error(series, center, coefs_arr,
+                                           CENTER, DELAY, DIM, STEP,
+                                           varianz, INSAMPLE, LENGTH)
+                            / sigma;
+        }
 
-              // old fprintf(stdout,"#out of sample error= %e\n",
-              //             forecast_error(INSAMPLE,LENGTH)/sigma);
-              sample_error.resize (dim_vector(2,1));
-              sample_error(1) = forecast_error(series, center, coefs_arr,
-                                               CENTER, DELAY, DIM, STEP,
-                                               varianz, INSAMPLE, LENGTH)
-                                / sigma;
-            }
+      // Create forecast if MAKECAST == true
+      NDArray forecast (dim_vector (0,0));
+      if (MAKECAST)
+        {
+        // old make_cast ();
+
+          forecast.resize(dim_vector(CLENGTH,1));
+
+          octave_idx_type dim=(DIM-1)*DELAY;
 
-          // Create forecast if MAKECAST == true
-          NDArray forecast (dim_vector (0,0));
-          if (MAKECAST)
-            {
-            // old make_cast ();
-
-              forecast.resize(dim_vector(CLENGTH,1));
-
-              octave_idx_type dim=(DIM-1)*DELAY;
-
-              OCTAVE_LOCAL_BUFFER (double, cast, dim + 1);
-              for (octave_idx_type i=0;i<=dim;i++)
-                cast[i]=series[LENGTH-1-dim+i];
+          OCTAVE_LOCAL_BUFFER (double, cast, dim + 1);
+          for (octave_idx_type i=0;i<=dim;i++)
+            cast[i]=series[LENGTH-1-dim+i];
 
-              for (octave_idx_type n=0;n<CLENGTH;n++)
-                {
-                  double new_el=coefs_arr[0];
-                  for (octave_idx_type i=1;i<=CENTER;i++)
-                    new_el += coefs_arr[i]*rbf(DELAY,DIM,varianz,&cast[dim],
-                                           center[i-1]);
-                //  old  fprintf(out,"%e\n",new_el*interval+min_val);
-                  forecast(n) = new_el*interval+min_val;
-                  for (octave_idx_type i=0;i<dim;i++)
-                    cast[i]=cast[i+1];
-                  cast[dim]=new_el;
-                }
+          for (octave_idx_type n=0;n<CLENGTH;n++)
+            {
+              double new_el=coefs_arr[0];
+              for (octave_idx_type i=1;i<=CENTER;i++)
+                new_el += coefs_arr[i]*rbf(DELAY,DIM,varianz,&cast[dim],
+                                       center[i-1]);
+            //  old  fprintf(out,"%e\n",new_el*interval+min_val);
+              forecast(n) = new_el*interval+min_val;
+              for (octave_idx_type i=0;i<dim;i++)
+                cast[i]=cast[i+1];
+              cast[dim]=new_el;
             }
+        }
 
-          // Create output
-          retval(0) = centers; 
-          retval(1) = variance_val; // variance;
-          retval(2) = coeff; 
-          retval(3) = sample_error; // sample error [in sample; out of sample];
-          retval(4) = forecast; // forecast values;
-        }
+      // Create output
+      retval(0) = centers; 
+      retval(1) = variance_val; // variance;
+      retval(2) = coeff; 
+      retval(3) = sample_error; // sample error [in sample; out of sample];
+      retval(4) = forecast; // forecast values;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__surrogates__.cc
--- a/src/__surrogates__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__surrogates__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -70,34 +70,30 @@
       F77_INT ispec         = TO_F77_INT (args(3).idx_type_value ());
       double seed           = args(4).double_value ();
 
-      if (! error_state)
-        {
+      F77_INT nmaxp = TO_F77_INT (input.rows ());
+      F77_INT mcmax = TO_F77_INT (input.columns ());
 
-          F77_INT nmaxp = TO_F77_INT (input.rows ());
-          F77_INT mcmax = TO_F77_INT (input.columns ());
+      Cell surro_data (dim_vector (nsur,1));
+      Matrix surro_tmp (input.dims ());
+      Matrix pars (nsur, 2);
 
-          Cell surro_data (dim_vector (nsur,1));
-          Matrix surro_tmp (input.dims ());
-          Matrix pars (nsur, 2);
+      for (octave_idx_type i = 0; i < nsur; i++)
+        {
+          F77_INT it_tmp;
+          double rel_discrepency_tmp;
 
-          for (octave_idx_type i = 0; i < nsur; i++)
-            {
-              F77_INT it_tmp;
-              double rel_discrepency_tmp;
+          F77_XFCN (ts_surrogates, TS_SURROGATES,
+                    (input.fortran_vec (), nmaxp, mcmax, imax, ispec,
+                     seed, surro_tmp.fortran_vec (), it_tmp,
+                     rel_discrepency_tmp));
 
-              F77_XFCN (ts_surrogates, TS_SURROGATES,
-                        (input.fortran_vec (), nmaxp, mcmax, imax, ispec,
-                         seed, surro_tmp.fortran_vec (), it_tmp,
-                         rel_discrepency_tmp));
+          surro_data(i) = surro_tmp;
+          pars(i,0)     = it_tmp;
+          pars(i,1)     = rel_discrepency_tmp;
+        }
 
-              surro_data(i) = surro_tmp;
-              pars(i,0)     = it_tmp;
-              pars(i,1)     = rel_discrepency_tmp;
-            }
-
-          retval(0) = surro_data;
-          retval(1) = pars;
-        }
+      retval(0) = surro_data;
+      retval(1) = pars;
 
     }
   return retval;
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__upo__.cc
--- a/src/__upo__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__upo__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -80,31 +80,26 @@
       F77_INT icen    = TO_F77_INT (args(9).int_value());
 
 
-
-      if (! error_state)
-        {
+      F77_INT lines_read = TO_F77_INT (in_out1.numel());
+      // Generating output vectors with estimated lengths
+      // The extra length (+1) is to store the actual lengths
+      NDArray olens (dim_vector (icen+1,1));
+      NDArray orbit_data (dim_vector ( (icen*20<lines_read?\
+                                                icen*20:lines_read)+1, 1));
+      NDArray acc (dim_vector (icen+1,1));
+      NDArray stability (dim_vector (icen+1,1));
 
-          F77_INT lines_read = TO_F77_INT (in_out1.numel());
-          // Generating output vectors with estimated lengths
-          // The extra length (+1) is to store the actual lengths
-          NDArray olens (dim_vector (icen+1,1));
-          NDArray orbit_data (dim_vector ( (icen*20<lines_read?\
-                                                    icen*20:lines_read)+1, 1));
-          NDArray acc (dim_vector (icen+1,1));
-          NDArray stability (dim_vector (icen+1,1));
-
-          F77_XFCN (ts_upo, TS_UPO,
-                    (m, eps, frac,teq, tdis, h, tacc, iper,icen, lines_read,
-                     in_out1.fortran_vec(), olens.fortran_vec(), 
-                     orbit_data.fortran_vec(), orbit_data.numel(), 
-                     acc.fortran_vec(), stability.fortran_vec()));
+      F77_XFCN (ts_upo, TS_UPO,
+                (m, eps, frac,teq, tdis, h, tacc, iper,icen, lines_read,
+                 in_out1.fortran_vec(), olens.fortran_vec(), 
+                 orbit_data.fortran_vec(), orbit_data.numel(), 
+                 acc.fortran_vec(), stability.fortran_vec()));
 
 
-          retval(0) = olens;
-          retval(1) = orbit_data;
-          retval(2) = acc;
-          retval(3) = stability;
-        }
+      retval(0) = olens;
+      retval(1) = orbit_data;
+      retval(2) = acc;
+      retval(3) = stability;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/__xzero__.cc
--- a/src/__xzero__.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/__xzero__.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -111,38 +111,35 @@
       double epsilon=EPS0/EPSF;
       octave_idx_type clength=(CLENGTH <= LENGTH) ? CLENGTH-STEP : LENGTH-STEP;
 
-      if (! error_state)
+      // Calculate fit
+      while (!alldone) {
+        alldone=1;
+        epsilon*=EPSF;
+        make_box(series1,box,list,LENGTH-STEP,NMAX,DIM,DELAY,epsilon);
+        for (octave_idx_type i=(DIM-1)*DELAY;i<clength;i++)
+          if (!done[i]) {
+            octave_idx_type actfound;
+            actfound=find_neighbors(series1,box,list,series2+i,LENGTH,NMAX,
+                                    DIM,DELAY,epsilon,found);
+            if (actfound >= MINN) {
+              for (octave_idx_type j=1;j<=STEP;j++)
+                error_array[j-1] += make_fit (series1, series2, found,
+                                              i,actfound,j);
+              done[i]=1;
+            }
+            alldone &= done[i];
+          }
+      }
+
+      // Create output
+      Matrix output (STEP,1);
+      for (octave_idx_type i=0;i<STEP;i++)
         {
-          // Calculate fit
-          while (!alldone) {
-            alldone=1;
-            epsilon*=EPSF;
-            make_box(series1,box,list,LENGTH-STEP,NMAX,DIM,DELAY,epsilon);
-            for (octave_idx_type i=(DIM-1)*DELAY;i<clength;i++)
-              if (!done[i]) {
-                octave_idx_type actfound;
-                actfound=find_neighbors(series1,box,list,series2+i,LENGTH,NMAX,
-                                        DIM,DELAY,epsilon,found);
-                if (actfound >= MINN) {
-                  for (octave_idx_type j=1;j<=STEP;j++)
-                    error_array[j-1] += make_fit (series1, series2, found,
-                                                  i,actfound,j);
-                  done[i]=1;
-                }
-                alldone &= done[i];
-              }
-          }
-
-          // Create output
-          Matrix output (STEP,1);
-          for (octave_idx_type i=0;i<STEP;i++)
-            {
-//            old fprintf(stdout,"%lu %e\n",i+1,
-//                      sqrt(error_array[i]/(clength-(DIM-1)*DELAY))/rms2);
-              output(i,0) = sqrt(error_array[i]/(clength-(DIM-1)*DELAY))/rms2;
-            }
-          retval(0) = output;
+          // old fprintf(stdout,"%lu %e\n",i+1,
+          //           sqrt(error_array[i]/(clength-(DIM-1)*DELAY))/rms2);
+          output(i,0) = sqrt(error_array[i]/(clength-(DIM-1)*DELAY))/rms2;
         }
+      retval(0) = output;
 
     }
   return retval;
diff -r 71f2c8fde0c5 -r e40a599d68cf src/lazy.cc
--- a/src/lazy.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/lazy.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -140,34 +140,31 @@
                        "Number of iterations (IMAX) must be a positive "
                        "integer");
 
-      if (! error_state)
-        {
-        // If vector is in 1 row: transpose (we will transpose the output to fit)
-            transposed = 0;
+      // If vector is in 1 row: transpose (we will transpose the output to fit)
+      transposed = 0;
 
-          if ((rows == 1) && (cols > 1))
-            {
-              transposed = 1;
-              in_out1    = in_out1.transpose();
-            }
+      if ((rows == 1) && (cols > 1))
+        {
+          transposed = 1;
+          in_out1    = in_out1.transpose();
+        }
 
-          F77_INT lines_read = TO_F77_INT (in_out1.numel());
-          NDArray in_out2 (Matrix (lines_read, 1));
+      F77_INT lines_read = TO_F77_INT (in_out1.numel());
+      NDArray in_out2 (Matrix (lines_read, 1));
 
-          F77_XFCN (ts_lazy, TS_LAZY,
-                    (m, rv, imax, lines_read,
-                      in_out1.fortran_vec(), in_out2.fortran_vec()));
+      F77_XFCN (ts_lazy, TS_LAZY,
+                (m, rv, imax, lines_read,
+                  in_out1.fortran_vec(), in_out2.fortran_vec()));
 
-          // Transpose the output to resemble the input
-          if (transposed)
-          {
-            in_out1 = in_out1.transpose();
-            in_out2 = in_out2.transpose();
-          }
+      // Transpose the output to resemble the input
+      if (transposed)
+      {
+        in_out1 = in_out1.transpose();
+        in_out2 = in_out2.transpose();
+      }
 
-          retval(0) = in_out1;
-          retval(1) = in_out2;
-        }
+      retval(0) = in_out1;
+      retval(1) = in_out2;
     }
   return retval;
 }
diff -r 71f2c8fde0c5 -r e40a599d68cf src/mutual.cc
--- a/src/mutual.cc	Mon Aug 26 12:51:20 2019 -0400
+++ b/src/mutual.cc	Mon Nov 29 13:43:29 2021 +0100
@@ -193,60 +193,57 @@
       int32NDArray h2 (dim_vector(partitions, partitions));
       OCTAVE_LOCAL_BUFFER (long, array, length);
 
-      if (! error_state)
-        {
-          // Load array
+      // Load array
 
-          // Rescale data and load array
-          // NOTE: currently supports only vectors so (dim == 1) always
-          if (dim == 1){
-          double mint, interval;
-          rescale_data(input,0,length,&mint,&interval);
+      // Rescale data and load array
+      // NOTE: currently supports only vectors so (dim == 1) always
+      if (dim == 1){
+      double mint, interval;
+      rescale_data(input,0,length,&mint,&interval);
 
-          for (octave_idx_type i=0;i<length;i++)
-            if (input(i,0) < 1.0)
-              array[i]=(long)(input(i,0)*(double)partitions);
-            else
-              array[i]=partitions-1;
-          }
+      for (octave_idx_type i=0;i<length;i++)
+        if (input(i,0) < 1.0)
+          array[i]=(long)(input(i,0)*(double)partitions);
+        else
+          array[i]=partitions-1;
+      }
 
-          double shannon = make_cond_entropy (0, h1, h11, h2, array, length,
-                                              partitions);
-          if (corrlength >= (long)length)
-            corrlength=length-1;
+      double shannon = make_cond_entropy (0, h1, h11, h2, array, length,
+                                          partitions);
+      if (corrlength >= (long)length)
+        corrlength=length-1;
 
-          // Construct the output
-          Matrix delay (corrlength+1,1);
-          // To save memory
-          int minf_len = 1;
+      // Construct the output
+      Matrix delay (corrlength+1,1);
+      // To save memory
+      int minf_len = 1;
 
-          if (nargout > 1)
-            minf_len = corrlength+1;
-          Matrix mutual_inf (minf_len,1);
+      if (nargout > 1)
+        minf_len = corrlength+1;
+      Matrix mutual_inf (minf_len,1);
 
-          // Assign output
-          delay(0,0)      = 0;
-          mutual_inf(0,0) = shannon;
-          for (octave_idx_type tau=1;tau<=corrlength;tau++) {
+      // Assign output
+      delay(0,0)      = 0;
+      mutual_inf(0,0) = shannon;
+      for (octave_idx_type tau=1;tau<=corrlength;tau++) {
 
-  //          fprintf(stdout,"%ld %e %e\n",tau,condent,condent/log((double)partitions));
-            delay(tau,0) = tau;
-            if (nargout > 1)
-              {
-                mutual_inf(tau,0) = make_cond_entropy(tau, h1, h11, h2, array,
-                                                      length, partitions);
-              }
+        // fprintf(stdout,"%ld %e %e\n",tau,condent,condent/log((double)partitions));
+        delay(tau,0) = tau;
+        if (nargout > 1)
+          {
+            mutual_inf(tau,0) = make_cond_entropy(tau, h1, h11, h2, array,
+                                                  length, partitions);
           }
+      }
 
-          if (transposed)
-            {
-              delay          = delay.transpose();
-              if (nargout > 1)
-                mutual_inf     = mutual_inf.transpose();
-            }
-          retval(0) = delay;
-          retval(1) = mutual_inf;
+      if (transposed)
+        {
+          delay          = delay.transpose();
+          if (nargout > 1)
+            mutual_inf     = mutual_inf.transpose();
         }
+      retval(0) = delay;
+      retval(1) = mutual_inf;
     }
   return retval;
 }
