[Scipy-svn] r4608 - in branches/Interpolate1D: . docs tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Aug 7 15:31:47 EDT 2008
Author: fcady
Date: 2008-08-07 14:31:44 -0500 (Thu, 07 Aug 2008)
New Revision: 4608
Modified:
branches/Interpolate1D/docs/tutorial.rst
branches/Interpolate1D/ndimage_wrapper.py
branches/Interpolate1D/tests/test_ndimage.py
Log:
more exhaustive tests (which resulted in fixing a few bugs that had been there) and more-or-less complete documentation
Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst 2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/docs/tutorial.rst 2008-08-07 19:31:44 UTC (rev 4608)
@@ -645,8 +645,57 @@
ND Interpolation
================================================
+1D and 2D interpolation are analogous in their user interface. However,
+the ND interpolation breaks from their patterns in certain important respects.
+The biggest difference is that the known data must be on a uniformly spaced
+grid; scattered points are not acceptable. Rather than a list of points, the
+user passes in an N-dimensional array of data. By default, the (i, j, k) element
+of the array is taken to give the value of the function at the point (i, j, k), but
+the user can specify a different set of starting points and spacings. The
+starting points and spacings can be different for each axis. Starting points
+and spacings are specified by the keywords 'starting_coords' and 'spacings',
+which can be passed as arrays or lists.
+Second, when calling the object / function, the points to be interpolated are
+passed as an nxL array, where n is the dimensionality of the data; each column
+of the array specifies a point at which to interpolate a value, and the returned
+array of length L gives the value at each of the points.
+
+Here is a basic example which illustrates these points: ::
+
+ In []: from interpolate import interpNd
+ In []: X, Y = meshgrid(arange(10.), arange(10.))
+ In []: Z = X+Y # function just adds coordinates
+ In []: coordinates = array([ [1.1, 2.2, 4.6],
+ [1.1, 2.2, 4.0 ])
+ In []: interpNd(Z, coordinates)
+ Out []: array([2.2, 4.4, 8.6])
+ # say the data start at the point (2, 1) and
+ # points are spaced 2 apart
+ In []: interpNd(Z, coordinates, starting_coords = array([1, 1], spacings=[2, 2])
+ Out []: array([ .1, 1.2, 3.3])
+
+By default, the interpolation is linear for points in-range and returns NaN for
+out-of-bounds; alternate types can be specified by the keywords
+'kind' and 'out', as in 2D interpolation. However,
+
+#) 'kind' must be a string ('linear', 'block', 'cubic', etc) indicating a type of
+ spline interpolation, or else an integers specifying the spline order.
+#) 'out' must be either NaN (the default), 'nearest', 'wrap', 'reflect' or 'constant'
+
+The user cannot pass in specially-tailored interpolation methods.
+
+There
+
+The second point is that all interpolation is done using splines. The keyword is still
+"kind", but only keywords specifying splines ('spline', 'cubic', 'quadratic', 'quintic', etc)
+are acceptable. If kind is an integer, that integer is taken to be the order of the spline.
+Note also that 'linear' denotes a spline of order 1, and 'block' denotes a spline of order
+zero.
+
+
+
================================================
ND Scattered Interpolation
================================================
Modified: branches/Interpolate1D/ndimage_wrapper.py
===================================================================
--- branches/Interpolate1D/ndimage_wrapper.py 2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/ndimage_wrapper.py 2008-08-07 19:31:44 UTC (rev 4608)
@@ -4,21 +4,28 @@
import numpy as np
import _nd_image
+def interpNd(data, coordinates, starting_coords=None, spacings=None, kind='linear',out=NaN):
+ return Interpolate1d(data = data,
+ starting_coords = starting_coords,
+ spacings = spacings,
+ kind = kind,
+ out = out
+ )(coordinates)
+
class InterpolateNd:
def __init__(self, data, starting_coords =None, spacings = None,
- order=3, out=NaN):
+ kind='linear', out=NaN):
""" data = array or list of lists
starting_coords = None, list, 1D array or 2D (nx1) array
spacings = None, list, 1D array or 2D (nx1) array
+ kind = string or integer
+ 0 = block extrapolation between midpoints
out = string in 'nearest', 'wrap', 'reflect', 'mirror', 'constant'
or just NaN
"""
# FIXME : include spline filtering
- if order < 0 or order > 5:
- raise RuntimeError, 'spline order not supported'
-
# checking format of input
data = array(data)
@@ -29,15 +36,54 @@
starting_coords = array(starting_coords)
assert starting_coords.size == data.ndim, "There must be one element of \
starting_coords per data dimension. Size mismatch."
- starting_coords = reshape(starting_coords, (data.ndim, 1))
+ starting_coords = np.reshape(starting_coords, (data.ndim, 1))
if spacings == None:
- spacings = np.zeros(( data.ndim, 1 ))
+ spacings = np.ones(( data.ndim, 1 ))
else:
spacings = array(spacings)
assert starting_coords.size == data.ndim, "There must be one element of \
starting_coords per data dimension"
- spacings = reshape(spacings, (data.ndim, 1))
+ spacings = np.reshape(spacings, (data.ndim, 1))
+ # determining the order
+ order_dict = \
+ { 0:0,
+ '0':0,
+ 'block':0,
+ 1:1,
+ '1':1,
+ 'linear':1,
+ 'Linear':1,
+ 2:2,
+ '2':2,
+ 'quadratic':2,
+ 'quad':2,
+ 'Quadratic':2,
+ 'Quad':2,
+ 3:3,
+ '3':3,
+ 'spline':3,
+ 'Spline':3,
+ 'cubic':3,
+ 'Cubic':3,
+ 4:4,
+ '4':4,
+ 'quartic':4,
+ 'Quartic':4,
+ 5:5,
+ '5':5,
+ 'quintic':5,
+ 'quint':5,
+ 'Quintic':5,
+ 'Quint':5
+ }
+ if order_dict.has_key(kind):
+ self.order = order_dict[kind]
+ elif isinstance(kind, int):
+ raise ValueError, "Only spline orders 0, 1, ..., 5 are supported"
+ else:
+ raise ValueError, "argument kind = %s not recognized" % str(kind)
+
# storing relevant data
self._data_array = data
self.ndim = data.ndim
@@ -46,7 +92,6 @@
self._min_coords = starting_coords
self._max_coords = self._min_coords + self._shape*self._spacings
self.out = out
- self.order = order
def __call__(self, coordinates):
""" coordinates is an n x L array, where n is the dimensionality of the data
Modified: branches/Interpolate1D/tests/test_ndimage.py
===================================================================
--- branches/Interpolate1D/tests/test_ndimage.py 2008-08-06 23:43:17 UTC (rev 4607)
+++ branches/Interpolate1D/tests/test_ndimage.py 2008-08-07 19:31:44 UTC (rev 4608)
@@ -7,7 +7,7 @@
import unittest
import time
-from numpy import arange, allclose, ones
+from numpy import arange, allclose, ones, array
import numpy as np
import ndimage_wrapper as nd
@@ -17,40 +17,97 @@
self.assert_(np.allclose(x, y))
def test_linear(self):
+ """ Make sure : basic linear works
+ """
boring_data = np.ones((5,5,5))
- interp = nd.InterpolateNd(boring_data, order = 1)
+ interp = nd.InterpolateNd(boring_data, kind = 'linear')
self.assertAllclose( interp(np.array([[2.3], [1.0], [3.9]])) , 1.0 )
+ def test_linear_not_1(self):
+ """ Make sure : linear interpolation works on a general dataset
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, kind = 'linear')
+ self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 )
+
def test_data_is_list(self):
+ """ Make sure : data can be entered as a list
+ """
boring_data = [ [1.0, 1.0, 1.0],
[1.0, 1.0, 1.0],
[1.0, 1.0, 1.0]]
- interp = nd.InterpolateNd(boring_data, order = 1)
+ interp = nd.InterpolateNd(boring_data)
self.assertAllclose( interp(np.array([[1.3], [1.0]])) , 1.0 )
def test_coords_is_1d(self):
+ """ Make sure : coordinates for a single point can be entered as a 1D array
+ """
boring_data = np.ones((5,5,5))
- interp = nd.InterpolateNd(boring_data, order = 1)
+ interp = nd.InterpolateNd(boring_data)
self.assertAllclose( interp(np.array([2.3, 1.0, 3.9])) , 1.0 )
def test_coords_is_list(self):
+ """ Make sure : coordinates for a single point can be entered as a list
+ """
boring_data = np.ones((5,5,5))
- interp = nd.InterpolateNd(boring_data, order = 1)
+ interp = nd.InterpolateNd(boring_data)
self.assertAllclose( interp([2.3, 1.0, 3.9]) , 1.0 )
+
+ def test_order2(self):
+ """ Make sure : quadratic interpolation works
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, kind = 2)
+ self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 3.3 )
+ def test_order0(self):
+ """ Make sure : block interpolation works
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, kind = 0)
+ self.assertAllclose( interp(np.array([[2.3], [1.1]])) , 3.0 )
+
+ def test_order3(self):
+ """ Make sure : cubic interpolation works
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, kind = 3)
+ self.assertAllclose( interp(np.array([[4.3], [4.1]])) , 8.4 )
+
+ def test_out(self):
+ """ Make sure : out-of-bounds returns NaN
+ """
+ boring_data = np.ones((5,5,5))
+ interp = nd.InterpolateNd(boring_data, kind = 'linear')
+ self.assert_( np.isnan(interp( np.array([[7.3], [1.0], [3.9]]) )))
+
+ def test_starting_coords(self):
+ """ Make sure : non-zero starting coordinates work correctly
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, starting_coords = array([2, 1]))
+ self.assertAllclose( interp(np.array([[2.3], [1.0]])) , 0.3 )
+
+ def test_spacings(self):
+ """ Make sure : spacings other than 1 work correctly
+ """
+ X, Y = np.meshgrid(arange(10.), arange(10.))
+ interesting_data = X+Y
+ interp = nd.InterpolateNd(interesting_data, spacings = array([2, 1]))
+ self.assertAllclose( interp(np.array([[2.4], [1.0]])) , 2.2 )
+
def runTest(self):
+ """ run all tests
+ """
test_list = [method_name for method_name in dir(self) if method_name.find('test')==0]
for test_name in test_list:
exec("self.%s()" % test_name)
-
- def test_order2(self):
- boring_data = np.ones((5,5,5))
- interp = nd.InterpolateNd(boring_data, order = 2)
- self.assertAllclose( interp(np.array([[2.3], [1.0], [3.9]])) , 1.0 )
- def test_out(self):
- pass
-
if __name__ == '__main__':
unittest.main()
\ No newline at end of file
More information about the Scipy-svn
mailing list