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

scipy-svn@scip... scipy-svn@scip...
Tue Aug 19 14:02:05 CDT 2008


Author: fcady
Date: 2008-08-19 14:02:03 -0500 (Tue, 19 Aug 2008)
New Revision: 4652

Modified:
   branches/Interpolate1D/__init__.py
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/interpolate1d.py
   branches/Interpolate1D/interpolate2d.py
   branches/Interpolate1D/setup.py
   branches/Interpolate1D/tests/test_interpolate2d.py
Log:
minor changes, some in preparation for ND scattered interpolation

Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/__init__.py	2008-08-19 19:02:03 UTC (rev 4652)
@@ -18,3 +18,5 @@
 # wrapper around Fortran implementing Tom's algorithm 526
 from algorithm526_wrapper import algorithm526
 
+from interpSNd import InterpolateSNd
+

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-19 19:02:03 UTC (rev 4652)
@@ -166,9 +166,9 @@
 interpolated.
 
 
---------------------------------------
+-----------------------------------------------------------
 User-defined Interpolation Methods
---------------------------------------
+-----------------------------------------------------------
 
 The string interface is designed to conveniently take care of most things a user would want
 to do in a way that is easy and, when something goes wrong, informative and helpful.
@@ -579,8 +579,8 @@
 
     newz = interp2d(x, y, z, newx, newy, kind='linear', out=NaN)
 
-where x, y, z, are 1D arrays or lists, and newx and newy may be either arrays, lists or scalars.  
-If they are scalars or zero-dimensional arrays, newz will be a scalar as well.  Otherwise
+where x, y, z, are arrays (1D or 2D) or lists, and newx and newy may be either arrays, lists or scalars.  
+If newx and newy are scalars or zero-dimensional arrays, newz will be a scalar as well.  Otherwise
 a vector is returned.  The only differences from intper1d are
 
 #) The known data points are specified by 3 arrays (x, y and z) rather than 2 (x and y).
@@ -737,8 +737,22 @@
 ND Scattered Interpolation
 ================================================
  
- Still in development.
+Still in development.
  
- Ideally the range of interpolation would be the convex hull of the known
- data points, and a Delaunay tesselation would be determined and stored
- at instantiation.  Then again, that would be very expensive.
\ No newline at end of file
+Ideally the range of interpolation would be the convex hull of the known
+data points, and a Delaunay tesselation would be determined and stored
+at instantiation.  Then again, that would be very expensive.
+ 
+InterpolateNd suffers from the requirement that data points be on a uniformly
+spaced grid.  This problem is solved by the callable class InterpolateSNd and
+the function interpSNd, which interpolate *scattered* N-dimensional data.
+
+First off a warning.  The code is in a rather preliminary stage which, while functioning,
+is VERY slow and not extensively tested.  This state of affairs is still being improved.
+
+The valid interpolation range is the convex hull of the data points.
+
+
+ 
+ 
+This functionality is still preliminary.
\ No newline at end of file

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/interpolate1d.py	2008-08-19 19:02:03 UTC (rev 4652)
@@ -164,7 +164,7 @@
                 interpolate from.  Note that 2-dimensional
                 y is not supported.
             
-            newx -- list of 1D numpy array
+            newx (when calling) -- list of 1D numpy array
                 x values at which to interpolate the value
                 of the function
                 

Modified: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/interpolate2d.py	2008-08-19 19:02:03 UTC (rev 4652)
@@ -36,26 +36,137 @@
 
 # functional interface: creates and calls an instance of objective interface
 def interp2d(x, y, z, newx, newy, kind='linear', out=NaN, bad_data=None):
+    """ A function for interpolation of 2D, real-valued data.
+        
+        Parameters
+        -----------
+            
+            x -- list or NumPy array
+                x includes the x-values for the data set to
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
+                    
+            y -- list or NumPy array
+                y includes the y-values for the data set  to
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
+                
+            z -- list or 1D NumPy array
+                z includes the z-values for the data set to
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
+                
+            newx -- list or NumPy array
+                newx includes the x-values for the points at
+                which to perform interpolation.  If it is in array, 
+                it must be 1 or 2-dimensional
+            
+            newy  -- list or NumPy array
+                newy includes the y-values for the points at
+                which to perform interpolation.  If it is in array, 
+                it must be 1 or 2-dimensional
+                
+        Optional Arguments
+        -------------------
+        
+            kind -- Usually a string.  But can be any type.
+                Specifies the type of interpolation to use for values within
+                the range of x.
+                
+                By default, linear interpolation is used.
+                
+                See below for details on other options.
+                
+            out  -- same as for kind
+                How to extrapolate values for outside the rectangle defined by
+                    min(x) <= newx[i] <= max(x)  ,  min(y) <= newy[i] <= max(y)
+                Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
+                a number') for all values below the region.
+                
+            bad_data -- list of numbers
+                List of numerical values (in x, y or z) which indicate unacceptable data. 
+                
+                If bad_data is not None (its default), all points whose x, y or z coordinate is in
+                bad_data, OR ones of whose coordinates is NaN, will be removed.  Note that
+                bad_data != None means NaNs will be removed even if they are not in
+                bad_data.
+                
+        Some Acceptable Input Strings
+        ------------------------
+        
+            "linear" -- linear interpolation : default
+            "spline" -- spline interpolation of default order
+            "quad", "quadratic" -- spline interpolation order 2
+            "cubic" -- spline interpolation order 3
+            "quartic" -- spline interpolation order 4
+            "quintic" -- spline interpolation order 5
+            "526" -- TOMS algorithm 526
+            
+        Other options for kind and out
+        ---------------------------------------------------
+        
+            If you choose to use a non-string argument, you must
+            be careful to use correct formatting.
+            
+            If a function is passed, it will be called when interpolating.
+            It is assumed to have the form 
+                newz = interp(x, y, z, newx, newy), 
+            where x, y, newx, and newy are all numpy arrays.
+            
+            If a callable class is passed, it is assumed to have format
+                instance = Class(x, y, z).
+            which can then be called by
+                newz = instance(newx, newy)
+            
+            If a callable object with method "init_xyz" or "set_xyz" is
+            passed, that method will be used to set x and y as follows
+                instance.set_xy(x, y)
+            and the object will be called during interpolation.
+                newz = instance(newx, newy)
+            If the "init_xyz" and "set_xyz" are not present, it will be called as
+                newz = argument(x, y, z, newx, newy)
+                
+            A primitive type which is not a string signifies a function
+            which is identically that value (e.g. val and 
+            lambda x, y, newx : val are equivalent).
+            
+        Example
+        ---------
+        
+            >>> import numpy
+            >>> from interpolate import Interpolate2d
+            >>> x = range(5)        # note list is permitted
+            >>> y = numpy.arange(5.)
+            >>> z = x+y
+            >>> newx = [.2, 2.3, 2.6, 7.0]
+            >>> newy = [1, 1, 1, 1]
+            >>> interp_func = Interpolate2d(x, y, z)
+            >>> interp_func(newx, newy)
+            array([1.2, 3.3, 3.6, NaN])
+    """
     return Interpolate2d(x, y, z, kind=kind, out=out, bad_data=bad_data)(newx, newy)
 
 # objective interface
 class Interpolate2d:
-    """ A callable class for interpolation of 1D, real-valued data.
+    """ A callable class for interpolation of 2D, real-valued data.
         
         Parameters
         -----------
             
-            x -- list or 1D NumPy array
+            x -- list or NumPy array
                 x includes the x-values for the data set to
-                interpolate from.
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
                     
-            y -- list or 1D NumPy array
+            y -- list or NumPy array
                 y includes the y-values for the data set  to
-                interpolate from.
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
                 
             z -- list or 1D NumPy array
                 z includes the z-values for the data set to
-                interpolate from.
+                interpolate from.  If it is in array, it must be
+                1 or 2-dimensional
                 
         Optional Arguments
         -------------------
@@ -91,6 +202,7 @@
             "cubic" -- spline interpolation order 3
             "quartic" -- spline interpolation order 4
             "quintic" -- spline interpolation order 5
+            "526" -- TOMS algorithm 526
             
         Other options for kind and out
         ---------------------------------------------------
@@ -131,9 +243,8 @@
             >>> newx = [.2, 2.3, 2.6, 7.0]
             >>> newy = [1, 1, 1, 1]
             >>> interp_func = Interpolate2d(x, y, z)
-            >>> interp_fuc(newx, newy)
+            >>> interp_func(newx, newy)
             array([1.2, 3.3, 3.6, NaN])
-            
     """
     def __init__(self, x, y, z, kind='linear', out=NaN, bad_data=None):
         
@@ -148,15 +259,15 @@
         # FIXME : perhaps allow 2D input if it is inthe form of meshgrid
          
         # check acceptable sizes and dimensions
+        # and make data 1D
         x = np.atleast_1d(x)
         y = np.atleast_1d(y)
         z = np.atleast_1d(z)
-        assert len(x) > 0 and len(y) > 0 and len(z)>0, "Arrays cannot be of zero length"
-        assert x.ndim == 1 , "x must be one-dimensional"
-        assert y.ndim == 1 , "y must be one-dimensional"
-        assert z.ndim == 1 , "z must be one-dimensional" 
-        assert len(x) == len(y) , "x and y must be of the same length"
-        assert len(x) == len(z) , "x and z must be of the same length"
+        assert x.shape == y.shape, "x and y must be the same shape"
+        assert x.shape == z.shape, "z must be the same shape as x and y"
+        assert x.ndim in [1, 2], "x and y must be at most 2-dimensional"
+        assert np.size(x) >= 0, "x, y and z cannot be of size zero"
+        x, y, z = map(np.ravel, [x, y, z])
         
         # remove bad data if applicable
         if bad_data is not None:
@@ -246,18 +357,22 @@
                                 isinstance(  newx , np.ndarray  ) and np.shape(newx) == () or \
                                 isinstance(  newy , np.ndarray  ) and np.shape(newy) == ()
         
-        # make input into a nice 1d, contiguous array
-        newx = atleast_1d_and_contiguous(newx, dtype=self._xdtype)
-        newy = atleast_1d_and_contiguous(newy, dtype=self._ydtype)
-        assert newx.ndim == 1, "newx can be at most 1-dimensional"
-        assert newy.ndim == 1, "newy can be at most 1-dimensional"
-        assert len(newx) == len(newy), "newx and newy must be the same length"
+        # check acceptable sizes and dimensions, record the shape of the input,
+        # and make data 1D
+        newx = np.atleast_1d(newx)
+        newy = np.atleast_1d(newy)
+        assert newx.shape == newy.shape, "newx and newy must be the same shape"
+        assert newx.ndim in [1, 2], "newx and newy must be at most 2-dimensional"
+        assert np.size(newx) >= 0, "newx and newy must be of size greater than zero"
+        shape_of_newx = newx.shape
+        newx, newy= map(np.ravel, [newx, newy])
         
+        # finding which points are in the valid interpolation range
         in_range_mask = (min(self._x) <= newx)  & (newx <= max(self._x)) & \
                                 (min(self._y) <= newy) & (newy <= max(self._y))        
         
-        # filling array of interpolated z-values
-        result = np.zeros(np.shape(newx), dtype = self._zdtype)
+        # making array of interpolated z-values
+        result = np.empty(np.shape(newx), dtype = self._zdtype)
         if sum(in_range_mask) > 0:  # if there are in-range values.  hack to deal
                                                # with behavior of np.vectorize on arrays of length 0
             result[in_range_mask] = self.kind(newx[in_range_mask], newy[in_range_mask])        
@@ -267,5 +382,7 @@
         # revert to scalar if applicable
         if input_is_scalar:
             result = result[0]
+        else:
+            result = np.reshape(result, shape_of_newx)
         
         return result
\ No newline at end of file

Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/setup.py	2008-08-19 19:02:03 UTC (rev 4652)
@@ -46,10 +46,12 @@
                                         'extensions/interp_526.f']
                         )
                         
+    # include tutorial on the module
     config.add_data_dir('docs')
 
     return config
 
 if __name__ == '__main__':
+    #from setuptools import setup
     from numpy.distutils.core import setup
     setup(**configuration().todict())
\ No newline at end of file

Modified: branches/Interpolate1D/tests/test_interpolate2d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-19 19:02:03 UTC (rev 4652)
@@ -17,7 +17,72 @@
     
     def assertAllclose(self, x, y):
         self.assert_(np.allclose(atleast_1d_and_contiguous(x), atleast_1d_and_contiguous(y)))
+     
+    def test_string_linear(self):
+        """ make sure : string 'linear' works
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
         
+        newx = np.arange(1, N)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind='linear', out='linear')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_string_quadratic(self):
+        """ make sure : string 'quadratic' works
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = np.arange(1, N)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind='quadratic', out='quad')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_string_cubic(self):
+        """make sure : string "cubic" works
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = np.arange(1, N)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind='cubic', out='cubic')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_string_526(self):
+        """ make sure : keyword '526' works
+            ie that TOMS algorithm 526 works
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = np.arange(1, N)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind='526', out='526')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+
     def test_callable_class_interpolation(self):
         """ make sure : instance of callable class with xyz not initiated works
         """
@@ -126,70 +191,21 @@
         self.assert_( np.isnan(newz[0]) )
         self.assert_( np.isnan(newz[-1]) )
         
-    def test_string_linear(self):
-        """ make sure : string 'linear' works
+    
+    def test_2D_input(self):
+        """ make sure : newx and newy can be 2D, and the result is also 2D
         """
         N = 7
         X, Y = meshgrid(arange(N), arange(N))
         Z = X + Y
         x, y, z = map(ravel, [X, Y, Z] )
         
-        newx = np.arange(1, N)-0.5
-        newy = newx
+        newx, newy = meshgrid(arange(1, N)-.2, arange(1,N)-.2)
         
         interp_func = Interpolate2d(x, y, z, kind='linear', out='linear')
         newz = interp_func(newx, newy)
         
         self.assertAllclose(newz, newx+newy)
         
-    def test_string_quadratic(self):
-        """ make sure : string 'quadratic' works
-        """
-        N = 7
-        X, Y = meshgrid(arange(N), arange(N))
-        Z = X + Y
-        x, y, z = map(ravel, [X, Y, Z] )
-        
-        newx = np.arange(1, N)-0.5
-        newy = newx
-        
-        interp_func = Interpolate2d(x, y, z, kind='quadratic', out='quad')
-        newz = interp_func(newx, newy)
-        
-        self.assertAllclose(newz, newx+newy)
-        
-    def test_string_cubic(self):
-        """make sure : string "cubic" works
-        """
-        N = 7
-        X, Y = meshgrid(arange(N), arange(N))
-        Z = X + Y
-        x, y, z = map(ravel, [X, Y, Z] )
-        
-        newx = np.arange(1, N)-0.5
-        newy = newx
-        
-        interp_func = Interpolate2d(x, y, z, kind='cubic', out='cubic')
-        newz = interp_func(newx, newy)
-        
-        self.assertAllclose(newz, newx+newy)
-        
-    def test_string_526(self):
-        """ make sure : keyword '526' works
-            ie that TOMS algorithm 526 works
-        """
-        N = 7
-        X, Y = meshgrid(arange(N), arange(N))
-        Z = X + Y
-        x, y, z = map(ravel, [X, Y, Z] )
-        
-        newx = np.arange(1, N)-0.5
-        newy = newx
-        
-        interp_func = Interpolate2d(x, y, z, kind='526', out='526')
-        newz = interp_func(newx, newy)
-        
-        self.assertAllclose(newz, newx+newy)
-        
 if __name__ == '__main__':
     unittest.main()                 
\ No newline at end of file



More information about the Scipy-svn mailing list