← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/source_f/neigh.f

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/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===========================================================================
 
27
c   Modified: Piotr Held <pjheld@gmail.com> (2015). 
 
28
c   This function is based on neigh.f of TISEAN 3.0.1 https://github.com/heggus/Tisean"
 
29
c===========================================================================
 
30
      subroutine base(nmax,y,id,m,jh,jpntr,eps)
 
31
      implicit none
 
32
c     -- parameters --
 
33
      integer*4 im, ii
 
34
      parameter (im = 100, ii = 100000000) 
 
35
 
 
36
c     -- input variables --
 
37
      integer*4 nmax, id, m
 
38
      real*8 eps 
 
39
 
 
40
c     -- input arrays --
 
41
      real*8 y(nmax)
 
42
      integer*4 jh(0:im*im),jpntr(nmax)
 
43
 
 
44
c     -- variables --
 
45
      integer*4 i, n
 
46
 
 
47
      do 10 i=0,im*im
 
48
 10      jh(i)=0
 
49
      do 20 n=(m-1)*id+1,nmax                                  ! make histogram
 
50
         i=mod(int(y(n)/eps)+ii,im)
 
51
         if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
 
52
 20      jh(i)=jh(i)+1
 
53
      do 30 i=1,im*im                                           ! accumulate it
 
54
 30      jh(i)=jh(i)+jh(i-1)
 
55
      do 40 n=(m-1)*id+1,nmax                           ! fill list of pointers
 
56
 
 
57
         i=mod(int(y(n)/eps)+ii,im)
 
58
         if(m.gt.1) i=im*i+mod(int(y(n-(m-1)*id)/eps)+ii,im)
 
59
         jpntr(jh(i))=n
 
60
 40      jh(i)=jh(i)-1
 
61
      end
 
62
c>-------------------------------------------------------------
 
63
      subroutine neigh(nmax,y,x,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
 
64
      implicit none
 
65
c     -- parameters --
 
66
      integer*4 im, ii
 
67
      parameter (im = 100, ii = 100000000) 
 
68
 
 
69
c     -- input variables --
 
70
      integer*4 nmax, id, m, nfound
 
71
      real*8 eps 
 
72
 
 
73
c     -- input arrays --
 
74
      real*8 y(nmax), x(nmax)
 
75
      integer*4 jh(0:im*im),jpntr(nmax), nlist(nmax)
 
76
 
 
77
c     -- variables --
 
78
      integer*4 i, j, n, np, k, kk, kloop, jk, ip, jj, nlast
 
79
 
 
80
      nfound=0
 
81
      kloop=1
 
82
      if(m.eq.1) kloop=0
 
83
      jj=int(y(n)/eps)
 
84
 
 
85
      kk=int(y(n-(m-1)*id)/eps)
 
86
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
 
87
         do 20 k=kk-kloop,kk+kloop
 
88
            jk=mod(j+ii,im)
 
89
            if(m.gt.1) jk=im*jk+mod(k+ii,im)
 
90
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
 
91
               np=jpntr(ip)
 
92
               if(np.gt.nlast) goto 20
 
93
               do 40 i=0,m-1
 
94
 40               if(abs(y(n-i*id)-x(np-i*id)).ge.eps) goto 30
 
95
               nfound=nfound+1
 
96
               nlist(nfound)=np                       ! make list of neighbours
 
97
 30            continue
 
98
 20         continue
 
99
 10      continue
 
100
      end
 
101
 
 
102
c versions for multivariate series
 
103
c author T. Schreiber (1999)
 
104
 
 
105
      subroutine mbase(nmax,mmax,nxx,y,id,m,jh,jpntr,eps)
 
106
      implicit none
 
107
c     -- parameters --
 
108
      integer*4 im, ii
 
109
      parameter (im = 100, ii = 100000000) 
 
110
 
 
111
c     -- input variables --
 
112
      integer*4 nmax, id, m, mmax, nxx
 
113
      real*8 eps 
 
114
 
 
115
c     -- input arrays --
 
116
      real*8 y(nxx, mmax)
 
117
      integer*4 jh(0:im*im), jpntr(nmax)
 
118
 
 
119
c     -- variables --
 
120
      integer*4 i, n, mt
 
121
 
 
122
      if(mmax.eq.1) then
 
123
         call base(nmax,y,id,m,jh,jpntr,eps)
 
124
         return
 
125
      endif
 
126
      mt=(m-1)/mmax+1
 
127
      do 10 i=0,im*im
 
128
 10      jh(i)=0
 
129
      do 20 n=(mt-1)*id+1,nmax                                 ! make histogram
 
130
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 
131
 20     jh(i)=jh(i)+1
 
132
      do 30 i=1,im*im                                          ! accumulate it
 
133
 30     jh(i)=jh(i)+jh(i-1)
 
134
      do 40 n=(mt-1)*id+1,nmax                          ! fill list of pointers
 
135
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 
136
        jpntr(jh(i))=n
 
137
 40     jh(i)=jh(i)-1
 
138
      end
 
139
c>--------------------------------------------------------------
 
140
      subroutine mneigh(nmax,mmax,nxx,y,n,nlast,id,m,jh,jpntr,eps,
 
141
     .   nlist,nfound)
 
142
      implicit none
 
143
c     -- parameters --
 
144
      integer*4 im, ii
 
145
      parameter (im = 100, ii = 100000000) 
 
146
 
 
147
c     -- input variables --
 
148
      integer*4 nmax, mmax, nxx, id, m, nfound
 
149
      real*8 eps 
 
150
 
 
151
c     -- input arrays --
 
152
      real*8 y(nxx,mmax)
 
153
      integer*4 jh(0:im*im),jpntr(nmax), nlist(nmax)
 
154
 
 
155
c     -- variables --
 
156
      integer*4 i, j, n, np, k, kk, jk, ip, jj, nlast, is, mcount
 
157
      integer*4 mt
 
158
 
 
159
      if(mmax.eq.1) then
 
160
         call neigh(nmax,y,y,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
 
161
         return
 
162
      endif
 
163
      mt=(m-1)/mmax+1
 
164
      nfound=0
 
165
      jj=int(y(n,1)/eps)
 
166
      kk=int(y(n,mmax)/eps)
 
167
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
 
168
         do 20 k=kk-1,kk+1
 
169
            jk=im*mod(j+ii,im)+mod(k+ii,im)
 
170
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
 
171
               np=jpntr(ip)
 
172
               if(np.gt.nlast) goto 20
 
173
               mcount=0
 
174
               do 40 i=mt-1,0,-1
 
175
                  do 40 is=1,mmax
 
176
                     mcount=mcount+1
 
177
                     if(mcount.gt.m) goto 1
 
178
 40                  if(abs(y(n-i*id,is)-y(np-i*id,is)).ge.eps) goto 30
 
179
 1             nfound=nfound+1
 
180
               nlist(nfound)=np                       ! make list of neighbours
 
181
 30            continue
 
182
 20         continue
 
183
 10      continue
 
184
      end
 
185
c>---------------------------------------------------------------------
 
186
c modified version for multivariate series
 
187
c author H. Kantz (2004)
 
188
 
 
189
      subroutine mneigh2(nmax,mdim,y,nx,vx,jh,jpntr,eps,
 
190
     .   nlist,nfound)
 
191
c
 
192
c     search neighbours for vx among the set of all y's
 
193
c     multivariate: mmax: spatial dimension
 
194
c     no additional delay!
 
195
      implicit none
 
196
c     -- parameters --
 
197
      integer*4 im, ii
 
198
      parameter (im = 100, ii = 100000000) 
 
199
 
 
200
c     -- input variables --
 
201
      integer*4 nmax, mdim, nx, nfound
 
202
      real*8 eps 
 
203
 
 
204
c     -- input arrays --
 
205
      real*8  y(nx,mdim), vx(mdim)
 
206
      integer*4 jh(0:im*im), jpntr(nmax), nlist(nmax)
 
207
 
 
208
c     -- variables --
 
209
      integer*4 j, np, k, kk, jk, ip, jj, is, mcount
 
210
 
 
211
      nfound=0
 
212
      jj=int(vx(1)/eps)
 
213
      kk=int(vx(mdim)/eps)
 
214
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
 
215
         do 20 k=kk-1,kk+1
 
216
            jk=im*mod(j+ii,im)+mod(k+ii,im)
 
217
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
 
218
               np=jpntr(ip)
 
219
c               if(np.gt.nlast) goto 20
 
220
               mcount=0
 
221
                  do 40 is=1,mdim
 
222
 40                  if(abs(vx(is)-y(np,is)).ge.eps) goto 30
 
223
               nfound=nfound+1
 
224
               nlist(nfound)=np                       ! make list of neighbours
 
225
 30            continue
 
226
 20         continue
 
227
 10      continue
 
228
      end
 
229
c>---------------------------------------------------------------------
 
230
      subroutine mbase2(nmax,mmax,nxx,y,jh,jpntr,eps)
 
231
 
 
232
      implicit none
 
233
c     -- parameters --
 
234
      integer*4 im, ii
 
235
      parameter (im = 100, ii = 100000000) 
 
236
 
 
237
c     -- input variables --
 
238
      integer*4 nmax, id, m, mmax, nxx
 
239
      real*8 eps 
 
240
 
 
241
c     -- input arrays --
 
242
      real*8 y(nxx, mmax)
 
243
      integer*4 jh(0:im*im), jpntr(nmax)
 
244
 
 
245
c     -- variables --
 
246
      integer*4 i, n
 
247
 
 
248
      if(mmax.eq.1) then
 
249
         call base(nmax,y,id,m,jh,jpntr,eps)
 
250
         return
 
251
      endif
 
252
      do 10 i=0,im*im
 
253
 10      jh(i)=0
 
254
      do 20 n=1,nmax                                 ! make histogram
 
255
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 
256
 20     jh(i)=jh(i)+1
 
257
      do 30 i=1,im*im                                          ! accumulate it
 
258
 30     jh(i)=jh(i)+jh(i-1)
 
259
      do 40 n=(mmax-1)*id+1,nmax                      ! fill list of pointers
 
260
        i=im*mod(int(y(n,1)/eps)+ii,im)+mod(int(y(n,mmax)/eps)+ii,im)
 
261
        jpntr(jh(i))=n
 
262
 40     jh(i)=jh(i)-1
 
263
      end