[Scipy-svn] r4577 - branches/Interpolate1D

scipy-svn@scip... scipy-svn@scip...
Tue Jul 29 12:56:53 CDT 2008


Author: fcady
Date: 2008-07-29 12:56:52 -0500 (Tue, 29 Jul 2008)
New Revision: 4577

Modified:
   branches/Interpolate1D/TODO.txt
   branches/Interpolate1D/__init__.py
   branches/Interpolate1D/erics_notes.txt
   branches/Interpolate1D/fitpack_wrapper.py
   branches/Interpolate1D/info.py
   branches/Interpolate1D/interpolate1d.py
   branches/Interpolate1D/interpolate_wrapper.py
Log:
incorporated many of Eric's suggestions

Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/TODO.txt	2008-07-29 17:56:52 UTC (rev 4577)
@@ -103,7 +103,7 @@
 
 **update for 2D and ND
     This will probably take the form of two additional
-    classes both based on interpolate1d.  Thus it probably
+    classes both modelled after interpolate1d.  Thus it probably
     shouldn't be done until interpolate1d is more settled.
     
     There is an interesting problem here.  Most of the extensions
@@ -114,7 +114,7 @@
     way to do it at first, just to get it working)
     
     We should probably use something other than FITPACK for this.
-    First off it's at most 2D.  But much worse, it doesn't evaluate at
+    But firstly it's at most 2D.  Much worse, it doesn't evaluate at
     a set of points; it evaluates over a grid, which requires inputs
     being in sorted order (both in x and y coordinates).  This makes
     input inconvenient and the code runs a lot slower than ndimage.
@@ -127,6 +127,13 @@
     based on FITPACK (or something else, depending on what
     happens with the smoothing parameter), and InterpolateNd
     which is based on ndimage.
+    
+    Another option is to have two classes: one for uniformly spaced data and
+    and one for scattered data.  Regularly spaced would use NDImage, and
+    scattered would start as an inefficient wrapper around FITPACK.  But longer
+    term the scattered data could be done with Delaunay triangulation,
+    perhaps something that would implicitly calculate the convex hull
+    of the points and interpolate within it.
 
 **high-level road map
     when the module is more established, there should be a page on

Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/__init__.py	2008-07-29 17:56:52 UTC (rev 4577)
@@ -9,4 +9,4 @@
 from fitpack_wrapper import Spline
 
 # wrapping for all supported interpolation types.
-from interpolate1d import Interpolate1d, interp1d
\ No newline at end of file
+from interpolate1d import Interpolate1d, interp1
\ No newline at end of file

Modified: branches/Interpolate1D/erics_notes.txt
===================================================================
--- branches/Interpolate1D/erics_notes.txt	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/erics_notes.txt	2008-07-29 17:56:52 UTC (rev 4577)
@@ -1,9 +1,39 @@
-*. I am glad to see your docstrings.
+*. Function signatures:
+    
+    def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
+                        kindkw={}, lowkw={}, highkw={}, \
+                        remove_bad_data = False, bad_data=[], interp_axis = 0):
 
-*. Move tests into a seprate test directory.
+     o. the trailing \ for new lines are not necessary for functions since line continuation
+        is implicit with the open/close parentheses.
+     o. [minor] I would just use NaN instead of np.NaN.  						
+     o. It is dangerous to initialize containers as default arguments because they are
+        effectively a singleton for the function.  We can discuss this if you don't know
+        what I am talking about.
+     o. kindkw, lowkw, and highkw aren't really necessary I don't think.  They should go.
+	 Xo. Do you need both remove_bad_data and bad_data?  If bad_data is None, then you 
+	    don't remove bad_data...
+	 Xo. I think I would change interp_axis to just be axis.  This is consistent with many
+	    other functions in numpy.
+	 Xo. The choice of whether axis=0 or axis=-1 by default is a reasonable question.
+	    fft defaults to axis=-1.  This is also the faster axis to operate across in the
+	    standard case.  It is, however, the opposite of how some people think about 
+	    things (columns vs. rows).  Talk to Travis O. for his take.  Mine is to use axis=-1.
+	
+	I think all of this might simplify the interface to the following:
+	
+    def interp1d(x, y, new_x, kind='linear', low=NaN, high=NaN, bad_data=None, axis=-1):
 
-*. Follow scipy/FORMAT_GUIDLINES.txt in the main scipy directory.
+*. _remove_bad_data should be left up to the interpolation method if it knows what to do.
+   otherwise, it is handled by this top level class.
+   And, we definitely don't want the list comprehension in the remove_bad_data class.
 
+X*. I am glad to see your docstrings.
+
+X*. Move tests into a seprate test directory.
+
+X*. Follow scipy/FORMAT_GUIDLINES.txt in the main scipy directory.
+
 For example:
 
 	test_callFormat -> test_call_format
@@ -14,14 +44,14 @@
 	describes them "here":http://www.python.org/doc/essays/styleguide.html.  A few
 	reminders follow:
 
-	   o Use 4 spaces for indentation levels.  Do not use tabs as they can result
+	   Xo Use 4 spaces for indentation levels.  Do not use tabs as they can result
 	     in indentation confusion.  Most editors have a feature that will insert 4
 	     spaces when the tab key is hit.  Also, many editors will automatically
 	     search/replace leading tabs with 4 spaces.
 
 	   o Only 80 characters on a line.
 
-	   o use all lowercase function names with underscore separated words:
+	   Xo use all lowercase function names with underscore separated words:
 
 	        def set_some_value()
 
@@ -30,7 +60,7 @@
 	        def setSomeValue()
 
 
-	   o use CamelCase class names:
+	   Xo use CamelCase class names:
 
 	        def BaseClass()
 
@@ -38,7 +68,7 @@
 
 	        def base_class()
 
-*. make_array_safe
+X*. make_array_safe
    I would prefer 'make_array_safe' named atleast_1d_and_contiguous().  This is more specific
    and it is immediately clear to other developers what the check does.  If you add other checks,
    then perhaps come up with a more generic name, but this explicit names, when possible, help
@@ -47,42 +77,11 @@
    Also, numpy has an ascontiguousrarray() function that would simplify the code a line or so.
 
    Also, this function lives in multiple places, interpolate1d and interpolate_wrapper.
-	
-*. Function signatures:
-    
-    def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
-                        kindkw={}, lowkw={}, highkw={}, \
-                        remove_bad_data = False, bad_data=[], interp_axis = 0):
 
-     o. the trailing \ for new lines are not necessary for functions since line continuation
-        is implicit with the open/close parentheses.
-     o. [minor] I would just use NaN instead of np.NaN.  						
-     o. It is dangerous to initialize containers as default arguments because they are
-        effectively a singleton for the function.  We can discuss this if you don't know
-        what I am talking about.
-     o. kindkw, lowkw, and highkw aren't really necessary I don't think.  They should go.
-	 o. Do you need both remove_bad_data and bad_data?  If bad_data is None, then you 
-	    don't remove bad_data...
-	 o. I think I would change interp_axis to just be axis.  This is consistent with many
-	    other functions in numpy.
-	 o. The choice of whether axis=0 or axis=-1 by default is a reasonable question.
-	    fft defaults to axis=-1.  This is also the faster axis to operate across in the
-	    standard case.  It is, however, the opposite of how some people think about 
-	    things (columns vs. rows).  Talk to Travis O. for his take.  Mine is to use axis=-1.
-	 o. 
-	
-	I think all of this might simplify the interface to the following:
-	
-    def interp1d(x, y, new_x, kind='linear', low=NaN, high=NaN, bad_data=None, axis=-1):
-
-*. isinstance(value, str) should be isinstance(value, basestring) so that we handle
+X*. isinstance(value, str) should be isinstance(value, basestring) so that we handle
    both strings and unicode correctly.
 
-*. _remove_bad_data should be left up to the interpolation method if it knows what to do.
-   otherwise, it is handled by this top level class.
-   And, we definitely don't want the list comprehension in the remove_bad_data class.
-
-*. If the input to interp1d is a scalar, the return value should be a scalar.
+X*. If the input to interp1d is a scalar, the return value should be a scalar.
    [add test to handle this.]
    This fails in the following:
 

Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/fitpack_wrapper.py	2008-07-29 17:56:52 UTC (rev 4577)
@@ -174,74 +174,4 @@
             assert ier==0,`ier`
             return z[:m]
         raise NotImplementedError,\
-              'finding roots unsupported for non-cubic splines'
-                    
-
-
-# testing
-import unittest
-import time
-from numpy import arange, allclose, ones
-
-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
+              'finding roots unsupported for non-cubic splines'
\ No newline at end of file

Modified: branches/Interpolate1D/info.py
===================================================================
--- branches/Interpolate1D/info.py	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/info.py	2008-07-29 17:56:52 UTC (rev 4577)
@@ -9,7 +9,7 @@
     and extrapolation of 1D (in both input and output) real-valued.  The
     primary function provided is:
 
-        interp1d(x, y, new_x) : from data points (x[i], y[i]), interpolates
+        interp1(x, y, new_x) :  from data points (x[i], y[i]), interpolates
                                         values for points in new_x and
                                         returns them as an array.  x and new_x
                                         must both be in sorted order.

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/interpolate1d.py	2008-07-29 17:56:52 UTC (rev 4577)
@@ -1,24 +1,14 @@
-
 # FIXME: information strings giving mathematical descriptions of the actions
 #     of the functions.
 
-from interpolate_wrapper import linear, logarithmic, block, block_average_above
+from interpolate_wrapper import linear, logarithmic, block, block_average_above, atleast_1d_and_contiguous
 from fitpack_wrapper import Spline
 import numpy as np
 from numpy import array, arange, empty, float64, NaN
-
-def make_array_safe(ary, typecode=np.float64):
-    """Used to make sure that inputs and outputs are
-    properly formatted.
-    """
-    ary = np.atleast_1d(np.asarray(ary, typecode))
-    if not ary.flags['CONTIGUOUS']:
-        ary = ary.copy()
-    return ary
     
-def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
-                    kindkw={}, lowkw={}, highkw={}, \
-                    remove_bad_data = False, bad_data=[], interp_axis = 0):
+def interp1d(x, y, new_x, interp = 'linear', low = NaN, high = NaN,
+                    interpkw = {}, lowkw={}, highkw={},
+                    bad_data = None):
     """ A function for interpolation of 1D data.
         
         Parameters
@@ -90,10 +80,15 @@
             >>> interp1d(x, y, new_x)
             array([.2, 2.3, 5.6, NaN])
     """
-    return Interpolate1d(x, y, kind=kind, low=low, high=high, \
-                                    kindkw=kindkw, lowkw=lowkw, highkw=highkw, \
-                                    remove_bad_data = remove_bad_data, bad_data=bad_data\
-                                    )(new_x)
+    return Interpolate1d(x, y, 
+                                interp = interp,
+                                low = low,
+                                high = high,
+                                interpkw = interpkw,
+                                lowkw = lowkw,
+                                highkw = highkw,
+                                bad_data = bad_data
+                                )(new_x)
 
 class Interpolate1d(object):
     """ A callable class for interpolation of 1D, real-valued data.
@@ -192,11 +187,12 @@
         ---------
         
             >>> import numpy
-            >>> from Interpolate1D import interp1d
+            >>> from interpolate1d import Interpolate1d
             >>> x = range(5)        # note list is permitted
             >>> y = numpy.arange(5.)
             >>> new_x = [.2, 2.3, 5.6]
-            >>> interp1d(x, y, new_x)
+            >>> interp_func = Interpolate1d(x, y)
+            >>> interp_fuc(new_x)
             array([.2, 2.3, 5.6, NaN])
     """
     # FIXME: more informative descriptions of sample arguments
@@ -204,19 +200,18 @@
     # FIXME : Allow copying or not of arrays.  non-copy + remove_bad_data should flash 
     #           a warning (esp if we interpolate missing values), but work anyway.
     
-    def __init__(self, x, y, kind='linear', low=np.NaN, high=np.NaN, \
-                        kindkw={}, lowkw={}, highkw={}, \
-                        remove_bad_data = False, bad_data=[]):
+    def __init__(self, x, y, interp = 'linear', low = NaN, high = NaN,
+                        interpkw={}, lowkw={}, highkw={}, bad_data = None):
         # FIXME: don't allow copying multiple times.
         # FIXME : allow no copying, in case user has huge dataset
         
         # remove bad data, is there is any
-        if remove_bad_data:
+        if bad_data is not None:
             x, y = self._remove_bad_data(x, y, bad_data)
         
         # check acceptable size and dimensions
-        x = np.array(x)
-        y = np.array(y)
+        x = np.atleast_1d(x)
+        y = np.atleast_1d(y)
         assert len(x) > 0 and len(y) > 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" 
@@ -227,7 +222,7 @@
         self._init_xy(x, y)
         
         # store interpolation functions for each range
-        self.kind = self._init_interp_method(kind, kindkw)
+        self.interp = self._init_interp_method(interp, interpkw)
         self.low = self._init_interp_method(low, lowkw)
         self.high = self._init_interp_method(high, highkw)
 
@@ -236,10 +231,10 @@
         # select proper dataypes and make arrays
         self._xdtype = {np.float32 : np.float32}.setdefault(type(x[0]), np.float64) # unless data is float32,  cast to float64
         self._ydtype = {np.float32 : np.float32}.setdefault(type(y[0]), np.float64)
-        self._x = make_array_safe(x, self._xdtype).copy()
-        self._y = make_array_safe(y, self._ydtype).copy()
+        self._x = atleast_1d_and_contiguous(x, self._xdtype).copy()
+        self._y = atleast_1d_and_contiguous(y, self._ydtype).copy()
 
-    def _remove_bad_data(self, x, y, bad_data = [None, np.NaN]):
+    def _remove_bad_data(self, x, y, bad_data = [None, NaN]):
         """ removes data points whose x or y coordinate is
             either in bad_data or is a NaN.
         """
@@ -272,7 +267,7 @@
             func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
                         'block_average_above':block_average_above}[interp_arg]
             result = lambda new_x : func(self._x, self._y, new_x, **kw)
-        elif interp_arg in ['Spline', Spline, 'spline']:
+        elif interp_arg in ['Spline', 'spline']:
             # use the Spline class from fitpack_wrapper
             # k = 3 unless otherwise specified
             result = Spline(self._x, self._y, **kw)
@@ -289,7 +284,7 @@
                 result = Spline(self._x, self._y, k=4)
             elif interp_arg in ['Quintic', 'quintic', 'Quin', 'quin']:
                 result = Spline(self._x, self._y, k=5)
-        elif isinstance(interp_arg, str):
+        elif isinstance(interp_arg, basestring):
             raise TypeError, "input string %s not valid" % interp_arg
         
         # secondary usage : user passes a callable class
@@ -314,9 +309,7 @@
                 
         # user passes a function to be called
         # Assume function has form of f(x, y, newx, **kw)
-        # FIXME : should other function forms be allowed?
         elif isfunction(interp_arg):
-            # assume x, y and newx are all passed to interp_arg
             result = lambda new_x : interp_arg(self._x, self._y, new_x, **kw)
         
         # default : user has passed a default value to always be returned
@@ -332,171 +325,39 @@
             Breaks x into pieces in-range, below-range, and above range.
             Performs appropriate operation on each and concatenates results.
         """
-        # FIXME : make_array_safe may also be called within the interpolation technique.
+        # FIXME : atleast_1d_and_contiguous may also be called within the interpolation technique.
         #   waste of time, but ok for the time being.
-        newx = make_array_safe(newx)
         
+        # if input is scalar or 0-dimemsional array, output will be scalar
+        input_is_scalar = np.isscalar(newx) or (isinstance(newx, type(np.array([1.0]))) and np.shape(newx) == ())
+        
+        newx_array = atleast_1d_and_contiguous(newx)
+        
         # masks indicate which elements fall into which interpolation region
-        low_mask = newx<self._x[0]
-        high_mask = newx>self._x[-1]
+        low_mask = newx_array<self._x[0]
+        high_mask = newx_array>self._x[-1]
         interp_mask = (~low_mask) & (~high_mask)
-        
+                
         # use correct function for x values in each region
-        if len(newx[low_mask]) == 0: new_low=np.array([])  # FIXME : remove need for if/else.
+        if len(newx_array[low_mask]) == 0: new_low=np.array([])  # FIXME : remove need for if/else.
                                                                             # if/else is a hack, since vectorize is failing
                                                                             # to work on lists/arrays of length 0
                                                                             # on the computer where this is being
                                                                             # developed
-        else: new_low = self.low(newx[low_mask])
-        if len(newx[interp_mask])==0: new_interp=np.array([])
-        else: new_interp = self.kind(newx[interp_mask])
-        if len(newx[high_mask]) == 0: new_high = np.array([])
-        else: new_high = self.high(newx[high_mask])
+        else: new_low = self.low(newx_array[low_mask])
+        if len(newx_array[interp_mask])==0: new_interp=np.array([])
+        else: new_interp = self.interp(newx_array[interp_mask])
+        if len(newx_array[high_mask]) == 0: new_high = np.array([])
+        else: new_high = self.high(newx_array[high_mask])
         
-        result = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
+        result_array = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
                                                                                           # Would be nice to say result = zeros(dtype=?)
                                                                                           # and fill in
         
-        return result
+        if input_is_scalar:
+            result = float(result_array)
+        else:
+            result = result_array
         
-# unit testing
-import unittest, time
-class Test(unittest.TestCase):
-    
-    def assertAllclose(self, x, y):
-        self.assert_(np.allclose(make_array_safe(x), make_array_safe(y)))
-        
-    def test_interpolate_wrapper(self):
-        """ run unit test contained in interpolate_wrapper.py
-        """
-        #print "\n\nTESTING _interpolate_wrapper MODULE"
-        from 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 fitpack_wrapper import Test
-        T = Test()
-        T.runTest()
-        
-    def test_instantiationFormat(self):
-        """ make sure : all allowed instantiation formats are supported
-        """
-        
-        # 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(k=2), high=Spline(k=2))
-        new_x = np.arange(N+1)-0.5
-        new_y = interp_func(new_x)
-        self.assertAllclose(new_x, new_y)
-        
-    def test_callFormat(self):
-        """ make sure : all allowed calling formats are supported
-        """
-        # 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_removeBad(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] = None; y[0]=np.NaN
-        new_x = arange(N+1)-0.5
-        new_y = interp1d(x, y, new_x, kind='linear', low='linear', high='linear', \
-                                    remove_bad_data = True, bad_data = [None])
-        self.assertAllclose(new_x, new_y)
-        
-    def test_intper1d(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_defaultExt(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', kindkw={'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', kindkw={'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', kindkw={'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)
-        
-        
-if __name__ == '__main__':
-    unittest.main()                 
\ No newline at end of file
+        return result
+  
\ No newline at end of file

Modified: branches/Interpolate1D/interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/interpolate_wrapper.py	2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/interpolate_wrapper.py	2008-07-29 17:56:52 UTC (rev 4577)
@@ -6,13 +6,9 @@
 import sys
 import _interpolate # C extension.  Does all the real work.
 
-def make_array_safe(ary, typecode = np.float64):
-    ary = np.atleast_1d(np.asarray(ary, typecode))
-    if not ary.flags['CONTIGUOUS']:
-        ary = ary.copy()
-    return ary
+def atleast_1d_and_contiguous(ary, typecode = np.float64):
+    return np.atleast_1d( np.ascontiguousarray(ary, typecode) )
 
-
 def linear(x, y, new_x):
     """ Linearly interpolates values in new_x based on the values in x and y
 
@@ -25,9 +21,9 @@
         new_x
             1-D array
     """
-    x = make_array_safe(x, np.float64)
-    y = make_array_safe(y, np.float64)
-    new_x = make_array_safe(new_x, np.float64)
+    x = atleast_1d_and_contiguous(x, np.float64)
+    y = atleast_1d_and_contiguous(y, np.float64)
+    new_x = atleast_1d_and_contiguous(new_x, np.float64)
 
     assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
     if len(y.shape) == 2:
@@ -52,9 +48,9 @@
         new_x
             1-D array
     """
-    x = make_array_safe(x, np.float64)
-    y = make_array_safe(y, np.float64)
-    new_x = make_array_safe(new_x, np.float64)
+    x = atleast_1d_and_contiguous(x, np.float64)
+    y = atleast_1d_and_contiguous(y, np.float64)
+    new_x = atleast_1d_and_contiguous(new_x, np.float64)
 
     assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
     if len(y.shape) == 2:
@@ -80,9 +76,9 @@
             1-D array
     """
     bad_index = None
-    x = make_array_safe(x, np.float64)
-    y = make_array_safe(y, np.float64)
-    new_x = make_array_safe(new_x, np.float64)
+    x = atleast_1d_and_contiguous(x, np.float64)
+    y = atleast_1d_and_contiguous(y, np.float64)
+    new_x = atleast_1d_and_contiguous(new_x, np.float64)
 
     assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
     if len(y.shape) == 2:
@@ -122,77 +118,4 @@
         # take requires the index array to be an Int
         indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(np.int))
         new_y = np.take(y, indices, axis=-1)
-        return new_y
-
-
-# Unit Test
-import unittest
-import time
-from numpy import arange, allclose, ones, NaN, isnan
-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_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
+        return new_y
\ No newline at end of file



More information about the Scipy-svn mailing list