← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/interpolate/fitpack/parsur.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/parsur.f
 
1
      subroutine parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest,
 
2
     * nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
 
3
c  given the set of ordered points f(i,j) in the idim-dimensional space,
 
4
c  corresponding to grid values (u(i),v(j)) ,i=1,...,mu ; j=1,...,mv,
 
5
c  parsur determines a smooth approximating spline surface s(u,v) , i.e.
 
6
c    f1 = s1(u,v)
 
7
c      ...                u(1) <= u <= u(mu) ; v(1) <= v <= v(mv)
 
8
c    fidim = sidim(u,v)
 
9
c  with sl(u,v), l=1,2,...,idim bicubic spline functions with common
 
10
c  knots tu(i),i=1,...,nu in the u-variable and tv(j),j=1,...,nv in the
 
11
c  v-variable.
 
12
c  in addition, these splines will be periodic in the variable u if
 
13
c  ipar(1) = 1 and periodic in the variable v if ipar(2) = 1.
 
14
c  if iopt=-1, parsur determines the least-squares bicubic spline
 
15
c  surface according to a given set of knots.
 
16
c  if iopt>=0, the number of knots of s(u,v) and their position
 
17
c  is chosen automatically by the routine. the smoothness of s(u,v) is
 
18
c  achieved by minimalizing the discontinuity jumps of the derivatives
 
19
c  of the splines at the knots. the amount of smoothness of s(u,v) is
 
20
c  determined by the condition that
 
21
c  fp=sumi=1,mu(sumj=1,mv(dist(f(i,j)-s(u(i),v(j)))**2))<=s,
 
22
c  with s a given non-negative constant.
 
23
c  the fit s(u,v) is given in its b-spline representation and can be
 
24
c  evaluated by means of routine surev.
 
25
c
 
26
c calling sequence:
 
27
c     call parsur(iopt,ipar,idim,mu,u,mv,v,f,s,nuest,nvest,nu,tu,
 
28
c    *  nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
 
29
c
 
30
c parameters:
 
31
c  iopt  : integer flag. unchanged on exit.
 
32
c          on entry iopt must specify whether a least-squares surface
 
33
c          (iopt=-1) or a smoothing surface (iopt=0 or 1)must be
 
34
c          determined.
 
35
c          if iopt=0 the routine will start with the initial set of
 
36
c          knots needed for determining the least-squares polynomial
 
37
c          surface.
 
38
c          if iopt=1 the routine will continue with the set of knots
 
39
c          found at the last call of the routine.
 
40
c          attention: a call with iopt=1 must always be immediately
 
41
c          preceded by another call with iopt = 1 or iopt = 0.
 
42
c  ipar  : integer array of dimension 2. unchanged on exit.
 
43
c          on entry ipar(1) must specify whether (ipar(1)=1) or not
 
44
c          (ipar(1)=0) the splines must be periodic in the variable u.
 
45
c          on entry ipar(2) must specify whether (ipar(2)=1) or not
 
46
c          (ipar(2)=0) the splines must be periodic in the variable v.
 
47
c  idim  : integer. on entry idim must specify the dimension of the
 
48
c          surface. 1 <= idim <= 3. unchanged on exit.
 
49
c  mu    : integer. on entry mu must specify the number of grid points
 
50
c          along the u-axis. unchanged on exit.
 
51
c          mu >= mumin where mumin=4-2*ipar(1)
 
52
c  u     : real array of dimension at least (mu). before entry, u(i)
 
53
c          must be set to the u-co-ordinate of the i-th grid point
 
54
c          along the u-axis, for i=1,2,...,mu. these values must be
 
55
c          supplied in strictly ascending order. unchanged on exit.
 
56
c  mv    : integer. on entry mv must specify the number of grid points
 
57
c          along the v-axis. unchanged on exit.
 
58
c          mv >= mvmin where mvmin=4-2*ipar(2)
 
59
c  v     : real array of dimension at least (mv). before entry, v(j)
 
60
c          must be set to the v-co-ordinate of the j-th grid point
 
61
c          along the v-axis, for j=1,2,...,mv. these values must be
 
62
c          supplied in strictly ascending order. unchanged on exit.
 
63
c  f     : real array of dimension at least (mu*mv*idim).
 
64
c          before entry, f(mu*mv*(l-1)+mv*(i-1)+j) must be set to the
 
65
c          l-th co-ordinate of the data point corresponding to the
 
66
c          the grid point (u(i),v(j)) for l=1,...,idim ,i=1,...,mu
 
67
c          and j=1,...,mv. unchanged on exit.
 
68
c          if ipar(1)=1 it is expected that f(mu*mv*(l-1)+mv*(mu-1)+j)
 
69
c          = f(mu*mv*(l-1)+j), l=1,...,idim ; j=1,...,mv
 
70
c          if ipar(2)=1 it is expected that f(mu*mv*(l-1)+mv*(i-1)+mv)
 
71
c          = f(mu*mv*(l-1)+mv*(i-1)+1), l=1,...,idim ; i=1,...,mu
 
72
c  s     : real. on entry (if iopt>=0) s must specify the smoothing
 
73
c          factor. s >=0. unchanged on exit.
 
74
c          for advice on the choice of s see further comments
 
75
c  nuest : integer. unchanged on exit.
 
76
c  nvest : integer. unchanged on exit.
 
77
c          on entry, nuest and nvest must specify an upper bound for the
 
78
c          number of knots required in the u- and v-directions respect.
 
79
c          these numbers will also determine the storage space needed by
 
80
c          the routine. nuest >= 8, nvest >= 8.
 
81
c          in most practical situation nuest = mu/2, nvest=mv/2, will
 
82
c          be sufficient. always large enough are nuest=mu+4+2*ipar(1),
 
83
c          nvest = mv+4+2*ipar(2), the number of knots needed for
 
84
c          interpolation (s=0). see also further comments.
 
85
c  nu    : integer.
 
86
c          unless ier=10 (in case iopt>=0), nu will contain the total
 
87
c          number of knots with respect to the u-variable, of the spline
 
88
c          surface returned. if the computation mode iopt=1 is used,
 
89
c          the value of nu should be left unchanged between subsequent
 
90
c          calls. in case iopt=-1, the value of nu should be specified
 
91
c          on entry.
 
92
c  tu    : real array of dimension at least (nuest).
 
93
c          on succesful exit, this array will contain the knots of the
 
94
c          splines with respect to the u-variable, i.e. the position of
 
95
c          the interior knots tu(5),...,tu(nu-4) as well as the position
 
96
c          of the additional knots tu(1),...,tu(4) and tu(nu-3),...,
 
97
c          tu(nu) needed for the b-spline representation.
 
98
c          if the computation mode iopt=1 is used,the values of tu(1)
 
99
c          ...,tu(nu) should be left unchanged between subsequent calls.
 
100
c          if the computation mode iopt=-1 is used, the values tu(5),
 
101
c          ...tu(nu-4) must be supplied by the user, before entry.
 
102
c          see also the restrictions (ier=10).
 
103
c  nv    : integer.
 
104
c          unless ier=10 (in case iopt>=0), nv will contain the total
 
105
c          number of knots with respect to the v-variable, of the spline
 
106
c          surface returned. if the computation mode iopt=1 is used,
 
107
c          the value of nv should be left unchanged between subsequent
 
108
c          calls. in case iopt=-1, the value of nv should be specified
 
109
c          on entry.
 
110
c  tv    : real array of dimension at least (nvest).
 
111
c          on succesful exit, this array will contain the knots of the
 
112
c          splines with respect to the v-variable, i.e. the position of
 
113
c          the interior knots tv(5),...,tv(nv-4) as well as the position
 
114
c          of the additional knots tv(1),...,tv(4) and tv(nv-3),...,
 
115
c          tv(nv) needed for the b-spline representation.
 
116
c          if the computation mode iopt=1 is used,the values of tv(1)
 
117
c          ...,tv(nv) should be left unchanged between subsequent calls.
 
118
c          if the computation mode iopt=-1 is used, the values tv(5),
 
119
c          ...tv(nv-4) must be supplied by the user, before entry.
 
120
c          see also the restrictions (ier=10).
 
121
c  c     : real array of dimension at least (nuest-4)*(nvest-4)*idim.
 
122
c          on succesful exit, c contains the coefficients of the spline
 
123
c          approximation s(u,v)
 
124
c  fp    : real. unless ier=10, fp contains the sum of squared
 
125
c          residuals of the spline surface returned.
 
126
c  wrk   : real array of dimension (lwrk). used as workspace.
 
127
c          if the computation mode iopt=1 is used the values of
 
128
c          wrk(1),...,wrk(4) should be left unchanged between subsequent
 
129
c          calls.
 
130
c  lwrk  : integer. on entry lwrk must specify the actual dimension of
 
131
c          the array wrk as declared in the calling (sub)program.
 
132
c          lwrk must not be too small.
 
133
c           lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+
 
134
c           4*(mu+mv)+q*idim where q is the larger of mv and nuest.
 
135
c  iwrk  : integer array of dimension (kwrk). used as workspace.
 
136
c          if the computation mode iopt=1 is used the values of
 
137
c          iwrk(1),.,iwrk(3) should be left unchanged between subsequent
 
138
c          calls.
 
139
c  kwrk  : integer. on entry kwrk must specify the actual dimension of
 
140
c          the array iwrk as declared in the calling (sub)program.
 
141
c          kwrk >= 3+mu+mv+nuest+nvest.
 
142
c  ier   : integer. unless the routine detects an error, ier contains a
 
143
c          non-positive value on exit, i.e.
 
144
c   ier=0  : normal return. the surface returned has a residual sum of
 
145
c            squares fp such that abs(fp-s)/s <= tol with tol a relat-
 
146
c            ive tolerance set to 0.001 by the program.
 
147
c   ier=-1 : normal return. the spline surface returned is an
 
148
c            interpolating surface (fp=0).
 
149
c   ier=-2 : normal return. the surface returned is the least-squares
 
150
c            polynomial surface. in this extreme case fp gives the
 
151
c            upper bound for the smoothing factor s.
 
152
c   ier=1  : error. the required storage space exceeds the available
 
153
c            storage space, as specified by the parameters nuest and
 
154
c            nvest.
 
155
c            probably causes : nuest or nvest too small. if these param-
 
156
c            eters are already large, it may also indicate that s is
 
157
c            too small
 
158
c            the approximation returned is the least-squares surface
 
159
c            according to the current set of knots. the parameter fp
 
160
c            gives the corresponding sum of squared residuals (fp>s).
 
161
c   ier=2  : error. a theoretically impossible result was found during
 
162
c            the iteration proces for finding a smoothing surface with
 
163
c            fp = s. probably causes : s too small.
 
164
c            there is an approximation returned but the corresponding
 
165
c            sum of squared residuals does not satisfy the condition
 
166
c            abs(fp-s)/s < tol.
 
167
c   ier=3  : error. the maximal number of iterations maxit (set to 20
 
168
c            by the program) allowed for finding a smoothing surface
 
169
c            with fp=s has been reached. probably causes : s too small
 
170
c            there is an approximation returned but the corresponding
 
171
c            sum of squared residuals does not satisfy the condition
 
172
c            abs(fp-s)/s < tol.
 
173
c   ier=10 : error. on entry, the input data are controlled on validity
 
174
c            the following restrictions must be satisfied.
 
175
c            -1<=iopt<=1, 0<=ipar(1)<=1, 0<=ipar(2)<=1, 1 <=idim<=3
 
176
c            mu >= 4-2*ipar(1),mv >= 4-2*ipar(2), nuest >=8, nvest >= 8,
 
177
c            kwrk>=3+mu+mv+nuest+nvest,
 
178
c            lwrk >= 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))
 
179
c             +4*(mu+mv)+max(nuest,mv)*idim
 
180
c            u(i-1)<u(i),i=2,..,mu, v(i-1)<v(i),i=2,...,mv
 
181
c            if iopt=-1: 8<=nu<=min(nuest,mu+4+2*ipar(1))
 
182
c                        u(1)<tu(5)<tu(6)<...<tu(nu-4)<u(mu)
 
183
c                        8<=nv<=min(nvest,mv+4+2*ipar(2))
 
184
c                        v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(mv)
 
185
c                    the schoenberg-whitney conditions, i.e. there must
 
186
c                    be subset of grid co-ordinates uu(p) and vv(q) such
 
187
c                    that   tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
 
188
c                           tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
 
189
c                     (see fpchec or fpchep)
 
190
c            if iopt>=0: s>=0
 
191
c                       if s=0: nuest>=mu+4+2*ipar(1)
 
192
c                               nvest>=mv+4+2*ipar(2)
 
193
c            if one of these conditions is found to be violated,control
 
194
c            is immediately repassed to the calling program. in that
 
195
c            case there is no approximation returned.
 
196
c
 
197
c further comments:
 
198
c   by means of the parameter s, the user can control the tradeoff
 
199
c   between closeness of fit and smoothness of fit of the approximation.
 
200
c   if s is too large, the surface will be too smooth and signal will be
 
201
c   lost ; if s is too small the surface will pick up too much noise. in
 
202
c   the extreme cases the program will return an interpolating surface
 
203
c   if s=0 and the constrained least-squares polynomial surface if s is
 
204
c   very large. between these extremes, a properly chosen s will result
 
205
c   in a good compromise between closeness of fit and smoothness of fit.
 
206
c   to decide whether an approximation, corresponding to a certain s is
 
207
c   satisfactory the user is highly recommended to inspect the fits
 
208
c   graphically.
 
209
c   recommended values for s depend on the accuracy of the data values.
 
210
c   if the user has an idea of the statistical errors on the data, he
 
211
c   can also find a proper estimate for s. for, by assuming that, if he
 
212
c   specifies the right s, parsur will return a surface s(u,v) which
 
213
c   exactly reproduces the surface underlying the data he can evaluate
 
214
c   the sum(dist(f(i,j)-s(u(i),v(j)))**2) to find a good estimate for s.
 
215
c   for example, if he knows that the statistical errors on his f(i,j)-
 
216
c   values is not greater than 0.1, he may expect that a good s should
 
217
c   have a value not larger than mu*mv*(0.1)**2.
 
218
c   if nothing is known about the statistical error in f(i,j), s must
 
219
c   be determined by trial and error, taking account of the comments
 
220
c   above. the best is then to start with a very large value of s (to
 
221
c   determine the le-sq polynomial surface and the corresponding upper
 
222
c   bound fp0 for s) and then to progressively decrease the value of s
 
223
c   ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
 
224
c   and more carefully as the approximation shows more detail) to
 
225
c   obtain closer fits.
 
226
c   to economize the search for a good s-value the program provides with
 
227
c   different modes of computation. at the first call of the routine, or
 
228
c   whenever he wants to restart with the initial set of knots the user
 
229
c   must set iopt=0.
 
230
c   if iopt = 1 the program will continue with the knots found at
 
231
c   the last call of the routine. this will save a lot of computation
 
232
c   time if parsur is called repeatedly for different values of s.
 
233
c   the number of knots of the surface returned and their location will
 
234
c   depend on the value of s and on the complexity of the shape of the
 
235
c   surface underlying the data. if the computation mode iopt = 1
 
236
c   is used, the knots returned may also depend on the s-values at
 
237
c   previous calls (if these were smaller). therefore, if after a number
 
238
c   of trials with different s-values and iopt=1,the user can finally
 
239
c   accept a fit as satisfactory, it may be worthwhile for him to call
 
240
c   parsur once more with the chosen value for s but now with iopt=0.
 
241
c   indeed, parsur may then return an approximation of the same quality
 
242
c   of fit but with fewer knots and therefore better if data reduction
 
243
c   is also an important objective for the user.
 
244
c   the number of knots may also depend on the upper bounds nuest and
 
245
c   nvest. indeed, if at a certain stage in parsur the number of knots
 
246
c   in one direction (say nu) has reached the value of its upper bound
 
247
c   (nuest), then from that moment on all subsequent knots are added
 
248
c   in the other (v) direction. this may indicate that the value of
 
249
c   nuest is too small. on the other hand, it gives the user the option
 
250
c   of limiting the number of knots the routine locates in any direction
 
251
c   for example, by setting nuest=8 (the lowest allowable value for
 
252
c   nuest), the user can indicate that he wants an approximation with
 
253
c   splines which are simple cubic polynomials in the variable u.
 
254
c
 
255
c  other subroutines required:
 
256
c    fppasu,fpchec,fpchep,fpknot,fprati,fpgrpa,fptrnp,fpback,
 
257
c    fpbacp,fpbspl,fptrpe,fpdisc,fpgivs,fprota
 
258
c
 
259
c  author:
 
260
c    p.dierckx
 
261
c    dept. computer science, k.u. leuven
 
262
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
 
263
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
 
264
c
 
265
c  latest update : march 1989
 
266
c
 
267
c  ..
 
268
c  ..scalar arguments..
 
269
      real*8 s,fp
 
270
      integer iopt,idim,mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier
 
271
c  ..array arguments..
 
272
      real*8 u(mu),v(mv),f(mu*mv*idim),tu(nuest),tv(nvest),
 
273
     * c((nuest-4)*(nvest-4)*idim),wrk(lwrk)
 
274
      integer ipar(2),iwrk(kwrk)
 
275
c  ..local scalars..
 
276
      real*8 tol,ub,ue,vb,ve,peru,perv
 
277
      integer i,j,jwrk,kndu,kndv,knru,knrv,kwest,l1,l2,l3,l4,
 
278
     * lfpu,lfpv,lwest,lww,maxit,nc,mf,mumin,mvmin
 
279
c  ..function references..
 
280
      integer max0
 
281
c  ..subroutine references..
 
282
c    fppasu,fpchec,fpchep
 
283
c  ..
 
284
c  we set up the parameters tol and maxit.
 
285
      maxit = 20
 
286
      tol = 0.1e-02
 
287
c  before starting computations a data check is made. if the input data
 
288
c  are invalid, control is immediately repassed to the calling program.
 
289
      ier = 10
 
290
      if(iopt.lt.(-1) .or. iopt.gt.1) go to 200
 
291
      if(ipar(1).lt.0 .or. ipar(1).gt.1) go to 200
 
292
      if(ipar(2).lt.0 .or. ipar(2).gt.1) go to 200
 
293
      if(idim.le.0 .or. idim.gt.3) go to 200
 
294
      mumin = 4-2*ipar(1)
 
295
      if(mu.lt.mumin .or. nuest.lt.8) go to 200
 
296
      mvmin = 4-2*ipar(2)
 
297
      if(mv.lt.mvmin .or. nvest.lt.8) go to 200
 
298
      mf = mu*mv
 
299
      nc = (nuest-4)*(nvest-4)
 
300
      lwest = 4+nuest*(mv*idim+11+4*ipar(1))+nvest*(11+4*ipar(2))+
 
301
     * 4*(mu+mv)+max0(nuest,mv)*idim
 
302
      kwest = 3+mu+mv+nuest+nvest
 
303
      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200
 
304
      do 10 i=2,mu
 
305
        if(u(i-1).ge.u(i)) go to 200
 
306
  10  continue
 
307
      do 20 i=2,mv
 
308
        if(v(i-1).ge.v(i)) go to 200
 
309
  20  continue
 
310
      if(iopt.ge.0) go to 100
 
311
      if(nu.lt.8 .or. nu.gt.nuest) go to 200
 
312
      ub = u(1)
 
313
      ue = u(mu)
 
314
      if (ipar(1).ne.0) go to 40
 
315
      j = nu
 
316
      do 30 i=1,4
 
317
        tu(i) = ub
 
318
        tu(j) = ue
 
319
        j = j-1
 
320
  30  continue
 
321
      call fpchec(u,mu,tu,nu,3,ier)
 
322
      if(ier.ne.0) go to 200
 
323
      go to 60
 
324
  40  l1 = 4
 
325
      l2 = l1
 
326
      l3 = nu-3
 
327
      l4 = l3
 
328
      peru = ue-ub
 
329
      tu(l2) = ub
 
330
      tu(l3) = ue
 
331
      do 50 j=1,3
 
332
        l1 = l1+1
 
333
        l2 = l2-1
 
334
        l3 = l3+1
 
335
        l4 = l4-1
 
336
        tu(l2) = tu(l4)-peru
 
337
        tu(l3) = tu(l1)+peru
 
338
  50  continue
 
339
      call fpchep(u,mu,tu,nu,3,ier)
 
340
      if(ier.ne.0) go to 200
 
341
  60  if(nv.lt.8 .or. nv.gt.nvest) go to 200
 
342
      vb = v(1)
 
343
      ve = v(mv)
 
344
      if (ipar(2).ne.0) go to 80
 
345
      j = nv
 
346
      do 70 i=1,4
 
347
        tv(i) = vb
 
348
        tv(j) = ve
 
349
        j = j-1
 
350
  70  continue
 
351
      call fpchec(v,mv,tv,nv,3,ier)
 
352
      if(ier.ne.0) go to 200
 
353
      go to 150
 
354
  80  l1 = 4
 
355
      l2 = l1
 
356
      l3 = nv-3
 
357
      l4 = l3
 
358
      perv = ve-vb
 
359
      tv(l2) = vb
 
360
      tv(l3) = ve
 
361
      do 90 j=1,3
 
362
        l1 = l1+1
 
363
        l2 = l2-1
 
364
        l3 = l3+1
 
365
        l4 = l4-1
 
366
        tv(l2) = tv(l4)-perv
 
367
        tv(l3) = tv(l1)+perv
 
368
  90  continue
 
369
      call fpchep(v,mv,tv,nv,3,ier)
 
370
      if (ier.eq.0) go to 150
 
371
      go to 200
 
372
 100  if(s.lt.0.) go to 200
 
373
      if(s.eq.0. .and. (nuest.lt.(mu+4+2*ipar(1)) .or.
 
374
     * nvest.lt.(mv+4+2*ipar(2))) )go to 200
 
375
      ier = 0
 
376
c  we partition the working space and determine the spline approximation
 
377
 150  lfpu = 5
 
378
      lfpv = lfpu+nuest
 
379
      lww = lfpv+nvest
 
380
      jwrk = lwrk-4-nuest-nvest
 
381
      knru = 4
 
382
      knrv = knru+mu
 
383
      kndu = knrv+mv
 
384
      kndv = kndu+nuest
 
385
      call fppasu(iopt,ipar,idim,u,mu,v,mv,f,mf,s,nuest,nvest,
 
386
     * tol,maxit,nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),
 
387
     * wrk(lfpu),wrk(lfpv),iwrk(1),iwrk(2),iwrk(3),iwrk(knru),
 
388
     * iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),jwrk,ier)
 
389
 200  return
 
390
      end
 
391