binaries/src/disembl/Tisean_3.0.1/source_f/randomize/cost/spikeauto.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 spike train 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(nhist=100000)
  32       dimension ihist0(nhist), ihist(nhist)
  33       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
  34
  35       iweight=ican('W',0)
  36       bininv=1./fmust("d")
  37       totbin=fmust("D")
  38       nbin=min(int(totbin*bininv)+1,nhist)
  39       inter=lopt("i",1)
  40       ncol=1
  41       end
  42
  43 c-------------------------------------------------------------------
  44 c print version information on cost function
  45 c
  46       subroutine what_cost()
  47       call ptext("Cost function: spike train 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# [-i -W#]")
  55       call popt("d","time span of one bin")
  56       call popt("D","total time spanned")
  57       call popt("i","expect intervals rather than times")
  58       call popt("W",
  59      .   "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 (0)")
  60       end
  61
  62 c-------------------------------------------------------------------
  63 c initialise all that is needed for cost function
  64 c
  65       subroutine cost_init()
  66       parameter(nhist=100000)
  67       dimension ihist0(nhist), ihist(nhist)
  68       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
  69
  70       call sauto(nbin,bininv,ihist0)
  71       end
  72
  73 c-------------------------------------------------------------------
  74 c initial transformation on time series and its inverse
  75 c here: series internally stored as intervals
  76 c
  77       subroutine cost_transform(nmax,mcmax,nxdum,x)
  78       parameter(nx=100000)
  79       parameter(nhist=100000)
  80       dimension ihist0(nhist), ihist(nhist)
  81       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
  82       dimension nxclu(nx)
  83       common /permutecom/ mxclu, nxclu
  84       dimension x(*), lx(nx)
  85
  86       if(inter.eq.1) return
  87       call sort(nmax,x,lx)
  88       do 10 n=nmax,2,-1
  89  10      x(n)=x(n)-x(n-1)
  90       mxclu=mxclu+1
  91       nxclu(mxclu)=1
  92       end
  93
  94       subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
  95       parameter(nhist=100000)
  96       dimension ihist0(nhist), ihist(nhist)
  97       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
  98       dimension x(*), y(*)
  99
 100       do 10 n=1,nmax
 101  10      y(n)=x(n)
 102       if(inter.eq.1) return
 103       do 20 n=2,nmax
 104  20      y(n)=y(n)+y(n-1)
 105       end
 106
 107 c-------------------------------------------------------------------
 108 c compute full cost function from scratch
 109 c
 110       subroutine cost_full(iv)
 111       parameter(nhist=100000)
 112       dimension ihist0(nhist), ihist(nhist)
 113       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
 114       common nmax,cost
 115
 116       call sauto(nbin,bininv,ihist)
 117       cost=aver(ihist0,ihist)
 118       if(iv.ne.0) call dump()
 119       end
 120
 121 c-------------------------------------------------------------------
 122 c compute changed cost function on exchange of n1 and n2
 123 c
 124       subroutine cost_update(nn1,nn2,cmax,iaccept,iv)
 125       parameter(nx=100000)
 126       parameter(nhist=100000)
 127       dimension ihist0(nhist), ihist(nhist)
 128       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
 129       dimension ihcop(nhist), x(nx)
 130       common nmax,cost,temp,cmin,rate,x
 131
 132       n1=min(nn1,nn2)
 133       n2=max(nn1,nn2)
 134       comp=0
 135       iaccept=0
 136       do 10 i=1,nbin
 137  10      ihcop(i)=ihist(i)
 138       dx=0
 139       do 20 nn=n1,1,-1
 140          if(nn.lt.n1) dx=dx+x(nn)
 141          if(int(dx*bininv)+1.gt.nbin) goto 1
 142          dxx=dx
 143          do 30 nnn=n1,n2-1
 144             dxx=dxx+x(nnn)
 145             il=int(dxx*bininv)+1
 146             if(il.gt.nbin) goto 20
 147  30         ihcop(il)=ihcop(il)-1
 148  20      continue
 149  1    dx=0
 150       do 40 nn=n2,1,-1
 151          if(nn.lt.n2) dx=dx+x(nn)
 152          if(int(dx*bininv)+1.gt.nbin) goto 2
 153          dxx=dx
 154          do 50 nnn=n2,nmax
 155             dxx=dxx+x(nnn)
 156             il=int(dxx*bininv)+1
 157             if(il.gt.nbin) goto 40
 158  50         ihcop(il)=ihcop(il)-1
 159  40      continue
 160  2    call exch(n1,n2)
 161       dx=0
 162       do 60 nn=n1,1,-1
 163          if(nn.lt.n1) dx=dx+x(nn)
 164          if(int(dx*bininv)+1.gt.nbin) goto 3
 165          dxx=dx
 166          do 70 nnn=n1,n2-1
 167             dxx=dxx+x(nnn)
 168             il=int(dxx*bininv)+1
 169             if(il.gt.nbin) goto 60
 170  70         ihcop(il)=ihcop(il)+1
 171  60      continue
 172  3    dx=0
 173       do 80 nn=n2,1,-1
 174          if(nn.lt.n2) dx=dx+x(nn)
 175          if(int(dx*bininv)+1.gt.nbin) goto 4
 176          dxx=dx
 177          do 90 nnn=n2,nmax
 178             dxx=dxx+x(nnn)
 179             il=int(dxx*bininv)+1
 180             if(il.gt.nbin) goto 80
 181  90         ihcop(il)=ihcop(il)+1
 182  80      continue
 183  4    comp=aver(ihist0,ihcop)
 184       if(comp.ge.cmax) then
 185          call exch(n1,n2)
 186          return
 187       endif
 188       cost=comp  ! if got here: accept
 189       iaccept=1
 190       if(iv.ne.0) call panic(ihcop)
 191       do 100 i=1,nbin
 192  100     ihist(i)=ihcop(i)
 193       end
 194
 195 c-------------------------------------------------------------------
 196 c compute autocorrealtion from scratch
 197 c
 198       subroutine sauto(nbin,bininv,ihist)
 199       parameter(nx=100000)
 200       dimension ihist(*)
 201       common nmax,cost,temp,cmin,rate,x
 202       dimension x(nx)
 203
 204       do 10 i=1,nbin
 205  10      ihist(i)=0
 206       do 20 n1=1,nmax
 207          dx=0
 208          do 30 n2=n1,nmax
 209             dx=dx+x(n2)
 210             il=int(dx*bininv)+1
 211             if(il.gt.nbin) goto 20
 212  30         ihist(il)=ihist(il)+1
 213  20      continue
 214       end
 215
 216 c-------------------------------------------------------------------
 217 c weighted average of autocorrelation
 218 c
 219       function aver(ih1,ih2)
 220       parameter(nhist=100000)
 221       dimension ih1(nhist), ih2(nhist)
 222       dimension ihist0(nhist), ihist(nhist)
 223       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
 224
 225       aver=0
 226       if(iweight.eq.0) then
 227          do 10 i=1,nbin
 228  10         aver=max(aver,real(abs(ih1(i)-ih2(i))))
 229       else if(iweight.eq.1) then
 230          do 20 i=1,nbin
 231  20         aver=aver+real(abs(ih1(i)-ih2(i)))/real(i)
 232       else if(iweight.eq.2) then
 233          do 30 i=1,nbin
 234  30         aver=aver+(ih1(i)-ih2(i))**2/real(i)
 235       endif
 236       end
 237
 238 c-------------------------------------------------------------------
 239 c diagnostic output
 240 c
 241       subroutine dump()
 242       parameter(nhist=100000)
 243       dimension ihist0(nhist), ihist(nhist)
 244       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
 245
 246       write(istderr(),'(5hgoal ,4i12)') (ihist0(n),n=1,min(4,nbin))
 247       write(istderr(),'(5his   ,4i12)') (ihist(n),n=1,min(4,nbin))
 248       write(istderr(),'(5hmiss ,4i12)')
 249      .   (abs(ihist0(n)-ihist(n)),n=1,min(4,nbin))
 250       write(istderr(),'()')
 251       end
 252
 253       subroutine panic(ihcop)
 254       parameter(nhist=100000)
 255       dimension ihist0(nhist), ihist(nhist)
 256       common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
 257       dimension ihcop(*)
 258
 259       call cost_full(0)
 260       write(istderr(),'(7hupdate ,4i12)') (ihcop(n),n=1,min(4,nbin))
 261       write(istderr(),'(7hfresh  ,4i12)') (ihist(n),n=1,min(4,nbin))
 262       write(istderr(),'(7hdiscr  ,4i12)')
 263      .   (abs(ihcop(n)-ihist(n)),n=1,min(4,nbin))
 264       write(istderr(),'()')
 265       end