[Scipy-svn] r4595 - in branches/Interpolate1D: . docs tests

scipy-svn@scip... scipy-svn@scip...
Fri Aug 1 18:34:04 CDT 2008


Author: fcady
Date: 2008-08-01 18:34:03 -0500 (Fri, 01 Aug 2008)
New Revision: 4595

Added:
   branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
Modified:
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/fitpack_wrapper2d.py
   branches/Interpolate1D/tests/test_fitpack_wrapper.py
Log:
more work and documentation for 2D interpolation

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-01 23:34:03 UTC (rev 4595)
@@ -471,14 +471,12 @@
 
 At instantiation:
 
-#) bbox
+#. bbox
 This is a 2-element list specifying the endpoints of the approximation interval.
 It default to [x[0],x[-1]]
-    
-#) w
+#. w
 a 1D sequence of weights which defaults to all ones.
-
-#) s 
+#. s 
 If s is zero, the interpolation is exact.  If s is not 0, the curve is smoothe subject to
 the constraint that sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
     
@@ -487,27 +485,30 @@
 
 At calling:
 
-#) nu
+#. nu
 Spline returns, not the spline function S, but the (nu)th derivative of S.  nu defaults
 to 0, so Spline usually returns the zeroth derivative of S, ie S.
 
+
 -----------------
 Special Methods
 -----------------
 
-#) set_smoothing_factor(s)
-#) get_knots
+The following special methods are also available, which are not wrapped by Interpolate1d.
+
+#. set_smoothing_factor(s = 0.0)
+#. get_knots
 returns the positions of the knots of the spline
-#) get_coeffs
+#. get_coeffs
 returns the coefficients of the 
-#) get_residual
+#. get_residual
 returns the weighted sum of the errors (due to smoothing) at the data points
 sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0)
-#) integral(a, b)
+#. integral(a, b)
 returns the integral from a to b
-#) derivatives(x)
+#. derivatives(x)
 returns all the derivatives of the spline at point x
-#) roots
+#. roots
 This only works for cubic splines.  But it returns the places where the spline
 is identically zero.
 
@@ -550,8 +551,43 @@
 The Objective Interface
 ------------------------------------------
 
-The objective interface for 2D is slightly more complicated.  
+The objective interface for 2D is similarly analogous to 1D.
 
+------------------------------------------
+The Spline2d Class
+------------------------------------------
+
+Just as with Spline, Spline2d is mostly intended to be wrapper by
+Interpolate2d, but it also functions as a stand-alone class featuring
+functionality not accessible through Interpolate2d.
+
+It is instantiated in virtually the same way ::
+
+    instance = Spline2d(x, y, z, kx = 3, ky = 3, s=0.0)
+    
+where x, y and z are 1D arrays.  It is called with arrays
+newx and newy, returning an array newz of the same length.
+
+Beyond basic usage, Spline2 also has the methods
+
+#) get_grid(self, x, y)
+x and y are treated as the coordinates of a grid, and all
+points on the grid are interpolated and returned in an array.
+    
+That is, if z = S.get_grid(x,y), z[i,j] is the interpolated value
+at the point (xi, yj)
+
+#) integral(xa, xb, ya, yb)
+Integrate the interpolated function over the indicated rectangle
+
+#) get_residual
+
+#) get_knots
+
+#) get_coeffs
+    
+
+
 ================================================
 ND Interpolation
 ================================================

Modified: branches/Interpolate1D/fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-01 23:34:03 UTC (rev 4595)
@@ -50,14 +50,11 @@
         [xb,xe] x [yb, ye] calculated from a given set of data points
         (x,y,z).
 
-        See also:
-
-        bisplrep, bisplev - an older wrapping of FITPACK
-        UnivariateSpline - a similar class for univariate spline interpolation
-        SmoothUnivariateSpline - to create a BivariateSpline through the
-                                 given points
-        LSQUnivariateSpline - to create a BivariateSpline using weighted
-                              least-squares fitting
+        More commenting needed
+        
+        If (xi, yi) is outside the interpolation range, it is
+        assigned the value of the nearest point that is within the
+        interpolation range.
     """
     def __init__(self, x=None, y=None, z=None, w=None, bbox=[None]*4, kx=3, ky=3, s=0.0, eps=None):
         """
@@ -118,8 +115,11 @@
     def __call__(self, x, y):
         """ Evaluate spline at positions x[i],y[i].
             x and y should be 1d arrays.
+            
+            If (xi, yi) is outside the interpolation range, it will be
+            assigned the value of the nearest point that is within the
+            interpolation range.
         """
-        # what happens when x contains duplicate values?
         
         if self._is_initialized is not True:
             raise Error, "x, y and z must be initialized before interpolating"
@@ -131,26 +131,24 @@
         data_grid = self.get_grid(sorted_x, sorted_y)
         
         # fixme : no list comprehension
-        z = np.array([ data_grid[np.searchsorted(sorted(x), x[i]), np.searchsorted(sorted(y),y[i])] \
+        z = np.array([ data_grid[np.searchsorted(sorted_x, x[i]), np.searchsorted(sorted_y,y[i])] \
                                     for i,xi in enumerate(x) ])
-            
+        
         return z
         
         
-    def get_grid(self, x, y, mth='array'):
+    def get_grid(self, x, y):
         """ Evaluate spline at positions x[i],y[j]."""
         
         if self._is_initialized is not True:
             raise Error, "x, y and z must be initialized before interpolating"
         
-        if mth=='array':
-            tx,ty,c = self.tck[:3]
-            kx,ky = self.degrees
-            z,ier = _dfitpack.bispev(tx,ty,c,kx,ky,x,y)
-            assert ier==0,'Invalid input: ier='+`ier`
-            return z
-        raise NotImplementedError
-
+        tx,ty,c = self.tck[:3]
+        kx,ky = self.degrees
+        z,ier = _dfitpack.bispev(tx,ty,c,kx,ky,x,y)
+        assert ier==0,'Invalid input: ier='+`ier`
+        return z
+        
     def get_residual(self):
         """ Return weighted sum of squared residuals of the spline
         approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0)

Modified: branches/Interpolate1D/tests/test_fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper.py	2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper.py	2008-08-01 23:34:03 UTC (rev 4595)
@@ -17,7 +17,7 @@
     def assertAllclose(self, x, y):
         self.assert_(np.allclose(x, y))
         
-    def test_linearInterp(self):
+    def test_k_1(self):
         """ make sure : linear interpolation (spline with order = 1, s = 0)works
         """
         N = 3000.
@@ -34,7 +34,7 @@
         #print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
         self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
         
-    def test_quadInterp(self):
+    def test_k_2(self):
         """ make sure : quadratic interpolation (spline with order = 2, s = 0)works
         """
         N = 3000.
@@ -49,7 +49,19 @@
         #print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
         self.assertAllclose(new_y, y)
         
+    def test_extrap(self):
+        """ make sure 1D extrapolation works
+        """
+        N = 3000.
+        x = np.arange(N)
+        y = np.arange(N)
         
+        interp_func = Spline(x, y, k=1)
+        newx = np.arange(-2, N)+0.5
+        newy = interp_func(newx)
+        
+        self.assertAllclose(newx, newy)
+        
     def test_inputFormat(self):
         """ make sure : it's possible to instantiate Spline without x and y
         """

Added: branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper2d.py	2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper2d.py	2008-08-01 23:34:03 UTC (rev 4595)
@@ -0,0 +1,106 @@
+""" test fitpack_wrapper2d.py
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+# testing
+import unittest
+import time
+from numpy import arange, allclose, ones, meshgrid, ravel, array
+import numpy as np
+from fitpack_wrapper2d import Spline2d
+
+class Test(unittest.TestCase):
+    
+    def assertAllclose(self, x, y):
+        self.assert_(np.allclose(x, y))
+        
+    def test_k_1(self):
+        """ make sure : linear interpolation (spline with order = 1, s = 0)works
+        """
+        N = 10.
+        X, Y = meshgrid( arange(N), arange(N) )
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z])
+        
+        newx = arange(N-1) +.5
+        newy = newx
+        
+        interp_func = Spline2d(x, y, z, kx=1, ky=1)
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_k_2(self):
+        """ make sure : quadratic interpolation (spline with order = 2, s = 0)works
+        """
+        N = 10.
+        X, Y = meshgrid( arange(N), arange(N) )
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z])
+        
+        newx = arange(N-1)+0.5
+        newy = newx
+        
+        interp_func = Spline2d(x, y, z, kx=2, ky=2)
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+        
+    def test_instantiate_without_init(self):
+        """ make sure : it's possible to instantiate Spline2d without setting x and y
+        """
+        N = 10.
+        X, Y = meshgrid( arange(N), arange(N) )
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z])
+        
+        newx = np.arange(N-1)+0.5
+        newy = newx
+        
+        interp_func = Spline2d(kx=1, ky=3)
+        interp_func.init_xyz(x, y, z)
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_extrap(self):
+        """ make sure : linear extrapolation works
+        """
+        N = 10.
+        X, Y = meshgrid( arange(N), arange(N) )
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z])
+        
+        newx = arange(N+8) +.5
+        newy = 2*newx
+        
+        interp_func = Spline2d(x, y, z, kx=1, ky=1)
+        newz = interp_func(newx, newy)
+        
+        print "newx: ", newx
+        print "newy: ", newy
+        print "sum : ", newx+newy
+        print "newz: ", newz
+        
+        print "Homer Simpson"
+        print interp_func(array([-2.0]),array([3.5]))
+        print interp_func(array([-7.0]),array([3.5]))
+        print interp_func(array([-2.0]),array([7]))
+        print "Bartman"
+        
+        self.assertAllclose(newz, newx+newy)
+    
+    def runTest(self):
+        test_list = [name for name in dir(self) if name.find('test_')==0]
+        for test_name in test_list:
+            exec("self.%s()" % test_name)
+        
+        
+if __name__ == '__main__':
+    unittest.main()
+    
+    
\ No newline at end of file



More information about the Scipy-svn mailing list