← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/__lyap_spec__.cc

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/__lyap_spec__.cc
 
1
/*
 
2
 *   This file is part of TISEAN
 
3
 *
 
4
 *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
 
5
 *
 
6
 *   TISEAN is free software; you can redistribute it and/or modify
 
7
 *   it under the terms of the GNU General Public License as published by
 
8
 *   the Free Software Foundation; either version 2 of the License, or
 
9
 *   (at your option) any later version.
 
10
 *
 
11
 *   TISEAN is distributed in the hope that it will be useful,
 
12
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
 *   GNU General Public License for more details.
 
15
 *
 
16
 *   You should have received a copy of the GNU General Public License
 
17
 *   along with TISEAN; if not, write to the Free Software
 
18
 *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
19
 */
 
20
/* Author: Rainer Hegger 
 
21
 * Modified: Piotr Held <pjheld@gmail.com> (2015).
 
22
 * This function is based on lyap_spec of TISEAN 3.0.1 
 
23
 * https://github.com/heggus/Tisean"
 
24
 */
 
25
 
 
26
#define HELPTEXT "Part of tisean package\n\
 
27
No argument checking\n\
 
28
FOR INTERNAL USE ONLY"
 
29
 
 
30
#include <octave/oct.h>
 
31
#include <octave/oct-map.h>
 
32
#include "routines_c/tsa.h"
 
33
#include <limits>
 
34
 
 
35
#define OUT 10
 
36
 
 
37
#define BOX 512
 
38
#define EPSMAX 1.0
 
39
// the delay vector is calculated in the associated m-file
 
40
// according to authors of TISEAN the DELAY > 1 did not work for some reason
 
41
#define DELAY 1
 
42
 
 
43
double sort(const double **series, unsigned long *found,
 
44
            octave_idx_type **indexes, double epsmin, bool epsset,
 
45
            octave_idx_type alldim, octave_idx_type MINNEIGHBORS,
 
46
            octave_idx_type LENGTH, long act,
 
47
            octave_idx_type *nfound, bool *enough)
 
48
{
 
49
  double maxeps=0.0,dx,dswap,maxdx;
 
50
  long self=0,i,j,del,hf,iswap,n1;
 
51
  octave_idx_type imax = *nfound;
 
52
 
 
53
  *enough=0;
 
54
 
 
55
  OCTAVE_LOCAL_BUFFER (double, abstand, LENGTH);
 
56
 
 
57
  for (i=0;i<imax;i++) {
 
58
    hf=found[i];
 
59
    if (hf != act) {
 
60
      maxdx=fabs(series[0][act]-series[0][hf]);
 
61
      for (j=1;j<alldim;j++) {
 
62
        n1=indexes[0][j];
 
63
        del=indexes[1][j];
 
64
        dx=fabs(series[n1][act-del]-series[n1][hf-del]);
 
65
        if (dx > maxdx) maxdx=dx;
 
66
      }
 
67
      abstand[i]=maxdx;
 
68
    }
 
69
    else {
 
70
      self=i;
 
71
    }
 
72
  }
 
73
 
 
74
  if (self != (imax-1)) {
 
75
    abstand[self]=abstand[imax-1];
 
76
    found[self]=found[imax-1];
 
77
  }
 
78
 
 
79
  for (i=0;i<MINNEIGHBORS;i++) {
 
80
    for (j=i+1;j<imax-1;j++) {
 
81
      if (abstand[j]<abstand[i]) {
 
82
        dswap=abstand[i];
 
83
        abstand[i]=abstand[j];
 
84
        abstand[j]=dswap;
 
85
        iswap=found[i];
 
86
        found[i]=found[j];
 
87
        found[j]=iswap;
 
88
      }
 
89
    }
 
90
  }
 
91
 
 
92
  if (!epsset || (abstand[MINNEIGHBORS-1] >= epsmin)) {
 
93
    *nfound=MINNEIGHBORS;
 
94
    *enough=1;
 
95
    maxeps=abstand[MINNEIGHBORS-1];
 
96
 
 
97
    return maxeps;
 
98
  }
 
99
 
 
100
  for (i=MINNEIGHBORS;i<imax-2;i++) {
 
101
    for (j=i+1;j<imax-1;j++) {
 
102
      if (abstand[j]<abstand[i]) {
 
103
        dswap=abstand[i];
 
104
        abstand[i]=abstand[j];
 
105
        abstand[j]=dswap;
 
106
        iswap=found[i];
 
107
        found[i]=found[j];
 
108
        found[j]=iswap;
 
109
      }
 
110
    }
 
111
    if (abstand[i] > epsmin) {
 
112
      (*nfound)=i+1;
 
113
      *enough=1;
 
114
      maxeps=abstand[i];
 
115
 
 
116
      return maxeps;
 
117
    }
 
118
  }
 
119
 
 
120
  maxeps=abstand[imax-2];
 
121
 
 
122
  return maxeps;
 
123
}
 
124
 
 
125
void make_dynamics(const double **series, octave_idx_type **box,
 
126
                   octave_idx_type **indexes, double &epsmin, bool epsset,
 
127
                   double EPSSTEP, octave_idx_type EMBED,
 
128
                   octave_idx_type MINNEIGHBORS, octave_idx_type LENGTH,
 
129
                   octave_idx_type DIMENSION,
 
130
                   octave_idx_type count, double &avneig, double &aveps,
 
131
                   double **dynamics, double *averr, octave_idx_type act)
 
132
{
 
133
 
 
134
  OCTAVE_LOCAL_BUFFER (unsigned long, found, LENGTH);
 
135
 
 
136
  OCTAVE_LOCAL_BUFFER (const double*, hser, DIMENSION);
 
137
  for (octave_idx_type i=0;i<DIMENSION;i++)
 
138
    hser[i]=series[i]+act;
 
139
 
 
140
  octave_idx_type alldim = DIMENSION * EMBED;
 
141
 
 
142
  double epsilon=epsmin/EPSSTEP;
 
143
  double foundeps = 0.0;
 
144
  octave_idx_type nfound=0;
 
145
  do {
 
146
    epsilon *= EPSSTEP;
 
147
    if (epsilon > EPSMAX)
 
148
      epsilon=EPSMAX;
 
149
    OCTAVE_LOCAL_BUFFER (long, list, LENGTH);
 
150
    make_multi_box(series,box,list,LENGTH-DELAY,BOX,DIMENSION,EMBED,
 
151
                   DELAY,epsilon);
 
152
    nfound=find_multi_neighbors(series,box,list,hser,BOX,
 
153
                                DIMENSION,EMBED,DELAY,epsilon,found);
 
154
    if (nfound > MINNEIGHBORS) {
 
155
      bool got_enough;
 
156
      foundeps=sort(series, found, indexes, epsmin, epsset, alldim,
 
157
                    MINNEIGHBORS, LENGTH, act, &nfound, &got_enough);
 
158
      if (got_enough)
 
159
        break;
 
160
    }
 
161
  } while (epsilon < EPSMAX);
 
162
 
 
163
 
 
164
  avneig += nfound;
 
165
  aveps += foundeps;
 
166
  if (!epsset)
 
167
    epsmin=aveps/count;
 
168
 
 
169
  // If less neighbors found than the minimum number
 
170
  // No sense continuing
 
171
  if (nfound < MINNEIGHBORS)
 
172
    {
 
173
      error_with_id ("Octave:tisean","Not enough neighbors found");
 
174
      return ;
 
175
    }
 
176
 
 
177
  // *_arr signifies pointer to data in *
 
178
  // This is done for optimization
 
179
  Matrix vec (alldim + 1,1);
 
180
  vec.fill (0.0);
 
181
  double *vec_arr = vec.fortran_vec ();
 
182
 
 
183
  Matrix mat (alldim+1, alldim+1);
 
184
  mat.fill (0.0);
 
185
  OCTAVE_LOCAL_BUFFER (double *, mat_arr, alldim + 1);
 
186
  for (octave_idx_type i = 0; i < alldim + 1; i++)
 
187
    mat_arr[i] = mat.fortran_vec () + (alldim+1) * i;
 
188
  
 
189
  for (octave_idx_type i=0;i<nfound;i++) {
 
190
    octave_idx_type act=found[i];
 
191
    mat_arr[0][0] += 1.0;
 
192
    for (octave_idx_type j=0;j<alldim;j++)
 
193
      mat_arr[0][j+1] += series[indexes[0][j]][act-indexes[1][j]];
 
194
    for (octave_idx_type j=0;j<alldim;j++) {
 
195
      double hv1=series[indexes[0][j]][act-indexes[1][j]];
 
196
      octave_idx_type hj=j+1;
 
197
      for (octave_idx_type k=j;k<alldim;k++)
 
198
        mat_arr[hj][k+1] += series[indexes[0][k]][act-indexes[1][k]]*hv1;
 
199
    }
 
200
  }
 
201
 
 
202
  for (octave_idx_type i=0;i<=alldim;i++)
 
203
    for (octave_idx_type j=i;j<=alldim;j++)
 
204
      mat_arr[j][i]=(mat_arr[i][j]/=(double)nfound);
 
205
 
 
206
  for (octave_idx_type d=0;d<DIMENSION;d++)
 
207
    {
 
208
      for (octave_idx_type i=0;i<=alldim;i++)
 
209
        vec_arr[i]=0.0;
 
210
      for (octave_idx_type i=0;i<nfound;i++) {
 
211
        octave_idx_type act=found[i];
 
212
        double hv=series[d][act+DELAY];
 
213
        vec_arr[0] += hv;
 
214
        for (octave_idx_type j=0;j<alldim;j++)
 
215
          vec_arr[j+1] += hv*series[indexes[0][j]][act-indexes[1][j]];
 
216
      }
 
217
      for (octave_idx_type i=0;i<=alldim;i++)
 
218
        vec_arr[i] /= (double)nfound;
 
219
 
 
220
      Matrix solved      = mat.solve (vec);
 
221
      double *solved_arr = solved.fortran_vec ();
 
222
 
 
223
    // If errors were raised (a singular matrix was encountered), 
 
224
    // there is no sense in countinuing
 
225
    if (error_state)
 
226
      {
 
227
        return ;
 
228
      }
 
229
 
 
230
      double new_vec = solved_arr[0];
 
231
      for (octave_idx_type i=1;i<=alldim;i++)
 
232
        dynamics[d][i-1] = solved_arr[i];
 
233
 
 
234
      for (octave_idx_type i=0;i<alldim;i++)
 
235
        new_vec += dynamics[d][i]*series[indexes[0][i]][act-indexes[1][i]];
 
236
      averr[d]  += (new_vec-series[d][act+DELAY])
 
237
                   *(new_vec-series[d][act+DELAY]);
 
238
    }
 
239
 
 
240
}
 
241
 
 
242
void gram_schmidt(octave_idx_type alldim, double **delta,
 
243
                  double *stretch)
 
244
{
 
245
 
 
246
  OCTAVE_LOCAL_BUFFER (double, diff, alldim);
 
247
  OCTAVE_LOCAL_BUFFER (double *, dnew, alldim);
 
248
  OCTAVE_LOCAL_BUFFER (double, dnew_data, alldim * alldim);
 
249
  for (octave_idx_type i=0;i<alldim;i++)
 
250
    dnew[i]=dnew_data + alldim * i;
 
251
 
 
252
  for (octave_idx_type i=0;i<alldim;i++) {
 
253
    for (octave_idx_type j=0;j<alldim;j++) 
 
254
      diff[j]=0.0;
 
255
    for (octave_idx_type j=0;j<i;j++) {
 
256
      double norm=0.0;
 
257
      for (octave_idx_type k=0;k<alldim;k++)
 
258
        norm += delta[i][k]*dnew[j][k];
 
259
      for (octave_idx_type k=0;k<alldim;k++)
 
260
        diff[k] -= norm*dnew[j][k];
 
261
    }
 
262
    double norm=0.0;
 
263
    for (octave_idx_type j=0;j<alldim;j++)
 
264
      norm += sqr(delta[i][j]+diff[j]);
 
265
    stretch[i]=(norm=sqrt(norm));
 
266
    for (octave_idx_type j=0;j<alldim;j++)
 
267
      dnew[i][j]=(delta[i][j]+diff[j])/norm;
 
268
  }
 
269
  for (octave_idx_type i=0;i<alldim;i++)
 
270
    for (octave_idx_type j=0;j<alldim;j++)
 
271
      delta[i][j]=dnew[i][j];
 
272
 
 
273
}
 
274
 
 
275
void make_iteration(octave_idx_type DIMENSION, octave_idx_type alldim,
 
276
                    double **dynamics, double **delta)
 
277
{
 
278
 
 
279
  OCTAVE_LOCAL_BUFFER (double *, dnew, alldim);
 
280
  OCTAVE_LOCAL_BUFFER (double, dnew_data, alldim * alldim);
 
281
  for (octave_idx_type i=0;i<alldim;i++)
 
282
    dnew[i]=dnew_data + alldim * i;
 
283
 
 
284
  for (octave_idx_type i=0;i<alldim;i++) {
 
285
    for (octave_idx_type j=0;j<DIMENSION;j++) {
 
286
      dnew[i][j]=dynamics[j][0]*delta[i][0];
 
287
      for (octave_idx_type k=1;k<alldim;k++)
 
288
        dnew[i][j] += dynamics[j][k]*delta[i][k];
 
289
    }
 
290
    for (octave_idx_type j=DIMENSION;j<alldim;j++)
 
291
      dnew[i][j]=delta[i][j-1];
 
292
  }
 
293
 
 
294
  for (octave_idx_type i=0;i<alldim;i++)
 
295
    for (octave_idx_type j=0;j<alldim;j++)
 
296
      delta[i][j]=dnew[i][j];
 
297
 
 
298
}
 
299
 
 
300
DEFUN_DLD (__lyap_spec__, args, nargout, HELPTEXT)
 
301
{
 
302
 
 
303
  int nargin = args.length ();
 
304
  octave_value_list retval;
 
305
 
 
306
  if ((nargin != 8 && nargin != 13) || nargout != 2)
 
307
  {
 
308
    print_usage ();
 
309
  }
 
310
  else
 
311
    {
 
312
 
 
313
      // Assign input
 
314
      Matrix input                 = args(0).matrix_value ();
 
315
      octave_idx_type EMBED        = args(1).idx_type_value ();
 
316
      octave_idx_type ITERATIONS   = args(2).idx_type_value ();
 
317
      double epsmin                = args(3).double_value ();
 
318
      bool epsset                  = args(4).bool_value ();
 
319
      double EPSSTEP               = args(5).double_value ();
 
320
      octave_idx_type MINNEIGHBORS = args(6).idx_type_value ();
 
321
      octave_idx_type it_pause     = args(7).idx_type_value ();
 
322
 
 
323
      octave_idx_type count = 0;
 
324
      if (nargin == 13)
 
325
        {
 
326
          count        = args(8).idx_type_value ();
 
327
          if (count > (ITERATIONS - (EMBED-1)*DELAY))
 
328
            error_with_id ("Octave:tisean", "The variable 'count' grew too "
 
329
                                            "large");
 
330
          if (count < 0)
 
331
            error_with_id ("Octave:tisean", "The variable 'count' is too "
 
332
                                            "small");
 
333
        }
 
334
 
 
335
      octave_idx_type LENGTH    = input.rows ();
 
336
      octave_idx_type DIMENSION = input.columns ();
 
337
      octave_idx_type alldim    = DIMENSION*EMBED;
 
338
 
 
339
      // The following Matrix need to passed when
 
340
      // performing an interrupted execution
 
341
      Matrix averr_matrix;
 
342
      if (count == 0)
 
343
        averr_matrix.resize (DIMENSION,1);
 
344
      else if (nargin == 13)
 
345
        averr_matrix = args(9).matrix_value ();
 
346
      else
 
347
        error_with_id ("Octave:tisean", "Wrong number of arguments (%d) "
 
348
                                        "for current count = %d",
 
349
                       nargin, count);
 
350
      double *averr = averr_matrix.fortran_vec ();
 
351
 
 
352
      Matrix delta_matrix;
 
353
      if (count == 0)
 
354
        delta_matrix.resize (alldim, alldim);
 
355
      else if (nargin == 13)
 
356
        delta_matrix = args(10).matrix_value ();
 
357
      else
 
358
        error_with_id ("Octave:tisean", "Wrong number of arguments (%d) "
 
359
                                        "for current count = %d",
 
360
                       nargin, count);
 
361
      OCTAVE_LOCAL_BUFFER (double *, delta, alldim);
 
362
      for (octave_idx_type i=0;i<alldim;i++)
 
363
        delta[i] = delta_matrix.fortran_vec () + alldim * i;
 
364
 
 
365
      double avneig;
 
366
      if (nargin == 13)
 
367
        avneig = args(11).double_value ();
 
368
      double aveps;
 
369
      if (nargin == 13)
 
370
        aveps  = args(12).double_value ();
 
371
 
 
372
      // Analyze and rescale input
 
373
      OCTAVE_LOCAL_BUFFER (double, interval, DIMENSION);
 
374
      OCTAVE_LOCAL_BUFFER (double, var, DIMENSION);
 
375
      double maxinterval=0.0;
 
376
      for (octave_idx_type i=0;i<DIMENSION;i++) {
 
377
        double min_val;
 
378
        rescale_data(input,i,LENGTH,&min_val,&interval[i]);
 
379
        if (interval[i] > maxinterval) 
 
380
          maxinterval=interval[i];
 
381
        double av;
 
382
        variance(input.column(i),LENGTH,&av,&var[i]);
 
383
      }
 
384
      
 
385
      // Series is a pointer to data stored in input
 
386
      // so input(i,j) == series[j][i]
 
387
      // This is done for optimization purposes
 
388
      OCTAVE_LOCAL_BUFFER (const double *, series, DIMENSION);
 
389
      for (octave_idx_type j = 0; j < DIMENSION; j++)
 
390
        {
 
391
          const double *ptr = input.fortran_vec ();
 
392
          series[j] = ptr + LENGTH * j;
 
393
        }
 
394
 
 
395
      // Allocate memory
 
396
      OCTAVE_LOCAL_BUFFER (octave_idx_type *, box, BOX);
 
397
      OCTAVE_LOCAL_BUFFER (octave_idx_type, box_data, BOX * BOX);
 
398
      for (octave_idx_type i=0;i<BOX;i++)
 
399
        box[i]=box_data + BOX * i;
 
400
 
 
401
      OCTAVE_LOCAL_BUFFER (double, factor, alldim);
 
402
      OCTAVE_LOCAL_BUFFER (double, lfactor, alldim);
 
403
 
 
404
      OCTAVE_LOCAL_BUFFER (double *, dynamics, DIMENSION);
 
405
      OCTAVE_LOCAL_BUFFER (double, dynamics_data, DIMENSION * alldim);
 
406
      for (octave_idx_type i=0;i<DIMENSION;i++)
 
407
        dynamics[i]=dynamics_data + alldim * i;
 
408
 
 
409
    // old  indexes=make_multi_index(dim,embed,DELAY);
 
410
      OCTAVE_LOCAL_BUFFER (octave_idx_type *, indexes, 2);
 
411
      OCTAVE_LOCAL_BUFFER (octave_idx_type, indexes_data, 2*DIMENSION * EMBED);
 
412
      indexes[0] = indexes_data;
 
413
      indexes[1] = indexes_data + (DIMENSION * EMBED);
 
414
 
 
415
      for (octave_idx_type i=0;i<DIMENSION * EMBED;i++) 
 
416
        {
 
417
          indexes[0][i]=i%DIMENSION;
 
418
          indexes[1][i]=(i/DIMENSION)*DELAY;
 
419
        }
 
420
     // end old indexes = make_multi_index();
 
421
 
 
422
      if (!error_state)
 
423
        {
 
424
 
 
425
          // Promote warnings connected with singular matrixes to errors
 
426
          set_warning_state ("Octave:nearly-singular-matrix","error");
 
427
          set_warning_state ("Octave:singular-matrix","error");
 
428
 
 
429
          // Prepare data for running first time
 
430
          if (count == 0)
 
431
            {
 
432
              averr_matrix.fill (0.0);
 
433
 
 
434
              if (epsset)
 
435
                epsmin /= maxinterval;
 
436
 
 
437
              // old rnd_init(0x098342L);
 
438
              TISEAN_rand generator (0x098342L);
 
439
              for (octave_idx_type i=0;i<10000;i++)
 
440
                generator.rnd_long();
 
441
              for (octave_idx_type i=0;i<alldim;i++) {
 
442
                factor[i]=0.0;
 
443
                for (octave_idx_type j=0;j<alldim;j++)
 
444
                  delta[i][j] = (double)generator.rnd_long()
 
445
                                / (double)std::numeric_limits<octave_idx_type>
 
446
                                          ::max ();
 
447
              }
 
448
              gram_schmidt(alldim, delta,lfactor);
 
449
 
 
450
              avneig = 0.0;
 
451
              aveps  = 0.0;
 
452
            }
 
453
 
 
454
          // Create output
 
455
          Matrix lyap_exp (1, 1 + alldim);
 
456
          time_t lasttime;
 
457
          time(&lasttime);
 
458
          bool pause_calc = false;
 
459
          for (octave_idx_type i = count + (EMBED-1) * DELAY;
 
460
               i < ITERATIONS && !pause_calc; i++)
 
461
            {
 
462
              count++;
 
463
              make_dynamics(series, box, indexes, epsmin, epsset, EPSSTEP,
 
464
                            EMBED, MINNEIGHBORS, LENGTH, DIMENSION, count,
 
465
                            avneig, aveps, dynamics, averr, i);
 
466
              // If there was an error
 
467
              // (matrix singularity or not enough neighbors)
 
468
              // No sense continuing
 
469
              if (error_state)
 
470
                {
 
471
                  return retval;
 
472
                }
 
473
              make_iteration(DIMENSION, alldim, dynamics, delta);
 
474
              gram_schmidt(alldim, delta,lfactor);
 
475
              for (octave_idx_type j=0;j<alldim;j++) {
 
476
                factor[j] += log(lfactor[j])/(double)DELAY;
 
477
              }
 
478
 
 
479
              if (((time(NULL)-lasttime) > OUT) || (i == (ITERATIONS-1))
 
480
                  || (count % it_pause == 0) )
 
481
                {
 
482
                  time(&lasttime);
 
483
 
 
484
                  // Create spectrum output
 
485
                  // old fprintf(stdout,"%ld ",count);
 
486
                  lyap_exp(0,0) = count;
 
487
                  for (octave_idx_type j=0;j<alldim;j++)
 
488
                    {
 
489
                      // old fprintf(stdout,"%e ",factor[j]/count);
 
490
                      lyap_exp(0, 1+j) = factor[j]/count;
 
491
                    }
 
492
 
 
493
                  pause_calc = true;
 
494
                }
 
495
            }
 
496
 
 
497
          // Create pause output
 
498
          if (count < (ITERATIONS - (EMBED-1)*DELAY))
 
499
            {
 
500
              octave_scalar_map pause_vars;
 
501
              pause_vars.setfield ("averr", averr_matrix);
 
502
              pause_vars.setfield ("delta", delta_matrix);
 
503
              pause_vars.setfield ("count", count);
 
504
              pause_vars.setfield ("avneig", avneig);
 
505
              pause_vars.setfield ("aveps", aveps);
 
506
              pause_vars.setfield ("epsmin", epsmin);
 
507
              retval(0) = lyap_exp; 
 
508
              retval(1) = pause_vars;
 
509
            }
 
510
 
 
511
          // Create final output
 
512
          if (count == (ITERATIONS - (EMBED-1)*DELAY))
 
513
            {
 
514
              double dim=0.0;
 
515
              octave_idx_type i;
 
516
              for (i=0;i<alldim;i++) {
 
517
                dim += factor[i];
 
518
                if (dim < 0.0)
 
519
                  break;
 
520
              }
 
521
              if (i < alldim)
 
522
                dim=i+(dim-factor[i])/fabs(factor[i]);
 
523
              else
 
524
                dim=alldim;
 
525
 
 
526
              // Create output pars
 
527
              octave_scalar_map pars;
 
528
 
 
529
              // Create rel_err
 
530
              Matrix rel_err (1, DIMENSION);
 
531
              // old fprintf(stdout,"#Average relative forecast errors:= ");
 
532
              for (octave_idx_type i=0;i<DIMENSION;i++)
 
533
                {
 
534
                  // old fprintf(stdout,"%e ",sqrt(averr[i]/count)/var[i]);
 
535
                  rel_err(0,i) = sqrt(averr[i]/count)/var[i];
 
536
                }
 
537
              pars.setfield ("rel_err", rel_err);
 
538
 
 
539
              // Create abs_err
 
540
              Matrix abs_err (1, DIMENSION);
 
541
              // old fprintf(stdout,"#Average absolute forecast errors:= ");
 
542
              for (octave_idx_type i=0;i<DIMENSION;i++)
 
543
                {
 
544
                  // old fprintf(stdout,"%e ",sqrt(averr[i]/count)
 
545
                  //                          *interval[i]);
 
546
                  abs_err(0,i) = sqrt(averr[i]/count)*interval[i];
 
547
                }
 
548
              pars.setfield ("abs_err", abs_err);
 
549
 
 
550
              // old fprintf(stdout,"#Average Neighborhood Size= %e\n",
 
551
              //             aveps*maxinterval/count);
 
552
              pars.setfield ("nsize", aveps*maxinterval/count);
 
553
 
 
554
              // old fprintf(stdout,"#Average num. of neighbors= %e\n",
 
555
              //             avneig/count);
 
556
              pars.setfield ("nno", avneig/count);
 
557
 
 
558
              // old fprintf(stdout,"#estimated KY-Dimension= %f\n",dim);
 
559
              pars.setfield ("ky_dim", dim);
 
560
 
 
561
              // Assign output
 
562
              retval(0) = lyap_exp;
 
563
              retval(1) = pars;
 
564
            }
 
565
        }
 
566
    }
 
567
  return retval;
 
568
}