1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
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.
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.
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
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/
33 call whatido("Create Multivariate Surrogate data",iverb)
36 nsur=min(999,ican("n",nsur))
39 r=rand(sqrt(abs(fcan("I",0.0))))
41 call columns(mc,mx,icol)
42 if(mcmax.eq.0) mcmax=max(1,mc)
43 isout=igetout(fout,iverb)
45 call nthstring(1,file)
46 call xreadfile(nmax,mcmax,nx,xx,nexcl,icol,file,iverb)
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")
55 if(nsur.gt.1.and.isout.eq.1)
56 . write(fout(index(fout," ")-3:72),'(i3.3)') isur
63 call store_spec(nmaxp,xamp(1,m),0)
64 call sort(nmaxp,xsort(1,m),list)
67 call rank(nmaxp,rwork,list)
68 20 call index2sort(nmaxp,x(1,m),list)
76 dspec=toxspec(nmaxp,mcmax,nx,xamp,y)
77 if(imax.ge.0.and.it.ge.imax) goto 2
79 60 call todist(nmaxp,xsort(1,m),y(1,m),x(1,m))
80 if(dspec.lt.ds0) goto 1
83 call xwritefile(nmaxp,mcmax,nx,y,fout,iverb)
85 call xwritefile(nmaxp,mcmax,nx,x,fout,iverb)
87 10 if(iv_surr(iverb).eq.1) write(istderr(),*)
88 . fout(1:index(fout," ")), ' (', it,
89 . ' iterations, relative discrepancy ', dspec, ')'
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)")
107 call ptext("Verbosity levels (add what you want):")
108 call ptext(" 1 = input/output" )
109 call ptext(" 2 = iterations / discrepancy")
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)
118 if(nmax.gt.nx.or.mmax.gt.mx) stop "toxspec: make nx/mx larger."
119 call rffti1(nmax,w2,iw)
123 call rfftf1(nmax,x(1,m),w1,w2,iw)
125 30 x(n,m)=x(n,m)/real(nmax)
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))
131 10 if(mod(nmax,2).eq.0) x(nmax,m)=sqrt(a(nmax,m))
135 60 goal(m)=x(2*n-1,m)-a(2*n-1,m)
138 50 x(2*n-1,m)=alpha+a(2*n-1,m)
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))
146 70 call rfftb1(nmax,x(1,m),w1,w2,iw)
150 90 toxspec=toxspec+(x(n,m)-w(n,m))**2
151 toxspec=sqrt((toxspec/nmax)/mmax)
154 function alp(mmax,goal)
162 10 f2=f2+sin(goal(m))
166 20 scos=scos+cos(alp-goal(m))
167 if(scos.lt.0) alp=alp+pi
170 subroutine todist(nmax,dist,x,y)
172 dimension x(nmax), dist(nmax), y(nmax), list(nx)
174 if(nmax.gt.nx) stop "todist: make nx larger."
175 call rank(nmax,x,list)
177 10 y(n)=dist(list(n))