[Scipy-svn] r4595 - in branches/Interpolate1D: . docs tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Aug 1 19:34:04 EDT 2008
Author: fcady
Date: 2008-08-01 18:34:03 -0500 (Fri, 01 Aug 2008)
New Revision: 4595
Added:
branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
Modified:
branches/Interpolate1D/docs/tutorial.rst
branches/Interpolate1D/fitpack_wrapper2d.py
branches/Interpolate1D/tests/test_fitpack_wrapper.py
Log:
more work and documentation for 2D interpolation
Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst 2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/docs/tutorial.rst 2008-08-01 23:34:03 UTC (rev 4595)
@@ -471,14 +471,12 @@
At instantiation:
-#) bbox
+#. bbox
This is a 2-element list specifying the endpoints of the approximation interval.
It default to [x[0],x[-1]]
-
-#) w
+#. w
a 1D sequence of weights which defaults to all ones.
-
-#) s
+#. s
If s is zero, the interpolation is exact. If s is not 0, the curve is smoothe subject to
the constraint that sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
@@ -487,27 +485,30 @@
At calling:
-#) nu
+#. nu
Spline returns, not the spline function S, but the (nu)th derivative of S. nu defaults
to 0, so Spline usually returns the zeroth derivative of S, ie S.
+
-----------------
Special Methods
-----------------
-#) set_smoothing_factor(s)
-#) get_knots
+The following special methods are also available, which are not wrapped by Interpolate1d.
+
+#. set_smoothing_factor(s = 0.0)
+#. get_knots
returns the positions of the knots of the spline
-#) get_coeffs
+#. get_coeffs
returns the coefficients of the
-#) get_residual
+#. get_residual
returns the weighted sum of the errors (due to smoothing) at the data points
sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0)
-#) integral(a, b)
+#. integral(a, b)
returns the integral from a to b
-#) derivatives(x)
+#. derivatives(x)
returns all the derivatives of the spline at point x
-#) roots
+#. roots
This only works for cubic splines. But it returns the places where the spline
is identically zero.
@@ -550,8 +551,43 @@
The Objective Interface
------------------------------------------
-The objective interface for 2D is slightly more complicated.
+The objective interface for 2D is similarly analogous to 1D.
+------------------------------------------
+The Spline2d Class
+------------------------------------------
+
+Just as with Spline, Spline2d is mostly intended to be wrapper by
+Interpolate2d, but it also functions as a stand-alone class featuring
+functionality not accessible through Interpolate2d.
+
+It is instantiated in virtually the same way ::
+
+ instance = Spline2d(x, y, z, kx = 3, ky = 3, s=0.0)
+
+where x, y and z are 1D arrays. It is called with arrays
+newx and newy, returning an array newz of the same length.
+
+Beyond basic usage, Spline2 also has the methods
+
+#) get_grid(self, x, y)
+x and y are treated as the coordinates of a grid, and all
+points on the grid are interpolated and returned in an array.
+
+That is, if z = S.get_grid(x,y), z[i,j] is the interpolated value
+at the point (xi, yj)
+
+#) integral(xa, xb, ya, yb)
+Integrate the interpolated function over the indicated rectangle
+
+#) get_residual
+
+#) get_knots
+
+#) get_coeffs
+
+
+
================================================
ND Interpolation
================================================
Modified: branches/Interpolate1D/fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper2d.py 2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/fitpack_wrapper2d.py 2008-08-01 23:34:03 UTC (rev 4595)
@@ -50,14 +50,11 @@
[xb,xe] x [yb, ye] calculated from a given set of data points
(x,y,z).
- See also:
-
- bisplrep, bisplev - an older wrapping of FITPACK
- UnivariateSpline - a similar class for univariate spline interpolation
- SmoothUnivariateSpline - to create a BivariateSpline through the
- given points
- LSQUnivariateSpline - to create a BivariateSpline using weighted
- least-squares fitting
+ More commenting needed
+
+ If (xi, yi) is outside the interpolation range, it is
+ assigned the value of the nearest point that is within the
+ interpolation range.
"""
def __init__(self, x=None, y=None, z=None, w=None, bbox=[None]*4, kx=3, ky=3, s=0.0, eps=None):
"""
@@ -118,8 +115,11 @@
def __call__(self, x, y):
""" Evaluate spline at positions x[i],y[i].
x and y should be 1d arrays.
+
+ If (xi, yi) is outside the interpolation range, it will be
+ assigned the value of the nearest point that is within the
+ interpolation range.
"""
- # what happens when x contains duplicate values?
if self._is_initialized is not True:
raise Error, "x, y and z must be initialized before interpolating"
@@ -131,26 +131,24 @@
data_grid = self.get_grid(sorted_x, sorted_y)
# fixme : no list comprehension
- z = np.array([ data_grid[np.searchsorted(sorted(x), x[i]), np.searchsorted(sorted(y),y[i])] \
+ z = np.array([ data_grid[np.searchsorted(sorted_x, x[i]), np.searchsorted(sorted_y,y[i])] \
for i,xi in enumerate(x) ])
-
+
return z
- def get_grid(self, x, y, mth='array'):
+ def get_grid(self, x, y):
""" Evaluate spline at positions x[i],y[j]."""
if self._is_initialized is not True:
raise Error, "x, y and z must be initialized before interpolating"
- if mth=='array':
- tx,ty,c = self.tck[:3]
- kx,ky = self.degrees
- z,ier = _dfitpack.bispev(tx,ty,c,kx,ky,x,y)
- assert ier==0,'Invalid input: ier='+`ier`
- return z
- raise NotImplementedError
-
+ tx,ty,c = self.tck[:3]
+ kx,ky = self.degrees
+ z,ier = _dfitpack.bispev(tx,ty,c,kx,ky,x,y)
+ assert ier==0,'Invalid input: ier='+`ier`
+ return z
+
def get_residual(self):
""" Return weighted sum of squared residuals of the spline
approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0)
Modified: branches/Interpolate1D/tests/test_fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper.py 2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper.py 2008-08-01 23:34:03 UTC (rev 4595)
@@ -17,7 +17,7 @@
def assertAllclose(self, x, y):
self.assert_(np.allclose(x, y))
- def test_linearInterp(self):
+ def test_k_1(self):
""" make sure : linear interpolation (spline with order = 1, s = 0)works
"""
N = 3000.
@@ -34,7 +34,7 @@
#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):
+ def test_k_2(self):
""" make sure : quadratic interpolation (spline with order = 2, s = 0)works
"""
N = 3000.
@@ -49,7 +49,19 @@
#print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
self.assertAllclose(new_y, y)
+ def test_extrap(self):
+ """ make sure 1D extrapolation works
+ """
+ N = 3000.
+ x = np.arange(N)
+ y = np.arange(N)
+ interp_func = Spline(x, y, k=1)
+ newx = np.arange(-2, N)+0.5
+ newy = interp_func(newx)
+
+ self.assertAllclose(newx, newy)
+
def test_inputFormat(self):
""" make sure : it's possible to instantiate Spline without x and y
"""
Added: branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper2d.py 2008-08-01 21:48:32 UTC (rev 4594)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper2d.py 2008-08-01 23:34:03 UTC (rev 4595)
@@ -0,0 +1,106 @@
+""" test fitpack_wrapper2d.py
+"""
+
+# hack to test on Field's computer
+import sys
+sys.path.append('c:/home/python/Interpolate1d')
+
+# testing
+import unittest
+import time
+from numpy import arange, allclose, ones, meshgrid, ravel, array
+import numpy as np
+from fitpack_wrapper2d import Spline2d
+
+class Test(unittest.TestCase):
+
+ def assertAllclose(self, x, y):
+ self.assert_(np.allclose(x, y))
+
+ def test_k_1(self):
+ """ make sure : linear interpolation (spline with order = 1, s = 0)works
+ """
+ N = 10.
+ X, Y = meshgrid( arange(N), arange(N) )
+ Z = X + Y
+ x, y, z = map(ravel, [X, Y, Z])
+
+ newx = arange(N-1) +.5
+ newy = newx
+
+ interp_func = Spline2d(x, y, z, kx=1, ky=1)
+ newz = interp_func(newx, newy)
+
+ self.assertAllclose(newz, newx+newy)
+
+ def test_k_2(self):
+ """ make sure : quadratic interpolation (spline with order = 2, s = 0)works
+ """
+ N = 10.
+ X, Y = meshgrid( arange(N), arange(N) )
+ Z = X + Y
+ x, y, z = map(ravel, [X, Y, Z])
+
+ newx = arange(N-1)+0.5
+ newy = newx
+
+ interp_func = Spline2d(x, y, z, kx=2, ky=2)
+ newz = interp_func(newx, newy)
+
+ self.assertAllclose(newz, newx+newy)
+
+
+ def test_instantiate_without_init(self):
+ """ make sure : it's possible to instantiate Spline2d without setting x and y
+ """
+ N = 10.
+ X, Y = meshgrid( arange(N), arange(N) )
+ Z = X + Y
+ x, y, z = map(ravel, [X, Y, Z])
+
+ newx = np.arange(N-1)+0.5
+ newy = newx
+
+ interp_func = Spline2d(kx=1, ky=3)
+ interp_func.init_xyz(x, y, z)
+ newz = interp_func(newx, newy)
+
+ self.assertAllclose(newz, newx+newy)
+
+ def test_extrap(self):
+ """ make sure : linear extrapolation works
+ """
+ N = 10.
+ X, Y = meshgrid( arange(N), arange(N) )
+ Z = X + Y
+ x, y, z = map(ravel, [X, Y, Z])
+
+ newx = arange(N+8) +.5
+ newy = 2*newx
+
+ interp_func = Spline2d(x, y, z, kx=1, ky=1)
+ newz = interp_func(newx, newy)
+
+ print "newx: ", newx
+ print "newy: ", newy
+ print "sum : ", newx+newy
+ print "newz: ", newz
+
+ print "Homer Simpson"
+ print interp_func(array([-2.0]),array([3.5]))
+ print interp_func(array([-7.0]),array([3.5]))
+ print interp_func(array([-2.0]),array([7]))
+ print "Bartman"
+
+ self.assertAllclose(newz, newx+newy)
+
+ 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
More information about the Scipy-svn
mailing list