Index: CHANGELOG
===================================================================
--- CHANGELOG	(revision 5885)
+++ CHANGELOG	(working copy)
@@ -1,3 +1,7 @@
+2008-07-26 Reorganized cbook and mlab methods related to numerical
+	   calculations that have little to do with the goals of those two
+	   modules into a separate module numerical_methods.py - DMK
+
 2008-07-24 Deprecated (raise NotImplementedError) all the mlab2
            functions from matplotlib.mlab out of concern that some of
            them were not clean room implementations. JDH
Index: lib/matplotlib/cbook.py
===================================================================
--- lib/matplotlib/cbook.py	(revision 5885)
+++ lib/matplotlib/cbook.py	(working copy)
@@ -1156,47 +1156,6 @@
 
     return result
 
-def less_simple_linear_interpolation( x, y, xi, extrap=False ):
-    """
-    This function provides simple (but somewhat less so than
-    simple_linear_interpolation) linear interpolation.
-    simple_linear_interpolation will give a list of point between a
-    start and an end, while this does true linear interpolation at an
-    arbitrary set of points.
-
-    This is very inefficient linear interpolation meant to be used
-    only for a small number of points in relatively non-intensive use
-    cases.
-    """
-    if is_scalar(xi): xi = [xi]
-
-    x = np.asarray(x)
-    y = np.asarray(y)
-    xi = np.asarray(xi)
-
-    s = list(y.shape)
-    s[0] = len(xi)
-    yi = np.tile( np.nan, s )
-
-    for ii,xx in enumerate(xi):
-        bb = x == xx
-        if np.any(bb):
-            jj, = np.nonzero(bb)
-            yi[ii] = y[jj[0]]
-        elif xx<x[0]:
-            if extrap:
-                yi[ii] = y[0]
-        elif xx>x[-1]:
-            if extrap:
-                yi[ii] = y[-1]
-        else:
-            jj, = np.nonzero(x<xx)
-            jj = max(jj)
-
-            yi[ii] = y[jj] + (xx-x[jj])/(x[jj+1]-x[jj]) * (y[jj+1]-y[jj])
-
-    return yi
-
 def recursive_remove(path):
     if os.path.isdir(path):
         for fname in glob.glob(os.path.join(path, '*')) + glob.glob(os.path.join(path, '.*')):
@@ -1335,60 +1294,6 @@
     ic1 = breakpoints
     return np.concatenate((ic0[:, np.newaxis], ic1[:, np.newaxis]), axis=1)
 
-def isvector(X):
-    """
-    Like the Matlab (TM) function with the same name, returns true if
-    the supplied numpy array or matrix looks like a vector, meaning it
-    has a one non-singleton axis (i.e., it can have multiple axes, but
-    all must have length 1, except for one of them).
-
-    If you just want to see if the array has 1 axis, use X.ndim==1
-
-    """
-    return np.prod(X.shape)==np.max(X.shape)
-
-def vector_lengths( X, P=2., axis=None ):
-    """
-    Finds the length of a set of vectors in n dimensions.  This is
-    like the numpy norm function for vectors, but has the ability to
-    work over a particular axis of the supplied array or matrix.
-
-    Computes (sum((x_i)^P))^(1/P) for each {x_i} being the elements of X along
-    the given axis.  If *axis* is *None*, compute over all elements of X.
-    """
-    X = np.asarray(X)
-    return (np.sum(X**(P),axis=axis))**(1./P)
-
-def distances_along_curve( X ):
-    """
-    Computes the distance between a set of successive points in N dimensions.
-
-    where X is an MxN array or matrix.  The distances between successive rows
-    is computed.  Distance is the standard Euclidean distance.
-    """
-    X = np.diff( X, axis=0 )
-    return vector_lengths(X,axis=1)
-
-def path_length(X):
-    """
-    Computes the distance travelled along a polygonal curve in N dimensions.
-
-
-    where X is an MxN array or matrix.  Returns an array of length M consisting
-    of the distance along the curve at each point (i.e., the rows of X).
-    """
-    X = distances_along_curve(X)
-    return np.concatenate( (np.zeros(1), np.cumsum(X)) )
-
-def is_closed_polygon(X):
-    """
-    Tests whether first and last object in a sequence are the same.  These are
-    presumably coordinates on a polygonal curve, in which case this function
-    tests if that curve is closed.
-
-    """
-    return np.all(X[0] == X[-1])
-
 # a dict to cross-map linestyle arguments
 _linestyles = [('-', 'solid'),
     ('--', 'dashed'),
Index: lib/matplotlib/numerical_methods.py
===================================================================
--- lib/matplotlib/numerical_methods.py	(revision 0)
+++ lib/matplotlib/numerical_methods.py	(revision 0)
@@ -0,0 +1,285 @@
+"""
+A collection of utility functions that do various numerical or geometrical
+manipulations.  
+"""
+import numpy as np
+from numpy import ma
+import matplotlib.cbook as cbook
+
+##################################################
+# Linear interpolation algorithms
+##################################################
+def less_simple_linear_interpolation( x, y, xi, extrap=False ):
+    """
+    This function provides simple (but somewhat less so than
+    cbook.simple_linear_interpolation) linear interpolation.
+    simple_linear_interpolation will give a list of point between a
+    start and an end, while this does true linear interpolation at an
+    arbitrary set of points.
+
+    This is very inefficient linear interpolation meant to be used
+    only for a small number of points in relatively non-intensive use
+    cases.  For real linear interpolation, use scipy.
+    """
+    if cbook.is_scalar(xi): xi = [xi]
+
+    x = np.asarray(x)
+    y = np.asarray(y)
+    xi = np.asarray(xi)
+
+    s = list(y.shape)
+    s[0] = len(xi)
+    yi = np.tile( np.nan, s )
+
+    for ii,xx in enumerate(xi):
+        bb = x == xx
+        if np.any(bb):
+            jj, = np.nonzero(bb)
+            yi[ii] = y[jj[0]]
+        elif xx<x[0]:
+            if extrap:
+                yi[ii] = y[0]
+        elif xx>x[-1]:
+            if extrap:
+                yi[ii] = y[-1]
+        else:
+            jj, = np.nonzero(x<xx)
+            jj = max(jj)
+
+            yi[ii] = y[jj] + (xx-x[jj])/(x[jj+1]-x[jj]) * (y[jj+1]-y[jj])
+
+    return yi
+
+def slopes(x,y):
+    """
+    SLOPES calculate the slope y'(x) Given data vectors X and Y SLOPES
+    calculates Y'(X), i.e the slope of a curve Y(X). The slope is
+    estimated using the slope obtained from that of a parabola through
+    any three consecutive points.
+
+    This method should be superior to that described in the appendix
+    of A CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russel
+    W. Stineman (Creative Computing July 1980) in at least one aspect:
+
+    Circles for interpolation demand a known aspect ratio between x-
+    and y-values.  For many functions, however, the abscissa are given
+    in different dimensions, so an aspect ratio is completely
+    arbitrary.
+
+    The parabola method gives very similar results to the circle
+    method for most regular cases but behaves much better in special
+    cases
+
+    Norbert Nemec, Institute of Theoretical Physics, University or
+    Regensburg, April 2006 Norbert.Nemec at physik.uni-regensburg.de
+
+    (inspired by a original implementation by Halldor Bjornsson,
+    Icelandic Meteorological Office, March 2006 halldor at vedur.is)
+    """
+    # Cast key variables as float.
+    x=np.asarray(x, np.float_)
+    y=np.asarray(y, np.float_)
+
+    yp=np.zeros(y.shape, np.float_)
+
+    dx=x[1:] - x[:-1]
+    dy=y[1:] - y[:-1]
+    dydx = dy/dx
+    yp[1:-1] = (dydx[:-1] * dx[1:] + dydx[1:] * dx[:-1])/(dx[1:] + dx[:-1])
+    yp[0] = 2.0 * dy[0]/dx[0] - yp[1]
+    yp[-1] = 2.0 * dy[-1]/dx[-1] - yp[-2]
+    return yp
+
+
+def stineman_interp(xi,x,y,yp=None):
+    """
+    STINEMAN_INTERP Well behaved data interpolation.  Given data
+    vectors X and Y, the slope vector YP and a new abscissa vector XI
+    the function stineman_interp(xi,x,y,yp) uses Stineman
+    interpolation to calculate a vector YI corresponding to XI.
+
+    Here's an example that generates a coarse sine curve, then
+    interpolates over a finer abscissa:
+
+      x = linspace(0,2*pi,20);  y = sin(x); yp = cos(x)
+      xi = linspace(0,2*pi,40);
+      yi = stineman_interp(xi,x,y,yp);
+      plot(x,y,'o',xi,yi)
+
+    The interpolation method is described in the article A
+    CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell
+    W. Stineman. The article appeared in the July 1980 issue of
+    Creative Computing with a note from the editor stating that while
+    they were
+
+      not an academic journal but once in a while something serious
+      and original comes in adding that this was
+      "apparently a real solution" to a well known problem.
+
+    For yp=None, the routine automatically determines the slopes using
+    the "slopes" routine.
+
+    X is assumed to be sorted in increasing order
+
+    For values xi[j] < x[0] or xi[j] > x[-1], the routine tries a
+    extrapolation.  The relevance of the data obtained from this, of
+    course, questionable...
+
+    original implementation by Halldor Bjornsson, Icelandic
+    Meteorolocial Office, March 2006 halldor at vedur.is
+
+    completely reworked and optimized for Python by Norbert Nemec,
+    Institute of Theoretical Physics, University or Regensburg, April
+    2006 Norbert.Nemec at physik.uni-regensburg.de
+
+    """
+
+    # Cast key variables as float.
+    x=np.asarray(x, np.float_)
+    y=np.asarray(y, np.float_)
+    assert x.shape == y.shape
+    N=len(y)
+
+    if yp is None:
+        yp = slopes(x,y)
+    else:
+        yp=np.asarray(yp, np.float_)
+
+    xi=np.asarray(xi, np.float_)
+    yi=np.zeros(xi.shape, np.float_)
+
+    # calculate linear slopes
+    dx = x[1:] - x[:-1]
+    dy = y[1:] - y[:-1]
+    s = dy/dx  #note length of s is N-1 so last element is #N-2
+
+    # find the segment each xi is in
+    # this line actually is the key to the efficiency of this implementation
+    idx = np.searchsorted(x[1:-1], xi)
+
+    # now we have generally: x[idx[j]] <= xi[j] <= x[idx[j]+1]
+    # except at the boundaries, where it may be that xi[j] < x[0] or xi[j] > x[-1]
+
+    # the y-values that would come out from a linear interpolation:
+    sidx = s.take(idx)
+    xidx = x.take(idx)
+    yidx = y.take(idx)
+    xidxp1 = x.take(idx+1)
+    yo = yidx + sidx * (xi - xidx)
+
+    # the difference that comes when using the slopes given in yp
+    dy1 = (yp.take(idx)- sidx) * (xi - xidx)       # using the yp slope of the left point
+    dy2 = (yp.take(idx+1)-sidx) * (xi - xidxp1) # using the yp slope of the right point
+
+    dy1dy2 = dy1*dy2
+    # The following is optimized for Python. The solution actually
+    # does more calculations than necessary but exploiting the power
+    # of numpy, this is far more efficient than coding a loop by hand
+    # in Python
+    yi = yo + dy1dy2 * np.choose(np.array(np.sign(dy1dy2), np.int32)+1,
+                                 ((2*xi-xidx-xidxp1)/((dy1-dy2)*(xidxp1-xidx)),
+                                  0.0,
+                                  1/(dy1+dy2),))
+    return yi
+
+##################################################
+# Code related to things in and around polygons
+##################################################
+def inside_poly(points, verts):
+    """
+    points is a sequence of x,y points
+    verts is a sequence of x,y vertices of a poygon
+
+    return value is a sequence of indices into points for the points
+    that are inside the polygon
+    """
+    res, =  np.nonzero(nxutils.points_inside_poly(points, verts))
+    return res
+
+def poly_below(xmin, xs, ys):
+    """
+    given a sequence of xs and ys, return the vertices of a polygon
+    that has a horzontal base at xmin and an upper bound at the ys.
+    xmin is a scalar.
+
+    intended for use with Axes.fill, eg
+    xv, yv = poly_below(0, x, y)
+    ax.fill(xv, yv)
+    """
+    xs = np.asarray(xs)
+    ys = np.asarray(ys)
+    Nx = len(xs)
+    Ny = len(ys)
+    assert(Nx==Ny)
+    x = xmin*np.ones(2*Nx)
+    y = np.ones(2*Nx)
+    x[:Nx] = xs
+    y[:Nx] = ys
+    y[Nx:] = ys[::-1]
+    return x, y
+
+
+def poly_between(x, ylower, yupper):
+    """
+    given a sequence of x, ylower and yupper, return the polygon that
+    fills the regions between them.  ylower or yupper can be scalar or
+    iterable.  If they are iterable, they must be equal in length to x
+
+    return value is x, y arrays for use with Axes.fill
+    """
+    Nx = len(x)
+    if not cbook.iterable(ylower):
+        ylower = ylower*np.ones(Nx)
+
+    if not cbook.iterable(yupper):
+        yupper = yupper*np.ones(Nx)
+
+    x = np.concatenate( (x, x[::-1]) )
+    y = np.concatenate( (yupper, ylower[::-1]) )
+    return x,y
+
+def is_closed_polygon(X):
+    """
+    Tests whether first and last object in a sequence are the same.  These are
+    presumably coordinates on a polygonal curve, in which case this function
+    tests if that curve is closed.
+
+    """
+    return np.all(X[0] == X[-1])
+
+##################################################
+# Vector and path length geometry calculations
+##################################################
+def vector_lengths( X, P=2., axis=None ):
+    """
+    Finds the length of a set of vectors in n dimensions.  This is
+    like the numpy norm function for vectors, but has the ability to
+    work over a particular axis of the supplied array or matrix.
+
+    Computes (sum((x_i)^P))^(1/P) for each {x_i} being the elements of X along
+    the given axis.  If *axis* is *None*, compute over all elements of X.
+    """
+    X = np.asarray(X)
+    return (np.sum(X**(P),axis=axis))**(1./P)
+
+def distances_along_curve( X ):
+    """
+    Computes the distance between a set of successive points in N dimensions.
+
+    where X is an MxN array or matrix.  The distances between successive rows
+    is computed.  Distance is the standard Euclidean distance.
+    """
+    X = np.diff( X, axis=0 )
+    return vector_lengths(X,axis=1)
+
+def path_length(X):
+    """
+    Computes the distance travelled along a polygonal curve in N dimensions.
+
+
+    where X is an MxN array or matrix.  Returns an array of length M consisting
+    of the distance along the curve at each point (i.e., the rows of X).
+    """
+    X = distances_along_curve(X)
+    return np.concatenate( (np.zeros(1), np.cumsum(X)) )
+
Index: lib/matplotlib/contour.py
===================================================================
--- lib/matplotlib/contour.py	(revision 5885)
+++ lib/matplotlib/contour.py	(working copy)
@@ -16,6 +16,7 @@
 import matplotlib.font_manager as font_manager
 import matplotlib.text as text
 import matplotlib.cbook as cbook
+import matplotlib.numerical_methods as numerical_methods
 
 # Import needed for adding manual selection capability to clabel
 from matplotlib.blocking_input import BlockingContourLabeler
@@ -339,7 +340,7 @@
         not empty (lc defaults to the empty list if None).  *spacing*
         is the space around the label in pixels to leave empty.
 
-        Do both of these tasks at once to avoid calling cbook.path_length
+        Do both of these tasks at once to avoid calling numerical_methods.path_length
         multiple times, which is relatively costly.
 
         The method used here involves calculating the path length
@@ -353,7 +354,7 @@
         hlw = lw/2.0
 
         # Check if closed and, if so, rotate contour so label is at edge
-        closed = cbook.is_closed_polygon(slc)
+        closed = numerical_methods.is_closed_polygon(slc)
         if closed:
             slc = np.r_[ slc[ind:-1], slc[:ind+1] ]
 
@@ -363,7 +364,7 @@
             ind = 0
 
         # Path length in pixel space
-        pl = cbook.path_length(slc)
+        pl = numerical_methods.path_length(slc)
         pl = pl-pl[ind]
 
         # Use linear interpolation to get points around label
@@ -373,7 +374,7 @@
         else:
             dp = np.zeros_like(xi)
 
-        ll = cbook.less_simple_linear_interpolation( pl, slc, dp+xi,
+        ll = numerical_methods.less_simple_linear_interpolation( pl, slc, dp+xi,
                                                      extrap=True )
 
         # get vector in pixel space coordinates from one point to other
@@ -399,16 +400,16 @@
             xi = dp + xi + np.array([-spacing,spacing])
 
             # Get indices near points of interest
-            I = cbook.less_simple_linear_interpolation(
+            I = numerical_methods.less_simple_linear_interpolation(
                 pl, np.arange(len(pl)), xi, extrap=False )
 
             # If those indices aren't beyond contour edge, find x,y
             if (not np.isnan(I[0])) and int(I[0])<>I[0]:
-                xy1 = cbook.less_simple_linear_interpolation(
+                xy1 = numerical_methods.less_simple_linear_interpolation(
                     pl, lc, [ xi[0] ] )
 
             if (not np.isnan(I[1])) and int(I[1])<>I[1]:
-                xy2 = cbook.less_simple_linear_interpolation(
+                xy2 = numerical_methods.less_simple_linear_interpolation(
                     pl, lc, [ xi[1] ] )
 
             # Make integer
@@ -476,7 +477,7 @@
                 # zero in print_label and locate_label.  Other than these
                 # functions, this is not necessary and should probably be
                 # eventually removed.
-                if cbook.is_closed_polygon( lc ):
+                if numerical_methods.is_closed_polygon( lc ):
                     slc = np.r_[ slc0, slc0[1:2,:] ]
                 else:
                     slc = slc0
Index: lib/matplotlib/pylab.py
===================================================================
--- lib/matplotlib/pylab.py	(revision 5885)
+++ lib/matplotlib/pylab.py	(working copy)
@@ -225,8 +225,12 @@
         sqrtm, prctile, center_matrix, rk4, exp_safe, amap,\
         sum_flat, mean_flat, rms_flat, l1norm, l2norm, norm, frange,\
         diagonal_matrix, base_repr, binary_repr, log2, ispower2,\
-        bivariate_normal, load, save, stineman_interp
+        bivariate_normal, load, save
 
+from matplotlib.numerical_methods import stineman_interp, slopes, \
+    stineman_interp, inside_poly, poly_below, poly_between, \
+    is_closed_polygon, path_length, distances_along_curve, vector_lengths
+
 from numpy import *
 from numpy.fft import *
 from numpy.random import *
@@ -237,10 +241,10 @@
      find, longest_contiguous_ones, longest_ones, prepca, prctile, prctile_rank, \
      center_matrix, rk4, bivariate_normal, get_xyz_where, get_sparse_matrix, dist, \
      dist_point_to_segment, segments_intersect, fftsurr, liaupunov, movavg, \
-     save, load, slopes, stineman_interp, inside_poly, poly_below, poly_between, exp_safe, \
+     save, load, exp_safe, \
      amap, rms_flat, l1norm, l2norm, norm_flat, frange, diagonal_matrix, identity, \
      base_repr, binary_repr, log2, ispower2, fromfunction_kw, rem, norm, orth, rank, sqrtm,\
-     mfuncC, approx_real, rec_append_field, rec_drop_fields, rec_join, csv2rec, rec2csv
+     mfuncC, approx_real, rec_append_field, rec_drop_fields, rec_join, csv2rec, rec2csv, isvector
 
 
 
Index: lib/matplotlib/mlab.py
===================================================================
--- lib/matplotlib/mlab.py	(revision 5885)
+++ lib/matplotlib/mlab.py	(working copy)
@@ -1292,191 +1292,6 @@
     else: return X
 
 
-def slopes(x,y):
-    """
-    SLOPES calculate the slope y'(x) Given data vectors X and Y SLOPES
-    calculates Y'(X), i.e the slope of a curve Y(X). The slope is
-    estimated using the slope obtained from that of a parabola through
-    any three consecutive points.
-
-    This method should be superior to that described in the appendix
-    of A CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russel
-    W. Stineman (Creative Computing July 1980) in at least one aspect:
-
-    Circles for interpolation demand a known aspect ratio between x-
-    and y-values.  For many functions, however, the abscissa are given
-    in different dimensions, so an aspect ratio is completely
-    arbitrary.
-
-    The parabola method gives very similar results to the circle
-    method for most regular cases but behaves much better in special
-    cases
-
-    Norbert Nemec, Institute of Theoretical Physics, University or
-    Regensburg, April 2006 Norbert.Nemec at physik.uni-regensburg.de
-
-    (inspired by a original implementation by Halldor Bjornsson,
-    Icelandic Meteorological Office, March 2006 halldor at vedur.is)
-    """
-    # Cast key variables as float.
-    x=np.asarray(x, np.float_)
-    y=np.asarray(y, np.float_)
-
-    yp=np.zeros(y.shape, np.float_)
-
-    dx=x[1:] - x[:-1]
-    dy=y[1:] - y[:-1]
-    dydx = dy/dx
-    yp[1:-1] = (dydx[:-1] * dx[1:] + dydx[1:] * dx[:-1])/(dx[1:] + dx[:-1])
-    yp[0] = 2.0 * dy[0]/dx[0] - yp[1]
-    yp[-1] = 2.0 * dy[-1]/dx[-1] - yp[-2]
-    return yp
-
-
-def stineman_interp(xi,x,y,yp=None):
-    """
-    STINEMAN_INTERP Well behaved data interpolation.  Given data
-    vectors X and Y, the slope vector YP and a new abscissa vector XI
-    the function stineman_interp(xi,x,y,yp) uses Stineman
-    interpolation to calculate a vector YI corresponding to XI.
-
-    Here's an example that generates a coarse sine curve, then
-    interpolates over a finer abscissa:
-
-      x = linspace(0,2*pi,20);  y = sin(x); yp = cos(x)
-      xi = linspace(0,2*pi,40);
-      yi = stineman_interp(xi,x,y,yp);
-      plot(x,y,'o',xi,yi)
-
-    The interpolation method is described in the article A
-    CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell
-    W. Stineman. The article appeared in the July 1980 issue of
-    Creative Computing with a note from the editor stating that while
-    they were
-
-      not an academic journal but once in a while something serious
-      and original comes in adding that this was
-      "apparently a real solution" to a well known problem.
-
-    For yp=None, the routine automatically determines the slopes using
-    the "slopes" routine.
-
-    X is assumed to be sorted in increasing order
-
-    For values xi[j] < x[0] or xi[j] > x[-1], the routine tries a
-    extrapolation.  The relevance of the data obtained from this, of
-    course, questionable...
-
-    original implementation by Halldor Bjornsson, Icelandic
-    Meteorolocial Office, March 2006 halldor at vedur.is
-
-    completely reworked and optimized for Python by Norbert Nemec,
-    Institute of Theoretical Physics, University or Regensburg, April
-    2006 Norbert.Nemec at physik.uni-regensburg.de
-
-    """
-
-    # Cast key variables as float.
-    x=np.asarray(x, np.float_)
-    y=np.asarray(y, np.float_)
-    assert x.shape == y.shape
-    N=len(y)
-
-    if yp is None:
-        yp = slopes(x,y)
-    else:
-        yp=np.asarray(yp, np.float_)
-
-    xi=np.asarray(xi, np.float_)
-    yi=np.zeros(xi.shape, np.float_)
-
-    # calculate linear slopes
-    dx = x[1:] - x[:-1]
-    dy = y[1:] - y[:-1]
-    s = dy/dx  #note length of s is N-1 so last element is #N-2
-
-    # find the segment each xi is in
-    # this line actually is the key to the efficiency of this implementation
-    idx = np.searchsorted(x[1:-1], xi)
-
-    # now we have generally: x[idx[j]] <= xi[j] <= x[idx[j]+1]
-    # except at the boundaries, where it may be that xi[j] < x[0] or xi[j] > x[-1]
-
-    # the y-values that would come out from a linear interpolation:
-    sidx = s.take(idx)
-    xidx = x.take(idx)
-    yidx = y.take(idx)
-    xidxp1 = x.take(idx+1)
-    yo = yidx + sidx * (xi - xidx)
-
-    # the difference that comes when using the slopes given in yp
-    dy1 = (yp.take(idx)- sidx) * (xi - xidx)       # using the yp slope of the left point
-    dy2 = (yp.take(idx+1)-sidx) * (xi - xidxp1) # using the yp slope of the right point
-
-    dy1dy2 = dy1*dy2
-    # The following is optimized for Python. The solution actually
-    # does more calculations than necessary but exploiting the power
-    # of numpy, this is far more efficient than coding a loop by hand
-    # in Python
-    yi = yo + dy1dy2 * np.choose(np.array(np.sign(dy1dy2), np.int32)+1,
-                                 ((2*xi-xidx-xidxp1)/((dy1-dy2)*(xidxp1-xidx)),
-                                  0.0,
-                                  1/(dy1+dy2),))
-    return yi
-
-def inside_poly(points, verts):
-    """
-    points is a sequence of x,y points
-    verts is a sequence of x,y vertices of a poygon
-
-    return value is a sequence of indices into points for the points
-    that are inside the polygon
-    """
-    res, =  np.nonzero(nxutils.points_inside_poly(points, verts))
-    return res
-
-def poly_below(xmin, xs, ys):
-    """
-    given a sequence of xs and ys, return the vertices of a polygon
-    that has a horzontal base at xmin and an upper bound at the ys.
-    xmin is a scalar.
-
-    intended for use with Axes.fill, eg
-    xv, yv = poly_below(0, x, y)
-    ax.fill(xv, yv)
-    """
-    xs = np.asarray(xs)
-    ys = np.asarray(ys)
-    Nx = len(xs)
-    Ny = len(ys)
-    assert(Nx==Ny)
-    x = xmin*np.ones(2*Nx)
-    y = np.ones(2*Nx)
-    x[:Nx] = xs
-    y[:Nx] = ys
-    y[Nx:] = ys[::-1]
-    return x, y
-
-
-def poly_between(x, ylower, yupper):
-    """
-    given a sequence of x, ylower and yupper, return the polygon that
-    fills the regions between them.  ylower or yupper can be scalar or
-    iterable.  If they are iterable, they must be equal in length to x
-
-    return value is x, y arrays for use with Axes.fill
-    """
-    Nx = len(x)
-    if not cbook.iterable(ylower):
-        ylower = ylower*np.ones(Nx)
-
-    if not cbook.iterable(yupper):
-        yupper = yupper*np.ones(Nx)
-
-    x = np.concatenate( (x, x[::-1]) )
-    y = np.concatenate( (yupper, ylower[::-1]) )
-    return x,y
-
 ### the following code was written and submitted by Fernando Perez
 ### from the ipython numutils package under a BSD license
 # begin fperez functions
@@ -1755,6 +1570,18 @@
     else:
         return len(bin_n)
 
+def isvector(X):
+    """
+    Like the Matlab (TM) function with the same name, returns true if
+    the supplied numpy array or matrix looks like a vector, meaning it
+    has a one non-singleton axis (i.e., it can have multiple axes, but
+    all must have length 1, except for one of them).
+
+    If you just want to see if the array has 1 axis, use X.ndim==1
+
+    """
+    return np.prod(X.shape)==np.max(X.shape)
+
 #from numpy import fromfunction as fromfunction_kw
 def fromfunction_kw(function, dimensions, **kwargs):
     """Drop-in replacement for fromfunction() from numpy
Index: API_CHANGES
===================================================================
--- API_CHANGES	(revision 5885)
+++ API_CHANGES	(working copy)
@@ -19,6 +19,10 @@
   maintained for the moment (in addition to their renamed versions), but they
   are depricated and will eventually be removed.
 
+* Moved several function in mlab.py into a separate module 
+  numerical_methods.py because they were unrelated to the initial purpose of
+  mlab and appeared more coherent elsewhere.
+
 Changes for 0.98.1
 ==================
 
Index: examples/pylab_examples/interp_demo.py
===================================================================
--- examples/pylab_examples/interp_demo.py	(revision 5885)
+++ examples/pylab_examples/interp_demo.py	(working copy)
@@ -1,6 +1,6 @@
 from matplotlib.pyplot import figure, show
 from numpy import pi, sin, linspace
-from matplotlib.mlab import stineman_interp
+from matplotlib.numerical_methods import stineman_interp
 
 x = linspace(0,2*pi,20);
 y = sin(x); yp = None
Index: examples/pylab_examples/fill_between.py
===================================================================
--- examples/pylab_examples/fill_between.py	(revision 5885)
+++ examples/pylab_examples/fill_between.py	(working copy)
@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-import matplotlib.mlab as mlab
+import matplotlib.numerical_methods as numerical_methods
 from pylab import figure, show
 import numpy as np
 
@@ -13,15 +13,15 @@
 ax3 = fig.add_subplot(313)
 
 
-xs, ys = mlab.poly_between(x, 0, y1)
+xs, ys = numerical_methods.poly_between(x, 0, y1)
 ax.fill(xs, ys)
 ax.set_ylabel('between y1 and 0')
 
-xs, ys = mlab.poly_between(x, y1, 1)
+xs, ys = numerical_methods.poly_between(x, y1, 1)
 ax2.fill(xs, ys)
 ax2.set_ylabel('between y1 and 1')
 
-xs, ys = mlab.poly_between(x, y1, y2)
+xs, ys = numerical_methods.poly_between(x, y1, y2)
 ax3.fill(xs, ys)
 ax3.set_ylabel('between y1 and y2')
 ax3.set_xlabel('x')
