website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/randomize/cost/spikespec.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   spike train power spectrum
  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(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
  36
  37       iweight=ican('W',0)
  38       fmax=fcan("F",0)
  39       nfreq=ican("#",0)
  40       inter=lopt("i",1)
  41       ncol=1
  42       end
  43
  44 c-------------------------------------------------------------------
  45 c print version information on cost function
  46 c
  47       subroutine what_cost()
  48       call ptext("Cost function: spike train power spectrum")
  49       end
  50
  51 c-------------------------------------------------------------------
  52 c print cost function specific usage message
  53 c
  54       subroutine usage_cost()
  55       call ptext("Cost function options: [-F# -## -w# -i]")
  56       call popt("W",
  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")
  62       end
  63
  64 c-------------------------------------------------------------------
  65 c initialise all that is needed for cost function
  66 c
  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
  73       dimension x(nx)
  74       common nmax,cost,temp,cmin,rate,x
  75
  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: ",
  82      .   x(nmax)-x(1)
  83       write(istderr(),*) "randomize_spikespec: computing ", nfreq,
  84      .   " frequencies up to ", fmax
  85       call sspect(nfreq,fmax/nfreq,sp0r,sp0i,sp0)
  86       end
  87
  88 c-------------------------------------------------------------------
  89 c initial transformation on time series and its inverse
  90 c
  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)
  98
  99       if(inter.eq.0) goto 1
 100       do 10 n=2,nmax
 101  10      x(n)=x(n)+x(n-1)
 102  1    call sort(nmax,x,lx)
 103       end
 104
 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)
 112
 113       do 10 n=1,nmax
 114  10      y(n)=x(n)
 115       if(inter.ne.1) return
 116       do 20 n=nmax,2,-1
 117  20      y(n)=y(n)-y(n-1)
 118       end
 119
 120 c-------------------------------------------------------------------
 121 c compute full cost function from scratch
 122 c
 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
 129       common nmax,cost
 130
 131       call sspect(nfreq,fmax/nfreq,spr,spi,sp)
 132
 133       cc=0
 134       do 10 i=1,nfreq
 135  10      call aver(cc,sp(i)-sp0(i),i)
 136       cost=cc
 137       if(iv.ne.0) call dump()
 138       end
 139
 140 c-------------------------------------------------------------------
 141 c compute changed cost function on exchange of n1 and n2
 142 c
 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
 151       data pi/3.1415926/
 152
 153       comp=0
 154       iaccept=0
 155       do 10 i=1,nfreq
 156          f=i*(fmax/nfreq)
 157          omega=2*pi*f
 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
 165       iaccept=1
 166       call exch(n1,n2)
 167       if(iv.ne.0) call panic(spcop)
 168       do 20 i=1,nfreq
 169          spr(i)=sprcop(i)
 170          spi(i)=spicop(i)
 171  20      sp(i)=spcop(i)
 172       end
 173
 174 c-------------------------------------------------------------------
 175 c compute spectrum from scratch
 176 c
 177       subroutine sspect(nfreq,fres,spr,spi,sp)
 178       parameter(nx=100000)
 179       dimension spr(*), spi(*), sp(*), x(nx)
 180       common nmax,cost,temp,cmin,rate,x
 181       data pi/3.1415926/
 182
 183       do 10 i=1,nfreq
 184          f=i*fres
 185          omega=2*pi*f
 186          sr=0
 187          si=0
 188          do 20 n=1,nmax
 189             sr=sr+cos(omega*x(n))
 190  20         si=si+sin(omega*x(n))
 191          spr(i)=sr
 192          spi(i)=si
 193  10      sp(i)=sr**2+si**2
 194       end
 195
 196 c-------------------------------------------------------------------
 197 c weighted average of autocorrelation
 198 c
 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
 205
 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)
 214       endif
 215       end
 216
 217 c-------------------------------------------------------------------
 218 c diagnostic output
 219 c
 220       subroutine dump()
 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
 226       common nmax
 227
 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(),'()')
 233       end
 234
 235       subroutine panic(spcop)
 236       parameter(mfreq=100000)
 237       dimension spcop(*)
 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
 242       common nmax
 243
 244       call cost_full(0)
 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(),'()')
 251       end