website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/xc2.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
  23 c   cross-correlation integral xc2
  24 c   see  H. Kantz, Phys.Rev.E49, 5091 (1994)
  25 c
  26 c   authors: T. Schreiber & H. Kantz (1998)
  27 c   multivariate version: H. Kantz (Jan 2007)
  28 c
  29 c=========================================================================
  30 c
  31       parameter(nx=1000000,me=30,meps=1000,mx=10)
  32       dimension x(nx,mx), c(me,meps), eps(meps), mdeps(meps), icol(mx)
  33       dimension y(nx,mx)
  34       integer mlist(2)
  35       character*72 file1, file2, fout
  36       data ipmin/1000/, res/2./, eps0/1e-30/, epsm/1e30/, id0/1/
  37       data iverb/1/
  38
  39 c=======================================================================
  40 c assume two input files (file1, file2)
  41 c - with identical structure in terms of colums
  42 c - with identical embedding spaces
  43 c - with individual length and exclusions: l, L, x, X
  44 c - multivariate data: maximum length nx, maximum dimension mx
  45 c
  46 c   norm: max-norm
  47 c
  48 c   no rescaling, since xc2 does not make sense if datasets are
  49 c                           not used in their original scalings
  50 c=================================================================
  51
  52       call whatido("cross correlation sum of two data sets",iverb)
  53       id=ican("d",id0)
  54       mmax=2
  55       mdim=1
  56
  57       call imcan("M",2,mc,mlist)
  58       if (mc.ge.2) then
  59        mmax=mlist(2)
  60        mdim=mlist(1)
  61        if (mmax*mdim.lt.2) stop 'Increase embedding dimension'
  62        if (mc.gt.2) print*, 'extra arguments of -m ignored'
  63       endif
  64
  65       ntmin=0
  66       ncmin=ican("n",1000)
  67       ipmin=ican("N",ipmin)
  68       res=fcan("#",res)
  69       feps=2**(1./res)
  70       eps0=fcan("r",eps0)
  71       epsm=fcan("R",epsm)
  72       nmaxx=ican("l",nx)
  73       nmaxy=ican("L",nx)
  74       nexcl1=ican("x",0)
  75       nexcl2=ican("X",0)
  76
  77       call columns(mc,mx,icol)
  78
  79       if (mc.gt.0.and.mc.ne.mdim) stop 'improper number of columns'
  80       isout=igetout(fout,0)
  81       if(fout.eq." ") isout=1
  82
  83       call nthstring(1,file1)
  84       if(file1.eq."-") stop "first input file name missing"
  85       call xreadfile(nmaxx,mdim,nx,x,nexcl1,icol,file1,iverb)
  86       call nthstring(2,file2)
  87       if(file2.eq."-") stop "second input file name missing"
  88       call xreadfile(nmaxy,mdim,nx,y,nexcl2,icol,file2,iverb)
  89
  90       if(isout.eq.1) then
  91         call addsuff(fout,file1,"_")
  92         call addsuff(fout,fout,file2)
  93         call addsuff(fout,fout,"_xc2")
  94       endif
  95
  96       epsmax=0.
  97       do imx=1,mmax
  98        call minmax(nmaxx,x(1,imx),xmin,xmax)
  99        epsmax=1.001*max(xmax-xmin,epsmax)
 100       enddo
 101       do imx=1,mmax
 102        call minmax(nmaxy,y(1,imx),xmin,xmax)
 103        epsmax=1.001*max(xmax-xmin,epsmax)
 104       enddo
 105
 106       neps=0
 107
 108       do 10 epsl=log(min(epsm,epsmax)),log(eps0),-log(feps)
 109          neps=neps+1
 110          if(neps.gt.meps) stop "xc2: make meps larger"
 111          eps(neps)=exp(epsl)
 112          do 20 m=1,mmax*mdim
 113  20         c(m,neps)=0
 114          if (mdim.eq.1) then
 115           call crosscor(nmaxx,x,nmaxy,y,eps(neps)
 116      .                 ,id,mmax,c(1,neps),ncmin,ipmin)
 117          else
 118           call mcrosscor(nmaxx,x,nmaxy,y,eps(neps),
 119      .      id,mmax,mdim,c(1,neps),ncmin,ipmin)
 120          endif
 121          mdd=mmax*mdim
 122          mdd1=max(2,mdim)
 123          do 30 m=mdd1,mdd
 124  30         if(c(m,neps).eq.0.) goto 1
 125          m=mdd+1
 126  1       mdd=m-1
 127          if(mdd.eq.mdim-1) stop
 128          mdeps(neps)=mdd
 129          call outfile(fout,iunit,iverb)
 130          do 40 m=mdd1,mdeps(1)
 131             write(iunit,'(4h#m= ,i5)') m
 132             do 50 nn=1,neps
 133                if(mdeps(nn).lt.m) goto 2
 134  50            write(iunit,*) eps(nn), c(m,nn)
 135  2          write(iunit,'()')
 136  40         write(iunit,'()')
 137          close(iunit)
 138  10      write(istderr(),*) eps(neps), mdd, c(mdd,neps)
 139       stop
 140       end
 141 c>--------------------------------------------------------------------
 142       subroutine usage()
 143 c usage message
 144
 145       call whatineed(
 146      .   "-M#,# [-d# -n# -N# -## -r# -R#"//
 147      .   " -o outfile -l# -x# -L# -X# -c#[,#] -V# -h] file1 file2")
 148       call popt("M",
 149      ."# of components, maximal embedding dimension [1,2]")
 150       call popt("d","delay [1]")
 151       call popt("n","minimal number of center points [1000]")
 152       call popt("N","maximal number of pairs [1000]")
 153       call popt("#","resolution, values per octave [2]")
 154       call popt("r",
 155      .   "minimal scale to be probed (as long as pairs found)")
 156       call popt("R","maximal scale to be probed [xmax-xmin]")
 157       call popt("l","length of time series 1 to be read [all data]")
 158       call popt("x","# of initial lines of 1 to be skipped [0]")
 159       call popt("L","length of time series 2 to be read [all data]")
 160       call popt("X","# of initial lines of 2 to be skipped [0]")
 161       call popt("c",
 162      ."columns to be read [1,2,3,.., # of components]")
 163       call pout("file1_file2_xc2")
 164       call pall()
 165       stop
 166       end
 167 c>--------------------------------------------------------------------
 168       subroutine crosscor(nmaxx,x,nmaxy,y,eps,id,m,c,ncmin,ipmin)
 169       parameter(im=100,ii=100000000,nx=1000000,mm=30)
 170       dimension y(nmaxy),x(nmaxx)
 171       dimension jh(0:im*im),ipairs(mm),c(m),jpntr(nx),nlist(nx)
 172
 173       if(nmaxx.gt.nx.or.m.gt.mm) stop "crosscor: make mm/nx larger."
 174       if(nmaxy.gt.nx.or.m.gt.mm) stop "crosscor: make mm/nx larger."
 175
 176       do 10 i=1,m-1
 177  10      ipairs(i)=0
 178       mb=min(m,2)
 179       call base(nmaxx,x,id,mb,jh,jpntr,eps)
 180       do 20 n=(m-1)*id+1,nmaxy
 181          call neigh(nx,y,x,n,nmaxx,id,mb,jh,jpntr,eps,nlist,nfound)
 182          do 30 nn=1,nfound                   ! all neighbours in two dimensions
 183             np=nlist(nn)
 184             if(np.lt.(m-1)*id+1) goto 30
 185             ipairs(1)=ipairs(1)+1
 186             do 40 i=mb,m-1
 187                if(abs(y(n-i*id)-x(np-i*id)).ge.eps) goto 30
 188  40            ipairs(i)=ipairs(i)+1            ! neighbours in 3..m dimensions
 189  30         continue
 190  20      if(n-(m-1)*id.ge.ncmin.and.ipairs(m-1).ge.ipmin) goto 1
 191       n=n-1
 192  1    s=real(n-(m-1)*id)*real(nmaxx-(m-1)*id) ! normalisation
 193       do 50 i=1,m-1
 194  50      if(s.gt.0.) c(i+1)=ipairs(i)/s
 195       end
 196 c>--------------------------------------------------------------------
 197       subroutine mcrosscor(nmaxx,x,nmaxy,y,eps,id,m,mdim,c,ncmin,ipmin)
 198
 199       parameter(im=100,nx=1000000,mm=30,mx=10)
 200
 201       dimension y(nx,mx),x(nx,mx)
 202       dimension jh(0:im*im),ipairs(mm),c(m*mdim),jpntr(nx),nlist(nx)
 203       dimension vx(mm)
 204
 205       if(nmaxx.gt.nx.or.m.gt.mm) stop "mcrosscor: make mm/nx larger."
 206       if(nmaxy.gt.nx.or.m.gt.mm) stop "mcrosscor: make mm/nx larger."
 207
 208       if (m*mdim.gt.mm) stop 'embedding x spatial dimension < 30 !'
 209       do 10 i=1,m*mdim
 210  10      ipairs(i)=0
 211
 212       call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
 213       do 20 n=(m-1)*id+1,nmaxx
 214          do ii=1,mdim
 215           vx(ii)=x(n,ii)
 216          enddo
 217          call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
 218          do 30 nn=1,nfound               ! all neighbours in mdim dimensions
 219             np=nlist(nn)
 220             if(np.lt.(m-1)*id+1) goto 30
 221             ipairs(mdim)=ipairs(mdim)+1
 222             do 40 i=1,m-1
 223              do 41 iim=1,mdim
 224                if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
 225                idim=i*mdim+iim
 226                ipairs(idim)=ipairs(idim)+1  ! neighbours mdim+1..m dimensions
 227  41           continue
 228  40          continue
 229  30         continue
 230  20      if(n-(m-1)*id.ge.ncmin.and.ipairs(m*mdim).ge.ipmin) goto 1
 231       n=n-1
 232  1    s=real(n-(m-1)*id)*real(nmaxy-(m-1)*id) ! normalisation
 233       do 50 i=mdim,mdim*m
 234  50      if(s.gt.0.) c(i)=ipairs(i)/s
 235
 236       return
 237       end
 238 c>---------------------------------------------------------------------
 239
 240
 241
 242
 243
 244