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=========================================================================
23 c cross-correlation integral xc2
24 c see H. Kantz, Phys.Rev.E49, 5091 (1994)
26 c authors: T. Schreiber & H. Kantz (1998)
27 c multivariate version: H. Kantz (Jan 2007)
29 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)
35 character*72 file1, file2, fout
36 data ipmin/1000/, res/2./, eps0/1e-30/, epsm/1e30/, id0/1/
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
48 c no rescaling, since xc2 does not make sense if datasets are
49 c not used in their original scalings
50 c=================================================================
52 call whatido("cross correlation sum of two data sets",iverb)
57 call imcan("M",2,mc,mlist)
61 if (mmax*mdim.lt.2) stop 'Increase embedding dimension'
62 if (mc.gt.2) print*, 'extra arguments of -m ignored'
77 call columns(mc,mx,icol)
79 if (mc.gt.0.and.mc.ne.mdim) stop 'improper number of columns'
81 if(fout.eq." ") isout=1
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)
91 call addsuff(fout,file1,"_")
92 call addsuff(fout,fout,file2)
93 call addsuff(fout,fout,"_xc2")
98 call minmax(nmaxx,x(1,imx),xmin,xmax)
99 epsmax=1.001*max(xmax-xmin,epsmax)
102 call minmax(nmaxy,y(1,imx),xmin,xmax)
103 epsmax=1.001*max(xmax-xmin,epsmax)
108 do 10 epsl=log(min(epsm,epsmax)),log(eps0),-log(feps)
110 if(neps.gt.meps) stop "xc2: make meps larger"
115 call crosscor(nmaxx,x,nmaxy,y,eps(neps)
116 . ,id,mmax,c(1,neps),ncmin,ipmin)
118 call mcrosscor(nmaxx,x,nmaxy,y,eps(neps),
119 . id,mmax,mdim,c(1,neps),ncmin,ipmin)
124 30 if(c(m,neps).eq.0.) goto 1
127 if(mdd.eq.mdim-1) stop
129 call outfile(fout,iunit,iverb)
130 do 40 m=mdd1,mdeps(1)
131 write(iunit,'(4h#m= ,i5)') m
133 if(mdeps(nn).lt.m) goto 2
134 50 write(iunit,*) eps(nn), c(m,nn)
138 10 write(istderr(),*) eps(neps), mdd, c(mdd,neps)
141 c>--------------------------------------------------------------------
146 . "-M#,# [-d# -n# -N# -## -r# -R#"//
147 . " -o outfile -l# -x# -L# -X# -c#[,#] -V# -h] file1 file2")
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]")
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]")
162 ."columns to be read [1,2,3,.., # of components]")
163 call pout("file1_file2_xc2")
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)
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."
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
184 if(np.lt.(m-1)*id+1) goto 30
185 ipairs(1)=ipairs(1)+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
190 20 if(n-(m-1)*id.ge.ncmin.and.ipairs(m-1).ge.ipmin) goto 1
192 1 s=real(n-(m-1)*id)*real(nmaxx-(m-1)*id) ! normalisation
194 50 if(s.gt.0.) c(i+1)=ipairs(i)/s
196 c>--------------------------------------------------------------------
197 subroutine mcrosscor(nmaxx,x,nmaxy,y,eps,id,m,mdim,c,ncmin,ipmin)
199 parameter(im=100,nx=1000000,mm=30,mx=10)
201 dimension y(nx,mx),x(nx,mx)
202 dimension jh(0:im*im),ipairs(mm),c(m*mdim),jpntr(nx),nlist(nx)
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."
208 if (m*mdim.gt.mm) stop 'embedding x spatial dimension < 30 !'
212 call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
213 do 20 n=(m-1)*id+1,nmaxx
217 call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
218 do 30 nn=1,nfound ! all neighbours in mdim dimensions
220 if(np.lt.(m-1)*id+1) goto 30
221 ipairs(mdim)=ipairs(mdim)+1
224 if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
226 ipairs(idim)=ipairs(idim)+1 ! neighbours mdim+1..m dimensions
230 20 if(n-(m-1)*id.ge.ncmin.and.ipairs(m*mdim).ge.ipmin) goto 1
232 1 s=real(n-(m-1)*id)*real(nmaxy-(m-1)*id) ! normalisation
234 50 if(s.gt.0.) c(i)=ipairs(i)/s
238 c>---------------------------------------------------------------------