X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=binaries%2Fsrc%2Fdisembl%2FTisean_3.0.1%2Fsource_c%2Fpca.c;fp=binaries%2Fsrc%2Fdisembl%2FTisean_3.0.1%2Fsource_c%2Fpca.c;h=cf3c42de3c0eeccff8f396497330b054d2e835f0;hb=a17c780665c109829426e062df4d75ff950725e0;hp=0000000000000000000000000000000000000000;hpb=f47da0247a9f9a8ac55571234064a0d3ded06b6c;p=jabaws.git

diff --git a/binaries/src/disembl/Tisean_3.0.1/source_c/pca.c b/binaries/src/disembl/Tisean_3.0.1/source_c/pca.c
new file mode 100644
index 0000000..cf3c42d
--- /dev/null
+++ b/binaries/src/disembl/Tisean_3.0.1/source_c/pca.c
@@ -0,0 +1,329 @@
+/*
+ *   This file is part of TISEAN
+ *
+ *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+ *
+ *   TISEAN is free software; you can redistribute it and/or modify
+ *   it under the terms of the GNU General Public License as published by
+ *   the Free Software Foundation; either version 2 of the License, or
+ *   (at your option) any later version.
+ *
+ *   TISEAN is distributed in the hope that it will be useful,
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *   GNU General Public License for more details.
+ *
+ *   You should have received a copy of the GNU General Public License
+ *   along with TISEAN; if not, write to the Free Software
+ *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+/*Author: Rainer Hegger Last modified: Jul 26, 2004 */
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include "routines/tsa.h"
+
+#define WID_STR "Performs a global PCA"
+
+unsigned long LENGTH=ULONG_MAX,exclude=0;
+unsigned int DIM=2,EMB=1,dimemb,LDIM=2,DELAY=1;
+unsigned int verbosity=0xff;
+char *outfile=NULL,stout=1,dim_set=0;
+unsigned int what_to_write=0,write_values=1,write_vectors=0;
+unsigned int write_comp=0,write_proj=0;
+unsigned int projection_set=0;
+char *infile=NULL,dimset=0,*column=NULL;
+double **series;
+
+void show_options(char *progname)
+{
+  what_i_do(progname,WID_STR);
+  fprintf(stderr," Usage: %s [options]\n",progname);
+  fprintf(stderr," Options:\n");
+  fprintf(stderr,"Everything not being a valid option will be interpreted"
+          " as a possible"
+          " datafile.\nIf no datafile is given stdin is read. Just - also"
+          " means stdin\n");
+  fprintf(stderr,"\t-l # of data to use [Default: whole file]\n");
+  fprintf(stderr,"\t-x # of lines to be ignore [Default: 0]\n");
+  fprintf(stderr,"\t-c columns to read [Default: 2]\n");
+  fprintf(stderr,"\t-m columns,embedding dim. to use [Default: 2,1]\n");
+  fprintf(stderr,"\t-d delay to use [Default: 1]\n");
+  fprintf(stderr,"\t-q projection dimension [Default: no projection]\n");
+  fprintf(stderr,"\t-W # what to write: [Default: 0]\n"
+	  "\t\t0 write eigenvalues only\n"
+	  "\t\t1 write eigenvectors\n"
+	  "\t\t2 write (projected) pca components\n"
+	  "\t\t3 write projected data\n");
+  fprintf(stderr,"\t-o output file name \n\t\t[Default: stdout; -o without "
+	  "value means 'datafile'.pca]\n");
+  fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
+          "0='only panic messages'\n\t\t"
+          "1='+ input/output messages'\n");
+  fprintf(stderr,"\t-h show these options\n");
+  exit(0);
+}
+
+void scan_options(int n,char **in)
+{
+  char *out;
+  
+  if ((out=check_option(in,n,'l','u')) != NULL)
+    sscanf(out,"%lu",&LENGTH);
+  if ((out=check_option(in,n,'x','u')) != NULL)
+    sscanf(out,"%lu",&exclude);
+  if ((out=check_option(in,n,'c','s')) != NULL)
+    column=out;
+  if ((out=check_option(in,n,'m','2')) != NULL) {
+    sscanf(out,"%u,%u",&DIM,&EMB);
+    dimset=1;
+  }
+  if ((out=check_option(in,n,'d','u')) != NULL)
+    sscanf(out,"%u",&DELAY);
+  if ((out=check_option(in,n,'q','u')) != NULL) {
+    sscanf(out,"%u",&LDIM);
+    projection_set=1;
+  }
+  if ((out=check_option(in,n,'V','u')) != NULL)
+    sscanf(out,"%u",&verbosity);
+  if ((out=check_option(in,n,'W','u')) != NULL) {
+    sscanf(out,"%u",&what_to_write);
+    switch(what_to_write) {
+    case 0: write_values=1;break;
+    case 1: write_values=0;write_vectors=1;break;
+    case 2: write_values=0;write_comp=1;break;
+    case 3: write_values=0;write_proj=1;break;
+    default: {
+      fprintf(stderr,"Wrong value for the -W flag. Exiting!\n");
+      exit(127);
+    }
+    }
+  }
+  if ((out=check_option(in,n,'o','o')) != NULL) {
+    stout=0;
+    if (strlen(out) > 0)
+      outfile=out;
+  }
+}
+
+void ordne(double *lyap,int *ord)
+{
+  long i,j,maxi;
+  double max;
+  
+  for (i=0;i<dimemb;i++)
+    ord[i]=i;
+
+  for (i=0;i<dimemb-1;i++)
+    for (j=i+1;j<dimemb;j++)
+      if (lyap[i] < lyap[j]) {
+	max=lyap[i];
+	lyap[i]=lyap[j];
+	lyap[j]=max;
+	maxi=ord[i];
+	ord[i]=ord[j];
+	ord[j]=maxi;
+      }
+}
+
+void make_pca(double *av)
+{
+  unsigned int i,j,k,i1,i2,j1,j2,k1,k2;
+  int *ord;
+  double **mat,*matarray,*eig,*sp,hsp=0.0;
+  FILE *fout=NULL;
+
+  check_alloc(ord=(int*)malloc(sizeof(int)*dimemb));
+  check_alloc(eig=(double*)malloc(sizeof(double)*dimemb));
+  check_alloc(matarray=(double*)malloc(sizeof(double)*dimemb*dimemb));
+  check_alloc(mat=(double**)malloc(sizeof(double*)*dimemb));
+  for (i=0;i<dimemb;i++)
+    mat[i]=(double*)(matarray+i*dimemb);
+
+  
+  for (i=0;i<dimemb;i++) {
+    i1=i/EMB;
+    i2=(i%EMB)*DELAY;
+    for (j=i;j<dimemb;j++) {
+      j1=j/EMB;
+      j2=(j%EMB)*DELAY;
+      mat[i][j]=0.0;
+      for (k=(EMB-1)*DELAY;k<LENGTH;k++)
+	mat[i][j] += series[i1][k-i2]*series[j1][k-j2];
+      mat[j][i]=(mat[i][j] /= (double)(LENGTH-(EMB-1)*DELAY));
+    }
+  }
+
+  eigen(mat,(unsigned long)dimemb,eig);
+  ordne(eig,ord);
+  
+  if (!stout) {
+    fout=fopen(outfile,"w");
+    if (verbosity&VER_INPUT)
+      fprintf(stderr,"Opened %s for writing\n",outfile);
+  }
+  else {
+    if (verbosity&VER_INPUT)
+      fprintf(stderr,"Writing to stdout\n");
+  }
+
+  for (i=0;i<dimemb;i++)
+    if (write_values) {
+      if (stout)
+	fprintf(stdout,"%d %e\n",i,eig[i]);
+      else
+	fprintf(fout,"%d %e\n",i,eig[i]);
+    }
+    else {
+      if (verbosity) {
+	if (stout)
+	  fprintf(stdout,"#%d %e\n",i,eig[i]);
+	else
+	  fprintf(fout,"#%d %e\n",i,eig[i]);
+      }
+    }
+  if (write_vectors) {
+    for (i=0;i<dimemb;i++) {
+      for (j=0;j<dimemb;j++) {
+	j1=ord[j];
+	if (stout)
+	  fprintf(stdout,"%e ",mat[i][j1]);
+	else
+	  fprintf(fout,"%e ",mat[i][j1]);
+      }
+      if (stout)
+	fprintf(stdout,"\n");
+      else
+	fprintf(fout,"\n");
+    }
+  }
+
+  if (write_comp) {
+    for (i=(EMB-1)*DELAY;i<LENGTH;i++) {
+      for (j=0;j<LDIM;j++) {
+	j1=ord[j];
+	hsp=0.0;
+	for (k=0;k<dimemb;k++) {
+	  k1=k/EMB;
+	  k2=(k%EMB)*DELAY;
+	  hsp += mat[k][j1]*(series[k1][i-k2]+av[k1]);
+	}
+	if (stout)
+	  fprintf(stdout,"%e ",hsp);
+	else
+	  fprintf(fout,"%e ",hsp);
+      }
+      if (stout)
+	fprintf(stdout,"\n");
+      else
+	fprintf(fout,"\n");
+    }
+  }
+
+  if (write_proj) {
+    check_alloc(sp=(double*)malloc(sizeof(double)*LDIM));
+    for (i=0;i<(EMB-1)*DELAY;i++) {
+      for (j=0;j<DIM;j++)
+	if (stout)
+	  fprintf(stdout,"%e ",series[j][i]+av[j]);
+	else
+	  fprintf(fout,"%e ",series[j][i]+av[j]);
+      if (stout)
+	fprintf(stdout,"\n");
+      else
+	fprintf(fout,"\n");
+    }
+    for (i=(EMB-1)*DELAY;i<LENGTH;i++) {
+      for (j=0;j<LDIM;j++) {
+	j1=ord[j];
+	sp[j]=0.0;
+	for (k=0;k<dimemb;k++) {
+	  k1=k/EMB;
+	  k2=(k%EMB)*DELAY;
+	  sp[j] += mat[k][j1]*series[k1][i-k2];
+	}
+      }
+      for (j=0;j<DIM;j++) {
+	hsp=0.0;
+	for (k=0;k<LDIM;k++) {
+	  k1=ord[k];
+	  hsp += mat[j*EMB][k1]*sp[k];
+	}
+	if (stout)
+	  fprintf(stdout,"%e ",hsp+av[j]);
+	else
+	  fprintf(fout,"%e ",hsp+av[j]);
+      }
+      if (stout)
+	fprintf(stdout,"\n");
+      else
+	fprintf(fout,"\n");
+    }
+    free(sp);
+  }
+
+  if (!stout)
+    fclose(fout);
+}
+
+int main(int argc,char **argv)
+{
+  char stdi=0;
+  unsigned int i,j;
+  double rms,*av;
+
+  if (scan_help(argc,argv))
+    show_options(argv[0]);
+
+  scan_options(argc,argv);
+
+#ifndef OMIT_WHAT_I_DO
+  if (verbosity&VER_INPUT)
+    what_i_do(argv[0],WID_STR);
+#endif
+
+  infile=search_datafile(argc,argv,NULL,verbosity);
+  if (infile == NULL)
+    stdi=1;
+
+  if (outfile == NULL) {
+    if (!stdi) {
+      check_alloc(outfile=(char*)calloc(strlen(infile)+5,(size_t)1));
+      strcpy(outfile,infile);
+      strcat(outfile,".pca");
+    }
+    else {
+      check_alloc(outfile=(char*)calloc((size_t)10,(size_t)1));
+      strcpy(outfile,"stdin.pca");
+    }
+  }
+  if (!stout)
+    test_outfile(outfile);
+
+  if (column == NULL)
+    series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIM,"",dimset,
+                                      verbosity);
+  else
+    series=(double**)get_multi_series(infile,&LENGTH,exclude,&DIM,column,
+                                      dimset,verbosity);
+  dimemb=DIM*EMB;
+  if (!projection_set)
+    LDIM=dimemb;
+  else {
+    if (LDIM < 1) LDIM=1;
+    if (LDIM > dimemb) LDIM=dimemb;
+  }
+
+  check_alloc(av=(double*)malloc(sizeof(double)*DIM));
+  for (j=0;j<DIM;j++) {
+    av[j]=rms=0.0;
+    variance(series[j],LENGTH,&av[j],&rms);
+    for (i=0;i<LENGTH;i++)
+      series[j][i] -= av[j];
+  }
+  make_pca(av);
+
+  return 0;
+}