--- /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 spike train power spectrum
+c author T. Schreiber (1999)
+c
+c-------------------------------------------------------------------
+c get cost function specific options
+c
+ subroutine opts_cost(ncol)
+ parameter(mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+
+ iweight=ican('W',0)
+ fmax=fcan("F",0)
+ nfreq=ican("#",0)
+ 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 power spectrum")
+ end
+
+c-------------------------------------------------------------------
+c print cost function specific usage message
+c
+ subroutine usage_cost()
+ call ptext("Cost function options: [-F# -## -w# -i]")
+ call popt("W",
+ . "average: 0=max(s) 1=|s|/f 2=(s/f)**2 3=|s| (0)")
+ call popt("F","maximal frequency (2*l / total time)")
+ call popt("#","number of frequencies (F* total time /2)")
+ call popt("w","frequency resolution (0)")
+ call popt("i","expect intervals rather than times")
+ end
+
+c-------------------------------------------------------------------
+c initialise all that is needed for cost function
+c
+ subroutine cost_init()
+ parameter(nx=100000,mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ dimension x(nx)
+ common nmax,cost,temp,cmin,rate,x
+
+ if(fmax.le.0.) fmax=2*nmax/(x(nmax)-x(1))
+ if(nfreq.le.0) nfreq=fmax*(x(nmax)-x(1))/2
+ if(nfreq.gt.mfreq) write(istderr(),'(a)')
+ . "truncated to ", mfreq," frequencies"
+ nfreq=min(mfreq,nfreq)
+ write(istderr(),*) "randomize_spikespec: total time covered: ",
+ . x(nmax)-x(1)
+ write(istderr(),*) "randomize_spikespec: computing ", nfreq,
+ . " frequencies up to ", fmax
+ call sspect(nfreq,fmax/nfreq,sp0r,sp0i,sp0)
+ end
+
+c-------------------------------------------------------------------
+c initial transformation on time series and its inverse
+c
+ subroutine cost_transform(nmax,mcmax,nxdum,x)
+ parameter(mfreq=100000,nx=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ dimension x(nx), lx(nx)
+
+ if(inter.eq.0) goto 1
+ do 10 n=2,nmax
+ 10 x(n)=x(n)+x(n-1)
+ 1 call sort(nmax,x,lx)
+ end
+
+ subroutine cost_inverse(nmax,mcmax,nxdum,x,y)
+ parameter(mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ dimension x(nmax), y(nmax)
+
+ do 10 n=1,nmax
+ 10 y(n)=x(n)
+ if(inter.ne.1) return
+ do 20 n=nmax,2,-1
+ 20 y(n)=y(n)-y(n-1)
+ end
+
+c-------------------------------------------------------------------
+c compute full cost function from scratch
+c
+ subroutine cost_full(iv)
+ parameter(mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ common nmax,cost
+
+ call sspect(nfreq,fmax/nfreq,spr,spi,sp)
+
+ cc=0
+ do 10 i=1,nfreq
+ 10 call aver(cc,sp(i)-sp0(i),i)
+ cost=cc
+ if(iv.ne.0) call dump()
+ end
+
+c-------------------------------------------------------------------
+c compute changed cost function on exchange of n1 and n2
+c
+ subroutine cost_update(n1,n2,cmax,iaccept,iv)
+ parameter(mfreq=100000,nx=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ dimension sprcop(mfreq), spicop(mfreq), spcop(mfreq), x(nx)
+ common nmax,cost,temp,cmin,rate,x
+ data pi/3.1415926/
+
+ comp=0
+ iaccept=0
+ do 10 i=1,nfreq
+ f=i*(fmax/nfreq)
+ omega=2*pi*f
+ xx=x(n1-1)+x(n1+1)-x(n1)
+ sprcop(i)=spr(i)-cos(omega*x(n1))+cos(omega*xx)
+ spicop(i)=spi(i)-sin(omega*x(n1))+sin(omega*xx)
+ spcop(i)=sprcop(i)**2+spicop(i)**2
+ call aver(comp,sp0(i)-spcop(i),i)
+ 10 if(comp.ge.cmax) return
+ cost=comp ! if got here: accept
+ iaccept=1
+ call exch(n1,n2)
+ if(iv.ne.0) call panic(spcop)
+ do 20 i=1,nfreq
+ spr(i)=sprcop(i)
+ spi(i)=spicop(i)
+ 20 sp(i)=spcop(i)
+ end
+
+c-------------------------------------------------------------------
+c compute spectrum from scratch
+c
+ subroutine sspect(nfreq,fres,spr,spi,sp)
+ parameter(nx=100000)
+ dimension spr(*), spi(*), sp(*), x(nx)
+ common nmax,cost,temp,cmin,rate,x
+ data pi/3.1415926/
+
+ do 10 i=1,nfreq
+ f=i*fres
+ omega=2*pi*f
+ sr=0
+ si=0
+ do 20 n=1,nmax
+ sr=sr+cos(omega*x(n))
+ 20 si=si+sin(omega*x(n))
+ spr(i)=sr
+ spi(i)=si
+ 10 sp(i)=sr**2+si**2
+ end
+
+c-------------------------------------------------------------------
+c weighted average of autocorrelation
+c
+ subroutine aver(cav,dc,n)
+ parameter(mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+
+ if(iweight.eq.0) then
+ cav=max(cav,abs(dc)) ! max (L-infinity) norm
+ else if(iweight.eq.1) then
+ cav=cav+abs(dc)/n ! L-1 norm (sum of moduli), weighted by freq.
+ else if(iweight.eq.2) then
+ cav=cav+(dc/n)**2 ! L-2 norm (sum of squares), weighted by freq.
+ else if(iweight.eq.2) then
+ cav=cav+abs(dc) ! L-1 norm (sum of moduli)
+ endif
+ end
+
+c-------------------------------------------------------------------
+c diagnostic output
+c
+ subroutine dump()
+ parameter(mfreq=100000)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ common nmax
+
+ write(istderr(),'(5hgoal ,4f9.3)') (sp0(n),n=1,min(4,nfreq))
+ write(istderr(),'(5his ,4f9.3)') (sp(n),n=1,min(4,nfreq))
+ write(istderr(),'(5hmiss ,4f9.3)')
+ . ((sp0(n)-sp(n)),n=1,min(4,nfreq))
+ write(istderr(),'()')
+ end
+
+ subroutine panic(spcop)
+ parameter(mfreq=100000)
+ dimension spcop(*)
+ dimension sp0r(mfreq), sp0i(mfreq), spr(mfreq), spi(mfreq),
+ . sp0(mfreq), sp(mfreq)
+ common /costcom/ nfreq, fmax, inter,
+ . sp0r, sp0i, sp0, spr, spi, sp, sd, sc, iweight
+ common nmax
+
+ call cost_full(0)
+ write(istderr(),'(7hupdate ,4f9.3)')
+ . (spcop(n),n=1,min(4,nfreq))
+ write(istderr(),'(7hfresh ,4f9.3)') (sp(n),n=1,min(4,nfreq))
+ write(istderr(),'(7hdiscr ,4f9.3)')
+ . ((spcop(n)-sp(n)),n=1,min(4,nfreq))
+ write(istderr(),'()')
+ end
git://source.jalview.org / jabaws.git / blobdiff