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 part of the TISEAN randomize package for constraint surrogates
24 c binned autocorrelation function of unevenly sampled data
25 c author T. Schreiber (1999)
27 c-------------------------------------------------------------------
28 c get cost function specific options
30 subroutine opts_cost(ncol)
31 parameter(nhist=100000,nx=100000)
32 dimension hnorm(nhist), h0(nhist), h(nhist)
34 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
39 nbin=min(int(totbin*bininv)+1,nhist)
43 c-------------------------------------------------------------------
44 c print version information on cost function
46 subroutine what_cost()
47 call ptext("Cost function: binned autocorrelation function")
50 c-------------------------------------------------------------------
51 c print cost function specific usage message
53 subroutine usage_cost()
54 call ptext("Cost function options: -d# -D# [-W#]")
55 call popt("d","time span of one bin")
56 call popt("D","total time spanned")
58 . "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 (0)")
61 c-------------------------------------------------------------------
62 c initialise all that is needed for cost function
64 subroutine cost_init()
65 parameter(nhist=100000,nx=100000)
66 dimension hnorm(nhist), h0(nhist), h(nhist), x(nx,2)
67 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
68 common nmax,cost,temp,cmin,rate,x
74 il=int((x(n2,2)-x(n1,2))*bininv)+1
75 if(il.gt.nbin) goto 20
76 30 hnorm(il)=hnorm(il)+1.
79 40 if(hnorm(i).gt.0.) hnorm(i)=1./hnorm(i)
80 call sauto(nbin,bininv,h0)
83 c-------------------------------------------------------------------
84 c initial transformation on time series and its inverse
85 c here: sort by increasing sample times, no inversion necessary
86 c also normalise to unit variance, zero mean
88 subroutine cost_transform(nmax,mcmax,nxdum,x)
89 parameter(nhist=100000,nx=100000)
90 dimension hnorm(nhist), h0(nhist), h(nhist)
91 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
92 dimension x(nxdum,2), lx(nx)
94 call indexx(nmax,x(1,2),lx)
95 call index2sort(nmax,x(1,2),lx)
96 call index2sort(nmax,x,lx)
97 call normal1(nmax,x,sc,sd)
100 subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
101 dimension x(nxdum,2), y(nxdum,2)
102 parameter(nhist=100000,nx=100000)
103 dimension hnorm(nhist), h0(nhist), h(nhist)
104 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
108 10 y(n,1)=x(n,1)*sd+sc
111 c-------------------------------------------------------------------
112 c compute full cost function from scratch
114 subroutine cost_full(iv)
115 parameter(nhist=100000,nx=100000)
116 dimension hnorm(nhist), h0(nhist), h(nhist)
117 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
120 call sauto(nbin,bininv,h)
122 if(iv.ne.0) call dump()
125 c-------------------------------------------------------------------
126 c compute changed cost function on exchange of n1 and n2
128 subroutine cost_update(nn1,nn2,cmax,iaccept,iv)
129 parameter(nhist=100000,nx=100000)
130 dimension hnorm(nhist), h0(nhist), h(nhist)
131 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
132 dimension hcop(nhist), x(nx,2)
133 common nmax,cost,temp,cmin,rate,x
142 il=int((x(n1,2)-x(nn,2))*bininv)+1
143 if(il.gt.nbin) goto 1
144 20 hcop(il)=hcop(il)-x(n1,1)*x(nn,1)
147 il=int((x(nn,2)-x(n1,2))*bininv)+1
148 if(il.gt.nbin) goto 2
149 30 if(nn.ne.n2) hcop(il)=hcop(il)-x(nn,1)*x(n1,1)
152 il=int((x(n2,2)-x(nn,2))*bininv)+1
153 if(il.gt.nbin) goto 3
154 40 hcop(il)=hcop(il)-x(n2,1)*x(nn,1)
157 il=int((x(nn,2)-x(n2,2))*bininv)+1
158 if(il.gt.nbin) goto 4
159 50 hcop(il)=hcop(il)-x(nn,1)*x(n2,1)
162 il=int((x(n1,2)-x(nn,2))*bininv)+1
163 if(il.gt.nbin) goto 5
164 60 hcop(il)=hcop(il)+x(n1,1)*x(nn,1)
167 il=int((x(nn,2)-x(n1,2))*bininv)+1
168 if(il.gt.nbin) goto 6
169 70 if(nn.ne.n2) hcop(il)=hcop(il)+x(nn,1)*x(n1,1)
172 il=int((x(n2,2)-x(nn,2))*bininv)+1
173 if(il.gt.nbin) goto 7
174 80 hcop(il)=hcop(il)+x(n2,1)*x(nn,1)
177 il=int((x(nn,2)-x(n2,2))*bininv)+1
178 if(il.gt.nbin) goto 8
179 90 hcop(il)=hcop(il)+x(nn,1)*x(n2,1)
181 if(comp.ge.cmax) then
185 cost=comp ! if got here: accept
187 if(iv.ne.0) call panic(hcop)
192 c-------------------------------------------------------------------
193 c compute autocorrealtion from scratch
195 subroutine sauto(nbin,bininv,h)
198 common nmax,cost,temp,cmin,rate,x
205 il=int((x(n2,2)-x(n1,2))*bininv)+1
206 if(il.gt.nbin) goto 20
207 30 h(il)=h(il)+x(n2,1)*x(n1,1)
211 c-------------------------------------------------------------------
212 c weighted average of autocorrelation
215 parameter(nhist=100000,nx=100000)
216 dimension hnorm(nhist), h0(nhist), h(nhist),
218 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
221 if(iweight.eq.0) then
223 10 aver=max(aver,abs((h1(i)-h2(i))*hnorm(i)))
224 else if(iweight.eq.1) then
226 20 aver=aver+abs((h1(i)-h2(i))*hnorm(i))/real(i)
227 else if(iweight.eq.2) then
229 30 aver=aver+((h1(i)-h2(i))*hnorm(i))**2/real(i)
233 c-------------------------------------------------------------------
237 parameter(nhist=100000,nx=100000)
238 dimension hnorm(nhist), h0(nhist), h(nhist)
239 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
241 write(istderr(),'(5hgoal ,4g15.5)') (h0(n),n=1,min(4,nbin))
242 write(istderr(),'(5his ,4g15.5)') (h(n),n=1,min(4,nbin))
243 write(istderr(),'(5hmiss ,4g15.5)')
244 . (abs(h0(n)-h(n)),n=1,min(4,nbin))
245 write(istderr(),'()')
248 subroutine panic(hcop)
249 parameter(nhist=100000,nx=100000)
250 dimension hnorm(nhist), h0(nhist), h(nhist)
251 common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
255 write(istderr(),'(7hupdate ,4g15.5)') (hcop(n),n=1,min(4,nbin))
256 write(istderr(),'(7hfresh ,4g15.5)') (h(n),n=1,min(4,nbin))
257 write(istderr(),'(7hdiscr ,4g15.5)')
258 . (abs(hcop(n)-h(n)),n=1,min(4,nbin))
259 write(istderr(),'()')