← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/__lzo_test__.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/__lzo_test__.cc
 
1
/*
 
2
 *   This file is part of TISEAN
 
3
 *
 
4
 *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber,
 
5
 *                           Piotr Held
 
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 lzo-test 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 <cmath>
 
32
 
 
33
#include "routines_c/tsa.h"
 
34
 
 
35
/*number of boxes for the neighbor search algorithm*/
 
36
#define NMAX 512
 
37
 
 
38
void sort(const Matrix &series, octave_idx_type *found, 
 
39
          double *abstand, octave_idx_type embed, octave_idx_type DELAY, 
 
40
          octave_idx_type MINN, octave_idx_type nfound,
 
41
          double **hser)
 
42
{
 
43
 
 
44
  octave_idx_type hdim = (embed-1) * DELAY;
 
45
  octave_idx_type dim = series.columns ();
 
46
 
 
47
  for (octave_idx_type i=0;i<nfound;i++) {
 
48
    octave_idx_type hf=found[i];
 
49
    abstand[i]=0.0;
 
50
    for (octave_idx_type j=0;j<dim;j++) {
 
51
      for (octave_idx_type k=0;k<=hdim;k += DELAY) {
 
52
      double dx=fabs(series(hf-k,j)-hser[hdim-k][j]);
 
53
      if (dx > abstand[i]) abstand[i]=dx;
 
54
      }
 
55
    }
 
56
  }
 
57
 
 
58
  for (octave_idx_type  i=0;i<MINN;i++)
 
59
    for (octave_idx_type  j=i+1;j<nfound;j++)
 
60
      if (abstand[j]<abstand[i]) {
 
61
      double dswap=abstand[i];
 
62
      abstand[i]=abstand[j];
 
63
      abstand[j]=dswap;
 
64
      octave_idx_type iswap=found[i];
 
65
      found[i]=found[j];
 
66
      found[j]=iswap;
 
67
      }
 
68
}
 
69
 
 
70
void make_fit(const Matrix &series, octave_idx_type dim,
 
71
              octave_idx_type act, octave_idx_type number,
 
72
              octave_idx_type istep, octave_idx_type *found,
 
73
              Matrix &error_array, Matrix &diffs)
 
74
{
 
75
  double casted;
 
76
  const double *help;
 
77
  octave_idx_type h;
 
78
  h=istep-1;
 
79
  octave_idx_type len = series.rows ();
 
80
  for (octave_idx_type j=0;j<dim;j++) 
 
81
    {
 
82
      casted=0.0;
 
83
  //  old  help=series[j]+istep;
 
84
      help=series.fortran_vec()+j*len+istep;
 
85
      for (octave_idx_type i=0;i<number;i++)
 
86
        casted += help[found[i]];
 
87
      casted /= (double)number;
 
88
      diffs(act,j) = casted-help[act];
 
89
      error_array(h,j)  += sqr(casted-help[act]);
 
90
    }
 
91
}
 
92
 
 
93
DEFUN_DLD (__lzo_test__, args, nargout, HELPTEXT)
 
94
{
 
95
 
 
96
  int nargin = args.length ();
 
97
  octave_value_list retval;
 
98
 
 
99
  if ((nargin != 13) || (nargout > 2))
 
100
  {
 
101
    print_usage ();
 
102
  }
 
103
  else
 
104
    {
 
105
 
 
106
      // Assign input
 
107
      Matrix input            = args(0).matrix_value ();
 
108
      octave_idx_type embed   = args(1).idx_type_value ();
 
109
      octave_idx_type DELAY   = args(2).idx_type_value ();
 
110
      octave_idx_type CLENGTH = args(3).idx_type_value ();
 
111
      bool clengthset         = args(4).bool_value ();
 
112
      octave_idx_type refstep = args(5).idx_type_value ();
 
113
      octave_idx_type MINN    = args(6).idx_type_value ();
 
114
      double EPS0             = args(7).double_value ();
 
115
      bool epsset             = args(8).bool_value ();
 
116
      double EPSF             = args(9).double_value ();
 
117
      octave_idx_type STEP    = args(10).idx_type_value ();
 
118
      octave_idx_type causal  = args(11).idx_type_value ();
 
119
      bool setsort            = args(12).bool_value ();
 
120
 
 
121
      octave_idx_type LENGTH  = input.rows ();
 
122
      octave_idx_type dim     = input.columns ();
 
123
 
 
124
      // Allocate memory and analyze input
 
125
      OCTAVE_LOCAL_BUFFER(double*, hser, dim);
 
126
      OCTAVE_LOCAL_BUFFER(double, av, dim);
 
127
      OCTAVE_LOCAL_BUFFER(double, rms, dim);
 
128
      OCTAVE_LOCAL_BUFFER(double, hinter, dim);
 
129
 
 
130
      double mind;
 
131
      double interval=0.0;
 
132
 
 
133
      for (octave_idx_type i=0;i<dim;i++) {
 
134
        rescale_data(input,i,LENGTH,&mind,&hinter[i]);
 
135
        variance(input.column(i),LENGTH,&av[i],&rms[i]);
 
136
        interval += hinter[i];
 
137
      }
 
138
      interval /= (double)dim;
 
139
 
 
140
      OCTAVE_LOCAL_BUFFER (long, list, LENGTH);
 
141
      OCTAVE_LOCAL_BUFFER (octave_idx_type, found, LENGTH);
 
142
      OCTAVE_LOCAL_BUFFER (unsigned long, hfound, LENGTH);
 
143
      OCTAVE_LOCAL_BUFFER (bool, done, LENGTH);
 
144
      OCTAVE_LOCAL_BUFFER (double, abstand, LENGTH);
 
145
 
 
146
      Matrix error_array (STEP, dim);
 
147
      Matrix diffs (LENGTH, dim);
 
148
      error_array.fill (0.0);
 
149
 
 
150
      MArray<octave_idx_type> box (dim_vector(NMAX, NMAX));
 
151
 
 
152
      // Compute forecast error
 
153
      if ( ! error_state)
 
154
        {
 
155
          for (octave_idx_type i=0;i<LENGTH;i++)
 
156
            done[i]=0;
 
157
 
 
158
          bool alldone = false;
 
159
          if (epsset)
 
160
            EPS0 /= interval;
 
161
 
 
162
          double epsilon = EPS0 / EPSF;
 
163
 
 
164
          if (!clengthset)
 
165
            CLENGTH=LENGTH;
 
166
          octave_idx_type clength = ((CLENGTH*refstep+STEP) <= LENGTH) 
 
167
                                    ? CLENGTH : (LENGTH-STEP)/refstep;
 
168
 
 
169
          // Compute estimates
 
170
          octave_idx_type actfound;
 
171
          octave_idx_type hi;
 
172
          while (!alldone) {
 
173
            alldone=1;
 
174
            epsilon*=EPSF;
 
175
            make_multi_box(input,box,list,LENGTH-(long)STEP,NMAX,dim,
 
176
                           embed,DELAY,epsilon);
 
177
            for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++)
 
178
              if (!done[i]) {
 
179
              hi=i*refstep;
 
180
              for (octave_idx_type j=0;j<dim;j++)
 
181
                {
 
182
//              old  hser[j]=series[j]+hi;
 
183
                  hser[j] = input.fortran_vec() + j * LENGTH + hi;
 
184
                }
 
185
              actfound=find_multi_neighbors(input,box,list,hser,NMAX,
 
186
                                            dim,embed,DELAY,epsilon,hfound);
 
187
              actfound=exclude_interval(actfound,hi-(long)causal+1,
 
188
                                  hi+causal+(embed-1)*DELAY-1,hfound,found);
 
189
              if (actfound >= MINN)
 
190
                {
 
191
                  if (setsort)
 
192
                    {
 
193
                      sort(input, found, abstand, embed, DELAY, MINN,
 
194
                           actfound, hser);
 
195
                      actfound=MINN;
 
196
                    }
 
197
                  for (octave_idx_type j=1;j<=STEP;j++) {
 
198
                    make_fit(input,dim,hi,actfound,j,found,error_array,diffs);
 
199
                  }
 
200
                  done[i]=1;
 
201
                }
 
202
              alldone &= done[i];
 
203
              }
 
204
          }
 
205
 
 
206
          // Create relative forecast error output
 
207
          Matrix rel_forecast_err (STEP, dim + 1);
 
208
          for (octave_idx_type i=0;i<STEP;i++) 
 
209
            {
 
210
              rel_forecast_err(i,0) = i + 1;
 
211
              for (octave_idx_type j=0;j<dim;j++) 
 
212
                {
 
213
//            old  fprintf(stdout,"%e ",
 
214
//                    sqrt(error[j][i]/(clength-(embed-1)*DELAY))/rms[j]);
 
215
                  rel_forecast_err(i,j+1) = sqrt(error_array(i,j)
 
216
                                          /(clength-(embed-1)*DELAY))/rms[j];
 
217
                }
 
218
            }
 
219
 
 
220
          // Create individual forecast error output
 
221
          Matrix ind_forecast_err (1,1);
 
222
          if (nargout > 1)
 
223
            {
 
224
              ind_forecast_err.resize(clength - (embed-1)*DELAY, dim);
 
225
              for (octave_idx_type i=(embed-1)*DELAY;i<clength;i++) 
 
226
                {
 
227
                  hi=i*refstep;
 
228
                  for (octave_idx_type j=0;j<dim;j++)
 
229
                    {
 
230
//                   old fprintf(stdout,"%e ",diffs[j][hi]*hinter[j]);
 
231
                      ind_forecast_err(i-(embed-1)*DELAY,j) = \
 
232
                      diffs(hi,j)*hinter[j];
 
233
                    }
 
234
                }
 
235
            }
 
236
 
 
237
          retval(0) = rel_forecast_err;
 
238
          retval(1) = ind_forecast_err;
 
239
        }
 
240
    }
 
241
  return retval;
 
242
}