[Scipy-svn] r2555 - in trunk/Lib/sandbox/spline: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jan 14 19:24:06 EST 2007
Author: jtravs
Date: 2007-01-14 18:23:52 -0600 (Sun, 14 Jan 2007)
New Revision: 2555
Added:
trunk/Lib/sandbox/spline/tests/dierckx_test_data.py
trunk/Lib/sandbox/spline/tests/test_fitpack.py
trunk/Lib/sandbox/spline/tests/test_spline.py
Removed:
trunk/Lib/sandbox/spline/tests/demos_xplt.py
trunk/Lib/sandbox/spline/tests/test_fitpack.py
trunk/Lib/sandbox/spline/tests/test_interpolate.py
Modified:
trunk/Lib/sandbox/spline/README.txt
trunk/Lib/sandbox/spline/fitpack.py
trunk/Lib/sandbox/spline/fitpack.pyf
trunk/Lib/sandbox/spline/info.py
trunk/Lib/sandbox/spline/spline.py
Log:
Added unit test suite. Fully converted to f2py.
Modified: trunk/Lib/sandbox/spline/README.txt
===================================================================
--- trunk/Lib/sandbox/spline/README.txt 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/README.txt 2007-01-15 00:23:52 UTC (rev 2555)
@@ -1,16 +1,19 @@
This module is meant to replace the spline functionality of scipy.interpolate.
-At the moment the code base is identical with a few minor cosmetic changes.
-It does not compile yet!
+Changes so far:
+ 1. removed c wrappers - only usig f2py interface now
+ 2. several bugs fixed (interface to percur, size of c array in tck)
+ 3. added 11 unit tests - now fully covered
+ 3. basic cleanup of the module files and docs
+
Planned changes are:
+ 1. possibly rationalise the interface to spline
+ 2. add further functions from fitpack fortran routines
+ 3. broaden documentation
+ 4. more tests
-1. remove dependence of f2c generated interface (use only f2py)
-2. cleanup!
-3. rationalise the interface (particularly the new spline/fitpack2 interface)
-4. build comprehensive unit test suite
-
Hopefully this module will end up in scipy, with scipy.interpolate only then
containing the interpolation functions interp1d, interp2d etc. which may depend
on the spline module and the new delaunay module (for scattered data).
-John Travers
+John Travers
\ No newline at end of file
Modified: trunk/Lib/sandbox/spline/fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/fitpack.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -18,7 +18,7 @@
NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK.
Pearu Peterson
-Modified by John Travers <jtravs at gmail.com>
+Modified by John Travers January 2007
Running test programs:
$ python fitpack.py 1 3 # run test programs 1, and 3
@@ -32,39 +32,44 @@
__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
'bisplrep', 'bisplev']
__version__ = "$Revision$"[10:-1]
-import _fitpack
+
from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
- dot, sin, cos, pi, arange, empty, int32
+ dot, sin, cos, pi, arange, empty, int32, where
myasarray = atleast_1d
-# Try to replace _fitpack interface with
-# f2py-generated version
+# f2py-generated interface to fitpack
import dfitpack
+# new class based interface to fitpack
+import spline
_iermess = {0:["""\
- The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
- -1:["""\
+ The spline has a residual sum of squares fp such that abs(fp-s)/
+ s<=0.001""",None],
+ -1:["""\
The spline is an interpolating spline (fp=0)""",None],
- -2:["""\
+ -2:["""\
The spline is weighted least-squares polynomial of degree k.
fp gives the upper bound fp0 for the smoothing factor s""",None],
- 1:["""\
+ 1:["""\
The required storage space exceeds the available storage space.
- Probable causes: data (x,y) size is too small or smoothing parameter s is too small (fp>s).""",ValueError],
- 2:["""\
+ Probable causes: data (x,y) size is too small or smoothing parameter
+ s is too small (fp>s).""",ValueError],
+ 2:["""\
A theoretically impossible results when finding a smoothin spline
- with fp = s. Probably causes: s too small. (abs(fp-s)/s>0.001)""",ValueError],
- 3:["""\
+ with fp = s. Probably causes: s too small.
+ (abs(fp-s)/s>0.001)""",ValueError],
+ 3:["""\
The maximal number of iterations (20) allowed for finding smoothing
spline with fp=s has been reached. Probably causes: s too small.
(abs(fp-s)/s>0.001)""",ValueError],
- 10:["""\
+ 10:["""\
Error on input data""",ValueError],
'unknown':["""\
An error occured""",TypeError]}
_iermess2 = {0:["""\
- The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
+ The spline has a residual sum of squares fp such that
+ abs(fp-s)/s<=0.001""",None],
-1:["""\
The spline is an interpolating spline (fp=0)""",None],
-2:["""\
@@ -75,7 +80,8 @@
norm least-squares solution of a rank deficient system.""",None],
1:["""\
The required storage space exceeds the available storage space.
- Probably causes: nxest or nyest too small or s is too small. (fp>s)""",ValueError],
+ Probably causes: nxest or nyest too small or s is too small.
+ (fp>s)""",ValueError],
2:["""\
A theoretically impossible results when finding a smoothin spline
with fp = s. Probably causes: s too small or badly chosen eps.
@@ -101,21 +107,19 @@
'unknown':["""\
An error occured""",TypeError]}
-_parcur_cache = {'t': array([],float), 'wrk': array([],float),
- 'iwrk':array([],int32), 'u': array([],float),'ub':0,'ue':1}
+_parcur_cache = {}
def splprep(x,w=None,u=None,ub=None,ue=None,k=3,task=0,s=None,t=None,
full_output=0,nest=None,per=0,quiet=1):
"""Find the B-spline representation of an N-dimensional curve.
Description:
-
- Given a list of N rank-1 arrays, x, which represent a curve in N-dimensional
- space parametrized by u, find a smooth approximating spline curve g(u).
+ Given a list of N rank-1 arrays, x, which represent a curve in
+ N-dimensional space parametrized by u, find a smooth approximating
+ spline curve g(u).
Uses the FORTRAN routine parcur from FITPACK
Inputs:
-
x -- A list of sample vector arrays representing the curve.
u -- An array of parameter values. If not given, these values are
calculated automatically as (M = len(x[0])):
@@ -124,15 +128,15 @@
u[i] = v[i] / v[M-1]
ub, ue -- The end-points of the parameters interval. Defaults to
u[0] and u[-1].
- k -- Degree of the spline. Cubic splines are recommended. Even values of
- k should be avoided especially with a small s-value.
+ k -- Degree of the spline. Cubic splines are recommended. Even values
+ of k should be avoided especially with a small s-value.
1 <= k <= 5.
task -- If task==0 find t and c for a given smoothing factor, s.
- If task==1 find t and c for another value of the smoothing factor,
- s. There must have been a previous call with task=0 or task=1
- for the same set of data.
- If task=-1 find the weighted least square spline for a given set of
- knots, t.
+ If task==1 find t and c for another value of the smoothing
+ factor, s. There must have been a previous call with task=0
+ or task=1 for the same set of data.
+ If task=-1 find the weighted least square spline for a given
+ set of knots, t.
s -- A smoothing condition. The amount of smoothness is determined by
satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
g(x) is the smoothed interpolation of (x,y). The user can use s to
@@ -148,12 +152,12 @@
help in determining the storage space. By default nest=m/2.
Always large enough is nest=m+k+1.
per -- If non-zero, data points are considered periodic with period
- x[m-1] - x[0] and a smooth periodic spline approximation is returned.
+ x[m-1] - x[0] and a smooth periodic spline approximation is
+ returned.
Values of y[m-1] and w[m-1] are not used.
quiet -- Non-zero to suppress messages.
Outputs: (tck, u, {fp, ier, msg})
-
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
coefficients, and the degree of the spline.
u -- An array of the values of the parameter.
@@ -165,19 +169,14 @@
msg -- A message corresponding to the integer flag, ier.
Remarks:
-
SEE splev for evaluation of the spline and its derivatives.
See also:
splrep, splev, sproot, spalde, splint - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
- if task<=0:
- _parcur_cache = {'t': array([],float), 'wrk': array([],float),
- 'iwrk':array([],int32),'u': array([],float),
- 'ub':0,'ue':1}
x=myasarray(x)
idim,m=x.shape
if per:
@@ -190,48 +189,56 @@
else: w=myasarray(w)
ipar=(u is not None)
if ipar:
- _parcur_cache['u']=u
- if ub is None: _parcur_cache['ub']=u[0]
- else: _parcur_cache['ub']=ub
- if ue is None: _parcur_cache['ue']=u[-1]
- else: _parcur_cache['ue']=ue
- else: _parcur_cache['u']=zeros(m,float)
+ if ub is None: ub=u[0]
+ if ue is None: ue=u[-1]
+ else: u=zeros(m,float); ub=ue=0.0
if not (1<=k<=5): raise TypeError, '1<=k=%d<=5 must hold'%(k)
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
if (not len(w)==m) or (ipar==1 and (not len(u)==m)):
raise TypeError,'Mismatch of input dimensions'
if s is None: s=m-sqrt(2*m)
- if t is None and task==-1: raise TypeError, 'Knots must be given for task=-1'
- if t is not None:
- _parcur_cache['t']=myasarray(t)
- n=len(_parcur_cache['t'])
- if task==-1 and n<2*k+2:
- raise TypeError, 'There must be at least 2*k+2 knots for task=-1'
+ if task==-1:
+ if t is None: raise TypeError, 'Knots must be given for task=-1'
+ t=myasarray(t)
+ n=len(t)
+ if n<2*k+2:
+ raise TypeError, 'There must be at least 2*k+2 knots for task=-1'
if m<=k: raise TypeError, 'm>k must hold'
if nest is None: nest=m+2*k
-
if (task>=0 and s==0) or (nest<0):
if per: nest=m+2*k
else: nest=m+k+1
nest=max(nest,2*k+3)
- u=_parcur_cache['u']
- ub=_parcur_cache['ub']
- ue=_parcur_cache['ue']
- t=_parcur_cache['t']
- wrk=_parcur_cache['wrk']
- iwrk=_parcur_cache['iwrk']
- t,c,o=_fitpack._parcur(ravel(transpose(x)),w,u,ub,ue,k,task,ipar,s,t,
- nest,wrk,iwrk,per)
- _parcur_cache['u']=o['u']
- _parcur_cache['ub']=o['ub']
- _parcur_cache['ue']=o['ue']
- _parcur_cache['t']=t
- _parcur_cache['wrk']=o['wrk']
- _parcur_cache['iwrk']=o['iwrk']
- ier,fp,n=o['ier'],o['fp'],len(t)
- u=o['u']
- c.shape=idim,n-k-1
- tcku = [t,list(c),k],u
+ if task==0:
+ u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth0(ipar,idim,u,x,
+ w,ub,ue,nest,k=k,s=s)
+ if task==1:
+ try:
+ u=_parcur_cache['u']
+ ub=_parcur_cache['ub']
+ ue=_parcur_cache['ue']
+ t=_parcur_cache['t']
+ wrk=_parcur_cache['wrk']
+ iwrk=_parcur_cache['iwrk']
+ n=_parcur_cache['n']
+ except KeyError:
+ raise ValueError, 'task=1 can only be called after task=0'
+ u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth1(ipar,idim,u,x,w,
+ ub,ue,nest,n,t,wrk,iwrk,k=k,s=s)
+ if task==-1:
+ u,ub,ue,n,t,c,fp,ier=dfitpack.parcur_lsq(ipar,idim,u,x,w,ub,ue,
+ nest,n,t,k=k)
+ if task>=0:
+ _parcur_cache['n']=n
+ _parcur_cache['u']=u
+ _parcur_cache['ub']=ub
+ _parcur_cache['ue']=ue
+ _parcur_cache['t']=t
+ _parcur_cache['wrk']=wrk
+ _parcur_cache['iwrk']=iwrk
+ c.shape=idim,nest
+ c = c[:][:n-k-1]
+ tcku = [t[:n],list(c),k],u
if ier<=0 and not quiet:
print _iermess[ier][0]
print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
@@ -251,21 +258,20 @@
else:
return tcku
-_curfit_cache = {'t': array([],float), 'wrk': array([],float),
- 'iwrk':array([],int32)}
+_splrep_cache = {}
+_percur_cache = {}
+
def splrep(x,y,w=None,xb=None,xe=None,k=3,task=0,s=1e-3,t=None,
full_output=0,per=0,quiet=1):
"""Find the B-spline representation of 1-D curve.
Description:
-
Given the set of data points (x[i], y[i]) determine a smooth spline
approximation of degree k on the interval xb <= x <= xe. The coefficients,
c, and the knot points, t, are returned. Uses the FORTRAN routine
curfit from FITPACK.
Inputs:
-
x, y -- The data points defining a curve y = f(x).
w -- Strictly positive rank-1 array of weights the same length as x and y.
The weights are used in computing the weighted least-squares spline
@@ -303,7 +309,6 @@
quiet -- Non-zero to suppress messages.
Outputs: (tck, {fp, ier, msg})
-
tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
coefficients, and the degree of the spline.
@@ -314,36 +319,35 @@
msg -- A message corresponding to the integer flag, ier.
Remarks:
+ See splev for evaluation of the spline and its derivatives.
- See splev for evaluation of the spline and its derivatives.
-
Example:
-
x = linspace(0, 10, 10)
y = sin(x)
tck = splrep(x, y)
x2 = linspace(0, 10, 200)
y2 = splev(x2, tck)
plot(x, y, 'o', x2, y2)
-
+
See also:
splprep, splev, sproot, spalde, splint - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
- if task<=0:
- _curfit_cache = {}
x,y=map(myasarray,[x,y])
m=len(x)
if w is None: w=ones(m,float)
else: w=myasarray(w)
- if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
+ if not len(w) == m:
+ raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
if (m != len(y)) or (m != len(w)):
- raise TypeError, 'Lengths of the first three arguments (x,y,w) must be equal'
+ raise TypeError, 'Lengths of the first three arguments'\
+ ' (x,y,w) must be equal'
if not (1<=k<=5):
- raise TypeError, 'Given degree of the spline (k=%d) is not supported. (1<=k<=5)'%(k)
- if m<=k: raise TypeError, 'm>k must hold'
+ raise TypeError, 'Given degree of the spline (k=%d) is not supported.'\
+ ' (1<=k<=5)'%(k)
+ if m<=k: raise TypeError, 'm>k must hold'
if xb is None: xb=x[0]
if xe is None: xe=x[-1]
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
@@ -352,31 +356,41 @@
task = -1
if task == -1:
if t is None: raise TypeError, 'Knots must be given for task=-1'
- numknots = len(t)
- _curfit_cache['t'] = empty((numknots + 2*k+2,),float)
- _curfit_cache['t'][k+1:-k-1] = t
- nest = len(_curfit_cache['t'])
- elif task == 0:
- if per:
- nest = max(m+2*k,2*k+3)
- else:
- nest = max(m+k+1,2*k+3)
- t = empty((nest,),float)
- _curfit_cache['t'] = t
- if task <= 0:
- _curfit_cache['wrk'] = empty((m*(k+1)+nest*(7+3*k),),float)
- _curfit_cache['iwrk'] = empty((nest,),int32)
- try:
- t=_curfit_cache['t']
- wrk=_curfit_cache['wrk']
- iwrk=_curfit_cache['iwrk']
- except KeyError:
- raise TypeError, "must call with task=1 only after"\
- " call with task=0,-1"
- if not per:
- n,c,fp,ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk, xb, xe, k, s)
- else:
- n,c,fp,ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s)
+ if per:
+ if task == 0:
+ n,t,c,fp,wrk,iwrk,ier = dfitpack.percur_smth0(x,y,w,k=k,s=s)
+ elif task ==1:
+ try:
+ t=_percur_cache['t']
+ wrk=_percur_cache['wrk']
+ iwrk=_percur_cache['iwrk']
+ n=_percur_cache['n']
+ except KeyError:
+ raise ValueError, 'task=1 can only be called after task=0'
+ n,t,c,fp,wrk,iwrk,ier = dfitpack.percur_smth1(x,y,w,n,t,wrk,iwrk,
+ k=k,s=s)
+ elif task==-1:
+ n,t,c,fp,ier = dfitpack.percur_lsq(x,y,w,n,t,k=k)
+ if task>=0:
+ _percur_cache['t']=t
+ _percur_cache['wrk']=wrk
+ _percur_cache['iwrk']=iwrk
+ _percur_cache['n']=n
+ if task == 0:
+ spl = spline.UnivariateSpline(x,y,w,[xb,xe],k=k,s=s)
+ elif task == 1:
+ try:
+ spl = _splrep_cache['spl']
+ except KeyError:
+ raise ValueError, 'task=1 can only be called after task=0'
+ spl.set_smoothing_factor(s)
+ elif task == -1:
+ t = t[where(t>xb)]
+ t = t[where(t<xe)]
+ spl = spline.LSQUnivariateSpline(x,y,t,w,[xb,xe],k=k)
+ if task>=0:
+ _splrep_cache['spl'] = spl
+ x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = spl._data
tck = (t[:n],c[:n-k-1],k)
if ier<=0 and not quiet:
print _iermess[ier][0]
@@ -397,22 +411,15 @@
else:
return tck
-def _ntlist(l): # return non-trivial list
- return l
- #if len(l)>1: return l
- #return l[0]
-
def splev(x,tck,der=0):
"""Evaulate a B-spline and its derivatives.
Description:
-
Given the knots and coefficients of a B-spline representation, evaluate
the value of the smoothing polynomial and it's derivatives.
This is a wrapper around the FORTRAN routines splev and splder of FITPACK.
Inputs:
-
x (u) -- a 1-D array of points at which to return the value of the
smoothed spline or its derivatives. If tck was returned from
splprep, then the parameter values, u should be given.
@@ -422,7 +429,6 @@
or equal to k).
Outputs: (y, )
-
y -- an array of values representing the spline function or curve.
If tck was returned from splrep, then this is a list of arrays
representing the curve in N-dimensional space.
@@ -430,54 +436,59 @@
See also:
splprep, splrep, sproot, spalde, splint - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
t,c,k=tck
try:
- c[0][0]
+ c[0][0] # check if c is >1-d
return map(lambda c,x=x,t=t,k=k,der=der:splev(x,[t,c,k],der),c)
- except: pass
+ except TypeError: pass # c is 1-d
if not (0<=der<=k):
raise ValueError,"0<=der=%d<=k=%d must hold"%(der,k)
x=myasarray(x)
- y,ier=_fitpack._spl_(x,der,t,c,k)
+ c = c.copy()
+ c.resize(len(t))
+ if (der>0):
+ y,ier=dfitpack.splder(t,c,k,x,der)
+ else:
+ y,ier=dfitpack.splev(t,c,k,x)
if ier==10: raise ValueError,"Invalid input data"
- if ier: raise TypeError,"An error occurred"
+ if ier: raise TypeError,"An unkown error occurred"
if len(y)>1: return y
return y[0]
-def splint(a,b,tck,full_output=0):
+def splint(a,b,tck,full_output=False):
"""Evaluate the definite integral of a B-spline.
Description:
-
Given the knots and coefficients of a B-spline, evaluate the definite
integral of the smoothing polynomial between two given points.
Inputs:
-
a, b -- The end-points of the integration interval.
tck -- A length 3 sequence describing the given spline (See splev).
- full_output -- Non-zero to return optional output.
+ full_output -- True to return optional output.
Outputs: (integral, {wrk})
-
integral -- The resulting integral.
- wrk -- An array containing the integrals of the normalized B-splines defined
- on the set of knots.
+ wrk -- An array containing the integrals of the normalized B-splines
+ defined on the set of knots.
-
See also:
splprep, splrep, sproot, spalde, splev - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
t,c,k=tck
- try: c[0][0];return _ntlist(map(lambda c,a=a,b=b,t=t,k=k:splint(a,b,[t,c,k]),c))
- except: pass
- aint,wrk=_fitpack._splint(t,c,k,a,b)
+ try:
+ c[0][0] # check if c is >1-d
+ return map(lambda c,a=a,b=b,t=t,k=k:splint(a,b,[t,c,k]),c)
+ except TypeError: pass # c is 1-d
+ c = c.copy()
+ c.resize(len(t))
+ aint,wrk = dfitpack.splint(t,c,k,a,b)
if full_output: return aint,wrk
else: return aint
@@ -485,38 +496,40 @@
"""Find the roots of a cubic B-spline.
Description:
-
Given the knots (>=8) and coefficients of a cubic B-spline return the
roots of the spline.
Inputs:
-
tck -- A length 3 sequence describing the given spline (See splev).
- The number of knots must be >= 8. The knots must be a montonically
- increasing sequence.
+ The number of knots must be >= 8. The knots must be a
+ montonically increasing sequence.
mest -- An estimate of the number of zeros (Default is 10).
Outputs: (zeros, )
-
zeros -- An array giving the roots of the spline.
See also:
splprep, splrep, splint, spalde, splev - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
t,c,k=tck
- if k==4: t=t[1:-1]
- if k==5: t=t[2:-2]
- try: c[0][0];return _ntlist(map(lambda c,t=t,k=k,mest=mest:sproot([t,c,k],mest),c))
- except: pass
+ if k != 3:
+ raise ValueError, "Sproot only valid for cubic bsplines"
+ try:
+ c[0][0] # check if c is >1-d
+ return map(lambda c,t=t,k=k,mest=mest:sproot([t,c,k],mest),c)
+ except TypeError: pass # c is 1-d
if len(t)<8:
- raise TypeError,"The number of knots %d>=8"%(len(t))
- z,ier=_fitpack._sproot(t,c,k,mest)
+ raise TypeError,"The number of knots be >=8"
+ c = c.copy()
+ c.resize(len(t))
+ z,m,ier=dfitpack.sproot(t,c,mest)
if ier==10:
- raise TypeError,"Invalid input data. t1<=..<=t4<t5<..<tn-3<=..<=tn must hold."
- if ier==0: return z
+ raise TypeError,"Invalid input data. t1<=..<=t4<t5<..<tn-3<=..<=tn" \
+ "must hold."
+ if ier==0: return z[0:m]
if ier==1:
print "Warning: the number of zeros exceeds mest"
return z
@@ -526,61 +539,56 @@
"""Evaluate all derivatives of a B-spline.
Description:
-
Given the knots and coefficients of a cubic B-spline compute all
derivatives up to order k at a point (or set of points).
Inputs:
-
tck -- A length 3 sequence describing the given spline (See splev).
x -- A point or a set of points at which to evaluate the derivatives.
Note that t(k) <= x <= t(n-k+1) must hold for each x.
Outputs: (results, )
-
results -- An array (or a list of arrays) containing all derivatives
up to order k inclusive for each point x.
See also:
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
bisplrep, bisplev - bivariate splines
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
t,c,k=tck
try:
- c[0][0]
- return _ntlist(map(lambda c,x=x,t=t,k=k:spalde(x,[t,c,k]),c))
- except: pass
+ c[0][0] # check if c is >1-d
+ return map(lambda c,x=x,t=t,k=k:spalde(x,[t,c,k]),c)
+ except TypeError: pass # c is 1-d
try: x=x.tolist()
except:
try: x=list(x)
except: x=[x]
if len(x)>1:
return map(lambda x,tck=tck:spalde(x,tck),x)
- d,ier=_fitpack._spalde(t,c,k,x[0])
+ c = c.copy()
+ c.resize(len(t))
+ d,ier=dfitpack.spalde(t,c,k,x[0])
if ier==0: return d
if ier==10:
raise TypeError,"Invalid input data. t(k)<=x<=t(n-k+1) must hold."
raise TypeError,"Unknown error"
-#def _curfit(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
-# full_output=0,nest=None,per=0,quiet=1):
-_surfit_cache = {'tx': array([],float),'ty': array([],float),
- 'wrk': array([],float), 'iwrk':array([],int32)}
+_surfit_cache = {}
+
def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,
s=None,eps=1e-16,tx=None,ty=None,full_output=0,
nxest=None,nyest=None,quiet=1):
"""Find a bivariate B-spline representation of a surface.
Description:
-
Given a set of data points (x[i], y[i], z[i]) representing a surface
z=f(x,y), compute a B-spline representation of the surface.
Inputs:
-
x, y, z -- Rank-1 arrays of data points.
w -- Rank-1 array of weights. By default w=ones(len(x)).
xb, xe -- End points of approximation interval in x.
@@ -609,7 +617,6 @@
quiet -- Non-zero to suppress printing of messages.
Outputs: (tck, {fp, ier, msg})
-
tck -- A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and
coefficients (c) of the bivariate B-spline representation of the
surface along with the degree of the spline.
@@ -621,40 +628,41 @@
msg -- A message corresponding to the integer flag, ier.
Remarks:
-
SEE bisplev to evaluate the value of the B-spline given its tck
representation.
See also:
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
x,y,z=map(myasarray,[x,y,z])
x,y,z=map(ravel,[x,y,z]) # ensure 1-d arrays.
m=len(x)
- if not (m==len(y)==len(z)): raise TypeError, 'len(x)==len(y)==len(z) must hold.'
+ if not (m==len(y)==len(z)):
+ raise TypeError, 'len(x)==len(y)==len(z) must hold.'
if w is None: w=ones(m,float)
else: w=myasarray(w)
- if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
+ if not len(w) == m:
+ raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
if xb is None: xb=x.min()
if xe is None: xe=x.max()
if yb is None: yb=y.min()
if ye is None: ye=y.max()
if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
if s is None: s=m-sqrt(2*m)
- if tx is None and task==-1: raise TypeError, 'Knots_x must be given for task=-1'
- if tx is not None: _surfit_cache['tx']=myasarray(tx)
- nx=len(_surfit_cache['tx'])
- if ty is None and task==-1: raise TypeError, 'Knots_y must be given for task=-1'
- if ty is not None: _surfit_cache['ty']=myasarray(ty)
- ny=len(_surfit_cache['ty'])
- if task==-1 and nx<2*kx+2:
- raise TypeError, 'There must be at least 2*kx+2 knots_x for task=-1'
- if task==-1 and ny<2*ky+2:
- raise TypeError, 'There must be at least 2*ky+2 knots_x for task=-1'
+ if task==-1:
+ if tx is None: raise TypeError, 'Knots_x must be given for task=-1'
+ _surfit_cache['tx']=myasarray(tx)
+ if ty is None: raise TypeError, 'Knots_y must be given for task=-1'
+ _surfit_cache['ty']=myasarray(ty)
+ if nx<2*kx+2:
+ raise TypeError,'There must be at least 2*kx+2 knots_x for task=-1'
+ if ny<2*ky+2:
+ raise TypeError,'There must be at least 2*ky+2 knots_y for task=-1'
if not ((1<=kx<=5) and (1<=ky<=5)):
- raise TypeError, 'Given degree of the spline (kx,ky=%d,%d) is not supported. (1<=k<=5)'%(kx,ky)
+ raise TypeError, 'Given degree of the spline (kx,ky=%d,%d) is not \
+ supported. (1<=k<=5)'%(kx,ky)
if m<(kx+1)*(ky+1): raise TypeError, 'm>=(kx+1)(ky+1) must hold'
if nxest is None: nxest=kx+sqrt(m/2)
if nyest is None: nyest=ky+sqrt(m/2)
@@ -662,25 +670,39 @@
if task>=0 and s==0:
nxest=int(kx+sqrt(3*m))
nyest=int(ky+sqrt(3*m))
- if task==-1:
- _surfit_cache['tx']=myasarray(tx)
- _surfit_cache['ty']=myasarray(ty)
- tx,ty=_surfit_cache['tx'],_surfit_cache['ty']
- wrk=_surfit_cache['wrk']
- iwrk=_surfit_cache['iwrk']
- u,v,km,ne=nxest-kx-1,nyest-ky-1,max(kx,ky)+1,max(nxest,nyest)
- bx,by=kx*v+ky+1,ky*u+kx+1
- b1,b2=bx,bx+v-ky
- if bx>by: b1,b2=by,by+u-kx
- lwrk1=u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
- lwrk2=u*v*(b2+1)+b2
- tx,ty,c,o = _fitpack._surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,task,s,eps,
- tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
- _curfit_cache['tx']=tx
- _curfit_cache['ty']=ty
- _curfit_cache['wrk']=o['wrk']
- ier,fp=o['ier'],o['fp']
- tck=[tx,ty,c,kx,ky]
+ if task == 0:
+ nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth0(x,y,z,w,xb,xe,yb,ye,
+ kx,ky,s=s,eps=eps,
+ nxest=nxest,nyest=nyest)
+ elif task==1:
+ try:
+ tx,ty=_surfit_cache['tx'],_surfit_cache['ty']
+ nx,ny= _surfit_cache['nx'],_surfit_cache['ny']
+ wrk1 = _surfit_cache['wrk1']
+ except KeyError:
+ raise ValueError, 'task=1 can only be called after task=0'
+ nx,tx,ny,ty,c,fp,wrk1,ier=dfitpack.surfit_smth1(x,y,z,nx,tx,ny,ty,wrk1,
+ w,xb,xe,yb,ye,kx,ky,s=s,eps=eps,
+ nxest=nxest,nyest=nyest)
+ elif task==-1:
+ tx = myasarray(tx)
+ ty = myasarray(ty)
+ tx,ty,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx,ty,w,
+ xb,xe,yb,ye,
+ kx,ky,eps,
+ nxest=nxest,nyest=nyest)
+ if ier>10:
+ tx,ty,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx,ty,w,\
+ xb,xe,yb,ye,\
+ kx,ky,eps,\
+ nxest=nxest,nyest=nyest)
+ if task>=0:
+ _surfit_cache['tx']=tx
+ _surfit_cache['ty']=ty
+ _surfit_cache['nx']=nx
+ _surfit_cache['ny']=ny
+ _surfit_cache['wrk1']=wrk1
+ tck=[tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)],kx,ky]
if ier<=0 and not quiet:
print _iermess2[ier][0]
print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
@@ -708,14 +730,12 @@
"""Evaluate a bivariate B-spline and its derivatives.
Description:
-
Return a rank-2 array of spline function values (or spline derivative
values) at points given by the cross-product of the rank-1 arrays x and y.
In special cases, return an array or just a float if either x or y or
both are floats.
Inputs:
-
x, y -- Rank-1 arrays specifying the domain over which to evaluate the
spline or its derivative.
tck -- A sequence of length 5 returned by bisplrep containing the knot
@@ -724,17 +744,15 @@
dx, dy -- The orders of the partial derivatives in x and y respectively.
Outputs: (vals, )
-
vals -- The B-pline or its derivative evaluated over the set formed by
the cross-product of x and y.
Remarks:
-
SEE bisprep to generate the tck representation.
See also:
splprep, splrep, splint, sproot, splev - evaluation, roots, integral
- UnivariateSpline, BivariateSpline - an alternative wrapping
+ UnivariateSpline, BivariateSpline - an alternative wrapping
of the FITPACK functions
"""
tx,ty,c,kx,ky=tck
@@ -743,7 +761,10 @@
x,y=map(myasarray,[x,y])
if (len(x.shape) != 1) or (len(y.shape) != 1):
raise ValueError, "First two entries should be rank-1 arrays."
- z,ier=_fitpack._bispev(tx,ty,c,kx,ky,x,y,dx,dy)
+ if (dx>0 or dy>0):
+ z,ier=dfitpack.parder(tx,ty,c,kx,ky,dx,dy,x,y)
+ else:
+ z,ier=dfitpack.bispev(tx,ty,c,kx,ky,x,y)
if ier==10: raise ValueError,"Invalid input data"
if ier: raise TypeError,"An error occurred"
z.shape=len(x),len(y)
@@ -751,7 +772,6 @@
if len(z[0])>1: return z[0]
return z[0][0]
-
if __name__ == "__main__":
import sys,string
runtest=range(10)
@@ -782,7 +802,7 @@
v,v1=f(x),f(x1)
nk=[]
for k in range(1,6):
- tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
+ tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
if at:t=tck[0][k:-k]
else: t=x1
nd=[]
@@ -811,7 +831,7 @@
v=f(x)
nk=[]
for k in range(1,6):
- tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
+ tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
nk.append([splint(ia,ib,tck),spalde(dx,tck)])
print "\nf = %s s=S_k(x;t,c) x in [%s, %s] > [%s, %s]"%(f(None),
`round(xb,3)`,`round(xe,3)`,
@@ -839,9 +859,9 @@
nk=[]
print " k : Roots of s(x) approx %s x in [%s,%s]:"%\
(f(None),`round(a,3)`,`round(b,3)`)
- for k in range(1,6):
- tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
- print ' %d : %s'%(k,`sproot(tck).tolist()`)
+ k=3
+ tck=splrep(x,v,s=s,per=per,k=k,xe=xe)
+ print ' %d : %s'%(k,`sproot(tck).tolist()`)
def test4(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
ia=0,ib=2*pi,dx=0.2*pi):
if xb is None: xb=a
@@ -853,8 +873,8 @@
print " u = %s N = %d"%(`round(dx,3)`,N)
print " k : [x(u), %s(x(u))] Error of splprep Error of splrep "%(f(0,None))
for k in range(1,6):
- tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
- tck=splrep(x,v,s=s,per=per,k=k,nest=-1)
+ tckp,u=splprep([x,v],s=s,per=per,k=k)
+ tck=splrep(x,v,s=s,per=per,k=k)
uv=splev(dx,tckp)
print " %d : %s %.1e %.1e"%\
(k,`map(lambda x:round(x,3),uv)`,
@@ -862,7 +882,7 @@
abs(splev(uv[0],tck)-f(uv[0])))
print "Derivatives of parametric cubic spline at u (first function):"
k=3
- tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
+ tckp,u=splprep([x,v],s=s,per=per,k=k)
for d in range(1,k+1):
uv=splev(dx,tckp,d)
put(" %s "%(`uv[0]`))
Modified: trunk/Lib/sandbox/spline/fitpack.pyf
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.pyf 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/fitpack.pyf 2007-01-15 00:23:52 UTC (rev 2555)
@@ -43,7 +43,6 @@
static int imax(int i1,int i2) {
return MAX(i1,i2);
}
-
static int calc_surfit_lwrk1(int m, int kx, int ky, int nxest, int nyest) {
int u = nxest-kx-1;
int v = nyest-ky-1;
@@ -64,8 +63,7 @@
int b2 = (bx<=by?bx+v-ky:by+u-kx);
return u*v*(b2+1)+b2;
}
-
-static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
+static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
int nxest, int nyest) {
int u = MAX(my, nxest);
return 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+my*(ky+1)+u;
@@ -77,7 +75,7 @@
!!!!!!!!!! Univariate spline !!!!!!!!!!!
subroutine splev(t,n,c,k,x,y,m,ier)
- ! y = splev(t,c,k,x)
+ ! y,ier = splev(t,c,k,x)
real*8 dimension(n),intent(in) :: t
integer intent(hide),depend(t) :: n=len(t)
real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
@@ -85,32 +83,32 @@
real*8 dimension(m),intent(in) :: x
real*8 dimension(m),depend(m),intent(out) :: y
integer intent(hide),depend(x) :: m=len(x)
- integer intent(hide) :: ier
+ integer intent(out) :: ier
end subroutine splev
subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
- ! dy = splder(t,c,k,x,[nu])
+ ! dy,ier = splder(t,c,k,x,[nu])
real*8 dimension(n) :: t
integer depend(t),intent(hide) :: n=len(t)
real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
integer :: k
- integer depend(k),check(0<=nu && nu<=k) :: nu = 1
+ integer optional,depend(k),check(0<=nu && nu<=k) :: nu = 1
real*8 dimension(m) :: x
real*8 dimension(m),depend(m),intent(out) :: y
integer depend(x),intent(hide) :: m=len(x)
real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
- integer intent(hide) :: ier
+ integer intent(out) :: ier
end subroutine splder
function splint(t,n,c,k,a,b,wrk)
- ! iy = splint(t,c,k,a,b)
+ ! iy,wrk = splint(t,c,k,a,b,wrk)
real*8 dimension(n),intent(in) :: t
integer intent(hide),depend(t) :: n=len(t)
real*8 dimension(n),depend(n),check(len(c)==n) :: c
integer intent(in) :: k
real*8 intent(in) :: a
real*8 intent(in) :: b
- real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
+ real*8 dimension(n),intent(out) :: wrk
real*8 :: splint
end function splint
@@ -127,10 +125,8 @@
subroutine spalde(t,n,c,k,x,d,ier)
! d,ier = spalde(t,c,k,x)
-
callprotoargument double*,int*,double*,int*,double*,double*,int*
callstatement {int k1=k+1; (*f2py_func)(t,&n,c,&k1,&x,d,&ier); }
-
real*8 dimension(n) :: t
integer intent(hide),depend(t) :: n=len(t)
real*8 dimension(n),depend(n),check(len(c)==n) :: c
@@ -140,84 +136,162 @@
integer intent(out) :: ier
end subroutine spalde
- subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in curfit.f
- integer :: iopt
+ subroutine percur_lsq(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
+ ! n,t,c,fp,ier = percur_lsq(x,y,w,n,t,[k])
+ fortranname percur
+ integer intent(hide):: iopt=-1
integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
real*8 dimension(m) :: x
real*8 dimension(m),depend(m),check(len(y)==m) :: y
real*8 dimension(m),depend(m),check(len(w)==m) :: w
- real*8 optional,depend(x),check(xb<=x[0]) :: xb = x[0]
- real*8 optional,depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
integer optional,check(1<=k && k <=5),intent(in) :: k=3
+ real*8 intent(hide),check(s>=0.0) :: s = 0.0
+ integer intent(hide),depend(m,k) :: nest=m+2*k
+ integer intent(in,out) :: n
+ real*8 dimension(n),intent(in,out) :: t
+ real*8 dimension(n),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk),intent(hide,cache) :: wrk
+ integer intent(hide),depend(m,k,nest) :: lwrk=m*(k+1)+nest*(8+5*k)
+ integer dimension(nest),intent(hide,cache) :: iwrk
+ integer intent(out) :: ier
+ end subroutine percur_lsq
+
+ subroutine percur_smth1(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
+ ! n,t,c,fp,wrk,iwrk,ier = percur_smth1(x,y,w,n,t,wrk,iwrk,[k,s])
+ fortranname percur
+ integer intent(hide):: iopt=1
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m) :: w
+ integer optional,check(1<=k && k <=5),intent(in) :: k=3
real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
+ integer intent(hide),depend(m,k) :: nest=m+2*k
+ integer intent(in,out),depend(nest) :: n=nest
+ real*8 dimension(nest),intent(in,out) :: t
real*8 dimension(n),intent(out) :: c
real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
+ real*8 dimension(lwrk),intent(in,out) :: wrk
+ integer intent(hide),depend(m,k,nest) :: lwrk=m*(k+1)+nest*(8+5*k)
+ integer dimension(nest),intent(in,out) :: iwrk
integer intent(out) :: ier
- end subroutine curfit
+ end subroutine percur_smth1
- subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
- ! in percur.f
- integer :: iopt
+ subroutine percur_smth0(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
+ ! n,t,c,fp,wrk,iwrk,ier = percur_smth0(x,y,w,[k,s])
+ fortranname percur
+ integer intent(hide):: iopt=0
integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
real*8 dimension(m) :: x
real*8 dimension(m),depend(m),check(len(y)==m) :: y
real*8 dimension(m),depend(m),check(len(w)==m) :: w
integer optional,check(1<=k && k <=5),intent(in) :: k=3
real*8 optional,check(s>=0.0) :: s = 0.0
- integer intent(hide),depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest),intent(inout) :: t
+ integer intent(hide),depend(m,k) :: nest=m+2*k
+ integer intent(out),depend(nest) :: n=nest
+ real*8 dimension(nest),intent(out) :: t
real*8 dimension(n),intent(out) :: c
real*8 intent(out) :: fp
- real*8 dimension(lwrk),intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest),intent(inout) :: iwrk
+ real*8 dimension(lwrk),intent(out) :: wrk
+ integer intent(hide),depend(m,k,nest) :: lwrk=m*(k+1)+nest*(8+5*k)
+ integer dimension(nest),intent(out) :: iwrk
integer intent(out) :: ier
- end subroutine percur
-
+ end subroutine percur_smth0
- subroutine parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,fp,wrk,lwrk,iwrk,ier)
- ! in parcur.f
- integer check(iopt>=-1 && iopt <= 1):: iopt
+ subroutine parcur_lsq(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,&
+ fp,wrk,lwrk,iwrk,ier)
+ !u,ub,ue,n,t,c,fp,ier=parcur_lsq(ipar,idim,u,x,w,ub,ue,nest,n,t,[k])
+ fortranname parcur
+ integer intent(hide) :: iopt = -1
integer check(ipar == 1 || ipar == 0) :: ipar
integer check(idim > 0 && idim < 11) :: idim
integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
- real*8 dimension(m), intent(inout) :: u
+ real*8 dimension(m), intent(in,out) :: u
integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
real*8 dimension(mx) :: x
real*8 dimension(m) :: w
- real*8 :: ub
- real*8 :: ue
- integer optional, check(1<=k && k<=5) :: k=3.0
- real*8 optional, check(s>=0.0) :: s = 0.0
- integer intent(hide), depend(t) :: nest=len(t)
- integer intent(out), depend(nest) :: n=nest
- real*8 dimension(nest), intent(inout) :: t
- integer intent(hide), depend(c,nest,idim), check(nc>=idim*nest) :: nc=len(c)
+ real*8 intent(in,out) :: ub
+ real*8 intent(in,out) :: ue
+ integer check(1<=k && k<=5) :: k=3.0
+ real*8 intent(hide),check(s>=0.0) :: s = 0.0
+ integer intent(in) :: nest
+ integer intent(in,out) :: n
+ real*8 dimension(nest), intent(in,out) :: t
+ integer intent(hide), depend(nest,idim) :: nc=idim*nest
real*8 dimension(nc), intent(out) :: c
real*8 intent(out) :: fp
- real*8 dimension(lwrk), intent(inout) :: wrk
- integer intent(hide),depend(wrk) :: lwrk=len(wrk)
- integer dimension(nest), intent(inout) :: iwrk
+ real*8 dimension(lwrk), intent(cache) :: wrk
+ integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
+ integer dimension(nest), intent(cache) :: iwrk
integer intent(out) :: ier
- end subroutine parcur
+ end subroutine parcur_lsq
+ subroutine parcur_smth1(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,&
+ c,fp,wrk,lwrk,iwrk,ier)
+ !u,ub,ue,n,t,c,fp,wrk,iwrk,ier=parcur_smth1(ipar,idim,u,x,w,ub,ue,nest,
+ ! n,t,wrk,iwrk,[k,s])
+ fortranname parcur
+ integer intent(hide) :: iopt = 1
+ integer check(ipar == 1 || ipar == 0) :: ipar
+ integer check(idim > 0 && idim < 11) :: idim
+ integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+ real*8 dimension(m), intent(in,out) :: u
+ integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+ real*8 dimension(mx) :: x
+ real*8 dimension(m) :: w
+ real*8 intent(in,out) :: ub
+ real*8 intent(in,out) :: ue
+ integer check(1<=k && k<=5) :: k=3.0
+ real*8 check(s>=0.0) :: s = 0.0
+ integer intent(in) :: nest
+ integer intent(in,out) :: n
+ real*8 dimension(nest), intent(in,out) :: t
+ integer intent(hide), depend(nest,idim) :: nc=idim*nest
+ real*8 dimension(nc), intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk), intent(in,out) :: wrk
+ integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
+ integer dimension(nest), intent(in,out) :: iwrk
+ integer intent(out) :: ier
+ end subroutine parcur_smth1
- subroutine fpcurf0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ subroutine parcur_smth0(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,&
+ c,fp,wrk,lwrk,iwrk,ier)
+ !u,ub,ue,n,t,c,fp,wrk,iwrk,ier=parcur_smth0(ipar,idim,u,x,w,
+ ! ub,ue,nest,[k,s])
+ fortranname parcur
+ integer intent(hide) :: iopt = 0
+ integer check(ipar == 1 || ipar == 0) :: ipar
+ integer check(idim > 0 && idim < 11) :: idim
+ integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+ real*8 dimension(m), intent(in,out) :: u
+ integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+ real*8 dimension(mx) :: x
+ real*8 dimension(m) :: w
+ real*8 intent(in,out) :: ub
+ real*8 intent(in,out) :: ue
+ integer check(1<=k && k<=5) :: k=3.0
+ real*8 check(s>=0.0) :: s = 0.0
+ integer intent(in) :: nest
+ integer intent(out) :: n
+ real*8 dimension(nest), intent(out) :: t
+ integer intent(hide), depend(nest,idim) :: nc=idim*nest
+ real*8 dimension(nc), intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk), intent(out) :: wrk
+ integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
+ integer dimension(nest), intent(out) :: iwrk
+ integer intent(out) :: ier
+ end subroutine parcur_smth0
+
+ subroutine fpcurf_smth0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,&
+ c,fp,fpint,wrk,nrdata,ier)
! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf0(x,y,k,[w,xb,xe,s,nest])
-
+ ! fpcurf_smth0(x,y,k,[w,xb,xe,s,nest])
fortranname fpcurf
callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
integer intent(hide) :: iopt = 0
real*8 dimension(m),intent(in,out) :: x
real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
@@ -240,16 +314,16 @@
real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
integer intent(out) :: ier
- end subroutine fpcurf0
+ end subroutine fpcurf_smth0
- subroutine fpcurf1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ subroutine fpcurf_smth1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,&
+ c,fp,fpint,wrk,nrdata,ier)
! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurf1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
+ ! fpcurf_smth1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
fortranname fpcurf
callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
integer intent(hide) :: iopt = 1
real*8 dimension(m),intent(in,out,overwrite) :: x
real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out,overwrite) :: y
@@ -272,16 +346,15 @@
real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
integer dimension(nest),depend(nest),check(len(nrdata)==nest),intent(in,out,cache,overwrite) :: nrdata
integer intent(in,out) :: ier
- end subroutine fpcurf1
+ end subroutine fpcurf_smth1
- subroutine fpcurfm1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ subroutine fpcurf_lsq(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,&
+ fp,fpint,wrk,nrdata,ier)
! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
- ! fpcurfm1(x,y,k,t,[w,xb,xe])
-
+ ! fpcurf_lsq(x,y,k,t,[w,xb,xe])
fortranname fpcurf
callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
integer intent(hide) :: iopt = -1
real*8 dimension(m),intent(in,out) :: x
real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
@@ -304,7 +377,7 @@
real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
integer intent(out) :: ier
- end subroutine fpcurfm1
+ end subroutine fpcurfm_lsq
!!!!!!!!!! Bivariate spline !!!!!!!!!!!
@@ -330,13 +403,81 @@
integer intent(out) :: ier
end subroutine bispev
- subroutine surfit_smth(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
+ subroutine parder(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,iwrk,kwrk,ier)
+ ! z,ier = parder(tx,ty,c,kx,ky,nux,nuy,x,y)
+ real*8 dimension(nx),intent(in) :: tx
+ integer intent(hide),depend(tx) :: nx=len(tx)
+ real*8 dimension(ny),intent(in) :: ty
+ integer intent(hide),depend(ty) :: ny=len(ty)
+ real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
+ check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
+ integer :: kx
+ integer :: ky
+ integer intent(in),depend(kx),check(0<=nux && nux<kx) :: nux
+ integer intent(in),depend(ky),check(0<=nuy && nuy<ky) :: nuy
+ real*8 intent(in),dimension(mx) :: x
+ integer intent(hide),depend(x) :: mx=len(x)
+ real*8 intent(in),dimension(my) :: y
+ integer intent(hide),depend(y) :: my=len(y)
+ real*8 dimension(mx,my),depend(mx,my),intent(out,c) :: z
+ real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
+ integer intent(hide),depend(mx,kx,my,ky) :: lwrk=mx*(kx+1)+my*(ky+1)
+ integer dimension(kwrk),depend(kwrk),intent(hide,cache) :: iwrk
+ integer intent(hide),depend(mx,my) :: kwrk=mx+my
+ integer intent(out) :: ier
+ end subroutine parder
+
+ subroutine surfit_smth1(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = surfit_smth(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
+ ! nx,tx,ny,ty,c,fp,wrk1,ier = surfit_smth1(x,y,z,nx,tx,ny,ty,wrk1,
+ ![w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
+ fortranname surfit
+ integer intent(hide) :: iopt=1
+ integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
+ :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(z)==m) :: z
+ real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
+ real*8 optional,depend(x,m) :: xb=dmin(x,m)
+ real*8 optional,depend(x,m) :: xe=dmax(x,m)
+ real*8 optional,depend(y,m) :: yb=dmin(y,m)
+ real*8 optional,depend(y,m) :: ye=dmax(y,m)
+ integer check(1<=kx && kx<=5) :: kx = 3
+ integer check(1<=ky && ky<=5) :: ky = 3
+ real*8 optional,check(0.0<=s) :: s = m
+ integer optional,depend(kx,m),check(nxest>=2*(kx+1)) &
+ :: nxest = imax(kx+1+sqrt(m/2),2*(kx+1))
+ integer optional,depend(ky,m),check(nyest>=2*(ky+1)) &
+ :: nyest = imax(ky+1+sqrt(m/2),2*(ky+1))
+ integer intent(hide),depend(nxest,nyest) :: nmax=MAX(nxest,nyest)
+ real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
+ integer intent(in,out) :: nx
+ real*8 dimension(nmax),intent(in,out),depend(nmax) :: tx
+ integer intent(in,out) :: ny
+ real*8 dimension(nmax),intent(in,out),depend(nmax) :: ty
+ real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
+ depend(kx,ky,nxest,nyest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk1),intent(in,out),depend(lwrk1) :: wrk1
+ integer intent(hide),depend(m,kx,ky,nxest,nyest) &
+ :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
+ real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
+ integer optional,intent(in),depend(kx,ky,nxest,nyest) &
+ :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
+ integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
+ integer intent(hide),depend(m,nxest,nyest,kx,ky) &
+ :: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
+ integer intent(out) :: ier
+ end subroutine surfit_smth1
+
+ subroutine surfit_smth0(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
+ nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
+ iwrk,kwrk,ier)
+ ! nx,tx,ny,ty,c,fp,ier = surfit_smth0(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
fortranname surfit
-
integer intent(hide) :: iopt=0
integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
:: m=len(x)
@@ -374,15 +515,13 @@
integer intent(hide),depend(m,nxest,nyest,kx,ky) &
:: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
integer intent(out) :: ier
- end subroutine surfit_smth
+ end subroutine surfit_smth0
subroutine surfit_lsq(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
iwrk,kwrk,ier)
! tx,ty,c,fp,ier = surfit_lsq(x,y,z,tx,ty,[w,xb,xe,yb,ye,kx,ky,eps,lwrk2])
-
fortranname surfit
-
integer intent(hide) :: iopt=-1
integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
:: m=len(x)
@@ -418,17 +557,15 @@
:: kwrk=m+(nx-2*kx-1)*(ny-2*ky-1)
integer intent(out) :: ier
end subroutine surfit_lsq
-
- subroutine regrid_smth(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,&
+
+ subroutine regrid_smth0(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,&
nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
- ! nx,tx,ny,ty,c,fp,ier = regrid_smth(x,y,z,[xb,xe,yb,ye,kx,ky,s])
-
+ ! nx,tx,ny,ty,c,fp,ier = regrid_smth0(x,y,z,[xb,xe,yb,ye,kx,ky,s])
fortranname regrid
-
integer intent(hide) :: iopt=0
integer intent(hide),depend(x,kx),check(mx>kx) :: mx=len(x)
real*8 dimension(mx) :: x
- integer intent(hide),depend(y,ky),check(my>ky) :: my=len(y)
+ integer intent(hide),depend(y,ky),check(my>ky) :: my=len(y)
real*8 dimension(my) :: y
real*8 dimension(mx*my),depend(mx,my),check(len(z)==mx*my) :: z
real*8 optional,depend(x,mx) :: xb=dmin(x,mx)
@@ -456,7 +593,7 @@
integer intent(hide),depend(mx,my,nxest,nyest) &
:: kwrk=3+mx+my+nxest+nyest
integer intent(out) :: ier
- end subroutine regrid_smth
+ end subroutine regrid_smth0
end interface
end python module dfitpack
Modified: trunk/Lib/sandbox/spline/info.py
===================================================================
--- trunk/Lib/sandbox/spline/info.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/info.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -2,27 +2,31 @@
Spline Tools
===================
-Wrappers around FITPACK functions:
+This package contains wrappers around the netlib Dierckx library for curve
+and surface fitting with splines.
+Two interfaces are provided. The prefered interface is an object orientated
+interface defined in spline.spline:
+
+ UnivariateSpline -- 1 dimensional, smoothing or interpolating spline
+ InterpolatedUnivariateSpline -- i dimensional interpolating spline
+ LSQUnivariateSpline -- 1 dimensional least squares fitted spline
+ SmoothBivariateSpline -- 2 dimensional smoothing spline on scattered data
+ LSQBivariateSpline -- a least scquares fitted spline for scattered data
+ RectBivariateSpline -- a 2 dimensional smoothing or interpolating spline
+ on regular gridded rectangular data
+
+An alternative interface is a more direct wrapping of the fortran rountines:
+
splrep -- find smoothing spline given (x,y) points on curve.
splprep -- find smoothing spline given parametrically defined curve.
splev -- evaluate the spline or its derivatives.
splint -- compute definite integral of a spline.
sproot -- find the roots of a cubic spline.
spalde -- compute all derivatives of a spline at given points.
- bisplrep -- find bivariate smoothing spline representation.
+ bisplrep -- find bivariate smoothing spline representation.
bisplev -- evaluate bivariate smoothing spline.
- UnivariateSpline -- A more recent, object-oriented wrapper;
- finds a (possibly smoothed) interpolating
- spline.
- InterpolatedUnivariateSpline
- LSQUnivariateSpline
- BivariateSpline -- A more recent, object-oriented wrapper;
- finds a interpolating spline for a
- bivariate function.
-
- SmoothBivariateSpline
"""
postpone_import = 1
Modified: trunk/Lib/sandbox/spline/spline.py
===================================================================
--- trunk/Lib/sandbox/spline/spline.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/spline.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -6,7 +6,7 @@
to double routines by Pearu Peterson.
"""
# Created by Pearu Peterson, June,August 2003
-# Modified by John Travers, October 2006
+# Modified by John Travers, October-January 2006
__all__ = [
'UnivariateSpline',
@@ -84,7 +84,7 @@
deviation of y[i].
"""
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
- data = dfitpack.fpcurf0(x,y,k,w=w,
+ data = dfitpack.fpcurf_smth0(x,y,k,w=w,
xb=bbox[0],xe=bbox[1],s=s)
if data[-1]==1:
# nest too small, setting to maximum bound
@@ -130,7 +130,7 @@
fpint.resize(nest)
nrdata.resize(nest)
args = data[:8] + (t,c,n,fpint,nrdata,data[13])
- data = dfitpack.fpcurf1(*args)
+ data = dfitpack.fpcurf_smth1(*args)
return data
def set_smoothing_factor(self, s):
@@ -144,7 +144,7 @@
'LSQ spline with fixed knots')
return
args = data[:6] + (s,) + data[7:]
- data = dfitpack.fpcurf1(*args)
+ data = dfitpack.fpcurf_smth1(*args)
if data[-1]==1:
# nest too small, setting to maximum bound
data = self._reset_nest(data)
@@ -158,8 +158,10 @@
"""
if nu is None:
- return dfitpack.splev(*(self._eval_args+(x,)))
- return dfitpack.splder(nu=nu,*(self._eval_args+(x,)))
+ sp,ier = dfitpack.splev(*(self._eval_args+(x,)))
+ return sp
+ sp,ier = dfitpack.splder(nu=nu,*(self._eval_args+(x,)))
+ return sp
def get_knots(self):
""" Return the positions of (boundary and interior)
@@ -186,7 +188,8 @@
""" Return definite integral of the spline between two
given points.
"""
- return dfitpack.splint(*(self._eval_args+(a,b)))
+ iy, wrk = dfitpack.splint(*(self._eval_args+(a,b)))
+ return iy
def derivatives(self, x):
""" Return all derivatives of the spline at the point x."""
@@ -205,11 +208,10 @@
assert ier==0,`ier`
return z[:m]
raise NotImplementedError,\
- 'finding roots unsupported for non-cubic splines'
+ 'finding roots not supported for non-cubic splines'
class InterpolatedUnivariateSpline(UnivariateSpline):
- """ Interpolated univariate spline approximation. Identical to
- UnivariateSpline with less error checking.
+ """ Interpolated univariate spline approximation.
"""
@@ -227,14 +229,12 @@
k=3 - degree of the univariate spline.
"""
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
- self._data = dfitpack.fpcurf0(x,y,k,w=w,
+ self._data = dfitpack.fpcurf_smth0(x,y,k,w=w,
xb=bbox[0],xe=bbox[1],s=0)
self._reset_class()
class LSQUnivariateSpline(UnivariateSpline):
- """ Weighted least-squares univariate spline
- approximation. Appears to be identical to UnivariateSpline with
- more error checking.
+ """ Weighted least-squares univariate spline approximation.
"""
@@ -264,7 +264,7 @@
if not alltrue(t[k+1:n-k]-t[k:n-k-1] > 0,axis=0):
raise ValueError,\
'Interior knots t must satisfy Schoenberg-Whitney conditions'
- data = dfitpack.fpcurfm1(x,y,k,t,w=w,xb=xb,xe=xe)
+ data = dfitpack.fpcurf_lsq(x,y,k,t,w=w,xb=xb,xe=xe)
self._data = data[:-3] + (None,None,data[-1])
self._reset_class()
@@ -363,7 +363,6 @@
LSQUnivariateSpline - to create a BivariateSpline using weighted
least-squares fitting
"""
-
def __init__(self, x, y, z, w=None,
bbox = [None]*4, kx=3, ky=3, s=None, eps=None):
"""
@@ -388,7 +387,7 @@
equations. 0 < eps < 1, default is 1e-16.
"""
xb,xe,yb,ye = bbox
- nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth(x,y,z,w,
+ nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth0(x,y,z,w,
xb,xe,yb,ye,
kx,ky,s=s,
eps=eps,lwrk2=1)
@@ -411,7 +410,6 @@
SmoothUnivariateSpline - to create a BivariateSpline through the
given points
"""
-
def __init__(self, x, y, z, tx, ty, w=None,
bbox = [None]*4,
kx=3, ky=3, eps=None):
@@ -471,7 +469,6 @@
bisplrep, bisplev - an older wrapping of FITPACK
UnivariateSpline - a similar class for univariate spline interpolation
"""
-
def __init__(self, x, y, z,
bbox = [None]*4, kx=3, ky=3, s=0):
"""
@@ -506,7 +503,7 @@
'y dimension of z must have same number of elements as y'
z = ravel(z)
xb,xe,yb,ye = bbox
- nx,tx,ny,ty,c,fp,ier = dfitpack.regrid_smth(x,y,z,
+ nx,tx,ny,ty,c,fp,ier = dfitpack.regrid_smth0(x,y,z,
xb,xe,yb,ye,
kx,ky,s)
if ier in [0,-1,-2]: # normal return
@@ -518,4 +515,3 @@
self.fp = fp
self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
self.degrees = kx,ky
-
Deleted: trunk/Lib/sandbox/spline/tests/demos_xplt.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/demos_xplt.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/tests/demos_xplt.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -1,54 +0,0 @@
-#!/usr/bin/env python
-# Created by Pearu Peterson, Aug 2003
-""" Test xplt based demos for interpolate.fitpack2 module
-"""
-__usage__ = """
-Build interpolate:
- python setup_interpolate.py build
-Run demos (assumes that scipy is installed):
- python -i tests/demos_xplt.py
-"""
-
-import sys
-from numpy.test.testing import set_package_path
-set_package_path()
-from interpolate.fitpack2 import UnivariateSpline,LSQUnivariateSpline,\
- InterpolatedUnivariateSpline
-from interpolate.fitpack2 import LSQBivariateSpline, SmoothBivariateSpline
-del sys.path[0]
-
-from scipy import *
-
-def demo1():
- x = arange(0,2*pi+pi/4,2*pi/8)
- xnew = arange(-pi/10,2*pi+pi/4+pi/10,pi/50)
- y = sin(x)
-
-
- def make_plot():
- xplt.plot(x,y,'x',xnew,spline(xnew),x,y,'b',xnew,sin(xnew),
- spline.get_knots(),spline(spline.get_knots()),'o')
-
- spline = UnivariateSpline(x,y,k=1)
- assert isinstance(spline,LSQUnivariateSpline)
- print 'Linear LSQ approximation of sin(x):',spline.__class__.__name__
- make_plot()
- print 'Residual=',spline.get_residual()
- raw_input('Press any key to continue..')
-
- spline.set_smoothing_factor(0)
- assert isinstance(spline,InterpolatedUnivariateSpline)
- print 'Linear interpolation of sin(x):',spline.__class__.__name__
- make_plot()
- print 'Residual=',spline.get_residual()
- raw_input('Press any key to continue..')
-
- spline = UnivariateSpline(x,y,k=1,s=0.1)
- print 'Linear smooth approximation of sin(x):',spline.__class__.__name__
- assert isinstance(spline,UnivariateSpline)
- make_plot()
- print 'Residual=',spline.get_residual()
- raw_input('Press any key to continue..')
-
-if __name__ == "__main__":
- demo1()
Added: trunk/Lib/sandbox/spline/tests/dierckx_test_data.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/dierckx_test_data.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/tests/dierckx_test_data.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -0,0 +1,184 @@
+from numpy import arange, array, sin, cos, atleast_1d, transpose, dot
+
+curfit_test = {
+#inputs:
+'y' : array([1.0,1.0,1.4,1.1,1.0,1.0,4.0,9.0,13.0,
+ 13.4,12.8,13.1,13.0,14.0,13.0,13.5,10.0,
+ 2.0,3.0,2.5,2.5,2.5,3.0,4.0,3.5]),
+'x' : arange(25.0)}
+
+curfit_test_smth = {
+#inputs:
+'k' : [3, 3, 3, 5, 5, 5],
+'s' : [1000, 60, 0, 60, 10, 0],
+'iopt' : [0, 1, 1, 0, 1, 1],
+# results:
+'fp' : array([0.265483E+03, 0.600380E+02, 0.000000E+00, 0.600201E+02,
+ 0.100002E+02, 0.000000E+00]),
+'ier' : array([-2, 0, -1, 0, 0, -1]),
+'n' : array([8, 11, 29, 13, 21, 31]),
+'t' : [array([0.0, 0.0, 0.0, 0.0, 24.0, 24.0, 24.0, 24.0]),
+ array([0.0, 0.0, 0.0, 0.0, 12.0, 15.0, 18.0, 24.0,
+ 24.0, 24.0, 24.0]),
+ array([0.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0, 5.0,
+ 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0,
+ 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0,
+ 22.0, 24.0, 24.0, 24.0, 24.0]),
+ array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.0, 24.0,
+ 24.0, 24.0, 24.0, 24.0, 24.0]),
+ array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 8.0,
+ 9.0, 12.0, 14.0, 15.0, 16.0, 17.0, 18.0, 24.0,
+ 24.0, 24.0, 24.0, 24.0, 24.0]),
+ array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 4.0,
+ 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
+ 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0,
+ 21.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0])],
+'c' : [array([-3.5307, 19.6679, 10.0310, 0.5226]),
+ array([1.7450, -6.1367, 21.0990, 11.1564, -0.1060, 2.0998,
+ 4.3288]),
+ array([1.0000, 0.4704, 1.8592, 0.9446, 1.1911, 0.2910,
+ 3.6448, 9.1297, 13.8362, 13.5254, 12.4622, 13.4258,
+ 12.4346, 14.8358, 12.2222, 14.2755, 11.6757, -0.9783,
+ 4.2375, 2.0285, 2.6486, 2.3770, 2.9986, 4.8757,
+ 3.5000]),
+ array([1.9049, -2.6240, 0.5773, 45.0508, -10.0998, 2.7788,
+ 4.2385]),
+ array([1.0179, -1.2191, 9.2752, -12.6640, 17.3909, 13.0655,
+ 11.7046, 14.2824, 14.9330, 8.7780, -2.6436, 5.4611,
+ 0.5460, 4.9279, 3.5774]),
+ array([1.0000, -0.7624, 3.5842, -0.4712, 2.4346, -0.3988,
+ 3.6292, 9.0283, 14.4832, 13.6084, 11.9251, 14.2042,
+ 11.3359, 16.1291, 11.2375, 14.3636, 14.1766, -4.8942,
+ 6.9048, 0.5376, 4.5146, 0.4259, 4.7049, 4.4036,
+ 3.5000])],
+'sp' : [array([-3.5, -0.8, 1.6, 3.7, 5.5, 7.0, 8.2, 9.2, 9.9, 10.5,
+ 10.8, 10.8, 10.8, 10.5, 10.1, 9.5, 8.9, 8.1, 7.2, 6.2,
+ 5.1, 4.0, 2.9, 1.7, 0.5]),
+ array([1.7, 0.4, 0.0, 0.6, 1.8, 3.6, 5.6, 7.8, 10.0, 12.0,
+ 13.5, 14.5, 14.7, 14.0, 12.6, 10.6, 8.2, 5.9, 3.9, 2.6,
+ 2.1, 2.1, 2.5, 3.3, 4.3]),
+ array([1.0, 1.0, 1.4, 1.1, 1.0, 1.0, 4.0, 9.0, 13.0, 13.4,
+ 12.8, 13.1, 13.0, 14.0, 13.0, 13.5, 10.0, 2.0, 3.0, 2.5,
+ 2.5, 2.5, 3.0, 4.0, 3.5]),
+ array([1.9, 0.5, -0.1, 0.3, 1.5, 3.4, 5.6, 8.0, 10.3, 12.3,
+ 13.8, 14.5, 14.5, 13.7, 12.3, 10.3, 8.2, 6.0, 4.1, 2.8,
+ 2.0, 2.0, 2.5, 3.4, 4.2]),
+ array([1.0, 0.9, 1.5, 1.2, 0.6, 1.2, 4.2, 8.9, 12.6, 13.6,
+ 13.2, 12.7, 12.9, 13.7, 14.2, 12.9, 8.6, 4.1, 2.3, 2.2,
+ 2.4, 2.7, 3.2, 3.7, 3.6]),
+ array([1.0, 1.0, 1.4, 1.1, 1.0, 1.0, 4.0, 9.0, 13.0, 13.4,
+ 12.8, 13.1, 13.0, 14.0, 13.0, 13.5, 10.0, 2.0, 3.0, 2.5,
+ 2.5, 2.5, 3.0, 4.0, 3.5])]
+}
+
+curfit_test_lsq = {
+#inputs:
+'k' : [3, 5],
+# results:
+'fp' : array([0.234182E+02, 0.155243E+02]),
+'ier' : array([0, 0]),
+'t' : [array([3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0]),
+ array([3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0])],
+'c' : [array([0.8505, 2.4795, -0.9480, 2.8341, 17.4216, 10.3854,
+ 17.1927, -2.9343, 5.6395, 2.4292, 3.7180]),
+ array([1.0199, -1.6630, 6.3167, -4.5892, 1.6163, 22.6478,
+ 3.5130, 27.5439, -16.2085, 15.0884, -5.3102, 8.6487,
+ 3.4724])],
+'sp' : [array([0.9, 1.5, 1.1, 0.5, 0.8, 2.2, 4.6, 8.1, 11.6, 13.8,
+ 14.1, 13.3, 12.7, 13.2, 13.7, 12.7, 9.3, 5.0, 1.8, 1.3,
+ 2.3, 3.4, 3.5, 3.2, 3.7]),
+ array([1.0, 0.8, 1.8, 0.9, 0.3, 1.6, 4.7, 8.6, 12.0, 13.7,
+ 13.7, 12.8, 12.7, 13.5, 14.2, 12.8, 8.9, 4.4, 1.8, 1.9,
+ 3.0, 3.0, 2.4, 4.2, 3.5])]
+}
+
+def f1(x,d=0):
+ if d is None: return "sin"
+ if x is None: return "sin(x)"
+ if d%4 == 0: return sin(x)
+ if d%4 == 1: return cos(x)
+ if d%4 == 2: return -sin(x)
+ if d%4 == 3: return -cos(x)
+
+def f2(x,y=0,dx=0,dy=0):
+ if x is None: return "sin(x+y)"
+ d=dx+dy
+ if d%4 == 0: return sin(x+y)
+ if d%4 == 1: return cos(x+y)
+ if d%4 == 2: return -sin(x+y)
+ if d%4 == 3: return -cos(x+y)
+
+myasarray = atleast_1d
+
+def norm2(x):
+ return dot(transpose(x),x)
+def makepairs(x,y):
+ x,y=map(myasarray,[x,y])
+ xy=array(map(lambda x,y:map(None,len(y)*[x],y),x,len(x)*[y]))
+ sh=xy.shape
+ xy.shape=sh[0]*sh[1],sh[2]
+ return transpose(xy)
+
+# very simple script to test interpolation with regrid and surfit
+# based on example data from netlib->dierckx->regrid
+
+#from numpy import *
+#from scipy.interpolate.fitpack2 import SmoothBivariateSpline, \
+ #RectBivariateSpline
+#import matplotlib
+#matplotlib.use('Agg')
+#import pylab
+
+## x,y coordinates
+#x = linspace(-1.5,1.5,11)
+#y = x
+## data taken from daregr
+#z = array([
+#[-0.0325, 0.0784, 0.0432, 0.0092, 0.1523, 0.0802, 0.0925, -0.0098, \
+ #0.0810, -0.0146, -0.0019],
+#[0.1276, 0.0223, 0.0357, 0.1858, 0.2818, 0.1675, 0.2239, 0.1671, \
+ #0.0843, 0.0151, 0.0427],
+#[0.0860, 0.1267, 0.1839, 0.3010, 0.5002, 0.4683, 0.4562, 0.2688, \
+ #0.1276, 0.1244, 0.0377],
+#[0.0802, 0.1803, 0.3055, 0.4403, 0.6116, 0.7178, 0.6797, 0.5218, \
+ #0.2624, 0.1341, -0.0233],
+#[0.1321, 0.2023, 0.4446, 0.7123, 0.7944, 0.9871, 0.8430, 0.6440, \
+ #0.4682, 0.1319, 0.1075],
+#[0.2561, 0.1900, 0.4614, 0.7322, 0.9777, 1.0463, 0.9481, 0.6649, \
+ #0.4491, 0.2442, 0.1341],
+#[0.0981, 0.2009, 0.4616, 0.5514, 0.7692, 0.9831, 0.7972, 0.5937, \
+ #0.4190, 0.1436, 0.0995],
+#[0.0991, 0.1545, 0.3399, 0.4940, 0.6328, 0.7168, 0.6886, 0.3925, \
+ #0.3015, 0.1758, 0.0928],
+#[-0.0197, 0.1479, 0.1225, 0.3254, 0.3847, 0.4767, 0.4324, 0.2827, \
+ #0.2287, 0.0999, 0.0785],
+#[0.0032, 0.0917, 0.0246, 0.1780, 0.2394, 0.1765, 0.1642, 0.2081, \
+ #0.1049, 0.0493, -0.0502],
+#[0.0101, 0.0297, 0.0468, 0.0221, 0.1074, 0.0433, 0.0626, 0.1436, \
+ #0.1092, -0.0232, 0.0132]])
+
+## plot original data
+#pylab.subplot(1,3,1)
+#pylab.imshow(z)
+#pylab.title('orig')
+
+## check regrid
+#mrs = RectBivariateSpline(x,y,z)
+#zr = mrs(x,y)
+#print sum(abs(zr-z))
+#pylab.subplot(1,3,2)
+#pylab.imshow(zr)
+#pylab.title('regrid')
+
+## check surfit
+## x increases in columns, y increases along rows
+#ym,xm = meshgrid(y,x) # deal with meshgrid badness
+#mrs = SmoothBivariateSpline(ravel(xm),ravel(ym),ravel(z),kx=3,ky=3,s=0)
+#zr = mrs(x,y)
+#print sum(abs(zr-z))
+#pylab.subplot(1,3,3)
+#pylab.imshow(zr)
+#pylab.title('surfit')
+
+#pylab.savefig('intplot.png')
+#pylab.close()
Deleted: trunk/Lib/sandbox/spline/tests/test_fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/test_fitpack.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/tests/test_fitpack.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -1,87 +0,0 @@
-#!/usr/bin/env python
-# Created by Pearu Peterson, June 2003
-""" Test functions for interpolate.fitpack2 module
-"""
-__usage__ = """
-Build interpolate:
- python setup_interpolate.py build
-Run tests if scipy is installed:
- python -c 'import scipy;scipy.interpolate.test(<level>)'
-Run tests if interpolate is not installed:
- python tests/test_fitpack.py [<level>]
-"""
-#import libwadpy
-
-import sys
-from numpy.testing import *
-from numpy import array
-set_package_path()
-from interpolate.fitpack2 import UnivariateSpline,LSQUnivariateSpline,\
- InterpolatedUnivariateSpline
-from interpolate.fitpack2 import LSQBivariateSpline, SmoothBivariateSpline,\
- RectBivariateSpline
-restore_path()
-
-class test_UnivariateSpline(NumpyTestCase):
- def check_linear_constant(self):
- x = [1,2,3]
- y = [3,3,3]
- lut = UnivariateSpline(x,y,k=1)
- assert_array_almost_equal(lut.get_knots(),[1,3])
- assert_array_almost_equal(lut.get_coeffs(),[3,3])
- assert_almost_equal(lut.get_residual(),0.0)
- assert_array_almost_equal(lut([1,1.5,2]),[3,3,3])
-
- def check_linear_1d(self):
- x = [1,2,3]
- y = [0,2,4]
- lut = UnivariateSpline(x,y,k=1)
- assert_array_almost_equal(lut.get_knots(),[1,3])
- assert_array_almost_equal(lut.get_coeffs(),[0,4])
- assert_almost_equal(lut.get_residual(),0.0)
- assert_array_almost_equal(lut([1,1.5,2]),[0,1,2])
-
-class test_LSQBivariateSpline(NumpyTestCase):
- def check_linear_constant(self):
- x = [1,1,1,2,2,2,3,3,3]
- y = [1,2,3,1,2,3,1,2,3]
- z = [3,3,3,3,3,3,3,3,3]
- s = 0.1
- tx = [1+s,3-s]
- ty = [1+s,3-s]
- lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1)
- #print lut.get_knots()
- #print lut.get_coeffs()
- #print lut.get_residual()
-
-class test_SmoothBivariateSpline(NumpyTestCase):
- def check_linear_constant(self):
- x = [1,1,1,2,2,2,3,3,3]
- y = [1,2,3,1,2,3,1,2,3]
- z = [3,3,3,3,3,3,3,3,3]
- lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1)
- assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3]))
- assert_array_almost_equal(lut.get_coeffs(),[3,3,3,3])
- assert_almost_equal(lut.get_residual(),0.0)
- assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[3,3],[3,3],[3,3]])
-
- def check_linear_1d(self):
- x = [1,1,1,2,2,2,3,3,3]
- y = [1,2,3,1,2,3,1,2,3]
- z = [0,0,0,2,2,2,4,4,4]
- lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1)
- assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3]))
- assert_array_almost_equal(lut.get_coeffs(),[0,0,4,4])
- assert_almost_equal(lut.get_residual(),0.0)
- assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[0,0],[1,1],[2,2]])
-
-class test_RectBivariateSpline(NumpyTestCase):
- def check_defaults(self):
- x = array([1,2,3,4,5])
- y = array([1,2,3,4,5])
- z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]])
- lut = RectBivariateSpline(x,y,z)
- assert_array_almost_equal(lut(x,y),z)
-
-if __name__ == "__main__":
- NumpyTest().run()
Added: trunk/Lib/sandbox/spline/tests/test_fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/test_fitpack.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/tests/test_fitpack.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+# Created by John Travers, January 2007
+""" Test functions for spline.fitpack module
+"""
+__usage__ = """
+Build spline:
+ python setup_spline.py build
+Run tests if scipy is installed:
+ python -c 'import scipy;scipy.spline.test(<level>)'
+Run tests if spline is not installed:
+ python tests/test_fitpack.py [<level>]
+"""
+
+import sys
+from numpy.testing import *
+from numpy import array, arange, around, pi, sin, ravel
+
+set_package_path()
+from spline.fitpack import splrep, splev, sproot, splint, spalde
+from spline.fitpack import bisplev, bisplrep, splprep
+restore_path()
+
+set_local_path()
+from dierckx_test_data import *
+restore_path()
+
+class test_splrep_slev(NumpyTestCase):
+ def check_curfit_against_dierckx_smth(self):
+ x,y = curfit_test['x'],curfit_test['y']
+ k,s = curfit_test_smth['k'],curfit_test_smth['s']
+ iopt = curfit_test_smth['iopt']
+ for i in range(len(k)):
+ tck = splrep(x,y,k=k[i],task=iopt[i],s=s[i])
+ out = splrep(x,y,k=k[i],task=iopt[i],s=s[i],full_output=1)
+ tck,fp=out[0],out[1]
+ sp = splev(x,tck)
+ assert_almost_equal(fp, curfit_test_smth['fp'][i], decimal=2)
+ assert_array_almost_equal(around(tck[0],1),
+ curfit_test_smth['t'][i])
+ assert_array_almost_equal(around(tck[1],4),
+ curfit_test_smth['c'][i], decimal=3)
+ assert_array_almost_equal(around(sp,1),
+ curfit_test_smth['sp'][i])
+
+ def check_curfit_against_dierckx_lsq(self):
+ """ Test against results obtined from the pure fortran routines.
+
+ Here we check simple spline creation and evaluation.
+ """
+ x,y = curfit_test['x'],curfit_test['y']
+ k = curfit_test_lsq['k']
+ for i in range(len(k)):
+ t = curfit_test_lsq['t'][i]
+ out = splrep(x,y,t=t,k=k[i],full_output=1)
+ tck,fp=out[0],out[1]
+ sp = splev(x,tck)
+ assert_almost_equal(fp, curfit_test_lsq['fp'][i], decimal=2)
+ assert_array_almost_equal(around(tck[1],4),
+ curfit_test_lsq['c'][i], decimal=3)
+ assert_array_almost_equal(around(sp,1),
+ curfit_test_lsq['sp'][i])
+
+class test_splint_spalde(NumpyTestCase):
+ def check_splint_spalde(self):
+ per = [0, 1, 0]
+ N = [20, 20, 50]
+ ia = [0, 0, 0.2*pi]
+ ib = [0, 0, pi]
+ a,b = 0,2*pi
+ dx = 0.2*pi
+ k = range(1,6)
+ for i in range(len(per)):
+ x=a+(b-a)*arange(N[i]+1,dtype=float)/float(N[i])
+ v=f1(x)
+ for j in range(len(k)):
+ tck = splrep(x,v,k=k[j],s=0,per=per[i])
+ ir = splint(ia[i],ib[i],tck)
+ dr = spalde(dx,tck)
+ assert_almost_equal(ir, f1(ib[i],-1)-f1(ia[i],-1), decimal=2)
+ d=0
+ for ddr in dr:
+ if d<k[j]-1:
+ assert_almost_equal(1, ddr/f1(dx,d), decimal=2)
+ d=d+1
+
+class test_splder(NumpyTestCase):
+ def check_splder(self):
+ N = 50
+ a,b = 0,2*pi
+ dx = 0.2*pi
+ k = range(3,6)
+ x=a+(b-a)*arange(N+1,dtype=float)/float(N)
+ v=f1(x)
+ for j in range(len(k)):
+ tck = splrep(x,v,k=k[j],s=0)
+ dr = spalde(dx,tck)
+ d2 = splev(dx,tck,der=1)
+ assert_almost_equal(1, dr[1]/f1(dx,1.0), decimal=2)
+
+class test_sproot(NumpyTestCase):
+ def check_sproot(self):
+ a=0
+ b=15
+ N=20
+ x=a+(b-a)*arange(N+1,dtype=float)/float(N)
+ v=f1(x)
+ k=3
+ tck = splrep(x,v,k=k,s=0)
+ ex = array([0.0, pi, 2.0*pi, 3.0*pi, 4.0*pi])
+ assert_array_almost_equal(sproot(tck),ex, decimal=3)
+
+class test_bisplev_bisplrep(NumpyTestCase):
+ def test_bisplev_bisplrep(self):
+ f=f2; kx=3; ky=3; xb=0; xe=2*pi
+ yb=0; ye=2*pi; Nx=20; Ny=20; s=0
+ x=xb+(xe-xb)*arange(Nx+1,dtype=float)/float(Nx)
+ y=yb+(ye-yb)*arange(Ny+1,dtype=float)/float(Ny)
+ xy=makepairs(x,y)
+ tck=bisplrep(xy[0],xy[1],f(xy[0],xy[1]),s=s,kx=kx,ky=ky)
+ tt=[tck[0][kx:-kx],tck[1][ky:-ky]]
+ t2=makepairs(tt[0],tt[1])
+ v1=bisplev(tt[0],tt[1],tck)
+ v2=f2(t2[0],t2[1])
+ v2.shape=len(tt[0]),len(tt[1])
+ assert_almost_equal(0.0, norm2(ravel(v1-v2)), decimal=5)
+
+class test_parcur(NumpyTestCase):
+ def check_parcur(self):
+ f=f1; per=0; s=0; a=0; b=2*pi;
+ N=[20,50]
+ ia=0; ib=2*pi; dx=0.2*pi
+ xb,xe=a,b
+ for i in range(len(N)):
+ x=a+(b-a)*arange(N[i]+1,dtype=float)/float(N[i])
+ x1=a+(b-a)*arange(1,N[i],dtype=float)/float(N[i]-1) # mid points
+ v,v1=f(x),f(x1)
+ for k in range(1,6):
+ tckp,u=splprep([x,v],s=s,per=per,k=k)
+ tck=splrep(x,v,s=s,per=per,k=k)
+ uv=splev(dx,tckp)
+ assert_almost_equal(0.0, around(abs(uv[1]-f(uv[0])),2),
+ decimal=1)
+ assert_almost_equal(0.0,
+ around(abs(splev(uv[0],tck)-f(uv[0])),2),decimal=1)
Deleted: trunk/Lib/sandbox/spline/tests/test_interpolate.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/test_interpolate.py 2007-01-14 22:48:12 UTC (rev 2554)
+++ trunk/Lib/sandbox/spline/tests/test_interpolate.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -1,198 +0,0 @@
-from numpy.testing import *
-from numpy import mgrid, pi, sin, ogrid
-import numpy as np
-
-set_package_path()
-from interpolate import interp1d, interp2d
-restore_path()
-
-
-class test_interp2d(NumpyTestCase):
- def test_interp2d(self):
- y, x = mgrid[0:pi:20j, 0:pi:21j]
- z = sin(x+y)
- I = interp2d(x, y, z)
- assert_almost_equal(I(1.0, 1.0), sin(2.0), decimal=2)
-
- v,u = ogrid[0:pi:24j, 0:pi:25j]
- assert_almost_equal(I(u.ravel(), v.ravel()), sin(v+u), decimal=2)
-
-
-class test_interp1d(NumpyTestCase):
-
- def setUp(self):
- self.x10 = np.arange(10.)
- self.y10 = np.arange(10.)
- self.x25 = self.x10.reshape((2,5))
- self.x2 = np.arange(2.)
- self.y2 = np.arange(2.)
- self.x1 = np.array([0.])
- self.y1 = np.array([0.])
-
- self.y210 = np.arange(20.).reshape((2, 10))
- self.y102 = np.arange(20.).reshape((10, 2))
-
- self.fill_value = -100.0
-
- def test_validation(self):
- """ Make sure that appropriate exceptions are raised when invalid values
- are given to the constructor.
- """
-
- # Only kind='linear' is implemented.
- self.assertRaises(NotImplementedError, interp1d, self.x10, self.y10, kind='cubic')
- interp1d(self.x10, self.y10, kind='linear')
-
- # x array must be 1D.
- self.assertRaises(ValueError, interp1d, self.x25, self.y10)
-
- # y array cannot be a scalar.
- self.assertRaises(ValueError, interp1d, self.x10, np.array(0))
-
- # Check for x and y arrays having the same length.
- self.assertRaises(ValueError, interp1d, self.x10, self.y2)
- self.assertRaises(ValueError, interp1d, self.x2, self.y10)
- self.assertRaises(ValueError, interp1d, self.x10, self.y102)
- interp1d(self.x10, self.y210)
- interp1d(self.x10, self.y102, axis=0)
-
- # Check for x and y having at least 1 element.
- self.assertRaises(ValueError, interp1d, self.x1, self.y10)
- self.assertRaises(ValueError, interp1d, self.x10, self.y1)
- self.assertRaises(ValueError, interp1d, self.x1, self.y1)
-
-
- def test_init(self):
- """ Check that the attributes are initialized appropriately by the
- constructor.
- """
-
- self.assert_(interp1d(self.x10, self.y10).copy)
- self.assert_(not interp1d(self.x10, self.y10, copy=False).copy)
- self.assert_(interp1d(self.x10, self.y10).bounds_error)
- self.assert_(not interp1d(self.x10, self.y10, bounds_error=False).bounds_error)
- self.assert_(np.isnan(interp1d(self.x10, self.y10).fill_value))
- self.assertEqual(
- interp1d(self.x10, self.y10, fill_value=3.0).fill_value,
- 3.0,
- )
- self.assertEqual(
- interp1d(self.x10, self.y10).axis,
- 0,
- )
- self.assertEqual(
- interp1d(self.x10, self.y210).axis,
- 1,
- )
- self.assertEqual(
- interp1d(self.x10, self.y102, axis=0).axis,
- 0,
- )
- assert_array_equal(
- interp1d(self.x10, self.y10).x,
- self.x10,
- )
- assert_array_equal(
- interp1d(self.x10, self.y10).y,
- self.y10,
- )
- assert_array_equal(
- interp1d(self.x10, self.y210).y,
- self.y210,
- )
-
-
- def test_linear(self):
- """ Check the actual implementation of linear interpolation.
- """
-
- interp10 = interp1d(self.x10, self.y10)
- assert_array_almost_equal(
- interp10(self.x10),
- self.y10,
- )
- assert_array_almost_equal(
- interp10(1.2),
- np.array([1.2]),
- )
- assert_array_almost_equal(
- interp10([2.4, 5.6, 6.0]),
- np.array([2.4, 5.6, 6.0]),
- )
-
-
- def test_bounds(self):
- """ Test that our handling of out-of-bounds input is correct.
- """
-
- extrap10 = interp1d(self.x10, self.y10, fill_value=self.fill_value,
- bounds_error=False)
- assert_array_equal(
- extrap10(11.2),
- np.array([self.fill_value]),
- )
- assert_array_equal(
- extrap10(-3.4),
- np.array([self.fill_value]),
- )
- assert_array_equal(
- extrap10._check_bounds(np.array([-1.0, 0.0, 5.0, 9.0, 11.0])),
- np.array([True, False, False, False, True]),
- )
-
- raises_bounds_error = interp1d(self.x10, self.y10, bounds_error=True)
- self.assertRaises(ValueError, raises_bounds_error, -1.0)
- self.assertRaises(ValueError, raises_bounds_error, 11.0)
- raises_bounds_error([0.0, 5.0, 9.0])
-
-
- def test_nd(self):
- """ Check the behavior when the inputs and outputs are multidimensional.
- """
-
- # Multidimensional input.
- interp10 = interp1d(self.x10, self.y10)
- assert_array_almost_equal(
- interp10(np.array([[3.4, 5.6], [2.4, 7.8]])),
- np.array([[3.4, 5.6], [2.4, 7.8]]),
- )
-
- # Multidimensional outputs.
- interp210 = interp1d(self.x10, self.y210)
- assert_array_almost_equal(
- interp210(1.5),
- np.array([[1.5], [11.5]]),
- )
- assert_array_almost_equal(
- interp210(np.array([1.5, 2.4])),
- np.array([[1.5, 2.4],
- [11.5, 12.4]]),
- )
-
- interp102 = interp1d(self.x10, self.y102, axis=0)
- assert_array_almost_equal(
- interp102(1.5),
- np.array([[3.0, 4.0]]),
- )
- assert_array_almost_equal(
- interp102(np.array([1.5, 2.4])),
- np.array([[3.0, 4.0],
- [4.8, 5.8]]),
- )
-
- # Both at the same time!
- x_new = np.array([[3.4, 5.6], [2.4, 7.8]])
- assert_array_almost_equal(
- interp210(x_new),
- np.array([[[3.4, 5.6], [2.4, 7.8]],
- [[13.4, 15.6], [12.4, 17.8]]]),
- )
- assert_array_almost_equal(
- interp102(x_new),
- np.array([[[6.8, 7.8], [11.2, 12.2]],
- [[4.8, 5.8], [15.6, 16.6]]]),
- )
-
-
-if __name__ == "__main__":
- NumpyTest().run()
Copied: trunk/Lib/sandbox/spline/tests/test_spline.py (from rev 2535, trunk/Lib/sandbox/spline/tests/test_fitpack.py)
===================================================================
--- trunk/Lib/sandbox/spline/tests/test_fitpack.py 2007-01-11 09:43:18 UTC (rev 2535)
+++ trunk/Lib/sandbox/spline/tests/test_spline.py 2007-01-15 00:23:52 UTC (rev 2555)
@@ -0,0 +1,167 @@
+#!/usr/bin/env python
+# Created by Pearu Peterson, June 2003
+# Modified by John Travers, October 2006
+""" Test functions for spline.spline module
+"""
+__usage__ = """
+Build spline:
+ python setup_spline.py build
+Run tests if scipy is installed:
+ python -c 'import scipy;scipy.spline.test(<level>)'
+Run tests if spline is not installed:
+ python tests/test_spline.py [<level>]
+"""
+
+import sys
+from numpy.testing import *
+from numpy import array, arange, around, pi, sin, cos
+
+set_package_path()
+from spline.spline import UnivariateSpline,LSQUnivariateSpline,\
+ InterpolatedUnivariateSpline
+from spline.spline import LSQBivariateSpline, SmoothBivariateSpline,\
+ RectBivariateSpline
+restore_path()
+
+set_local_path()
+from dierckx_test_data import *
+restore_path()
+
+class test_UnivariateSpline(NumpyTestCase):
+ def check_linear_constant(self):
+ x = [1,2,3]
+ y = [3,3,3]
+ lut = UnivariateSpline(x,y,k=1)
+ assert_array_almost_equal(lut.get_knots(),[1,3])
+ assert_array_almost_equal(lut.get_coeffs(),[3,3])
+ assert_almost_equal(lut.get_residual(),0.0)
+ assert_array_almost_equal(lut([1,1.5,2]),[3,3,3])
+
+ def check_linear_1d(self):
+ x = [1,2,3]
+ y = [0,2,4]
+ lut = UnivariateSpline(x,y,k=1)
+ assert_array_almost_equal(lut.get_knots(),[1,3])
+ assert_array_almost_equal(lut.get_coeffs(),[0,4])
+ assert_almost_equal(lut.get_residual(),0.0)
+ assert_array_almost_equal(lut([1,1.5,2]),[0,1,2])
+
+ def check_curfit_against_dierckx(self):
+ """ Test against results obtined from the pure fortran routines.
+
+ Here we check simple spline creation and evaluation.
+ """
+ x,y = curfit_test['x'],curfit_test['y']
+ k,s = curfit_test_smth['k'],curfit_test_smth['s']
+ iopt = curfit_test_smth['iopt']
+ for i in range(len(k)):
+ if iopt[i] == 0:
+ uspl = UnivariateSpline(x,y,k=k[i],s=s[i])
+ elif iopt[i] == 1:
+ uspl.set_smoothing_factor(s[i])
+ assert_almost_equal(uspl.get_residual(),
+ curfit_test_smth['fp'][i], decimal=2)
+ n = uspl._data[7]
+ assert_equal(n,curfit_test_smth['n'][i])
+ assert_array_almost_equal(around(uspl.get_knots(),1),
+ curfit_test_smth['t'][i][k[i]:n-k[i]])
+ assert_array_almost_equal(around(uspl.get_coeffs(),4),
+ curfit_test_smth['c'][i], decimal=3)
+ assert_array_almost_equal(around(uspl(x),1),
+ curfit_test_smth['sp'][i])
+
+ def check_spint_spalde(self):
+ per = [0, 0, 0]
+ N = [20, 20, 50]
+ ia = [0, 0, 0.2*pi]
+ ib = [0, 0, pi]
+ a,b = 0,2*pi
+ dx = 0.2*pi
+ k = range(1,6)
+ for i in range(len(per)):
+ x=a+(b-a)*arange(N[i]+1,dtype=float)/float(N[i])
+ v=f1(x)
+ for j in range(len(k)):
+ uspl = UnivariateSpline(x,v,k=k[j],s=0)
+ ir = uspl.integral(ia[i],ib[i])
+ dr = uspl.derivatives(dx)
+ assert_almost_equal(ir, f1(ib[i],-1)-f1(ia[i],-1), decimal=2)
+ d=0
+ for ddr in dr:
+ if d<k[j]-1:
+ assert_almost_equal(1, ddr/f1(dx,d), decimal=2)
+ d=d+1
+
+ def check_sproot(self):
+ a=0
+ b=15
+ N=20
+ x=a+(b-a)*arange(N+1,dtype=float)/float(N)
+ v=f1(x)
+ k=3
+ uspl = UnivariateSpline(x,v,k=k,s=0)
+ ex = array([0.0, pi, 2.0*pi, 3.0*pi, 4.0*pi])
+ assert_array_almost_equal(uspl.roots(),ex, decimal=3)
+
+class test_LSQUnivariateSpline(NumpyTestCase):
+ def check_curfit_against_dierckx(self):
+ """ Test against results obtined from the pure fortran routines.
+
+ Here we check simple spline creation and evaluation.
+ """
+ x,y = curfit_test['x'],curfit_test['y']
+ k = curfit_test_lsq['k']
+ for i in range(len(k)):
+ t = curfit_test_lsq['t'][i]
+ lsquspl = LSQUnivariateSpline(x,y,t,k=k[i])
+ assert_almost_equal(lsquspl.get_residual(),
+ curfit_test_lsq['fp'][i], decimal=2)
+ assert_array_almost_equal(around(lsquspl.get_coeffs(),4),
+ curfit_test_lsq['c'][i], decimal=3)
+ assert_array_almost_equal(around(lsquspl(x),1),
+ curfit_test_lsq['sp'][i])
+
+class test_LSQBivariateSpline(NumpyTestCase):
+ def check_linear_constant(self):
+ x = [1,1,1,2,2,2,3,3,3]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = [3,3,3,3,3,3,3,3,3]
+ s = 0.1
+ tx = [1+s,3-s]
+ ty = [1+s,3-s]
+ lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1)
+ #print lut.get_knots()
+ #print lut.get_coeffs()
+ #print lut.get_residual()
+
+class test_SmoothBivariateSpline(NumpyTestCase):
+ def check_linear_constant(self):
+ x = [1,1,1,2,2,2,3,3,3]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = [3,3,3,3,3,3,3,3,3]
+ lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1)
+ assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3]))
+ assert_array_almost_equal(lut.get_coeffs(),[3,3,3,3])
+ assert_almost_equal(lut.get_residual(),0.0)
+ assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[3,3],[3,3],[3,3]])
+
+ def check_linear_1d(self):
+ x = [1,1,1,2,2,2,3,3,3]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = [0,0,0,2,2,2,4,4,4]
+ lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1)
+ assert_array_almost_equal(lut.get_knots(),([1,1,3,3],[1,1,3,3]))
+ assert_array_almost_equal(lut.get_coeffs(),[0,0,4,4])
+ assert_almost_equal(lut.get_residual(),0.0)
+ assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[0,0],[1,1],[2,2]])
+
+class test_RectBivariateSpline(NumpyTestCase):
+ def check_defaults(self):
+ x = array([1,2,3,4,5])
+ y = array([1,2,3,4,5])
+ z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]])
+ lut = RectBivariateSpline(x,y,z)
+ assert_array_almost_equal(lut(x,y),z)
+
+if __name__ == "__main__":
+ NumpyTest().run()
More information about the Scipy-svn
mailing list