From 94be20ad1e12a394ad2e6b01a9272fd2d705fc3e Mon Sep 17 00:00:00 2001 From: Antony Lee Date: Sat, 25 Nov 2017 22:15:43 -0800 Subject: [PATCH] Improve less_simple_linear_interpolation. - Rewrite less_simple_linear_interpolation to use np.interp and np.polyfit (for the extrapolation part), document it properly and note the differences with np.interp. - Deprecate extrapolation-less mode (which is well covered by np.interp). - Simplify its usage in contour: only one place needs extrapolation, so the others can use np.interp; -1 is a perfectly reasonable marker for out of bounds which avoids having to dance with the errstate. --- lib/matplotlib/contour.py | 57 +++++++++++--------------------- lib/matplotlib/mlab.py | 68 ++++++++++++++++++++------------------- 2 files changed, 54 insertions(+), 71 deletions(-) diff --git a/lib/matplotlib/contour.py b/lib/matplotlib/contour.py index ecbde9ae91b5..db72b31331db 100644 --- a/lib/matplotlib/contour.py +++ b/lib/matplotlib/contour.py @@ -417,25 +417,17 @@ def calc_label_rot_and_inline(self, slc, ind, lw, lc=None, spacing=5): else: dp = np.zeros_like(xi) - ll = mlab.less_simple_linear_interpolation(pl, slc, dp + xi, - extrap=True) - + ll = mlab.less_simple_linear_interpolation(pl, slc, dp + xi) # get vector in pixel space coordinates from one point to other dd = np.diff(ll, axis=0).ravel() - # Get angle of vector - must be calculated in pixel space for - # text rotation to work correctly - if np.all(dd == 0): # Must deal with case of zero length label - rotation = 0.0 - else: - rotation = np.rad2deg(np.arctan2(dd[1], dd[0])) + # text rotation to work correctly (for zero-sized labels, atan2(0, 0) + # is still well-defined). + rotation = np.rad2deg(np.arctan2(dd[1], dd[0])) if self.rightside_up: # Fix angle so text is never upside-down - if rotation > 90: - rotation = rotation - 180.0 - if rotation < -90: - rotation = 180.0 + rotation + rotation = (rotation + 90) % 180 - 90 # Break contour if desired nlc = [] @@ -443,37 +435,26 @@ def calc_label_rot_and_inline(self, slc, ind, lw, lc=None, spacing=5): # Expand range by spacing xi = dp + xi + np.array([-spacing, spacing]) - # Get indices near points of interest - I = mlab.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 = mlab.less_simple_linear_interpolation( - pl, lc, [xi[0]]) - - if (not np.isnan(I[1])) and int(I[1]) != I[1]: - xy2 = mlab.less_simple_linear_interpolation( - pl, lc, [xi[1]]) - - # Round to integer values but keep as float - # To allow check against nan below - # Ignore nans here to avoid throwing an error on Appveyor build - # (can possibly be removed when build uses numpy 1.13) - with np.errstate(invalid='ignore'): - I = [np.floor(I[0]), np.ceil(I[1])] + # Get (integer) indices near points of interest; use -1 as marker + # for out of bounds. + I = np.interp(xi, pl, np.arange(len(pl)), left=-1, right=-1) + I = [np.floor(I[0]).astype(int), np.ceil(I[1]).astype(int)] + if I[0] != -1: + xy1 = mlab.less_simple_linear_interpolation(pl, lc, [xi[0]]) + if I[1] != -1: + xy2 = mlab.less_simple_linear_interpolation(pl, lc, [xi[1]]) # Actually break contours if closed: # This will remove contour if shorter than label - if np.all(~np.isnan(I)): - nlc.append(np.r_[xy2, lc[int(I[1]):int(I[0]) + 1], xy1]) + if np.all(I != -1): + nlc.append(np.row_stack([xy2, lc[I[1]:I[0]+1], xy1])) else: # These will remove pieces of contour if they have length zero - if not np.isnan(I[0]): - nlc.append(np.r_[lc[:int(I[0]) + 1], xy1]) - if not np.isnan(I[1]): - nlc.append(np.r_[xy2, lc[int(I[1]):]]) + if I[0] != -1: + nlc.append(np.row_stack([lc[:I[0]+1], xy1])) + if I[1] != -1: + nlc.append(np.row_stack([xy2, lc[I[1]:]])) # The current implementation removes contours completely # covered by labels. Uncomment line below to keep diff --git a/lib/matplotlib/mlab.py b/lib/matplotlib/mlab.py index 3664e869814e..78614820a102 100644 --- a/lib/matplotlib/mlab.py +++ b/lib/matplotlib/mlab.py @@ -3404,44 +3404,46 @@ def griddata(x, y, z, xi, yi, interp='nn'): ################################################## # Linear interpolation algorithms ################################################## -def less_simple_linear_interpolation(x, y, xi, extrap=False): +_sentinel = object() +def less_simple_linear_interpolation(x, y, xi, extrap=_sentinel): """ - This function provides simple (but somewhat less so than - :func:`cbook.simple_linear_interpolation`) linear interpolation. - :func:`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. + Linearly interpolate and extrapolate a function. - 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. + This function differs from `np.interp` by its ability to linearly + extrapolate beyond the lowest and highest x, and to handle multidimensional + output. + + Parameters + ---------- + x : a 1D array-like + The x-values where the interpolated function is known. + y : a 2D array-like + The values of the interpolated function. + xi : a 1D array-like + The x-values where to interpolate or extrapolate the function. + + Returns + ------- + The linearly interpolated and extrapolated values as a 2D array. """ x = np.asarray(x) y = np.asarray(y) - xi = np.atleast_1d(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 + xi = np.asarray(xi) + out = np.column_stack([np.interp(xi, x, col, left=np.nan, right=np.nan) + for col in y.T]) + if extrap is not _sentinel: + cbook.warn_deprecated( + "2.2", "The 'extrap' keyword argument to " + "'less_simple_linear_interpolation' is deprecated " + "(use np.interp instead).") + if extrap: + mask = xi < x[0] + if np.any(mask): + out[mask] = np.polyval(np.polyfit(x[:2], y[:2], 1), xi[mask]) + mask = xi > x[-1] + if np.any(mask): + out[mask] = np.polyval(np.polyfit(x[-2:], y[-2:], 1), xi[mask]) + return out def slopes(x, y):