Skip to content

Improve less_simple_linear_interpolation. #9859

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

Improve less_simple_linear_interpolation. #9859

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 19 additions & 38 deletions lib/matplotlib/contour.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -417,63 +417,44 @@ 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 = []
if len(lc):
# 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
Expand Down
68 changes: 35 additions & 33 deletions lib/matplotlib/mlab.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The original is using the same extrapolation method as np.interp uses by default: setting the past-the-end values to the end values of the input array. Are you changing the behavior because what mpl actually needs is linear extrapolation?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wow, I actually missed that point because my understanding of the contour code is that it indeed needs extrapolation (but I need to look at it again and run some examples). If we don't actually need extrapolation then we can definitely get rid of this function in favor of using interp.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I convinced myself that linear extrapolation was not needed because the call to calc_label_rot_and_inline is protected by print_label which is explicitly there to check that the contour line is long enough to host the label text, so we don't expect overflow anyways (so flat extrapolation is as good as anything).

Rewriting in terms of np.interp...

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):
Expand Down