2
fitpack --- curve and surface fitting with splines
4
fitpack is based on a collection of Fortran routines DIERCKX
5
by P. Dierckx (see http://www.netlib.org/dierckx/) transformed
6
to double routines by Pearu Peterson.
8
# Created by Pearu Peterson, June,August 2003
12
'InterpolatedUnivariateSpline',
13
'LSQUnivariateSpline',
16
'SmoothBivariateSpline']
19
from scipy_base import zeros, Float, concatenate, alltrue
20
from scipy_distutils.misc_util import PostponedException
23
except: dfitpack = PostponedException()
25
################ Univariate spline ####################
27
_curfit_messages = {1:"""
28
The required storage space exceeds the available storage space, as
29
specified by the parameter nest: nest too small. If nest is already
30
large (say nest > m/2), it may also indicate that s is too small.
31
The approximation returned is the weighted least-squares spline
32
according to the knots t[0],t[1],...,t[n-1]. (n=nest) the parameter fp
33
gives the corresponding weighted sum of squared residuals (fp>s).
36
A theoretically impossible result was found during the iteration
37
proces for finding a smoothing spline with fp = s: s too small.
38
There is an approximation returned but the corresponding weighted sum
39
of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
41
The maximal number of iterations maxit (set to 20 by the program)
42
allowed for finding a smoothing spline with fp=s has been reached: s
44
There is an approximation returned but the corresponding weighted sum
45
of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
47
Error on entry, no approximation returned. The following conditions
49
xb<=x[0]<x[1]<...<x[m-1]<=xe, w[i]>0, i=0..m-1
51
xb<t[k+1]<t[k+2]<...<t[n-k-2]<xe"""
54
class UnivariateSpline:
55
""" Univariate spline s(x) of degree k on the interval
56
[xb,xe] calculated from a given set of data points
60
def __init__(self, x, y, w=None, bbox = [None]*2, k=3, s=None):
63
x,y - 1-d sequences of data points (x must be
64
in strictly ascending order)
67
w - positive 1-d sequence of weights
68
bbox - 2-sequence specifying the boundary of
69
the approximation interval.
70
By default, bbox=[x[0],x[-1]]
71
k=3 - degree of the univariate spline.
72
s - positive smoothing factor defined for
74
sum((w[i]*(y[i]-s(x[i])))**2) <= s
75
Default s=len(w) which should be a good value
76
if 1/w[i] is an estimate of the standard
79
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
80
data = dfitpack.fpcurf0(x,y,k,w=w,
81
xb=bbox[0],xe=bbox[1],s=s)
83
# nest too small, setting to maximum bound
84
data = self._reset_nest(data)
88
def _reset_class(self):
90
n,t,c,k,ier = data[7],data[8],data[9],data[5],data[-1]
91
self._eval_args = t[:n],c[:n-k-1],k
93
# the spline returned has a residual sum of squares fp
94
# such that abs(fp-s)/s <= tol with tol a relative
95
# tolerance set to 0.001 by the program
98
# the spline returned is an interpolating spline
99
self.__class__ = InterpolatedUnivariateSpline
101
# the spline returned is the weighted least-squares
102
# polynomial of degree k. In this extreme case fp gives
103
# the upper bound fp0 for the smoothing factor s.
104
self.__class__ = LSQUnivariateSpline
108
self.__class__ = LSQUnivariateSpline
109
message = _curfit_messages.get(ier,'ier=%s' % (ier))
110
warnings.warn(message)
112
def _reset_nest(self, data, nest=None):
115
k,m = data[5],len(data[0])
116
nest = m+k+1 # this is the maximum bound for nest
118
assert n<=nest,"nest can only be increased"
119
t,c,fpint,nrdata = data[8].copy(),data[9].copy(),\
120
data[11].copy(),data[12].copy()
125
args = data[:8] + (t,c,n,fpint,nrdata,data[13])
126
data = dfitpack.fpcurf1(*args)
129
def set_smoothing_factor(self, s):
130
""" Continue spline computation with the given smoothing
131
factor s and with the knots found at the last call.
135
warnings.warn('smoothing factor unchanged for'
136
'LSQ spline with fixed knots')
138
args = data[:6] + (s,) + data[7:]
139
data = dfitpack.fpcurf1(*args)
141
# nest too small, setting to maximum bound
142
data = self._reset_nest(data)
146
def __call__(self, x, nu=None):
147
""" Evaluate spline (or its nu-th derivative) at positions x.
148
Note: x can be unordered but the evaluation is
149
more efficient if x is (partially) ordered.
152
return dfitpack.splev(*(self._eval_args+(x,)))
153
return dfitpack.splder(nu=nu,*(self._eval_args+(x,)))
156
""" Return the positions of (boundary and interior)
160
k,n = data[5],data[7]
161
return data[8][k:n-k]
163
def get_coeffs(self):
164
""" Return spline coefficients."""
166
k,n = data[5],data[7]
167
return data[9][:n-k-1]
169
def get_residual(self):
170
""" Return weighted sum of squared residuals of the spline
171
approximation: sum ((w[i]*(y[i]-s(x[i])))**2)
173
return self._data[10]
175
def integral(self, a, b):
176
""" Return definite integral of the spline between two
179
return dfitpack.splint(*(self._eval_args+(a,b)))
181
def derivatives(self, x):
182
""" Return all derivatives of the spline at the point x."""
183
d,ier = dfitpack.spalde(*(self._eval_args+(x,)))
188
""" Return the zeros of the spline.
190
Restriction: only cubic splines are supported by fitpack.
194
z,m,ier = dfitpack.sproot(*self._eval_args[:2])
197
raise NotImplementedError,'finding roots unsupported for non-cubic splines'
199
class InterpolatedUnivariateSpline(UnivariateSpline):
200
""" Interpolated univariate spline approximation."""
201
def __init__(self, x, y, w=None, bbox = [None]*2, k=3):
204
x,y - 1-d sequences of data points (x must be
205
in strictly ascending order)
208
w - positive 1-d sequence of weights
209
bbox - 2-sequence specifying the boundary of
210
the approximation interval.
211
By default, bbox=[x[0],x[-1]]
212
k=3 - degree of the univariate spline.
214
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
215
self._data = dfitpack.fpcurf0(x,y,k,w=w,
216
xb=bbox[0],xe=bbox[1],s=0)
219
class LSQUnivariateSpline(UnivariateSpline):
220
""" Weighted least-squares univariate spline approximation."""
222
def __init__(self, x, y, t, w=None, bbox = [None]*2, k=3):
225
x,y - 1-d sequences of data points (x must be
226
in strictly ascending order)
227
t - 1-d sequence of the positions of user-defined
228
interior knots of the spline (t must be in strictly
229
ascending order and bbox[0]<t[0]<...<t[-1]<bbox[-1])
232
w - positive 1-d sequence of weights
233
bbox - 2-sequence specifying the boundary of
234
the approximation interval.
235
By default, bbox=[x[0],x[-1]]
236
k=3 - degree of the univariate spline.
238
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
241
if xb is None: xb = x[0]
242
if xe is None: xe = x[-1]
243
t = concatenate(([xb]*(k+1),t,[xe]*(k+1)))
245
if not alltrue(t[k+1:n-k]-t[k:n-k-1] > 0):
247
'Interior knots t must satisfy Schoenberg-Whitney conditions'
248
data = dfitpack.fpcurfm1(x,y,k,t,w=w,xb=xb,xe=xe)
249
self._data = data[:-3] + (None,None,data[-1])
253
################ Bivariate spline ####################
255
_surfit_messages = {1:"""
256
The required storage space exceeds the available storage space: nxest
257
or nyest too small, or s too small.
258
The weighted least-squares spline corresponds to the current set of
261
A theoretically impossible result was found during the iteration
262
process for finding a smoothing spline with fp = s: s too small or
264
Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
266
the maximal number of iterations maxit (set to 20 by the program)
267
allowed for finding a smoothing spline with fp=s has been reached:
269
Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
271
No more knots can be added because the number of b-spline coefficients
272
(nx-kx-1)*(ny-ky-1) already exceeds the number of data points m:
273
either s or m too small.
274
The weighted least-squares spline corresponds to the current set of
277
No more knots can be added because the additional knot would (quasi)
278
coincide with an old one: s too small or too large a weight to an
279
inaccurate data point.
280
The weighted least-squares spline corresponds to the current set of
283
Error on entry, no approximation returned. The following conditions
285
xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
287
xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
288
yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""",
290
The coefficients of the spline returned have been computed as the
291
minimal norm least-squares solution of a (numerically) rank deficient
292
system (deficiency=%i). If deficiency is large, the results may be
293
inaccurate. Deficiency may strongly depend on the value of eps."""
297
class BivariateSpline:
298
""" Bivariate spline s(x,y) of degrees kx and ky on the rectangle
299
[xb,xe] x [yb, ye] calculated from a given set of data points
303
def get_residual(self):
304
""" Return weighted sum of squared residuals of the spline
305
approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2)
309
""" Return a tuple (tx,ty) where tx,ty contain knots positions
310
of the spline with respect to x-, y-variable, respectively.
311
The position of interior and additional knots are given as
312
t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.
315
def get_coeffs(self):
316
""" Return spline coefficients."""
318
def __call__(self,x,y,mth='array'):
319
""" Evaluate spline at positions x,y."""
321
tx,ty,c = self.tck[:3]
323
z,ier = dfitpack.bispev(tx,ty,c,kx,ky,x,y)
324
assert ier==0,'Invalid input: ier='+`ier`
326
raise NotImplementedError
328
class SmoothBivariateSpline(BivariateSpline):
329
""" Smooth bivariate spline approximation."""
331
def __init__(self, x, y, z, w=None,
332
bbox = [None]*4, kx=3, ky=3, s=None, eps=None):
335
x,y,z - 1-d sequences of data points (order is not
338
w - positive 1-d sequence of weights
339
bbox - 4-sequence specifying the boundary of
340
the rectangular approximation domain.
341
By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
342
kx,ky=3,3 - degrees of the bivariate spline.
343
s - positive smoothing factor defined for
344
estimation condition:
345
sum((w[i]*(z[i]-s(x[i],y[i])))**2) <= s
346
Default s=len(w) which should be a good value
347
if 1/w[i] is an estimate of the standard
349
eps - a threshold for determining the effective rank
350
of an over-determined linear system of
351
equations. 0 < eps < 1, default is 1e-16.
354
nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth(x,y,z,w,\
356
kx,ky,s=s,eps=eps,lwrk2=1)
357
if ier in [0,-1,-2]: # normal return
360
message = _surfit_messages.get(ier,'ier=%s' % (ier))
361
warnings.warn(message)
364
self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
367
class LSQBivariateSpline(BivariateSpline):
368
""" Weighted least-squares spline approximation.
370
def __init__(self, x, y, z, tx, ty, w=None,
372
kx=3, ky=3, eps=None):
375
x,y,z - 1-d sequences of data points (order is not
377
tx,ty - strictly ordered 1-d sequences of knots
380
w - positive 1-d sequence of weights
381
bbox - 4-sequence specifying the boundary of
382
the rectangular approximation domain.
383
By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
384
kx,ky=3,3 - degrees of the bivariate spline.
385
eps - a threshold for determining the effective rank
386
of an over-determined linear system of
387
equations. 0 < eps < 1, default is 1e-16.
391
tx1 = zeros((nx,),Float)
392
ty1 = zeros((ny,),Float)
393
tx1[kx+1:nx-kx-1] = tx
394
ty1[ky+1:ny-ky-1] = ty
397
tx1,ty1,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx1,ty1,w,\
401
tx1,ty1,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx1,ty1,w,\
404
if ier in [0,-1,-2]: # normal return
408
deficiency = (nx-kx-1)*(ny-ky-1)+ier
409
message = _surfit_messages.get(-3) % (deficiency)
411
message = _surfit_messages.get(ier,'ier=%s' % (ier))
412
warnings.warn(message)
added
removed
Lib/interpolate/fitpack2.py
