website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/xrecur.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 cross-recurrence plot
  23 c authors H. Kantz & T. Schreiber (2004)
  24 c modified by H. Kantz Feb 2007 (multivariate version)
  25       program xrecur
  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./
  33
  34       call whatido(
  35      ."cross recurrence plot of two scalar or vector valued data sets"
  36      .             ,iverb)
  37
  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
  45 c
  46 c either: fixed epsilon, plot certain percentage of all
  47 c         neighbours found
  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
  52 c
  53 c=================================================================
  54
  55       id=ican("d",id0)
  56
  57       mmax=2
  58       mdim=1
  59
  60       call imcan("m",2,mc,mlist)
  61       if (mc.ge.2) then
  62        mmax=mlist(2)
  63        mdim=mlist(1)
  64        if (mc.gt.2) print*,'extra arguments of -m ignored'
  65       endif
  66       ntmin=0
  67       kmin=ican("k",0)
  68       idown1=ican("s",1)
  69       idown2=ican("S",1)
  70       eps=fcan("r",eps)
  71       xperc=fcan("%",xperc)
  72       nmaxx=ican("l",nx)
  73       nmaxy=ican("L",nx)
  74       nexcl1=ican("x",0)
  75       nexcl2=ican("X",0)
  76       iscal=0
  77       iscal=lopt("n",1)
  78
  79       call columns(mc,mx,icol)
  80       if (mc.gt.0.and.mc.ne.mdim) stop 'improper number of columns'
  81       isout=igetout(fout,0)
  82       if(fout.eq." ") isout=1
  83 c
  84       call nthstring(1,file1)
  85       if(file1.eq."-") stop 'missing input filename'
  86       if (mdim.gt.1)
  87      .  call xreadfile(nmaxx,mdim,nx,x,nexcl1,icol,file1,iverb)
  88       if (mdim.eq.1)
  89      .  call xreadfile(nmaxx,1,nx1,x1,nexcl1,icol,file1,iverb)
  90 c
  91       call nthstring(2,file2)
  92       if(file2.eq."-") stop 'missing second input filename'
  93       if (mdim.gt.1)
  94      . call xreadfile(nmaxy,mdim,nx,y,nexcl2,icol,file2,iverb)
  95       if (mdim.eq.1)
  96      . call xreadfile(nmaxy,1,nx1,y1,nexcl2,icol,file2,iverb)
  97
  98       if(isout.eq.1) then
  99         call addsuff(fout,file1,"_")
 100         call addsuff(fout,fout,file2)
 101         call addsuff(fout,fout,"_xrec")
 102       endif
 103
 104 c     rescale data if flag -n is not set (iscal=1)
 105       if (iscal.ne.1) then
 106        print*,'normalizing data to unit interval'
 107
 108        if (mdim.gt.1) then
 109
 110 c      rescaling each component of file 1
 111         do imx=1,mdim
 112          xmin=x(1,imx)
 113          xmax=x(1,imx)
 114          do i1=2,nmaxx
 115           xmax=max(xmax,x(i1,imx))
 116           xmin=min(xmin,x(i1,imx))
 117          enddo
 118          scal=.9999d0/(xmax-xmin)
 119          do i1=1,nmaxx
 120           x(i1,imx)=(x(i1,imx)-xmin)*scal
 121          enddo
 122         enddo
 123
 124 c     rescaling each component of file 2
 125
 126         do imx=1,mdim
 127          xmin=y(1,imx)
 128          xmax=y(1,imx)
 129          do i1=2,nmaxy
 130           xmax=max(xmax,y(i1,imx))
 131           xmin=min(xmin,y(i1,imx))
 132          enddo
 133          scal=.9999d0/(xmax-xmin)
 134          do i1=1,nmaxy
 135           y(i1,imx)=(y(i1,imx)-xmin)*scal
 136          enddo
 137         enddo
 138
 139        else
 140
 141 c     rescaling the single component of file 1
 142          xmin=x1(1)
 143          xmax=x1(1)
 144          do i1=2,nmaxx
 145           xmax=max(xmax,x1(i1))
 146           xmin=min(xmin,x1(i1))
 147          enddo
 148          scal=.9999d0/(xmax-xmin)
 149          do i1=1,nmaxx
 150           x1(i1)=(x1(i1)-xmin)*scal
 151          enddo
 152
 153 c     rescaling the single component of file 2
 154
 155          xmin=y1(1)
 156          xmax=y1(1)
 157          do i1=2,nmaxy
 158           xmax=max(xmax,y1(i1))
 159           xmin=min(xmin,y1(i1))
 160          enddo
 161          scal=.9999d0/(xmax-xmin)
 162          do i1=1,nmaxy
 163           y1(i1)=(y1(i1)-xmin)*scal
 164          enddo
 165
 166        endif
 167       endif
 168
 169       call outfile(fout,iunit,1)
 170       ntot=0
 171
 172       if (kmin.eq.0) then
 173 c     search all neighbours with distance < eps
 174
 175         xperc=xperc/100.
 176
 177         if (mdim.eq.1) then
 178           call crossrec(nmaxx,x1,nmaxy,y1,eps,
 179      .                   id,mmax,iunit,xperc,ntot,idown1,idown2)
 180          else
 181           call mcrossrec(nmaxx,x,nmaxy,y,eps,
 182      .              id,mmax,mdim,iunit,xperc,ntot,idown1,idown2)
 183         endif
 184
 185 c>-----------------------------------------------------------
 186       else
 187
 188        do i=(mmax-1)*id+1,nmaxx
 189         inot(i)=1
 190        enddo
 191        epsfac=1.1
 192        eps=eps/epsfac
 193
 194        do 10 io=1,100
 195          eps=eps*epsfac
 196          if (mdim.eq.1) then
 197           call crossrec1(nmaxx,x1,nmaxy,y1,eps,
 198      .               id,mmax,kmin,iunit,inot,ntodo,ntot,idown1,idown2)
 199          else
 200           call mcrossrec1(nmaxx,x,nmaxy,y,eps,id,mmax,mdim,kmin,iunit,
 201      .                    inot,ntodo,ntot,idown1,idown2)
 202          endif
 203          if (iverb.eq.1) print*,eps,ntodo,ntot
 204          if (ntodo.eq.0) goto 7
 205  10    continue
 206 c>--------------------------------------------------------
 207       endif
 208
 209  7    close(iunit)
 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
 215       end
 216 c========================================================
 217       subroutine usage()
 218 c usage message
 219
 220       call whatineed(
 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]")
 226       call popt("r",
 227      ."diameter of the neighbourhood as absolute value [.001]")
 228       call popt("k",
 229      ."find the # closest points, starting with diameter r")
 230       call popt("%",
 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")
 240       call pall()
 241       stop
 242       end
 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)
 249       nseed=13413241
 250
 251       if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec: make nx larger."
 252       mb=min(m,2)
 253
 254       call base(nmaxy,y,id,mb,jh,jpntr,eps)
 255       nnull=(m-1)*id+1
 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
 259             np=nlist(nn)
 260             if (np.lt.nnull) goto 30
 261             if (mod(np-nnull,idown2).ne.0) goto 30
 262             do 40 i=mb,m-1
 263                if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
 264  40         continue
 265             call random(nseed,rr)
 266             if (rr.le.xperc) write(iunit,*)n,np
 267             if (rr.le.xperc) ntot=ntot+1
 268  30         continue
 269  20      continue
 270       end
 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)
 277       dimension vx(mx)
 278       nseed=134512331
 279
 280       if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec: make nx larger."
 281
 282       call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
 283       nnull=(m-1)*id+1
 284       do 20 n=nnull,nmaxx,idown1
 285          do ii=1,mdim
 286           vx(ii)=x(n,ii)
 287          enddo
 288          call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
 289          do 30 nn=1,nfound               ! all neighbours in mdim dimensions
 290             np=nlist(nn)
 291             if (np.lt.nnull) goto 30
 292             if (mod(np-nnull,idown2).ne.0) goto 30
 293             do 40 i=1,m-1
 294              do 41 iim=1,mdim
 295                if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
 296  41          continue
 297  40         continue
 298              call random(nseed,rr)
 299              if (rr.le.xperc) write(iunit,*)n,np
 300              if (rr.le.xperc) ntot=ntot+1
 301  30      continue
 302  20   continue
 303       return
 304       end
 305
 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)
 311       ntodo=0
 312
 313       if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "crossrec1: make nx larger."
 314       mb=min(m,2)
 315
 316       call base(nmaxy,y,id,mb,jh,jpntr,eps)
 317       nnull=(m-1)*id+1
 318       do 20 n=nnull,nmaxx,idown1
 319          if (inot(n).eq.0) goto 20
 320          ntodo=ntodo+1
 321          call neigh(nx,x,y,n,nx,id,mb,jh,jpntr,eps,nlist,nfound)
 322          if (nfound.lt.kmin) goto 20
 323          nreal=0
 324          do 30 nn=1,nfound                   ! all neighbours in two dimensions
 325             np=nlist(nn)
 326             if(np.lt.nnull) goto 30
 327             if (mod(np-nnull,idown2).ne.0) goto 30
 328             do 40 i=mb,m-1
 329                if(abs(x(n-i*id)-y(np-i*id)).ge.eps) goto 30
 330  40         continue
 331             nreal=nreal+1
 332             nlist(nreal)=np
 333  30      continue
 334          if (nreal.lt.kmin) goto 20
 335            ntodo=ntodo-1
 336            inot(n)=0
 337            ntot=ntot+nreal
 338            do in=1,nreal
 339             write(iunit,*)n,nlist(in)
 340            enddo
 341  20      continue
 342       end
 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)
 350       dimension vx(mx)
 351       ntodo=0
 352
 353       if(nmaxx.gt.nx.or.nmaxy.gt.nx) stop "mcrossrec1: make nx larger."
 354
 355       call mbase(nmaxy,mdim,nx,y,id,1,jh,jpntr,eps)
 356       nnull=(m-1)*id+1
 357       do 20 n=nnull,nmaxx,idown1
 358          if (inot(n).eq.0) goto 20
 359          ntodo=ntodo+1
 360          do ii=1,mdim
 361           vx(ii)=x(n,ii)
 362          enddo
 363          call mneigh2(nmaxy,mdim,y,nx,vx,jh,jpntr,eps,nlist,nfound)
 364          if (nfound.lt.kmin) goto 20
 365          nreal=0
 366          do 30 nn=1,nfound               ! all neighbours in mdim dimensions
 367             np=nlist(nn)
 368             if(np.lt.nnull) goto 30
 369             if (mod(np-nnull,idown2).ne.0) goto 30
 370             do 40 i=0,m-1
 371              do 41 iim=1,mdim
 372                if(abs(x(n-i*id,iim)-y(np-i*id,iim)).ge.eps) goto 30
 373  41          continue
 374  40         continue
 375             nreal=nreal+1
 376             nlist(nreal)=np
 377  30      continue
 378          if (nreal.lt.kmin) goto 20
 379            inot(n)=0
 380            ntot=ntot+nreal
 381            do in=1,nreal
 382             write(iunit,*)n,nlist(in)
 383            enddo
 384            ntodo=ntodo-1
 385  20      continue
 386       return
 387       end
 388
 389       subroutine random(iseed,s)
 390 c
 391 c     random number generator of Park & Miller
 392       integer*8 ifac,ibase,iargument
 393       ifac=7**5
 394       ibase=2**30-1
 395       im=im+2**30
 396       iargument=iseed
 397       iargument=mod(iargument*ifac,ibase)
 398       s=float(iargument)/float(ibase)
 399       iseed=iargument
 400       return
 401       end
 402 c>------------------------------------