[Scipy-svn] r2254 - in trunk/Lib/interpolate: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Oct 11 19:29:20 EDT 2006
Author: stefan
Date: 2006-10-11 18:28:59 -0500 (Wed, 11 Oct 2006)
New Revision: 2254
Modified:
trunk/Lib/interpolate/fitpack.pyf
trunk/Lib/interpolate/fitpack2.py
trunk/Lib/interpolate/tests/test_fitpack.py
Log:
Add RectBivariateSpline [for John Travers]. Fix long lines.
Modified: trunk/Lib/interpolate/fitpack.pyf
===================================================================
--- trunk/Lib/interpolate/fitpack.pyf 2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/fitpack.pyf 2006-10-11 23:28:59 UTC (rev 2254)
@@ -63,6 +63,12 @@
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,
+ 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;
+}
'''
interface
@@ -411,7 +417,46 @@
:: 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,&
+ 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])
+ 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)
+ 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)
+ real*8 optional,depend(x,mx) :: xe=dmax(x,mx)
+ real*8 optional,depend(y,my) :: yb=dmin(y,my)
+ real*8 optional,depend(y,my) :: ye=dmax(y,my)
+ integer optional,check(1<=kx && kx<=5) :: kx = 3
+ integer optional,check(1<=ky && ky<=5) :: ky = 3
+ real*8 optional,check(0.0<=s) :: s = 0.0
+ integer intent(hide),depend(kx,mx),check(nxest>=2*(kx+1)) &
+ :: nxest = mx+kx+1
+ integer intent(hide),depend(ky,my),check(nyest>=2*(ky+1)) &
+ :: nyest = my+ky+1
+ integer intent(out) :: nx
+ real*8 dimension(nxest),intent(out),depend(nxest) :: tx
+ integer intent(out) :: ny
+ real*8 dimension(nyest),intent(out),depend(nyest) :: 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(lwrk),intent(cache,hide),depend(lwrk) :: wrk
+ integer intent(hide),depend(mx,my,kx,ky,nxest,nyest) &
+ :: lwrk=calc_regrid_lwrk(mx,my,kx,ky,nxest,nyest)
+ integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
+ integer intent(hide),depend(mx,my,nxest,nyest) &
+ :: kwrk=3+mx+my+nxest+nyest
+ integer intent(out) :: ier
+ end subroutine regrid_smth
+
end interface
end python module dfitpack
Modified: trunk/Lib/interpolate/fitpack2.py
===================================================================
--- trunk/Lib/interpolate/fitpack2.py 2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/fitpack2.py 2006-10-11 23:28:59 UTC (rev 2254)
@@ -13,10 +13,11 @@
'LSQUnivariateSpline',
'LSQBivariateSpline',
- 'SmoothBivariateSpline']
+ 'SmoothBivariateSpline',
+ 'RectBivariateSpline']
import warnings
-from numpy import zeros, concatenate, alltrue
+from numpy import zeros, concatenate, alltrue, ravel, all, diff
import dfitpack
@@ -134,6 +135,7 @@
def set_smoothing_factor(self, s):
""" Continue spline computation with the given smoothing
factor s and with the knots found at the last call.
+
"""
data = self._data
if data[6]==-1:
@@ -150,8 +152,9 @@
def __call__(self, x, nu=None):
""" Evaluate spline (or its nu-th derivative) at positions x.
- Note: x can be unordered but the evaluation is
- more efficient if x is (partially) ordered.
+ Note: x can be unordered but the evaluation is more efficient
+ if x is (partially) ordered.
+
"""
if nu is None:
return dfitpack.splev(*(self._eval_args+(x,)))
@@ -166,14 +169,15 @@
return data[8][k:n-k]
def get_coeffs(self):
- """ Return spline coefficients."""
+ """Return spline coefficients."""
data = self._data
k,n = data[5],data[7]
return data[9][:n-k-1]
def get_residual(self):
- """ Return weighted sum of squared residuals of the spline
+ """Return weighted sum of squared residuals of the spline
approximation: sum ((w[i]*(y[i]-s(x[i])))**2,axis=0)
+
"""
return self._data[10]
@@ -199,11 +203,15 @@
z,m,ier = dfitpack.sproot(*self._eval_args[:2])
assert ier==0,`ier`
return z[:m]
- raise NotImplementedError,'finding roots unsupported for non-cubic splines'
+ raise NotImplementedError,\
+ 'finding roots unsupported for non-cubic splines'
class InterpolatedUnivariateSpline(UnivariateSpline):
- """ Interpolated univariate spline approximation. Identical to UnivariateSpline with less error checking."""
+ """ Interpolated univariate spline approximation. Identical to
+ UnivariateSpline with less error checking.
+ """
+
def __init__(self, x, y, w=None, bbox = [None]*2, k=3):
"""
Input:
@@ -223,8 +231,12 @@
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. Appears to be identical to UnivariateSpline with
+ more error checking.
+ """
+
def __init__(self, x, y, t, w=None, bbox = [None]*2, k=3):
"""
Input:
@@ -361,7 +373,8 @@
w - positive 1-d sequence of weights
bbox - 4-sequence specifying the boundary of
the rectangular approximation domain.
- By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
+ By default, bbox=[min(x,tx),max(x,tx),
+ min(y,ty),max(y,ty)]
kx,ky=3,3 - degrees of the bivariate spline.
s - positive smoothing factor defined for
estimation condition:
@@ -374,9 +387,10 @@
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,\
- xb,xe,yb,ye,\
- kx,ky,s=s,eps=eps,lwrk2=1)
+ nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth(x,y,z,w,
+ xb,xe,yb,ye,
+ kx,ky,s=s,
+ eps=eps,lwrk2=1)
if ier in [0,-1,-2]: # normal return
pass
else:
@@ -410,7 +424,8 @@
w - positive 1-d sequence of weights
bbox - 4-sequence specifying the boundary of
the rectangular approximation domain.
- By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
+ By default, bbox=[min(x,tx),max(x,tx),
+ min(y,ty),max(y,ty)]
kx,ky=3,3 - degrees of the bivariate spline.
eps - a threshold for determining the effective rank
of an over-determined linear system of
@@ -443,3 +458,63 @@
self.fp = fp
self.tck = tx1,ty1,c
self.degrees = kx,ky
+
+class RectBivariateSpline(BivariateSpline):
+ """ Bivariate spline approximation over a rectangular mesh.
+
+ Can be used for both smoothing or interpolating data.
+
+ See also:
+
+ SmoothBivariateSpline - a smoothing bivariate spline for scattered data
+ 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):
+ """
+ Input:
+ x,y - 1-d sequences of coordinates in strictly ascending order
+ z - 2-d array of data with shape (x.size,y.size)
+ Optional input:
+ bbox - 4-sequence specifying the boundary of
+ the rectangular approximation domain.
+ By default, bbox=[min(x,tx),max(x,tx),
+ min(y,ty),max(y,ty)]
+ kx,ky=3,3 - degrees of the bivariate spline.
+ s - positive smoothing factor defined for
+ estimation condition:
+ sum((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0) <= s
+ Default s=0 which is for interpolation
+ """
+ x,y = ravel(x),ravel(y)
+ if not all(diff(x) > 0.0):
+ raise TypeError,'x must be strictly increasing'
+ if not all(diff(y) > 0.0):
+ raise TypeError,'y must be strictly increasing'
+ if not ((x.min() == x[0]) and (x.max() == x[-1])):
+ raise TypeError, 'x must be strictly ascending'
+ if not ((y.min() == y[0]) and (y.max() == y[-1])):
+ raise TypeError, 'y must be strictly ascending'
+ if not x.size == z.shape[0]:
+ raise TypeError,\
+ 'x dimension of z must have same number of elements as x'
+ if not y.size == z.shape[1]:
+ raise TypeError,\
+ '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,
+ xb,xe,yb,ye,
+ kx,ky,s)
+ if ier in [0,-1,-2]: # normal return
+ pass
+ else:
+ message = _surfit_messages.get(ier,'ier=%s' % (ier))
+ warnings.warn(message)
+
+ self.fp = fp
+ self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
+ self.degrees = kx,ky
+
Modified: trunk/Lib/interpolate/tests/test_fitpack.py
===================================================================
--- trunk/Lib/interpolate/tests/test_fitpack.py 2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/tests/test_fitpack.py 2006-10-11 23:28:59 UTC (rev 2254)
@@ -14,10 +14,12 @@
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
+from interpolate.fitpack2 import LSQBivariateSpline, SmoothBivariateSpline,\
+ RectBivariateSpline
restore_path()
class test_UnivariateSpline(ScipyTestCase):
@@ -73,5 +75,13 @@
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(ScipyTestCase):
+ 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__":
ScipyTest().run()
More information about the Scipy-svn
mailing list