← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/sandbox/spline/fitpack/clocur.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/spline/fitpack/clocur.f
 
1
      subroutine clocur(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,
 
2
     * wrk,lwrk,iwrk,ier)
 
3
c  given the ordered set of m points x(i) in the idim-dimensional space
 
4
c  with x(1)=x(m), and given also a corresponding set of strictly in-
 
5
c  creasing values u(i) and the set of positive numbers w(i),i=1,2,...,m
 
6
c  subroutine clocur determines a smooth approximating closed spline
 
7
c  curve s(u), i.e.
 
8
c      x1 = s1(u)
 
9
c      x2 = s2(u)       u(1) <= u <= u(m)
 
10
c      .........
 
11
c      xidim = sidim(u)
 
12
c  with sj(u),j=1,2,...,idim periodic spline functions of degree k with
 
13
c  common knots t(j),j=1,2,...,n.
 
14
c  if ipar=1 the values u(i),i=1,2,...,m must be supplied by the user.
 
15
c  if ipar=0 these values are chosen automatically by clocur as
 
16
c      v(1) = 0
 
17
c      v(i) = v(i-1) + dist(x(i),x(i-1)) ,i=2,3,...,m
 
18
c      u(i) = v(i)/v(m) ,i=1,2,...,m
 
19
c  if iopt=-1 clocur calculates the weighted least-squares closed spline
 
20
c  curve according to a given set of knots.
 
21
c  if iopt>=0 the number of knots of the splines sj(u) and the position
 
22
c  t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
 
23
c  ness of s(u) is then achieved by minimalizing the discontinuity
 
24
c  jumps of the k-th derivative of s(u) at the knots t(j),j=k+2,k+3,...,
 
25
c  n-k-1. the amount of smoothness is determined by the condition that
 
26
c  f(p)=sum((w(i)*dist(x(i),s(u(i))))**2) be <= s, with s a given non-
 
27
c  negative constant, called the smoothing factor.
 
28
c  the fit s(u) is given in the b-spline representation and can be
 
29
c  evaluated by means of subroutine curev.
 
30
c
 
31
c  calling sequence:
 
32
c     call clocur(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,
 
33
c    * fp,wrk,lwrk,iwrk,ier)
 
34
c
 
35
c  parameters:
 
36
c   iopt  : integer flag. on entry iopt must specify whether a weighted
 
37
c           least-squares closed spline curve (iopt=-1) or a smoothing
 
38
c           closed spline curve (iopt=0 or 1) must be determined. if
 
39
c           iopt=0 the routine will start with an initial set of knots
 
40
c           t(i)=u(1)+(u(m)-u(1))*(i-k-1),i=1,2,...,2*k+2. if iopt=1 the
 
41
c           routine will continue with the knots found at the last call.
 
42
c           attention: a call with iopt=1 must always be immediately
 
43
c           preceded by another call with iopt=1 or iopt=0.
 
44
c           unchanged on exit.
 
45
c   ipar  : integer flag. on entry ipar must specify whether (ipar=1)
 
46
c           the user will supply the parameter values u(i),or whether
 
47
c           (ipar=0) these values are to be calculated by clocur.
 
48
c           unchanged on exit.
 
49
c   idim  : integer. on entry idim must specify the dimension of the
 
50
c           curve. 0 < idim < 11.
 
51
c           unchanged on exit.
 
52
c   m     : integer. on entry m must specify the number of data points.
 
53
c           m > 1. unchanged on exit.
 
54
c   u     : real array of dimension at least (m). in case ipar=1,before
 
55
c           entry, u(i) must be set to the i-th value of the parameter
 
56
c           variable u for i=1,2,...,m. these values must then be
 
57
c           supplied in strictly ascending order and will be unchanged
 
58
c           on exit. in case ipar=0, on exit,the array will contain the
 
59
c           values u(i) as determined by clocur.
 
60
c   mx    : integer. on entry mx must specify the actual dimension of
 
61
c           the array x as declared in the calling (sub)program. mx must
 
62
c           not be too small (see x). unchanged on exit.
 
63
c   x     : real array of dimension at least idim*m.
 
64
c           before entry, x(idim*(i-1)+j) must contain the j-th coord-
 
65
c           inate of the i-th data point for i=1,2,...,m and j=1,2,...,
 
66
c           idim. since first and last data point must coincide it
 
67
c           means that x(j)=x(idim*(m-1)+j),j=1,2,...,idim.
 
68
c           unchanged on exit.
 
69
c   w     : real array of dimension at least (m). before entry, w(i)
 
70
c           must be set to the i-th value in the set of weights. the
 
71
c           w(i) must be strictly positive. w(m) is not used.
 
72
c           unchanged on exit. see also further comments.
 
73
c   k     : integer. on entry k must specify the degree of the splines.
 
74
c           1<=k<=5. it is recommended to use cubic splines (k=3).
 
75
c           the user is strongly dissuaded from choosing k even,together
 
76
c           with a small s-value. unchanged on exit.
 
77
c   s     : real.on entry (in case iopt>=0) s must specify the smoothing
 
78
c           factor. s >=0. unchanged on exit.
 
79
c           for advice on the choice of s see further comments.
 
80
c   nest  : integer. on entry nest must contain an over-estimate of the
 
81
c           total number of knots of the splines returned, to indicate
 
82
c           the storage space available to the routine. nest >=2*k+2.
 
83
c           in most practical situation nest=m/2 will be sufficient.
 
84
c           always large enough is nest=m+2*k, the number of knots
 
85
c           needed for interpolation (s=0). unchanged on exit.
 
86
c   n     : integer.
 
87
c           unless ier = 10 (in case iopt >=0), n will contain the
 
88
c           total number of knots of the smoothing spline curve returned
 
89
c           if the computation mode iopt=1 is used this value of n
 
90
c           should be left unchanged between subsequent calls.
 
91
c           in case iopt=-1, the value of n must be specified on entry.
 
92
c   t     : real array of dimension at least (nest).
 
93
c           on succesful exit, this array will contain the knots of the
 
94
c           spline curve,i.e. the position of the interior knots t(k+2),
 
95
c           t(k+3),..,t(n-k-1) as well as the position of the additional
 
96
c           t(1),t(2),..,t(k+1)=u(1) and u(m)=t(n-k),...,t(n) needed for
 
97
c           the b-spline representation.
 
98
c           if the computation mode iopt=1 is used, the values of t(1),
 
99
c           t(2),...,t(n) should be left unchanged between subsequent
 
100
c           calls. if the computation mode iopt=-1 is used, the values
 
101
c           t(k+2),...,t(n-k-1) must be supplied by the user, before
 
102
c           entry. see also the restrictions (ier=10).
 
103
c   nc    : integer. on entry nc must specify the actual dimension of
 
104
c           the array c as declared in the calling (sub)program. nc
 
105
c           must not be too small (see c). unchanged on exit.
 
106
c   c     : real array of dimension at least (nest*idim).
 
107
c           on succesful exit, this array will contain the coefficients
 
108
c           in the b-spline representation of the spline curve s(u),i.e.
 
109
c           the b-spline coefficients of the spline sj(u) will be given
 
110
c           in c(n*(j-1)+i),i=1,2,...,n-k-1 for j=1,2,...,idim.
 
111
c   fp    : real. unless ier = 10, fp contains the weighted sum of
 
112
c           squared residuals of the spline curve returned.
 
113
c   wrk   : real array of dimension at least m*(k+1)+nest*(7+idim+5*k).
 
114
c           used as working space. if the computation mode iopt=1 is
 
115
c           used, the values wrk(1),...,wrk(n) should be left unchanged
 
116
c           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. lwrk
 
119
c           must not be too small (see wrk). unchanged on exit.
 
120
c   iwrk  : integer array of dimension at least (nest).
 
121
c           used as working space. if the computation mode iopt=1 is
 
122
c           used,the values iwrk(1),...,iwrk(n) should be left unchanged
 
123
c           between subsequent calls.
 
124
c   ier   : integer. unless the routine detects an error, ier contains a
 
125
c           non-positive value on exit, i.e.
 
126
c    ier=0  : normal return. the close curve returned has a residual
 
127
c             sum of squares fp such that abs(fp-s)/s <= tol with tol a
 
128
c             relative tolerance set to 0.001 by the program.
 
129
c    ier=-1 : normal return. the curve returned is an interpolating
 
130
c             spline curve (fp=0).
 
131
c    ier=-2 : normal return. the curve returned is the weighted least-
 
132
c             squares point,i.e. each spline sj(u) is a constant. in
 
133
c             this extreme case fp gives the upper bound fp0 for the
 
134
c             smoothing factor s.
 
135
c    ier=1  : error. the required storage space exceeds the available
 
136
c             storage space, as specified by the parameter nest.
 
137
c             probably causes : nest too small. if nest is already
 
138
c             large (say nest > m/2), it may also indicate that s is
 
139
c             too small
 
140
c             the approximation returned is the least-squares closed
 
141
c             curve according to the knots t(1),t(2),...,t(n). (n=nest)
 
142
c             the parameter fp gives the corresponding weighted sum of
 
143
c             squared residuals (fp>s).
 
144
c    ier=2  : error. a theoretically impossible result was found during
 
145
c             the iteration proces for finding a smoothing curve with
 
146
c             fp = s. probably causes : s too small.
 
147
c             there is an approximation returned but the corresponding
 
148
c             weighted sum of squared residuals does not satisfy the
 
149
c             condition abs(fp-s)/s < tol.
 
150
c    ier=3  : error. the maximal number of iterations maxit (set to 20
 
151
c             by the program) allowed for finding a smoothing curve
 
152
c             with fp=s has been reached. probably causes : s too small
 
153
c             there is an approximation returned but the corresponding
 
154
c             weighted sum of squared residuals does not satisfy the
 
155
c             condition abs(fp-s)/s < tol.
 
156
c    ier=10 : error. on entry, the input data are controlled on validity
 
157
c             the following restrictions must be satisfied.
 
158
c             -1<=iopt<=1, 1<=k<=5, m>1, nest>2*k+2, w(i)>0,i=1,2,...,m
 
159
c             0<=ipar<=1, 0<idim<=10, lwrk>=(k+1)*m+nest*(7+idim+5*k),
 
160
c             nc>=nest*idim, x(j)=x(idim*(m-1)+j), j=1,2,...,idim
 
161
c             if ipar=0: sum j=1,idim (x(i*idim+j)-x((i-1)*idim+j))**2>0
 
162
c                        i=1,2,...,m-1.
 
163
c             if ipar=1: u(1)<u(2)<...<u(m)
 
164
c             if iopt=-1: 2*k+2<=n<=min(nest,m+2*k)
 
165
c                         u(1)<t(k+2)<t(k+3)<...<t(n-k-1)<u(m)
 
166
c                            (u(1)=0 and u(m)=1 in case ipar=0)
 
167
c                       the schoenberg-whitney conditions, i.e. there
 
168
c                       must be a subset of data points uu(j) with
 
169
c                       uu(j) = u(i) or u(i)+(u(m)-u(1)) such that
 
170
c                         t(j) < uu(j) < t(j+k+1), j=k+1,...,n-k-1
 
171
c             if iopt>=0: s>=0
 
172
c                         if s=0 : nest >= m+2*k
 
173
c             if one of these conditions is found to be violated,control
 
174
c             is immediately repassed to the calling program. in that
 
175
c             case there is no approximation returned.
 
176
c
 
177
c  further comments:
 
178
c   by means of the parameter s, the user can control the tradeoff
 
179
c   between closeness of fit and smoothness of fit of the approximation.
 
180
c   if s is too large, the curve will be too smooth and signal will be
 
181
c   lost ; if s is too small the curve will pick up too much noise. in
 
182
c   the extreme cases the program will return an interpolating curve if
 
183
c   s=0 and the weighted least-squares point if s is very large.
 
184
c   between these extremes, a properly chosen s will result in a good
 
185
c   compromise between closeness of fit and smoothness of fit.
 
186
c   to decide whether an approximation, corresponding to a certain s is
 
187
c   satisfactory the user is highly recommended to inspect the fits
 
188
c   graphically.
 
189
c   recommended values for s depend on the weights w(i). if these are
 
190
c   taken as 1/d(i) with d(i) an estimate of the standard deviation of
 
191
c   x(i), a good s-value should be found in the range (m-sqrt(2*m),m+
 
192
c   sqrt(2*m)). if nothing is known about the statistical error in x(i)
 
193
c   each w(i) can be set equal to one and s determined by trial and
 
194
c   error, taking account of the comments above. the best is then to
 
195
c   start with a very large value of s ( to determine the weighted
 
196
c   least-squares point and the upper bound fp0 for s) and then to
 
197
c   progressively decrease the value of s ( say by a factor 10 in the
 
198
c   beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
 
199
c   approximating curve shows more detail) to obtain closer fits.
 
200
c   to economize the search for a good s-value the program provides with
 
201
c   different modes of computation. at the first call of the routine, or
 
202
c   whenever he wants to restart with the initial set of knots the user
 
203
c   must set iopt=0.
 
204
c   if iopt=1 the program will continue with the set of knots found at
 
205
c   the last call of the routine. this will save a lot of computation
 
206
c   time if clocur is called repeatedly for different values of s.
 
207
c   the number of knots of the spline returned and their location will
 
208
c   depend on the value of s and on the complexity of the shape of the
 
209
c   curve underlying the data. but, if the computation mode iopt=1 is
 
210
c   used, the knots returned may also depend on the s-values at previous
 
211
c   calls (if these were smaller). therefore, if after a number of
 
212
c   trials with different s-values and iopt=1, the user can finally
 
213
c   accept a fit as satisfactory, it may be worthwhile for him to call
 
214
c   clocur once more with the selected value for s but now with iopt=0.
 
215
c   indeed, clocur may then return an approximation of the same quality
 
216
c   of fit but with fewer knots and therefore better if data reduction
 
217
c   is also an important objective for the user.
 
218
c
 
219
c   the form of the approximating curve can strongly be affected  by
 
220
c   the choice of the parameter values u(i). if there is no physical
 
221
c   reason for choosing a particular parameter u, often good results
 
222
c   will be obtained with the choice of clocur(in case ipar=0), i.e.
 
223
c        v(1)=0, v(i)=v(i-1)+q(i), i=2,...,m, u(i)=v(i)/v(m), i=1,..,m
 
224
c   where
 
225
c        q(i)= sqrt(sum j=1,idim (xj(i)-xj(i-1))**2 )
 
226
c   other possibilities for q(i) are
 
227
c        q(i)= sum j=1,idim (xj(i)-xj(i-1))**2
 
228
c        q(i)= sum j=1,idim abs(xj(i)-xj(i-1))
 
229
c        q(i)= max j=1,idim abs(xj(i)-xj(i-1))
 
230
c        q(i)= 1
 
231
c
 
232
c
 
233
c  other subroutines required:
 
234
c    fpbacp,fpbspl,fpchep,fpclos,fpdisc,fpgivs,fpknot,fprati,fprota
 
235
c
 
236
c  references:
 
237
c   dierckx p. : algorithms for smoothing data with periodic and
 
238
c                parametric splines, computer graphics and image
 
239
c                processing 20 (1982) 171-184.
 
240
c   dierckx p. : algorithms for smoothing data with periodic and param-
 
241
c                etric splines, report tw55, dept. computer science,
 
242
c                k.u.leuven, 1981.
 
243
c   dierckx p. : curve and surface fitting with splines, monographs on
 
244
c                numerical analysis, oxford university press, 1993.
 
245
c
 
246
c  author:
 
247
c    p.dierckx
 
248
c    dept. computer science, k.u. leuven
 
249
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
 
250
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
 
251
c
 
252
c  creation date : may 1979
 
253
c  latest update : march 1987
 
254
c
 
255
c  ..
 
256
c  ..scalar arguments..
 
257
      real*8 s,fp
 
258
      integer iopt,ipar,idim,m,mx,k,nest,n,nc,lwrk,ier
 
259
c  ..array arguments..
 
260
      real*8 u(m),x(mx),w(m),t(nest),c(nc),wrk(lwrk)
 
261
      integer iwrk(nest)
 
262
c  ..local scalars..
 
263
      real*8 per,tol,dist
 
264
      integer i,ia1,ia2,ib,ifp,ig1,ig2,iq,iz,i1,i2,j1,j2,k1,k2,lwest,
 
265
     * maxit,m1,nmin,ncc,j
 
266
c  ..function references..
 
267
      real*8 sqrt
 
268
c  we set up the parameters tol and maxit
 
269
      maxit = 20
 
270
      tol = 0.1e-02
 
271
c  before starting computations a data check is made. if the input data
 
272
c  are invalid, control is immediately repassed to the calling program.
 
273
      ier = 10
 
274
      if(iopt.lt.(-1) .or. iopt.gt.1) go to 90
 
275
      if(ipar.lt.0 .or. ipar.gt.1) go to 90
 
276
      if(idim.le.0 .or. idim.gt.10) go to 90
 
277
      if(k.le.0 .or. k.gt.5) go to 90
 
278
      k1 = k+1
 
279
      k2 = k1+1
 
280
      nmin = 2*k1
 
281
      if(m.lt.2 .or. nest.lt.nmin) go to 90
 
282
      ncc = nest*idim
 
283
      if(mx.lt.m*idim .or. nc.lt.ncc) go to 90
 
284
      lwest = m*k1+nest*(7+idim+5*k)
 
285
      if(lwrk.lt.lwest) go to 90
 
286
      i1 = idim
 
287
      i2 = m*idim
 
288
      do 5 j=1,idim
 
289
         if(x(i1).ne.x(i2)) go to 90
 
290
         i1 = i1-1
 
291
         i2 = i2-1
 
292
   5  continue
 
293
      if(ipar.ne.0 .or. iopt.gt.0) go to 40
 
294
      i1 = 0
 
295
      i2 = idim
 
296
      u(1) = 0.
 
297
      do 20 i=2,m
 
298
         dist = 0.
 
299
         do 10 j1=1,idim
 
300
            i1 = i1+1
 
301
            i2 = i2+1
 
302
            dist = dist+(x(i2)-x(i1))**2
 
303
  10     continue
 
304
         u(i) = u(i-1)+sqrt(dist)
 
305
  20  continue
 
306
      if(u(m).le.0.) go to 90
 
307
      do 30 i=2,m
 
308
         u(i) = u(i)/u(m)
 
309
  30  continue
 
310
      u(m) = 0.1e+01
 
311
  40  if(w(1).le.0.) go to 90
 
312
      m1 = m-1
 
313
      do 50 i=1,m1
 
314
         if(u(i).ge.u(i+1) .or. w(i).le.0.) go to 90
 
315
  50  continue
 
316
      if(iopt.ge.0) go to 70
 
317
      if(n.le.nmin .or. n.gt.nest) go to 90
 
318
      per = u(m)-u(1)
 
319
      j1 = k1
 
320
      t(j1) = u(1)
 
321
      i1 = n-k
 
322
      t(i1) = u(m)
 
323
      j2 = j1
 
324
      i2 = i1
 
325
      do 60 i=1,k
 
326
         i1 = i1+1
 
327
         i2 = i2-1
 
328
         j1 = j1+1
 
329
         j2 = j2-1
 
330
         t(j2) = t(i2)-per
 
331
         t(i1) = t(j1)+per
 
332
  60  continue
 
333
      call fpchep(u,m,t,n,k,ier)
 
334
      if (ier.eq.0) go to 80
 
335
      go to 90
 
336
  70  if(s.lt.0.) go to 90
 
337
      if(s.eq.0. .and. nest.lt.(m+2*k)) go to 90
 
338
      ier = 0
 
339
c we partition the working space and determine the spline approximation.
 
340
  80  ifp = 1
 
341
      iz = ifp+nest
 
342
      ia1 = iz+ncc
 
343
      ia2 = ia1+nest*k1
 
344
      ib = ia2+nest*k
 
345
      ig1 = ib+nest*k2
 
346
      ig2 = ig1+nest*k2
 
347
      iq = ig2+nest*k1
 
348
      call fpclos(iopt,idim,m,u,mx,x,w,k,s,nest,tol,maxit,k1,k2,n,t,
 
349
     * ncc,c,fp,wrk(ifp),wrk(iz),wrk(ia1),wrk(ia2),wrk(ib),wrk(ig1),
 
350
     * wrk(ig2),wrk(iq),iwrk,ier)
 
351
  90  return
 
352
      end