website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/pc.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   embed using principal components
  23 C   author Thomas Schreiber (1998)
  24 c===========================================================================
  25       parameter(nx=1000000, me=500)
  26       dimension x(nx), c(me,me), d(me), xc(me), z(me,me)
  27       character*72 file, fout
  28       data id/1/, isvd/2/
  29       data iverb/1/
  30
  31       call whatido("embed using principal components",iverb)
  32       m=imust("m")
  33       if(m.gt.me) stop "svd: make me larger."
  34       id=ican("d",id)
  35       isvd=min(ican("q",isvd),m)
  36       nmax=ican("l",nx)
  37       nexcl=ican("x",0)
  38       jcol=ican("c",0)
  39       isout=igetout(fout,iverb)
  40
  41       call nthstring(1,file)
  42       call readfile(nmax,x,nexcl,jcol,file,iverb)
  43       call normal(nmax,x,sc,sd)
  44       call svd_vectors(nmax,m,id,x,c,z,d)
  45       if(iv_io(iverb).eq.1) write(istderr(),*)
  46      .   "#, fraction of variance, accumulative fraction"
  47       ctot=0.
  48       do 10 i=1,m
  49  10      ctot=ctot+d(m+1-i)
  50       cacc=0.
  51       do 20 i=1,m
  52          cacc=cacc+d(m+1-i)
  53  20      if(iv_io(iverb).eq.1)
  54      .      write(istderr(),*) i, d(m+1-i)/ctot, cacc/ctot
  55       if(file.eq."-") file="stdin"
  56       if(isout.eq.1) call addsuff(fout,file,"_pc")
  57       call outfile(fout,iunit,iverb)
  58       do 30 n=(m-1)*id+1,nmax
  59          do 40 i=1,isvd
  60             s=0
  61             do 50 j=1,m
  62  50            s=s+z(j,m+1-i)*x(n-(j-1)*id)
  63  40         xc(i)=s
  64  30      write(iunit,*) (xc(i),i=1,isvd)
  65       end
  66
  67       subroutine usage()
  68 c usage message
  69
  70       call whatineed(
  71      .   "-m# [-d# -q# -o outfile -l# -x# -c# -V# -h] file")
  72       call popt("m","initial embedding dimension")
  73       call popt("d","delay for initial embedding (1)")
  74       call popt("q","number of principal components (2)")
  75       call popt("l","number of values to be read (all)")
  76       call popt("x","number of values to be skipped (0)")
  77       call popt("c","column to be read (1 or file,#)")
  78       call pout("file_pc")
  79       call pall()
  80       stop
  81       end
  82
  83       subroutine svd_vectors(nmax,m,id,x,c,z,d)
  84       parameter(me=500)
  85       dimension x(nmax), c(me,*), d(m), w1(me), w2(me), z(me,*)
  86
  87       if(m.gt.me) stop "svd_vectors: make me larger."
  88       do 10 i=1,m
  89          do 10 j=i,m
  90             s=0.
  91             do 20 n=(m-1)*id+1,nmax
  92  20            s=s+x(n-(i-1)*id)*x(n-(j-1)*id)
  93             c(i,j)=s/(nmax-(m-1)*id)
  94  10         c(j,i)=c(i,j)
  95       call rs(me,m,c,d,1,z,w1,w2,ierr)
  96       end