website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/lazy.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   lazy.f
  23 c   simple nonlinear noise reduction
  24 c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
  25 c      University Press (1997,2004)
  26 c   author T. Schreiber (1998)
  27 c===========================================================================
  28       parameter(nx=1000000)
  29       dimension x(nx), x0(nx), xc(nx)
  30       character*72 file, fout
  31       data eps/0./, frac/0./, imax/1/
  32       data iverb/1/
  33
  34       call whatido("simple nonlinear noise reduction",iverb)
  35       m=imust("m")
  36       eps=fcan("r",eps)
  37       frac=fcan("v",frac)
  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       if(eps.eq.0.and.frac.eq.0.) call usage()
  44
  45       call nthstring(1,file)
  46       call readfile(nmax,x,nexcl,jcol,file,iverb)
  47       if(file.eq."-") file="stdin"
  48       if(isout.eq.1) call addsuff(fout,file,"_l")
  49       call rms(nmax,x,sc,sd)
  50       if(frac.gt.0) eps=sd*frac
  51       do 10 n=1,nmax
  52  10      x0(n)=x(n)
  53       do 20 it=1,imax
  54          call nrlazy(nmax,x,xc,m,eps)
  55          if(fout.ne." ".or.isout.eq.1.or.it.eq.imax) then
  56             if(isout.eq.1) call suffix(fout,"c")
  57             call outfile(fout,iunit,iverb)
  58             do 30 n=1,nmax
  59  30            write(iunit,*) xc(n), x0(n)-xc(n)
  60             if(iunit.ne.istdout()) close(iunit)
  61             if(iv_io(iverb).eq.1) call writereport(nmax,fout)
  62          endif
  63          eps=0
  64          do 40 n=1,nmax
  65             eps=eps+(xc(n)-x(n))**2
  66  40         x(n)=xc(n)
  67          eps=sqrt(eps/nmax)
  68          if(eps.eq.0.) then
  69             if(iv_io(iverb).eq.1) write(istderr(),*)
  70      .      'Zero correction, finished'
  71             stop
  72          endif
  73  20      if(iv_io(iverb).eq.1) write(istderr(),*)
  74      .      'New diameter of neighbourhoods is ', eps
  75       end
  76
  77       subroutine usage()
  78 c usage message
  79
  80       call whatineed(
  81      .   "-m# [-r# | -v#] [-i# -o outfile -l# -x# -c# -V# -h] file")
  82       call ptext("either -r or -v must be present")
  83       call popt("m","embedding dimension")
  84       call popt("r","absolut radius of neighbourhoods")
  85       call popt("v","same as fraction of standard deviation")
  86       call popt("i","number of iterations (1)")
  87       call popt("l","number of values to be read (all)")
  88       call popt("x","number of values to be skipped (0)")
  89       call popt("c","column to be read (1 or file,#)")
  90       call pout("file_lc,  file_lcc (etc.)")
  91       call pall()
  92       stop
  93       end
  94
  95       subroutine nrlazy(nmax,y,yc,m,eps)
  96       parameter(im=100,ii=100000000,nx=1000000)
  97       dimension y(nmax),yc(nmax),jh(0:im*im),jpntr(nx),nlist(nx)
  98
  99       if(nmax.gt.nx) stop "nrlazy: make nx larger."
 100       call base(nmax,y,1,m,jh,jpntr,eps)
 101       do 10 n=1,nmax
 102  10      yc(n)=y(n)
 103       do 20 n=m,nmax
 104          call neigh(nmax,y,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
 105          av=0
 106          do 30 nn=1,nfound
 107  30         av=av+y(nlist(nn)-(m-1)/2)              ! average middle coordinate
 108  20      yc(n-(m-1)/2)=av/nfound
 109       end