binaries/src/disembl/Tisean_3.0.1/source_f/surrogates.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   Create multivariate surrogate data
  23 c   author T. Schreiber (1999)
  24 c===========================================================================
  25       parameter(nx=100000,mx=20)
  26       dimension xx(nx,mx), x(nx,mx), y(nx,mx), xamp(nx,mx),
  27      .   xsort(nx,mx), list(nx), icol(mx), rwork(nx)
  28       character*72 file, fout
  29       data nsur/1/, imax/-1/
  30       external rand
  31       data iverb/15/
  32
  33       call whatido("Create Multivariate Surrogate data",iverb)
  34       nmax=ican("l",nx)
  35       nexcl=ican("x",0)
  36       nsur=min(999,ican("n",nsur))
  37       imax=ican("i",imax)
  38       ispec=lopt("S",1)
  39       r=rand(sqrt(abs(fcan("I",0.0))))
  40       mcmax=ican("m",0)
  41       call columns(mc,mx,icol)
  42       if(mcmax.eq.0) mcmax=max(1,mc)
  43       isout=igetout(fout,iverb)
  44
  45       call nthstring(1,file)
  46       call xreadfile(nmax,mcmax,nx,xx,nexcl,icol,file,iverb)
  47       nmaxp=nless(nmax)
  48       if(nmaxp.ne.nmax.and.iv_io(iverb).eq.1)
  49      .   write(istderr(),*) "surrogates: using first ", nmaxp
  50       if(file.eq."-") file="stdin"
  51       if(isout.eq.1) call addsuff(fout,file,"_surr")
  52       if(nsur.gt.1.and.isout.eq.1) call suffix(fout,"_000")
  53
  54       do 10 isur=1,nsur
  55          if(nsur.gt.1.and.isout.eq.1)
  56      .      write(fout(index(fout," ")-3:72),'(i3.3)') isur
  57          do 20 m=1,mcmax
  58             do 30 n=1,nmaxp
  59                x(n,m)=xx(n,m)
  60                y(n,m)=x(n,m)
  61                xamp(n,m)=x(n,m)
  62  30            xsort(n,m)=x(n,m)
  63             call store_spec(nmaxp,xamp(1,m),0)
  64             call sort(nmaxp,xsort(1,m),list)
  65             do 40 n=1,nmaxp
  66  40            rwork(n)=rand(0.0)
  67             call rank(nmaxp,rwork,list)
  68  20         call index2sort(nmaxp,x(1,m),list)
  69          it=-1
  70          dspec=r1mach(2)
  71  1       it=it+1
  72          do 50 m=1,mcmax
  73             do 50 n=1,nmaxp
  74  50            y(n,m)=x(n,m)
  75          ds0=dspec
  76          dspec=toxspec(nmaxp,mcmax,nx,xamp,y)
  77          if(imax.ge.0.and.it.ge.imax) goto 2
  78          do 60 m=1,mcmax
  79  60         call todist(nmaxp,xsort(1,m),y(1,m),x(1,m))
  80          if(dspec.lt.ds0) goto 1
  81  2       continue
  82          if(ispec.gt.0) then
  83             call xwritefile(nmaxp,mcmax,nx,y,fout,iverb)
  84          else
  85             call xwritefile(nmaxp,mcmax,nx,x,fout,iverb)
  86          endif
  87  10      if(iv_surr(iverb).eq.1) write(istderr(),*)
  88      .      fout(1:index(fout," ")), ' (', it,
  89      .      ' iterations, relative discrepancy ', dspec,   ')'
  90       end
  91
  92       subroutine usage()
  93 c usage message
  94
  95       call whatineed(
  96      .   "[-n# -i# -S -I# -o outfile -l# -x# -m# -c#[,#] -V# -h] file")
  97       call popt("n","number of surrogates (1)")
  98       call popt("i","number of iterations (until no change)")
  99       call popt("S","make spectrum exact rather than distribution")
 100       call popt("I","seed for random numbers")
 101       call popt("l","number of values to be read (all)")
 102       call popt("x","number of values to be skipped (0)")
 103       call popt("m","number of columns to be read (1)")
 104       call popt("c","columns to be read (1)")
 105       call pout("file_surr(_nnn)")
 106       call pall()
 107       call ptext("Verbosity levels (add what you want):")
 108       call ptext("          1 = input/output" )
 109       call ptext("          2 = iterations / discrepancy")
 110       stop
 111       end
 112
 113       function toxspec(nmax,mmax,nxx,a,x)
 114       parameter(nx=100000,mx=20,tol=1e-5)
 115       dimension x(nxx,mmax), a(nxx,mmax), w(nx,mx), w1(nx),
 116      .   w2(nx), iw(15), goal(mx)
 117
 118       if(nmax.gt.nx.or.mmax.gt.mx) stop "toxspec: make nx/mx larger."
 119       call rffti1(nmax,w2,iw)
 120       do 10 m=1,mmax
 121          do 20 n=1,nmax
 122  20         w(n,m)=x(n,m)
 123          call rfftf1(nmax,x(1,m),w1,w2,iw)
 124          do 30 n=1,nmax
 125  30         x(n,m)=x(n,m)/real(nmax)
 126          x(1,m)=sqrt(a(1,m))
 127          do 40 n=2,(nmax+1)/2
 128             pha=atan2(x(2*n-1,m),x(2*n-2,m))
 129             x(2*n-2,m)=sqrt(a(2*n-2,m))
 130  40         x(2*n-1,m)=pha
 131  10      if(mod(nmax,2).eq.0) x(nmax,m)=sqrt(a(nmax,m))
 132       if(mmax.gt.1) then
 133          do 50 n=2,(nmax+1)/2
 134             do 60 m=1,mmax
 135  60            goal(m)=x(2*n-1,m)-a(2*n-1,m)
 136             alpha=alp(mmax,goal)
 137             do 50 m=1,mmax
 138  50            x(2*n-1,m)=alpha+a(2*n-1,m)
 139       endif
 140       do 70 m=1,mmax
 141          do 80 n=2,(nmax+1)/2
 142             c=x(2*n-2,m)*cos(x(2*n-1,m))
 143             s=x(2*n-2,m)*sin(x(2*n-1,m))
 144             x(2*n-1,m)=s
 145  80         x(2*n-2,m)=c
 146  70      call rfftb1(nmax,x(1,m),w1,w2,iw)
 147       toxspec=0
 148       do 90 m=1,mmax
 149          do 90 n=1,nmax
 150  90         toxspec=toxspec+(x(n,m)-w(n,m))**2
 151       toxspec=sqrt((toxspec/nmax)/mmax)
 152       end
 153
 154       function alp(mmax,goal)
 155       dimension goal(mmax)
 156       data pi/3.1415926/
 157
 158       f1=0
 159       f2=0
 160       do 10 m=1,mmax
 161          f1=f1+cos(goal(m))
 162  10      f2=f2+sin(goal(m))
 163       alp=atan2(f2,f1)
 164       scos=0
 165       do 20 m=1,mmax
 166  20      scos=scos+cos(alp-goal(m))
 167       if(scos.lt.0) alp=alp+pi
 168       end
 169
 170       subroutine todist(nmax,dist,x,y)
 171       parameter(nx=100000)
 172       dimension x(nmax), dist(nmax), y(nmax), list(nx)
 173
 174       if(nmax.gt.nx) stop "todist: make nx larger."
 175       call rank(nmax,x,list)
 176       do 10 n=1,nmax
 177  10      y(n)=dist(list(n))
 178       end
 179