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 spike train power spectrum
25 c author T. Schreiber (1999)
27 c-------------------------------------------------------------------
28 c get cost function specific options
30 subroutine opts_cost(ncol)
31 parameter(mfreq=100000)
32 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
33 . sp0(mfreq), sp(mfreq)
34 common /costcom/ nfreq, fmax, inter,
35 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
44 c-------------------------------------------------------------------
45 c print version information on cost function
47 subroutine what_cost()
48 call ptext("Cost function: spike train power spectrum")
51 c-------------------------------------------------------------------
52 c print cost function specific usage message
54 subroutine usage_cost()
55 call ptext("Cost function options: [-F# -## -w# -i]")
57 . "average: 0=max(s) 1=|s|/f 2=(s/f)**2 3=|s| (0)")
58 call popt("F","maximal frequency (2*l / total time)")
59 call popt("#","number of frequencies (F* total time /2)")
60 call popt("w","frequency resolution (0)")
61 call popt("i","expect intervals rather than times")
64 c-------------------------------------------------------------------
65 c initialise all that is needed for cost function
67 subroutine cost_init()
68 parameter(nx=100000,mfreq=100000)
69 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
70 . sp0(mfreq), sp(mfreq)
71 common /costcom/ nfreq, fmax, inter,
72 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
74 common nmax,cost,temp,cmin,rate,x
76 if(fmax.le.0.) fmax=2*nmax/(x(nmax)-x(1))
77 if(nfreq.le.0) nfreq=fmax*(x(nmax)-x(1))/2
78 if(nfreq.gt.mfreq) write(istderr(),'(a)')
79 . "truncated to ", mfreq," frequencies"
80 nfreq=min(mfreq,nfreq)
81 write(istderr(),*) "randomize_spikespec: total time covered: ",
83 write(istderr(),*) "randomize_spikespec: computing ", nfreq,
84 . " frequencies up to ", fmax
85 call sspect(nfreq,fmax/nfreq,sp0r,sp0i,sp0)
88 c-------------------------------------------------------------------
89 c initial transformation on time series and its inverse
91 subroutine cost_transform(nmax,mcmax,nxdum,x)
92 parameter(mfreq=100000,nx=100000)
93 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
94 . sp0(mfreq), sp(mfreq)
95 common /costcom/ nfreq, fmax, inter,
96 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
97 dimension x(nx), lx(nx)
102 1 call sort(nmax,x,lx)
105 subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
106 parameter(mfreq=100000)
107 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
108 . sp0(mfreq), sp(mfreq)
109 common /costcom/ nfreq, fmax, inter,
110 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
111 dimension x(nmax), y(nmax)
115 if(inter.ne.1) return
120 c-------------------------------------------------------------------
121 c compute full cost function from scratch
123 subroutine cost_full(iv)
124 parameter(mfreq=100000)
125 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
126 . sp0(mfreq), sp(mfreq)
127 common /costcom/ nfreq, fmax, inter,
128 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
131 call sspect(nfreq,fmax/nfreq,spr,spi,sp)
135 10 call aver(cc,sp(i)-sp0(i),i)
137 if(iv.ne.0) call dump()
140 c-------------------------------------------------------------------
141 c compute changed cost function on exchange of n1 and n2
143 subroutine cost_update(n1,n2,cmax,iaccept,iv)
144 parameter(mfreq=100000,nx=100000)
145 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
146 . sp0(mfreq), sp(mfreq)
147 common /costcom/ nfreq, fmax, inter,
148 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
149 dimension sprcop(mfreq), spicop(mfreq), spcop(mfreq), x(nx)
150 common nmax,cost,temp,cmin,rate,x
158 xx=x(n1-1)+x(n1+1)-x(n1)
159 sprcop(i)=spr(i)-cos(omega*x(n1))+cos(omega*xx)
160 spicop(i)=spi(i)-sin(omega*x(n1))+sin(omega*xx)
161 spcop(i)=sprcop(i)**2+spicop(i)**2
162 call aver(comp,sp0(i)-spcop(i),i)
163 10 if(comp.ge.cmax) return
164 cost=comp ! if got here: accept
167 if(iv.ne.0) call panic(spcop)
174 c-------------------------------------------------------------------
175 c compute spectrum from scratch
177 subroutine sspect(nfreq,fres,spr,spi,sp)
179 dimension spr(*), spi(*), sp(*), x(nx)
180 common nmax,cost,temp,cmin,rate,x
189 sr=sr+cos(omega*x(n))
190 20 si=si+sin(omega*x(n))
196 c-------------------------------------------------------------------
197 c weighted average of autocorrelation
199 subroutine aver(cav,dc,n)
200 parameter(mfreq=100000)
201 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
202 . sp0(mfreq), sp(mfreq)
203 common /costcom/ nfreq, fmax, inter,
204 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
206 if(iweight.eq.0) then
207 cav=max(cav,abs(dc)) ! max (L-infinity) norm
208 else if(iweight.eq.1) then
209 cav=cav+abs(dc)/n ! L-1 norm (sum of moduli), weighted by freq.
210 else if(iweight.eq.2) then
211 cav=cav+(dc/n)**2 ! L-2 norm (sum of squares), weighted by freq.
212 else if(iweight.eq.2) then
213 cav=cav+abs(dc) ! L-1 norm (sum of moduli)
217 c-------------------------------------------------------------------
221 parameter(mfreq=100000)
222 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
223 . sp0(mfreq), sp(mfreq)
224 common /costcom/ nfreq, fmax, inter,
225 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
228 write(istderr(),'(5hgoal ,4f9.3)') (sp0(n),n=1,min(4,nfreq))
229 write(istderr(),'(5his ,4f9.3)') (sp(n),n=1,min(4,nfreq))
230 write(istderr(),'(5hmiss ,4f9.3)')
231 . ((sp0(n)-sp(n)),n=1,min(4,nfreq))
232 write(istderr(),'()')
235 subroutine panic(spcop)
236 parameter(mfreq=100000)
238 dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
239 . sp0(mfreq), sp(mfreq)
240 common /costcom/ nfreq, fmax, inter,
241 . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
245 write(istderr(),'(7hupdate ,4f9.3)')
246 . (spcop(n),n=1,min(4,nfreq))
247 write(istderr(),'(7hfresh ,4f9.3)') (sp(n),n=1,min(4,nfreq))
248 write(istderr(),'(7hdiscr ,4f9.3)')
249 . ((spcop(n)-sp(n)),n=1,min(4,nfreq))
250 write(istderr(),'()')