website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/randomize/cost/auto.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
  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: 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(nn1,nn2,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       n1=min(nn1,nn2)
 119       n2=max(nn1,nn2)
 120       comp=0
 121       iaccept=0
 122       do 10 n=1,nlag
 123          cc=c(n)
 124          dx=x(n2)-x(n1)
 125          if(n1-n.ge.1) cc=cc+dx*x(n1-n)
 126          if(n2+n.le.nmax) cc=cc-dx*x(n2+n)
 127          if(n2-n1.eq.n) goto 1
 128          if(n1+n.le.nmax) cc=cc+dx*x(n1+n)
 129          if(n2-n.ge.1) cc=cc-dx*x(n2-n)
 130  1       call aver(comp,c0(n)-cc,n)
 131          if(comp.ge.cmax) return
 132  10      ccop(n)=cc
 133       cost=comp  ! if got here: accept
 134       iaccept=1
 135       call exch(n1,n2)
 136       do 20 n=1,nlag
 137  20      c(n)=ccop(n)
 138       end
 139
 140 c-------------------------------------------------------------------
 141 c compute autocorrelation from scratch
 142 c
 143       subroutine auto(nlag,c)
 144       parameter(nx=100000)
 145       dimension c(*), x(nx)
 146       common nmax,cost,temp,cmin,rate,x
 147
 148       do 10 n=1,nlag
 149          cc=0
 150          do 20 i=n+1,nmax
 151  20         cc=cc+x(i-n)*x(i)
 152  10      c(n)=cc
 153       end
 154
 155 c-------------------------------------------------------------------
 156 c weighted average of autocorrelation
 157 c
 158       subroutine aver(cav,dc,n)
 159       parameter(mlag=100000)
 160       dimension c0(mlag), c(mlag)
 161       common /costcom/ nlag, c0, c, sd, sc, iweight
 162       common nmax
 163
 164       if(iweight.eq.0) then
 165          cav=max(cav,abs(dc)/real(nmax-n))
 166       else if(iweight.eq.1) then
 167          cav=cav+abs(dc)/real((nmax-n)*n)
 168       else if(iweight.eq.2) then
 169          cav=cav+(dc/real((nmax-n)*n))**2
 170       else
 171          cav=max(cav,abs(dc)/real((nmax-n)*n))
 172       endif
 173       end
 174