website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/project.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   nonlinear noise reduction
  23 c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
  24 c      University Press (1997,2004)
  25 c   authors T. Schreiber, H. Kantz, R. Hegger (1998) based on earlier versions
  26 c===========================================================================
  27       parameter(nx=100000)
  28       dimension x(nx), x0(nx), xc(nx)
  29       character*72 file, fout
  30       data imax/1/
  31       data iverb/1/
  32
  33       call whatido("nonlinear noise reduction (see also: noise)",iverb)
  34       m=imust("m")
  35       nq=m-imust("q")
  36       eps=fmust("r",eps)
  37       kmin=imust("k")
  38       imax=ican("i",imax)
  39       nmax=ican("l",nx)
  40       nexcl=ican("x",0)
  41       jcol=ican("c",0)
  42       isout=igetout(fout,iverb)
  43
  44       call nthstring(1,file)
  45       call readfile(nmax,x,nexcl,jcol,file,iverb)
  46       if(file.eq."-") file="stdin"
  47       if(isout.eq.1) call addsuff(fout,file,"_")
  48
  49       do 10 n=1,nmax
  50  10      x0(n)=x(n)
  51       do 20 it=1,imax
  52          call clean(nmax,x,xc,m,kmin,nq,eps,iverb)
  53          if(fout.ne." ".or.isout.eq.1.or.it.eq.imax) then
  54             if(isout.eq.1) call suffix(fout,"c")
  55             call outfile(fout,iunit,iverb)
  56             do 30 n=1,nmax
  57  30            write(iunit,*) xc(n), x0(n)-xc(n)
  58             if(iunit.ne.istdout()) close(iunit)
  59             if(iv_io(iverb).eq.1) call writereport(nmax,fout)
  60          endif
  61          eps=0
  62          do 40 n=1,nmax
  63             eps=eps+(xc(n)-x(n))**2
  64  40         x(n)=xc(n)
  65          eps=sqrt(eps/nmax)
  66  20      if(iv_io(iverb).eq.1)
  67      .      write(istderr(),*) 'New diameter of neighbourhoods is ', eps
  68       end
  69
  70       subroutine usage()
  71 c usage message
  72
  73       call whatineed(
  74      .   "-m# -q# -r# -k# [-i# -o outfile -l# -x# -c# -V# -h] file")
  75       call popt("m","embedding dimension")
  76       call popt("q","dimension of manifold")
  77       call popt("r","radius of neighbourhoods")
  78       call popt("k","minimal number of neighbours")
  79       call popt("i","number of iterations (1)")
  80       call popt("l","number of values to be read (all)")
  81       call popt("x","number of values to be skipped (0)")
  82       call popt("c","column to be read (1 or file,#)")
  83       call pout("file_c, file_cc (etc.)")
  84       call pall()
  85       call ptext("Verbosity levels (add what you want):")
  86       call ptext("          1 = input/output" )
  87       call ptext("          2 = state of neighbour search")
  88       write(istderr(),'()')
  89       stop
  90       end
  91
  92       subroutine clean(nmax,y,yc,m,kmin,nq,d,iverb)
  93       parameter(im=100,ii=100000000,nx=100000,mm=15,small=0.0001)
  94       dimension y(nmax),yc(nmax),r(mm),ju(nx),c(mm,mm),cm(mm),
  95      .  jh(0:im*im),jpntr(nx),nlist(nx), zcm(mm,nx)
  96
  97       if(nmax.gt.nx.or.m.gt.mm) stop "clean: make mm/nx larger."
  98       sr=2*small+m-2                                        ! ${\rm tr}(1/r)=1$
  99       do 10 i=1,m
 100          r(i)=sr
 101  10      if(i.eq.m.or.i.eq.1) r(i)=sr/small
 102       do 20 i=1,nmax
 103  20      yc(i)=y(i)
 104       do 30 istep=1,2
 105          eps=d
 106          iu=nmax-m+1
 107          do 40 i=1,iu
 108  40         ju(i)=i+m-1
 109  1       call base(nmax,y,1,m,jh,jpntr,eps)
 110          iunp=0
 111          do 50 nn=1,iu                                        ! find neighbours
 112             n=ju(nn)
 113             call neigh(nmax,y,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
 114             if(nfound.lt.kmin) then               ! not enough neighbours found
 115                iunp=iunp+1                                ! mark for next sweep
 116                ju(iunp)=n
 117             else                                      ! fine: enough neighbours
 118                do 90 i=1,m                              ! centre of mass vector
 119                   s=0
 120                   do 100 np=1,nfound
 121  100                 s=s+y(nlist(np)-m+i)
 122  90               cm(i)=s/nfound
 123                if(istep.eq.1) then                  ! just store centre of mass
 124                   do 110 i=1,m
 125  110                 zcm(i,n)=cm(i)
 126                else
 127                   do 120 i=1,m                ! corrected centre of mass vector
 128                      s=0
 129                      do 130 np=1,nfound
 130  130                    s=s+zcm(i,nlist(np))
 131  120                 cm(i)=2*cm(i)-s/nfound
 132                   do 140 i=1,m                      ! compute covariance matrix
 133                      do 140 j=i,m
 134                         s=0
 135                         do 150 np=1,nfound
 136                            jm=nlist(np)-m
 137  150                       s=s+(y(jm+i)-cm(i))*(y(jm+j)-cm(j))
 138                         c(i,j)=r(i)*r(j)*s/nfound
 139  140                    c(j,i)=c(i,j)
 140                  call eigen(c,m)               ! find eigenvectors (decreasing)
 141                  do 160 i=1,m
 142                     s=0
 143                     do 170 iq=m-nq+1,m
 144                        do 170 j=1,m
 145  170                      s=s+(y(n-m+j)-cm(j))*c(i,iq)*c(j,iq)*r(j)
 146  160                yc(n-m+i)=yc(n-m+i)-s/r(i)/r(i)
 147                endif
 148             endif
 149  50         continue
 150          iu=iunp
 151          if(iv_uncorr(iverb).eq.1)
 152      .      write(istderr(),*) "With ", eps, iunp, " uncorrected"
 153          eps=eps*sqrt(2.)
 154  30      if(iunp.ne.0) goto 1
 155       end
 156
 157 c driver for diagonalisation routines
 158 c see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
 159 c      University Press (1997)
 160 c Copyright (C) T. Schreiber (1997)
 161
 162       subroutine eigen(c,kk)
 163       parameter(md=15)
 164       dimension c(md,md),d(md),w1(md),w2(md),z(md,md)
 165       if(kk.gt.md) stop "eigen: make md larger."
 166
 167       call rs(md,kk,c,d,1,z,w1,w2,ierr)
 168       do 10 i=1,kk
 169          do 10 j=1,kk
 170  10         c(i,j)=z(i,kk+1-j)
 171       end
 172