[Scipy-svn] r4690 - in trunk/scipy/interpolate: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Sep 5 07:55:42 CDT 2008


Author: oliphant
Date: 2008-09-05 07:55:38 -0500 (Fri, 05 Sep 2008)
New Revision: 4690

Added:
   trunk/scipy/interpolate/interpolate_wrapper.py
   trunk/scipy/interpolate/tests/test_interpolate_wrapper.py
Log:
Add interpolate wrapper and tests.

Added: trunk/scipy/interpolate/interpolate_wrapper.py
===================================================================
--- trunk/scipy/interpolate/interpolate_wrapper.py	2008-09-05 12:50:12 UTC (rev 4689)
+++ trunk/scipy/interpolate/interpolate_wrapper.py	2008-09-05 12:55:38 UTC (rev 4690)
@@ -0,0 +1,138 @@
+""" helper_funcs.py.
+    scavenged from enthought,interpolate
+"""
+
+import numpy as np
+import sys
+import _interpolate # C extension.  Does all the real work.
+
+def atleast_1d_and_contiguous(ary, dtype = np.float64):
+    return np.atleast_1d( np.ascontiguousarray(ary, dtype) )
+
+def nearest(x, y, new_x):
+    """ Rounds each new_x[i] to the closest value in x
+        and returns corresponding y.
+    """
+    shifted_x = np.concatenate(( np.array([x[0]-1]) , x[0:-1] ))
+    
+    midpoints_of_x = atleast_1d_and_contiguous( .5*(x + shifted_x) )
+    new_x = atleast_1d_and_contiguous(new_x)
+    
+    TINY = 1e-10
+    indices = np.searchsorted(midpoints_of_x, new_x+TINY)-1
+    indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(np.int))
+    new_y = np.take(y, indices, axis=-1)
+    
+    return new_y
+    
+    
+
+def linear(x, y, new_x):
+    """ Linearly interpolates values in new_x based on the values in x and y
+
+        Parameters
+        ----------
+        x
+            1-D array
+        y
+            1-D or 2-D array
+        new_x
+            1-D array
+    """
+    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:
+        new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
+        for i in range(len(new_y)): # for each row
+            _interpolate.linear_dddd(x, y[i], new_x, new_y[i])
+    else:
+        new_y = np.zeros(len(new_x), np.float64)
+        _interpolate.linear_dddd(x, y, new_x, new_y)
+
+    return new_y
+
+def logarithmic(x, y, new_x):
+    """ Linearly interpolates values in new_x based in the log space of y.
+
+        Parameters
+        ----------
+        x
+            1-D array
+        y
+            1-D or 2-D array
+        new_x
+            1-D array
+    """
+    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:
+        new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
+        for i in range(len(new_y)):
+            _interpolate.loginterp_dddd(x, y[i], new_x, new_y[i])
+    else:
+        new_y = np.zeros(len(new_x), np.float64)
+        _interpolate.loginterp_dddd(x, y, new_x, new_y)
+
+    return new_y
+    
+def block_average_above(x, y, new_x):
+    """ Linearly interpolates values in new_x based on the values in x and y
+
+        Parameters
+        ----------
+        x
+            1-D array
+        y
+            1-D or 2-D array
+        new_x
+            1-D array
+    """
+    bad_index = None
+    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:
+        new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
+        for i in range(len(new_y)):
+            bad_index = _interpolate.block_averave_above_dddd(x, y[i], 
+                                                            new_x, new_y[i])
+            if bad_index is not None:
+                break                                                
+    else:
+        new_y = np.zeros(len(new_x), np.float64)
+        bad_index = _interpolate.block_average_above_dddd(x, y, new_x, new_y)
+
+    if bad_index is not None:
+        msg = "block_average_above cannot extrapolate and new_x[%d]=%f "\
+              "is out of the x range (%f, %f)" % \
+              (bad_index, new_x[bad_index], x[0], x[-1])
+        raise ValueError, msg
+              
+    return new_y
+
+def block(x, y, new_x):
+        """ Essentially a step function.
+        
+            For each new_x[i], finds largest j such that
+            x[j] < new_x[j], and returns y[j].
+        """
+        # find index of values in x that preceed values in x
+        # This code is a little strange -- we really want a routine that
+        # returns the index of values where x[j] < x[index]
+        TINY = 1e-10
+        indices = np.searchsorted(x, new_x+TINY)-1
+
+        # If the value is at the front of the list, it'll have -1.
+        # In this case, we will use the first (0), element in the array.
+        # 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
\ No newline at end of file

Added: trunk/scipy/interpolate/tests/test_interpolate_wrapper.py
===================================================================
--- trunk/scipy/interpolate/tests/test_interpolate_wrapper.py	2008-09-05 12:50:12 UTC (rev 4689)
+++ trunk/scipy/interpolate/tests/test_interpolate_wrapper.py	2008-09-05 12:55:38 UTC (rev 4690)
@@ -0,0 +1,86 @@
+""" module to test interpolate_wrapper.py
+"""
+
+# Unit Test
+import unittest
+import time
+from numpy import arange, allclose, ones, NaN, isnan
+import numpy as np
+
+# functionality to be tested
+from scipy.interpolate.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()
+    



More information about the Scipy-svn mailing list