3
fitpack (dierckx in netlib) --- A Python-C wrapper to FITPACK (by P. Dierckx).
4
FITPACK is a collection of FORTRAN programs for curve and surface
5
fitting with splines and tensor product splines.
8
http://www.cs.kuleuven.ac.be/cwis/research/nalag/research/topics/fitpack.html
10
http://www.netlib.org/dierckx/index.html
12
Copyright 2002 Pearu Peterson all rights reserved,
13
Pearu Peterson <pearu@cens.ioc.ee>
14
Permission to use, modify, and distribute this software is given under the
15
terms of the SciPy (BSD style) license. See LICENSE.txt that came with
16
this distribution for specifics.
18
NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
21
Modified by John Travers January 2007
23
Running test programs:
24
$ python fitpack.py 1 3 # run test programs 1, and 3
25
$ python fitpack.py # run all available test programs
27
TODO: Make interfaces to the following fitpack functions:
28
For univariate splines: cocosp, concon, fourco
29
For bivariate splines: profil, parsur, surev
32
__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
33
'bisplrep', 'bisplev', 'insert']
34
__version__ = "$Revision: 2902 $"[10:-1]
36
from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
37
dot, sin, cos, pi, arange, empty, int32, where, resize
38
myasarray = atleast_1d
40
# f2py-generated interface to fitpack
42
# new class based interface to fitpack
46
The spline has a residual sum of squares fp such that abs(fp-s)/
49
The spline is an interpolating spline (fp=0)""",None],
51
The spline is weighted least-squares polynomial of degree k.
52
fp gives the upper bound fp0 for the smoothing factor s""",None],
54
The required storage space exceeds the available storage space.
55
Probable causes: data (x,y) size is too small or smoothing parameter
56
s is too small (fp>s).""",ValueError],
58
A theoretically impossible results when finding a smoothin spline
59
with fp = s. Probably causes: s too small.
60
(abs(fp-s)/s>0.001)""",ValueError],
62
The maximal number of iterations (20) allowed for finding smoothing
63
spline with fp=s has been reached. Probably causes: s too small.
64
(abs(fp-s)/s>0.001)""",ValueError],
66
Error on input data""",ValueError],
68
An error occured""",TypeError]}
71
The spline has a residual sum of squares fp such that
72
abs(fp-s)/s<=0.001""",None],
74
The spline is an interpolating spline (fp=0)""",None],
76
The spline is weighted least-squares polynomial of degree kx and ky.
77
fp gives the upper bound fp0 for the smoothing factor s""",None],
79
Warning. The coefficients of the spline have been computed as the minimal
80
norm least-squares solution of a rank deficient system.""",None],
82
The required storage space exceeds the available storage space.
83
Probably causes: nxest or nyest too small or s is too small.
84
(fp>s)""",ValueError],
86
A theoretically impossible results when finding a smoothin spline
87
with fp = s. Probably causes: s too small or badly chosen eps.
88
(abs(fp-s)/s>0.001)""",ValueError],
90
The maximal number of iterations (20) allowed for finding smoothing
91
spline with fp=s has been reached. Probably causes: s too small.
92
(abs(fp-s)/s>0.001)""",ValueError],
94
No more knots can be added because the number of B-spline coefficients
95
already exceeds the number of data points m. Probably causes: either
96
s or m too small. (fp>s)""",ValueError],
98
No more knots can be added because the additional knot would coincide
99
with an old one. Probably cause: s too small or too large a weight
100
to an inaccurate data point. (fp>s)""",ValueError],
102
Error on input data""",ValueError],
104
rwrk2 too small, i.e. there is not enough workspace for computing
105
the minimal least-squares solution of a rank deficient system of linear
106
equations.""",ValueError],
108
An error occured""",TypeError]}
112
def splprep(x,w=None,u=None,ub=None,ue=None,k=3,task=0,s=None,t=None,
113
full_output=0,nest=None,per=0,quiet=1):
114
"""Find the B-spline representation of an N-dimensional curve.
117
Given a list of N rank-1 arrays, x, which represent a curve in
118
N-dimensional space parametrized by u, find a smooth approximating
120
Uses the FORTRAN routine parcur from FITPACK
123
x -- A list of sample vector arrays representing the curve.
124
u -- An array of parameter values. If not given, these values are
125
calculated automatically as (M = len(x[0])):
127
v[i] = v[i-1] + distance(x[i],x[i-1])
129
ub, ue -- The end-points of the parameters interval. Defaults to
131
k -- Degree of the spline. Cubic splines are recommended. Even values
132
of k should be avoided especially with a small s-value.
134
task -- If task==0 find t and c for a given smoothing factor, s.
135
If task==1 find t and c for another value of the smoothing
136
factor, s. There must have been a previous call with task=0
137
or task=1 for the same set of data.
138
If task=-1 find the weighted least square spline for a given
140
s -- A smoothing condition. The amount of smoothness is determined by
141
satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
142
g(x) is the smoothed interpolation of (x,y). The user can use s to
143
control the tradeoff between closeness and smoothness of fit. Larger
144
s means more smoothing while smaller values of s indicate less
145
smoothing. Recommended values of s depend on the weights, w. If the
146
weights represent the inverse of the standard-deviation of y, then a
147
good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
148
where m is the number of datapoints in x, y, and w.
149
t -- The knots needed for task=-1.
150
full_output -- If non-zero, then return optional outputs.
151
nest -- An over-estimate of the total number of knots of the spline to
152
help in determining the storage space. By default nest=m/2.
153
Always large enough is nest=m+k+1.
154
per -- If non-zero, data points are considered periodic with period
155
x[m-1] - x[0] and a smooth periodic spline approximation is
157
Values of y[m-1] and w[m-1] are not used.
158
quiet -- Non-zero to suppress messages.
160
Outputs: (tck, u, {fp, ier, msg})
161
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
162
coefficients, and the degree of the spline.
163
u -- An array of the values of the parameter.
165
fp -- The weighted sum of squared residuals of the spline approximation.
166
ier -- An integer flag about splrep success. Success is indicated
167
if ier<=0. If ier in [1,2,3] an error occurred but was not raised.
168
Otherwise an error is raised.
169
msg -- A message corresponding to the integer flag, ier.
172
SEE splev for evaluation of the spline and its derivatives.
175
splrep, splev, sproot, spalde, splint - evaluation, roots, integral
176
bisplrep, bisplev - bivariate splines
177
UnivariateSpline, BivariateSpline - an alternative wrapping
178
of the FITPACK functions
183
for i in range(idim):
184
if x[i][0]!=x[i][-1]:
185
if quiet<2:print 'Warning: Setting x[%d][%d]=x[%d][0]'%(i,m,i)
187
if not 0<idim<11: raise TypeError,'0<idim<11 must hold'
188
if w is None: w=ones(m,float)
192
if ub is None: ub=u[0]
193
if ue is None: ue=u[-1]
194
else: u=zeros(m,float); ub=ue=0.0
195
if not (1<=k<=5): raise TypeError, '1<=k=%d<=5 must hold'%(k)
196
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
197
if (not len(w)==m) or (ipar==1 and (not len(u)==m)):
198
raise TypeError,'Mismatch of input dimensions'
199
if s is None: s=m-sqrt(2*m)
201
if t is None: raise TypeError, 'Knots must be given for task=-1'
205
raise TypeError, 'There must be at least 2*k+2 knots for task=-1'
206
if m<=k: raise TypeError, 'm>k must hold'
207
if nest is None: nest=m+2*k
208
if (task>=0 and s==0) or (nest<0):
214
u,n,t,c,fp,wrk,iwrk,ier=dfitpack.clocur_smth0(ipar,idim,u,x,w,nest,
217
u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth0(ipar,idim,u,x,
218
w,ub,ue,nest,k=k,s=s)
223
wrk=_parcur_cache['wrk']
224
iwrk=_parcur_cache['iwrk']
227
ub=_parcur_cache['ub']
228
ue=_parcur_cache['ue']
230
raise ValueError, 'task=1 can only be called after task=0'
232
u,n,t,c,fp,wrk,iwrk,ier=dfitpack.clocur_smth1(ipar,idim,u,x,
233
w,n,t,wrk,iwrk,k=k,s=s)
235
u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth1(ipar,idim,u,x,
236
w,ub,ue,n,t,wrk,iwrk,k=k,s=s)
239
u,n,t,c,fp,ier=dfitpack.clocur_lsq(ipar,idim,u,x,w,n,t,k=k)
241
u,ub,ue,n,t,c,fp,ier=dfitpack.parcur_lsq(ipar,idim,u,x,w,ub,ue,
247
_parcur_cache['wrk']=wrk
248
_parcur_cache['iwrk']=iwrk
250
_parcur_cache['ub']=ub
251
_parcur_cache['ue']=ue
255
tcku = [t[:n],list(c),k],u
256
if ier<=0 and not quiet:
257
print _iermess[ier][0]
258
print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
259
if ier>0 and not full_output:
261
print "Warning: "+_iermess[ier][0]
264
raise _iermess[ier][1],_iermess[ier][0]
266
raise _iermess['unknown'][1],_iermess['unknown'][0]
269
return tcku,fp,ier,_iermess[ier][0]
271
return tcku,fp,ier,_iermess['unknown'][0]
278
def splrep(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
279
full_output=0,per=0,quiet=1):
280
"""Find the B-spline representation of 1-D curve.
283
Given the set of data points (x[i], y[i]) determine a smooth spline
284
approximation of degree k on the interval xb <= x <= xe. The coefficients,
285
c, and the knot points, t, are returned. Uses the FORTRAN routine
289
x, y -- The data points defining a curve y = f(x).
290
w -- Strictly positive rank-1 array of weights the same length as x and y.
291
The weights are used in computing the weighted least-squares spline
292
fit. If the errors in the y values have standard-deviation given by the
293
vector d, then w should be 1/d. Default is ones(len(x)).
294
xb, xe -- The interval to fit. If None, these default to x[0] and x[-1]
296
k -- The order of the spline fit. It is recommended to use cubic splines.
297
Even order splines should be avoided especially with small s values.
299
task -- If task==0 find t and c for a given smoothing factor, s.
300
If task==1 find t and c for another value of the
301
smoothing factor, s. There must have been a previous
302
call with task=0 or task=1 for the same set of data
303
(t will be stored an used internally)
304
If task=-1 find the weighted least square spline for
305
a given set of knots, t. These should be interior knots
306
as knots on the ends will be added automatically.
307
s -- A smoothing condition. The amount of smoothness is determined by
308
satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
309
g(x) is the smoothed interpolation of (x,y). The user can use s to
310
control the tradeoff between closeness and smoothness of fit. Larger
311
s means more smoothing while smaller values of s indicate less
312
smoothing. Recommended values of s depend on the weights, w. If the
313
weights represent the inverse of the standard-deviation of y, then a
314
good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
315
where m is the number of datapoints in x, y, and w.
316
default : s = m-sqrt(2*m) if weights are supplied.
317
s = 0.0 (interpolating) if no weights are supplied.
318
t -- The knots needed for task=-1. If given then task is automatically
320
full_output -- If non-zero, then return optional outputs.
321
per -- If non-zero, data points are considered periodic with period
322
x[m-1] - x[0] and a smooth periodic spline approximation is returned.
323
Values of y[m-1] and w[m-1] are not used.
324
quiet -- Non-zero to suppress messages.
326
Outputs: (tck, {fp, ier, msg})
327
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
328
coefficients, and the degree of the spline.
330
fp -- The weighted sum of squared residuals of the spline approximation.
331
ier -- An integer flag about splrep success. Success is indicated if
332
ier<=0. If ier in [1,2,3] an error occurred but was not raised.
333
Otherwise an error is raised.
334
msg -- A message corresponding to the integer flag, ier.
337
See splev for evaluation of the spline and its derivatives.
340
x = linspace(0, 10, 10)
343
x2 = linspace(0, 10, 200)
345
plot(x, y, 'o', x2, y2)
348
splprep, splev, sproot, spalde, splint - evaluation, roots, integral
349
bisplrep, bisplev - bivariate splines
350
UnivariateSpline, BivariateSpline - an alternative wrapping
351
of the FITPACK functions
353
x,y=map(myasarray,[x,y])
357
if s is None: s = 0.0
360
if s is None: s = m-sqrt(2*m)
362
raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
363
if (m != len(y)) or (m != len(w)):
364
raise TypeError, 'Lengths of the first three arguments'\
365
' (x,y,w) must be equal'
367
raise TypeError, 'Given degree of the spline (k=%d) is not supported.'\
369
if m<=k: raise TypeError, 'm>k must hold'
370
if xb is None: xb=x[0]
371
if xe is None: xe=x[-1]
372
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
376
if t is None: raise TypeError, 'Knots must be given for task=-1'
379
n,t,c,fp,wrk,iwrk,ier = dfitpack.percur_smth0(x,y,w,k=k,s=s)
383
wrk=_percur_cache['wrk']
384
iwrk=_percur_cache['iwrk']
387
raise ValueError, 'task=1 can only be called after task=0'
388
n,t,c,fp,wrk,iwrk,ier = dfitpack.percur_smth1(x,y,w,n,t,wrk,iwrk,
391
n,t,c,fp,ier = dfitpack.percur_lsq(x,y,w,t,k=k)
394
_percur_cache['wrk']=wrk
395
_percur_cache['iwrk']=iwrk
399
spl = spline.UnivariateSpline(x,y,w,[xb,xe],k=k,s=s)
402
spl = _splrep_cache['spl']
404
raise ValueError, 'task=1 can only be called after task=0'
405
spl.set_smoothing_factor(s)
409
spl = spline.LSQUnivariateSpline(x,y,t,w,[xb,xe],k=k)
411
_splrep_cache['spl'] = spl
412
x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = spl._data
413
tck = (t[:n],c[:n-k-1],k)
414
if ier<=0 and not quiet:
415
print _iermess[ier][0]
416
print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
417
if ier>0 and not full_output:
419
print "Warning: "+_iermess[ier][0]
422
raise _iermess[ier][1],_iermess[ier][0]
424
raise _iermess['unknown'][1],_iermess['unknown'][0]
427
return tck,fp,ier,_iermess[ier][0]
429
return tck,fp,ier,_iermess['unknown'][0]
433
def splev(x,tck,der=0):
434
"""Evaulate a B-spline and its derivatives.
437
Given the knots and coefficients of a B-spline representation, evaluate
438
the value of the smoothing polynomial and it's derivatives.
439
This is a wrapper around the FORTRAN routines splev and splder of FITPACK.
442
x (u) -- a 1-D array of points at which to return the value of the
443
smoothed spline or its derivatives. If tck was returned from
444
splprep, then the parameter values, u should be given.
445
tck -- A sequence of length 3 returned by splrep or splprep containg the
446
knots, coefficients, and degree of the spline.
447
der -- The order of derivative of the spline to compute (must be less than
451
y -- an array of values representing the spline function or curve.
452
If tck was returned from splrep, then this is a list of arrays
453
representing the curve in N-dimensional space.
456
splprep, splrep, sproot, spalde, splint - evaluation, roots, integral
457
bisplrep, bisplev - bivariate splines
458
UnivariateSpline, BivariateSpline - an alternative wrapping
459
of the FITPACK functions
463
c[0][0] # check if c is >1-d
464
parametric = True # it is
465
except TypeError: parametric = False # c is 1-d
467
return map(lambda c,x=x,t=t,k=k,der=der:splev(x,[t,c,k],der),c)
470
raise ValueError,"0<=der=%d<=k=%d must hold"%(der,k)
475
y,ier=dfitpack.splder(t,c,k,x,der)
477
y,ier=dfitpack.splev(t,c,k,x)
478
if ier==10: raise ValueError,"Invalid input data"
479
if ier: raise TypeError,"An unkown error occurred"
480
if len(y)>1: return y
483
def splint(a,b,tck,full_output=False):
484
"""Evaluate the definite integral of a B-spline.
487
Given the knots and coefficients of a B-spline, evaluate the definite
488
integral of the smoothing polynomial between two given points.
491
a, b -- The end-points of the integration interval.
492
tck -- A length 3 sequence describing the given spline (See splev).
493
full_output -- True to return optional output.
495
Outputs: (integral, {wrk})
496
integral -- The resulting integral.
497
wrk -- An array containing the integrals of the normalized B-splines
498
defined on the set of knots.
501
splprep, splrep, sproot, spalde, splev - evaluation, roots, integral
502
bisplrep, bisplev - bivariate splines
503
UnivariateSpline, BivariateSpline - an alternative wrapping
504
of the FITPACK functions
508
c[0][0] # check if c is >1-d
509
parametric = True # it is
510
except TypeError: parametric = False # c is 1-d
512
return map(lambda c,a=a,b=b,t=t,k=k:splint(a,b,[t,c,k]),c)
516
aint,wrk = dfitpack.splint(t,c,k,a,b)
517
if full_output: return aint,wrk
520
def sproot(tck,mest=10):
521
"""Find the roots of a cubic B-spline.
524
Given the knots (>=8) and coefficients of a cubic B-spline return the
528
tck -- A length 3 sequence describing the given spline (See splev).
529
The number of knots must be >= 8. The knots must be a
530
montonically increasing sequence.
531
mest -- An estimate of the number of zeros (Default is 10).
534
zeros -- An array giving the roots of the spline.
537
splprep, splrep, splint, spalde, splev - evaluation, roots, integral
538
bisplrep, bisplev - bivariate splines
539
UnivariateSpline, BivariateSpline - an alternative wrapping
540
of the FITPACK functions
544
raise ValueError, "Sproot only valid for cubic bsplines"
546
c[0][0] # check if c is >1-d
547
parametric = True # it is
548
except TypeError: parametric = False # c is 1-d
550
return map(lambda c,t=t,k=k,mest=mest:sproot([t,c,k],mest),c)
553
raise TypeError,"The number of knots be >=8"
556
z,m,ier=dfitpack.sproot(t,c,mest)
558
raise TypeError,"Invalid input data. t1<=..<=t4<t5<..<tn-3<=..<=tn" \
560
if ier==0: return z[0:m]
562
print "Warning: the number of zeros exceeds mest"
564
raise TypeError,"Unknown error"
567
"""Evaluate all derivatives of a B-spline.
570
Given the knots and coefficients of a cubic B-spline compute all
571
derivatives up to order k at a point (or set of points).
574
tck -- A length 3 sequence describing the given spline (See splev).
575
x -- A point or a set of points at which to evaluate the derivatives.
576
Note that t(k) <= x <= t(n-k+1) must hold for each x.
579
results -- An array (or a list of arrays) containing all derivatives
580
up to order k inclusive for each point x.
583
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
584
bisplrep, bisplev - bivariate splines
585
UnivariateSpline, BivariateSpline - an alternative wrapping
586
of the FITPACK functions
590
c[0][0] # check if c is >1-d
591
parametric = True # it is
592
except TypeError: parametric = False # c is 1-d
594
return map(lambda c,x=x,t=t,k=k:spalde(x,[t,c,k]),c)
601
return map(lambda x,tck=tck:spalde(x,tck),x)
604
d,ier=dfitpack.spalde(t,c,k,x[0])
607
raise TypeError,"Invalid input data. t(k)<=x<=t(n-k+1) must hold."
608
raise TypeError,"Unknown error"
610
def insert(x,tck,m=1,per=0):
611
"""Insert knots into a B-spline.
615
Given the knots and coefficients of a B-spline representation, create a
616
new B-spline with a knot inserted m times at point x.
617
This is a wrapper around the FORTRAN routine insert of FITPACK.
621
x (u) -- A 1-D point at which to insert a new knot(s). If tck was returned
622
from splprep, then the parameter values, u should be given.
623
tck -- A sequence of length 3 returned by splrep or splprep containg the
624
knots, coefficients, and degree of the spline.
625
m -- The number of times to insert the given knot (its multiplicity).
626
per -- If non-zero, input spline is considered periodic.
630
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
631
coefficients, and the degree of the new spline.
634
t(k+1) <= x <= t(n-k), where k is the degree of the spline.
635
In case of a periodic spline (per != 0) there must be
636
either at least k interior knots t(j) satisfying t(k+1)<t(j)<=x
637
or at least k interior knots t(j) satisfying x<=t(j)<t(n-k).
648
tt, cc_val, kk = insert(x, [t, c_vals, k], m)
659
tt, nn, cc, ier = dfitpack.insert(per, t, n, c, k, x)
660
if ier==10: raise ValueError,"Invalid input data"
661
if ier: raise TypeError,"An error occurred"
669
def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,
670
s=None,eps=1e-16,tx=None,ty=None,full_output=0,
671
nxest=None,nyest=None,quiet=1):
672
"""Find a bivariate B-spline representation of a surface.
675
Given a set of data points (x[i], y[i], z[i]) representing a surface
676
z=f(x,y), compute a B-spline representation of the surface.
679
x, y, z -- Rank-1 arrays of data points.
680
w -- Rank-1 array of weights. By default w=ones(len(x)).
681
xb, xe -- End points of approximation interval in x.
682
yb, ye -- End points of approximation interval in y.
683
By default xb, xe, yb, ye = x.min(), x.max(), y.min(), y.max()
684
kx, ky -- The degrees of the spline (1 <= kx, ky <= 5). Third order
685
(kx=ky=3) is recommended.
686
task -- If task=0, find knots in x and y and coefficients for a given
688
If task=1, find knots and coefficients for another value of the
689
smoothing factor, s. bisplrep must have been previously called
690
with task=0 or task=1.
691
If task=-1, find coefficients for a given set of knots tx, ty.
692
s -- A non-negative smoothing factor. If weights correspond
693
to the inverse of the standard-deviation of the errors in z,
694
then a good s-value should be found in the range
695
(m-sqrt(2*m),m+sqrt(2*m)) where m=len(x)
696
eps -- A threshold for determining the effective rank of an
697
over-determined linear system of equations (0 < eps < 1)
698
--- not likely to need changing.
699
tx, ty -- Rank-1 arrays of the knots of the spline for task=-1
700
full_output -- Non-zero to return optional outputs.
701
nxest, nyest -- Over-estimates of the total number of knots.
702
If None then nxest = max(kx+sqrt(m/2),2*kx+3),
703
nyest = max(ky+sqrt(m/2),2*ky+3)
704
quiet -- Non-zero to suppress printing of messages.
706
Outputs: (tck, {fp, ier, msg})
707
tck -- A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and
708
coefficients (c) of the bivariate B-spline representation of the
709
surface along with the degree of the spline.
711
fp -- The weighted sum of squared residuals of the spline approximation.
712
ier -- An integer flag about splrep success. Success is indicated if
713
ier<=0. If ier in [1,2,3] an error occurred but was not raised.
714
Otherwise an error is raised.
715
msg -- A message corresponding to the integer flag, ier.
718
SEE bisplev to evaluate the value of the B-spline given its tck
722
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
723
UnivariateSpline, BivariateSpline - an alternative wrapping
724
of the FITPACK functions
726
x,y,z=map(myasarray,[x,y,z])
727
x,y,z=map(ravel,[x,y,z]) # ensure 1-d arrays.
729
if not (m==len(y)==len(z)):
730
raise TypeError, 'len(x)==len(y)==len(z) must hold.'
731
if w is None: w=ones(m,float)
734
raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
735
if xb is None: xb=x.min()
736
if xe is None: xe=x.max()
737
if yb is None: yb=y.min()
738
if ye is None: ye=y.max()
739
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
740
if s is None: s=m-sqrt(2*m)
742
if tx is None: raise TypeError, 'Knots_x must be given for task=-1'
743
_surfit_cache['tx']=myasarray(tx)
744
if ty is None: raise TypeError, 'Knots_y must be given for task=-1'
745
_surfit_cache['ty']=myasarray(ty)
747
raise TypeError,'There must be at least 2*kx+2 knots_x for task=-1'
749
raise TypeError,'There must be at least 2*ky+2 knots_y for task=-1'
750
if not ((1<=kx<=5) and (1<=ky<=5)):
751
raise TypeError, 'Given degree of the spline (kx,ky=%d,%d) is not \
752
supported. (1<=k<=5)'%(kx,ky)
753
if m<(kx+1)*(ky+1): raise TypeError, 'm>=(kx+1)(ky+1) must hold'
754
if nxest is None: nxest=kx+sqrt(m/2)
755
if nyest is None: nyest=ky+sqrt(m/2)
756
nxest,nyest=max(nxest,2*kx+3),max(nyest,2*ky+3)
758
nxest=int(kx+sqrt(3*m))
759
nyest=int(ky+sqrt(3*m))
761
nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth0(x,y,z,w,xb,xe,yb,ye,
763
nxest=nxest,nyest=nyest)
766
tx,ty=_surfit_cache['tx'],_surfit_cache['ty']
767
nx,ny= _surfit_cache['nx'],_surfit_cache['ny']
768
wrk1 = _surfit_cache['wrk1']
770
raise ValueError, 'task=1 can only be called after task=0'
771
nx,tx,ny,ty,c,fp,wrk1,ier=dfitpack.surfit_smth1(x,y,z,nx,tx,ny,ty,wrk1,
772
w,xb,xe,yb,ye,kx,ky,s=s,eps=eps,
773
nxest=nxest,nyest=nyest)
777
tx,ty,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx,ty,w,
780
nxest=nxest,nyest=nyest)
782
tx,ty,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx,ty,w,\
785
nxest=nxest,nyest=nyest)
787
_surfit_cache['tx']=tx
788
_surfit_cache['ty']=ty
789
_surfit_cache['nx']=nx
790
_surfit_cache['ny']=ny
791
_surfit_cache['wrk1']=wrk1
792
tck=[tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)],kx,ky]
793
if ier<=0 and not quiet:
794
print _iermess2[ier][0]
795
print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
797
ierm=min(11,max(-3,ier))
798
if ierm>0 and not full_output:
799
if ier in [1,2,3,4,5]:
800
print "Warning: "+_iermess2[ierm][0]
801
print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
805
raise _iermess2[ierm][1],_iermess2[ierm][0]
807
raise _iermess2['unknown'][1],_iermess2['unknown'][0]
810
return tck,fp,ier,_iermess2[ierm][0]
812
return tck,fp,ier,_iermess2['unknown'][0]
816
def bisplev(x,y,tck,dx=0,dy=0):
817
"""Evaluate a bivariate B-spline and its derivatives.
820
Return a rank-2 array of spline function values (or spline derivative
821
values) at points given by the cross-product of the rank-1 arrays x and y.
822
In special cases, return an array or just a float if either x or y or
826
x, y -- Rank-1 arrays specifying the domain over which to evaluate the
827
spline or its derivative.
828
tck -- A sequence of length 5 returned by bisplrep containing the knot
829
locations, the coefficients, and the degree of the spline:
831
dx, dy -- The orders of the partial derivatives in x and y respectively.
834
vals -- The B-pline or its derivative evaluated over the set formed by
835
the cross-product of x and y.
838
SEE bisprep to generate the tck representation.
841
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
842
UnivariateSpline, BivariateSpline - an alternative wrapping
843
of the FITPACK functions
846
if not (0<=dx<kx): raise ValueError,"0<=dx=%d<kx=%d must hold"%(dx,kx)
847
if not (0<=dy<ky): raise ValueError,"0<=dy=%d<ky=%d must hold"%(dy,ky)
848
x,y=map(myasarray,[x,y])
849
if (len(x.shape) != 1) or (len(y.shape) != 1):
850
raise ValueError, "First two entries should be rank-1 arrays."
852
z,ier=dfitpack.parder(tx,ty,c,kx,ky,dx,dy,x,y)
854
z,ier=dfitpack.bispev(tx,ty,c,kx,ky,x,y)
855
if ier==10: raise ValueError,"Invalid input data"
856
if ier: raise TypeError,"An error occurred"
857
z.shape=len(x),len(y)
858
if len(z)>1: return z
859
if len(z[0])>1: return z[0]
862
if __name__ == "__main__":
865
if len(sys.argv[1:])>0:
866
runtest=map(string.atoi,sys.argv[1:])
869
return dot(transpose(x),x)
871
if d is None: return "sin"
872
if x is None: return "sin(x)"
873
if d%4 == 0: return sin(x)
874
if d%4 == 1: return cos(x)
875
if d%4 == 2: return -sin(x)
876
if d%4 == 3: return -cos(x)
877
def f2(x,y=0,dx=0,dy=0):
878
if x is None: return "sin(x+y)"
880
if d%4 == 0: return sin(x+y)
881
if d%4 == 1: return cos(x+y)
882
if d%4 == 2: return -sin(x+y)
883
if d%4 == 3: return -cos(x+y)
884
def test1(f=f1,per=0,s=0,a=0,b=2*pi,N=20,at=0,xb=None,xe=None):
887
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
888
x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
892
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
897
nd.append(norm2(f(t,d)-splev(t,tck,d)))
899
print "\nf = %s s=S_k(x;t,c) x in [%s, %s] > [%s, %s]"%(f(None),
900
`round(xb,3)`,`round(xe,3)`,
901
`round(a,3)`,`round(b,3)`)
902
if at: str="at knots"
903
else: str="at the middle of nodes"
904
print " per=%d s=%s Evaluation %s"%(per,`s`,str)
905
print " k : |f-s|^2 |f'-s'| |f''-.. |f'''-. |f''''- |f'''''"
913
def test2(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
914
ia=0,ib=2*pi,dx=0.2*pi):
917
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
921
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
922
nk.append([splint(ia,ib,tck),spalde(dx,tck)])
923
print "\nf = %s s=S_k(x;t,c) x in [%s, %s] > [%s, %s]"%(f(None),
924
`round(xb,3)`,`round(xe,3)`,
925
`round(a,3)`,`round(b,3)`)
926
print " per=%d s=%s N=%d [a, b] = [%s, %s] dx=%s"%(per,`s`,N,`round(ia,3)`,`round(ib,3)`,`round(dx,3)`)
927
print " k : int(s,[a,b]) Int.Error Rel. error of s^(d)(dx) d = 0, .., k"
932
put(" %d %s%.8f %.1e "%(k,sr,abs(r[0]),
933
abs(r[0]-(f(ib,-1)-f(ia,-1)))))
936
put(" %.1e "%(abs(1-dr/f(dx,d))))
940
def test3(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
941
ia=0,ib=2*pi,dx=0.2*pi):
944
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
947
print " k : Roots of s(x) approx %s x in [%s,%s]:"%\
948
(f(None),`round(a,3)`,`round(b,3)`)
950
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
951
print ' %d : %s'%(k,`sproot(tck).tolist()`)
952
def test4(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
953
ia=0,ib=2*pi,dx=0.2*pi):
956
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
957
x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
960
print " u = %s N = %d"%(`round(dx,3)`,N)
961
print " k : [x(u), %s(x(u))] Error of splprep Error of splrep "%(f(0,None))
963
tckp,u=splprep([x,v],s=s,per=per,k=k)
964
tck=splrep(x,v,s=s,per=per,k=k)
966
print " %d : %s %.1e %.1e"%\
967
(k,`map(lambda x:round(x,3),uv)`,
969
abs(splev(uv[0],tck)-f(uv[0])))
970
print "Derivatives of parametric cubic spline at u (first function):"
972
tckp,u=splprep([x,v],s=s,per=per,k=k)
973
for d in range(1,k+1):
975
put(" %s "%(`uv[0]`))
978
x,y=map(myasarray,[x,y])
979
xy=array(map(lambda x,y:map(None,len(y)*[x],y),x,len(x)*[y]))
981
xy.shape=sh[0]*sh[1],sh[2]
983
def test5(f=f2,kx=3,ky=3,xb=0,xe=2*pi,yb=0,ye=2*pi,Nx=20,Ny=20,s=0):
984
x=xb+(xe-xb)*arange(Nx+1,dtype=float)/float(Nx)
985
y=yb+(ye-yb)*arange(Ny+1,dtype=float)/float(Ny)
987
tck=bisplrep(xy[0],xy[1],f(xy[0],xy[1]),s=s,kx=kx,ky=ky)
988
tt=[tck[0][kx:-kx],tck[1][ky:-ky]]
989
t2=makepairs(tt[0],tt[1])
990
v1=bisplev(tt[0],tt[1],tck)
992
v2.shape=len(tt[0]),len(tt[1])
993
print norm2(ravel(v1-v2))
996
******************************************
997
\tTests of splrep and splev
998
******************************************"""
1005
test1(b=1.5*pi,xe=2*pi,per=1,s=1e-1)
1008
******************************************
1009
\tTests of splint and spalde
1010
******************************************"""
1013
test2(ia=0.2*pi,ib=pi)
1014
test2(ia=0.2*pi,ib=pi,N=50)
1017
******************************************
1019
******************************************"""
1021
print "Note that if k is not 3, some roots are missed or incorrect"
1024
******************************************
1025
\tTests of splprep, splrep, and splev
1026
******************************************"""
1031
******************************************
1032
\tTests of bisplrep, bisplev
1033
******************************************"""
added
removed
scipy/sandbox/spline/fitpack.py
