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   utilities for neighbour search
c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
c      University Press (1997)
c   author T. Schreiber (1999)
c   last modified H. Kantz Feb.2007
c===========================================================================
      subroutine base(nmax,y,id,m,jh,jpntr,eps)
      parameter(im=100,ii=100000000) 
      dimension y(nmax),jh(0:im*im),jpntr(nmax)

      do 10 i=0,im*im
 10      jh(i)=0
      do 20 n=(m-1)*id+1,nmax                                  ! make histogram
         i=mod(int(y(n)/eps)+ii,im)
         if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
 20      jh(i)=jh(i)+1
      do 30 i=1,im*im                                           ! accumulate it
 30      jh(i)=jh(i)+jh(i-1)
      do 40 n=(m-1)*id+1,nmax                           ! fill list of pointers

         i=mod(int(y(n)/eps)+ii,im)
         if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
         jpntr(jh(i))=n
 40      jh(i)=jh(i)-1
      end

      subroutine neigh(nmax,y,x,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
      parameter(im=100,ii=100000000) 
      dimension y(nmax),x(nmax),jh(0:im*im),jpntr(nmax),nlist(nmax)

      nfound=0
      kloop=1
      if(m.eq.1) kloop=0
      jj=int(y(n)/eps)

      kk=int(y(n-(m-1)*id)/eps)
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
         do 20 k=kk-kloop,kk+kloop
            jk=mod(j+ii,im)
            if(m.gt.1) jk=im*jk+mod(k+ii,im)
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
               np=jpntr(ip)
               if(np.gt.nlast) goto 20
               do 40 i=0,m-1
 40               if(abs(y(n-i*id)-x(np-i*id)).ge.eps) goto 30
               nfound=nfound+1
               nlist(nfound)=np                       ! make list of neighbours
 30            continue
 20         continue
 10      continue
      end

c versions for multivariate series
c author T. Schreiber (1999)

      subroutine mbase(nmax,mmax,nxx,y,id,m,jh,jpntr,eps)
      parameter(im=100,ii=100000000) 
      dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax)

      if(mmax.eq.1) then
         call base(nmax,y,id,m,jh,jpntr,eps)
         return
      endif
      mt=(m-1)/mmax+1
      do 10 i=0,im*im
 10      jh(i)=0
      do 20 n=(mt-1)*id+1,nmax                                 ! make histogram
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 20     jh(i)=jh(i)+1
      do 30 i=1,im*im                                          ! accumulate it
 30     jh(i)=jh(i)+jh(i-1)
      do 40 n=(mt-1)*id+1,nmax                          ! fill list of pointers
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
        jpntr(jh(i))=n
 40     jh(i)=jh(i)-1
      end

      subroutine mneigh(nmax,mmax,nxx,y,n,nlast,id,m,jh,jpntr,eps,
     .   nlist,nfound)
      parameter(im=100,ii=100000000) 
      dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax),nlist(nmax)

      if(mmax.eq.1) then
         call neigh(nmax,y,y,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
         return
      endif
      mt=(m-1)/mmax+1
      nfound=0
      jj=int(y(n,1)/eps)
      kk=int(y(n,mmax)/eps)
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
         do 20 k=kk-1,kk+1
            jk=im*mod(j+ii,im)+mod(k+ii,im)
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
               np=jpntr(ip)
               if(np.gt.nlast) goto 20
               mcount=0
               do 40 i=mt-1,0,-1
                  do 40 is=1,mmax
                     mcount=mcount+1
                     if(mcount.gt.m) goto 1
 40                  if(abs(y(n-i*id,is)-y(np-i*id,is)).ge.eps) goto 30
 1             nfound=nfound+1
               nlist(nfound)=np                       ! make list of neighbours
 30            continue
 20         continue
 10      continue
      end
c>---------------------------------------------------------------------
c modified version for multivariate series
c author H. Kantz (2004)

      subroutine mneigh2(nmax,mdim,y,nx,vx,jh,jpntr,eps,
     .   nlist,nfound)
c
c     search neighbours for vx among the set of all y's
c     multivariate: mmax: spatial dimension
c     no additional delay!
      parameter(im=100,ii=100000000) 
      dimension y(nx,mdim),jh(0:im*im),jpntr(nmax),nlist(nmax)
      dimension vx(mdim)

      nfound=0
      jj=int(vx(1)/eps)
      kk=int(vx(mdim)/eps)
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
         do 20 k=kk-1,kk+1
            jk=im*mod(j+ii,im)+mod(k+ii,im)
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
               np=jpntr(ip)
c               if(np.gt.nlast) goto 20
               mcount=0
                  do 40 is=1,mdim
 40                  if(abs(vx(is)-y(np,is)).ge.eps) goto 30
 1             nfound=nfound+1
               nlist(nfound)=np                       ! make list of neighbours
 30            continue
 20         continue
 10      continue
      end

      subroutine mbase2(nmax,mmax,nxx,y,jh,jpntr,eps)
      parameter(im=100,ii=100000000) 
      dimension y(nxx,mmax),jh(0:im*im),jpntr(nmax)

      if(mmax.eq.1) then
         call base(nmax,y,id,m,jh,jpntr,eps)
         return
      endif
      do 10 i=0,im*im
 10      jh(i)=0
      do 20 n=1,nmax                                 ! make histogram
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 20     jh(i)=jh(i)+1
      do 30 i=1,im*im                                          ! accumulate it
 30     jh(i)=jh(i)+jh(i-1)
      do 40 n=(mmax-1)*id+1,nmax                      ! fill list of pointers
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
        jpntr(jh(i))=n
 40     jh(i)=jh(i)-1
      end