← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/sandbox/pyloess/src/stl.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sandbox/pyloess/src/stl.f
 
1
      subroutine stl(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,
 
2
     &nljump,ni,no,rw,season,trend,work)
 
3
      integer n, np, ns, nt, nl, isdeg, itdeg, ildeg, nsjump, ntjump, 
 
4
     &nljump, ni, no, k
 
5
      integer newns, newnt, newnl, newnp
 
6
      double precision y(n), rw(n), season(n), trend(n), work(n+2*np,5)
 
7
      logical userw
 
8
      userw = .false.
 
9
      k = 0
 
10
      do 23000 i = 1,n
 
11
          trend(i) = 0.D0
 
12
23000 continue
 
13
      newns = max(3,ns)
 
14
      newnt = max(3,nt)
 
15
      newnl = max(3,nl)
 
16
      newnp = max(2,np)
 
17
      if(.not.(mod(newns,2) .eq. 0))goto 23002
 
18
      newns = newns + 1
 
19
23002 continue
 
20
      if(.not.(mod(newnt,2) .eq. 0))goto 23004
 
21
      newnt = newnt + 1
 
22
23004 continue
 
23
      if(.not.(mod(newnl,2) .eq. 0))goto 23006
 
24
      newnl = newnl + 1
 
25
23006 continue
 
26
23008 continue
 
27
      call onestp(y,n,newnp,newns,newnt,newnl,isdeg,itdeg,ildeg,nsjump,
 
28
     &ntjump,nljump,ni,userw,rw,season, trend, work)
 
29
      k = k+1
 
30
      if(.not.(k .gt. no))goto 23011
 
31
      goto 23010
 
32
23011 continue
 
33
      do 23013 i = 1,n
 
34
          work(i,1) = trend(i)+season(i)
 
35
23013 continue
 
36
      call rwts(y,n,work(1,1),rw)
 
37
      userw = .true.
 
38
23009 goto 23008
 
39
23010 continue
 
40
      if(.not.(no .le. 0))goto 23015
 
41
      do 23017 i = 1,n
 
42
          rw(i) = 1.D0
 
43
23017 continue
 
44
23015 continue
 
45
      return
 
46
      end
 
47
      
 
48
      
 
49
      subroutine ess(y,n,len,ideg,njump,userw,rw,ys,res)
 
50
      integer n, len, ideg, njump, newnj, nleft, nright, nsh, k, i, j
 
51
      double precision y(n), rw(n), ys(n), res(n), delta
 
52
      logical ok, userw
 
53
      if(.not.(n .lt. 2))goto 23019
 
54
      ys(1) = y(1)
 
55
      return
 
56
23019 continue
 
57
      newnj = min(njump, n-1)
 
58
      if(.not.(len .ge. n))goto 23021
 
59
      nleft = 1
 
60
      nright = n
 
61
      do 23023 i = 1,n,newnj 
 
62
      call est(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,userw,rw,ok)
 
63
      if(.not.( .not. ok))goto 23025
 
64
      ys(i) = y(i)
 
65
23025 continue
 
66
23023 continue
 
67
      goto 23022
 
68
23021 continue
 
69
      if(.not.(newnj .eq. 1))goto 23027
 
70
      nsh = (len+1)/2
 
71
      nleft = 1
 
72
      nright = len
 
73
      do 23029 i = 1,n 
 
74
      if(.not.(i .gt. nsh  .and.  nright .ne. n))goto 23031
 
75
      nleft = nleft+1
 
76
      nright = nright+1
 
77
23031 continue
 
78
      call est(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,userw,rw,ok)
 
79
      if(.not.( .not. ok))goto 23033
 
80
      ys(i) = y(i)
 
81
23033 continue
 
82
23029 continue
 
83
      goto 23028
 
84
23027 continue
 
85
      nsh = (len+1)/2
 
86
      do 23035 i = 1,n,newnj 
 
87
      if(.not.(i .lt. nsh))goto 23037
 
88
      nleft = 1
 
89
      nright = len
 
90
      goto 23038
 
91
23037 continue
 
92
      if(.not.(i .ge. n-nsh+1))goto 23039
 
93
      nleft = n-len+1
 
94
      nright = n
 
95
      goto 23040
 
96
23039 continue
 
97
      nleft = i-nsh+1
 
98
      nright = len+i-nsh
 
99
23040 continue
 
100
23038 continue
 
101
      call est(y,n,len,ideg,dble(i),ys(i),nleft,nright,res,userw,rw,ok)
 
102
      if(.not.( .not. ok))goto 23041
 
103
      ys(i) = y(i)
 
104
23041 continue
 
105
23035 continue
 
106
23028 continue
 
107
23022 continue
 
108
      if(.not.(newnj .ne. 1))goto 23043
 
109
      do 23045 i = 1,n-newnj,newnj 
 
110
      delta = (ys(i+newnj)-ys(i))/dble(newnj)
 
111
      do 23047 j = i+1,i+newnj-1
 
112
      ys(j) = ys(i)+delta*dble(j-i)
 
113
23047 continue
 
114
23045 continue
 
115
      k = ((n-1)/newnj)*newnj+1
 
116
      if(.not.(k .ne. n))goto 23049
 
117
      call est(y,n,len,ideg,dble(n),ys(n),nleft,nright,res,userw,rw,ok)
 
118
      if(.not.( .not. ok))goto 23051
 
119
      ys(n) = y(n)
 
120
23051 continue
 
121
      if(.not.(k .ne. n-1))goto 23053
 
122
      delta = (ys(n)-ys(k))/dble(n-k)
 
123
      do 23055 j = k+1,n-1
 
124
      ys(j) = ys(k)+delta*dble(j-k)
 
125
23055 continue
 
126
23053 continue
 
127
23049 continue
 
128
23043 continue
 
129
      return
 
130
      end
 
131
      
 
132
      
 
133
      subroutine est(y,n,len,ideg,xs,ys,nleft,nright,w,userw,rw,ok)
 
134
      integer n, len, ideg, nleft, nright, j
 
135
      double precision y(n), w(n), rw(n), xs, ys, range, h, h1, h9, a, 
 
136
     &b, c, r
 
137
      logical userw,ok
 
138
      range = dble(n)-dble(1)
 
139
      h = max(xs-dble(nleft),dble(nright)-xs)
 
140
      if(.not.(len .gt. n))goto 23057
 
141
      h = h+dble((len-n)/2)
 
142
23057 continue
 
143
      h9 = .999999D0*h
 
144
      h1 = .000001D0*h
 
145
      a = 0.D0
 
146
      do 23059 j = nleft,nright 
 
147
      w(j) = 0.D0
 
148
      r = dabs(dble(j)-xs)
 
149
      if(.not.(r .le. h9))goto 23061
 
150
      if(.not.(r .le. h1))goto 23063
 
151
      w(j) = 1.D0
 
152
      goto 23064
 
153
23063 continue
 
154
      w(j) = (1.D0-(r/h)**3.D0)**3.D0
 
155
23064 continue
 
156
      if(.not.(userw))goto 23065
 
157
      w(j) = rw(j)*w(j)
 
158
23065 continue
 
159
      a = a+w(j)
 
160
23061 continue
 
161
23059 continue
 
162
      if(.not.(a .le. 0.D0))goto 23067
 
163
      ok = .false.
 
164
      goto 23068
 
165
23067 continue
 
166
      ok = .true.
 
167
      do 23069 j = nleft,nright
 
168
      w(j) = w(j)/a
 
169
23069 continue
 
170
      if(.not.((h .gt. 0.D0) .and. (ideg .gt. 0)))goto 23071
 
171
      a = 0.D0
 
172
      do 23073 j = nleft,nright
 
173
      a = a+w(j)*dble(j)
 
174
23073 continue
 
175
      b = xs-a
 
176
      c = 0.D0
 
177
      do 23075 j = nleft,nright
 
178
      c = c+w(j)*(dble(j)-a)**2.D0
 
179
23075 continue
 
180
      if(.not.(sqrt(c) .gt. .001*range))goto 23077
 
181
      b = b/c
 
182
      do 23079 j = nleft,nright
 
183
      w(j) = w(j)*(b*(dble(j)-a)+1.D0)
 
184
23079 continue
 
185
23077 continue
 
186
23071 continue
 
187
      ys = 0.D0
 
188
      do 23081 j = nleft,nright
 
189
      ys = ys+w(j)*y(j)
 
190
23081 continue
 
191
23068 continue
 
192
      return
 
193
      end
 
194
      
 
195
      
 
196
      subroutine fts(x,n,np,trend,work)
 
197
      integer n, np
 
198
      double precision x(n), trend(n), work(n)
 
199
      call ma(x,n,np,trend)
 
200
      call ma(trend,n-np+1,np,work)
 
201
      call ma(work,n-2*np+2,3,trend)
 
202
      return
 
203
      end
 
204
      
 
205
      
 
206
      subroutine ma(x, n, len, ave)
 
207
      integer n, len, i, j, k, m, newn
 
208
      double precision x(n), ave(n), flen, v
 
209
      newn = n-len+1
 
210
      flen = dble(len)
 
211
      v = 0.D0
 
212
      do 23083 i = 1,len
 
213
      v = v+x(i)
 
214
23083 continue
 
215
      ave(1) = v/flen
 
216
      if(.not.(newn .gt. 1))goto 23085
 
217
      k = len
 
218
      m = 0
 
219
      do 23087 j = 2, newn 
 
220
      k = k+1
 
221
      m = m+1
 
222
      v = v-x(m)+x(k)
 
223
      ave(j) = v/flen
 
224
23087 continue
 
225
23085 continue
 
226
      return
 
227
      end
 
228
      
 
229
      
 
230
      subroutine onestp(y,n,np,ns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,
 
231
     &nljump,ni,userw,rw,season,trend,work)
 
232
      integer n,ni,np,ns,nt,nsjump,ntjump,nl,nljump,isdeg,itdeg,ildeg
 
233
      double precision y(n),rw(n),season(n),trend(n),work(n+2*np,5)
 
234
      logical userw
 
235
      do 23089 j = 1,ni 
 
236
      do 23091 i = 1,n
 
237
      work(i,1) = y(i)-trend(i)
 
238
23091 continue
 
239
      call ss(work(1,1),n,np,ns,isdeg,nsjump,userw,rw,work(1,2),work(1,
 
240
     &3),work(1,4),work(1,5),season)
 
241
      call fts(work(1,2),n+2*np,np,work(1,3),work(1,1))
 
242
      call ess(work(1,3),n,nl,ildeg,nljump,.false.,work(1,4),work(1,1),
 
243
     &work(1,5))
 
244
      do 23093 i = 1,n
 
245
      season(i) = work(np+i,2)-work(i,1)
 
246
23093 continue
 
247
      do 23095 i = 1,n
 
248
      work(i,1) = y(i)-season(i)
 
249
23095 continue
 
250
      call ess(work(1,1),n,nt,itdeg,ntjump,userw,rw,trend,work(1,3))
 
251
23089 continue
 
252
      return
 
253
      end
 
254
      
 
255
      
 
256
      subroutine rwts(y,n,fit,rw)
 
257
      integer mid(2), n
 
258
      double precision y(n), fit(n), rw(n), cmad, c9, c1, r
 
259
      do 23097 i = 1,n
 
260
      rw(i) = dabs(y(i)-fit(i))
 
261
23097 continue
 
262
      mid(1) = n/2+1
 
263
      mid(2) = n-mid(1)+1
 
264
      call psort(rw,n,mid,2)
 
265
      cmad = 3.D0*(rw(mid(1))+rw(mid(2)))
 
266
      c9 = .999999D0*cmad
 
267
      c1 = .000001D0*cmad
 
268
      do 23099 i = 1,n 
 
269
      r = dabs(y(i)-fit(i))
 
270
      if(.not.(r .le. c1))goto 23101
 
271
      rw(i) = 1.
 
272
      goto 23102
 
273
23101 continue
 
274
      if(.not.(r .le. c9))goto 23103
 
275
      rw(i) = (1.D0-(r/cmad)**2.D0)**2.D0
 
276
      goto 23104
 
277
23103 continue
 
278
      rw(i) = 0.D0
 
279
23104 continue
 
280
23102 continue
 
281
23099 continue
 
282
      return
 
283
      end
 
284
      
 
285
      
 
286
      subroutine ss(y,n,np,ns,isdeg,nsjump,userw,rw,season,work1,work2,
 
287
     &work3,work4)
 
288
      integer n, np, ns, isdeg, nsjump, nright, nleft, i, j, k
 
289
      double precision y(n), rw(n), season(n+2*np), work1(n), work2(n), 
 
290
     &work3(n), work4(n), xs
 
291
      logical userw,ok
 
292
      j=1
 
293
23105 if(.not.(j .le. np))goto 23107
 
294
      k = (n-j)/np+1
 
295
      do 23108 i = 1,k
 
296
      work1(i) = y((i-1)*np+j)
 
297
23108 continue
 
298
      if(.not.(userw))goto 23110
 
299
      do 23112 i = 1,k
 
300
      work3(i) = rw((i-1)*np+j)
 
301
23112 continue
 
302
23110 continue
 
303
      call ess(work1,k,ns,isdeg,nsjump,userw,work3,work2(2),work4)
 
304
      xs = 0.D0
 
305
      nright = min0(ns,k)
 
306
      call est(work1,k,ns,isdeg,xs,work2(1),1,nright,work4,userw,work3,
 
307
     &ok)
 
308
      if(.not.( .not. ok))goto 23114
 
309
      work2(1) = work2(2)
 
310
23114 continue
 
311
      xs = dble(k+1)
 
312
      nleft = max0(1,k-ns+1)
 
313
      call est(work1,k,ns,isdeg,xs,work2(k+2),nleft,k,work4,userw,work3,
 
314
     &ok)
 
315
      if(.not.( .not. ok))goto 23116
 
316
      work2(k+2) = work2(k+1)
 
317
23116 continue
 
318
      do 23118 m = 1,k+2
 
319
      season((m-1)*np+j) = work2(m)
 
320
23118 continue
 
321
       j=j+1
 
322
      goto 23105
 
323
23107 continue
 
324
      return
 
325
      end
 
326
      
 
327
      
 
328
      subroutine stlez(y, n, np, ns, isdeg, itdeg, robust, no, rw, 
 
329
     &season, trend, work)
 
330
      logical robust
 
331
      integer n, i, j, np, ns, no, nt, nl, ni, nsjump, ntjump, nljump, 
 
332
     &newns, newnp
 
333
      integer isdeg, itdeg, ildeg
 
334
      double precision y(n), rw(n), season(n), trend(n), work(n+2*np,7)
 
335
      double precision maxs, mins, maxt, mint, maxds, maxdt, difs, dift
 
336
      ildeg = itdeg
 
337
      newns = max(3,ns)
 
338
      if(.not.(mod(newns,2) .eq. 0))goto 23120
 
339
      newns = newns+1
 
340
23120 continue
 
341
      newnp = max(2,np)
 
342
      nt = (1.5*newnp)/(1 - 1.5/newns) + 0.5
 
343
      nt = max(3,nt)
 
344
      if(.not.(mod(nt,2) .eq. 0))goto 23122
 
345
      nt = nt+1
 
346
23122 continue
 
347
      nl = newnp
 
348
      if(.not.(mod(nl,2) .eq. 0))goto 23124
 
349
      nl = nl+1
 
350
23124 continue
 
351
      if(.not.(robust))goto 23126
 
352
      ni = 1
 
353
      goto 23127
 
354
23126 continue
 
355
      ni = 2
 
356
23127 continue
 
357
      nsjump = max(1,int(dble(newns)/10.D0 + 0.9D0))
 
358
      ntjump = max(1,int(dble(nt)/10.D0 + 0.9D0))
 
359
      nljump = max(1,int(dble(nl)/10.D0 + 0.9D0))
 
360
      do 23128 i = 1,n
 
361
      trend(i) = 0.D0
 
362
23128 continue
 
363
      call onestp(y,n,newnp,newns,nt,nl,isdeg,itdeg,ildeg,nsjump,ntjump,
 
364
     &nljump,ni,.false.,rw,season,trend,work)
 
365
      no = 0
 
366
      if(.not.(robust))goto 23130
 
367
      j=1
 
368
23132 if(.not.(j .le. 15))goto 23134
 
369
      do 23135 i = 1,n
 
370
      work(i,6) = season(i)
 
371
      work(i,7) = trend(i)
 
372
      work(i,1) = trend(i)+season(i)
 
373
23135 continue
 
374
      call rwts(y,n,work(1,1),rw)
 
375
      call onestp(y, n, newnp, newns, nt, nl, isdeg, itdeg, ildeg, 
 
376
     &nsjump,ntjump, nljump, ni, .true., rw, season, trend, work)
 
377
      no = no+1
 
378
      maxs = work(1,6)
 
379
      mins = work(1,6)
 
380
      maxt = work(1,7)
 
381
      mint = work(1,7)
 
382
      maxds = dabs(work(1,6) - season(1))
 
383
      maxdt = dabs(work(1,7) - trend(1))
 
384
      do 23137 i = 2,n
 
385
      if(.not.(maxs .lt. work(i,6)))goto 23139
 
386
      maxs = work(i,6)
 
387
23139 continue
 
388
      if(.not.(maxt .lt. work(i,7)))goto 23141
 
389
      maxt = work(i,7)
 
390
23141 continue
 
391
      if(.not.(mins .gt. work(i,6)))goto 23143
 
392
      mins = work(i,6)
 
393
23143 continue
 
394
      if(.not.(mint .gt. work(i,7)))goto 23145
 
395
      mint = work(i,7)
 
396
23145 continue
 
397
      difs = dabs(work(i,6) - season(i))
 
398
      dift = dabs(work(i,7) - trend(i))
 
399
      if(.not.(maxds .lt. difs))goto 23147
 
400
      maxds = difs
 
401
23147 continue
 
402
      if(.not.(maxdt .lt. dift))goto 23149
 
403
      maxdt = dift
 
404
23149 continue
 
405
23137 continue
 
406
      if(.not.((maxds/(maxs-mins) .lt. .01)  .and.  
 
407
     &    (maxdt/(maxt-mint) .lt. .01)))goto 23151
 
408
      goto 23134
 
409
23151 continue
 
410
       j=j+1
 
411
      goto 23132
 
412
23134 continue
 
413
23130 continue
 
414
      if(.not.( .not. robust))goto 23153
 
415
      do 23155 i = 1,n
 
416
      rw(i) = 1.D0
 
417
23155 continue
 
418
23153 continue
 
419
      return
 
420
      end
 
421
      
 
422
      
 
423
      subroutine psort(a,n,ind,ni)
 
424
      double precision a(n)
 
425
      integer n,ind(ni),ni
 
426
      integer indu(16),indl(16),iu(16),il(16),p,jl,ju,i,j,m,k,ij,l
 
427
      double precision t,tt
 
428
      if(.not.(n .lt. 0 .or. ni .lt. 0))goto 23157
 
429
      return
 
430
23157 continue
 
431
      if(.not.(n .lt. 2  .or.  ni .eq. 0))goto 23159
 
432
      return
 
433
23159 continue
 
434
      jl = 1
 
435
      ju = ni
 
436
      indl(1) = 1
 
437
      indu(1) = ni
 
438
      i = 1
 
439
      j = n
 
440
      m = 1
 
441
23161 continue
 
442
      if(.not.(i .lt. j))goto 23164
 
443
      go to 10
 
444
23164 continue
 
445
23166 continue
 
446
      m = m-1
 
447
      if(.not.(m .eq. 0))goto 23169
 
448
      goto 23163
 
449
23169 continue
 
450
      i = il(m)
 
451
      j = iu(m)
 
452
      jl = indl(m)
 
453
      ju = indu(m)
 
454
      if(.not.(jl .le. ju))goto 23171
 
455
23173 if(.not.(j-i .gt. 10))goto 23174
 
456
10    k = i
 
457
      ij = (i+j)/2
 
458
      t = a(ij)
 
459
      if(.not.(a(i) .gt. t))goto 23175
 
460
      a(ij) = a(i)
 
461
      a(i) = t
 
462
      t = a(ij)
 
463
23175 continue
 
464
      l = j
 
465
      if(.not.(a(j) .lt. t))goto 23177
 
466
      a(ij) = a(j)
 
467
      a(j) = t
 
468
      t = a(ij)
 
469
      if(.not.(a(i) .gt. t))goto 23179
 
470
      a(ij) = a(i)
 
471
      a(i) = t
 
472
      t = a(ij)
 
473
23179 continue
 
474
23177 continue
 
475
23181 continue
 
476
      l = l-1
 
477
      if(.not.(a(l) .le. t))goto 23184
 
478
      tt = a(l)
 
479
23186 continue
 
480
      k = k+1
 
481
23187 if(.not.(a(k) .ge. t))goto 23186
 
482
      if(.not.(k .gt. l))goto 23189
 
483
      goto 23183
 
484
23189 continue
 
485
      a(l) = a(k)
 
486
      a(k) = tt
 
487
23184 continue
 
488
23182 goto 23181
 
489
23183 continue
 
490
      indl(m) = jl
 
491
      indu(m) = ju
 
492
      p = m
 
493
      m = m+1
 
494
      if(.not.(l-i .le. j-k))goto 23191
 
495
      il(p) = k
 
496
      iu(p) = j
 
497
      j = l
 
498
23193 continue
 
499
      if(.not.(jl .gt. ju))goto 23196
 
500
      goto 23167
 
501
23196 continue
 
502
      if(.not.(ind(ju) .le. j))goto 23198
 
503
      goto 23195
 
504
23198 continue
 
505
      ju = ju-1
 
506
23194 goto 23193
 
507
23195 continue
 
508
      indl(p) = ju+1
 
509
      goto 23192
 
510
23191 continue
 
511
      il(p) = i
 
512
      iu(p) = l
 
513
      i = k
 
514
23200 continue
 
515
      if(.not.(jl .gt. ju))goto 23203
 
516
      goto 23167
 
517
23203 continue
 
518
      if(.not.(ind(jl) .ge. i))goto 23205
 
519
      goto 23202
 
520
23205 continue
 
521
      jl = jl+1
 
522
23201 goto 23200
 
523
23202 continue
 
524
      indu(p) = jl-1
 
525
23192 continue
 
526
      goto 23173
 
527
23174 continue
 
528
      if(.not.(i .eq. 1))goto 23207
 
529
      goto 23168
 
530
23207 continue
 
531
      i = i-1
 
532
23209 continue
 
533
      i = i+1
 
534
      if(.not.(i .eq. j))goto 23212
 
535
      goto 23211
 
536
23212 continue
 
537
      t = a(i+1)
 
538
      if(.not.(a(i) .gt. t))goto 23214
 
539
      k = i
 
540
23216 continue
 
541
      a(k+1) = a(k)
 
542
      k = k-1
 
543
23217 if(.not.(t .ge. a(k)))goto 23216
 
544
      a(k+1) = t
 
545
23214 continue
 
546
23210 goto 23209
 
547
23211 continue
 
548
23171 continue
 
549
23167 goto 23166
 
550
23168 continue
 
551
23162 goto 23161
 
552
23163 continue
 
553
      return
 
554
      end