--- /dev/null
+c===========================================================================
+c
+c This file is part of TISEAN
+c
+c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
+c
+c TISEAN is free software; you can redistribute it and/or modify
+c it under the terms of the GNU General Public License as published by
+c the Free Software Foundation; either version 2 of the License, or
+c (at your option) any later version.
+c
+c TISEAN is distributed in the hope that it will be useful,
+c but WITHOUT ANY WARRANTY; without even the implied warranty of
+c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+c GNU General Public License for more details.
+c
+c You should have received a copy of the GNU General Public License
+c along with TISEAN; if not, write to the Free Software
+c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+c
+c===========================================================================
+c part of the TISEAN randomize package for constraint surrogates
+c cost function
+c binned spike train autocorrelation function
+c author T. Schreiber (1999)
+c
+c-------------------------------------------------------------------
+c get cost function specific options
+c
+ subroutine opts_cost(ncol)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+
+ iweight=ican('W',0)
+ bininv=1./fmust("d")
+ totbin=fmust("D")
+ nbin=min(int(totbin*bininv)+1,nhist)
+ inter=lopt("i",1)
+ ncol=1
+ end
+
+c-------------------------------------------------------------------
+c print version information on cost function
+c
+ subroutine what_cost()
+ call ptext("Cost function: spike train autocorrelation function")
+ end
+
+c-------------------------------------------------------------------
+c print cost function specific usage message
+c
+ subroutine usage_cost()
+ call ptext("Cost function options: -d# -D# [-i -W#]")
+ call popt("d","time span of one bin")
+ call popt("D","total time spanned")
+ call popt("i","expect intervals rather than times")
+ call popt("W",
+ . "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 (0)")
+ end
+
+c-------------------------------------------------------------------
+c initialise all that is needed for cost function
+c
+ subroutine cost_init()
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+
+ call sauto(nbin,bininv,ihist0)
+ end
+
+c-------------------------------------------------------------------
+c initial transformation on time series and its inverse
+c here: series internally stored as intervals
+c
+ subroutine cost_transform(nmax,mcmax,nxdum,x)
+ parameter(nx=100000)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+ dimension nxclu(nx)
+ common /permutecom/ mxclu, nxclu
+ dimension x(*), lx(nx)
+
+ if(inter.eq.1) return
+ call sort(nmax,x,lx)
+ do 10 n=nmax,2,-1
+ 10 x(n)=x(n)-x(n-1)
+ mxclu=mxclu+1
+ nxclu(mxclu)=1
+ end
+
+ subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+ dimension x(*), y(*)
+
+ do 10 n=1,nmax
+ 10 y(n)=x(n)
+ if(inter.eq.1) return
+ do 20 n=2,nmax
+ 20 y(n)=y(n)+y(n-1)
+ end
+
+c-------------------------------------------------------------------
+c compute full cost function from scratch
+c
+ subroutine cost_full(iv)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+ common nmax,cost
+
+ call sauto(nbin,bininv,ihist)
+ cost=aver(ihist0,ihist)
+ if(iv.ne.0) call dump()
+ end
+
+c-------------------------------------------------------------------
+c compute changed cost function on exchange of n1 and n2
+c
+ subroutine cost_update(nn1,nn2,cmax,iaccept,iv)
+ parameter(nx=100000)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+ dimension ihcop(nhist), x(nx)
+ common nmax,cost,temp,cmin,rate,x
+
+ n1=min(nn1,nn2)
+ n2=max(nn1,nn2)
+ comp=0
+ iaccept=0
+ do 10 i=1,nbin
+ 10 ihcop(i)=ihist(i)
+ dx=0
+ do 20 nn=n1,1,-1
+ if(nn.lt.n1) dx=dx+x(nn)
+ if(int(dx*bininv)+1.gt.nbin) goto 1
+ dxx=dx
+ do 30 nnn=n1,n2-1
+ dxx=dxx+x(nnn)
+ il=int(dxx*bininv)+1
+ if(il.gt.nbin) goto 20
+ 30 ihcop(il)=ihcop(il)-1
+ 20 continue
+ 1 dx=0
+ do 40 nn=n2,1,-1
+ if(nn.lt.n2) dx=dx+x(nn)
+ if(int(dx*bininv)+1.gt.nbin) goto 2
+ dxx=dx
+ do 50 nnn=n2,nmax
+ dxx=dxx+x(nnn)
+ il=int(dxx*bininv)+1
+ if(il.gt.nbin) goto 40
+ 50 ihcop(il)=ihcop(il)-1
+ 40 continue
+ 2 call exch(n1,n2)
+ dx=0
+ do 60 nn=n1,1,-1
+ if(nn.lt.n1) dx=dx+x(nn)
+ if(int(dx*bininv)+1.gt.nbin) goto 3
+ dxx=dx
+ do 70 nnn=n1,n2-1
+ dxx=dxx+x(nnn)
+ il=int(dxx*bininv)+1
+ if(il.gt.nbin) goto 60
+ 70 ihcop(il)=ihcop(il)+1
+ 60 continue
+ 3 dx=0
+ do 80 nn=n2,1,-1
+ if(nn.lt.n2) dx=dx+x(nn)
+ if(int(dx*bininv)+1.gt.nbin) goto 4
+ dxx=dx
+ do 90 nnn=n2,nmax
+ dxx=dxx+x(nnn)
+ il=int(dxx*bininv)+1
+ if(il.gt.nbin) goto 80
+ 90 ihcop(il)=ihcop(il)+1
+ 80 continue
+ 4 comp=aver(ihist0,ihcop)
+ if(comp.ge.cmax) then
+ call exch(n1,n2)
+ return
+ endif
+ cost=comp ! if got here: accept
+ iaccept=1
+ if(iv.ne.0) call panic(ihcop)
+ do 100 i=1,nbin
+ 100 ihist(i)=ihcop(i)
+ end
+
+c-------------------------------------------------------------------
+c compute autocorrealtion from scratch
+c
+ subroutine sauto(nbin,bininv,ihist)
+ parameter(nx=100000)
+ dimension ihist(*)
+ common nmax,cost,temp,cmin,rate,x
+ dimension x(nx)
+
+ do 10 i=1,nbin
+ 10 ihist(i)=0
+ do 20 n1=1,nmax
+ dx=0
+ do 30 n2=n1,nmax
+ dx=dx+x(n2)
+ il=int(dx*bininv)+1
+ if(il.gt.nbin) goto 20
+ 30 ihist(il)=ihist(il)+1
+ 20 continue
+ end
+
+c-------------------------------------------------------------------
+c weighted average of autocorrelation
+c
+ function aver(ih1,ih2)
+ parameter(nhist=100000)
+ dimension ih1(nhist), ih2(nhist)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+
+ aver=0
+ if(iweight.eq.0) then
+ do 10 i=1,nbin
+ 10 aver=max(aver,real(abs(ih1(i)-ih2(i))))
+ else if(iweight.eq.1) then
+ do 20 i=1,nbin
+ 20 aver=aver+real(abs(ih1(i)-ih2(i)))/real(i)
+ else if(iweight.eq.2) then
+ do 30 i=1,nbin
+ 30 aver=aver+(ih1(i)-ih2(i))**2/real(i)
+ endif
+ end
+
+c-------------------------------------------------------------------
+c diagnostic output
+c
+ subroutine dump()
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+
+ write(istderr(),'(5hgoal ,4i12)') (ihist0(n),n=1,min(4,nbin))
+ write(istderr(),'(5his ,4i12)') (ihist(n),n=1,min(4,nbin))
+ write(istderr(),'(5hmiss ,4i12)')
+ . (abs(ihist0(n)-ihist(n)),n=1,min(4,nbin))
+ write(istderr(),'()')
+ end
+
+ subroutine panic(ihcop)
+ parameter(nhist=100000)
+ dimension ihist0(nhist), ihist(nhist)
+ common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
+ dimension ihcop(*)
+
+ call cost_full(0)
+ write(istderr(),'(7hupdate ,4i12)') (ihcop(n),n=1,min(4,nbin))
+ write(istderr(),'(7hfresh ,4i12)') (ihist(n),n=1,min(4,nbin))
+ write(istderr(),'(7hdiscr ,4i12)')
+ . (abs(ihcop(n)-ihist(n)),n=1,min(4,nbin))
+ write(istderr(),'()')
+ end
git://source.jalview.org / jabaws.git / blobdiff