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

scipy-svn@scip... scipy-svn@scip...
Fri Aug 1 16:47:40 CDT 2008


Author: fcady
Date: 2008-08-01 16:47:36 -0500 (Fri, 01 Aug 2008)
New Revision: 4593

Added:
   branches/Interpolate1D/interpolate2d.py
   branches/Interpolate1D/tests/
   branches/Interpolate1D/tests/regression_test.py
   branches/Interpolate1D/tests/test_fitpack_wrapper.py
   branches/Interpolate1D/tests/test_interpolate1d.py
   branches/Interpolate1D/tests/test_interpolate2d.py
   branches/Interpolate1D/tests/test_interpolate_wrapper.py
Modified:
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/fitpack_wrapper.py
   branches/Interpolate1D/interpolate1d.py
Log:
interpolate2d.py added.  Maybe not the perfect api, but it is functional.  Also added unit tests

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-01 21:47:36 UTC (rev 4593)
@@ -7,7 +7,7 @@
 the most basic mathematical tools available to a researcher.
 
 The interpolate package provides tools for interpolating and extrapolating new data points from a known set of data points.  
-Interpolate provides both a functional interface that is flexible and easy to use as well as an object oriented interface that 
+It provides a functional interface that is flexible and easy to use as well as an object oriented interface that 
 can be more efficient and flexible for some cases.  It is able to interpolate and extrapolate in 1D, 2D, and even N 
 dimensions. *[FIXME : 1D only right now]*
 
@@ -472,24 +472,24 @@
 At instantiation:
 
 #) bbox
-    This is a 2-element list specifying the endpoints of the approximation interval.
-    It default to [x[0],x[-1]]
+This is a 2-element list specifying the endpoints of the approximation interval.
+It default to [x[0],x[-1]]
     
 #) w
-    a 1D sequence of weights which defaults to all ones.
+a 1D sequence of weights which defaults to all ones.
 
 #) 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
+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
     
-    BEWARE : in the current implementation of the code, if s is small but not zero,
-            instantiating Spline can become painfully slow.
+BEWARE : in the current implementation of the code, if s is small but not zero,
+    instantiating Spline can become painfully slow.
 
 At calling:
 
 #) 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.
+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
@@ -497,16 +497,73 @@
 
 #) set_smoothing_factor(s)
 #) get_knots
-    returns the positions of the knots of the spline
+returns the positions of the knots of the spline
 #) get_coeffs
-    returns the coefficients of the 
+returns the coefficients of the 
 #) 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)
+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)
-    returns the integral from a to b
+returns the integral from a to b
 #) derivatives(x)
-    returns all the derivatives of the spline at point x
+returns all the derivatives of the spline at point x
 #) roots
-    This only works for cubic splines.  But it returns the places where the spline
-    is identically zero.
+This only works for cubic splines.  But it returns the places where the spline
+is identically zero.
+
+
+================================================
+2D Interpolation
+================================================
+
+*[This is being written preemptively]*
+
+In 2D interpolation, known data are of the form (x, y, z), and we interpolate
+z_new from (x_new, y_new).  
+
+As in the case of 1D interpolation, there is a convenient functional interface
+for 2D interpolation as well as a callable object which can be more efficient.
+In analogy to 1D interpolation, the function is interp2d and the class is Interpolate2d.
+
+------------------------------------------
+The Functional Interface
+------------------------------------------
+
+The functional interface is virtually identical to that for interp1d: ::
+
+    new_z = interp2d(x, y, z, new_x, new_y)
+
+The range of interpolation is the rectangle given by the largest and
+smallest values in x and y.
+As in the case of 1D, string arguments can be passed
+as keywords to specify particular types of interpolation.  The keywords
+in this case are kind (for in-range interpolation) and out (for out-of-bounds).
+
+By default, out of bounds returns NaN, and in-bounds returns a linear
+interpolation.
+
+new_x and new_y may be either arrays, lists or scalars.  If they are
+scalars or zero-dimensional arrays, new_z will be a scalar as well.  Otherwise
+a vector is returned.
+
+------------------------------------------
+The Objective Interface
+------------------------------------------
+
+The objective interface for 2D is slightly more complicated.  
+
+================================================
+ND Interpolation
+================================================
+
+
+
+================================================
+ND Scattered Interpolation
+================================================
+ 
+ Still in development.
+ 
+ Ideally the range of interpolation would be the convex hull of the known
+ data points, and a Delaunay triangulation would be determined and stored
+ at instantiation.  Then again, that would be VERY expensive.
\ No newline at end of file

Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/fitpack_wrapper.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -1,26 +1,20 @@
 """
-This module is used for spline interpolation, and functions
-as a wrapper around the FITPACK Fortran interpolation
-package.  Its functionality is contained in the Spline
-class.
+    This module is used for spline interpolation, and functions
+    as a wrapper around the FITPACK Fortran interpolation
+    package.  Its functionality is contained in the Spline
+    class.
 
-Spline is primarily meant to be called by Interpolate1d
-or inter1d, but is a stand-alone class in its own right.
+    Spline is primarily meant to be called by Interpolate1d
+    or interp1d, but is a stand-alone class in its own right
+    that is not available through these interfaces.
 
+    The code has been modified from an older version of
+    scipy.interpolate, where it was directly called by the
+    user.  As such, it includes functionality not available through
+    Interpolate1d.  For this reason, users may wish to get
+    under the hood.
 
-The code has been modified from an older version of
-scipy.interpolate, where it was directly called by the
-user.  As such, it includes functionality not available through
-Interpolate1d.  For this reason, users may wish to get
-under the hood.
-
 """
-# FIXME : CLEAN UP THIS FILE!  scipy.interpolate contained a lot of
-#       nice functionality that is only partially in this file.
-#       The question is whether to copy over the full functionality
-#       to the point where we may as well include fitting.py from
-#       scipy.interpolate, or whether we should strip this down some.
-#       Until that's decided, cleaning is premature.
 
 import numpy as np
 
@@ -29,26 +23,26 @@
 
 class Spline(object):
     """ Univariate spline s(x) of degree k on the interval
-    [xb,xe] calculated from a given set of data points
-    (x,y).
+        [xb,xe] calculated from a given set of data points
+        (x,y).
 
-    Can include least-squares fitting.
+        Can include least-squares fitting.
 
     """
 
     def __init__(self, x=None, y=None, w=None, bbox = [None]*2, k=3, s=0.0):
         """
         Input:
-          x,y   - 1-d sequences of data points (x must be
+            x,y   - 1-d sequences of data points (x must be
                   in strictly ascending order)
 
         Optional input:
-          w          - positive 1-d sequence of weights
-          bbox     - 2-sequence specifying the boundary of
+            k=3     - degree of the univariate spline.
+            w       - positive 1-d sequence of weights
+            bbox   - 2-sequence specifying the boundary of
                        the approximation interval.
                        By default, bbox=[x[0],x[-1]]
-          k=3     - degree of the univariate spline.
-          s          - positive smoothing factor defined for
+            s        - positive smoothing factor defined for
                        estimation condition:
                          sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
                         Default s=0

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/interpolate1d.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -10,18 +10,19 @@
 # dictionary of interpolation functions/classes/objects
 method_register = \
                 { # functions
-                    'linear' : linear, 
-                    'logarithmic' : logarithmic, 
-                    'block' : block, 
+                    'linear' : linear,  'Linear' : linear, 
+                    'logarithmic' : logarithmic, 'Logarithmic' : logarithmic, 
+                    'block' : block, 'Block' : block, 
                     'block_average_above' : block_average_above, 
-                    'nearest' : nearest,
+                    'Block_average_above' : block_average_above, 
+                    'nearest' : nearest, 'Nearest' : nearest,
                     
                     # Splines
                     'Spline' : Spline, 'spline' : Spline,
                     'Quadratic' : Spline(k=2), 'quadratic' : Spline(k=2),
                     'Quad' : Spline(k=2), 'quad' : Spline(k=2),
                     'Cubic' : Spline(k=3), 'cubic' : Spline(k=3),
-                    'Quartic' : Spline(k=3), 'quartic' : Spline(k=3),
+                    'Quartic' : Spline(k=4), 'quartic' : Spline(k=4),
                     'Quar' : Spline(k=4), 'quar' : Spline(k=4),
                     'Quintic' : Spline(k=5), 'quintic' : Spline(k=5),
                     'Quin' : Spline(k=5), 'quin' : Spline(k=5)
@@ -76,20 +77,14 @@
                 Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
                 a number') for all values above the range of x.
                 
-            bad_data -- list
+            bad_data -- list of numbers
                 List of numerical values (in x or y) which indicate unacceptable data. 
                 
                 If bad_data is not None (its default), all points whose x or y coordinate is in
-                bad_data, OR ones of whose coordinates is NaN, will be removed.
+                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.
                 
-            interpkw -- dictionary
-                If interp is set to a function, class or callable object, this contains
-                additional keywords.
-                
-            lowkw (highkw) -- dictionary
-                like interpkw, but for extrap_low and extrap_high
-                
-            
         Some Acceptable Input Strings
         ------------------------
         
@@ -97,14 +92,13 @@
             "logarithmic" -- logarithmic interpolation : linear in log space?
             "block" --
             "block_average_above' -- block average above
-            "Spline" -- spline interpolation.  keyword k (defaults to 3) 
-                indicates order of spline
+            "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
             
-        Other options for interp, extrap_low, and extrap_high
+        Other options for kind, low, and high
         ---------------------------------------------------
         
             If you choose to use a non-string argument, you must
@@ -112,21 +106,21 @@
             
             If a function is passed, it will be called when interpolating.
             It is assumed to have the form 
-                newy = interp(x, y, newx, **kw), 
+                newy = interp(x, y, newx), 
             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, **kw).
+                instance = Class(x, y).
             which can then be called by
                 new_y = instance(new_x)
             
             If a callable object with method "init_xy" or "set_xy" is
             passed, that method will be used to set x and y as follows
-                instance.set_xy(x, y, **kw)
+                instance.set_xy(x, y)
             and the object will be called during interpolation.
                 new_y = instance(new_x)
             If the "init_xy" and "set_xy" are not present, it will be called as
-                new_y = argument(new_x)
+                new_y = argument(x, y, new_x)
                 
             A primitive type which is not a string signifies a function
             which is identically that value (e.g. val and 
@@ -187,11 +181,13 @@
                 Same options as for 'kind'.  Defaults to returning numpy.NaN ('not 
                 a number') for all values above the range of x.
                 
-            bad_data -- list
+            bad_data -- list of numbers
                 List of numerical values (in x or y) which indicate unacceptable data. 
                 
                 If bad_data is not None (its default), all points whose x or y coordinate is in
-                bad_data, OR ones of whose coordinates is NaN, will be removed.
+                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
         ------------------------
@@ -200,8 +196,7 @@
             "logarithmic" -- logarithmic interpolation : linear in log space?
             "block" --
             "block_average_above' -- block average above
-            "Spline" -- spline interpolation.  keyword k (defaults to 3) 
-                indicates order of spline
+            "Spline" -- spline interpolation of default order
             "quad", "quadratic" -- spline interpolation order 2
             "cubic" -- spline interpolation order 3
             "quartic" -- spline interpolation order 4
@@ -229,7 +224,7 @@
             and the object will be called during interpolation.
                 new_y = instance(new_x)
             If the "init_xy" and "set_xy" are not present, it will be called as
-                new_y = argument(new_x)
+                new_y = argument(x, y, new_x)
                 
             A primitive type which is not a string signifies a function
             which is identically that value (e.g. val and 
@@ -239,7 +234,7 @@
         ---------
         
             >>> import numpy
-            >>> from interpolate1d import Interpolate1d
+            >>> from interpolate import Interpolate1d
             >>> x = range(5)        # note list is permitted
             >>> y = numpy.arange(5.)
             >>> new_x = [.2, 2.3, 5.6, 7.0]

Added: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/interpolate2d.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,261 @@
+
+from numpy import NaN, array
+import numpy as np
+from fitpack_wrapper2d import Spline2d
+
+def atleast_1d_and_contiguous(ary, dtype = np.float64):
+    # FIXME : don't have in 2 places
+    return np.atleast_1d( np.ascontiguousarray(ary, dtype) )
+
+# dictionary of interpolation functions/classes/objects
+method_register = \
+                { 
+                    'linear' : Spline2d(kx=1, ky=1),
+                    'Linear' :  Spline2d(kx=1, ky=1),
+                    'Spline' : Spline2d(), 'spline' : Spline2d(),
+                    'Quadratic' : Spline2d(kx=2, ky=2), 'quadratic' : Spline2d(kx=2, ky=2),
+                    'Quad' : Spline2d(kx=2, ky=2), 'quad' : Spline2d(kx=2, ky=2),
+                    'Cubic' : Spline2d(kx=3, ky=3), 'cubic' : Spline2d(kx=3, ky=3),
+                    'Quartic' : Spline2d(kx=4, ky=4), 'quartic' : Spline2d(kx=4, ky=4),
+                    'Quar' : Spline2d(kx=4, ky=4), 'quar' : Spline2d(kx=4, ky=4),
+                    'Quintic' : Spline2d(kx=5, ky=5), 'quintic' : Spline2d(kx=5, ky=5),
+                    'Quin' : Spline2d(kx=5, ky=5), 'quin' : Spline2d(kx=5, ky=5)
+                }
+                
+# dictionary of types for casting.  key = possible datatype, value = datatype it is cast to
+# BEWARE : if you cast things to integers, you will lose interpolation ability
+dtype_register = {np.float32 : np.float32, 
+                            np.float64 : np.float64
+                            }
+dtype_default = np.float64
+
+def interp2d(x, y, z, newx, newy, kind='linear', out=NaN, bad_data=None):
+    return Interpolate2d(x, y, z, kind=kind, out=out, bad_data=bad_data)(newx, newy)
+    
+class Interpolate2d:
+    """ A callable class for interpolation of 1D, real-valued data.
+        
+        Parameters
+        -----------
+            
+            x -- list or 1D NumPy array
+                x includes the x-values for the data set to
+                interpolate from.
+                    
+            y -- list or 1D NumPy array
+                y includes the y-values for the data set  to
+                interpolate from.
+                
+            z -- list or 1D NumPy array
+                z includes the z-values for the data set to
+                interpolate from.
+                
+        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
+            
+        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_fuc(newx, newy)
+            array([1.2, 3.3, 3.6, NaN])
+            
+    """
+    def __init__(self, x, y, z, kind='linear', out=NaN, bad_data=None):
+        
+        self._init_xyz(x, y, z, bad_data)
+        
+        self.kind = self._init_interp_method(kind)
+        self.out = self._init_interp_method(out)
+        
+    def _init_xyz(self, x, y, z, bad_data):
+        # all must be vectors.  Should allow meshgrid, though?  Perhaps.
+        # Mention that in documentation.
+        
+        if bad_data is not None:
+            try: # check that bad_data contains only numerical values
+                sum_of_bad_data = sum(bad_data)
+            except:
+                raise TypeError, "bad_data must be either None \
+                        or a list of numbers"  
+            x, y, z = self._remove_bad_data(x, y, z, bad_data)
+         
+        # check acceptable sizes and dimensions
+        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"
+            
+        # select proper dataypes and make arrays
+        self._xdtype = dtype_register.setdefault(type(x[0]), dtype_default)
+        self._ydtype = dtype_register.setdefault(type(y[0]), dtype_default)
+        self._zdtype = dtype_register.setdefault(type(z[0]), dtype_default)
+        self._x = atleast_1d_and_contiguous(x, self._xdtype).copy()
+        self._y = atleast_1d_and_contiguous(y, self._ydtype).copy()
+        self._z = atleast_1d_and_contiguous(z, self._zdtype).copy()
+                
+    def _init_interp_method(self, method):
+        """ returns the interpolating function specified by interp_arg.
+        """        
+        from inspect import isclass, isfunction
+        
+        # primary usage : user passes a string indicating a known function
+        # pick interpolator accordingly
+        if isinstance(method, basestring):
+            interpolator = method_register.setdefault(method, None )
+            if interpolator is None: 
+                raise TypeError, "input string %s not valid" % method
+        else:
+            interpolator = method
+        
+        # interpolator is a callable : function, class, or instance of class
+        if hasattr(interpolator, '__call__'):
+            # function
+            if isfunction(interpolator):
+                result = lambda newx, newy : interpolator(self._x, self._y, self._z, newx, newy)
+                
+            # callable class 
+            elif isclass(interpolator):
+                if hasattr(interpolator, 'set_xyz'):
+                    result = interpolator()
+                    result.set_xyz(self._x, self._y, self._z)
+                if hasattr(interpolator, 'init_xyz'):
+                    result = interpolator()
+                    result.init_xyz(self._x, self._y, self._z)
+                else:
+                    result = interpolator(self._x, self._y, self._z)
+                
+            # instance of callable class
+            else:
+                if hasattr(interpolator, 'init_xyz'):
+                    result = interpolator
+                    result.init_xyz(self._x, self._y, self._z)
+                elif hasattr(interpolator, 'set_xyz'):
+                    result = interpolator
+                    result.set_xyz(self._x, self._y, self._z)
+                else:
+                    result = lambda newx, newy : interpolator(self._x, self._y, self._z, newx, newy)
+            
+        # non-callable : user has passed a default value to always be returned
+        else:
+            result = np.vectorize(lambda newx, newy : interpolator)
+        
+        return result
+        
+    def _remove_bad_data(self, x, y, z, bad_data):
+        """ removes data points whose x or y coordinate is
+            either in bad_data or is a NaN.
+        """
+
+        bad_data_mask = np.isnan(x) | np.isnan(y) | np.isnan(z)
+        for bad_num in bad_data:
+            bad_data_mask =  \
+                    bad_data_mask | (x==bad_num) | (y==bad_num) | (z==bad_num)
+        
+        x = x[~bad_data_mask]
+        y = y[~bad_data_mask]
+        z = z[~bad_data_mask]
+        
+        return x, y, z
+        
+    def __call__(self, newx, newy):
+        
+        # record if input is scalar or 0-dimemsional array, in which case output will be scalar
+        input_is_scalar = np.isscalar(newx) or  np.isscalar(newy) or \
+                                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)
+        assert newx.ndim == 1, "newx can be at most 1-dimensional"
+        newy = atleast_1d_and_contiguous(newy, dtype=self._ydtype)
+        assert newy.ndim == 1, "newy can be at most 1-dimensional"
+        assert len(newx) == len(newy), "newx and newy must be the same length"
+        
+        in_range_mask = (min(self._x) <= newx)  & (newx <= max(self._x)) & \
+                                (min(self._y) <= newy) & (newy <= max(self._y))        
+        
+        result = np.zeros(np.shape(newx))
+        if sum(in_range_mask) > 0:
+            # hack to deal with behavior of vectorize on arrays of length 0
+            result[in_range_mask] = self.kind(newx[in_range_mask], newy[in_range_mask])        
+        if sum(~in_range_mask) > 0:
+            # hack to deal with behavior of vectorize on arrays of length 0
+            result[~in_range_mask] = self.out(newx[~in_range_mask], newy[~in_range_mask])
+        
+        if input_is_scalar:
+            result = result[0]
+        
+        return result
\ No newline at end of file

Added: branches/Interpolate1D/tests/regression_test.py
===================================================================
--- branches/Interpolate1D/tests/regression_test.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/tests/regression_test.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,37 @@
+""" 
+    regression test:
+
+    This script runs a simple regression test on the functionality of
+    the interpolation module.  Currently, when run, it times each
+    unit test in interpolate1d.py and stores those times in a dict
+    of dicts; outer keys are time test was performed, and inner
+    keys are names of tests run.
+
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+import shelve, time
+from test_interpolate1d import Test
+
+# name of log file to which all data is stored.
+filename = 'regression_test.dbm'
+
+log_total = shelve.open(filename)
+current_time = str(time.localtime()[0:5]) # specified up to the minute
+
+# run all tests in interpolate1d's test class
+test_list = [name for name in dir(Test) if name.find('test_') == 0]
+log_now = {}
+
+# record time taken for each test
+for test_name in test_list:
+    t1 = time.clock()
+    eval('Test.%s' % test_name)
+    t2 = time.clock()
+    log_now[test_name] = t2-t1
+
+log_total[current_time] = log_now
+log_total.close()

Added: branches/Interpolate1D/tests/test_fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,76 @@
+""" module to test fitpack_wrapper.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
+import numpy as np
+from fitpack_wrapper import Spline
+
+class Test(unittest.TestCase):
+    
+    def assertAllclose(self, x, y):
+        self.assert_(np.allclose(x, y))
+        
+    def test_linearInterp(self):
+        """ make sure : linear interpolation (spline with order = 1, s = 0)works
+        """
+        N = 3000.
+        x = np.arange(N)
+        y = np.arange(N)
+        #T1 = time.clock()
+        interp_func = Spline(x, y, k=1)
+        #T2 = time.clock()
+        #print "time to create order 1 spline interpolation function with N = %i:" % N, T2 - T1
+        new_x = np.arange(N)+0.5
+        #t1 = time.clock()
+        new_y = interp_func(new_x)
+        #t2 = time.clock()
+        #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):
+        """ make sure : quadratic interpolation (spline with order = 2, s = 0)works
+        """
+        N = 3000.
+        x = np.arange(N)
+        y = x**2
+        interp_func = Spline(x, y, k=2)
+        #print "time to create order 1 spline interpolation function with N = %i:" % N, T2 - T1
+        new_x = np.arange(N)+0.5
+        #t1 = time.clock()
+        new_y = interp_func(x)
+        #t2 = time.clock()
+        #print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
+        self.assertAllclose(new_y, y)
+        
+        
+    def test_inputFormat(self):
+        """ make sure : it's possible to instantiate Spline without x and y
+        """
+        #print "testing input format"
+        N = 3000.
+        x = np.arange(N)
+        y = np.arange(N)
+        interp_func = Spline(k=1)
+        interp_func.init_xy(x, y)
+        new_x = np.arange(N)+0.5
+        new_y = interp_func(new_x)
+        self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
+    
+    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

Added: branches/Interpolate1D/tests/test_interpolate1d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate1d.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/tests/test_interpolate1d.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,170 @@
+""" Test module for interpolate1d.py
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+import numpy as np
+from numpy import arange
+
+from interpolate1d import interp1d, Interpolate1d, atleast_1d_and_contiguous
+from fitpack_wrapper import Spline
+
+# unit testing
+import unittest, time
+class Test(unittest.TestCase):
+    
+    def assertAllclose(self, x, y):
+        self.assert_(np.allclose(atleast_1d_and_contiguous(x), atleast_1d_and_contiguous(y)))
+        
+    def test_interpolate_wrapper(self):
+        """ run unit test contained in interpolate_wrapper.py
+        """
+        #print "\n\nTESTING _interpolate_wrapper MODULE"
+        from test_interpolate_wrapper import Test
+        T = Test()
+        T.runTest()
+        
+    def test_fitpack_wrapper(self):
+        """ run unit test contained in fitpack_wrapper.py
+        """
+        #print "\n\nTESTING _fitpack_wrapper MODULE"
+        from test_fitpack_wrapper import Test
+        T = Test()
+        T.runTest()
+        
+    def test_callable_class_interpolation(self):
+        """ make sure : an instance of a callable class in which
+            x and y haven't been initiated works
+        """
+        N = 7 #must be > 5
+        x = np.arange(N)
+        y = np.arange(N)
+        interp_func = Interpolate1d(x, y, kind=Spline(k=2), low=Spline, high=Spline)
+        new_x = np.arange(N+1)-0.5
+        new_y = interp_func(new_x)
+        self.assertAllclose(new_x, new_y)
+        
+    def test_no_out_of_range_args(self):
+        """ make sure : having no out-of-range elements in new_x is fine
+        """
+        # There was a bug with this earlier.        
+        N = 5
+        x = arange(N)
+        y = arange(N)
+        new_x = arange(1,N-1)+.2
+        interp_func = Interpolate1d(x, y, kind='linear', low='linear', high=np.NaN)
+        new_y = interp_func(new_x)
+        self.assertAllclose(new_x, new_y)
+        
+    def test_remove_bad_data(self):
+        """make sure : interp1d works with bad data
+        """
+        N = 7.0 # must be >=5
+        x = arange(N); x[2] = np.NaN
+        y = arange(N); y[4] = 52.3; y[0]=np.NaN
+        new_x = arange(N+1)-0.5
+        new_y = interp1d(x, y, new_x, kind='linear', low='linear', 
+                                high='linear', bad_data = [52.3])
+        self.assertAllclose(new_x, new_y)
+        
+    def test_interp1d(self):
+        """ make sure : interp1d works, at least in the linear case
+        """
+        N = 7
+        x = arange(N)
+        y = arange(N)
+        new_x = arange(N+1)-0.5
+        new_y = interp1d(x, y, new_x, kind='linear', low='linear', high='linear')        
+        self.assertAllclose(new_x, new_y)
+        
+        
+    def test_spline1_default_extrapolation(self):
+        """ make sure : spline order 1 (linear) interpolation works correctly
+            make sure : default extrapolation works
+        """
+        #print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
+        N = 7 # must be > 5
+        x = np.arange(N)
+        y = np.arange(N)
+        interp_func = Interpolate1d(x, y, kind=Spline(k=1), low=None, high=599.73)
+        new_x = np.arange(N+1)-0.5
+        new_y = interp_func(new_x)
+        
+        self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
+        self.assert_(new_y[0] == None)
+        self.assert_(new_y[-1] == 599.73)
+        
+    def test_spline2(self):
+        """ make sure : order-2 splines work on linear data
+            make sure : order-2 splines work on non-linear data
+            make sure : 'cubic' and 'quad' as arguments yield
+                                the desired spline
+        """
+        #print "\n\nTESTING 2nd ORDER SPLINE"
+        N = 7 #must be > 5
+        x = np.arange(N)
+        y = np.arange(N)
+        T1 = time.clock()
+        interp_func = Interpolate1d(x, y, kind=Spline(k=2), low='spline', high='spline')
+        T2 = time.clock()
+        #print "time to create 2nd order spline interp function with N = %i: " % N, T2 - T1
+        new_x = np.arange(N+1)-0.5
+        t1 = time.clock()
+        new_y = interp_func(new_x)
+        t2 = time.clock()
+        #print "time to evaluate 2nd order spline interp function with N = %i: " % N, t2 - t1
+        self.assertAllclose(new_x, new_y)
+        
+        # make sure for non-linear data
+        N = 7
+        x = np.arange(N)
+        y = x**2
+        interp_func = Interpolate1d(x, y, kind=Spline(k=2), low='quad', high='cubic')
+        new_x = np.arange(N+1)-0.5
+        new_y = interp_func(new_x)
+        self.assertAllclose(new_x**2, new_y)
+        
+        
+    def test_linear(self):
+        """ make sure : linear interpolation works 
+            make sure : linear extrapolation works
+        """
+        #print "\n\nTESTING LINEAR INTERPOLATION"
+        N = 7
+        x = arange(N)
+        y = arange(N)
+        new_x = arange(N+1)-0.5
+        T1 = time.clock()
+        interp_func = Interpolate1d(x, y, kind='linear', low='linear', high='linear')
+        T2 = time.clock()
+        #print "time to create linear interp function with N = %i: " % N, T2 - T1
+        t1 = time.clock()
+        new_y = interp_func(new_x)
+        t2 = time.clock()
+        #print "time to create linear interp function with N = %i: " % N, t2 - t1
+        
+        self.assertAllclose(new_x, new_y)
+        
+    def test_scalar_input(self):
+        """ make sure : newx being a scalar or a 0-degree array
+                            makes the output a scalar
+        """
+        N = 7
+        x = arange(N)
+        y = arange(N)
+        interp_func = Interpolate1d(x, y, kind='linear', low='linear', high='linear')
+        
+        # scalar input
+        newx1 = 0.5
+        newy1 = interp_func(newx1)
+        self.assert_( np.isscalar(newy1) )
+        
+        # zero-degree array
+        newx2 = np.array(0.5)
+        newy2 = interp_func(newx2)
+        self.assert_( np.isscalar(newy2) )
+        
+if __name__ == '__main__':
+    unittest.main()                 
\ No newline at end of file

Added: branches/Interpolate1D/tests/test_interpolate2d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,178 @@
+""" Test module for interpolate1d.py
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+import numpy as np
+from numpy import arange, meshgrid, ravel
+
+from interpolate2d import interp2d, Interpolate2d, atleast_1d_and_contiguous
+from fitpack_wrapper2d import Spline2d
+
+# unit testing
+import unittest, time
+class Test(unittest.TestCase):
+    
+    def assertAllclose(self, x, y):
+        self.assert_(np.allclose(atleast_1d_and_contiguous(x), atleast_1d_and_contiguous(y)))
+        
+    def test_callable_class_interpolation(self):
+        """ make sure : instance of callable class with xyz not initiated works
+        """
+        N = 7
+        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 = Interpolate2d(x, y, z, kind=Spline2d(kx=1, ky=1))
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_no_out_of_range_args(self):
+        """ make sure : having no out-of-range elements in new_x is fine
+        """
+        # There was a bug with this earlier.        
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = arange(1,N-1)+.2
+        newy = arange(1,N-1)+.3
+        
+        interp_func = Interpolate2d(x, y, z, kind='linear', out = 'linear')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx + newy)
+        
+    def test_remove_bad_data(self):
+        """make sure : works with bad data
+        """
+        N = 5.
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        x[2] = 55.0
+        z[4] = np.NaN
+        
+        newx = arange(1, N-1)-0.5
+        newy = newx
+        
+        newz = interp2d(x, y, z, newx, newy, kind='linear', bad_data = [55.0])
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_interp2d(self):
+        """ make sure : interp2d works, at least in the linear case
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = arange(1, N)-0.5
+        newy = newx
+        
+        newz = interp2d(x, y, z, newx, newy, kind='linear', out='linear')    
+        
+        self.assertAllclose(newz, newx+newy)
+        
+    def test_scalar_input(self):
+        """ make sure : scalar input or a 0-degree array works
+        """
+        
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+                
+        interp_func = Interpolate2d(x, y, z, kind='linear', out='linear')
+        
+        # scalar input
+        newx1 = 0.5
+        newy1 = 0.7
+        
+        newz1 = interp_func(newx1, newy1)
+        self.assert_( np.isscalar(newz1) )
+        
+        # zero-degree array
+        newx2 = 0.5
+        newy2 = 0.7
+        
+        newz2 = interp_func(newx2, newy2)
+        self.assert_( np.isscalar(newz2) )
+        
+    def test_extrapolation(self):
+        """ make sure : default extrapolation works
+        """
+        N = 7
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = np.arange(0, N+1)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind=Spline2d(kx=1, ky=1) )
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz[1:5], newx[1:5]+newy[1:5] )
+        self.assert_( np.isnan(newz[0]) )
+        self.assert_( np.isnan(newz[-1]) )
+        
+    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)
+        
+if __name__ == '__main__':
+    unittest.main()                 
\ No newline at end of file

Added: branches/Interpolate1D/tests/test_interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate_wrapper.py	2008-08-01 00:06:19 UTC (rev 4592)
+++ branches/Interpolate1D/tests/test_interpolate_wrapper.py	2008-08-01 21:47:36 UTC (rev 4593)
@@ -0,0 +1,90 @@
+""" module to test interpolate_wrapper.py
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+# Unit Test
+import unittest
+import time
+from numpy import arange, allclose, ones, NaN, isnan
+import numpy as np
+
+# functionality to be tested
+from interpolate_wrapper import atleast_1d_and_contiguous, \
+        linear, logarithmic, block_average_above, block, nearest
+
+class Test(unittest.TestCase):
+    
+    def assertAllclose(self, x, y, rtol=1.0e-5):
+        for i, xi in enumerate(x):
+            self.assert_(allclose(xi, y[i], rtol) or (isnan(xi) and isnan(y[i])))
+            
+    def test_nearest(self):
+        N = 5
+        x = arange(N)
+        y = arange(N)
+        self.assertAllclose(y, nearest(x, y, x+.1))
+        self.assertAllclose(y, nearest(x, y, x-.1))
+        
+    def test_linear(self):
+        N = 3000.
+        x = arange(N)
+        y = arange(N)
+        new_x = arange(N)+0.5
+        t1 = time.clock()
+        new_y = linear(x, y, new_x)
+        t2 = time.clock()
+        #print "time for linear interpolation with N = %i:" % N, t2 - t1
+        
+        self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
+        
+    def test_block_average_above(self):
+        N = 3000.
+        x = arange(N)
+        y = arange(N)
+        
+        new_x = arange(N/2)*2
+        t1 = time.clock()
+        new_y = block_average_above(x, y, new_x)
+        t2 = time.clock()
+        #print "time for block_avg_above interpolation with N = %i:" % N, t2 - t1
+        self.assertAllclose(new_y[:5], [0.0, 0.5, 2.5, 4.5, 6.5])
+
+    def test_linear2(self):
+        N = 3000.
+        x = arange(N)
+        y = ones((100,N)) * arange(N)
+        new_x = arange(N)+0.5
+        t1 = time.clock()
+        new_y = linear(x, y, new_x)
+        t2 = time.clock()
+        #print "time for 2D linear interpolation with N = %i:" % N, t2 - t1
+        self.assertAllclose(new_y[:5,:5],
+                            [[ 0.5, 1.5, 2.5, 3.5, 4.5],
+                             [ 0.5, 1.5, 2.5, 3.5, 4.5],
+                             [ 0.5, 1.5, 2.5, 3.5, 4.5],
+                             [ 0.5, 1.5, 2.5, 3.5, 4.5],
+                             [ 0.5, 1.5, 2.5, 3.5, 4.5]])
+                             
+    def test_logarithmic(self):
+        N = 4000.
+        x = arange(N)
+        y = arange(N)
+        new_x = arange(N)+0.5
+        t1 = time.clock()
+        new_y = logarithmic(x, y, new_x)
+        t2 = time.clock()
+        #print "time for logarithmic interpolation with N = %i:" % N, t2 - t1
+        correct_y = [np.NaN, 1.41421356, 2.44948974, 3.46410162, 4.47213595]
+        self.assertAllclose(new_y[:5], correct_y)
+        
+    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