binaries/src/disembl/Tisean_3.0.1/source_c/pca.c

   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 Last modified: Jul 26, 2004 */
  21 #include <stdio.h>
  22 #include <string.h>
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <limits.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Performs a global PCA"
  29
  30 unsigned long LENGTH=ULONG_MAX,exclude=0;
  31 unsigned int DIM=2,EMB=1,dimemb,LDIM=2,DELAY=1;
  32 unsigned int verbosity=0xff;
  33 char *outfile=NULL,stout=1,dim_set=0;
  34 unsigned int what_to_write=0,write_values=1,write_vectors=0;
  35 unsigned int write_comp=0,write_proj=0;
  36 unsigned int projection_set=0;
  37 char *infile=NULL,dimset=0,*column=NULL;
  38 double **series;
  39
  40 void show_options(char *progname)
  41 {
  42   what_i_do(progname,WID_STR);
  43   fprintf(stderr," Usage: %s [options]\n",progname);
  44   fprintf(stderr," Options:\n");
  45   fprintf(stderr,"Everything not being a valid option will be interpreted"
  46           " as a possible"
  47           " datafile.\nIf no datafile is given stdin is read. Just - also"
  48           " means stdin\n");
  49   fprintf(stderr,"\t-l # of data to use [Default: whole file]\n");
  50   fprintf(stderr,"\t-x # of lines to be ignore [Default: 0]\n");
  51   fprintf(stderr,"\t-c columns to read [Default: 2]\n");
  52   fprintf(stderr,"\t-m columns,embedding dim. to use [Default: 2,1]\n");
  53   fprintf(stderr,"\t-d delay to use [Default: 1]\n");
  54   fprintf(stderr,"\t-q projection dimension [Default: no projection]\n");
  55   fprintf(stderr,"\t-W # what to write: [Default: 0]\n"
  56           "\t\t0 write eigenvalues only\n"
  57           "\t\t1 write eigenvectors\n"
  58           "\t\t2 write (projected) pca components\n"
  59           "\t\t3 write projected data\n");
  60   fprintf(stderr,"\t-o output file name \n\t\t[Default: stdout; -o without "
  61           "value means 'datafile'.pca]\n");
  62   fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
  63           "0='only panic messages'\n\t\t"
  64           "1='+ input/output messages'\n");
  65   fprintf(stderr,"\t-h show these options\n");
  66   exit(0);
  67 }
  68
  69 void scan_options(int n,char **in)
  70 {
  71   char *out;
  72
  73   if ((out=check_option(in,n,'l','u')) != NULL)
  74     sscanf(out,"%lu",&LENGTH);
  75   if ((out=check_option(in,n,'x','u')) != NULL)
  76     sscanf(out,"%lu",&exclude);
  77   if ((out=check_option(in,n,'c','s')) != NULL)
  78     column=out;
  79   if ((out=check_option(in,n,'m','2')) != NULL) {
  80     sscanf(out,"%u,%u",&DIM,&EMB);
  81     dimset=1;
  82   }
  83   if ((out=check_option(in,n,'d','u')) != NULL)
  84     sscanf(out,"%u",&DELAY);
  85   if ((out=check_option(in,n,'q','u')) != NULL) {
  86     sscanf(out,"%u",&LDIM);
  87     projection_set=1;
  88   }
  89   if ((out=check_option(in,n,'V','u')) != NULL)
  90     sscanf(out,"%u",&verbosity);
  91   if ((out=check_option(in,n,'W','u')) != NULL) {
  92     sscanf(out,"%u",&what_to_write);
  93     switch(what_to_write) {
  94     case 0: write_values=1;break;
  95     case 1: write_values=0;write_vectors=1;break;
  96     case 2: write_values=0;write_comp=1;break;
  97     case 3: write_values=0;write_proj=1;break;
  98     default: {
  99       fprintf(stderr,"Wrong value for the -W flag. Exiting!\n");
 100       exit(127);
 101     }
 102     }
 103   }
 104   if ((out=check_option(in,n,'o','o')) != NULL) {
 105     stout=0;
 106     if (strlen(out) > 0)
 107       outfile=out;
 108   }
 109 }
 110
 111 void ordne(double *lyap,int *ord)
 112 {
 113   long i,j,maxi;
 114   double max;
 115
 116   for (i=0;i<dimemb;i++)
 117     ord[i]=i;
 118
 119   for (i=0;i<dimemb-1;i++)
 120     for (j=i+1;j<dimemb;j++)
 121       if (lyap[i] < lyap[j]) {
 122         max=lyap[i];
 123         lyap[i]=lyap[j];
 124         lyap[j]=max;
 125         maxi=ord[i];
 126         ord[i]=ord[j];
 127         ord[j]=maxi;
 128       }
 129 }
 130
 131 void make_pca(double *av)
 132 {
 133   unsigned int i,j,k,i1,i2,j1,j2,k1,k2;
 134   int *ord;
 135   double **mat,*matarray,*eig,*sp,hsp=0.0;
 136   FILE *fout=NULL;
 137
 138   check_alloc(ord=(int*)malloc(sizeof(int)*dimemb));
 139   check_alloc(eig=(double*)malloc(sizeof(double)*dimemb));
 140   check_alloc(matarray=(double*)malloc(sizeof(double)*dimemb*dimemb));
 141   check_alloc(mat=(double**)malloc(sizeof(double*)*dimemb));
 142   for (i=0;i<dimemb;i++)
 143     mat[i]=(double*)(matarray+i*dimemb);
 144
 145
 146   for (i=0;i<dimemb;i++) {
 147     i1=i/EMB;
 148     i2=(i%EMB)*DELAY;
 149     for (j=i;j<dimemb;j++) {
 150       j1=j/EMB;
 151       j2=(j%EMB)*DELAY;
 152       mat[i][j]=0.0;
 153       for (k=(EMB-1)*DELAY;k<LENGTH;k++)
 154         mat[i][j] += series[i1][k-i2]*series[j1][k-j2];
 155       mat[j][i]=(mat[i][j] /= (double)(LENGTH-(EMB-1)*DELAY));
 156     }
 157   }
 158
 159   eigen(mat,(unsigned long)dimemb,eig);
 160   ordne(eig,ord);
 161
 162   if (!stout) {
 163     fout=fopen(outfile,"w");
 164     if (verbosity&VER_INPUT)
 165       fprintf(stderr,"Opened %s for writing\n",outfile);
 166   }
 167   else {
 168     if (verbosity&VER_INPUT)
 169       fprintf(stderr,"Writing to stdout\n");
 170   }
 171
 172   for (i=0;i<dimemb;i++)
 173     if (write_values) {
 174       if (stout)
 175         fprintf(stdout,"%d %e\n",i,eig[i]);
 176       else
 177         fprintf(fout,"%d %e\n",i,eig[i]);
 178     }
 179     else {
 180       if (verbosity) {
 181         if (stout)
 182           fprintf(stdout,"#%d %e\n",i,eig[i]);
 183         else
 184           fprintf(fout,"#%d %e\n",i,eig[i]);
 185       }
 186     }
 187   if (write_vectors) {
 188     for (i=0;i<dimemb;i++) {
 189       for (j=0;j<dimemb;j++) {
 190         j1=ord[j];
 191         if (stout)
 192           fprintf(stdout,"%e ",mat[i][j1]);
 193         else
 194           fprintf(fout,"%e ",mat[i][j1]);
 195       }
 196       if (stout)
 197         fprintf(stdout,"\n");
 198       else
 199         fprintf(fout,"\n");
 200     }
 201   }
 202
 203   if (write_comp) {
 204     for (i=(EMB-1)*DELAY;i<LENGTH;i++) {
 205       for (j=0;j<LDIM;j++) {
 206         j1=ord[j];
 207         hsp=0.0;
 208         for (k=0;k<dimemb;k++) {
 209           k1=k/EMB;
 210           k2=(k%EMB)*DELAY;
 211           hsp += mat[k][j1]*(series[k1][i-k2]+av[k1]);
 212         }
 213         if (stout)
 214           fprintf(stdout,"%e ",hsp);
 215         else
 216           fprintf(fout,"%e ",hsp);
 217       }
 218       if (stout)
 219         fprintf(stdout,"\n");
 220       else
 221         fprintf(fout,"\n");
 222     }
 223   }
 224
 225   if (write_proj) {
 226     check_alloc(sp=(double*)malloc(sizeof(double)*LDIM));
 227     for (i=0;i<(EMB-1)*DELAY;i++) {
 228       for (j=0;j<DIM;j++)
 229         if (stout)
 230           fprintf(stdout,"%e ",series[j][i]+av[j]);
 231         else
 232           fprintf(fout,"%e ",series[j][i]+av[j]);
 233       if (stout)
 234         fprintf(stdout,"\n");
 235       else
 236         fprintf(fout,"\n");
 237     }
 238     for (i=(EMB-1)*DELAY;i<LENGTH;i++) {
 239       for (j=0;j<LDIM;j++) {
 240         j1=ord[j];
 241         sp[j]=0.0;
 242         for (k=0;k<dimemb;k++) {
 243           k1=k/EMB;
 244           k2=(k%EMB)*DELAY;
 245           sp[j] += mat[k][j1]*series[k1][i-k2];
 246         }
 247       }
 248       for (j=0;j<DIM;j++) {
 249         hsp=0.0;
 250         for (k=0;k<LDIM;k++) {
 251           k1=ord[k];
 252           hsp += mat[j*EMB][k1]*sp[k];
 253         }
 254         if (stout)
 255           fprintf(stdout,"%e ",hsp+av[j]);
 256         else
 257           fprintf(fout,"%e ",hsp+av[j]);
 258       }
 259       if (stout)
 260         fprintf(stdout,"\n");
 261       else
 262         fprintf(fout,"\n");
 263     }
 264     free(sp);
 265   }
 266
 267   if (!stout)
 268     fclose(fout);
 269 }
 270
 271 int main(int argc,char **argv)
 272 {
 273   char stdi=0;
 274   unsigned int i,j;
 275   double rms,*av;
 276
 277   if (scan_help(argc,argv))
 278     show_options(argv[0]);
 279
 280   scan_options(argc,argv);
 281
 282 #ifndef OMIT_WHAT_I_DO
 283   if (verbosity&VER_INPUT)
 284     what_i_do(argv[0],WID_STR);
 285 #endif
 286
 287   infile=search_datafile(argc,argv,NULL,verbosity);
 288   if (infile == NULL)
 289     stdi=1;
 290
 291   if (outfile == NULL) {
 292     if (!stdi) {
 293       check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
 294       strcpy(outfile,infile);
 295       strcat(outfile,".pca");
 296     }
 297     else {
 298       check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
 299       strcpy(outfile,"stdin.pca");
 300     }
 301   }
 302   if (!stout)
 303     test_outfile(outfile);
 304
 305   if (column == NULL)
 306     series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIM,"",dimset,
 307                                       verbosity);
 308   else
 309     series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIM,column,
 310                                       dimset,verbosity);
 311   dimemb=DIM*EMB;
 312   if (!projection_set)
 313     LDIM=dimemb;
 314   else {
 315     if (LDIM < 1) LDIM=1;
 316     if (LDIM > dimemb) LDIM=dimemb;
 317   }
 318
 319   check_alloc(av=(double*)malloc(sizeof(double)*DIM));
 320   for (j=0;j<DIM;j++) {
 321     av[j]=rms=0.0;
 322     variance(series[j],LENGTH,&av[j],&rms);
 323     for (i=0;i<LENGTH;i++)
 324       series[j][i] -= av[j];
 325   }
 326   make_pca(av);
 327
 328   return 0;
 329 }