# HG changeset patch # User Markus Mützel # Date 1638189809 -3600 # Mon Nov 29 13:43:29 2021 +0100 # Node ID e40a599d68cf3061d04a9dac30e36751ed20acb2 # Parent 71f2c8fde0c59f4942dd22da89bba162e8ae6c97 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 --- src/__boxcount__.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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 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::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 histo (maxembed * dimension, 0.0); - start_box(series, which_dims, histo, maxembed, dimension, DELAY, - length, epsi, Q); - - if (Q != 1.0) - for (std::vector::iterator it = histo.begin (); - it != histo.end (); it++) - *it = log(*it)/(1.0-Q); - - histo_list.push_back (histo); + std::list>::const_iterator it_hist; + it_hist = histo_list.cbegin (); + for (octave_idx_type j=0;jhist[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>::const_iterator it_hist; - it_hist = histo_list.cbegin (); - for (octave_idx_type j=0;jhist[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 --- src//__c1__.cc.orig 2015-08-14 17:25:52.000000000 -0500 +++ 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 --- src/__d2__.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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 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=0)?sn-MINDIST:0; - octave_idx_type n2=(sn+MINDIST 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 WHEN) || (n == (nmax-1)) || - (imin > (HOWOFTEN-1)) || (counter % it_pause == 0)) - { - time(&lasttime); - - if (imin > (HOWOFTEN-1)) + for (octave_idx_type i1=0;i1 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 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 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 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 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 --- src/__delay__.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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 --- src/__false_nearest__.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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= 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= 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 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 1) - { - double sumerror=0.0; - for (octave_idx_type j=0;j 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 1) + { + double sumerror=0.0; + for (octave_idx_type j=0;j 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 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= MINN) + done=1; + for (octave_idx_type j=0;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 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 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 --- src/__lyap_spec__.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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 - ::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 + ::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 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 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 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 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 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 1) + { + double sumerror=0.0; + for (octave_idx_type j=0;j box (dim_vector(NMAX,NMAX)); - if ( ! error_state) - { - for (octave_idx_type j=0;j indexes (dim_vector (alldim, 2)); - for (octave_idx_type i=0;i= 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 indexes (dim_vector (alldim, 2)); + for (octave_idx_type i=0;i= 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 box (dim_vector(NMAX, NMAX)); // Compute forecast error - if ( ! error_state) - { - for (octave_idx_type i=0;i= 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 1) + { + ind_forecast_err.resize(clength - (embed-1)*DELAY, dim); + for (octave_idx_type i=(embed-1)*DELAY;i= 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 1) - { - ind_forecast_err.resize(clength - (embed-1)*DELAY, dim); - for (octave_idx_type i=(embed-1)*DELAY;i 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 --- src/__rbf__.cc Mon Aug 26 12:51:20 2019 -0400 +++ src/__rbf__.cc Mon Nov 29 13:43:29 2021 +0100 @@ -152,169 +152,161 @@ for (octave_idx_type j=0;j= 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= 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 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 --- src/mutual.cc Mon Aug 26 12:51:20 2019 -0400 +++ 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= (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; } --- src/__lfo_test__.cc.orig 2015-08-14 17:25:52.000000000 -0500 +++ 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 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 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 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 MINN) + { + make_fit(series, indexes, found, STEP, DIM, COMP, + actfound,i,newcast); + for (octave_idx_type j=0;j