[Scipy-svn] r4467 - in trunk/scipy/interpolate: . tests
scipy-svn@scip...
scipy-svn@scip...
Mon Jun 23 16:14:21 CDT 2008
Author: ptvirtan
Date: 2008-06-23 16:14:11 -0500 (Mon, 23 Jun 2008)
New Revision: 4467
Modified:
trunk/scipy/interpolate/fitpack.py
trunk/scipy/interpolate/fitpack.pyf
trunk/scipy/interpolate/fitpack2.py
trunk/scipy/interpolate/tests/test_fitpack.py
Log:
Wrap and expose dblint from dfitpack. (Implements #206). Add corresponding tests.
Modified: trunk/scipy/interpolate/fitpack.py
===================================================================
--- trunk/scipy/interpolate/fitpack.py 2008-06-23 20:53:22 UTC (rev 4466)
+++ trunk/scipy/interpolate/fitpack.py 2008-06-23 21:14:11 UTC (rev 4467)
@@ -842,6 +842,28 @@
if len(z[0])>1: return z[0]
return z[0][0]
+def dblint(xa,xb,ya,yb,tck):
+ """Evaluate the integral of a spline over area [xa,xb] x [ya,yb].
+
+ Parameters
+ ----------
+ xa, xb : float
+ The end-points of the x integration interval.
+ ya, yb : float
+ The end-points of the y integration interval.
+ tck : list [tx, ty, c, kx, ky]
+ A sequence of length 5 returned by bisplrep containing the knot
+ locations tx, ty, the coefficients c, and the degrees kx, ky
+ of the spline.
+
+ Returns
+ -------
+ integ : float
+ The value of the resulting integral.
+ """
+ tx,ty,c,kx,ky=tck
+ return dfitpack.dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
+
def insert(x,tck,m=1,per=0):
"""Insert knots into a B-spline.
Modified: trunk/scipy/interpolate/fitpack.pyf
===================================================================
--- trunk/scipy/interpolate/fitpack.pyf 2008-06-23 20:53:22 UTC (rev 4466)
+++ trunk/scipy/interpolate/fitpack.pyf 2008-06-23 21:14:11 UTC (rev 4467)
@@ -456,7 +456,24 @@
:: kwrk=3+mx+my+nxest+nyest
integer intent(out) :: ier
end subroutine regrid_smth
-
+
+ function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
+ ! iy = dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
+ 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
+ real*8 intent(in) :: xb
+ real*8 intent(in) :: xe
+ real*8 intent(in) :: yb
+ real*8 intent(in) :: ye
+ real*8 dimension(nx+ny-kx-ky-2),depend(nx,ny,kx,ky),intent(cache,hide) :: wrk
+ real*8 :: dblint
+ end function dblint
end interface
end python module dfitpack
Modified: trunk/scipy/interpolate/fitpack2.py
===================================================================
--- trunk/scipy/interpolate/fitpack2.py 2008-06-23 20:53:22 UTC (rev 4466)
+++ trunk/scipy/interpolate/fitpack2.py 2008-06-23 21:14:11 UTC (rev 4467)
@@ -352,6 +352,27 @@
assert ier==0,'Invalid input: ier='+`ier`
return z
raise NotImplementedError
+
+ def integral(self, xa, xb, ya, yb):
+ """
+ Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
+
+ Parameters
+ ----------
+ xa, xb : float
+ The end-points of the x integration interval.
+ ya, yb : float
+ The end-points of the y integration interval.
+
+ Returns
+ -------
+ integ : float
+ The value of the resulting integral.
+
+ """
+ tx,ty,c = self.tck[:3]
+ kx,ky = self.degrees
+ return dfitpack.dblint(tx,ty,c,kx,ky,xa,xb,ya,yb)
class SmoothBivariateSpline(BivariateSpline):
""" Smooth bivariate spline approximation.
Modified: trunk/scipy/interpolate/tests/test_fitpack.py
===================================================================
--- trunk/scipy/interpolate/tests/test_fitpack.py 2008-06-23 20:53:22 UTC (rev 4466)
+++ trunk/scipy/interpolate/tests/test_fitpack.py 2008-06-23 21:14:11 UTC (rev 4467)
@@ -14,7 +14,7 @@
import sys
from scipy.testing import *
-from numpy import array
+from numpy import array, diff
from scipy.interpolate.fitpack2 import UnivariateSpline,LSQUnivariateSpline,\
InterpolatedUnivariateSpline
from scipy.interpolate.fitpack2 import LSQBivariateSpline, \
@@ -48,10 +48,49 @@
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()
+ assert_almost_equal(lut(2,2), 3.)
+
+ def test_bilinearity(self):
+ x = [1,1,1,2,2,2,3,3,3]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = [0,7,8,3,4,7,1,3,4]
+ s = 0.1
+ tx = [1+s,3-s]
+ ty = [1+s,3-s]
+ lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1)
+
+ tx, ty = lut.get_knots()
+
+ for xa, xb in zip(tx[:-1], tx[1:]):
+ for ya, yb in zip(ty[:-1], ty[1:]):
+ for t in [0.1, 0.5, 0.9]:
+ for s in [0.3, 0.4, 0.7]:
+ xp = xa*(1-t) + xb*t
+ yp = ya*(1-s) + yb*s
+ zp = (+ lut(xa, ya)*(1-t)*(1-s)
+ + lut(xb, ya)*t*(1-s)
+ + lut(xa, yb)*(1-t)*s
+ + lut(xb, yb)*t*s)
+ assert_almost_equal(lut(xp,yp), zp)
+
+ def test_integral(self):
+ x = [1,1,1,2,2,2,8,8,8]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = array([0,7,8,3,4,7,1,3,4])
+
+ s = 0.1
+ tx = [1+s,3-s]
+ ty = [1+s,3-s]
+ lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1)
+ tx, ty = lut.get_knots()
+
+ tz = lut(tx, ty)
+ trpz = .25*(diff(tx)[:,None]*diff(ty)[None,:]
+ *(tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum()
+
+ assert_almost_equal(lut.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz)
+
class TestSmoothBivariateSpline(TestCase):
def test_linear_constant(self):
x = [1,1,1,2,2,2,3,3,3]
@@ -73,6 +112,29 @@
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]])
+ def test_integral(self):
+ x = [1,1,1,2,2,2,4,4,4]
+ y = [1,2,3,1,2,3,1,2,3]
+ z = array([0,7,8,3,4,7,1,3,4])
+
+ lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1,s=0)
+ tx = [1,2,4]
+ ty = [1,2,3]
+
+ tz = lut(tx, ty)
+ trpz = .25*(diff(tx)[:,None]*diff(ty)[None,:]
+ *(tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum()
+ assert_almost_equal(lut.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz)
+
+ lut2 = SmoothBivariateSpline(x,y,z,kx=2,ky=2,s=0)
+ assert_almost_equal(lut2.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz,
+ decimal=0) # the quadratures give 23.75 and 23.85
+
+ tz = lut(tx[:-1], ty[:-1])
+ trpz = .25*(diff(tx[:-1])[:,None]*diff(ty[:-1])[None,:]
+ *(tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum()
+ assert_almost_equal(lut.integral(tx[0], tx[-2], ty[0], ty[-2]), trpz)
+
class TestRectBivariateSpline(TestCase):
def test_defaults(self):
x = array([1,2,3,4,5])
More information about the Scipy-svn
mailing list