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

scipy-svn@scip... scipy-svn@scip...
Mon Aug 4 15:27:29 CDT 2008


Author: fcady
Date: 2008-08-04 15:27:24 -0500 (Mon, 04 Aug 2008)
New Revision: 4600

Added:
   branches/Interpolate1D/docs/XY-crossection.png
   branches/Interpolate1D/docs/XY-surface.png
Modified:
   branches/Interpolate1D/TODO.txt
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/fitpack_wrapper2d.py
   branches/Interpolate1D/interpolate2d.py
Log:
more work done on documentation, a format check on the input, and some bug fixes.  This should be rpetty good now.

Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt	2008-08-04 19:48:18 UTC (rev 4599)
+++ branches/Interpolate1D/TODO.txt	2008-08-04 20:27:24 UTC (rev 4600)
@@ -5,7 +5,7 @@
 at appropriate places in the code.
 
 
-************ PRESSING ISSUES ***********
+************ API and/or MAJOR ISSUES ***********
 
 **handle smoothing 
     There is a question of whether, if, and how much to allow
@@ -19,7 +19,11 @@
     This appears resolved, in that Spline can do smoothing (but defaults
     to not), and user-defined classes are allowed to smooth.
 
+**remove list comprehension in interpolate2d.Interpolate2d.__call__
 
+**possibly allow interp2d to return a 2d array?  Like in the case that
+    x and y are 2d arrays in meshgrid format.
+
 **pick best spline
     Under-the-hood machinery currently comes from _interpolate.cpp
     (used in enthought.interpolate) and FITPACK (Fortran, used in 
@@ -31,15 +35,6 @@
     regular grid.  I'm inclined to stay with FITPACK, except for the
     slow performance when x is small (we could add a hack to not
     let s be tiny > 0).
-
-
-**clean up fitpack_wrapper.py
-    Currently it is all hacked together from scipy.interpolate.
-    There is unused functionality, commented-out functionality,
-    and other undesirables.  Once it has been established what
-    we will include, and what should be accessible to the user, 
-    that stuff needs to be cleaned up, and the
-    rest needs to be either removed or set aside.
     
 **allow y to be 2-dimensional?
     It is not decided whether this feature should be supported, but
@@ -69,7 +64,11 @@
     
     See also "allow y to be 2-dimensional?" below
     
+**Put Spline2d into fitpack_wrapper.py
+    It's out now so that is can be worked on separately, but the
+    fitpack-based code should all be in one module.
     
+    
 
 *********** DOCUMENTATION-TYPE AND NON-URGENT TASKS *******
 
@@ -113,7 +112,7 @@
     
 ********* LONGER TERM ************
 
-**update for 2D and ND
+**update for ND
     This will take the form of two additional
     classes both modelled after interpolate1d.  Thus it probably
     shouldn't be done until interpolate1d is more settled.

Added: branches/Interpolate1D/docs/XY-crossection.png
===================================================================
(Binary files differ)


Property changes on: branches/Interpolate1D/docs/XY-crossection.png
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + application/octet-stream

Added: branches/Interpolate1D/docs/XY-surface.png
===================================================================
(Binary files differ)


Property changes on: branches/Interpolate1D/docs/XY-surface.png
___________________________________________________________________
Name: svn:executable
   + *
Name: svn:mime-type
   + application/octet-stream

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-04 19:48:18 UTC (rev 4599)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-04 20:27:24 UTC (rev 4600)
@@ -6,28 +6,30 @@
 a sample of a function into an approximation of that function for every value, and is one of
 the most basic mathematical tools available to a researcher.
 
-The interpolate package provides tools for interpolating and extrapolating new data points from a known set of data points.  
-It provides a functional interface that is flexible and easy to use as well as an object oriented interface that 
-can be more efficient and flexible for some cases.  It is able to interpolate and extrapolate in 1D, 2D, and even N 
-dimensions. *[FIXME : 1D only right now]*
+The interpolate package provides tools for interpolating and extrapolating new data points 
+from a known set of data points.  It provides a functional interface that is flexible and easy 
+to use as well as an object oriented interface that can be more efficient and flexible for some 
+cases.  It is able to interpolate and extrapolate in 1D, 2D, and even N dimensions. *[FIXME : 1D only right now]*
 
-For 1D interpolation, it handles linear and spline(cubic, quadratic, and quintic) for both uniformly and non-uniformly spaced 
-data points "out of the box."  Users can choose the behavior when new values fall outside the range of known data, and
-with a little more work, they can incorporate interpolation methods that are specially tailored to their needs.
+For 1D interpolation, it handles both uniformly and non-uniformly spaced data points, 
+and many popular interpolation methods (particularly splines) are available by simply passing in a keyword
+string ('linear' for linear interpolation, 'cubic' for spline order 3, etc).  With a little 
+more work, users can design and incorporate interpolation methods that are specially tailored to their needs.
 
-For 2D interpolation, *[FIXME : include this]*
+2D interpolation is analogous to 1D, with a nearly identical user interface (corrected for the fact
+that each data point has one more component) and functionality.
 
 This tutorial covers how to use the interpolate module, provides some basic examples, and shows
 them at work in realistic sample sessions.  These sessions demonstrate how to use the 
 interpolate module, but also highlight some of the uses of interpolation techniques.
 
-================================================
-1D Interpolation with the Functional Interface
-================================================
+======================
+1D Interpolation 
+======================
 
--------------
-Basic Usage
--------------
+---------------------------------------------------------
+Basic Usage with the Functional Interface
+---------------------------------------------------------
 
 For most users, the primary feature of interpolate is the function 'interp1d'.  It takes
 in known data points and the points at which to interpolate values, and returns
@@ -161,6 +163,7 @@
 using bad_data removed both those points, so the entire range of x is linearly
 interpolated.
 
+
 --------------------------------------
 User-defined Interpolation Methods
 --------------------------------------
@@ -244,12 +247,25 @@
     Out []: array([ 5.7, 4.0, 5.7 ])
 
 
+-----------------------------
+Smoothing data
+-----------------------------
+Interpolate1d does not explicitly support smoothing of noisy data.  Other packages
+can be used to pre-filter the data before calling Interpolate1d on it, but
+a more convenient option is provided by the Spline class described below.
 
+Spline supports data smoothing with the keyword s (which defaults to 0, i.e. no
+smoothing).  If you want smoothing, you can set kind equal to an instance of the
+Spline class which has s != 0.  For example: ::
 
+    In []: import interpolate as I
+    In []: newy = I.interp1d(x, y, newx, kind=I.Spline(s=4.0) )
 
-================================================
+
+
+--------------------------------------------------------------------
 1D Interpolation with the Object Interface
-================================================
+--------------------------------------------------------------------
 
 interp1d is built as a wrapper around the class Interpolate1d.  If you want to
 interpolate multiple times from the same dataset, it can be more efficient
@@ -281,17 +297,17 @@
     In []: new_y2 == new_Y2
     Out []: True
     
-==================================================
+--------------------------------------------------------------------
 Sample Data Analysis Sessions Using Interpolate
-==================================================
+--------------------------------------------------------------------
 
 Below are several sample sessions or code pieces from various applications
 showing uses for interpolation and how it can be done using the
 interpolate module.
 
------------------------------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Estimating Function Statistics and Displaying Data
------------------------------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 In this session, the geologist
 has a data set of data indicating the temperature at various
@@ -322,9 +338,9 @@
     # are not necessarily uniform, but it is easy to uniformly sample the interpolated function
     In []: average_temp = average( I.interp1d(depth, temp, linspace(0,20,100), 'cubic', bad_data=[1000]) )
     
----------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Modeling from a small dataset
----------------------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 This computational biologist wants to model the growth rate of 
 cancer cells in tissue.  For several levels of blood glucose, he has measurements 
@@ -366,9 +382,9 @@
         def __call__(self, glucose_level):
             return self.CO2_to_growth( self.glucose_to_CO2( glucose_level ))
 
---------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Optimization
---------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 This engineer is developing a piece of hardware, and needs to find the optimal
 thickness for a thin film it contains.  Because performance (by some metric) is at a premium,
@@ -404,9 +420,9 @@
 More sophisticated optimization tools are also available from the scipy.optimize
 module.
 
-===================
+-----------------------------------------
 The Spline Class
-===================
+-----------------------------------------
 
 Interpolate1d, with the string arguments 'spline', 'cubic', 'quad', 'quintic', etc, is
 actually a wrapper around the Spline class, which contains fast and powerful Fortran
@@ -415,9 +431,9 @@
 
 This section describes the operation of the Spline class.
 
-----------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Intro to Splines
-----------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Splines are a class of functions which 
 #) are easy and quick to evaluate, 
@@ -440,9 +456,9 @@
 through all the data points but is smoother than it would be if it had to. k=3 is the most 
 common order of spline used in interpolation, and is often called a cubic spline.
 
--------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Basic Usage
--------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 At instantiation, the user passes in x, y, and possibly the spline order k (which defaults to 3).
 Calls are then made with the new_x array. ::
@@ -465,9 +481,9 @@
     In []: interp2.init_xy(x, y)
 
 
----------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Optional Arguments
----------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 At instantiation:
 
@@ -490,27 +506,27 @@
 to 0, so Spline usually returns the zeroth derivative of S, ie S.
 
 
------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 Special Methods
------------------
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-The following special methods are also available, which are not wrapped by Interpolate1d.
+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
+    returns the positions of the knots of the spline
 #. get_coeffs
-returns the coefficients of the 
+    returns the coefficients of the 
 #. 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)
+    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)
-returns the integral from a to b
+    returns the integral from a to b
 #. derivatives(x)
-returns all the derivatives of the spline at point x
+    returns all the derivatives of the spline at point x
 #. roots
-This only works for cubic splines.  But it returns the places where the spline
-is identically zero.
+    This only works for cubic splines.  But it returns the places where the spline
+    is identically zero.
 
 
 ================================================
@@ -530,29 +546,68 @@
 The Functional Interface
 ------------------------------------------
 
-The functional interface is virtually identical to that for interp1d: ::
+The functional interface is completely analogous to that for interp1d: ::
 
-    new_z = interp2d(x, y, z, new_x, new_y)
+    newz = interp2d(x, y, z, newx, newy, kind='linear', out=NaN)
 
-The range of interpolation is the rectangle given by the largest and
-smallest values in x and y.
-As in the case of 1D, string arguments can be passed
-as keywords to specify particular types of interpolation.  The keywords
-in this case are kind (for in-range interpolation) and out (for out-of-bounds).
+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
+a vector is returned.  The only differences from intper1d are
 
-By default, out of bounds returns NaN, and in-bounds returns a linear
-interpolation.
+#. The known data points are specified by 3 arrays (x, y and z) rather than 2 (x and y).
+    z is the dependent variable, while x and y are independent variables.
+#. Where to interpolate values is specified by two arrays, newx and newy, rather
+    than only one array.
+#. The extrapolation keywords "low" and "high" are replaced by the single keyword "out"
+    for out-of-bounds.
+#. Not all of the same keyword arguments are available for 1D and 2D.  The main ones like 
+    'linear', 'cubic' and 'spline', however, work in both cases, and we try to give analogous
+    methods the same name.  But some methods are particular to, or have only been written
+    for, one praticular dimensionality.
 
-new_x and new_y may be either arrays, lists or scalars.  If they are
-scalars or zero-dimensional arrays, new_z will be a scalar as well.  Otherwise
-a vector is returned.
+As in 1D, linear interpolation is used by default, while out of bounds returns NaN.
 
+Here is an example session: ::
+
+    In []: from interpolate import interp2d
+    In []: x, y = meshgrid( arange(10.), arange(10.) )
+    In []: z = x*y
+    In []: x, y, z = map(ravel, [x, y, z] )
+    # calling with scalars to get the value at a point
+    In []: interp2d(x, y, z, 5.5, 7.0)
+    Out []: 38.5
+    # plotting a cross-section of the surface
+    In []: newx = arange(0, 5, .2)
+    In []: newy = arange(0, 5, .2)
+    In []: plot( interp2d(x, y, z, newx, newy ))
+    
+.. image :: XY-crossection.png
+
+We can also use imshow to see the whole surface, though we must use
+reshape to turn z back into a 2D array. ::
+
+    In []: X, Y = meshgrid( arange(0,9,.2), arange(0,9,.2) )
+    In []: Z = interp2d(x, y, z, ravel(X), ravel(Y))
+    In []: len(Z)
+    Out []: 2025
+    # sqrt(2025) = 45
+    In []: imshow( reshape(Z, (45,45) ) )
+
+.. image :: XY-surface.png
+
 ------------------------------------------
 The Objective Interface
 ------------------------------------------
 
-The objective interface for 2D is similarly analogous to 1D.
+The objective interface for 2D is similarly analogous to 1D.  The class is
+instantiated with x, y, z, kind and out.  It is called with newx and newy. ::
 
+    In []: interp_obj = Interpolate2d(x, y, z)
+    In []: newz = interp_obj(newx, newy)
+
+
+
+
 ------------------------------------------
 The Spline2d Class
 ------------------------------------------
@@ -571,20 +626,18 @@
 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)
-
+    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
-
+    Integrate the interpolated function over the indicated rectangle
 #) get_residual
-
+    Same as Spline
 #) get_knots
-
+    Same as Spline
 #) get_coeffs
+    Same as Spline
     
 
 

Modified: branches/Interpolate1D/fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-04 19:48:18 UTC (rev 4599)
+++ branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-04 20:27:24 UTC (rev 4600)
@@ -120,22 +120,26 @@
             assigned the value of the nearest point that is within the
             interpolation range.
         """
+        # FIXME : this function calls self.get_grid, which is extremely inefficient.  However,
+        # I don't believe Fitpack really provides functionality to interpolate at scattered values.
         
-        if self._is_initialized is not True:
+        if self._is_initialized is False:
             raise Error, "x, y and z must be initialized before interpolating"
             
         # check input format
-        assert isinstance(x, np.ndarray) and isinstance(y, np.ndarray), \
+        assert ( isinstance(x, np.ndarray) and isinstance(y, np.ndarray) ), \
                     "newx and newy must both be numpy arrays"
         assert len(x) == len(y), "newx and newy must be of the same length"
         
         # sort only once for efficiency
-        sorted_x = sorted(x)
-        sorted_y = sorted(y)
+        sorted_x = x.copy()
+        sorted_x.sort()
+        sorted_y = y.copy()
+        sorted_y.sort()
         
         data_grid = self.get_grid(sorted_x, sorted_y)
         
-        # fixme : no list comprehension
+        # FIXME : no list comprehension
         z = np.array([ data_grid[np.searchsorted(sorted_x, x[i]), np.searchsorted(sorted_y,y[i])] \
                                     for i,xi in enumerate(x) ])
         

Modified: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py	2008-08-04 19:48:18 UTC (rev 4599)
+++ branches/Interpolate1D/interpolate2d.py	2008-08-04 20:27:24 UTC (rev 4600)
@@ -239,22 +239,23 @@
         
         # 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"
-        newy = atleast_1d_and_contiguous(newy, dtype=self._ydtype)
         assert newy.ndim == 1, "newy can be at most 1-dimensional"
         assert len(newx) == len(newy), "newx and newy must be the same length"
         
         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))
-        if sum(in_range_mask) > 0:
-            # hack to deal with behavior of vectorize on arrays of length 0
+        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])        
         if sum(~in_range_mask) > 0:
-            # hack to deal with behavior of vectorize on arrays of length 0
             result[~in_range_mask] = self.out(newx[~in_range_mask], newy[~in_range_mask])
         
+        # revert to scalar if applicable
         if input_is_scalar:
             result = result[0]
         



More information about the Scipy-svn mailing list