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===========================================================================
22 c cross-recurrence plot
23 c authors H. Kantz & T. Schreiber (2004)
24 c modified by H. Kantz Feb 2007 (multivariate version)
26 parameter(nx=1000000,me=30,meps=1000,mx=10,nx1=10000000)
27 dimension x(nx,mx), y(nx,mx)
28 dimension x1(nx1),y1(nx1)
29 integer mlist(2), inot(nx1), icol(mx)
30 character*72 file1, file2, fout
31 data eps/1e-3/, id0/1/
32 data iverb/1/, xperc/100./
35 ."cross recurrence plot of two scalar or vector valued data sets"
38 c ================================================================
39 c assume two input files (file1, file2)
40 c - with identical structure in terms of colums
41 c - with identical embedding spaces
42 c - with individual length and exclusions: l1,l2, x1,x2
43 c - univariate data: maximum time series length nx1
44 c - multivariate data: maximum length nx, maximum dimension mx
46 c either: fixed epsilon, plot certain percentage of all
48 c or: fix the numbers of points of the y time series to be
49 c found as neighbours of the x time series: non-symmetric!
50 c norm: max-norm after rescaling of the individual components
51 c unless rescaling is switched off
53 c=================================================================
60 call imcan("m",2,mc,mlist)
64 if (mc.gt.2) print*,'extra arguments of -m ignored'
79 call columns(mc,mx,icol)
80 if (mc.gt.0.and.mc.ne.mdim) stop 'improper number of columns'
82 if(fout.eq." ") isout=1
84 call nthstring(1,file1)
85 if(file1.eq."-") stop 'missing input filename'
87 . call xreadfile(nmaxx,mdim,nx,x,nexcl1,icol,file1,iverb)
89 . call xreadfile(nmaxx,1,nx1,x1,nexcl1,icol,file1,iverb)
91 call nthstring(2,file2)
92 if(file2.eq."-") stop 'missing second input filename'
94 . call xreadfile(nmaxy,mdim,nx,y,nexcl2,icol,file2,iverb)
96 . call xreadfile(nmaxy,1,nx1,y1,nexcl2,icol,file2,iverb)
99 call addsuff(fout,file1,"_")
100 call addsuff(fout,fout,file2)
101 call addsuff(fout,fout,"_xrec")
104 c rescale data if flag -n is not set (iscal=1)
106 print*,'normalizing data to unit interval'
110 c rescaling each component of file 1
115 xmax=max(xmax,x(i1,imx))
116 xmin=min(xmin,x(i1,imx))
118 scal=.9999d0/(xmax-xmin)
120 x(i1,imx)=(x(i1,imx)-xmin)*scal
124 c rescaling each component of file 2
130 xmax=max(xmax,y(i1,imx))
131 xmin=min(xmin,y(i1,imx))
133 scal=.9999d0/(xmax-xmin)
135 y(i1,imx)=(y(i1,imx)-xmin)*scal
141 c rescaling the single component of file 1
145 xmax=max(xmax,x1(i1))
146 xmin=min(xmin,x1(i1))
148 scal=.9999d0/(xmax-xmin)
150 x1(i1)=(x1(i1)-xmin)*scal
153 c rescaling the single component of file 2
158 xmax=max(xmax,y1(i1))
159 xmin=min(xmin,y1(i1))
161 scal=.9999d0/(xmax-xmin)
163 y1(i1)=(y1(i1)-xmin)*scal
169 call outfile(fout,iunit,1)
173 c search all neighbours with distance < eps
178 call crossrec(nmaxx,x1,nmaxy,y1,eps,
179 . id,mmax,iunit,xperc,ntot,idown1,idown2)
181 call mcrossrec(nmaxx,x,nmaxy,y,eps,
182 . id,mmax,mdim,iunit,xperc,ntot,idown1,idown2)
185 c>-----------------------------------------------------------
188 do i=(mmax-1)*id+1,nmaxx
197 call crossrec1(nmaxx,x1,nmaxy,y1,eps,
198 . id,mmax,kmin,iunit,inot,ntodo,ntot,idown1,idown2)
200 call mcrossrec1(nmaxx,x,nmaxy,y,eps,id,mmax,mdim,kmin,iunit,
201 . inot,ntodo,ntot,idown1,idown2)
203 if (iverb.eq.1) print*,eps,ntodo,ntot
204 if (ntodo.eq.0) goto 7
206 c>--------------------------------------------------------
210 print*,ntot,' points contained in the recurrence plot'
211 if (kmin.gt.0) print*,'last epsilon:',eps
212 if (kmin.gt.0) print*,'average number of neighbours:',
213 . real(ntot)*idown1/nmaxx
214 if (kmin.eq.0) print*,'using eps=',eps
216 c========================================================
221 ."[ -m#,# -d# -r# -k# -o outfile -l# -x# -L# -X# -c#[,#] -%# -V#
222 . -n -h] file1 file2")
223 call popt("m","# of components, embedding dimension [1,2]")
224 call popt("c","columns to be read [1,2,3,...,# of components]")
225 call popt("d","delay [1]")
227 ."diameter of the neighbourhood as absolute value [.001]")
229 ."find the # closest points, starting with diameter r")
231 ."print only percentage of dots [100], no effect if -k is set")
232 call popt("l","length of time series 1 to be read [all data]")
233 call popt("x","# of initial lines in 1 to be skipped [0]")
234 call popt("s","use only every # delay vector of file 1 [1]")
235 call popt("L","length of time series 2 to be read [all data]")
236 call popt("X","# of initial lines in 2 to be skipped [0]")
237 call popt("S","use only every # delay vector of file 2 [1]")
238 call popt("n","if set: do NOT normalize data to unit interval")
239 call pout("file1_file2_xrec")
243 c>--------------------------------------------------------------------
244 subroutine crossrec(nmaxx,x,nmaxy,y,eps,
245 . id,m,iunit,xperc,ntot,idown1,idown2)
246 parameter(im=100,ii=100000000,nx=1000000)
247 dimension y(nx),x(nx)
248 dimension jh(0:im*im),jpntr(nx),nlist(nx)
251 if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec: make nx larger."
254 call base(nmaxy,y,id,mb,jh,jpntr,eps)
256 do 20 n=nnull,nmaxx,idown1
257 call neigh(nx,x,y,n,nx,id,mb,jh,jpntr,eps,nlist,nfound)
258 do 30 nn=1,nfound ! all neighbours in two dimensions
260 if (np.lt.nnull) goto 30
261 if (mod(np-nnull,idown2).ne.0) goto 30
263 if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
265 call random(nseed,rr)
266 if (rr.le.xperc) write(iunit,*)n,np
267 if (rr.le.xperc) ntot=ntot+1
271 c>--------------------------------------------------------------------
272 subroutine mcrossrec(nmaxx,x,nmaxy,y,eps,
273 . id,m,mdim,iunit,xperc,ntot,idown1,idown2)
274 parameter(im=100,nx=1000000,mx=10)
275 dimension y(nx,mx),x(nx,mx)
276 dimension jh(0:im*im),jpntr(nx),nlist(nx)
280 if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec: make nx larger."
282 call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
284 do 20 n=nnull,nmaxx,idown1
288 call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
289 do 30 nn=1,nfound ! all neighbours in mdim dimensions
291 if (np.lt.nnull) goto 30
292 if (mod(np-nnull,idown2).ne.0) goto 30
295 if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
298 call random(nseed,rr)
299 if (rr.le.xperc) write(iunit,*)n,np
300 if (rr.le.xperc) ntot=ntot+1
306 subroutine crossrec1(nmaxx,x,nmaxy,y,eps,id,m,kmin,iunit,
307 . inot,ntodo,ntot,idown1,idown2)
308 parameter(im=100,ii=100000000,nx=1000000)
309 dimension y(nmaxy),x(nmaxx),inot(nmaxx)
310 dimension jh(0:im*im),jpntr(nx),nlist(nx)
313 if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec1: make nx larger."
316 call base(nmaxy,y,id,mb,jh,jpntr,eps)
318 do 20 n=nnull,nmaxx,idown1
319 if (inot(n).eq.0) goto 20
321 call neigh(nx,x,y,n,nx,id,mb,jh,jpntr,eps,nlist,nfound)
322 if (nfound.lt.kmin) goto 20
324 do 30 nn=1,nfound ! all neighbours in two dimensions
326 if(np.lt.nnull) goto 30
327 if (mod(np-nnull,idown2).ne.0) goto 30
329 if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
334 if (nreal.lt.kmin) goto 20
339 write(iunit,*)n,nlist(in)
343 c>--------------------------------------------------------------------
344 subroutine mcrossrec1(nmaxx,x,nmaxy,y,eps,
345 . id,m,mdim,kmin,iunit,inot,
346 . ntodo,ntot,idown1,idown2)
347 parameter(im=100,nx=1000000,mx=10)
348 dimension y(nx,mx),x(nx,mx),inot(nmaxx)
349 dimension jh(0:im*im),jpntr(nx),nlist(nx)
353 if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec1: make nx larger."
355 call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
357 do 20 n=nnull,nmaxx,idown1
358 if (inot(n).eq.0) goto 20
363 call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
364 if (nfound.lt.kmin) goto 20
366 do 30 nn=1,nfound ! all neighbours in mdim dimensions
368 if(np.lt.nnull) goto 30
369 if (mod(np-nnull,idown2).ne.0) goto 30
372 if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
378 if (nreal.lt.kmin) goto 20
382 write(iunit,*)n,nlist(in)
389 subroutine random(iseed,s)
391 c random number generator of Park & Miller
392 integer*8 ifac,ibase,iargument
397 iargument=mod(iargument*ifac,ibase)
398 s=float(iargument)/float(ibase)
402 c>------------------------------------