← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/__ghkss__.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/__ghkss__.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 ghkss 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
 
 
32
#include "routines_c/tsa.h"
 
33
 
 
34
#define BOX 1024
 
35
 
 
36
int eps_set=0,euclidean=0,dimset=0;
 
37
 
 
38
octave_idx_type dim,qdim=2,minn=50,iterations=1,embed=5,length;
 
39
octave_idx_type comp=1;
 
40
unsigned int delay=1;
 
41
unsigned int verbosity=0;
 
42
double mineps,epsfac;
 
43
 
 
44
char resize_eps;
 
45
 
 
46
double *d_min,*d_max,d_max_max;
 
47
double **delta,**corr;
 
48
double *metric,trace;
 
49
long **box,*list;
 
50
unsigned long *flist;
 
51
int emb_offset;
 
52
unsigned int ibox=BOX-1;
 
53
unsigned int *index_comp,*index_embed;
 
54
 
 
55
/*these are global to save time*/
 
56
int *sorted;
 
57
double *av,**mat,*matarray,*eig;
 
58
 
 
59
 
 
60
void sort(double *x,int *n)
 
61
{
 
62
  octave_idx_type i,j;
 
63
  long iswap;
 
64
  double dswap;
 
65
  
 
66
  for (i=0;i<dim;i++)
 
67
    n[i]=i;
 
68
  
 
69
  for (i=0;i<dim-1;i++)
 
70
    for (j=i+1;j<dim;j++)
 
71
      if (x[j] > x[i]) {
 
72
      dswap=x[i];
 
73
      x[i]=x[j];
 
74
      x[j]=dswap;
 
75
      iswap=n[i];
 
76
      n[i]=n[j];
 
77
      n[j]=iswap;
 
78
      }
 
79
}
 
80
 
 
81
void mmb(const Matrix &series,double eps)
 
82
{
 
83
  octave_idx_type i,x,y;
 
84
  double ieps=1.0/eps;
 
85
 
 
86
  for (x=0;x<BOX;x++)
 
87
    for (y=0;y<BOX;y++)
 
88
      box[x][y] = -1;
 
89
   
 
90
  for (i=emb_offset;i<length;i++) {
 
91
    x=(int)(series(i,0)*ieps)&ibox;
 
92
    y=(int)(series(i-emb_offset,comp-1)*ieps)&ibox;
 
93
    list[i]=box[x][y];
 
94
    box[x][y]=i;
 
95
  }
 
96
}
 
97
 
 
98
unsigned long fmn(const Matrix &series, long which,double eps)
 
99
{
 
100
  unsigned long nf=0;
 
101
  octave_idx_type i,i1,i2,j,j1,k,li;
 
102
  long k1;
 
103
  long element;
 
104
  double dx=0.0;
 
105
  
 
106
  i=(int)(series(which,0)/eps)&ibox;
 
107
  j=(int)(series(which-emb_offset,comp-1)/eps)&ibox;
 
108
  
 
109
  for (i1=i-1;i1<=i+1;i1++) {
 
110
    i2=i1&ibox;
 
111
    for (j1=j-1;j1<=j+1;j1++) {
 
112
      element=box[i2][j1&ibox];
 
113
      while (element != -1) {
 
114
      for (k=0;k<embed;k++) {
 
115
        k1= -k*(int)delay;
 
116
        for (li=0;li<comp;li++) {
 
117
          dx=fabs(series(which+k1,li)-series(element+k1,li));
 
118
          if (dx > eps)
 
119
            break;
 
120
        }
 
121
        if (dx > eps)
 
122
          break;
 
123
      }
 
124
      if (dx <= eps)
 
125
        flist[nf++]=element;
 
126
      element=list[element];
 
127
      }
 
128
    }
 
129
  }
 
130
  return nf;
 
131
}
 
132
 
 
133
void make_correction(const Matrix &series, unsigned long n,octave_idx_type nf)
 
134
{
 
135
  octave_idx_type i,i1,i2,j,j1,j2,k,k1,k2,hs;
 
136
  double help;
 
137
  
 
138
  for (i=0;i<dim;i++) {
 
139
    i1=index_comp[i];
 
140
    i2=index_embed[i];
 
141
    help=0.0;
 
142
    for (j=0;j<nf;j++)
 
143
      help += series(flist[j]-i2,i1);
 
144
    av[i]=help/nf;
 
145
  }
 
146
 
 
147
  for (i=0;i<dim;i++) {
 
148
    i1=index_comp[i];
 
149
    i2=index_embed[i];
 
150
    for (j=i;j<dim;j++) {
 
151
      help=0.0;
 
152
      j1=index_comp[j];
 
153
      j2=index_embed[j];
 
154
      for (k=0;k<nf;k++) {
 
155
      hs=flist[k];
 
156
      help += series(hs-i2,i1)*series(hs-j2,j1);
 
157
      }
 
158
      mat[i][j]=(help/nf-av[i]*av[j])*metric[i]*metric[j];
 
159
      mat[j][i]=mat[i][j];
 
160
    }
 
161
  }
 
162
 
 
163
  eigen(mat,dim,eig);
 
164
  sort(eig,sorted);
 
165
 
 
166
  for (i=0;i<dim;i++) {
 
167
    help=0.0;
 
168
    for (j=qdim;j<dim;j++) {
 
169
      hs=sorted[j];
 
170
      for (k=0;k<dim;k++) {
 
171
      k1=index_comp[k];
 
172
      k2=index_embed[k];
 
173
      help += (series(n-k2,k1)-av[k])*mat[k][hs]*mat[i][hs]*metric[k];
 
174
      }
 
175
    }
 
176
    corr[n][i]=help/metric[i];
 
177
  }
 
178
}
 
179
 
 
180
void handle_trend(unsigned long n,octave_idx_type nf)
 
181
{
 
182
  octave_idx_type i,i1,i2,j;
 
183
  double help;
 
184
  
 
185
  for (i=0;i<dim;i++) {
 
186
    help=0.0;
 
187
    for (j=0;j<nf;j++)
 
188
      help += corr[flist[j]][i];
 
189
    av[i]=help/nf;
 
190
  }
 
191
 
 
192
  for (i=0;i<dim;i++) {
 
193
    i1=index_comp[i];
 
194
    i2=index_embed[i];
 
195
    delta[i1][n-i2] += (corr[n][i]-av[i])/(trace*metric[i]);
 
196
  }
 
197
}
 
198
 
 
199
void set_correction(Matrix &series)
 
200
{
 
201
  octave_idx_type i,j;
 
202
  double help;
 
203
 
 
204
  OCTAVE_LOCAL_BUFFER (double, hav, comp);
 
205
  OCTAVE_LOCAL_BUFFER (double, hsigma, comp);
 
206
  if ( ! error_state)
 
207
    {
 
208
      for (j=0;j<comp;j++)
 
209
        hav[j]=hsigma[j]=0.0;
 
210
 
 
211
      for (i=0;i<length;i++)
 
212
        for (j=0;j<comp;j++) {
 
213
          hav[j] += (help=delta[j][i]);
 
214
          hsigma[j] += help*help;
 
215
        }
 
216
 
 
217
      for (j=0;j<comp;j++) {
 
218
        hav[j] /= length;
 
219
        hsigma[j]=sqrt(fabs(hsigma[j]/length-hav[j]*hav[j]));
 
220
      }
 
221
      if (verbosity) {
 
222
        for (i=0;i<comp;i++) {
 
223
          octave_stdout << "Average shift of component " << i+1 << " = "
 
224
                        << hav[i]*d_max[i] << "\n";
 
225
          octave_stdout << "Average rms correction of comp. " << i+1 << " = "
 
226
                        << hsigma[i]*d_max[i] << "\n\n";
 
227
        }
 
228
      }
 
229
      for (i=0;i<length;i++)
 
230
        for (j=0;j<comp;j++)
 
231
          series(i,j) -= delta[j][i];
 
232
 
 
233
      if (resize_eps) {
 
234
        mineps /= epsfac;
 
235
        if (verbosity)
 
236
          octave_stdout << "Reset minimal neighbourhood size to "
 
237
                        << mineps*d_max_max << "\n";
 
238
      }
 
239
 
 
240
      resize_eps=0;
 
241
    }
 
242
}
 
243
 
 
244
DEFUN_DLD (__ghkss__, args, , HELPTEXT)
 
245
{
 
246
 
 
247
  int epscount,*ok;
 
248
  octave_idx_type iter,nfound,n;
 
249
  octave_idx_type i,j;
 
250
  char all_done;
 
251
  unsigned long allfound;
 
252
  double epsilon;
 
253
  double **hser;
 
254
 
 
255
  int nargin = args.length ();
 
256
  octave_value_list retval;
 
257
 
 
258
  if (nargin != 11)
 
259
  {
 
260
    print_usage ();
 
261
  }
 
262
  else
 
263
    {
 
264
 
 
265
      // Assign input
 
266
      Matrix input = args(0).matrix_value();
 
267
      embed        = args(1).int_value();
 
268
      comp         = args(2).int_value();
 
269
      delay        = args(3).int_value();
 
270
      qdim         = args(4).int_value();
 
271
      minn         = args(5).int_value();
 
272
      mineps       = args(6).double_value();
 
273
      eps_set      = args(7).int_value();
 
274
      iterations   = args(8).int_value();
 
275
      euclidean    = args(9).int_value();
 
276
      verbosity    = args(10).int_value();
 
277
 
 
278
      dim=comp*embed;
 
279
      emb_offset=(embed-1)*delay;
 
280
      length = input.rows();
 
281
 
 
282
      // Prepare for noise reduction and allocate memory
 
283
      check_alloc(d_min=(double*)malloc(sizeof(double)*comp));
 
284
      check_alloc(d_max=(double*)malloc(sizeof(double)*comp));
 
285
      d_max_max=0.0;
 
286
      for (i=0;i<comp;i++) {
 
287
        rescale_data(input,i,length,&d_min[i],&d_max[i]);
 
288
        if (d_max[i] > d_max_max)
 
289
          d_max_max=d_max[i];
 
290
      }
 
291
 
 
292
      if (!eps_set)
 
293
        mineps=1./1000.;
 
294
      else
 
295
        mineps /= d_max_max;
 
296
      epsfac=sqrt(2.0);
 
297
 
 
298
      check_alloc(box=(long**)malloc(sizeof(long*)*BOX));
 
299
      for (i=0;i<BOX;i++)
 
300
        check_alloc(box[i]=(long*)malloc(sizeof(long)*BOX));
 
301
 
 
302
      check_alloc(list=(long*)malloc(sizeof(long)*length));
 
303
      check_alloc(flist=(unsigned long*)malloc(sizeof(long)*length));
 
304
 
 
305
      check_alloc(metric=(double*)malloc(sizeof(double)*dim));
 
306
      trace=0.0;
 
307
      if (euclidean) {
 
308
        for (i=0;i<dim;i++) {
 
309
          metric[i]=1.0;
 
310
          trace += 1./metric[i];
 
311
        }
 
312
      }
 
313
      else {
 
314
        for (i=0;i<dim;i++) {
 
315
          if ((i >= comp) && (i < ((long)dim-(long)comp))) 
 
316
          metric[i]=1.0;
 
317
          else 
 
318
          metric[i]=1.0e3;
 
319
          trace += 1./metric[i];
 
320
        }
 
321
      }
 
322
 
 
323
      check_alloc(corr=(double**)malloc(sizeof(double*)*length));
 
324
      for (i=0;i<length;i++)
 
325
        check_alloc(corr[i]=(double*)malloc(sizeof(double)*dim));
 
326
      check_alloc(ok=(int*)malloc(sizeof(int)*length));
 
327
      check_alloc(delta=(double**)malloc(sizeof(double*)*comp));
 
328
      for (i=0;i<comp;i++)
 
329
        check_alloc(delta[i]=(double*)malloc(sizeof(double)*length));
 
330
      check_alloc(index_comp=(unsigned int*)malloc(sizeof(int)*dim));
 
331
      check_alloc(index_embed=(unsigned int*)malloc(sizeof(int)*dim));
 
332
      check_alloc(av=(double*)malloc(sizeof(double)*dim));
 
333
      check_alloc(sorted=(int*)malloc(sizeof(int)*dim));
 
334
      check_alloc(eig=(double*)malloc(sizeof(double)*dim));
 
335
      check_alloc(matarray=(double*)malloc(sizeof(double)*dim*dim));
 
336
      check_alloc(mat=(double**)malloc(sizeof(double*)*dim));
 
337
      for (i=0;i<dim;i++)
 
338
        mat[i]=(double*)(matarray+dim*i);
 
339
      check_alloc(hser=(double**)malloc(sizeof(double*)*comp));
 
340
 
 
341
      if ( ! error_state)
 
342
        {
 
343
          // Create output matrix
 
344
          Matrix output (length, comp);
 
345
 
 
346
          for (i=0;i<dim;i++) {
 
347
            index_comp[i]=i%comp;
 
348
            index_embed[i]=(i/comp)*delay;
 
349
          }
 
350
 
 
351
          // Calculate the noise reduction
 
352
          resize_eps=0;
 
353
          for (iter=1;iter<=iterations;iter++) 
 
354
            {
 
355
              for (i=0;i<length;i++) {
 
356
                ok[i]=0;
 
357
                for (j=0;j<dim;j++)
 
358
                corr[i][j]=0.0;
 
359
                for (j=0;j<comp;j++)
 
360
                delta[j][i]=0.0;
 
361
              }
 
362
              epsilon=mineps;
 
363
              all_done=0;
 
364
              epscount=1;
 
365
              allfound=0;
 
366
              if (verbosity)
 
367
                octave_stdout << "Starting iteration " << iter << "\n";
 
368
              while(!all_done) {
 
369
                mmb(input, epsilon);
 
370
                all_done=1;
 
371
                for (n=emb_offset;n<length;n++)
 
372
                if (!ok[n]) {
 
373
                  nfound=fmn(input,n,epsilon);
 
374
                  if (nfound >= minn) {
 
375
                    make_correction(input,n,nfound);
 
376
                    ok[n]=epscount;
 
377
                    if (epscount == 1)
 
378
                      resize_eps=1;
 
379
                    allfound++;
 
380
                  }
 
381
                  else
 
382
                    all_done=0;
 
383
                }
 
384
                if (verbosity)
 
385
                  octave_stdout << "Corrected " << allfound << " points with epsilon= "
 
386
                                << epsilon*d_max_max << "\n";
 
387
                if (std::isinf(epsilon*d_max_max))
 
388
                  {
 
389
                    error_with_id ("Octave:tisean", "cannot reduce noise on input data");
 
390
                    return retval;
 
391
                  }
 
392
                epsilon *= epsfac;
 
393
                epscount++;
 
394
              }
 
395
              if (verbosity)
 
396
                octave_stdout << "Start evaluating the trend\n";
 
397
 
 
398
              epsilon=mineps;
 
399
              allfound=0;
 
400
              for (i=1;i<epscount;i++) {
 
401
                mmb(input,epsilon);
 
402
                for (n=emb_offset;n<length;n++)
 
403
                if (ok[n] == i) {
 
404
                  nfound=fmn(input,n,epsilon);
 
405
                  handle_trend(n,nfound);
 
406
                  allfound++;
 
407
                }
 
408
                if (verbosity)
 
409
                octave_stdout << "Trend subtracted for " << allfound << " points with epsilon= " 
 
410
                              << epsilon*d_max_max << "\n";
 
411
                epsilon *= epsfac;
 
412
              }
 
413
              set_correction(input);
 
414
 
 
415
              if (iter == iterations)
 
416
                for (i=0;i<length;i++) 
 
417
                  {
 
418
                    for (j=0;j<comp;j++) 
 
419
                      {
 
420
                      // old  fprintf(file,"%e ",series[j][i]*d_max[j]+d_min[j]);
 
421
                        output(i,j) = input(i,j)*d_max[j]+d_min[j];
 
422
                      }
 
423
                  }
 
424
            }
 
425
          retval(0) = output;
 
426
        }
 
427
    // Deallocate of all the memory
 
428
      delete[] d_min;
 
429
      delete[] d_max;
 
430
      for (i=0;i<BOX;i++)
 
431
        delete[] box[i];
 
432
      delete[] box;
 
433
      delete[] list;
 
434
      delete[] flist;
 
435
      delete[] metric;
 
436
      for (i=0;i<length;i++)
 
437
        delete[] corr[i];
 
438
      delete[] corr;
 
439
      delete[] ok;
 
440
      for (i=0;i<comp;i++)
 
441
        delete[] delta[i];
 
442
      delete[] delta;
 
443
      delete[] index_comp;
 
444
      delete[] index_embed;
 
445
      delete[] av;
 
446
      delete[] sorted;
 
447
      delete[] eig;
 
448
      delete[] matarray;
 
449
      delete[] mat;
 
450
      delete[] hser;
 
451
    }
 
452
  return retval;
 
453
}