← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/interpolate/fitpack/regrid.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/interpolate/fitpack/regrid.f
 
1
      subroutine regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,
 
2
     * nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
 
3
c given the set of values z(i,j) on the rectangular grid (x(i),y(j)),
 
4
c i=1,...,mx;j=1,...,my, subroutine regrid determines a smooth bivar-
 
5
c iate spline approximation s(x,y) of degrees kx and ky on the rect-
 
6
c angle xb <= x <= xe, yb <= y <= ye.
 
7
c if iopt = -1 regrid calculates the least-squares spline according
 
8
c to a given set of knots.
 
9
c if iopt >= 0 the total numbers nx and ny of these knots and their
 
10
c position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
 
11
c ally by the routine. the smoothness of s(x,y) is then achieved by
 
12
c minimalizing the discontinuity jumps in the derivatives of s(x,y)
 
13
c across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
 
14
c the amounth of smoothness is determined by the condition that f(p) =
 
15
c sum ((z(i,j)-s(x(i),y(j))))**2) be <= s, with s a given non-negative
 
16
c constant, called the smoothing factor.
 
17
c the fit is given in the b-spline representation (b-spline coefficients
 
18
c c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
 
19
c uated by means of subroutine bispev.
 
20
c
 
21
c calling sequence:
 
22
c     call regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
 
23
c    *  nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
 
24
c
 
25
c parameters:
 
26
c  iopt  : integer flag. on entry iopt must specify whether a least-
 
27
c          squares spline (iopt=-1) or a smoothing spline (iopt=0 or 1)
 
28
c          must be determined.
 
29
c          if iopt=0 the routine will start with an initial set of knots
 
30
c          tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=
 
31
c          1,...,ky+1. if iopt=1 the routine will continue with the set
 
32
c          of knots found at the last call of the routine.
 
33
c          attention: a call with iopt=1 must always be immediately pre-
 
34
c                     ceded by another call with iopt=1 or iopt=0 and
 
35
c                     s.ne.0.
 
36
c          unchanged on exit.
 
37
c  mx    : integer. on entry mx must specify the number of grid points
 
38
c          along the x-axis. mx > kx . unchanged on exit.
 
39
c  x     : real array of dimension at least (mx). before entry, x(i)
 
40
c          must be set to the x-co-ordinate of the i-th grid point
 
41
c          along the x-axis, for i=1,2,...,mx. these values must be
 
42
c          supplied in strictly ascending order. unchanged on exit.
 
43
c  my    : integer. on entry my must specify the number of grid points
 
44
c          along the y-axis. my > ky . unchanged on exit.
 
45
c  y     : real array of dimension at least (my). before entry, y(j)
 
46
c          must be set to the y-co-ordinate of the j-th grid point
 
47
c          along the y-axis, for j=1,2,...,my. these values must be
 
48
c          supplied in strictly ascending order. unchanged on exit.
 
49
c  z     : real array of dimension at least (mx*my).
 
50
c          before entry, z(my*(i-1)+j) must be set to the data value at
 
51
c          the grid point (x(i),y(j)) for i=1,...,mx and j=1,...,my.
 
52
c          unchanged on exit.
 
53
c  xb,xe : real values. on entry xb,xe,yb and ye must specify the bound-
 
54
c  yb,ye   aries of the rectangular approximation domain.
 
55
c          xb<=x(i)<=xe,i=1,...,mx; yb<=y(j)<=ye,j=1,...,my.
 
56
c          unchanged on exit.
 
57
c  kx,ky : integer values. on entry kx and ky must specify the degrees
 
58
c          of the spline. 1<=kx,ky<=5. it is recommended to use bicubic
 
59
c          (kx=ky=3) splines. unchanged on exit.
 
60
c  s     : real. on entry (in case iopt>=0) s must specify the smoothing
 
61
c          factor. s >=0. unchanged on exit.
 
62
c          for advice on the choice of s see further comments
 
63
c  nxest : integer. unchanged on exit.
 
64
c  nyest : integer. unchanged on exit.
 
65
c          on entry, nxest and nyest must specify an upper bound for the
 
66
c          number of knots required in the x- and y-directions respect.
 
67
c          these numbers will also determine the storage space needed by
 
68
c          the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1).
 
69
c          in most practical situation nxest = mx/2, nyest=my/2, will
 
70
c          be sufficient. always large enough are nxest=mx+kx+1, nyest=
 
71
c          my+ky+1, the number of knots needed for interpolation (s=0).
 
72
c          see also further comments.
 
73
c  nx    : integer.
 
74
c          unless ier=10 (in case iopt >=0), nx will contain the total
 
75
c          number of knots with respect to the x-variable, of the spline
 
76
c          approximation returned. if the computation mode iopt=1 is
 
77
c          used, the value of nx should be left unchanged between sub-
 
78
c          sequent calls.
 
79
c          in case iopt=-1, the value of nx should be specified on entry
 
80
c  tx    : real array of dimension nmax.
 
81
c          on succesful exit, this array will contain the knots of the
 
82
c          spline with respect to the x-variable, i.e. the position of
 
83
c          the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the
 
84
c          position of the additional knots tx(1)=...=tx(kx+1)=xb and
 
85
c          tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat.
 
86
c          if the computation mode iopt=1 is used, the values of tx(1),
 
87
c          ...,tx(nx) should be left unchanged between subsequent calls.
 
88
c          if the computation mode iopt=-1 is used, the values tx(kx+2),
 
89
c          ...tx(nx-kx-1) must be supplied by the user, before entry.
 
90
c          see also the restrictions (ier=10).
 
91
c  ny    : integer.
 
92
c          unless ier=10 (in case iopt >=0), ny will contain the total
 
93
c          number of knots with respect to the y-variable, of the spline
 
94
c          approximation returned. if the computation mode iopt=1 is
 
95
c          used, the value of ny should be left unchanged between sub-
 
96
c          sequent calls.
 
97
c          in case iopt=-1, the value of ny should be specified on entry
 
98
c  ty    : real array of dimension nmax.
 
99
c          on succesful exit, this array will contain the knots of the
 
100
c          spline with respect to the y-variable, i.e. the position of
 
101
c          the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the
 
102
c          position of the additional knots ty(1)=...=ty(ky+1)=yb and
 
103
c          ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat.
 
104
c          if the computation mode iopt=1 is used, the values of ty(1),
 
105
c          ...,ty(ny) should be left unchanged between subsequent calls.
 
106
c          if the computation mode iopt=-1 is used, the values ty(ky+2),
 
107
c          ...ty(ny-ky-1) must be supplied by the user, before entry.
 
108
c          see also the restrictions (ier=10).
 
109
c  c     : real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
 
110
c          on succesful exit, c contains the coefficients of the spline
 
111
c          approximation s(x,y)
 
112
c  fp    : real. unless ier=10, fp contains the sum of squared
 
113
c          residuals of the spline approximation returned.
 
114
c  wrk   : real array of dimension (lwrk). used as workspace.
 
115
c          if the computation mode iopt=1 is used the values of wrk(1),
 
116
c          ...,wrk(4) should be left unchanged between subsequent calls.
 
117
c  lwrk  : integer. on entry lwrk must specify the actual dimension of
 
118
c          the array wrk as declared in the calling (sub)program.
 
119
c          lwrk must not be too small.
 
120
c           lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
 
121
c            my*(ky+1) +u
 
122
c           where u is the larger of my and nxest.
 
123
c  iwrk  : integer array of dimension (kwrk). used as workspace.
 
124
c          if the computation mode iopt=1 is used the values of iwrk(1),
 
125
c          ...,iwrk(3) should be left unchanged between subsequent calls
 
126
c  kwrk  : integer. on entry kwrk must specify the actual dimension of
 
127
c          the array iwrk as declared in the calling (sub)program.
 
128
c          kwrk >= 3+mx+my+nxest+nyest.
 
129
c  ier   : integer. unless the routine detects an error, ier contains a
 
130
c          non-positive value on exit, i.e.
 
131
c   ier=0  : normal return. the spline returned has a residual sum of
 
132
c            squares fp such that abs(fp-s)/s <= tol with tol a relat-
 
133
c            ive tolerance set to 0.001 by the program.
 
134
c   ier=-1 : normal return. the spline returned is an interpolating
 
135
c            spline (fp=0).
 
136
c   ier=-2 : normal return. the spline returned is the least-squares
 
137
c            polynomial of degrees kx and ky. in this extreme case fp
 
138
c            gives the upper bound for the smoothing factor s.
 
139
c   ier=1  : error. the required storage space exceeds the available
 
140
c            storage space, as specified by the parameters nxest and
 
141
c            nyest.
 
142
c            probably causes : nxest or nyest too small. if these param-
 
143
c            eters are already large, it may also indicate that s is
 
144
c            too small
 
145
c            the approximation returned is the least-squares spline
 
146
c            according to the current set of knots. the parameter fp
 
147
c            gives the corresponding sum of squared residuals (fp>s).
 
148
c   ier=2  : error. a theoretically impossible result was found during
 
149
c            the iteration proces for finding a smoothing spline with
 
150
c            fp = s. probably causes : s too small.
 
151
c            there is an approximation returned but the corresponding
 
152
c            sum of squared residuals does not satisfy the condition
 
153
c            abs(fp-s)/s < tol.
 
154
c   ier=3  : error. the maximal number of iterations maxit (set to 20
 
155
c            by the program) allowed for finding a smoothing spline
 
156
c            with fp=s has been reached. probably causes : s too small
 
157
c            there is an approximation returned but the corresponding
 
158
c            sum of squared residuals does not satisfy the condition
 
159
c            abs(fp-s)/s < tol.
 
160
c   ier=10 : error. on entry, the input data are controlled on validity
 
161
c            the following restrictions must be satisfied.
 
162
c            -1<=iopt<=1, 1<=kx,ky<=5, mx>kx, my>ky, nxest>=2*kx+2,
 
163
c            nyest>=2*ky+2, kwrk>=3+mx+my+nxest+nyest,
 
164
c            lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
 
165
c             my*(ky+1) +max(my,nxest),
 
166
c            xb<=x(i-1)<x(i)<=xe,i=2,..,mx,yb<=y(j-1)<y(j)<=ye,j=2,..,my
 
167
c            if iopt=-1: 2*kx+2<=nx<=min(nxest,mx+kx+1)
 
168
c                        xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
 
169
c                        2*ky+2<=ny<=min(nyest,my+ky+1)
 
170
c                        yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
 
171
c                    the schoenberg-whitney conditions, i.e. there must
 
172
c                    be subset of grid co-ordinates xx(p) and yy(q) such
 
173
c                    that   tx(p) < xx(p) < tx(p+kx+1) ,p=1,...,nx-kx-1
 
174
c                           ty(q) < yy(q) < ty(q+ky+1) ,q=1,...,ny-ky-1
 
175
c            if iopt>=0: s>=0
 
176
c                        if s=0 : nxest>=mx+kx+1, nyest>=my+ky+1
 
177
c            if one of these conditions is found to be violated,control
 
178
c            is immediately repassed to the calling program. in that
 
179
c            case there is no approximation returned.
 
180
c
 
181
c further comments:
 
182
c   regrid does not allow individual weighting of the data-values.
 
183
c   so, if these were determined to widely different accuracies, then
 
184
c   perhaps the general data set routine surfit should rather be used
 
185
c   in spite of efficiency.
 
186
c   by means of the parameter s, the user can control the tradeoff
 
187
c   between closeness of fit and smoothness of fit of the approximation.
 
188
c   if s is too large, the spline will be too smooth and signal will be
 
189
c   lost ; if s is too small the spline will pick up too much noise. in
 
190
c   the extreme cases the program will return an interpolating spline if
 
191
c   s=0 and the least-squares polynomial (degrees kx,ky) if s is
 
192
c   very large. between these extremes, a properly chosen s will result
 
193
c   in a good compromise between closeness of fit and smoothness of fit.
 
194
c   to decide whether an approximation, corresponding to a certain s is
 
195
c   satisfactory the user is highly recommended to inspect the fits
 
196
c   graphically.
 
197
c   recommended values for s depend on the accuracy of the data values.
 
198
c   if the user has an idea of the statistical errors on the data, he
 
199
c   can also find a proper estimate for s. for, by assuming that, if he
 
200
c   specifies the right s, regrid will return a spline s(x,y) which
 
201
c   exactly reproduces the function underlying the data he can evaluate
 
202
c   the sum((z(i,j)-s(x(i),y(j)))**2) to find a good estimate for this s
 
203
c   for example, if he knows that the statistical errors on his z(i,j)-
 
204
c   values is not greater than 0.1, he may expect that a good s should
 
205
c   have a value not larger than mx*my*(0.1)**2.
 
206
c   if nothing is known about the statistical error in z(i,j), s must
 
207
c   be determined by trial and error, taking account of the comments
 
208
c   above. the best is then to start with a very large value of s (to
 
209
c   determine the least-squares polynomial and the corresponding upper
 
210
c   bound fp0 for s) and then to progressively decrease the value of s
 
211
c   ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
 
212
c   and more carefully as the approximation shows more detail) to
 
213
c   obtain closer fits.
 
214
c   to economize the search for a good s-value the program provides with
 
215
c   different modes of computation. at the first call of the routine, or
 
216
c   whenever he wants to restart with the initial set of knots the user
 
217
c   must set iopt=0.
 
218
c   if iopt=1 the program will continue with the set of knots found at
 
219
c   the last call of the routine. this will save a lot of computation
 
220
c   time if regrid is called repeatedly for different values of s.
 
221
c   the number of knots of the spline returned and their location will
 
222
c   depend on the value of s and on the complexity of the shape of the
 
223
c   function underlying the data. if the computation mode iopt=1
 
224
c   is used, the knots returned may also depend on the s-values at
 
225
c   previous calls (if these were smaller). therefore, if after a number
 
226
c   of trials with different s-values and iopt=1, the user can finally
 
227
c   accept a fit as satisfactory, it may be worthwhile for him to call
 
228
c   regrid once more with the selected value for s but now with iopt=0.
 
229
c   indeed, regrid may then return an approximation of the same quality
 
230
c   of fit but with fewer knots and therefore better if data reduction
 
231
c   is also an important objective for the user.
 
232
c   the number of knots may also depend on the upper bounds nxest and
 
233
c   nyest. indeed, if at a certain stage in regrid the number of knots
 
234
c   in one direction (say nx) has reached the value of its upper bound
 
235
c   (nxest), then from that moment on all subsequent knots are added
 
236
c   in the other (y) direction. this may indicate that the value of
 
237
c   nxest is too small. on the other hand, it gives the user the option
 
238
c   of limiting the number of knots the routine locates in any direction
 
239
c   for example, by setting nxest=2*kx+2 (the lowest allowable value for
 
240
c   nxest), the user can indicate that he wants an approximation which
 
241
c   is a simple polynomial of degree kx in the variable x.
 
242
c
 
243
c  other subroutines required:
 
244
c    fpback,fpbspl,fpregr,fpdisc,fpgivs,fpgrre,fprati,fprota,fpchec,
 
245
c    fpknot
 
246
c
 
247
c  references:
 
248
c   dierckx p. : a fast algorithm for smoothing data on a rectangular
 
249
c                grid while using spline functions, siam j.numer.anal.
 
250
c                19 (1982) 1286-1304.
 
251
c   dierckx p. : a fast algorithm for smoothing data on a rectangular
 
252
c                grid while using spline functions, report tw53, dept.
 
253
c                computer science,k.u.leuven, 1980.
 
254
c   dierckx p. : curve and surface fitting with splines, monographs on
 
255
c                numerical analysis, oxford university press, 1993.
 
256
c
 
257
c  author:
 
258
c    p.dierckx
 
259
c    dept. computer science, k.u. leuven
 
260
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
 
261
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
 
262
c
 
263
c  creation date : may 1979
 
264
c  latest update : march 1989
 
265
c
 
266
c  ..
 
267
c  ..scalar arguments..
 
268
      real*8 xb,xe,yb,ye,s,fp
 
269
      integer iopt,mx,my,kx,ky,nxest,nyest,nx,ny,lwrk,kwrk,ier
 
270
c  ..array arguments..
 
271
      real*8 x(mx),y(my),z(mx*my),tx(nxest),ty(nyest),
 
272
     * c((nxest-kx-1)*(nyest-ky-1)),wrk(lwrk)
 
273
      integer iwrk(kwrk)
 
274
c  ..local scalars..
 
275
      real*8 tol
 
276
      integer i,j,jwrk,kndx,kndy,knrx,knry,kwest,kx1,kx2,ky1,ky2,
 
277
     * lfpx,lfpy,lwest,lww,maxit,nc,nminx,nminy,mz
 
278
c  ..function references..
 
279
      integer max0
 
280
c  ..subroutine references..
 
281
c    fpregr,fpchec
 
282
c  ..
 
283
c  we set up the parameters tol and maxit.
 
284
      maxit = 20
 
285
      tol = 0.1e-02
 
286
c  before starting computations a data check is made. if the input data
 
287
c  are invalid, control is immediately repassed to the calling program.
 
288
      ier = 10
 
289
      if(kx.le.0 .or. kx.gt.5) go to 70
 
290
      kx1 = kx+1
 
291
      kx2 = kx1+1
 
292
      if(ky.le.0 .or. ky.gt.5) go to 70
 
293
      ky1 = ky+1
 
294
      ky2 = ky1+1
 
295
      if(iopt.lt.(-1) .or. iopt.gt.1) go to 70
 
296
      nminx = 2*kx1
 
297
      if(mx.lt.kx1 .or. nxest.lt.nminx) go to 70
 
298
      nminy = 2*ky1
 
299
      if(my.lt.ky1 .or. nyest.lt.nminy) go to 70
 
300
      mz = mx*my
 
301
      nc = (nxest-kx1)*(nyest-ky1)
 
302
      lwest = 4+nxest*(my+2*kx2+1)+nyest*(2*ky2+1)+mx*kx1+
 
303
     * my*ky1+max0(nxest,my)
 
304
      kwest = 3+mx+my+nxest+nyest
 
305
      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 70
 
306
      if(xb.gt.x(1) .or. xe.lt.x(mx)) go to 70
 
307
      do 10 i=2,mx
 
308
        if(x(i-1).ge.x(i)) go to 70
 
309
  10  continue
 
310
      if(yb.gt.y(1) .or. ye.lt.y(my)) go to 70
 
311
      do 20 i=2,my
 
312
        if(y(i-1).ge.y(i)) go to 70
 
313
  20  continue
 
314
      if(iopt.ge.0) go to 50
 
315
      if(nx.lt.nminx .or. nx.gt.nxest) go to 70
 
316
      j = nx
 
317
      do 30 i=1,kx1
 
318
        tx(i) = xb
 
319
        tx(j) = xe
 
320
        j = j-1
 
321
  30  continue
 
322
      call fpchec(x,mx,tx,nx,kx,ier)
 
323
      if(ier.ne.0) go to 70
 
324
      if(ny.lt.nminy .or. ny.gt.nyest) go to 70
 
325
      j = ny
 
326
      do 40 i=1,ky1
 
327
        ty(i) = yb
 
328
        ty(j) = ye
 
329
        j = j-1
 
330
  40  continue
 
331
      call fpchec(y,my,ty,ny,ky,ier)
 
332
      if (ier.eq.0) go to 60
 
333
      go to 70
 
334
  50  if(s.lt.0.) go to 70
 
335
      if(s.eq.0. .and. (nxest.lt.(mx+kx1) .or. nyest.lt.(my+ky1)) )
 
336
     * go to 70
 
337
      ier = 0
 
338
c  we partition the working space and determine the spline approximation
 
339
  60  lfpx = 5
 
340
      lfpy = lfpx+nxest
 
341
      lww = lfpy+nyest
 
342
      jwrk = lwrk-4-nxest-nyest
 
343
      knrx = 4
 
344
      knry = knrx+mx
 
345
      kndx = knry+my
 
346
      kndy = kndx+nxest
 
347
      call fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
 
348
     * tol,maxit,nc,nx,tx,ny,ty,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),
 
349
     * wrk(lfpx),wrk(lfpy),iwrk(1),iwrk(2),iwrk(3),iwrk(knrx),
 
350
     * iwrk(knry),iwrk(kndx),iwrk(kndy),wrk(lww),jwrk,ier)
 
351
  70  return
 
352
      end
 
353