website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/randomize/cost/uneven.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   binned autocorrelation function of unevenly sampled data
  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(nhist=100000,nx=100000)
  32       dimension hnorm(nhist), h0(nhist), h(nhist)
  33       character*80 filet
  34       common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
  35
  36       iweight=ican('W',0)
  37       bininv=1./fmust("d")
  38       totbin=fmust("D")
  39       nbin=min(int(totbin*bininv)+1,nhist)
  40       ncol=2
  41       end
  42
  43 c-------------------------------------------------------------------
  44 c print version information on cost function
  45 c
  46       subroutine what_cost()
  47       call ptext("Cost function: binned autocorrelation function")
  48       end
  49
  50 c-------------------------------------------------------------------
  51 c print cost function specific usage message
  52 c
  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")
  57       call popt("W",
  58      .   "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 (0)")
  59       end
  60
  61 c-------------------------------------------------------------------
  62 c initialise all that is needed for cost function
  63 c
  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
  69
  70       do 10 i=1,nbin
  71  10      hnorm(i)=0.
  72       do 20 n1=1,nmax
  73          do 30 n2=n1,nmax
  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.
  77  20      continue
  78       do 40 i=1,nbin
  79  40      if(hnorm(i).gt.0.) hnorm(i)=1./hnorm(i)
  80       call sauto(nbin,bininv,h0)
  81       end
  82
  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
  87 c
  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)
  93
  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)
  98       end
  99
 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
 105
 106       do 10 n=1,nmax
 107          y(n,2)=x(n,2)
 108  10      y(n,1)=x(n,1)*sd+sc
 109       end
 110
 111 c-------------------------------------------------------------------
 112 c compute full cost function from scratch
 113 c
 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
 118       common nmax,cost
 119
 120       call sauto(nbin,bininv,h)
 121       cost=aver(h0,h)
 122       if(iv.ne.0) call dump()
 123       end
 124
 125 c-------------------------------------------------------------------
 126 c compute changed cost function on exchange of n1 and n2
 127 c
 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
 134
 135       n1=min(nn1,nn2)
 136       n2=max(nn1,nn2)
 137       comp=0
 138       iaccept=0
 139       do 10 i=1,nbin
 140  10      hcop(i)=h(i)
 141       do 20 nn=n1-1,1,-1
 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)
 145  1    continue
 146       do 30 nn=n1,nmax
 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)
 150  2    continue
 151       do 40 nn=n2-1,1,-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)
 155  3    continue
 156       do 50 nn=n2,nmax
 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)
 160  4    call exch(n1,n2)
 161       do 60 nn=n1-1,1,-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)
 165  5    continue
 166       do 70 nn=n1,nmax
 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)
 170  6    continue
 171       do 80 nn=n2-1,1,-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)
 175  7    continue
 176       do 90 nn=n2,nmax
 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)
 180  8    comp=aver(h0,hcop)
 181       if(comp.ge.cmax) then
 182          call exch(n1,n2)
 183          return
 184       endif
 185       cost=comp  ! if got here: accept
 186       iaccept=1
 187       if(iv.ne.0) call panic(hcop)
 188       do 100 i=1,nbin
 189  100     h(i)=hcop(i)
 190       end
 191
 192 c-------------------------------------------------------------------
 193 c compute autocorrealtion from scratch
 194 c
 195       subroutine sauto(nbin,bininv,h)
 196       parameter(nx=100000)
 197       dimension h(*)
 198       common nmax,cost,temp,cmin,rate,x
 199       dimension x(nx,2)
 200
 201       do 10 i=1,nbin
 202  10      h(i)=0
 203       do 20 n1=1,nmax
 204          do 30 n2=n1,nmax
 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)
 208  20      continue
 209       end
 210
 211 c-------------------------------------------------------------------
 212 c weighted average of autocorrelation
 213 c
 214       function aver(h1,h2)
 215       parameter(nhist=100000,nx=100000)
 216       dimension hnorm(nhist), h0(nhist), h(nhist),
 217      .   h1(*), h2(*)
 218       common /costcom/ bininv, nbin, hnorm, h0, h, sd, sc, iweight
 219
 220       aver=0
 221       if(iweight.eq.0) then
 222          do 10 i=1,nbin
 223  10         aver=max(aver,abs((h1(i)-h2(i))*hnorm(i)))
 224       else if(iweight.eq.1) then
 225          do 20 i=1,nbin
 226  20         aver=aver+abs((h1(i)-h2(i))*hnorm(i))/real(i)
 227       else if(iweight.eq.2) then
 228          do 30 i=1,nbin
 229  30         aver=aver+((h1(i)-h2(i))*hnorm(i))**2/real(i)
 230       endif
 231       end
 232
 233 c-------------------------------------------------------------------
 234 c diagnostic output
 235 c
 236       subroutine dump()
 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
 240
 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(),'()')
 246       end
 247
 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
 252       dimension hcop(*)
 253
 254       call cost_full(0)
 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(),'()')
 260       end