Hi,
On Fri, 2008-07-25 at 12:43 -0500, John Hunter wrote:
> Also, I would rather not put the geometry functions in cbook, eg
> distances_along_curve and path_length and friends. Perhaps we should
> have some sort of geometry module where all these things live (there
> are some in mlab as well) but barring that I would rather have
> math/geometry stuff in mlab and general purpose helpers in cbook.
> Let's move all those before the release so we don't have to worry
> about API breakage later.
>
I started looking at this reorganization, but it seems to me that mlab
already has a number of functions that don't seem to have much to do
with its initial purpose - matlab compatibility. This was pretty
confusing for me when I initially started using matplotlib as the mlab
namespace was a mix of familiar and unfamiliar.
I decided instead to group several functions into a "numerical_methods"
module that includes most of my new cbook functions as well as some
linear interpolation and polygon stuff from mlab. I moved isvector to
mlab since this really is an imitation of a matlab function of the same
name. I have propagated these changes through to other functions in
matplotlib that use them and added remarks in CHANGELOG and API_CHANGES.
As these changes may be debatable, I haven't committed them. Instead, I
am attaching a patchset. I probably won't check email this weekend, so
I will let the powers that be decide what to do with this.
Along the way, I noticed these is some duplication in the examples
directory between pylab_examples and mpl_examples.
Cheers,
David
> JDH
--
**********************************
David M. Kaplan
Charge de Recherche 1
Institut de Recherche pour le Developpement
Centre de Recherche Halieutique Mediterraneenne et Tropicale
av. Jean Monnet
B.P. 171
34203 Sete cedex
France
Phone: +33 (0)4 99 57 32 27
Fax: +33 (0)4 99 57 32 95
http://www.ur097.ird.fr/team/dkaplan/index.html
**********************************
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')
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Matplotlib-devel mailing list
Matplotlib-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-devel
