binaries/src/disembl/Tisean_3.0.1/source_f/randomize/cost/autop.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   part of the TISEAN randomize package for constraint surrogates
  23 c   cost function
  24 c   autocorrelation function with periodic continuation
  25 c   author T. Schreiber (1999)
  26 c
  27 c-------------------------------------------------------------------
  28 c get cost function specific options
  29 c
  30       subroutine opts_cost(ncol)
  31       parameter(mlag=100000)
  32       dimension c0(mlag), c(mlag)
  33       common /costcom/ nlag, c0, c, sd, sc, iweight
  34
  35       nlag=imust('D')
  36       iweight=ican('W',0)
  37       ncol=1
  38       end
  39
  40 c-------------------------------------------------------------------
  41 c print version information on cost function
  42 c
  43       subroutine what_cost()
  44       call ptext("Cost function: periodic autocorrelation")
  45       end
  46
  47 c-------------------------------------------------------------------
  48 c print cost function specific usage message
  49 c
  50       subroutine usage_cost()
  51       call ptext("Cost function options: -D# [-W#]")
  52       call popt("D","number of lags")
  53       call popt("W",
  54      .   "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 3=max(c)/lag (0)")
  55       end
  56
  57 c-------------------------------------------------------------------
  58 c initialise all that is needed for cost function
  59 c
  60       subroutine cost_init()
  61       parameter(mlag=100000)
  62       dimension c0(mlag), c(mlag)
  63       common /costcom/ nlag, c0, c, sd, sc, iweight
  64
  65       if(nlag.gt.mlag) write(istderr(),'(a)')
  66      .   "truncated to ", mlag," lags"
  67       nlag=min(mlag,nlag)
  68       call auto(nlag,c0)
  69       end
  70
  71 c-------------------------------------------------------------------
  72 c initial transformation on time series and its inverse
  73 c
  74       subroutine cost_transform(nmax,mcmax,nxdum,x)
  75       dimension x(nmax)
  76       parameter(mlag=100000)
  77       dimension c0(mlag), c(mlag)
  78       common /costcom/ nlag, c0, c, sd, sc, iweight
  79
  80       call normal1(nmax,x,sc,sd)
  81       end
  82
  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
  88
  89       do 10 n=1,nmax
  90  10      y(n)=x(n)*sd+sc
  91       end
  92
  93 c-------------------------------------------------------------------
  94 c compute full cost function from scratch
  95 c
  96       subroutine cost_full(iv)
  97       parameter(mlag=100000)
  98       dimension c0(mlag), c(mlag)
  99       common /costcom/ nlag, c0, c, sd, sc, iweight
 100       common nmax,cost
 101
 102       call auto(nlag,c)
 103       cc=0
 104       do 10 n=1,nlag
 105  10      call aver(cc,c0(n)-c(n),n)
 106       cost=cc
 107       end
 108
 109 c-------------------------------------------------------------------
 110 c compute changed cost function on exchange of n1 and n2
 111 c
 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
 117
 118       comp=0
 119       iaccept=0
 120       do 10 n=1,nlag
 121          cc=c(n)
 122          dx=x(n2)-x(n1)
 123          nd1=n1-n
 124          if(nd1.lt.1) nd1=nd1+nmax
 125          if(nd1.ne.n2) cc=cc+dx*x(nd1)
 126          nu1=n1+n
 127          if(nu1.gt.nmax) nu1=nu1-nmax
 128          if(nu1.ne.n2) cc=cc+dx*x(nu1)
 129          nd2=n2-n
 130          if(nd2.lt.1) nd2=nd2+nmax
 131          if(nd2.ne.n1) cc=cc-dx*x(nd2)
 132          nu2=n2+n
 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
 137  10      ccop(n)=cc
 138       cost=comp  ! if got here: accept
 139       iaccept=1
 140       call exch(n1,n2)
 141       do 20 n=1,nlag
 142  20      c(n)=ccop(n)
 143       end
 144
 145 c-------------------------------------------------------------------
 146 c compute autocorrelation from scratch
 147 c
 148       subroutine auto(nlag,c)
 149       parameter(nx=100000)
 150       dimension c(*), x(nx)
 151       common nmax,cost,temp,cmin,rate,x
 152
 153       do 10 n=1,nlag
 154          cc=0
 155          do 20 i=1,nmax
 156             ii=i-n
 157             if(ii.lt.1) ii=ii+nmax
 158  20         cc=cc+x(ii)*x(i)
 159  10      c(n)=cc
 160       end
 161
 162 c-------------------------------------------------------------------
 163 c weighted average of autocorrelation
 164 c
 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
 169       common nmax
 170
 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
 177       else
 178          cav=max(cav,abs(dc)/real(nmax*n))
 179       endif
 180       end
 181