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 autocorrelation function with periodic continuation
25 c author T. Schreiber (1999)
27 c-------------------------------------------------------------------
28 c get cost function specific options
30 subroutine opts_cost(ncol)
31 parameter(mlag=100000)
32 dimension c0(mlag), c(mlag)
33 common /costcom/ nlag, c0, c, sd, sc, iweight
40 c-------------------------------------------------------------------
41 c print version information on cost function
43 subroutine what_cost()
44 call ptext("Cost function: periodic autocorrelation")
47 c-------------------------------------------------------------------
48 c print cost function specific usage message
50 subroutine usage_cost()
51 call ptext("Cost function options: -D# [-W#]")
52 call popt("D","number of lags")
54 . "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 3=max(c)/lag (0)")
57 c-------------------------------------------------------------------
58 c initialise all that is needed for cost function
60 subroutine cost_init()
61 parameter(mlag=100000)
62 dimension c0(mlag), c(mlag)
63 common /costcom/ nlag, c0, c, sd, sc, iweight
65 if(nlag.gt.mlag) write(istderr(),'(a)')
66 . "truncated to ", mlag," lags"
71 c-------------------------------------------------------------------
72 c initial transformation on time series and its inverse
74 subroutine cost_transform(nmax,mcmax,nxdum,x)
76 parameter(mlag=100000)
77 dimension c0(mlag), c(mlag)
78 common /costcom/ nlag, c0, c, sd, sc, iweight
80 call normal1(nmax,x,sc,sd)
83 subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
84 dimension x(nmax), y(nmax)
85 parameter(mlag=100000)
86 dimension c0(mlag), c(mlag)
87 common /costcom/ nlag, c0, c, sd, sc, iweight
93 c-------------------------------------------------------------------
94 c compute full cost function from scratch
96 subroutine cost_full(iv)
97 parameter(mlag=100000)
98 dimension c0(mlag), c(mlag)
99 common /costcom/ nlag, c0, c, sd, sc, iweight
105 10 call aver(cc,c0(n)-c(n),n)
109 c-------------------------------------------------------------------
110 c compute changed cost function on exchange of n1 and n2
112 subroutine cost_update(n1,n2,cmax,iaccept,iv)
113 parameter(mlag=100000,nx=100000)
114 dimension c0(mlag), c(mlag), ccop(mlag), x(nx)
115 common /costcom/ nlag, c0, c, sd, sc, iweight
116 common nmax,cost,temp,cmin,rate,x
124 if(nd1.lt.1) nd1=nd1+nmax
125 if(nd1.ne.n2) cc=cc+dx*x(nd1)
127 if(nu1.gt.nmax) nu1=nu1-nmax
128 if(nu1.ne.n2) cc=cc+dx*x(nu1)
130 if(nd2.lt.1) nd2=nd2+nmax
131 if(nd2.ne.n1) cc=cc-dx*x(nd2)
133 if(nu2.gt.nmax) nu2=nu2-nmax
134 if(nu2.ne.n1) cc=cc-dx*x(nu2)
135 call aver(comp,c0(n)-cc,n)
136 if(comp.ge.cmax) return
138 cost=comp ! if got here: accept
145 c-------------------------------------------------------------------
146 c compute autocorrelation from scratch
148 subroutine auto(nlag,c)
150 dimension c(*), x(nx)
151 common nmax,cost,temp,cmin,rate,x
157 if(ii.lt.1) ii=ii+nmax
162 c-------------------------------------------------------------------
163 c weighted average of autocorrelation
165 subroutine aver(cav,dc,n)
166 parameter(mlag=100000)
167 dimension c0(mlag), c(mlag)
168 common /costcom/ nlag, c0, c, sd, sc, iweight
171 if(iweight.eq.0) then
172 cav=max(cav,abs(dc)/real(nmax))
173 else if(iweight.eq.1) then
174 cav=cav+abs(dc)/real(nmax*n)
175 else if(iweight.eq.2) then
176 cav=cav+(dc/real(nmax*n))**2
178 cav=max(cav,abs(dc)/real(nmax*n))