[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