1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
|
subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,
* wrk,lwrk,iwrk,ier)
c given the set of data points (x(i),y(i)) and the set of positive
c numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
c approximation of degree k on the interval xb <= x <= xe.
c if iopt=-1 curfit calculates the weighted least-squares spline
c according to a given set of knots.
c if iopt>=0 the number of knots of the spline s(x) and the position
c t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
c ness of s(x) is then achieved by minimalizing the discontinuity
c jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
c n-k-1. the amount of smoothness is determined by the condition that
c f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
c negative constant, called the smoothing factor.
c the fit s(x) is given in the b-spline representation (b-spline coef-
c ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
c subroutine splev.
c
c calling sequence:
c call curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,
c * lwrk,iwrk,ier)
c
c parameters:
c iopt : integer flag. on entry iopt must specify whether a weighted
c least-squares spline (iopt=-1) or a smoothing spline (iopt=
c 0 or 1) must be determined. if iopt=0 the routine will start
c with an initial set of knots t(i)=xb, t(i+k+1)=xe, i=1,2,...
c k+1. if iopt=1 the routine will continue with the knots
c found at the last call of the routine.
c attention: a call with iopt=1 must always be immediately
c preceded by another call with iopt=1 or iopt=0.
c unchanged on exit.
c m : integer. on entry m must specify the number of data points.
c m > k. unchanged on exit.
c x : real array of dimension at least (m). before entry, x(i)
c must be set to the i-th value of the independent variable x,
c for i=1,2,...,m. these values must be supplied in strictly
c ascending order. unchanged on exit.
c y : real array of dimension at least (m). before entry, y(i)
c must be set to the i-th value of the dependent variable y,
c for i=1,2,...,m. unchanged on exit.
c w : real array of dimension at least (m). before entry, w(i)
c must be set to the i-th value in the set of weights. the
c w(i) must be strictly positive. unchanged on exit.
c see also further comments.
c xb,xe : real values. on entry xb and xe must specify the boundaries
c of the approximation interval. xb<=x(1), xe>=x(m).
c unchanged on exit.
c k : integer. on entry k must specify the degree of the spline.
c 1<=k<=5. it is recommended to use cubic splines (k=3).
c the user is strongly dissuaded from choosing k even,together
c with a small s-value. unchanged on exit.
c s : real.on entry (in case iopt>=0) s must specify the smoothing
c factor. s >=0. unchanged on exit.
c for advice on the choice of s see further comments.
c nest : integer. on entry nest must contain an over-estimate of the
c total number of knots of the spline returned, to indicate
c the storage space available to the routine. nest >=2*k+2.
c in most practical situation nest=m/2 will be sufficient.
c always large enough is nest=m+k+1, the number of knots
c needed for interpolation (s=0). unchanged on exit.
c n : integer.
c unless ier =10 (in case iopt >=0), n will contain the
c total number of knots of the spline approximation returned.
c if the computation mode iopt=1 is used this value of n
c should be left unchanged between subsequent calls.
c in case iopt=-1, the value of n must be specified on entry.
c t : real array of dimension at least (nest).
c on succesful exit, this array will contain the knots of the
c spline,i.e. the position of the interior knots t(k+2),t(k+3)
c ...,t(n-k-1) as well as the position of the additional knots
c t(1)=t(2)=...=t(k+1)=xb and t(n-k)=...=t(n)=xe needed for
c the b-spline representation.
c if the computation mode iopt=1 is used, the values of t(1),
c t(2),...,t(n) should be left unchanged between subsequent
c calls. if the computation mode iopt=-1 is used, the values
c t(k+2),...,t(n-k-1) must be supplied by the user, before
c entry. see also the restrictions (ier=10).
c c : real array of dimension at least (nest).
c on succesful exit, this array will contain the coefficients
c c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
c fp : real. unless ier=10, fp contains the weighted sum of
c squared residuals of the spline approximation returned.
c wrk : real array of dimension at least (m*(k+1)+nest*(7+3*k)).
c used as working space. if the computation mode iopt=1 is
c used, the values wrk(1),...,wrk(n) should be left unchanged
c between subsequent calls.
c lwrk : integer. on entry,lwrk must specify the actual dimension of
c the array wrk as declared in the calling (sub)program.lwrk
c must not be too small (see wrk). unchanged on exit.
c iwrk : integer array of dimension at least (nest).
c used as working space. if the computation mode iopt=1 is
c used,the values iwrk(1),...,iwrk(n) should be left unchanged
c between subsequent calls.
c ier : integer. unless the routine detects an error, ier contains a
c non-positive value on exit, i.e.
c ier=0 : normal return. the spline returned has a residual sum of
c squares fp such that abs(fp-s)/s <= tol with tol a relat-
c ive tolerance set to 0.001 by the program.
c ier=-1 : normal return. the spline returned is an interpolating
c spline (fp=0).
c ier=-2 : normal return. the spline returned is the weighted least-
c squares polynomial of degree k. in this extreme case fp
c gives the upper bound fp0 for the smoothing factor s.
c ier=1 : error. the required storage space exceeds the available
c storage space, as specified by the parameter nest.
c probably causes : nest too small. if nest is already
c large (say nest > m/2), it may also indicate that s is
c too small
c the approximation returned is the weighted least-squares
c spline according to the knots t(1),t(2),...,t(n). (n=nest)
c the parameter fp gives the corresponding weighted sum of
c squared residuals (fp>s).
c ier=2 : error. a theoretically impossible result was found during
c the iteration proces for finding a smoothing spline with
c fp = s. probably causes : s too small.
c there is an approximation returned but the corresponding
c weighted sum of squared residuals does not satisfy the
c condition abs(fp-s)/s < tol.
c ier=3 : error. the maximal number of iterations maxit (set to 20
c by the program) allowed for finding a smoothing spline
c with fp=s has been reached. probably causes : s too small
c there is an approximation returned but the corresponding
c weighted sum of squared residuals does not satisfy the
c condition abs(fp-s)/s < tol.
c ier=10 : error. on entry, the input data are controlled on validity
c the following restrictions must be satisfied.
c -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
c xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
c if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
c xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
c the schoenberg-whitney conditions, i.e. there
c must be a subset of data points xx(j) such that
c t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
c if iopt>=0: s>=0
c if s=0 : nest >= m+k+1
c if one of these conditions is found to be violated,control
c is immediately repassed to the calling program. in that
c case there is no approximation returned.
c
c further comments:
c by means of the parameter s, the user can control the tradeoff
c between closeness of fit and smoothness of fit of the approximation.
c if s is too large, the spline will be too smooth and signal will be
c lost ; if s is too small the spline will pick up too much noise. in
c the extreme cases the program will return an interpolating spline if
c s=0 and the weighted least-squares polynomial of degree k if s is
c very large. between these extremes, a properly chosen s will result
c in a good compromise between closeness of fit and smoothness of fit.
c to decide whether an approximation, corresponding to a certain s is
c satisfactory the user is highly recommended to inspect the fits
c graphically.
c recommended values for s depend on the weights w(i). if these are
c taken as 1/d(i) with d(i) an estimate of the standard deviation of
c y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
c sqrt(2*m)). if nothing is known about the statistical error in y(i)
c each w(i) can be set equal to one and s determined by trial and
c error, taking account of the comments above. the best is then to
c start with a very large value of s ( to determine the least-squares
c polynomial and the corresponding upper bound fp0 for s) and then to
c progressively decrease the value of s ( say by a factor 10 in the
c beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
c approximation shows more detail) to obtain closer fits.
c to economize the search for a good s-value the program provides with
c different modes of computation. at the first call of the routine, or
c whenever he wants to restart with the initial set of knots the user
c must set iopt=0.
c if iopt=1 the program will continue with the set of knots found at
c the last call of the routine. this will save a lot of computation
c time if curfit is called repeatedly for different values of s.
c the number of knots of the spline returned and their location will
c depend on the value of s and on the complexity of the shape of the
c function underlying the data. but, if the computation mode iopt=1
c is used, the knots returned may also depend on the s-values at
c previous calls (if these were smaller). therefore, if after a number
c of trials with different s-values and iopt=1, the user can finally
c accept a fit as satisfactory, it may be worthwhile for him to call
c curfit once more with the selected value for s but now with iopt=0.
c indeed, curfit may then return an approximation of the same quality
c of fit but with fewer knots and therefore better if data reduction
c is also an important objective for the user.
c
c other subroutines required:
c fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota
c
c references:
c dierckx p. : an algorithm for smoothing, differentiation and integ-
c ration of experimental data using spline functions,
c j.comp.appl.maths 1 (1975) 165-184.
c dierckx p. : a fast algorithm for smoothing data on a rectangular
c grid while using spline functions, siam j.numer.anal.
c 19 (1982) 1286-1304.
c dierckx p. : an improved algorithm for curve fitting with spline
c functions, report tw54, dept. computer science,k.u.
c leuven, 1981.
c dierckx p. : curve and surface fitting with splines, monographs on
c numerical analysis, oxford university press, 1993.
c
c author:
c p.dierckx
c dept. computer science, k.u. leuven
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c creation date : may 1979
c latest update : march 1987
c
c ..
c ..scalar arguments..
real*8 xb,xe,s,fp
integer iopt,m,k,nest,n,lwrk,ier
c ..array arguments..
real*8 x(m),y(m),w(m),t(nest),c(nest),wrk(lwrk)
integer iwrk(nest)
c ..local scalars..
real*8 tol
integer i,ia,ib,ifp,ig,iq,iz,j,k1,k2,lwest,maxit,nmin
c ..
c we set up the parameters tol and maxit
maxit = 20
tol = 0.1d-02
c before starting computations a data check is made. if the input data
c are invalid, control is immediately repassed to the calling program.
ier = 10
if(k.le.0 .or. k.gt.5) go to 50
k1 = k+1
k2 = k1+1
if(iopt.lt.(-1) .or. iopt.gt.1) go to 50
nmin = 2*k1
if(m.lt.k1 .or. nest.lt.nmin) go to 50
lwest = m*k1+nest*(7+3*k)
if(lwrk.lt.lwest) go to 50
if(xb.gt.x(1) .or. xe.lt.x(m)) go to 50
do 10 i=2,m
if(x(i-1).gt.x(i)) go to 50
10 continue
if(iopt.ge.0) go to 30
if(n.lt.nmin .or. n.gt.nest) go to 50
j = n
do 20 i=1,k1
t(i) = xb
t(j) = xe
j = j-1
20 continue
call fpchec(x,m,t,n,k,ier)
if (ier.eq.0) go to 40
go to 50
30 if(s.lt.0.) go to 50
if(s.eq.0. .and. nest.lt.(m+k1)) go to 50
c we partition the working space and determine the spline approximation.
40 ifp = 1
iz = ifp+nest
ia = iz+nest
ib = ia+nest*k1
ig = ib+nest*k2
iq = ig+nest*k2
call fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,
* wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
50 return
end
|
