website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/neigh.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   utilities for neighbour search
  23 c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
  24 c      University Press (1997)
  25 c   author T. Schreiber (1999)
  26 c   last modified H. Kantz Feb.2007
  27 c===========================================================================
  28       subroutine base(nmax,y,id,m,jh,jpntr,eps)
  29       parameter(im=100,ii=100000000)
  30       dimension y(nmax),jh(0:im*im),jpntr(nmax)
  31
  32       do 10 i=0,im*im
  33  10      jh(i)=0
  34       do 20 n=(m-1)*id+1,nmax                                  ! make histogram
  35          i=mod(int(y(n)/eps)+ii,im)
  36          if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
  37  20      jh(i)=jh(i)+1
  38       do 30 i=1,im*im                                           ! accumulate it
  39  30      jh(i)=jh(i)+jh(i-1)
  40       do 40 n=(m-1)*id+1,nmax                           ! fill list of pointers
  41
  42          i=mod(int(y(n)/eps)+ii,im)
  43          if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
  44          jpntr(jh(i))=n
  45  40      jh(i)=jh(i)-1
  46       end
  47
  48       subroutine neigh(nmax,y,x,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
  49       parameter(im=100,ii=100000000)
  50       dimension y(nmax),x(nmax),jh(0:im*im),jpntr(nmax),nlist(nmax)
  51
  52       nfound=0
  53       kloop=1
  54       if(m.eq.1) kloop=0
  55       jj=int(y(n)/eps)
  56
  57       kk=int(y(n-(m-1)*id)/eps)
  58       do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
  59          do 20 k=kk-kloop,kk+kloop
  60             jk=mod(j+ii,im)
  61             if(m.gt.1) jk=im*jk+mod(k+ii,im)
  62             do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
  63                np=jpntr(ip)
  64                if(np.gt.nlast) goto 20
  65                do 40 i=0,m-1
  66  40               if(abs(y(n-i*id)-x(np-i*id)).ge.eps) goto 30
  67                nfound=nfound+1
  68                nlist(nfound)=np                       ! make list of neighbours
  69  30            continue
  70  20         continue
  71  10      continue
  72       end
  73
  74 c versions for multivariate series
  75 c author T. Schreiber (1999)
  76
  77       subroutine mbase(nmax,mmax,nxx,y,id,m,jh,jpntr,eps)
  78       parameter(im=100,ii=100000000)
  79       dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax)
  80
  81       if(mmax.eq.1) then
  82          call base(nmax,y,id,m,jh,jpntr,eps)
  83          return
  84       endif
  85       mt=(m-1)/mmax+1
  86       do 10 i=0,im*im
  87  10      jh(i)=0
  88       do 20 n=(mt-1)*id+1,nmax                                 ! make histogram
  89         i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
  90  20     jh(i)=jh(i)+1
  91       do 30 i=1,im*im                                          ! accumulate it
  92  30     jh(i)=jh(i)+jh(i-1)
  93       do 40 n=(mt-1)*id+1,nmax                          ! fill list of pointers
  94         i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
  95         jpntr(jh(i))=n
  96  40     jh(i)=jh(i)-1
  97       end
  98
  99       subroutine mneigh(nmax,mmax,nxx,y,n,nlast,id,m,jh,jpntr,eps,
 100      .   nlist,nfound)
 101       parameter(im=100,ii=100000000)
 102       dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax),nlist(nmax)
 103
 104       if(mmax.eq.1) then
 105          call neigh(nmax,y,y,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
 106          return
 107       endif
 108       mt=(m-1)/mmax+1
 109       nfound=0
 110       jj=int(y(n,1)/eps)
 111       kk=int(y(n,mmax)/eps)
 112       do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
 113          do 20 k=kk-1,kk+1
 114             jk=im*mod(j+ii,im)+mod(k+ii,im)
 115             do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
 116                np=jpntr(ip)
 117                if(np.gt.nlast) goto 20
 118                mcount=0
 119                do 40 i=mt-1,0,-1
 120                   do 40 is=1,mmax
 121                      mcount=mcount+1
 122                      if(mcount.gt.m) goto 1
 123  40                  if(abs(y(n-i*id,is)-y(np-i*id,is)).ge.eps) goto 30
 124  1             nfound=nfound+1
 125                nlist(nfound)=np                       ! make list of neighbours
 126  30            continue
 127  20         continue
 128  10      continue
 129       end
 130 c>---------------------------------------------------------------------
 131 c modified version for multivariate series
 132 c author H. Kantz (2004)
 133
 134       subroutine mneigh2(nmax,mdim,y,nx,vx,jh,jpntr,eps,
 135      .   nlist,nfound)
 136 c
 137 c     search neighbours for vx among the set of all y's
 138 c     multivariate: mmax: spatial dimension
 139 c     no additional delay!
 140       parameter(im=100,ii=100000000)
 141       dimension y(nx,mdim),jh(0:im*im),jpntr(nmax),nlist(nmax)
 142       dimension vx(mdim)
 143
 144       nfound=0
 145       jj=int(vx(1)/eps)
 146       kk=int(vx(mdim)/eps)
 147       do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
 148          do 20 k=kk-1,kk+1
 149             jk=im*mod(j+ii,im)+mod(k+ii,im)
 150             do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
 151                np=jpntr(ip)
 152 c               if(np.gt.nlast) goto 20
 153                mcount=0
 154                   do 40 is=1,mdim
 155  40                  if(abs(vx(is)-y(np,is)).ge.eps) goto 30
 156  1             nfound=nfound+1
 157                nlist(nfound)=np                       ! make list of neighbours
 158  30            continue
 159  20         continue
 160  10      continue
 161       end
 162
 163       subroutine mbase2(nmax,mmax,nxx,y,jh,jpntr,eps)
 164       parameter(im=100,ii=100000000)
 165       dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax)
 166
 167       if(mmax.eq.1) then
 168          call base(nmax,y,id,m,jh,jpntr,eps)
 169          return
 170       endif
 171       do 10 i=0,im*im
 172  10      jh(i)=0
 173       do 20 n=1,nmax                                 ! make histogram
 174         i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 175  20     jh(i)=jh(i)+1
 176       do 30 i=1,im*im                                          ! accumulate it
 177  30     jh(i)=jh(i)+jh(i-1)
 178       do 40 n=(mmax-1)*id+1,nmax                      ! fill list of pointers
 179         i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 180         jpntr(jh(i))=n
 181  40     jh(i)=jh(i)-1
 182       end