IUPred wrapper, tester and new binary for X64 linux systems
[jabaws.git] / 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 }