[Scipy-svn] r4652 - in branches/Interpolate1D: . docs tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Aug 19 15:02:05 EDT 2008
Author: fcady
Date: 2008-08-19 14:02:03 -0500 (Tue, 19 Aug 2008)
New Revision: 4652
Modified:
branches/Interpolate1D/__init__.py
branches/Interpolate1D/docs/tutorial.rst
branches/Interpolate1D/interpolate1d.py
branches/Interpolate1D/interpolate2d.py
branches/Interpolate1D/setup.py
branches/Interpolate1D/tests/test_interpolate2d.py
Log:
minor changes, some in preparation for ND scattered interpolation
Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/__init__.py 2008-08-19 19:02:03 UTC (rev 4652)
@@ -18,3 +18,5 @@
# wrapper around Fortran implementing Tom's algorithm 526
from algorithm526_wrapper import algorithm526
+from interpSNd import InterpolateSNd
+
Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/docs/tutorial.rst 2008-08-19 19:02:03 UTC (rev 4652)
@@ -166,9 +166,9 @@
interpolated.
---------------------------------------
+-----------------------------------------------------------
User-defined Interpolation Methods
---------------------------------------
+-----------------------------------------------------------
The string interface is designed to conveniently take care of most things a user would want
to do in a way that is easy and, when something goes wrong, informative and helpful.
@@ -579,8 +579,8 @@
newz = interp2d(x, y, z, newx, newy, kind='linear', out=NaN)
-where x, y, z, are 1D arrays or lists, and newx and newy may be either arrays, lists or scalars.
-If they are scalars or zero-dimensional arrays, newz will be a scalar as well. Otherwise
+where x, y, z, are arrays (1D or 2D) or lists, and newx and newy may be either arrays, lists or scalars.
+If newx and newy are scalars or zero-dimensional arrays, newz will be a scalar as well. Otherwise
a vector is returned. The only differences from intper1d are
#) The known data points are specified by 3 arrays (x, y and z) rather than 2 (x and y).
@@ -737,8 +737,22 @@
ND Scattered Interpolation
================================================
- Still in development.
+Still in development.
- Ideally the range of interpolation would be the convex hull of the known
- data points, and a Delaunay tesselation would be determined and stored
- at instantiation. Then again, that would be very expensive.
\ No newline at end of file
+Ideally the range of interpolation would be the convex hull of the known
+data points, and a Delaunay tesselation would be determined and stored
+at instantiation. Then again, that would be very expensive.
+
+InterpolateNd suffers from the requirement that data points be on a uniformly
+spaced grid. This problem is solved by the callable class InterpolateSNd and
+the function interpSNd, which interpolate *scattered* N-dimensional data.
+
+First off a warning. The code is in a rather preliminary stage which, while functioning,
+is VERY slow and not extensively tested. This state of affairs is still being improved.
+
+The valid interpolation range is the convex hull of the data points.
+
+
+
+
+This functionality is still preliminary.
\ No newline at end of file
Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/interpolate1d.py 2008-08-19 19:02:03 UTC (rev 4652)
@@ -164,7 +164,7 @@
interpolate from. Note that 2-dimensional
y is not supported.
- newx -- list of 1D numpy array
+ newx (when calling) -- list of 1D numpy array
x values at which to interpolate the value
of the function
Modified: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/interpolate2d.py 2008-08-19 19:02:03 UTC (rev 4652)
@@ -36,26 +36,137 @@
# functional interface: creates and calls an instance of objective interface
def interp2d(x, y, z, newx, newy, kind='linear', out=NaN, bad_data=None):
+ """ A function for interpolation of 2D, real-valued data.
+
+ Parameters
+ -----------
+
+ x -- list or NumPy array
+ x includes the x-values for the data set to
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
+
+ y -- list or NumPy array
+ y includes the y-values for the data set to
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
+
+ z -- list or 1D NumPy array
+ z includes the z-values for the data set to
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
+
+ newx -- list or NumPy array
+ newx includes the x-values for the points at
+ which to perform interpolation. If it is in array,
+ it must be 1 or 2-dimensional
+
+ newy -- list or NumPy array
+ newy includes the y-values for the points at
+ which to perform interpolation. If it is in array,
+ it must be 1 or 2-dimensional
+
+ 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
+ "526" -- TOMS algorithm 526
+
+ 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_func(newx, newy)
+ array([1.2, 3.3, 3.6, NaN])
+ """
return Interpolate2d(x, y, z, kind=kind, out=out, bad_data=bad_data)(newx, newy)
# objective interface
class Interpolate2d:
- """ A callable class for interpolation of 1D, real-valued data.
+ """ A callable class for interpolation of 2D, real-valued data.
Parameters
-----------
- x -- list or 1D NumPy array
+ x -- list or NumPy array
x includes the x-values for the data set to
- interpolate from.
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
- y -- list or 1D NumPy array
+ y -- list or NumPy array
y includes the y-values for the data set to
- interpolate from.
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
z -- list or 1D NumPy array
z includes the z-values for the data set to
- interpolate from.
+ interpolate from. If it is in array, it must be
+ 1 or 2-dimensional
Optional Arguments
-------------------
@@ -91,6 +202,7 @@
"cubic" -- spline interpolation order 3
"quartic" -- spline interpolation order 4
"quintic" -- spline interpolation order 5
+ "526" -- TOMS algorithm 526
Other options for kind and out
---------------------------------------------------
@@ -131,9 +243,8 @@
>>> newx = [.2, 2.3, 2.6, 7.0]
>>> newy = [1, 1, 1, 1]
>>> interp_func = Interpolate2d(x, y, z)
- >>> interp_fuc(newx, newy)
+ >>> interp_func(newx, newy)
array([1.2, 3.3, 3.6, NaN])
-
"""
def __init__(self, x, y, z, kind='linear', out=NaN, bad_data=None):
@@ -148,15 +259,15 @@
# FIXME : perhaps allow 2D input if it is inthe form of meshgrid
# check acceptable sizes and dimensions
+ # and make data 1D
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"
+ assert x.shape == y.shape, "x and y must be the same shape"
+ assert x.shape == z.shape, "z must be the same shape as x and y"
+ assert x.ndim in [1, 2], "x and y must be at most 2-dimensional"
+ assert np.size(x) >= 0, "x, y and z cannot be of size zero"
+ x, y, z = map(np.ravel, [x, y, z])
# remove bad data if applicable
if bad_data is not None:
@@ -246,18 +357,22 @@
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)
- newy = atleast_1d_and_contiguous(newy, dtype=self._ydtype)
- assert newx.ndim == 1, "newx can be at most 1-dimensional"
- assert newy.ndim == 1, "newy can be at most 1-dimensional"
- assert len(newx) == len(newy), "newx and newy must be the same length"
+ # check acceptable sizes and dimensions, record the shape of the input,
+ # and make data 1D
+ newx = np.atleast_1d(newx)
+ newy = np.atleast_1d(newy)
+ assert newx.shape == newy.shape, "newx and newy must be the same shape"
+ assert newx.ndim in [1, 2], "newx and newy must be at most 2-dimensional"
+ assert np.size(newx) >= 0, "newx and newy must be of size greater than zero"
+ shape_of_newx = newx.shape
+ newx, newy= map(np.ravel, [newx, newy])
+ # finding which points are in the valid interpolation range
in_range_mask = (min(self._x) <= newx) & (newx <= max(self._x)) & \
(min(self._y) <= newy) & (newy <= max(self._y))
- # filling array of interpolated z-values
- result = np.zeros(np.shape(newx), dtype = self._zdtype)
+ # making array of interpolated z-values
+ result = np.empty(np.shape(newx), dtype = self._zdtype)
if sum(in_range_mask) > 0: # if there are in-range values. hack to deal
# with behavior of np.vectorize on arrays of length 0
result[in_range_mask] = self.kind(newx[in_range_mask], newy[in_range_mask])
@@ -267,5 +382,7 @@
# revert to scalar if applicable
if input_is_scalar:
result = result[0]
+ else:
+ result = np.reshape(result, shape_of_newx)
return result
\ No newline at end of file
Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/setup.py 2008-08-19 19:02:03 UTC (rev 4652)
@@ -46,10 +46,12 @@
'extensions/interp_526.f']
)
+ # include tutorial on the module
config.add_data_dir('docs')
return config
if __name__ == '__main__':
+ #from setuptools import setup
from numpy.distutils.core import setup
setup(**configuration().todict())
\ No newline at end of file
Modified: branches/Interpolate1D/tests/test_interpolate2d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate2d.py 2008-08-18 20:58:28 UTC (rev 4651)
+++ branches/Interpolate1D/tests/test_interpolate2d.py 2008-08-19 19:02:03 UTC (rev 4652)
@@ -17,7 +17,72 @@
def assertAllclose(self, x, y):
self.assert_(np.allclose(atleast_1d_and_contiguous(x), atleast_1d_and_contiguous(y)))
+
+ 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)
+
+ def test_string_526(self):
+ """ make sure : keyword '526' works
+ ie that TOMS algorithm 526 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='526', out='526')
+ newz = interp_func(newx, newy)
+
+ self.assertAllclose(newz, newx+newy)
+
def test_callable_class_interpolation(self):
""" make sure : instance of callable class with xyz not initiated works
"""
@@ -126,70 +191,21 @@
self.assert_( np.isnan(newz[0]) )
self.assert_( np.isnan(newz[-1]) )
- def test_string_linear(self):
- """ make sure : string 'linear' works
+
+ def test_2D_input(self):
+ """ make sure : newx and newy can be 2D, and the result is also 2D
"""
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
+ newx, newy = meshgrid(arange(1, N)-.2, arange(1,N)-.2)
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)
-
- def test_string_526(self):
- """ make sure : keyword '526' works
- ie that TOMS algorithm 526 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='526', out='526')
- newz = interp_func(newx, newy)
-
- self.assertAllclose(newz, newx+newy)
-
if __name__ == '__main__':
unittest.main()
\ No newline at end of file
More information about the Scipy-svn
mailing list