forked from pool/octave-forge-tisean
3726 lines
134 KiB
Diff
3726 lines
134 KiB
Diff
|
# 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;
|
||
|
|
}
|