Sign Up
Log In
Log In
or
Sign Up
Places
All Projects
Status Monitor
Collapse sidebar
openSUSE:Factory:Rebuild
octave-forge-tisean
tisean-drop-error_state-use.patch
Overview
Repositories
Revisions
Requests
Users
Attributes
Meta
File tisean-drop-error_state-use.patch of Package octave-forge-tisean
# 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 # Updated by Marius Schamschula <mps@macports.org> # Date 2023-03-10 Remove usage of `error_state` (bug #61583). * 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.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src//__c1__.cc 2023-03-09 18:58:35.000000000 -0600 @@ -78,65 +78,61 @@ octave_idx_type iverb = verbose; - if (! error_state) - { + octave_idx_type lines_read = input.rows (); //nmax in d1() + octave_idx_type columns_read = input.columns (); - octave_idx_type lines_read = input.rows (); //nmax in d1() - octave_idx_type columns_read = 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); - - // Seed the rand() function for d1() - F77_XFCN (rand, RAND, (sqrt(seed))); - - for (octave_idx_type m = mindim; m <= maxdim; m++) - { - octave_scalar_map tmp (keys); - tmp.setfield ("dim", 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) - { - 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)); - - 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); - - output.assign (idx_vector(m-mindim), tmp); - } - - retval(0) = output; - } - } - return retval; + 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))); + + for (octave_idx_type m = mindim; m <= maxdim; m++) + { + octave_scalar_map tmp (keys); + tmp.setfield ("dim", 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) + { + 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)); + + 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); + + output.assign (idx_vector(m-mindim), tmp); + } + + 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.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/__lfo_ar__.cc 2023-03-10 05:42:42.000000000 -0600 @@ -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,100 @@ 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"); - - // 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; - - if (verbose) - printf("\nStarting new dataset\n\n"); - - for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) - { - if (verbose) - { - printf ("For epsilon = %e, the count = %d\n", - epsilon,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; - - 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 - // 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]; - } - } - } - 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++) - { - //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; - } - } - // 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); + // 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; + + // Create output + Matrix output (output_rows,dim+4); + octave_idx_type count = 0; + + if (verbose) + printf("\nStarting new dataset\n\n"); + + for (double epsilon=EPS0;epsilon<EPS1*EPSF;epsilon*=EPSF) + { + if (verbose) + { + printf ("For epsilon = %e, the count = %d\n", + epsilon,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; + + 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); + 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]; + } + } + } + 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++) + { + //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; + } + } + // 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/__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.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/__lzo_run__.cc 2023-03-09 19:11:14.000000000 -0600 @@ -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); - - // old indexes=make_multi_index(dim,embed,DELAY); - - 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; - } - - // end old index = make_multi_index(); - - // old rnd_init(seed); - TISEAN_rand generator (seed); - - double epsilon0=EPS0/EPSF; - - 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); - - 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]; - } - } - } - - if (row_count != 0) - { - output.resize (row_count, dim); - retval(0) = output; - } - else - retval(0) = Matrix (0,0); - } - } + 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); + + 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; + } + + // end old index = make_multi_index(); + + // old rnd_init(seed); + TISEAN_rand generator (seed); + + double epsilon0=EPS0/EPSF; + + 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); + + 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]; + } + } + } + + 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.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/__surrogates__.cc 2023-03-09 19:13:16.000000000 -0600 @@ -63,35 +63,30 @@ octave_idx_type ispec = args(3).idx_type_value (); double seed = args(4).double_value (); - if (! error_state) - { + octave_idx_type nmaxp = input.rows (); + octave_idx_type mcmax = input.columns (); - octave_idx_type nmaxp = input.rows (); - octave_idx_type mcmax = input.columns (); - - 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++) - { - octave_idx_type 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)); - - surro_data(i) = surro_tmp; - pars(i,0) = it_tmp; - pars(i,1) = rel_discrepency_tmp; - } - - retval(0) = surro_data; - retval(1) = pars; - } - - } + 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++) + { + octave_idx_type 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)); + + surro_data(i) = surro_tmp; + pars(i,0) = it_tmp; + pars(i,1) = rel_discrepency_tmp; + } + + retval(0) = surro_data; + retval(1) = pars; + } return retval; } diff -r 71f2c8fde0c5 -r e40a599d68cf src/__upo__.cc --- a/src/__upo__.cc.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/__upo__.cc 2023-03-09 19:15:03.000000000 -0600 @@ -74,31 +74,26 @@ int icen = args(9).int_value(); - - if (! error_state) - { - - int lines_read = 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())); - - - retval(0) = olens; - retval(1) = orbit_data; - retval(2) = acc; - retval(3) = stability; - } - } + int lines_read = 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())); + + + 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.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/lazy.cc 2023-03-09 19:16:23.000000000 -0600 @@ -133,36 +133,33 @@ "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 ((rows == 1) && (cols > 1)) - { - transposed = 1; - in_out1 = in_out1.transpose(); - } - - int lines_read = 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())); - - // 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; - } - } - return retval; + // 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(); + } + + int lines_read = 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())); + + // 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; + } + 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; } --- a/src/__lfo_test__.cc.orig 2015-08-14 17:25:52.000000000 -0500 +++ b/src/__lfo_test__.cc 2023-03-09 19:39:33.000000000 -0600 @@ -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,80 @@ 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"); - - // 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; - - // 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]); - - 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)); - - // 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++) - { - 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]; - } - } - - retval(0) = rel; - retval(1) = ind; - } - } + // 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 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]) + { + + 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); + 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]; + } + } + + 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 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++) + { + // old fprintf(stdout,"%e ",individual[j][i]); + ind(i-(EMBED-1)*DELAY, j) = individual[j][i]; + } + } + + retval(0) = rel; + retval(1) = ind; + } return retval; }
Locations
Projects
Search
Status Monitor
Help
OpenBuildService.org
Documentation
API Documentation
Code of Conduct
Contact
Support
@OBShq
Terms
openSUSE Build Service is sponsored by
The Open Build Service is an
openSUSE project
.
Sign Up
Log In
Places
Places
All Projects
Status Monitor
