← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/__polynom__.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/__polynom__.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 polynom 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 "routines_c/tsa.h"
 
32
 
 
33
double polynom(double *series, octave_idx_type DELAY, octave_idx_type N,
 
34
               octave_idx_type act,octave_idx_type dim,long cur,long fac)
 
35
{
 
36
  double ret=1.0;
 
37
 
 
38
  octave_idx_type n  = cur/fac;
 
39
  octave_idx_type hi = act-(dim-1)*DELAY;
 
40
  for (octave_idx_type j=1;j<=n;j++)
 
41
    ret *= series[hi];
 
42
  if (dim > 1) 
 
43
    ret *= polynom(series, DELAY, N,act,dim-1,cur-n*fac,fac/(N+1));
 
44
 
 
45
  return ret;
 
46
}
 
47
 
 
48
octave_idx_type number_pars(octave_idx_type DIM, octave_idx_type ord,
 
49
                            octave_idx_type start)
 
50
{
 
51
 
 
52
  octave_idx_type ret=0;
 
53
 
 
54
  if (ord == 1)
 
55
    for (octave_idx_type i=start;i<=DIM;i++)
 
56
      ret += 1;
 
57
  else
 
58
    for (octave_idx_type i=start;i<=DIM;i++)
 
59
      ret += number_pars(DIM, ord-1,i);
 
60
 
 
61
  return ret;
 
62
}
 
63
 
 
64
void make_coding(std::vector <octave_idx_type> &coding_vec, octave_idx_type N,
 
65
                 octave_idx_type ord, octave_idx_type d,
 
66
                 octave_idx_type fac, octave_idx_type cur)
 
67
{
 
68
  if (d == -1)
 
69
    coding_vec.push_back (cur);
 
70
  else
 
71
    for (octave_idx_type j=0;j<=ord;j++)
 
72
      make_coding(coding_vec, N, ord-j,d-1,fac*(N+1),cur+j*fac);
 
73
}
 
74
 
 
75
void make_fit(double *series, octave_idx_type *coding, double *results,
 
76
              octave_idx_type INSAMPLE, octave_idx_type N,
 
77
              octave_idx_type DIM, octave_idx_type DELAY,
 
78
              long maxencode, octave_idx_type pars)
 
79
{
 
80
 
 
81
  Matrix mat_b (pars,1);
 
82
  double *b = mat_b.fortran_vec ();
 
83
 
 
84
  Matrix mat (pars,pars);
 
85
  OCTAVE_LOCAL_BUFFER (double *, mat_arr, pars);
 
86
 
 
87
  for (octave_idx_type i=0;i<pars;i++)
 
88
    mat_arr[i]=mat.fortran_vec () + i * pars;
 
89
 
 
90
  for (octave_idx_type i=0;i<pars;i++) {
 
91
    b[i]=0.0;
 
92
    for (octave_idx_type j=0;j<pars;j++)
 
93
      mat_arr[i][j]=0.0;
 
94
  }
 
95
 
 
96
  for (octave_idx_type i=0;i<pars;i++)
 
97
    for (octave_idx_type j=i;j<pars;j++)
 
98
      for (octave_idx_type k=(DIM-1)*DELAY;k<INSAMPLE-1;k++)
 
99
        mat_arr[i][j] += polynom(series,DELAY,N,k,DIM,coding[i],maxencode)*
 
100
          polynom(series,DELAY,N,k,DIM,coding[j],maxencode);
 
101
  for (octave_idx_type i=0;i<pars;i++)
 
102
    for (octave_idx_type j=i;j<pars;j++)
 
103
      mat_arr[j][i]=(mat_arr[i][j] /= (INSAMPLE-1-(DIM-1)*DELAY));
 
104
 
 
105
  for (octave_idx_type i=0;i<pars;i++) {
 
106
    for (octave_idx_type j=(DIM-1)*DELAY;j<INSAMPLE-1;j++)
 
107
      b[i] += series[j+1]*polynom(series,DELAY,N,j,DIM,coding[i],maxencode);
 
108
    b[i] /= (INSAMPLE-1-(DIM-1)*DELAY);
 
109
  }
 
110
 
 
111
// old  solvele(mat_arr,b,pars);
 
112
  Matrix solved_vec = mat.solve (mat_b);
 
113
  double *solved_ptr = solved_vec.fortran_vec ();
 
114
 
 
115
  for (octave_idx_type i=0;i<pars;i++)
 
116
    results[i]=solved_ptr[i];
 
117
 
 
118
}
 
119
 
 
120
void decode(octave_idx_type N,octave_idx_type *out,int dim,long cur,long fac)
 
121
{
 
122
 
 
123
  octave_idx_type n=cur/fac;
 
124
  out[dim]=n;
 
125
  if (dim > 0) 
 
126
    decode(N,out,dim-1,cur-(long)n*fac,fac/(N+1));
 
127
}
 
128
 
 
129
double make_error(double *series, octave_idx_type *coding, double *results,
 
130
                  octave_idx_type N, octave_idx_type DIM,
 
131
                  octave_idx_type DELAY, long maxencode, octave_idx_type pars,
 
132
                  octave_idx_type i0, octave_idx_type i1)
 
133
{
 
134
 
 
135
  double err=0.0;
 
136
  for (octave_idx_type j=i0+(DIM-1)*DELAY;j<(long)i1-1;j++) {
 
137
    double h=0.0;
 
138
    for (octave_idx_type k=0;k<pars;k++) 
 
139
      h += results[k]*polynom(series,DELAY,N,j,DIM,coding[k],maxencode);
 
140
    err += (series[j+1]-h)*(series[j+1]-h);
 
141
  }
 
142
  return err /= (double)(i1-i0-(DIM-1)*DELAY);
 
143
}
 
144
 
 
145
void make_cast (NDArray &forecast, double *series, octave_idx_type *coding,
 
146
                double * results, octave_idx_type LENGTH,
 
147
                octave_idx_type CLENGTH, octave_idx_type N,
 
148
                octave_idx_type DIM, octave_idx_type DELAY, long maxencode,
 
149
                octave_idx_type pars,
 
150
                double std_dev)
 
151
{
 
152
 
 
153
  for (octave_idx_type i=0;i<=(DIM-1)*DELAY;i++)
 
154
    series[i]=series[LENGTH-(DIM-1)*DELAY-1+i];
 
155
 
 
156
  octave_idx_type hi=(DIM-1)*DELAY;
 
157
  forecast.resize (dim_vector (CLENGTH,1));
 
158
  for (octave_idx_type i=1;i<=CLENGTH;i++) 
 
159
    {
 
160
      double casted=0.0;
 
161
      for (octave_idx_type k=0;k<pars;k++)
 
162
        casted += results[k]*polynom(series,DELAY,N,(DIM-1)*DELAY,DIM,
 
163
                                     coding[k],maxencode);
 
164
 
 
165
    // old fprintf(fcast,"%e\n",casted*std_dev);
 
166
      forecast(i-1) = casted*std_dev;
 
167
      for (octave_idx_type j=0;j<(DIM-1)*DELAY;j++)
 
168
        series[j]=series[j+1];
 
169
      series[hi]=casted;
 
170
    }
 
171
}
 
172
 
 
173
DEFUN_DLD (__polynom__, args, nargout, HELPTEXT)
 
174
{
 
175
 
 
176
  int nargin = args.length ();
 
177
  octave_value_list retval;
 
178
 
 
179
  if (nargin != 6 || nargout != 5)
 
180
  {
 
181
    print_usage ();
 
182
  }
 
183
  else
 
184
    {
 
185
 
 
186
      // Assign input
 
187
      NDArray input = args(0).array_value ();
 
188
      octave_idx_type DIM      = args(1).idx_type_value ();
 
189
      octave_idx_type DELAY    = args(2).idx_type_value ();
 
190
      octave_idx_type N        = args(3).idx_type_value ();
 
191
      octave_idx_type INSAMPLE = args(4).idx_type_value ();
 
192
      octave_idx_type CLENGTH  = args(5).idx_type_value ();
 
193
 
 
194
      octave_idx_type LENGTH = input.numel ();
 
195
      double *series = input.fortran_vec ();
 
196
 
 
197
      // Analyze inputs
 
198
      double std_dev, av;
 
199
      variance(input,LENGTH,&av,&std_dev);
 
200
 
 
201
      // Rescale input
 
202
      for (octave_idx_type i=0;i<LENGTH;i++)
 
203
        series[i] /= std_dev;
 
204
 
 
205
      // Create help values for the fit
 
206
      long maxencode=1;
 
207
      for (octave_idx_type i=1;i<DIM;i++)
 
208
        maxencode *= (N+1);
 
209
 
 
210
      octave_idx_type pars = 1;
 
211
      for (octave_idx_type i=1;i<=N;i++) {
 
212
        pars += number_pars(DIM, i,1);
 
213
      }
 
214
 
 
215
      // Allocate memory
 
216
      OCTAVE_LOCAL_BUFFER (double, results, pars);
 
217
      OCTAVE_LOCAL_BUFFER (octave_idx_type, opar, DIM);
 
218
 
 
219
      // Create coding
 
220
      std::vector <octave_idx_type> coding_vec;
 
221
      coding_vec.reserve (pars);
 
222
      make_coding(coding_vec,N,N,DIM-1,1,0);
 
223
      octave_idx_type *coding = coding_vec.data();
 
224
 
 
225
      if (! error_state)
 
226
        {
 
227
          // Promote warnings connected with singular matrixes to errors
 
228
          set_warning_state ("Octave:nearly-singular-matrix","error");
 
229
          set_warning_state ("Octave:singular-matrix","error");
 
230
 
 
231
          // Make the fit
 
232
          make_fit (series, coding, results, INSAMPLE, N, DIM, DELAY,
 
233
                    maxencode, pars);
 
234
 
 
235
          // If error encountered during the fit there is no sense to continue
 
236
          if (error_state)
 
237
            {
 
238
              return retval;
 
239
            }
 
240
 
 
241
          // Create outputs
 
242
 
 
243
          // Create output that contains the number of free parameters
 
244
          // old fprintf(file,"#number of free parameters= %d\n\n",pars);
 
245
          octave_idx_type free_par = pars;
 
246
 
 
247
          // Create output that contains the norm used for the fit
 
248
          // old fprintf(file,"#used norm for the fit= %e\n",std_dev);
 
249
          double fit_norm = std_dev;
 
250
 
 
251
          // Create coefficients output
 
252
          Matrix coeffs (pars, DIM + 1);
 
253
          for (octave_idx_type j=0;j<pars;j++)
 
254
            {
 
255
              decode(N,opar,DIM-1,coding[j],maxencode);
 
256
              octave_idx_type sumpar=0;
 
257
              for (octave_idx_type k=0;k<DIM;k++)
 
258
                {
 
259
                  sumpar += opar[k];
 
260
                //  old fprintf(file,"%d ",opar[k]);
 
261
                  coeffs(j, k) = opar[k];
 
262
                }
 
263
          //  old fprintf(file,"%e\n",results[j]
 
264
          //                          /pow(std_dev,(double)(sumpar-1)));
 
265
              coeffs(j,DIM) = results[j]/pow(std_dev,(double)(sumpar-1));
 
266
 
 
267
            }
 
268
 
 
269
          // Create sample error
 
270
          // 1st element of sample_error is the insample error
 
271
          // 2nd element of sample_error is the out of sample error (if exists)
 
272
          NDArray sample_err (dim_vector (1,1));
 
273
 
 
274
    // old in_error = make_error((unsigned long)0,INSAMPLE)
 
275
    //     fprintf(file,"#average insample error= %e\n",
 
276
    //             sqrt(in_error)*std_dev);
 
277
          sample_err(0) = sqrt(make_error(series, coding, results, N, DIM,
 
278
                                          DELAY, maxencode, pars,
 
279
                                          0,INSAMPLE))
 
280
                          * std_dev;
 
281
 
 
282
          if (INSAMPLE < LENGTH) 
 
283
            {
 
284
            // old out_error=make_error(INSAMPLE,LENGTH);
 
285
            //     fprintf(file,"#average out of sample error= %e\n",
 
286
            //             sqrt(out_error)*std_dev);
 
287
              sample_err.resize (dim_vector (2,1));
 
288
              sample_err(1) = sqrt (make_error (series, coding, results, N,
 
289
                                                DIM, DELAY, maxencode, pars,
 
290
                                                INSAMPLE, LENGTH)) 
 
291
                              * std_dev;
 
292
            }
 
293
 
 
294
          // Create forecast
 
295
          NDArray forecast (dim_vector (0,0));
 
296
          if (CLENGTH > 0)
 
297
            make_cast(forecast, series, coding, results, LENGTH, CLENGTH, N,
 
298
                      DIM, DELAY, maxencode, pars, std_dev);
 
299
 
 
300
          // Assign outputs
 
301
          retval(0) = free_par;
 
302
          retval(1) = fit_norm;
 
303
          retval(2) = coeffs;
 
304
          retval(3) = sample_err;
 
305
          retval(4) = forecast;
 
306
        }
 
307
    }
 
308
  return retval;
 
309
}