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

scipy-svn@scip... scipy-svn@scip...
Thu Aug 7 14:31:47 CDT 2008


Author: fcady
Date: 2008-08-07 14:31:44 -0500 (Thu, 07 Aug 2008)
New Revision: 4608

Modified:
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/ndimage_wrapper.py
   branches/Interpolate1D/tests/test_ndimage.py
Log:
more exhaustive tests (which resulted in fixing a few bugs that had been there) and more-or-less complete documentation

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-07 19:31:44 UTC (rev 4608)
@@ -645,8 +645,57 @@
 ND Interpolation
 ================================================
 
+1D and 2D interpolation are analogous in their user interface.  However,
+the ND interpolation breaks from their patterns in certain important respects.
 
+The biggest difference is that the known data must be on a uniformly spaced
+grid; scattered points are not acceptable.  Rather than a list of points, the
+user passes in an N-dimensional array of data.  By default, the (i, j, k) element 
+of the array is taken to give the value of the function at the point (i, j, k), but
+the user can specify a different set of starting points and spacings.  The
+starting points and spacings can be different for each axis.  Starting points
+and spacings are specified by the keywords 'starting_coords' and 'spacings',
+which can be passed as arrays or lists.
 
+Second, when calling the object / function, the points to be interpolated are
+passed as an nxL array, where n is the dimensionality of the data; each column
+of the array specifies a point at which to interpolate a value, and the returned
+array of length L gives the value at each of the points.
+
+Here is a basic example which illustrates these points: ::
+
+    In []: from interpolate import interpNd
+    In []: X, Y = meshgrid(arange(10.), arange(10.))
+    In []: Z = X+Y # function just adds coordinates
+    In []: coordinates = array([ [1.1, 2.2, 4.6],
+                                          [1.1, 2.2, 4.0 ])
+    In []: interpNd(Z, coordinates)
+    Out []: array([2.2, 4.4, 8.6])
+    # say the data start at the point (2, 1) and
+    # points are spaced 2 apart
+    In []: interpNd(Z, coordinates, starting_coords = array([1, 1], spacings=[2, 2])
+    Out []: array([ .1, 1.2, 3.3])
+    
+By default, the interpolation is linear for points in-range and returns NaN for
+out-of-bounds; alternate types can be specified by the keywords
+'kind' and 'out', as in 2D interpolation.  However, 
+
+#) 'kind' must be a string ('linear', 'block', 'cubic', etc) indicating a type of
+    spline interpolation, or else an integers specifying the spline order.
+#) 'out' must be either NaN (the default), 'nearest', 'wrap', 'reflect' or 'constant'
+
+The user cannot pass in specially-tailored interpolation methods.
+
+There 
+
+The second point is that all interpolation is done using splines.  The keyword is still
+"kind", but only keywords specifying splines ('spline', 'cubic', 'quadratic', 'quintic', etc)
+are acceptable.  If kind is an integer, that integer is taken to be the order of the spline.
+Note also that 'linear' denotes a spline of order 1, and 'block' denotes a spline of order
+zero.  
+
+
+
 ================================================
 ND Scattered Interpolation
 ================================================

Modified: branches/Interpolate1D/ndimage_wrapper.py
===================================================================
--- branches/Interpolate1D/ndimage_wrapper.py	2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/ndimage_wrapper.py	2008-08-07 19:31:44 UTC (rev 4608)
@@ -4,21 +4,28 @@
 import numpy as np
 import _nd_image
 
+def interpNd(data, coordinates, starting_coords=None, spacings=None, kind='linear',out=NaN):
+    return Interpolate1d(data = data,
+                                    starting_coords = starting_coords,
+                                    spacings = spacings,
+                                    kind = kind,
+                                    out = out
+                                    )(coordinates)
+
 class InterpolateNd:
     def __init__(self, data, starting_coords =None, spacings = None, 
-                        order=3, out=NaN):
+                        kind='linear', out=NaN):
         """ data = array or list of lists
             starting_coords = None, list, 1D array or 2D (nx1) array
             spacings = None, list, 1D array or 2D (nx1) array
+            kind = string or integer
+                0 = block extrapolation between midpoints
             out = string in 'nearest', 'wrap', 'reflect', 'mirror', 'constant'
                         or just NaN
         """
         
         # FIXME : include spline filtering
         
-        if order < 0 or order > 5:
-            raise RuntimeError, 'spline order not supported'
-        
         # checking format of input
         data = array(data)
         
@@ -29,15 +36,54 @@
             starting_coords = array(starting_coords)
             assert starting_coords.size == data.ndim, "There must be one element of \
                             starting_coords per data dimension.  Size mismatch."
-            starting_coords = reshape(starting_coords, (data.ndim, 1))
+            starting_coords = np.reshape(starting_coords, (data.ndim, 1))
         if spacings == None:
-            spacings = np.zeros(( data.ndim, 1 ))
+            spacings = np.ones(( data.ndim, 1 ))
         else:
             spacings = array(spacings)
             assert starting_coords.size == data.ndim, "There must be one element of \
                             starting_coords per data dimension"
-            spacings = reshape(spacings, (data.ndim, 1))
+            spacings = np.reshape(spacings, (data.ndim, 1))
         
+        # determining the order
+        order_dict = \
+            { 0:0,
+                '0':0,
+                'block':0,
+                1:1,
+                '1':1,
+                'linear':1,
+                'Linear':1,
+                2:2,
+                '2':2,
+                'quadratic':2,
+                'quad':2,
+                'Quadratic':2,
+                'Quad':2,
+                3:3,
+                '3':3,
+                'spline':3,
+                'Spline':3,
+                'cubic':3,
+                'Cubic':3,
+                4:4,
+                '4':4,
+                'quartic':4,
+                'Quartic':4,
+                5:5,
+                '5':5,
+                'quintic':5,
+                'quint':5,
+                'Quintic':5,
+                'Quint':5
+                }
+        if order_dict.has_key(kind):
+            self.order = order_dict[kind]
+        elif isinstance(kind, int):
+            raise ValueError, "Only spline orders 0, 1, ..., 5 are supported"
+        else:
+            raise ValueError, "argument kind = %s not recognized" % str(kind)
+                
         # storing relevant data
         self._data_array = data
         self.ndim = data.ndim
@@ -46,7 +92,6 @@
         self._min_coords = starting_coords
         self._max_coords = self._min_coords + self._shape*self._spacings
         self.out = out
-        self.order = order
         
     def __call__(self, coordinates):
         """ coordinates is an n x L array, where n is the dimensionality of the data

Modified: branches/Interpolate1D/tests/test_ndimage.py
===================================================================
--- branches/Interpolate1D/tests/test_ndimage.py	2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/tests/test_ndimage.py	2008-08-07 19:31:44 UTC (rev 4608)
@@ -7,7 +7,7 @@
 
 import unittest
 import time
-from numpy import arange, allclose, ones
+from numpy import arange, allclose, ones, array
 import numpy as np
 import ndimage_wrapper as nd
 
@@ -17,40 +17,97 @@
         self.assert_(np.allclose(x, y))
     
     def test_linear(self):
+        """ Make sure : basic linear works
+        """
         boring_data = np.ones((5,5,5))
-        interp = nd.InterpolateNd(boring_data, order = 1)
+        interp = nd.InterpolateNd(boring_data, kind = 'linear')
         self.assertAllclose( interp(np.array([[2.3], [1.0], [3.9]])) , 1.0 )
         
+    def test_linear_not_1(self):
+        """ Make sure : linear interpolation works on a general dataset
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, kind = 'linear')
+        self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 )
+        
     def test_data_is_list(self):
+        """ Make sure : data can be entered as a list
+        """
         boring_data = [ [1.0, 1.0, 1.0],
                               [1.0, 1.0, 1.0],
                               [1.0, 1.0, 1.0]]
-        interp = nd.InterpolateNd(boring_data, order = 1)
+        interp = nd.InterpolateNd(boring_data)
         self.assertAllclose( interp(np.array([[1.3], [1.0]])) , 1.0 )
         
     def test_coords_is_1d(self):
+        """ Make sure : coordinates for a single point can be entered as a 1D array
+        """
         boring_data = np.ones((5,5,5))
-        interp = nd.InterpolateNd(boring_data, order = 1)
+        interp = nd.InterpolateNd(boring_data)
         self.assertAllclose( interp(np.array([2.3, 1.0, 3.9])) , 1.0 )
         
     def test_coords_is_list(self):
+        """ Make sure : coordinates for a single point can be entered as a list
+        """
         boring_data = np.ones((5,5,5))
-        interp = nd.InterpolateNd(boring_data, order = 1)
+        interp = nd.InterpolateNd(boring_data)
         self.assertAllclose( interp([2.3, 1.0, 3.9]) , 1.0 )
+            
+    def test_order2(self):
+        """ Make sure : quadratic interpolation works
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, kind = 2)
+        self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 )
         
+    def test_order0(self):
+        """ Make sure : block interpolation works
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, kind = 0)
+        self.assertAllclose( interp(np.array([[2.3], [1.1]])) , 3.0 )
+        
+    def test_order3(self):
+        """ Make sure : cubic interpolation works
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, kind = 3)
+        self.assertAllclose( interp(np.array([[4.3], [4.1]])) , 8.4 )
+        
+    def test_out(self):
+        """ Make sure : out-of-bounds returns NaN
+        """
+        boring_data = np.ones((5,5,5))
+        interp = nd.InterpolateNd(boring_data, kind = 'linear')
+        self.assert_( np.isnan(interp(  np.array([[7.3], [1.0], [3.9]])  )))
+        
+    def test_starting_coords(self):
+        """ Make sure : non-zero starting coordinates work correctly
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, starting_coords = array([2, 1]))
+        self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 0.3 )
+    
+    def test_spacings(self):
+        """ Make sure : spacings other than 1 work correctly
+        """
+        X, Y = np.meshgrid(arange(10.), arange(10.))
+        interesting_data = X+Y
+        interp = nd.InterpolateNd(interesting_data, spacings = array([2, 1]))
+        self.assertAllclose( interp(np.array([[2.4], [1.0]])) , 2.2 )
+        
     def runTest(self):
+        """ run all tests
+        """
         test_list = [method_name for method_name in dir(self) if method_name.find('test')==0]
         for test_name in test_list:
             exec("self.%s()" % test_name)
-            
-    def test_order2(self):
-        boring_data = np.ones((5,5,5))
-        interp = nd.InterpolateNd(boring_data, order = 2)
-        self.assertAllclose( interp(np.array([[2.3], [1.0], [3.9]])) , 1.0 )
         
-    def test_out(self):
-        pass
         
-        
 if __name__ == '__main__':
     unittest.main()
\ No newline at end of file



More information about the Scipy-svn mailing list