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

scipy-svn@scip... scipy-svn@scip...
Thu Aug 7 16:05:51 CDT 2008


Author: fcady
Date: 2008-08-07 16:05:37 -0500 (Thu, 07 Aug 2008)
New Revision: 4609

Added:
   branches/Interpolate1D/interpolateNd.py
Removed:
   branches/Interpolate1D/erics_notes.txt
   branches/Interpolate1D/fitpack_wrapper2d.py
   branches/Interpolate1D/ndimage_wrapper.py
   branches/Interpolate1D/regression_test.py
Modified:
   branches/Interpolate1D/TODO.txt
   branches/Interpolate1D/__init__.py
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/fitpack_wrapper.py
   branches/Interpolate1D/info.py
   branches/Interpolate1D/interpolate1d.py
   branches/Interpolate1D/interpolate2d.py
   branches/Interpolate1D/setup.py
   branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
   branches/Interpolate1D/tests/test_ndimage.py
Log:
msotly cleaning up and condensing the existing files.  The different projects and classes were scattered among a range of different files.  Also improvements to the documentation

Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/TODO.txt	2008-08-07 21:05:37 UTC (rev 4609)
@@ -4,37 +4,15 @@
 note with ideas.  More info is often contained in FIXMEs
 at appropriate places in the code.
 
+Note that I don't include much about docstrings, documentation, or
+comments here.  Improvement of those, esp as other stuff evolves,
+is a given.
 
+
 ************ API and/or MAJOR ISSUES ***********
 
-**handle smoothing 
-    There is a question of whether, if, and how much to allow
-    smoothing of the data.  If we smooth, we're not technically
-    interpolating, but users often want to smooth data.  
-
-    In fitpack_wrapper the smoothing parameter is s.  It now defaults
-    to 0.0 (exact interpolation).  Zero smoothing and moderate (s ~ 1)
-    are fast, but small non-zero s makes the algorithm VERY slow.
-    
-    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 
-    scipy.interpolate).  This isn't necessarily the best.  Other code
-    is used in scipy.ndimage and scipy.signal.  There is surely other
-    code out there too.  Figure out what is best and incorporate it.
-
-    Signal code is slower than FITPACK, and NDImage requires a
-    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).
     
 **allow y to be 2-dimensional?
     It is not decided whether this feature should be supported, but
@@ -53,99 +31,43 @@
     deleted for each column.
 
 **better handling of variable types
-    Currently everything is cast to a float64 if it is not already
-    a float32.  Is this the best way to do it?
-
-    There's also the question of Y being 2-dimensional and/or
-    incorporating strings for record arrays.  My instinct is to
-    have Interpolate1d only interpolate functions that are from
-    R1 -> R1.  That's how it is currently designed.  Other stuff
-    can be made as wrappers around Interpolate1d.
+    Currently only 1D interpolation deals with both float32
+    and float64 data.
     
-    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.
-    
-    
+**include spline pre-filtering in ND interpolation,
+    or understand why we don't need to do that.
 
-*********** DOCUMENTATION-TYPE AND NON-URGENT TASKS *******
+**README file that describes the architecture of the
+    package and also includes license information.
 
-**improve regression tests
-    desired for fitpack_wrapper and _interpolate_wrapper
-    as well as Interpolate1d.  What I have now is
-    really basic.
-
-
-**improve unit tests in tests directory
-
-
-**comment all files
-    There are comments there already, but they should be
-    made better.  Plus, features are changing, which requires
-    updating the documentation.
-
-
-**doc strings for interpolate1d and its members
-    There's docstrings there already, but they should be
-    made better.  In particular, it must be ensured that
-    they are of the proper format and include examples.
-
-    The doc strings for __init__.py, interpolate1d.py,
-    Interpolate1d, and interp1d are virtually identical
-    and very long; perhaps a master string can be stored
-    somewhere that they all reference.  This would make
-    updates of documentation easier.
-
-
 **figure out NumPy stuff with vectorize.
     In function Interpolate1d.__call__
     It would be nice to remove the hack I used.
     I believe vectorize is supposed to handle arrays of
     length 0, but it's not working on my computer.
-    
 
-**allow newx to be in non-sorted order
+**allow newx to be in non-sorted order for 1D
     This requires rethinking the partition of newx into
     low, high and mid
     
+**put all extension files into their own directory
+    
 ********* LONGER TERM ************
 
-**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.
+**allow for scattered ND data
+    A Delaunay triangulation may be a good way to do this. I
+    would really like to code it up, but that may be a bit preliminary.
     
-    There is an interesting problem here.  Most of the extensions
-    I have assume a regular grid.  First off, this is often not general.
-    Secondly, if I DO use a regular grid, how do I deal with bad
-    data?  The best way is probably a pre-processing where you
-    interpolate values for the bad points (linear would be a nice simple
-    way to do it at first, just to get it working)
-    
-    We should probably use something other than FITPACK for this.
-    But firstly it's at most 2D.  Much worse, it doesn't evaluate at
-    a set of points; it evaluates over a grid, which requires inputs
-    being in sorted order (both in x and y coordinates).  This makes
-    input inconvenient and the code runs a lot slower than ndimage.
-    ND image has 2 main downsides :1) it requires x and y be uniform
-    spacing (since its really interpolating entries of an array, rather
-    than values of a function on Rn), and 2) the C code is TERRIBLY
-    documented/commented.  But that can be fixed.
-    
-    So we may have two separate classes: Interpolate1d which is
-    based on FITPACK (or something else, depending on what
-    happens with the smoothing parameter), and InterpolateNd
-    which is based on ndimage.
-    
-    Another option is to have two classes: one for uniformly spaced data and
-    and one for scattered data.  Regularly spaced would use NDImage, and
-    scattered would start as an inefficient wrapper around FITPACK.  But longer
-    term the scattered data could be done with Delaunay triangulation,
-    perhaps something that would implicitly calculate the convex hull
-    of the points and interpolate within it.
+**rethink extension packages and capabilities
+    Under-the-hood machinery currently comes from _interpolate.cpp
+    (used in enthought.interpolate), FITPACK (Fortran, used in 
+    scipy.interpolate) and ndimage.  This isn't necessarily the best.  
+    Other code is used in scipy.signal.  There is surely other
+    code out there too.  Figure out what is best and incorporate it.
 
+    Currently, 1D and 2D both mostly wrap FITPACK through the Spline
+    and Spline2d classes, and ND wraps ndimage.
+
 **high-level road map
     when the module is more established, there should be a page on
     the wiki which describes the big-picture of the module; what

Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/__init__.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -2,11 +2,16 @@
 # information about included functions
 from info import __doc__
 
+# primary usage for 1, 2, and N-dimensional interpolation
+from interpolate1d import Interpolate1d, interp1d
+from interpolate2d import Interpolate2d, interp2d
+from interpolateNd import InterpolateNd, interpNd
+
 # basic interpolation routines
+# wrapped by interpolate*.py files
 from interpolate_wrapper import linear, logarithmic, block, block_average_above
 
 # support for spline interpolation
+# wrapped by interpolate*.py files
 from fitpack_wrapper import Spline
 
-# wrapping for all supported interpolation types.
-from interpolate1d import Interpolate1d, interp1
\ No newline at end of file

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-07 21:05:37 UTC (rev 4609)
@@ -683,16 +683,15 @@
 #) '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'
+    where the strings indicate extrapolation methods.
 
 The user cannot pass in specially-tailored interpolation methods.
 
-There 
+There is also an objective interface that interpNd wraps around.  The
+class InterpolateNd is instantiated with all arguments other than coordinates,
+and called with coordinates.
 
-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.  
+Finally, removal of bad data points is not supported for ND interpolation.
 
 
 

Deleted: branches/Interpolate1D/erics_notes.txt
===================================================================
--- branches/Interpolate1D/erics_notes.txt	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/erics_notes.txt	2008-08-07 21:05:37 UTC (rev 4609)
@@ -1,92 +0,0 @@
-*. _remove_bad_data should be left up to the interpolation method if it knows what to do.
-   otherwise, it is handled by this top level class.
-   And, we definitely don't want the list comprehension in the remove_bad_data class.
-
-*. Function signatures:
-    
-    def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
-                        kindkw={}, lowkw={}, highkw={}, \
-                        remove_bad_data = False, bad_data=[], interp_axis = 0):
-    
-     o. It is dangerous to initialize containers as default arguments because they are
-        effectively a singleton for the function.  We can discuss this if you don't know
-        what I am talking about.
-     Xo. kindkw, lowkw, and highkw aren't really necessary I don't think.  They should go.
-     Xo. [minor] I would just use NaN instead of np.NaN.  						
-     Xo. the trailing \ for new lines are not necessary for functions since line continuation
-        is implicit with the open/close parentheses.
-	 Xo. Do you need both remove_bad_data and bad_data?  If bad_data is None, then you 
-	    don't remove bad_data...
-	 Xo. I think I would change interp_axis to just be axis.  This is consistent with many
-	    other functions in numpy.
-	 Xo. The choice of whether axis=0 or axis=-1 by default is a reasonable question.
-	    fft defaults to axis=-1.  This is also the faster axis to operate across in the
-	    standard case.  It is, however, the opposite of how some people think about 
-	    things (columns vs. rows).  Talk to Travis O. for his take.  Mine is to use axis=-1.
-	
-	I think all of this might simplify the interface to the following:
-	
-    def interp1d(x, y, new_x, kind='linear', low=NaN, high=NaN, bad_data=None, axis=-1):
-
-X*. Follow scipy/FORMAT_GUIDLINES.txt in the main scipy directory.
-
-For example:
-
-	test_callFormat -> test_call_format
-
-	Here are the rules:
-
-	Follow the standard Python formatting rules when writing code for SciPy.  Guido
-	describes them "here":http://www.python.org/doc/essays/styleguide.html.  A few
-	reminders follow:
-
-	   Xo Use 4 spaces for indentation levels.  Do not use tabs as they can result
-	     in indentation confusion.  Most editors have a feature that will insert 4
-	     spaces when the tab key is hit.  Also, many editors will automatically
-	     search/replace leading tabs with 4 spaces.
-
-	   o Only 80 characters on a line.
-
-	   Xo use all lowercase function names with underscore separated words:
-
-	        def set_some_value()
-
-	     instead of:
-
-	        def setSomeValue()
-
-
-	   Xo use CamelCase class names:
-
-	        def BaseClass()
-
-	     instead of:
-
-	        def base_class()
-X*. I am glad to see your docstrings.
-
-X*. Move tests into a seprate test directory.
-
-X*. make_array_safe
-   I would prefer 'make_array_safe' named atleast_1d_and_contiguous().  This is more specific
-   and it is immediately clear to other developers what the check does.  If you add other checks,
-   then perhaps come up with a more generic name, but this explicit names, when possible, help
-   readability.
-
-   Also, numpy has an ascontiguousrarray() function that would simplify the code a line or so.
-
-   Also, this function lives in multiple places, interpolate1d and interpolate_wrapper.
-
-X*. isinstance(value, str) should be isinstance(value, basestring) so that we handle
-   both strings and unicode correctly.
-
-X*. If the input to interp1d is a scalar, the return value should be a scalar.
-   [add test to handle this.]
-   This fails in the following:
-
-    In [5]: x = arange(10)
-	In [6]: y = arange(10)*2.
-	In [7]: interpolate1d.interp1d(x,y,3.2)
-	Out[7]: array([ 6.4])
-
-  
\ No newline at end of file

Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/fitpack_wrapper.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -7,6 +7,10 @@
     Spline is primarily meant to be called by Interpolate1d
     or interp1d, but is a stand-alone class in its own right
     that is not available through these interfaces.
+    
+    Spline2d is primarily meant to be called by Interpolate2d
+    or interp2d, but is a stand-alone class in its own right
+    that is not available through these interfaces.
 
     The code has been modified from an older version of
     scipy.interpolate, where it was directly called by the
@@ -17,6 +21,7 @@
 """
 
 import numpy as np
+import warnings
 
 import _dfitpack # extension module containing FITPACK subroutines in Fortran
 
@@ -167,4 +172,206 @@
             assert ier==0,`ier`
             return z[:m]
         raise NotImplementedError,\
-              'finding roots unsupported for non-cubic splines'
\ No newline at end of file
+              'finding roots unsupported for non-cubic splines'
+              
+
+############################
+## BELOW THIS POINT IS CODE FOR 2D INTERPOLATION
+############################
+
+_surfit_messages = {1:"""
+        The required storage space exceeds the available storage space: nxest
+        or nyest too small, or s too small.
+        The weighted least-squares spline corresponds to the current set of
+        knots.""",
+                            2:"""
+        A theoretically impossible result was found during the iteration
+        process for finding a smoothing spline with fp = s: s too small or
+        badly chosen eps.
+        Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
+                            3:"""
+        the maximal number of iterations maxit (set to 20 by the program)
+        allowed for finding a smoothing spline with fp=s has been reached:
+        s too small.
+        Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
+                            4:"""
+        No more knots can be added because the number of b-spline coefficients
+        (nx-kx-1)*(ny-ky-1) already exceeds the number of data points m:
+        either s or m too small.
+        The weighted least-squares spline corresponds to the current set of
+        knots.""",
+                            5:"""
+        No more knots can be added because the additional knot would (quasi)
+        coincide with an old one: s too small or too large a weight to an
+        inaccurate data point.
+        The weighted least-squares spline corresponds to the current set of
+        knots.""",
+                            10:"""
+        Error on entry, no approximation returned. The following conditions
+        must hold:
+        xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
+        If iopt==-1, then
+          xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
+          yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""",
+                            -3:"""
+        The coefficients of the spline returned have been computed as the
+        minimal norm least-squares solution of a (numerically) rank deficient
+        system (deficiency=%i). If deficiency is large, the results may be
+        inaccurate. Deficiency may strongly depend on the value of eps."""
+                    }
+
+class Spline2d(object):
+    """ Bivariate spline s(x,y) of degrees kx and ky on the rectangle
+        [xb,xe] x [yb, ye] calculated from a given set of data points
+        (x,y,z).
+
+        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):
+        """
+            Input:
+              x,y,z  - 1-d sequences of data points (order is not
+                       important)
+            Optional input:
+              w          - positive 1-d sequence of weights
+              bbox       - 4-sequence specifying the boundary of
+                           the rectangular approximation domain.
+                           By default, bbox=[min(x,tx),max(x,tx),
+                                             min(y,ty),max(y,ty)]
+              kx,ky=3,3  - degrees of the bivariate spline.
+              s          - positive smoothing factor defined for
+                           estimation condition:
+                             sum((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0) <= s
+                           Default s=len(w) which should be a good value
+                           if 1/w[i] is an estimate of the standard
+                           deviation of z[i].
+              eps        - a threshold for determining the effective rank
+                           of an over-determined linear system of
+                           equations. 0 < eps < 1, default is 1e-16.
+        """
+        
+        self._w = w
+        self._bbox = bbox
+        self._kx = kx
+        self._ky = kx
+        self._s = s
+        self._eps = eps
+        
+        if x is not None and y is not None and z is not None:
+            self.init_xyz(x, y, z)
+            self._is_initialized = True
+        else:
+            self._is_initialized = False
+        
+    def init_xyz(self, x, y, z):
+        xb,xe,yb,ye = self._bbox
+        nx,tx,ny,ty,c,fp,wrk1,ier = _dfitpack.surfit_smth(x,y,z,
+                                                         self._w,
+                                                         xb, xe, yb, ye,
+                                                         self._kx, self._ky,
+                                                         s=self._s,
+                                                         eps=self._eps, lwrk2=1)
+        if ier in [0,-1,-2]: # normal return
+            pass
+        else:
+            message = _surfit_messages.get(ier,'ier=%s' % (ier))
+            warnings.warn(message)
+
+        self.fp = fp
+        self.tck = tx[:nx],ty[:ny],c[:(nx-self._kx-1)*(ny-self._ky-1)]
+        self.degrees = self._kx,self._ky
+        
+        self._is_initialized = True
+        
+    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.
+        """
+        # 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 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) ), \
+                    "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 = 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
+        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):
+        """ 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"
+        
+        # check input format
+        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"
+        
+        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)
+        """
+        return self.fp
+    def get_knots(self):
+        """ Return a tuple (tx,ty) where tx,ty contain knots positions
+            of the spline with respect to x-, y-variable, respectively.
+            The position of interior and additional knots are given as
+              t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.
+        """
+        return self.tck[:2]
+    def get_coeffs(self):
+        """ Return spline coefficients."""
+        return self.tck[2]
+    
+    
+    def integral(self, xa, xb, ya, yb):
+        """
+            Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
+            
+            Parameters
+            ----------
+            xa, xb : float
+                The end-points of the x integration interval.
+            ya, yb : float
+                The end-points of the y integration interval.
+            
+            Returns
+            -------
+            integ : float
+                The value of the resulting integral.
+            
+        """
+        tx,ty,c = self.tck[:3]
+        kx,ky = self.degrees
+        return _dfitpack.dblint(tx,ty,c,kx,ky,xa,xb,ya,yb)
\ No newline at end of file

Deleted: branches/Interpolate1D/fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/fitpack_wrapper2d.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -1,206 +0,0 @@
-import warnings
-from numpy import zeros, concatenate, alltrue, ravel, all, diff
-import numpy as np
-
-import _dfitpack
-
-_surfit_messages = {1:"""
-        The required storage space exceeds the available storage space: nxest
-        or nyest too small, or s too small.
-        The weighted least-squares spline corresponds to the current set of
-        knots.""",
-                            2:"""
-        A theoretically impossible result was found during the iteration
-        process for finding a smoothing spline with fp = s: s too small or
-        badly chosen eps.
-        Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
-                            3:"""
-        the maximal number of iterations maxit (set to 20 by the program)
-        allowed for finding a smoothing spline with fp=s has been reached:
-        s too small.
-        Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
-                            4:"""
-        No more knots can be added because the number of b-spline coefficients
-        (nx-kx-1)*(ny-ky-1) already exceeds the number of data points m:
-        either s or m too small.
-        The weighted least-squares spline corresponds to the current set of
-        knots.""",
-                            5:"""
-        No more knots can be added because the additional knot would (quasi)
-        coincide with an old one: s too small or too large a weight to an
-        inaccurate data point.
-        The weighted least-squares spline corresponds to the current set of
-        knots.""",
-                            10:"""
-        Error on entry, no approximation returned. The following conditions
-        must hold:
-        xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
-        If iopt==-1, then
-          xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
-          yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""",
-                            -3:"""
-        The coefficients of the spline returned have been computed as the
-        minimal norm least-squares solution of a (numerically) rank deficient
-        system (deficiency=%i). If deficiency is large, the results may be
-        inaccurate. Deficiency may strongly depend on the value of eps."""
-                    }
-
-class Spline2d(object):
-    """ Bivariate spline s(x,y) of degrees kx and ky on the rectangle
-        [xb,xe] x [yb, ye] calculated from a given set of data points
-        (x,y,z).
-
-        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):
-        """
-            Input:
-              x,y,z  - 1-d sequences of data points (order is not
-                       important)
-            Optional input:
-              w          - positive 1-d sequence of weights
-              bbox       - 4-sequence specifying the boundary of
-                           the rectangular approximation domain.
-                           By default, bbox=[min(x,tx),max(x,tx),
-                                             min(y,ty),max(y,ty)]
-              kx,ky=3,3  - degrees of the bivariate spline.
-              s          - positive smoothing factor defined for
-                           estimation condition:
-                             sum((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0) <= s
-                           Default s=len(w) which should be a good value
-                           if 1/w[i] is an estimate of the standard
-                           deviation of z[i].
-              eps        - a threshold for determining the effective rank
-                           of an over-determined linear system of
-                           equations. 0 < eps < 1, default is 1e-16.
-        """
-        
-        self._w = w
-        self._bbox = bbox
-        self._kx = kx
-        self._ky = kx
-        self._s = s
-        self._eps = eps
-        
-        if x is not None and y is not None and z is not None:
-            self.init_xyz(x, y, z)
-            self._is_initialized = True
-        else:
-            self._is_initialized = False
-        
-    def init_xyz(self, x, y, z):
-        xb,xe,yb,ye = self._bbox
-        nx,tx,ny,ty,c,fp,wrk1,ier = _dfitpack.surfit_smth(x,y,z,
-                                                         self._w,
-                                                         xb, xe, yb, ye,
-                                                         self._kx, self._ky,
-                                                         s=self._s,
-                                                         eps=self._eps, lwrk2=1)
-        if ier in [0,-1,-2]: # normal return
-            pass
-        else:
-            message = _surfit_messages.get(ier,'ier=%s' % (ier))
-            warnings.warn(message)
-
-        self.fp = fp
-        self.tck = tx[:nx],ty[:ny],c[:(nx-self._kx-1)*(ny-self._ky-1)]
-        self.degrees = self._kx,self._ky
-        
-        self._is_initialized = True
-        
-    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.
-        """
-        # 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 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) ), \
-                    "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 = 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
-        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):
-        """ 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"
-        
-        # check input format
-        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"
-        
-        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)
-        """
-        return self.fp
-    def get_knots(self):
-        """ Return a tuple (tx,ty) where tx,ty contain knots positions
-            of the spline with respect to x-, y-variable, respectively.
-            The position of interior and additional knots are given as
-              t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.
-        """
-        return self.tck[:2]
-    def get_coeffs(self):
-        """ Return spline coefficients."""
-        return self.tck[2]
-    
-    
-    def integral(self, xa, xb, ya, yb):
-        """
-            Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
-            
-            Parameters
-            ----------
-            xa, xb : float
-                The end-points of the x integration interval.
-            ya, yb : float
-                The end-points of the y integration interval.
-            
-            Returns
-            -------
-            integ : float
-                The value of the resulting integral.
-            
-        """
-        tx,ty,c = self.tck[:3]
-        kx,ky = self.degrees
-        return _dfitpack.dblint(tx,ty,c,kx,ky,xa,xb,ya,yb)
-        
-# RectBivariateSpline in scipy.interpolate is for a rectangular grid and presumably must faster.  There are 3 levels of niceness: scattered
-#       data, irregular grid, and regular grids.  Spline2d is for the first level and thus slow.  ndimage is for the 3rd level and thus fast.
-#       I vote to no explicitly treat the 3rd level, but RecBivariateSpline does that if we want to implement it in the future.
\ No newline at end of file

Modified: branches/Interpolate1D/info.py
===================================================================
--- branches/Interpolate1D/info.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/info.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -3,19 +3,30 @@
 
 __doc__ = \
 """
-    Interpolation of 1D data
-
     This module provides several functions and classes for interpolation
-    and extrapolation of 1D (in both input and output) real-valued.  The
-    primary function provided is:
+    and extrapolation of real-valued data.  The primary interface is through
+    the functions:
 
-        interp1d(x, y, new_x) :  from data points (x[i], y[i]), interpolates
+        interp1d(x, y, newx) :  from data points (x[i], y[i]), interpolates
                                         values for points in new_x and
-                                        returns them as an array.  x and new_x
+                                        returns them as an array.  x and newx
                                         must both be in sorted order.
+        
+        interp2d(x, y, z, newx, newy): from data points (x[i], y[i], z[i]), interpolates
+                                        values for points (newx[i], newy[i]) and
+                                        returns them as an array.
+                                        
+        interpNd(data, coordinates): data contains known values (which are
+                                        assume to be uniformly spaced), and coordinates
+                                        indicates where to interpolate new values.
+                                        
+    Each function defaults to linear interpolation for in-range points, and
+    returning a NaN when points are out of range.  However, by specifying keywords
+    a variety of techniques can be chosen or, for the first wo functions, even explicitly
+    written by the user.
+    
+    The following callable classes are also provided:
 
-    Classes provided include:
-
         Interpolate1d  :   an object for interpolation of
                                 various kinds.  interp1d is a wrapper
                                 around this class.
@@ -24,21 +35,21 @@
                                 wraps this class if spline interpolation
                                 is used.  However, not all functionality
                                 of Spline is available through Interpolate1d.
+                                
+        Interpolate2d  :
         
-    Functions provided include:
+        Spline 2d  :  
+        
+        InterpolateNd  :
+        
+    These functions and classes constitute the primary api.  However, several
+    additional functions are also provided (the primary api in many cases calls
+    these functions):
 
         linear : linear interpolation
         logarithmic :  logarithmic interpolation
         block : block interpolation
         block_average_above : block average above interpolation
-        
-    The dependency/interface architecture is as follows:
-        interpolate1d.py is viewed by the user (through interp1d and Interpolate1d)
-            It depends on fitpack_wrapper.py and interpolate_wrapper.py
-        fitpack_wrapper is viewed by the user (through Spline)
-            It depends on dfitpack.pyd, a Fortran extension module
-        interpolate_wrapper is viewed by he user (through functions linear, etc)
-            It depends on _interpolate.pyd, a C extension module.
 
 """
         

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/interpolate1d.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -130,7 +130,7 @@
         ---------
         
             >>> import numpy
-            >>> from Interpolate1D import interp1d
+            >>> from interpolate import interp1d
             >>> x = range(5)        # note list is permitted
             >>> y = numpy.arange(5.)
             >>> new_x = [.2, 2.3, 5.6]

Modified: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/interpolate2d.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -1,7 +1,7 @@
 
 from numpy import NaN, array
 import numpy as np
-from fitpack_wrapper2d import Spline2d
+from fitpack_wrapper import Spline2d
 
 def atleast_1d_and_contiguous(ary, dtype = np.float64):
     # FIXME : don't have in 2 places

Copied: branches/Interpolate1D/interpolateNd.py (from rev 4608, branches/Interpolate1D/ndimage_wrapper.py)
===================================================================
--- branches/Interpolate1D/ndimage_wrapper.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/interpolateNd.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -0,0 +1,298 @@
+""" ND Interpolation wrapping using code from NDImage"""
+
+from numpy import array, arange, NaN
+import numpy as np
+import _nd_image
+
+def interpNd(data, coordinates, starting_coords=None, spacings=None, kind='linear',out=NaN):
+    """ A function for interpolation of 1D, real-valued data.
+        
+        Parameters
+        -----------
+            
+            data -- NumPy array (N-dimensional) or list of lists
+                indicates the known values of the function.
+            
+            coordinates -- array or list
+                To interpolate at a set of L points, array must be NxL, where
+                each column denotes a point.  If only one point is desired, its
+                coordinates may be entered as either a list or an array.
+                
+        Optional Arguments
+        -------------------
+        
+            starting_coords -- array or list
+                indicates the point in space
+                whose value is given by data[0, ..., 0].
+                Defaults to being all zeros.
+                
+            spacings -- array or list
+                jth component gives spacing
+                of points along the jth axis.  Defaults
+                to being all ones.
+        
+            kind -- A string or integer
+                Indicates what interpolation method to perform on
+                points within the region of interpolation
+                
+                0, 'block' -- block interpolation based on interval midpoints
+                1, 'linear' -- linear interpolation
+                2, 'quadratic' -- spline order 2 interpolation
+                3, 'cubic' -- cubic spline interpolation
+                4, 'quartic' -- 4th order spline interpolation
+                5, 'quintic' -- 5th order spine interpolation
+                
+            out -- string or NaN
+                Indicates how to extrapolate values at points outside 
+                the region of interpolation.
+            
+                NaN -- return NaN for all points out of range
+                'nearest' -- return value at nearest valid point
+                'constant' -- returns 0
+                'wrap' -- points over one boundary wrap around to the other side
+                'reflect' -- out-of-bounds points are reflected into the valid region
+                
+        Example
+        --------
+        
+            >>> import numpy as np
+            >>> from interpolate import interpNd
+            >>> boring_data = np.ones((5,5,5))
+            >>> nd.interpNd(boring_data, np.array([[2.3], [1.0], [3.9]]))
+            1.0
+    """
+    return InterpolateNd(data = data,
+                                    starting_coords = starting_coords,
+                                    spacings = spacings,
+                                    kind = kind,
+                                    out = out
+                                    )(coordinates)
+
+class InterpolateNd:
+    """ A callable class for interpolation of 1D, real-valued data.
+        
+        Parameters
+        -----------
+            
+            data -- NumPy array (N-dimensional) or list of lists
+                indicates the known values of the function.
+                
+        Optional Arguments
+        -------------------
+        
+            starting_coords -- array or list
+                indicates the point in space
+                whose value is given by data[0, ..., 0].
+                Defaults to being all zeros.
+                
+            spacings -- array or list
+                jth component gives spacing
+                of points along the jth axis.  Defaults
+                to being all ones.
+        
+            kind -- A string or integer
+                Indicates what interpolation method to perform on
+                points within the region of interpolation
+                
+                0, 'block' -- block interpolation based on interval midpoints
+                1, 'linear' -- linear interpolation
+                2, 'quadratic' -- spline order 2 interpolation
+                3, 'cubic' -- cubic spline interpolation
+                4, 'quartic' -- 4th order spline interpolation
+                5, 'quintic' -- 5th order spine interpolation
+                
+            out -- string or NaN
+                Indicates how to extrapolate values at points outside 
+                the region of interpolation.
+            
+                NaN -- return NaN for all points out of range
+                'nearest' -- return value at nearest valid point
+                'constant' -- returns 0
+                'wrap' -- points over one boundary wrap around to the other side
+                'reflect' -- out-of-bounds points are reflected into the valid region
+                
+        Example
+        --------
+        
+            >>> import numpy as np
+            >>> from interpolate import InterpolateNd
+            >>> boring_data = np.ones((5,5,5))
+            >>> nd.InterpolateNd(boring_data)( np.array([[2.3], [1.0], [3.9]]) )
+            1.0
+    """
+    def __init__(self, data, starting_coords =None, spacings = None, 
+                        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
+        
+        # checking format of input
+        data = array(data)
+        
+        # for proper processing later, starting_coords and spacings must be of shape (data.ndim, 1)
+        if starting_coords == None:
+            starting_coords = np.zeros(( data.ndim, 1 ))
+        else:
+            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 = np.reshape(starting_coords, (data.ndim, 1))
+        if spacings == None:
+            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 = 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
+        self._shape = np.shape(data)
+        self._spacings = spacings
+        self._min_coords = starting_coords
+        self._max_coords = self._min_coords + self._shape*self._spacings
+        self.out = out
+        
+    def __call__(self, coordinates):
+        """ coordinates is an n x L array, where n is the dimensionality of the data
+            and L is number of points.  That is, each column of coordinates
+            indicates a point at which to interpolate.
+        """
+        
+        # format checking
+        coordinates = array(coordinates)
+        if coordinates.ndim == 1: # passed in a single point
+            coordinates = np.reshape(coordinates, ( self.ndim, 1))
+        assert coordinates.ndim == 2, "Coordinates must be 1 or 2 dimensional"
+        n, num_points = coordinates.shape
+        assert n == self.ndim, "The first dimension of the input \
+                must be as long as the dimensionality of the space"
+        
+        # converting from points in ND space to array indices
+        indices = (coordinates - self._min_coords)/self._spacings
+        
+        if self.out in ['nearest', 'wrap', 'reflect', 'mirror', 'constant']:
+            # out of bounds can be performed by _interpolate_array_entry
+            result = self._interpolate_array_entry(self._data_array, indices, self.order, out = self.out)
+        else:
+            # need to return NaN when entry is out of bounds
+            in_bounds_mask = self._index_in_bounds(indices)
+            in_bounds = indices[:, in_bounds_mask]
+            out_bounds = indices[:, ~in_bounds_mask]
+            
+            result = np.zeros(num_points)
+            result[in_bounds_mask] = \
+                self._interpolate_array_entry(self._data_array, indices[:,in_bounds_mask], self.order)
+            result[~in_bounds_mask] = NaN
+            
+        return result
+        
+    
+    def _interpolate_array_entry(self, data_array, indices, order, out='nearest'):
+        """ indices is nxL matrix, where n is data_array.ndim
+            returns array of length L giving interpolated entries.
+        """
+        
+        extrap_code_register = { 'nearest':0,
+                                        'wrap': 1,
+                                        'reflect':2,
+                                        'mirror':3,
+                                        'constant':4,
+                                        }
+        
+        n, L = np.shape(indices)
+        
+        output = np.zeros( L , dtype=np.float64 ) # place to store the data
+
+        # geometric transform takes data_array, interpolates its values at indices, and
+        # stores those values in output.  Other parameters give details of interpolation method.
+        _nd_image.geometric_transform(data_array, None, indices, None, None, \
+               output, order, extrap_code_register[out], 0.0, None, None)
+               
+        return output
+        
+    def _index_in_bounds(self, indices):
+        """ return an array of bools saying which
+            points are in interpolation bounds
+        """
+        shape_as_column_vec = np.reshape(self._shape, (self.ndim, 1))
+        
+        # entry is 1 if that coordinate of a point is in its bounds
+        index_in_bounds = (0 <= indices) & \
+                                    (indices <= shape_as_column_vec)
+        
+        # for each point, number of coordinates that are in bounds
+        num_indices_in_bounds = np.sum(index_in_bounds, axis=0)
+        
+        # True if each coordinate for the point is in bounds
+        return num_indices_in_bounds == self.ndim
+        
+    def _coord_in_bounds(self, coordinates):
+        """ return an array of bools saying which
+            points are in interpolation bounds
+        """
+        # entry is 1 if that coordinate of a point is in its bounds
+        coord_in_bounds = (self._min_coords <= coordinates) & \
+                                    (coordinates <= self._max_coords)
+        
+        # for each point, number of coordinates that are in bounds
+        num_coords_in_bounds = np.sum(coord_in_bounds, axis=0)
+        
+        # True if each coordinate for the point is in bounds
+        return num_coords_in_bounds == self.ndim
+        
+    
+        
+    
+        
+        
+        
+    
+    
\ No newline at end of file

Deleted: branches/Interpolate1D/ndimage_wrapper.py
===================================================================
--- branches/Interpolate1D/ndimage_wrapper.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/ndimage_wrapper.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -1,191 +0,0 @@
-""" ND Interpolation wrapping using code from NDImage"""
-
-from numpy import array, arange, NaN
-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, 
-                        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
-        
-        # checking format of input
-        data = array(data)
-        
-        # for proper processing later, starting_coords and spacings must be of shape (data.ndim, 1)
-        if starting_coords == None:
-            starting_coords = np.zeros(( data.ndim, 1 ))
-        else:
-            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 = np.reshape(starting_coords, (data.ndim, 1))
-        if spacings == None:
-            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 = 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
-        self._shape = np.shape(data)
-        self._spacings = spacings
-        self._min_coords = starting_coords
-        self._max_coords = self._min_coords + self._shape*self._spacings
-        self.out = out
-        
-    def __call__(self, coordinates):
-        """ coordinates is an n x L array, where n is the dimensionality of the data
-            and L is number of points.  That is, each column of coordinates
-            indicates a point at which to interpolate.
-        """
-        
-        # format checking
-        coordinates = array(coordinates)
-        if coordinates.ndim == 1: # passed in a single point
-            coordinates = np.reshape(coordinates, ( self.ndim, 1))
-        assert coordinates.ndim == 2, "Coordinates must be 1 or 2 dimensional"
-        n, num_points = coordinates.shape
-        assert n == self.ndim, "The first dimension of the input \
-                must be as long as the dimensionality of the space"
-        
-        # converting from points in ND space to array indices
-        indices = (coordinates - self._min_coords)/self._spacings
-        
-        if self.out in ['nearest', 'wrap', 'reflect', 'mirror', 'constant']:
-            # out of bounds can be performed by _interpolate_array_entry
-            result = self._interpolate_array_entry(self._data_array, indices, self.order, out = self.out)
-        else:
-            # need to return NaN when entry is out of bounds
-            in_bounds_mask = self._index_in_bounds(indices)
-            in_bounds = indices[:, in_bounds_mask]
-            out_bounds = indices[:, ~in_bounds_mask]
-            
-            result = np.zeros(num_points)
-            result[in_bounds_mask] = \
-                self._interpolate_array_entry(self._data_array, indices[:,in_bounds_mask], self.order)
-            result[~in_bounds_mask] = NaN
-            
-        return result
-        
-    
-    def _interpolate_array_entry(self, data_array, indices, order, out='nearest'):
-        """ indices is nxL matrix, where n is data_array.ndim
-            returns array of length L giving interpolated entries.
-        """
-        
-        extrap_code_register = { 'nearest':0,
-                                        'wrap': 1,
-                                        'reflect':2,
-                                        'mirror':3,
-                                        'constant':4,
-                                        }
-        
-        n, L = np.shape(indices)
-        
-        output = np.zeros( L , dtype=np.float64 ) # place to store the data
-
-        # geometric transform takes data_array, interpolates its values at indices, and
-        # stores those values in output.  Other parameters give details of interpolation method.
-        _nd_image.geometric_transform(data_array, None, indices, None, None, \
-               output, order, extrap_code_register[out], 0.0, None, None)
-               
-        return output
-        
-    def _index_in_bounds(self, indices):
-        """ return an array of bools saying which
-            points are in interpolation bounds
-        """
-        shape_as_column_vec = np.reshape(self._shape, (self.ndim, 1))
-        
-        # entry is 1 if that coordinate of a point is in its bounds
-        index_in_bounds = (0 <= indices) & \
-                                    (indices <= shape_as_column_vec)
-        
-        # for each point, number of coordinates that are in bounds
-        num_indices_in_bounds = np.sum(index_in_bounds, axis=0)
-        
-        # True if each coordinate for the point is in bounds
-        return num_indices_in_bounds == self.ndim
-        
-    def _coord_in_bounds(self, coordinates):
-        """ return an array of bools saying which
-            points are in interpolation bounds
-        """
-        # entry is 1 if that coordinate of a point is in its bounds
-        coord_in_bounds = (self._min_coords <= coordinates) & \
-                                    (coordinates <= self._max_coords)
-        
-        # for each point, number of coordinates that are in bounds
-        num_coords_in_bounds = np.sum(coord_in_bounds, axis=0)
-        
-        # True if each coordinate for the point is in bounds
-        return num_coords_in_bounds == self.ndim
-        
-    
-        
-    
-        
-        
-        
-    
-    
\ No newline at end of file

Deleted: branches/Interpolate1D/regression_test.py
===================================================================
--- branches/Interpolate1D/regression_test.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/regression_test.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -1,33 +0,0 @@
-""" 
-    regression test:
-
-    This script runs a simple regression test on the functionality of
-    the interpolation module.  Currently, when run, it times each
-    unit test in interpolate1d.py and stores those times in a dict
-    of dicts; outer keys are time test was performed, and inner
-    keys are names of tests run.
-
-"""
-
-import shelve, time
-#from tests/test_interpolate1d import Test
-
-# name of log file to which all data is stored.
-filename = 'regression_test.dbm'
-
-#log_total = shelve.open(filename)
-current_time = str(time.localtime()[0:5]) # specified up to the minute
-
-# run all tests in interpolate1d's test class
-#test_list = [name for name in dir(Test) if name.find('test_') == 0]
-log_now = {}
-
-# record time taken for each test
-#for test_name in test_list:
-  #  t1 = time.clock()
-   # eval('Test.%s' % test_name)
-    #t2 = time.clock()
-    #log_now[test_name] = t2-t1
-
-#log_total[current_time] = log_now
-#log_total.close()

Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/setup.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -40,7 +40,6 @@
                         include_dirs=['ndimage']+[get_include()],
                         )
                         
-    # FIXME : add documentation files
     config.add_data_dir('docs')
 
     return config

Modified: branches/Interpolate1D/tests/test_fitpack_wrapper2d.py
===================================================================
--- branches/Interpolate1D/tests/test_fitpack_wrapper2d.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/tests/test_fitpack_wrapper2d.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -10,7 +10,7 @@
 import time
 from numpy import arange, allclose, ones, meshgrid, ravel, array
 import numpy as np
-from fitpack_wrapper2d import Spline2d
+from fitpack_wrapper import Spline2d
 
 class Test(unittest.TestCase):
     

Modified: branches/Interpolate1D/tests/test_ndimage.py
===================================================================
--- branches/Interpolate1D/tests/test_ndimage.py	2008-08-07 19:31:44 UTC (rev 4608)
+++ branches/Interpolate1D/tests/test_ndimage.py	2008-08-07 21:05:37 UTC (rev 4609)
@@ -9,13 +9,20 @@
 import time
 from numpy import arange, allclose, ones, array
 import numpy as np
-import ndimage_wrapper as nd
+import interpolateNd as nd
 
 class Test (unittest.TestCase):
     
     def assertAllclose(self, x, y):
         self.assert_(np.allclose(x, y))
     
+    def test_interpNd(self):
+        """ Make sure : the function interpNd works
+        """
+        boring_data = np.ones((5,5,5))
+        answer = nd.interpNd(boring_data, np.array([[2.3], [1.0], [3.9]]))
+        self.assertAllclose( answer , 1.0 )
+        
     def test_linear(self):
         """ Make sure : basic linear works
         """



More information about the Scipy-svn mailing list