1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
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.
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.
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
21 c===========================================================================
22 c part of the TISEAN randomize package for constraint surrogates
24 c binned spike train autocorrelation function
25 c author T. Schreiber (1999)
27 c-------------------------------------------------------------------
28 c get cost function specific options
30 subroutine opts_cost(ncol)
31 parameter(nhist=100000)
32 dimension ihist0(nhist), ihist(nhist)
33 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
38 nbin=min(int(totbin*bininv)+1,nhist)
43 c-------------------------------------------------------------------
44 c print version information on cost function
46 subroutine what_cost()
47 call ptext("Cost function: spike train autocorrelation function")
50 c-------------------------------------------------------------------
51 c print cost function specific usage message
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")
59 . "average: 0=max(c) 1=|c|/lag 2=(c/lag)**2 (0)")
62 c-------------------------------------------------------------------
63 c initialise all that is needed for cost function
65 subroutine cost_init()
66 parameter(nhist=100000)
67 dimension ihist0(nhist), ihist(nhist)
68 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
70 call sauto(nbin,bininv,ihist0)
73 c-------------------------------------------------------------------
74 c initial transformation on time series and its inverse
75 c here: series internally stored as intervals
77 subroutine cost_transform(nmax,mcmax,nxdum,x)
79 parameter(nhist=100000)
80 dimension ihist0(nhist), ihist(nhist)
81 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
83 common /permutecom/ mxclu, nxclu
84 dimension x(*), lx(nx)
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
102 if(inter.eq.1) return
107 c-------------------------------------------------------------------
108 c compute full cost function from scratch
110 subroutine cost_full(iv)
111 parameter(nhist=100000)
112 dimension ihist0(nhist), ihist(nhist)
113 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
116 call sauto(nbin,bininv,ihist)
117 cost=aver(ihist0,ihist)
118 if(iv.ne.0) call dump()
121 c-------------------------------------------------------------------
122 c compute changed cost function on exchange of n1 and n2
124 subroutine cost_update(nn1,nn2,cmax,iaccept,iv)
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
140 if(nn.lt.n1) dx=dx+x(nn)
141 if(int(dx*bininv)+1.gt.nbin) goto 1
146 if(il.gt.nbin) goto 20
147 30 ihcop(il)=ihcop(il)-1
151 if(nn.lt.n2) dx=dx+x(nn)
152 if(int(dx*bininv)+1.gt.nbin) goto 2
157 if(il.gt.nbin) goto 40
158 50 ihcop(il)=ihcop(il)-1
163 if(nn.lt.n1) dx=dx+x(nn)
164 if(int(dx*bininv)+1.gt.nbin) goto 3
169 if(il.gt.nbin) goto 60
170 70 ihcop(il)=ihcop(il)+1
174 if(nn.lt.n2) dx=dx+x(nn)
175 if(int(dx*bininv)+1.gt.nbin) goto 4
180 if(il.gt.nbin) goto 80
181 90 ihcop(il)=ihcop(il)+1
183 4 comp=aver(ihist0,ihcop)
184 if(comp.ge.cmax) then
188 cost=comp ! if got here: accept
190 if(iv.ne.0) call panic(ihcop)
192 100 ihist(i)=ihcop(i)
195 c-------------------------------------------------------------------
196 c compute autocorrealtion from scratch
198 subroutine sauto(nbin,bininv,ihist)
201 common nmax,cost,temp,cmin,rate,x
211 if(il.gt.nbin) goto 20
212 30 ihist(il)=ihist(il)+1
216 c-------------------------------------------------------------------
217 c weighted average of autocorrelation
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
226 if(iweight.eq.0) then
228 10 aver=max(aver,real(abs(ih1(i)-ih2(i))))
229 else if(iweight.eq.1) then
231 20 aver=aver+real(abs(ih1(i)-ih2(i)))/real(i)
232 else if(iweight.eq.2) then
234 30 aver=aver+(ih1(i)-ih2(i))**2/real(i)
238 c-------------------------------------------------------------------
242 parameter(nhist=100000)
243 dimension ihist0(nhist), ihist(nhist)
244 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
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(),'()')
253 subroutine panic(ihcop)
254 parameter(nhist=100000)
255 dimension ihist0(nhist), ihist(nhist)
256 common /costcom/ inter, bininv, nbin, ihist0, ihist, iweight
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(),'()')