[Scipy-svn] r4595 - in branches/Interpolate1D: . docs tests

scipy-svn@scip... scipy-svn@scip...
Fri Aug 1 18:34:04 CDT 2008

```Author: fcady
Date: 2008-08-01 18:34:03 -0500 (Fri, 01 Aug 2008)
New Revision: 4595

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).

-
-        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_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
"""

===================================================================
--- 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

```