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 c cross-correlation integral xc2 c see H. Kantz, Phys.Rev.E49, 5091 (1994) c c authors: T. Schreiber & H. Kantz (1998) c multivariate version: H. Kantz (Jan 2007) c c========================================================================= c parameter(nx=1000000,me=30,meps=1000,mx=10) dimension x(nx,mx), c(me,meps), eps(meps), mdeps(meps), icol(mx) dimension y(nx,mx) integer mlist(2) character*72 file1, file2, fout data ipmin/1000/, res/2./, eps0/1e-30/, epsm/1e30/, id0/1/ data iverb/1/ c======================================================================= c assume two input files (file1, file2) c - with identical structure in terms of colums c - with identical embedding spaces c - with individual length and exclusions: l, L, x, X c - multivariate data: maximum length nx, maximum dimension mx c c norm: max-norm c c no rescaling, since xc2 does not make sense if datasets are c not used in their original scalings c================================================================= call whatido("cross correlation sum of two data sets",iverb) id=ican("d",id0) mmax=2 mdim=1 call imcan("M",2,mc,mlist) if (mc.ge.2) then mmax=mlist(2) mdim=mlist(1) if (mmax*mdim.lt.2) stop 'Increase embedding dimension' if (mc.gt.2) print*, 'extra arguments of -m ignored' endif ntmin=0 ncmin=ican("n",1000) ipmin=ican("N",ipmin) res=fcan("#",res) feps=2**(1./res) eps0=fcan("r",eps0) epsm=fcan("R",epsm) nmaxx=ican("l",nx) nmaxy=ican("L",nx) nexcl1=ican("x",0) nexcl2=ican("X",0) call columns(mc,mx,icol) if (mc.gt.0.and.mc.ne.mdim) stop 'improper number of columns' isout=igetout(fout,0) if(fout.eq." ") isout=1 call nthstring(1,file1) if(file1.eq."-") stop "first input file name missing" call xreadfile(nmaxx,mdim,nx,x,nexcl1,icol,file1,iverb) call nthstring(2,file2) if(file2.eq."-") stop "second input file name missing" call xreadfile(nmaxy,mdim,nx,y,nexcl2,icol,file2,iverb) if(isout.eq.1) then call addsuff(fout,file1,"_") call addsuff(fout,fout,file2) call addsuff(fout,fout,"_xc2") endif epsmax=0. do imx=1,mmax call minmax(nmaxx,x(1,imx),xmin,xmax) epsmax=1.001*max(xmax-xmin,epsmax) enddo do imx=1,mmax call minmax(nmaxy,y(1,imx),xmin,xmax) epsmax=1.001*max(xmax-xmin,epsmax) enddo neps=0 do 10 epsl=log(min(epsm,epsmax)),log(eps0),-log(feps) neps=neps+1 if(neps.gt.meps) stop "xc2: make meps larger" eps(neps)=exp(epsl) do 20 m=1,mmax*mdim 20 c(m,neps)=0 if (mdim.eq.1) then call crosscor(nmaxx,x,nmaxy,y,eps(neps) . ,id,mmax,c(1,neps),ncmin,ipmin) else call mcrosscor(nmaxx,x,nmaxy,y,eps(neps), . id,mmax,mdim,c(1,neps),ncmin,ipmin) endif mdd=mmax*mdim mdd1=max(2,mdim) do 30 m=mdd1,mdd 30 if(c(m,neps).eq.0.) goto 1 m=mdd+1 1 mdd=m-1 if(mdd.eq.mdim-1) stop mdeps(neps)=mdd call outfile(fout,iunit,iverb) do 40 m=mdd1,mdeps(1) write(iunit,'(4h#m= ,i5)') m do 50 nn=1,neps if(mdeps(nn).lt.m) goto 2 50 write(iunit,*) eps(nn), c(m,nn) 2 write(iunit,'()') 40 write(iunit,'()') close(iunit) 10 write(istderr(),*) eps(neps), mdd, c(mdd,neps) stop end c>-------------------------------------------------------------------- subroutine usage() c usage message call whatineed( . "-M#,# [-d# -n# -N# -## -r# -R#"// . " -o outfile -l# -x# -L# -X# -c#[,#] -V# -h] file1 file2") call popt("M", ."# of components, maximal embedding dimension [1,2]") call popt("d","delay [1]") call popt("n","minimal number of center points [1000]") call popt("N","maximal number of pairs [1000]") call popt("#","resolution, values per octave [2]") call popt("r", . "minimal scale to be probed (as long as pairs found)") call popt("R","maximal scale to be probed [xmax-xmin]") call popt("l","length of time series 1 to be read [all data]") call popt("x","# of initial lines of 1 to be skipped [0]") call popt("L","length of time series 2 to be read [all data]") call popt("X","# of initial lines of 2 to be skipped [0]") call popt("c", ."columns to be read [1,2,3,.., # of components]") call pout("file1_file2_xc2") call pall() stop end c>-------------------------------------------------------------------- subroutine crosscor(nmaxx,x,nmaxy,y,eps,id,m,c,ncmin,ipmin) parameter(im=100,ii=100000000,nx=1000000,mm=30) dimension y(nmaxy),x(nmaxx) dimension jh(0:im*im),ipairs(mm),c(m),jpntr(nx),nlist(nx) if(nmaxx.gt.nx.or.m.gt.mm) stop "crosscor: make mm/nx larger." if(nmaxy.gt.nx.or.m.gt.mm) stop "crosscor: make mm/nx larger." do 10 i=1,m-1 10 ipairs(i)=0 mb=min(m,2) call base(nmaxx,x,id,mb,jh,jpntr,eps) do 20 n=(m-1)*id+1,nmaxy call neigh(nx,y,x,n,nmaxx,id,mb,jh,jpntr,eps,nlist,nfound) do 30 nn=1,nfound ! all neighbours in two dimensions np=nlist(nn) if(np.lt.(m-1)*id+1) goto 30 ipairs(1)=ipairs(1)+1 do 40 i=mb,m-1 if(abs(y(n-i*id)-x(np-i*id)).ge.eps) goto 30 40 ipairs(i)=ipairs(i)+1 ! neighbours in 3..m dimensions 30 continue 20 if(n-(m-1)*id.ge.ncmin.and.ipairs(m-1).ge.ipmin) goto 1 n=n-1 1 s=real(n-(m-1)*id)*real(nmaxx-(m-1)*id) ! normalisation do 50 i=1,m-1 50 if(s.gt.0.) c(i+1)=ipairs(i)/s end c>-------------------------------------------------------------------- subroutine mcrosscor(nmaxx,x,nmaxy,y,eps,id,m,mdim,c,ncmin,ipmin) parameter(im=100,nx=1000000,mm=30,mx=10) dimension y(nx,mx),x(nx,mx) dimension jh(0:im*im),ipairs(mm),c(m*mdim),jpntr(nx),nlist(nx) dimension vx(mm) if(nmaxx.gt.nx.or.m.gt.mm) stop "mcrosscor: make mm/nx larger." if(nmaxy.gt.nx.or.m.gt.mm) stop "mcrosscor: make mm/nx larger." if (m*mdim.gt.mm) stop 'embedding x spatial dimension < 30 !' do 10 i=1,m*mdim 10 ipairs(i)=0 call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps) do 20 n=(m-1)*id+1,nmaxx do ii=1,mdim vx(ii)=x(n,ii) enddo call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound) do 30 nn=1,nfound ! all neighbours in mdim dimensions np=nlist(nn) if(np.lt.(m-1)*id+1) goto 30 ipairs(mdim)=ipairs(mdim)+1 do 40 i=1,m-1 do 41 iim=1,mdim if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30 idim=i*mdim+iim ipairs(idim)=ipairs(idim)+1 ! neighbours mdim+1..m dimensions 41 continue 40 continue 30 continue 20 if(n-(m-1)*id.ge.ncmin.and.ipairs(m*mdim).ge.ipmin) goto 1 n=n-1 1 s=real(n-(m-1)*id)*real(nmaxy-(m-1)*id) ! normalisation do 50 i=mdim,mdim*m 50 if(s.gt.0.) c(i)=ipairs(i)/s return end c>---------------------------------------------------------------------