--- /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 cross-recurrence plot
+c authors H. Kantz & T. Schreiber (2004)
+c modified by H. Kantz Feb 2007 (multivariate version)
+ program xrecur
+ parameter(nx=1000000,me=30,meps=1000,mx=10,nx1=10000000)
+ dimension x(nx,mx), y(nx,mx)
+ dimension x1(nx1),y1(nx1)
+ integer mlist(2), inot(nx1), icol(mx)
+ character*72 file1, file2, fout
+ data eps/1e-3/, id0/1/
+ data iverb/1/, xperc/100./
+
+ call whatido(
+ ."cross recurrence plot of two scalar or vector valued data sets"
+ . ,iverb)
+
+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: l1,l2, x1,x2
+c - univariate data: maximum time series length nx1
+c - multivariate data: maximum length nx, maximum dimension mx
+c
+c either: fixed epsilon, plot certain percentage of all
+c neighbours found
+c or: fix the numbers of points of the y time series to be
+c found as neighbours of the x time series: non-symmetric!
+c norm: max-norm after rescaling of the individual components
+c unless rescaling is switched off
+c
+c=================================================================
+
+ 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 (mc.gt.2) print*,'extra arguments of -m ignored'
+ endif
+ ntmin=0
+ kmin=ican("k",0)
+ idown1=ican("s",1)
+ idown2=ican("S",1)
+ eps=fcan("r",eps)
+ xperc=fcan("%",xperc)
+ nmaxx=ican("l",nx)
+ nmaxy=ican("L",nx)
+ nexcl1=ican("x",0)
+ nexcl2=ican("X",0)
+ iscal=0
+ iscal=lopt("n",1)
+
+ 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
+c
+ call nthstring(1,file1)
+ if(file1.eq."-") stop 'missing input filename'
+ if (mdim.gt.1)
+ . call xreadfile(nmaxx,mdim,nx,x,nexcl1,icol,file1,iverb)
+ if (mdim.eq.1)
+ . call xreadfile(nmaxx,1,nx1,x1,nexcl1,icol,file1,iverb)
+c
+ call nthstring(2,file2)
+ if(file2.eq."-") stop 'missing second input filename'
+ if (mdim.gt.1)
+ . call xreadfile(nmaxy,mdim,nx,y,nexcl2,icol,file2,iverb)
+ if (mdim.eq.1)
+ . call xreadfile(nmaxy,1,nx1,y1,nexcl2,icol,file2,iverb)
+
+ if(isout.eq.1) then
+ call addsuff(fout,file1,"_")
+ call addsuff(fout,fout,file2)
+ call addsuff(fout,fout,"_xrec")
+ endif
+
+c rescale data if flag -n is not set (iscal=1)
+ if (iscal.ne.1) then
+ print*,'normalizing data to unit interval'
+
+ if (mdim.gt.1) then
+
+c rescaling each component of file 1
+ do imx=1,mdim
+ xmin=x(1,imx)
+ xmax=x(1,imx)
+ do i1=2,nmaxx
+ xmax=max(xmax,x(i1,imx))
+ xmin=min(xmin,x(i1,imx))
+ enddo
+ scal=.9999d0/(xmax-xmin)
+ do i1=1,nmaxx
+ x(i1,imx)=(x(i1,imx)-xmin)*scal
+ enddo
+ enddo
+
+c rescaling each component of file 2
+
+ do imx=1,mdim
+ xmin=y(1,imx)
+ xmax=y(1,imx)
+ do i1=2,nmaxy
+ xmax=max(xmax,y(i1,imx))
+ xmin=min(xmin,y(i1,imx))
+ enddo
+ scal=.9999d0/(xmax-xmin)
+ do i1=1,nmaxy
+ y(i1,imx)=(y(i1,imx)-xmin)*scal
+ enddo
+ enddo
+
+ else
+
+c rescaling the single component of file 1
+ xmin=x1(1)
+ xmax=x1(1)
+ do i1=2,nmaxx
+ xmax=max(xmax,x1(i1))
+ xmin=min(xmin,x1(i1))
+ enddo
+ scal=.9999d0/(xmax-xmin)
+ do i1=1,nmaxx
+ x1(i1)=(x1(i1)-xmin)*scal
+ enddo
+
+c rescaling the single component of file 2
+
+ xmin=y1(1)
+ xmax=y1(1)
+ do i1=2,nmaxy
+ xmax=max(xmax,y1(i1))
+ xmin=min(xmin,y1(i1))
+ enddo
+ scal=.9999d0/(xmax-xmin)
+ do i1=1,nmaxy
+ y1(i1)=(y1(i1)-xmin)*scal
+ enddo
+
+ endif
+ endif
+
+ call outfile(fout,iunit,1)
+ ntot=0
+
+ if (kmin.eq.0) then
+c search all neighbours with distance < eps
+
+ xperc=xperc/100.
+
+ if (mdim.eq.1) then
+ call crossrec(nmaxx,x1,nmaxy,y1,eps,
+ . id,mmax,iunit,xperc,ntot,idown1,idown2)
+ else
+ call mcrossrec(nmaxx,x,nmaxy,y,eps,
+ . id,mmax,mdim,iunit,xperc,ntot,idown1,idown2)
+ endif
+
+c>-----------------------------------------------------------
+ else
+
+ do i=(mmax-1)*id+1,nmaxx
+ inot(i)=1
+ enddo
+ epsfac=1.1
+ eps=eps/epsfac
+
+ do 10 io=1,100
+ eps=eps*epsfac
+ if (mdim.eq.1) then
+ call crossrec1(nmaxx,x1,nmaxy,y1,eps,
+ . id,mmax,kmin,iunit,inot,ntodo,ntot,idown1,idown2)
+ else
+ call mcrossrec1(nmaxx,x,nmaxy,y,eps,id,mmax,mdim,kmin,iunit,
+ . inot,ntodo,ntot,idown1,idown2)
+ endif
+ if (iverb.eq.1) print*,eps,ntodo,ntot
+ if (ntodo.eq.0) goto 7
+ 10 continue
+c>--------------------------------------------------------
+ endif
+
+ 7 close(iunit)
+ print*,ntot,' points contained in the recurrence plot'
+ if (kmin.gt.0) print*,'last epsilon:',eps
+ if (kmin.gt.0) print*,'average number of neighbours:',
+ . real(ntot)*idown1/nmaxx
+ if (kmin.eq.0) print*,'using eps=',eps
+ end
+c========================================================
+ subroutine usage()
+c usage message
+
+ call whatineed(
+ ."[ -m#,# -d# -r# -k# -o outfile -l# -x# -L# -X# -c#[,#] -%# -V#
+ . -n -h] file1 file2")
+ call popt("m","# of components, embedding dimension [1,2]")
+ call popt("c","columns to be read [1,2,3,...,# of components]")
+ call popt("d","delay [1]")
+ call popt("r",
+ ."diameter of the neighbourhood as absolute value [.001]")
+ call popt("k",
+ ."find the # closest points, starting with diameter r")
+ call popt("%",
+ ."print only percentage of dots [100], no effect if -k is set")
+ call popt("l","length of time series 1 to be read [all data]")
+ call popt("x","# of initial lines in 1 to be skipped [0]")
+ call popt("s","use only every # delay vector of file 1 [1]")
+ call popt("L","length of time series 2 to be read [all data]")
+ call popt("X","# of initial lines in 2 to be skipped [0]")
+ call popt("S","use only every # delay vector of file 2 [1]")
+ call popt("n","if set: do NOT normalize data to unit interval")
+ call pout("file1_file2_xrec")
+ call pall()
+ stop
+ end
+c>--------------------------------------------------------------------
+ subroutine crossrec(nmaxx,x,nmaxy,y,eps,
+ . id,m,iunit,xperc,ntot,idown1,idown2)
+ parameter(im=100,ii=100000000,nx=1000000)
+ dimension y(nx),x(nx)
+ dimension jh(0:im*im),jpntr(nx),nlist(nx)
+ nseed=13413241
+
+ if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec: make nx larger."
+ mb=min(m,2)
+
+ call base(nmaxy,y,id,mb,jh,jpntr,eps)
+ nnull=(m-1)*id+1
+ do 20 n=nnull,nmaxx,idown1
+ call neigh(nx,x,y,n,nx,id,mb,jh,jpntr,eps,nlist,nfound)
+ do 30 nn=1,nfound ! all neighbours in two dimensions
+ np=nlist(nn)
+ if (np.lt.nnull) goto 30
+ if (mod(np-nnull,idown2).ne.0) goto 30
+ do 40 i=mb,m-1
+ if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
+ 40 continue
+ call random(nseed,rr)
+ if (rr.le.xperc) write(iunit,*)n,np
+ if (rr.le.xperc) ntot=ntot+1
+ 30 continue
+ 20 continue
+ end
+c>--------------------------------------------------------------------
+ subroutine mcrossrec(nmaxx,x,nmaxy,y,eps,
+ . id,m,mdim,iunit,xperc,ntot,idown1,idown2)
+ parameter(im=100,nx=1000000,mx=10)
+ dimension y(nx,mx),x(nx,mx)
+ dimension jh(0:im*im),jpntr(nx),nlist(nx)
+ dimension vx(mx)
+ nseed=134512331
+
+ if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec: make nx larger."
+
+ call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
+ nnull=(m-1)*id+1
+ do 20 n=nnull,nmaxx,idown1
+ 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.nnull) goto 30
+ if (mod(np-nnull,idown2).ne.0) goto 30
+ 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
+ 41 continue
+ 40 continue
+ call random(nseed,rr)
+ if (rr.le.xperc) write(iunit,*)n,np
+ if (rr.le.xperc) ntot=ntot+1
+ 30 continue
+ 20 continue
+ return
+ end
+
+ subroutine crossrec1(nmaxx,x,nmaxy,y,eps,id,m,kmin,iunit,
+ . inot,ntodo,ntot,idown1,idown2)
+ parameter(im=100,ii=100000000,nx=1000000)
+ dimension y(nmaxy),x(nmaxx),inot(nmaxx)
+ dimension jh(0:im*im),jpntr(nx),nlist(nx)
+ ntodo=0
+
+ if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec1: make nx larger."
+ mb=min(m,2)
+
+ call base(nmaxy,y,id,mb,jh,jpntr,eps)
+ nnull=(m-1)*id+1
+ do 20 n=nnull,nmaxx,idown1
+ if (inot(n).eq.0) goto 20
+ ntodo=ntodo+1
+ call neigh(nx,x,y,n,nx,id,mb,jh,jpntr,eps,nlist,nfound)
+ if (nfound.lt.kmin) goto 20
+ nreal=0
+ do 30 nn=1,nfound ! all neighbours in two dimensions
+ np=nlist(nn)
+ if(np.lt.nnull) goto 30
+ if (mod(np-nnull,idown2).ne.0) goto 30
+ do 40 i=mb,m-1
+ if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
+ 40 continue
+ nreal=nreal+1
+ nlist(nreal)=np
+ 30 continue
+ if (nreal.lt.kmin) goto 20
+ ntodo=ntodo-1
+ inot(n)=0
+ ntot=ntot+nreal
+ do in=1,nreal
+ write(iunit,*)n,nlist(in)
+ enddo
+ 20 continue
+ end
+c>--------------------------------------------------------------------
+ subroutine mcrossrec1(nmaxx,x,nmaxy,y,eps,
+ . id,m,mdim,kmin,iunit,inot,
+ . ntodo,ntot,idown1,idown2)
+ parameter(im=100,nx=1000000,mx=10)
+ dimension y(nx,mx),x(nx,mx),inot(nmaxx)
+ dimension jh(0:im*im),jpntr(nx),nlist(nx)
+ dimension vx(mx)
+ ntodo=0
+
+ if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec1: make nx larger."
+
+ call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
+ nnull=(m-1)*id+1
+ do 20 n=nnull,nmaxx,idown1
+ if (inot(n).eq.0) goto 20
+ ntodo=ntodo+1
+ do ii=1,mdim
+ vx(ii)=x(n,ii)
+ enddo
+ call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
+ if (nfound.lt.kmin) goto 20
+ nreal=0
+ do 30 nn=1,nfound ! all neighbours in mdim dimensions
+ np=nlist(nn)
+ if(np.lt.nnull) goto 30
+ if (mod(np-nnull,idown2).ne.0) goto 30
+ do 40 i=0,m-1
+ do 41 iim=1,mdim
+ if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
+ 41 continue
+ 40 continue
+ nreal=nreal+1
+ nlist(nreal)=np
+ 30 continue
+ if (nreal.lt.kmin) goto 20
+ inot(n)=0
+ ntot=ntot+nreal
+ do in=1,nreal
+ write(iunit,*)n,nlist(in)
+ enddo
+ ntodo=ntodo-1
+ 20 continue
+ return
+ end
+
+ subroutine random(iseed,s)
+c
+c random number generator of Park & Miller
+ integer*8 ifac,ibase,iargument
+ ifac=7**5
+ ibase=2**30-1
+ im=im+2**30
+ iargument=iseed
+ iargument=mod(iargument*ifac,ibase)
+ s=float(iargument)/float(ibase)
+ iseed=iargument
+ return
+ end
+c>------------------------------------
git://source.jalview.org / jabaws.git / blobdiff