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.
22
Running test programs:
23
$ python fitpack.py 1 3 # run test programs 1, and 3
24
$ python fitpack.py # run all available test programs
26
TODO: Make interfaces to the following fitpack functions:
27
For univariate splines: cocosp, concon, fourco, insert
28
For bivariate splines: profil, regrid, parsur, surev
31
__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
32
'bisplrep', 'bisplev', 'insert']
33
__version__ = "$Revision: 2762 $"[10:-1]
35
from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
36
dot, sin, cos, pi, arange, empty, int32
37
myasarray = atleast_1d
39
# Try to replace _fitpack interface with
40
# f2py-generated version
44
The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
46
The spline is an interpolating spline (fp=0)""",None],
48
The spline is weighted least-squares polynomial of degree k.
49
fp gives the upper bound fp0 for the smoothing factor s""",None],
51
The required storage space exceeds the available storage space.
52
Probable causes: data (x,y) size is too small or smoothing parameter s is too small (fp>s).""",ValueError],
54
A theoretically impossible results when finding a smoothin spline
55
with fp = s. Probably causes: s too small. (abs(fp-s)/s>0.001)""",ValueError],
57
The maximal number of iterations (20) allowed for finding smoothing
58
spline with fp=s has been reached. Probably causes: s too small.
59
(abs(fp-s)/s>0.001)""",ValueError],
61
Error on input data""",ValueError],
63
An error occured""",TypeError]}
66
The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
68
The spline is an interpolating spline (fp=0)""",None],
70
The spline is weighted least-squares polynomial of degree kx and ky.
71
fp gives the upper bound fp0 for the smoothing factor s""",None],
73
Warning. The coefficients of the spline have been computed as the minimal
74
norm least-squares solution of a rank deficient system.""",None],
76
The required storage space exceeds the available storage space.
77
Probably causes: nxest or nyest too small or s is too small. (fp>s)""",ValueError],
79
A theoretically impossible results when finding a smoothin spline
80
with fp = s. Probably causes: s too small or badly chosen eps.
81
(abs(fp-s)/s>0.001)""",ValueError],
83
The maximal number of iterations (20) allowed for finding smoothing
84
spline with fp=s has been reached. Probably causes: s too small.
85
(abs(fp-s)/s>0.001)""",ValueError],
87
No more knots can be added because the number of B-spline coefficients
88
already exceeds the number of data points m. Probably causes: either
89
s or m too small. (fp>s)""",ValueError],
91
No more knots can be added because the additional knot would coincide
92
with an old one. Probably cause: s too small or too large a weight
93
to an inaccurate data point. (fp>s)""",ValueError],
95
Error on input data""",ValueError],
97
rwrk2 too small, i.e. there is not enough workspace for computing
98
the minimal least-squares solution of a rank deficient system of linear
99
equations.""",ValueError],
101
An error occured""",TypeError]}
103
_parcur_cache = {'t': array([],float), 'wrk': array([],float),
104
'iwrk':array([],int32), 'u': array([],float),'ub':0,'ue':1}
106
def splprep(x,w=None,u=None,ub=None,ue=None,k=3,task=0,s=None,t=None,
107
full_output=0,nest=None,per=0,quiet=1):
108
"""Find the B-spline representation of an N-dimensional curve.
112
Given a list of N rank-1 arrays, x, which represent a curve in N-dimensional
113
space parametrized by u, find a smooth approximating spline curve g(u).
114
Uses the FORTRAN routine parcur from FITPACK
118
x -- A list of sample vector arrays representing the curve.
119
u -- An array of parameter values. If not given, these values are
120
calculated automatically as (M = len(x[0])):
122
v[i] = v[i-1] + distance(x[i],x[i-1])
124
ub, ue -- The end-points of the parameters interval. Defaults to
126
k -- Degree of the spline. Cubic splines are recommended. Even values of
127
k should be avoided especially with a small s-value.
129
task -- If task==0 find t and c for a given smoothing factor, s.
130
If task==1 find t and c for another value of the smoothing factor,
131
s. There must have been a previous call with task=0 or task=1
132
for the same set of data.
133
If task=-1 find the weighted least square spline for a given set of
135
s -- A smoothing condition. The amount of smoothness is determined by
136
satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
137
g(x) is the smoothed interpolation of (x,y). The user can use s to
138
control the tradeoff between closeness and smoothness of fit. Larger
139
s means more smoothing while smaller values of s indicate less
140
smoothing. Recommended values of s depend on the weights, w. If the
141
weights represent the inverse of the standard-deviation of y, then a
142
good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
143
where m is the number of datapoints in x, y, and w.
144
t -- The knots needed for task=-1.
145
full_output -- If non-zero, then return optional outputs.
146
nest -- An over-estimate of the total number of knots of the spline to
147
help in determining the storage space. By default nest=m/2.
148
Always large enough is nest=m+k+1.
149
per -- If non-zero, data points are considered periodic with period
150
x[m-1] - x[0] and a smooth periodic spline approximation is returned.
151
Values of y[m-1] and w[m-1] are not used.
152
quiet -- Non-zero to suppress messages.
154
Outputs: (tck, u, {fp, ier, msg})
156
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
157
coefficients, and the degree of the spline.
158
u -- An array of the values of the parameter.
160
fp -- The weighted sum of squared residuals of the spline approximation.
161
ier -- An integer flag about splrep success. Success is indicated
162
if ier<=0. If ier in [1,2,3] an error occurred but was not raised.
163
Otherwise an error is raised.
164
msg -- A message corresponding to the integer flag, ier.
168
SEE splev for evaluation of the spline and its derivatives.
171
splrep, splev, sproot, spalde, splint - evaluation, roots, integral
172
bisplrep, bisplev - bivariate splines
173
UnivariateSpline, BivariateSpline - an alternative wrapping
174
of the FITPACK functions
177
_parcur_cache = {'t': array([],float), 'wrk': array([],float),
178
'iwrk':array([],int32),'u': array([],float),
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)
193
if ub is None: _parcur_cache['ub']=u[0]
194
else: _parcur_cache['ub']=ub
195
if ue is None: _parcur_cache['ue']=u[-1]
196
else: _parcur_cache['ue']=ue
197
else: _parcur_cache['u']=zeros(m,float)
198
if not (1<=k<=5): raise TypeError, '1<=k=%d<=5 must hold'%(k)
199
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
200
if (not len(w)==m) or (ipar==1 and (not len(u)==m)):
201
raise TypeError,'Mismatch of input dimensions'
202
if s is None: s=m-sqrt(2*m)
203
if t is None and task==-1: raise TypeError, 'Knots must be given for task=-1'
205
_parcur_cache['t']=myasarray(t)
206
n=len(_parcur_cache['t'])
207
if task==-1 and n<2*k+2:
208
raise TypeError, 'There must be at least 2*k+2 knots for task=-1'
209
if m<=k: raise TypeError, 'm>k must hold'
210
if nest is None: nest=m+2*k
212
if (task>=0 and s==0) or (nest<0):
217
ub=_parcur_cache['ub']
218
ue=_parcur_cache['ue']
220
wrk=_parcur_cache['wrk']
221
iwrk=_parcur_cache['iwrk']
222
t,c,o=_fitpack._parcur(ravel(transpose(x)),w,u,ub,ue,k,task,ipar,s,t,
224
_parcur_cache['u']=o['u']
225
_parcur_cache['ub']=o['ub']
226
_parcur_cache['ue']=o['ue']
228
_parcur_cache['wrk']=o['wrk']
229
_parcur_cache['iwrk']=o['iwrk']
230
ier,fp,n=o['ier'],o['fp'],len(t)
233
tcku = [t,list(c),k],u
234
if ier<=0 and not quiet:
235
print _iermess[ier][0]
236
print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
237
if ier>0 and not full_output:
239
print "Warning: "+_iermess[ier][0]
242
raise _iermess[ier][1],_iermess[ier][0]
244
raise _iermess['unknown'][1],_iermess['unknown'][0]
247
return tcku,fp,ier,_iermess[ier][0]
249
return tcku,fp,ier,_iermess['unknown'][0]
253
_curfit_cache = {'t': array([],float), 'wrk': array([],float),
254
'iwrk':array([],int32)}
255
def splrep(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
256
full_output=0,per=0,quiet=1):
257
"""Find the B-spline representation of 1-D curve.
261
Given the set of data points (x[i], y[i]) determine a smooth spline
262
approximation of degree k on the interval xb <= x <= xe. The coefficients,
263
c, and the knot points, t, are returned. Uses the FORTRAN routine
268
x, y -- The data points defining a curve y = f(x).
269
w -- Strictly positive rank-1 array of weights the same length as x and y.
270
The weights are used in computing the weighted least-squares spline
271
fit. If the errors in the y values have standard-deviation given by the
272
vector d, then w should be 1/d. Default is ones(len(x)).
273
xb, xe -- The interval to fit. If None, these default to x[0] and x[-1]
275
k -- The order of the spline fit. It is recommended to use cubic splines.
276
Even order splines should be avoided especially with small s values.
278
task -- If task==0 find t and c for a given smoothing factor, s.
279
If task==1 find t and c for another value of the
280
smoothing factor, s. There must have been a previous
281
call with task=0 or task=1 for the same set of data
282
(t will be stored an used internally)
283
If task=-1 find the weighted least square spline for
284
a given set of knots, t. These should be interior knots
285
as knots on the ends will be added automatically.
286
s -- A smoothing condition. The amount of smoothness is determined by
287
satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
288
g(x) is the smoothed interpolation of (x,y). The user can use s to
289
control the tradeoff between closeness and smoothness of fit. Larger
290
s means more smoothing while smaller values of s indicate less
291
smoothing. Recommended values of s depend on the weights, w. If the
292
weights represent the inverse of the standard-deviation of y, then a
293
good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
294
where m is the number of datapoints in x, y, and w.
295
default : s=m-sqrt(2*m) if weights are supplied.
296
s = 0.0 (interpolating) if no weights are supplied.
297
t -- The knots needed for task=-1. If given then task is automatically
299
full_output -- If non-zero, then return optional outputs.
300
per -- If non-zero, data points are considered periodic with period
301
x[m-1] - x[0] and a smooth periodic spline approximation is returned.
302
Values of y[m-1] and w[m-1] are not used.
303
quiet -- Non-zero to suppress messages.
305
Outputs: (tck, {fp, ier, msg})
307
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
308
coefficients, and the degree of the spline.
310
fp -- The weighted sum of squared residuals of the spline approximation.
311
ier -- An integer flag about splrep success. Success is indicated if
312
ier<=0. If ier in [1,2,3] an error occurred but was not raised.
313
Otherwise an error is raised.
314
msg -- A message corresponding to the integer flag, ier.
318
See splev for evaluation of the spline and its derivatives.
322
x = linspace(0, 10, 10)
325
x2 = linspace(0, 10, 200)
327
plot(x, y, 'o', x2, y2)
330
splprep, splev, sproot, spalde, splint - evaluation, roots, integral
331
bisplrep, bisplev - bivariate splines
332
UnivariateSpline, BivariateSpline - an alternative wrapping
333
of the FITPACK functions
337
x,y=map(myasarray,[x,y])
341
if s is None: s = 0.0
344
if s is None: s = m-sqrt(2*m)
345
if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
346
if (m != len(y)) or (m != len(w)):
347
raise TypeError, 'Lengths of the first three arguments (x,y,w) must be equal'
349
raise TypeError, 'Given degree of the spline (k=%d) is not supported. (1<=k<=5)'%(k)
350
if m<=k: raise TypeError, 'm>k must hold'
351
if xb is None: xb=x[0]
352
if xe is None: xe=x[-1]
353
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
357
if t is None: raise TypeError, 'Knots must be given for task=-1'
359
_curfit_cache['t'] = empty((numknots + 2*k+2,),float)
360
_curfit_cache['t'][k+1:-k-1] = t
361
nest = len(_curfit_cache['t'])
364
nest = max(m+2*k,2*k+3)
366
nest = max(m+k+1,2*k+3)
367
t = empty((nest,),float)
368
_curfit_cache['t'] = t
370
if per: _curfit_cache['wrk'] = empty((m*(k+1)+nest*(8+5*k),),float)
371
else: _curfit_cache['wrk'] = empty((m*(k+1)+nest*(7+3*k),),float)
372
_curfit_cache['iwrk'] = empty((nest,),int32)
375
wrk=_curfit_cache['wrk']
376
iwrk=_curfit_cache['iwrk']
378
raise TypeError, "must call with task=1 only after"\
379
" call with task=0,-1"
381
n,c,fp,ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk, xb, xe, k, s)
383
n,c,fp,ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s)
384
tck = (t[:n],c[:n],k)
385
if ier<=0 and not quiet:
386
print _iermess[ier][0]
387
print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
388
if ier>0 and not full_output:
390
print "Warning: "+_iermess[ier][0]
393
raise _iermess[ier][1],_iermess[ier][0]
395
raise _iermess['unknown'][1],_iermess['unknown'][0]
398
return tck,fp,ier,_iermess[ier][0]
400
return tck,fp,ier,_iermess['unknown'][0]
404
def _ntlist(l): # return non-trivial list
406
#if len(l)>1: return l
409
def splev(x,tck,der=0):
410
"""Evaulate a B-spline and its derivatives.
414
Given the knots and coefficients of a B-spline representation, evaluate
415
the value of the smoothing polynomial and it's derivatives.
416
This is a wrapper around the FORTRAN routines splev and splder of FITPACK.
420
x (u) -- a 1-D array of points at which to return the value of the
421
smoothed spline or its derivatives. If tck was returned from
422
splprep, then the parameter values, u should be given.
423
tck -- A sequence of length 3 returned by splrep or splprep containg the
424
knots, coefficients, and degree of the spline.
425
der -- The order of derivative of the spline to compute (must be less than
430
y -- an array of values representing the spline function or curve.
431
If tck was returned from splrep, then this is a list of arrays
432
representing the curve in N-dimensional space.
435
splprep, splrep, sproot, spalde, splint - evaluation, roots, integral
436
bisplrep, bisplev - bivariate splines
437
UnivariateSpline, BivariateSpline - an alternative wrapping
438
of the FITPACK functions
447
return map(lambda c,x=x,t=t,k=k,der=der:splev(x,[t,c,k],der),c)
450
raise ValueError,"0<=der=%d<=k=%d must hold"%(der,k)
452
y,ier=_fitpack._spl_(x,der,t,c,k)
453
if ier==10: raise ValueError,"Invalid input data"
454
if ier: raise TypeError,"An error occurred"
455
if len(y)>1: return y
458
def splint(a,b,tck,full_output=0):
459
"""Evaluate the definite integral of a B-spline.
463
Given the knots and coefficients of a B-spline, evaluate the definite
464
integral of the smoothing polynomial between two given points.
468
a, b -- The end-points of the integration interval.
469
tck -- A length 3 sequence describing the given spline (See splev).
470
full_output -- Non-zero to return optional output.
472
Outputs: (integral, {wrk})
474
integral -- The resulting integral.
475
wrk -- An array containing the integrals of the normalized B-splines defined
480
splprep, splrep, sproot, spalde, splev - evaluation, roots, integral
481
bisplrep, bisplev - bivariate splines
482
UnivariateSpline, BivariateSpline - an alternative wrapping
483
of the FITPACK functions
492
return _ntlist(map(lambda c,a=a,b=b,t=t,k=k:splint(a,b,[t,c,k]),c))
494
aint,wrk=_fitpack._splint(t,c,k,a,b)
495
if full_output: return aint,wrk
498
def sproot(tck,mest=10):
499
"""Find the roots of a cubic B-spline.
503
Given the knots (>=8) and coefficients of a cubic B-spline return the
508
tck -- A length 3 sequence describing the given spline (See splev).
509
The number of knots must be >= 8. The knots must be a montonically
511
mest -- An estimate of the number of zeros (Default is 10).
515
zeros -- An array giving the roots of the spline.
518
splprep, splrep, splint, spalde, splev - evaluation, roots, integral
519
bisplrep, bisplev - bivariate splines
520
UnivariateSpline, BivariateSpline - an alternative wrapping
521
of the FITPACK functions
532
return _ntlist(map(lambda c,t=t,k=k,mest=mest:sproot([t,c,k],mest),c))
535
raise TypeError,"The number of knots %d>=8"%(len(t))
536
z,ier=_fitpack._sproot(t,c,k,mest)
538
raise TypeError,"Invalid input data. t1<=..<=t4<t5<..<tn-3<=..<=tn must hold."
541
print "Warning: the number of zeros exceeds mest"
543
raise TypeError,"Unknown error"
546
"""Evaluate all derivatives of a B-spline.
550
Given the knots and coefficients of a cubic B-spline compute all
551
derivatives up to order k at a point (or set of points).
555
tck -- A length 3 sequence describing the given spline (See splev).
556
x -- A point or a set of points at which to evaluate the derivatives.
557
Note that t(k) <= x <= t(n-k+1) must hold for each x.
561
results -- An array (or a list of arrays) containing all derivatives
562
up to order k inclusive for each point x.
565
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
566
bisplrep, bisplev - bivariate splines
567
UnivariateSpline, BivariateSpline - an alternative wrapping
568
of the FITPACK functions
577
return _ntlist(map(lambda c,x=x,t=t,k=k:spalde(x,[t,c,k]),c))
584
return map(lambda x,tck=tck:spalde(x,tck),x)
585
d,ier=_fitpack._spalde(t,c,k,x[0])
588
raise TypeError,"Invalid input data. t(k)<=x<=t(n-k+1) must hold."
589
raise TypeError,"Unknown error"
591
#def _curfit(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
592
# full_output=0,nest=None,per=0,quiet=1):
594
_surfit_cache = {'tx': array([],float),'ty': array([],float),
595
'wrk': array([],float), 'iwrk':array([],int32)}
596
def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,
597
s=None,eps=1e-16,tx=None,ty=None,full_output=0,
598
nxest=None,nyest=None,quiet=1):
599
"""Find a bivariate B-spline representation of a surface.
603
Given a set of data points (x[i], y[i], z[i]) representing a surface
604
z=f(x,y), compute a B-spline representation of the surface.
608
x, y, z -- Rank-1 arrays of data points.
609
w -- Rank-1 array of weights. By default w=ones(len(x)).
610
xb, xe -- End points of approximation interval in x.
611
yb, ye -- End points of approximation interval in y.
612
By default xb, xe, yb, ye = x.min(), x.max(), y.min(), y.max()
613
kx, ky -- The degrees of the spline (1 <= kx, ky <= 5). Third order
614
(kx=ky=3) is recommended.
615
task -- If task=0, find knots in x and y and coefficients for a given
617
If task=1, find knots and coefficients for another value of the
618
smoothing factor, s. bisplrep must have been previously called
619
with task=0 or task=1.
620
If task=-1, find coefficients for a given set of knots tx, ty.
621
s -- A non-negative smoothing factor. If weights correspond
622
to the inverse of the standard-deviation of the errors in z,
623
then a good s-value should be found in the range
624
(m-sqrt(2*m),m+sqrt(2*m)) where m=len(x)
625
eps -- A threshold for determining the effective rank of an
626
over-determined linear system of equations (0 < eps < 1)
627
--- not likely to need changing.
628
tx, ty -- Rank-1 arrays of the knots of the spline for task=-1
629
full_output -- Non-zero to return optional outputs.
630
nxest, nyest -- Over-estimates of the total number of knots.
631
If None then nxest = max(kx+sqrt(m/2),2*kx+3),
632
nyest = max(ky+sqrt(m/2),2*ky+3)
633
quiet -- Non-zero to suppress printing of messages.
635
Outputs: (tck, {fp, ier, msg})
637
tck -- A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and
638
coefficients (c) of the bivariate B-spline representation of the
639
surface along with the degree of the spline.
641
fp -- The weighted sum of squared residuals of the spline approximation.
642
ier -- An integer flag about splrep success. Success is indicated if
643
ier<=0. If ier in [1,2,3] an error occurred but was not raised.
644
Otherwise an error is raised.
645
msg -- A message corresponding to the integer flag, ier.
649
SEE bisplev to evaluate the value of the B-spline given its tck
653
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
654
UnivariateSpline, BivariateSpline - an alternative wrapping
655
of the FITPACK functions
657
x,y,z=map(myasarray,[x,y,z])
658
x,y,z=map(ravel,[x,y,z]) # ensure 1-d arrays.
660
if not (m==len(y)==len(z)): raise TypeError, 'len(x)==len(y)==len(z) must hold.'
661
if w is None: w=ones(m,float)
663
if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
664
if xb is None: xb=x.min()
665
if xe is None: xe=x.max()
666
if yb is None: yb=y.min()
667
if ye is None: ye=y.max()
668
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
669
if s is None: s=m-sqrt(2*m)
670
if tx is None and task==-1: raise TypeError, 'Knots_x must be given for task=-1'
671
if tx is not None: _surfit_cache['tx']=myasarray(tx)
672
nx=len(_surfit_cache['tx'])
673
if ty is None and task==-1: raise TypeError, 'Knots_y must be given for task=-1'
674
if ty is not None: _surfit_cache['ty']=myasarray(ty)
675
ny=len(_surfit_cache['ty'])
676
if task==-1 and nx<2*kx+2:
677
raise TypeError, 'There must be at least 2*kx+2 knots_x for task=-1'
678
if task==-1 and ny<2*ky+2:
679
raise TypeError, 'There must be at least 2*ky+2 knots_x for task=-1'
680
if not ((1<=kx<=5) and (1<=ky<=5)):
681
raise TypeError, 'Given degree of the spline (kx,ky=%d,%d) is not supported. (1<=k<=5)'%(kx,ky)
682
if m<(kx+1)*(ky+1): raise TypeError, 'm>=(kx+1)(ky+1) must hold'
683
if nxest is None: nxest=kx+sqrt(m/2)
684
if nyest is None: nyest=ky+sqrt(m/2)
685
nxest,nyest=max(nxest,2*kx+3),max(nyest,2*ky+3)
687
nxest=int(kx+sqrt(3*m))
688
nyest=int(ky+sqrt(3*m))
690
_surfit_cache['tx']=myasarray(tx)
691
_surfit_cache['ty']=myasarray(ty)
692
tx,ty=_surfit_cache['tx'],_surfit_cache['ty']
693
wrk=_surfit_cache['wrk']
694
iwrk=_surfit_cache['iwrk']
695
u,v,km,ne=nxest-kx-1,nyest-ky-1,max(kx,ky)+1,max(nxest,nyest)
696
bx,by=kx*v+ky+1,ky*u+kx+1
698
if bx>by: b1,b2=by,by+u-kx
699
lwrk1=u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
701
tx,ty,c,o = _fitpack._surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,task,s,eps,
702
tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
703
_curfit_cache['tx']=tx
704
_curfit_cache['ty']=ty
705
_curfit_cache['wrk']=o['wrk']
706
ier,fp=o['ier'],o['fp']
708
if ier<=0 and not quiet:
709
print _iermess2[ier][0]
710
print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
712
ierm=min(11,max(-3,ier))
713
if ierm>0 and not full_output:
714
if ier in [1,2,3,4,5]:
715
print "Warning: "+_iermess2[ierm][0]
716
print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
720
raise _iermess2[ierm][1],_iermess2[ierm][0]
722
raise _iermess2['unknown'][1],_iermess2['unknown'][0]
725
return tck,fp,ier,_iermess2[ierm][0]
727
return tck,fp,ier,_iermess2['unknown'][0]
731
def bisplev(x,y,tck,dx=0,dy=0):
732
"""Evaluate a bivariate B-spline and its derivatives.
736
Return a rank-2 array of spline function values (or spline derivative
737
values) at points given by the cross-product of the rank-1 arrays x and y.
738
In special cases, return an array or just a float if either x or y or
743
x, y -- Rank-1 arrays specifying the domain over which to evaluate the
744
spline or its derivative.
745
tck -- A sequence of length 5 returned by bisplrep containing the knot
746
locations, the coefficients, and the degree of the spline:
748
dx, dy -- The orders of the partial derivatives in x and y respectively.
752
vals -- The B-pline or its derivative evaluated over the set formed by
753
the cross-product of x and y.
757
SEE bisprep to generate the tck representation.
760
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
761
UnivariateSpline, BivariateSpline - an alternative wrapping
762
of the FITPACK functions
765
if not (0<=dx<kx): raise ValueError,"0<=dx=%d<kx=%d must hold"%(dx,kx)
766
if not (0<=dy<ky): raise ValueError,"0<=dy=%d<ky=%d must hold"%(dy,ky)
767
x,y=map(myasarray,[x,y])
768
if (len(x.shape) != 1) or (len(y.shape) != 1):
769
raise ValueError, "First two entries should be rank-1 arrays."
770
z,ier=_fitpack._bispev(tx,ty,c,kx,ky,x,y,dx,dy)
771
if ier==10: raise ValueError,"Invalid input data"
772
if ier: raise TypeError,"An error occurred"
773
z.shape=len(x),len(y)
774
if len(z)>1: return z
775
if len(z[0])>1: return z[0]
778
def insert(x,tck,m=1,per=0):
779
"""Insert knots into a B-spline.
783
Given the knots and coefficients of a B-spline representation, create a
784
new B-spline with a knot inserted m times at point x.
785
This is a wrapper around the FORTRAN routine insert of FITPACK.
789
x (u) -- A 1-D point at which to insert a new knot(s). If tck was returned
790
from splprep, then the parameter values, u should be given.
791
tck -- A sequence of length 3 returned by splrep or splprep containg the
792
knots, coefficients, and degree of the spline.
793
m -- The number of times to insert the given knot (its multiplicity).
794
per -- If non-zero, input spline is considered periodic.
798
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
799
coefficients, and the degree of the new spline.
802
t(k+1) <= x <= t(n-k), where k is the degree of the spline.
803
In case of a periodic spline (per != 0) there must be
804
either at least k interior knots t(j) satisfying t(k+1)<t(j)<=x
805
or at least k interior knots t(j) satisfying x<=t(j)<t(n-k).
816
tt, cc_val, kk = insert(x, [t, c_vals, k], m)
820
tt, cc, ier = _fitpack._insert(per, t, c, k, x, m)
821
if ier==10: raise ValueError,"Invalid input data"
822
if ier: raise TypeError,"An error occurred"
825
if __name__ == "__main__":
828
if len(sys.argv[1:])>0:
829
runtest=map(string.atoi,sys.argv[1:])
832
return dot(transpose(x),x)
834
if d is None: return "sin"
835
if x is None: return "sin(x)"
836
if d%4 == 0: return sin(x)
837
if d%4 == 1: return cos(x)
838
if d%4 == 2: return -sin(x)
839
if d%4 == 3: return -cos(x)
840
def f2(x,y=0,dx=0,dy=0):
841
if x is None: return "sin(x+y)"
843
if d%4 == 0: return sin(x+y)
844
if d%4 == 1: return cos(x+y)
845
if d%4 == 2: return -sin(x+y)
846
if d%4 == 3: return -cos(x+y)
847
def test1(f=f1,per=0,s=0,a=0,b=2*pi,N=20,at=0,xb=None,xe=None):
850
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
851
x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
855
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
860
nd.append(norm2(f(t,d)-splev(t,tck,d)))
862
print "\nf = %s s=S_k(x;t,c) x in [%s, %s] > [%s, %s]"%(f(None),
863
`round(xb,3)`,`round(xe,3)`,
864
`round(a,3)`,`round(b,3)`)
865
if at: str="at knots"
866
else: str="at the middle of nodes"
867
print " per=%d s=%s Evaluation %s"%(per,`s`,str)
868
print " k : |f-s|^2 |f'-s'| |f''-.. |f'''-. |f''''- |f'''''"
876
def test2(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
877
ia=0,ib=2*pi,dx=0.2*pi):
880
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
884
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
885
nk.append([splint(ia,ib,tck),spalde(dx,tck)])
886
print "\nf = %s s=S_k(x;t,c) x in [%s, %s] > [%s, %s]"%(f(None),
887
`round(xb,3)`,`round(xe,3)`,
888
`round(a,3)`,`round(b,3)`)
889
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)`)
890
print " k : int(s,[a,b]) Int.Error Rel. error of s^(d)(dx) d = 0, .., k"
895
put(" %d %s%.8f %.1e "%(k,sr,abs(r[0]),
896
abs(r[0]-(f(ib,-1)-f(ia,-1)))))
899
put(" %.1e "%(abs(1-dr/f(dx,d))))
903
def test3(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
904
ia=0,ib=2*pi,dx=0.2*pi):
907
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
910
print " k : Roots of s(x) approx %s x in [%s,%s]:"%\
911
(f(None),`round(a,3)`,`round(b,3)`)
913
tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
914
print ' %d : %s'%(k,`sproot(tck).tolist()`)
915
def test4(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
916
ia=0,ib=2*pi,dx=0.2*pi):
919
x=a+(b-a)*arange(N+1,dtype=float)/float(N) # nodes
920
x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
923
print " u = %s N = %d"%(`round(dx,3)`,N)
924
print " k : [x(u), %s(x(u))] Error of splprep Error of splrep "%(f(0,None))
926
tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
927
tck=splrep(x,v,s=s,per=per,k=k)
929
print " %d : %s %.1e %.1e"%\
930
(k,`map(lambda x:round(x,3),uv)`,
932
abs(splev(uv[0],tck)-f(uv[0])))
933
print "Derivatives of parametric cubic spline at u (first function):"
935
tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
936
for d in range(1,k+1):
938
put(" %s "%(`uv[0]`))
941
x,y=map(myasarray,[x,y])
942
xy=array(map(lambda x,y:map(None,len(y)*[x],y),x,len(x)*[y]))
944
xy.shape=sh[0]*sh[1],sh[2]
946
def test5(f=f2,kx=3,ky=3,xb=0,xe=2*pi,yb=0,ye=2*pi,Nx=20,Ny=20,s=0):
947
x=xb+(xe-xb)*arange(Nx+1,dtype=float)/float(Nx)
948
y=yb+(ye-yb)*arange(Ny+1,dtype=float)/float(Ny)
950
tck=bisplrep(xy[0],xy[1],f(xy[0],xy[1]),s=s,kx=kx,ky=ky)
951
tt=[tck[0][kx:-kx],tck[1][ky:-ky]]
952
t2=makepairs(tt[0],tt[1])
953
v1=bisplev(tt[0],tt[1],tck)
955
v2.shape=len(tt[0]),len(tt[1])
956
print norm2(ravel(v1-v2))
959
******************************************
960
\tTests of splrep and splev
961
******************************************"""
968
test1(b=1.5*pi,xe=2*pi,per=1,s=1e-1)
971
******************************************
972
\tTests of splint and spalde
973
******************************************"""
976
test2(ia=0.2*pi,ib=pi)
977
test2(ia=0.2*pi,ib=pi,N=50)
980
******************************************
982
******************************************"""
984
print "Note that if k is not 3, some roots are missed or incorrect"
987
******************************************
988
\tTests of splprep, splrep, and splev
989
******************************************"""
994
******************************************
995
\tTests of bisplrep, bisplev
996
******************************************"""
added
removed
scipy/interpolate/fitpack.py
