[Scipy-svn] r4690 - in trunk/scipy/interpolate: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Sep 5 08:55:42 EDT 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