10000 Merge pull request #24714 from anntzer/jacobian_angle · timhoffm/matplotlib@05fc1b3 · GitHub
[go: up one dir, main page]
Skip to content

Commit 05fc1b3

Browse files
QuLogicQuLogic
authored
Merge pull request matplotlib#24714 from anntzer/jacobian_angle
Improve handling of degenerate jacobians in non-rectilinear grids.
2 parents 27c5294 + 5ca9137 commit 05fc1b3

File tree

3 files changed

+102
-41
lines changed

3 files changed

+102
-41
lines changed

‎lib/mpl_toolkits/axisartist/floating_axes.py

Lines changed: 4 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -71,25 +71,19 @@ def trf_xy(x, y):
7171

7272
if self.nth_coord == 0:
7373
mask = (ymin <= yy0) & (yy0 <= ymax)
74-
(xx1, yy1), (dxx1, dyy1), (dxx2, dyy2) = \
75-
grid_helper_curvelinear._value_and_jacobian(
74+
(xx1, yy1), angle_normal, angle_tangent = \
75+
grid_helper_curvelinear._value_and_jac_angle(
7676
trf_xy, self.value, yy0[mask], (xmin, xmax), (ymin, ymax))
7777
labels = self._grid_info["lat_labels"]
7878

7979
elif self.nth_coord == 1:
8080
mask = (xmin <= xx0) & (xx0 <= xmax)
81-
(xx1, yy1), (dxx2, dyy2), (dxx1, dyy1) = \
82-
grid_helper_curvelinear._value_and_jacobian(
81+
(xx1, yy1), angle_tangent, angle_normal = \
82+
grid_helper_curvelinear._value_and_jac_angle(
8383
trf_xy, xx0[mask], self.value, (xmin, xmax), (ymin, ymax))
8484
labels = self._grid_info["lon_labels"]
8585

8686
labels = [l for l, m in zip(labels, mask) if m]
87-
88-
angle_normal = np.arctan2(dyy1, dxx1)
89-
angle_tangent = np.arctan2(dyy2, dxx2)
90-
mm = (dyy1 == 0) & (dxx1 == 0) # points with degenerate normal
91-
angle_normal[mm] = angle_tangent[mm] + np.pi / 2
92-
9387
tick_to_axes = self.get_tick_transform(axes) - axes.transAxes
9488
in_01 = functools.partial(
9589
mpl.transforms._interval_contains_close, (0, 1))

‎lib/mpl_toolkits/axisartist/grid_helper_curvelinear.py

Lines changed: 70 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -16,29 +16,75 @@
1616
from .grid_finder import GridFinder
1717

1818

19-
def _value_and_jacobian(func, xs, ys, xlims, ylims):
19+
def _value_and_jac_angle(func, xs, ys, xlim, ylim):
2020
"""
21-
Compute *func* and its derivatives along x and y at positions *xs*, *ys*,
22-
while ensuring that finite difference calculations don't try to evaluate
23-
values outside of *xlims*, *ylims*.
21+
Parameters
22+
----------
23+
func : callable
24+
A function that transforms the coordinates of a point (x, y) to a new coordinate
25+
system (u, v), and which can also take x and y as arrays of shape *shape* and
26+
returns (u, v) as a ``(2, shape)`` array.
27+
xs, ys : array-likes
28+
Points where *func* and its derivatives will be evaluated.
29+
xlim, ylim : pairs of floats
30+
(min, max) beyond which *func* should not be evaluated.
31+
32+
Returns
33+
-------
34+
val
35+
Value of *func* at each point of ``(xs, ys)``.
36+
thetas_dx
37+
Angles (in radians) defined by the (u, v) components of the numerically
38+
differentiated df/dx vector, at each point of ``(xs, ys)``. If needed, the
39+
differentiation step size is increased until at least one component of df/dx
40+
is nonzero, under the constraint of not going out of the *xlims*, *ylims*
41+
bounds. If the gridline at a point is actually null (and the angle is thus not
42+
well defined), the derivatives are evaluated after taking a small step along y;
43+
this ensures e.g. that the tick at r=0 on a radial axis of a polar plot is
44+
parallel with the ticks at r!=0.
45+
thetas_dy
46+
Like *thetas_dx*, but for df/dy.
2447
"""
25-
eps = np.finfo(float).eps ** (1/2) # see e.g. scipy.optimize.approx_fprime
48+
49+
shape = np.broadcast_shapes(np.shape(xs), np.shape(ys))
2650
val = func(xs, ys)
27-
# Take the finite difference step in the direction where the bound is the
28-
# furthest; the step size is min of epsilon and distance to that bound.
29-
xlo, xhi = sorted(xlims)
30-
dxlo = xs - xlo
31-
dxhi = xhi - xs
32-
xeps = (np.take([-1, 1], dxhi >= dxlo)
33-
* np.minimum(eps, np.maximum(dxlo, dxhi)))
34-
val_dx = func(xs + xeps, ys)
35-
ylo, yhi = sorted(ylims)
36-
dylo = ys - ylo
37-
dyhi = yhi - ys
38-
yeps = (np.take([-1, 1], dyhi >= dylo)
39-
* np.minimum(eps, np.maximum(dylo, dyhi)))
40-
val_dy = func(xs, ys + yeps)
41-
return (val, (val_dx - val) / xeps, (val_dy - val) / yeps)
51+
52+
# Take finite difference steps towards the furthest bound; the step size will be the
53+
# min of epsilon and the distance to that bound.
54+
eps0 = np.finfo(float).eps ** (1/2) # cf. scipy.optimize.approx_fprime
55+
56+
def calc_eps(vals, lim):
57+
lo, hi = sorted(lim)
58+
dlo = vals - lo
59+
dhi = hi - vals
60+
eps_max = np.maximum(dlo, dhi)
61+
eps = np.where(dhi >= dlo, 1, -1) * np.minimum(eps0, eps_max)
62+
return eps, eps_max
63+
64+
xeps, xeps_max = calc_eps(xs, xlim)
65+
yeps, yeps_max = calc_eps(ys, ylim)
66+
67+
def calc_thetas(dfunc, ps, eps_p0, eps_max, eps_q):
68+
thetas_dp = np.full(shape, np.nan)
69+
missing = np.full(shape, True)
70+
eps_p = eps_p0
71+
for it, eps_q in enumerate([0, eps_q]):
72+
while missing.any() and (abs(eps_p) < eps_max).any():
73+
if it == 0 and (eps_p > 1).any():
74+
break # Degenerate derivative, move a bit along the other coord.
75+
eps_p = np.minimum(eps_p, eps_max)
76+
df_x, df_y = (dfunc(eps_p, eps_q) - dfunc(0, eps_q)) / eps_p
77+
good = missing & ((df_x != 0) | (df_y != 0))
78+
thetas_dp[good] = np.arctan2(df_y, df_x)[good]
79+
missing &= ~good
80+
eps_p *= 2
81+
return thetas_dp
82+
83+
thetas_dx = calc_thetas(lambda eps_p, eps_q: func(xs + eps_p, ys + eps_q),
84+
xs, xeps, xeps_max, yeps)
85+
thetas_dy = calc_thetas(lambda eps_p, eps_q: func(xs + eps_q, ys + eps_p),
86+
ys, yeps, yeps_max, xeps)
87+
return (val, thetas_dx, thetas_dy)
4288

4389

4490
class FixedAxisArtistHelper(_FixedAxisArtistHelperBase):
@@ -167,12 +213,11 @@ def trf_xy(x, y):
167213
elif self.nth_coord == 1:
168214
xx0 = (xmin + xmax) / 2
169215
yy0 = self.value
170-
xy1, dxy1_dx, dxy1_dy = _value_and_jacobian(
216+
xy1, angle_dx, angle_dy = _value_and_jac_angle(
171217
trf_xy, xx0, yy0, (xmin, xmax), (ymin, ymax))
172218
p = axes.transAxes.inverted().transform(xy1)
173219
if 0 <= p[0] <= 1 and 0 <= p[1] <= 1:
174-
d = [dxy1_dy, dxy1_dx][self.nth_coord]
175-
return xy1, np.rad2deg(np.arctan2(*d[::-1]))
220+
return xy1, np.rad2deg([angle_dy, angle_dx][self.nth_coord])
176221
else:
177222
return None, None
178223

@@ -197,23 +242,17 @@ def trf_xy(x, y):
197242
# find angles
198243
if self.nth_coord == 0:
199244
mask = (e0 <= yy0) & (yy0 <= e1)
200-
(xx1, yy1), (dxx1, dyy1), (dxx2, dyy2) = _value_and_jacobian(
245+
(xx1, yy1), angle_normal, angle_tangent = _value_and_jac_angle(
201246
trf_xy, self.value, yy0[mask], (-np.inf, np.inf), (e0, e1))
202247
labels = self._grid_info["lat_labels"]
203248

204249
elif self.nth_coord == 1:
205250
mask = (e0 <= xx0) & (xx0 <= e1)
206-
(xx1, yy1), (dxx2, dyy2), (dxx1, dyy1) = _value_and_jacobian(
251+
(xx1, yy1), angle_tangent, angle_normal = _value_and_jac_angle(
207252
trf_xy, xx0[mask], self.value, (-np.inf, np.inf), (e0, e1))
208253
labels = self._grid_info["lon_labels"]
209254

210255
labels = [l for l, m in zip(labels, mask) if m]
211-
212-
angle_normal = np.arctan2(dyy1, dxx1)
213-
angle_tangent = np.arctan2(dyy2, dxx2)
214-
mm = (dyy1 == 0) & (dxx1 == 0) # points with degenerate normal
215-
angle_normal[mm] = angle_tangent[mm] + np.pi / 2
216-
217256
tick_to_axes = self.get_tick_transform(axes) - axes.transAxes
218257
in_01 = functools.partial(
219258
mpl.transforms._interval_contains_close, (0, 1))

‎lib/mpl_toolkits/axisartist/tests/test_floating_axes.py

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
import numpy as np
22

3+
import pytest
4+
35
import matplotlib.pyplot as plt
46
import matplotlib.projections as mprojections
57
import matplotlib.transforms as mtransforms
@@ -113,3 +115,29 @@ def test_axis_direction():
113115
ax.axis['y'] = ax.new_floating_axis(nth_coord=1, value=0,
114116
axis_direction='left')
115117
assert ax.axis['y']._axis_direction == 'left'
118+
119+
120+
def test_transform_with_zero_derivatives():
121+
# The transform is really a 45° rotation
122+
# tr(x, y) = x-y, x+y; inv_tr(u, v) = (u+v)/2, (u-v)/2
123+
# with an additional x->exp(-x**-2) on each coordinate.
124+
# Therefore all ticks should be at +/-45°, even the one at zero where the
125+
# transform derivatives are zero.
126+
127+
# at x=0, exp(-x* 9787 *-2)=0; div-by-zero can be ignored.
128+
@np.errstate(divide="ignore")
129+
def tr(x, y):
130+
return np.exp(-x**-2) - np.exp(-y**-2), np.exp(-x**-2) + np.exp(-y**-2)
131+
132+
def inv_tr(u, v):
133+
return (-np.log((u+v)/2))**(1/2), (-np.log((v-u)/2))**(1/2)
134+
135+
fig = plt.figure()
136+
ax = fig.add_subplot(
137+
axes_class=FloatingAxes, grid_helper=GridHelperCurveLinear(
138+
(tr, inv_tr), extremes=(0, 10, 0, 10)))
139+
fig.canvas.draw()
140+
141+
for k in ax.axis:
142+
for l, a in ax.axis[k].major_ticks.locs_angles:
143+
assert a % 90 == pytest.approx(45, abs=1e-3)

0 commit comments

Comments
 (0)
0