[ENH]: contour at discontinuities: add max_jump parameter

[nschloe]

Problem

Follow up on #21309.

The problem is contour on discontinuous functions. Mpl makes the (sane) assumption that, if a quad's corners have the values 0.0 and 0.78, the level contour of 0.51 (for example) has to pass somewhere through the quad, and draws it. This is incorrect where the function is actually discontinuous, but of course mpl has no way of knowing that.

See this stackoverflow post or here.

Proposed solution

A simple solution for this would be to allow for a user-supplied float parameter, max_jump or discontinuity_threshold or something like that, defaulting to np.inf. If within a quad the function values differ by more than this threshold, a discontinuity is assumed and the level contour is assumed not to pass through the quad.

Additional context and prior art

No response

[nschloe]

I've just implemented this functionality in an alternative contouring algorithm in mplx. It's just two-hundred lines of code, no C++, the path-finding is done by networkx. Perhaps that's something worth including into mpl one day; anyway it can serve as an inspiration.

import matplotlib.pyplot as plt
import numpy as np

import mplx

x = np.linspace(-2.0, 2.0, 100)
y = np.linspace(-2.0, 2.0, 100)

X, Y = np.meshgrid(x, y)
Z = X + 1j * Y

vals = np.imag(np.log(Z))

# plt.contour(X, Y, vals, levels=[-2.0, -1.0, 0.0, 1.0, 2.0])  # draws wrong lines
mplx.contour(X, Y, vals, levels=[-2.0, -1.0, 0.0, 1.0, 2.0], max_jump=1.0)
mplx.discontour(X, Y, vals, min_jump=1.0, linestyle=":", color="r")

plt.gca().set_aspect("equal")
plt.savefig("out.png", bbox_inches="tight")
plt.show()

[jklymak]

Thanks for the proof-of-concept. I think this is a reasonable API. Whether it is reasonable to implement in the contouring code is a different issue. @ianthomas23 may be able to give guidance.

[ianthomas23]

I am glad that you have found a solution to your problem, and thanks for packaging it up so that it is easy to check. Here are some observations.

Firstly, you are not dealing with saddle quads correctly. If you contour z = [[0, 1], [1, 0.5]] this is what you see

Secondly, the difference between code that you just use yourself and code that is in matplotlib is that for the latter you need to cover all the functionality that the majority of users expect. Hence you need:

1. Support for filled contours, not just contour lines.
2. Support for z being a masked array, and the two corner_mask possibilities.

You will find that this will take significantly more than 200 lines of code.

Thirdly, performance. Using the test code below with a problem size of 1000x1000 takes 1.2 seconds using mpl on my dev machine, whereas your code takes 7.6 minutes, which is ~370 times slower. This is why the mpl algorithm is written in C++.

Test code:

import matplotlib.pyplot as plt
import mplx
import numpy as np
import time

n = 100

x, y = np.meshgrid(np.linspace(0, 1, n), np.linspace(0, 1, n))
rng = np.random.default_rng(824761)
z = rng.random(size=(n, n))
levels = [0.1, 0.3, 0.5, 0.7, 0.9]
print(f"Problem size {n} x {n}")

fig, axes = plt.subplots(ncols=2)

start = time.perf_counter()
axes[0].contour(x, y, z, levels)
mpl_s = time.perf_counter() - start
axes[0].plot(x, y, 'k-', x.T, y.T, 'k-', alpha=0.2)
axes[0].set_title("mpl")
print(f"mpl: {mpl_s:.3f} s")

plt.sca(axes[1])
start = time.perf_counter()
mplx.contour(x, y, z, levels)
mplx_s = time.perf_counter() - start
axes[1].plot(x, y, 'k-', x.T, y.T, 'k-', alpha=0.2)
axes[1].set_title("mplx")
print(f"mplx: {mplx_s:.3f} s")
print(f"Ratio: {mplx_s/mpl_s:.3f}")

plt.show()

[jklymak]

I think its a non-starter to not use the C++ code! I think the question is if an extra discontinuity check is possible in the C++ code. I'm not sure how hard it would be to propagate the check down, or if that is a reasonable new feature.

[ianthomas23]

The real question is whether there is somebody interested enough in the problem to do the actual work!

There are already 3 workarounds for this problem:

1. Split the grid into two along y=0 and separately contour y <= 0 and y >= 0. Care is needed to correctly set the z levels in the region of interest (x <= 0, y == 0).
2. Create the grid and calculate the contours in (radius, angle) space and transform back to (x, y) space.
3. Use tricontour.

[nschloe]

Jody is right, contour in mplx isn't meant as a drop-in replacement for plt.contour(), more a technical demo. I thought it might be interesting since it covers my use case (that was my main motivation writing it) and also demonstrates what's possible.

Some comments, nonetheless:

Support for filled contours, not just contour lines.

That won't be possible with with max_jump. Note that contours could just suddently end in the middle of nowhere. Contour lines of discontinuous graphs don't separate the domain in connected pieces. So, that whole thing would be for contour only.

Support for z being a masked array, and the two corner_mask possibilities.

I imagine adding this would be easy. The algorithm basically treats all quads individually, checks if a line goes through it and how, and only in the end forms the paths via a depth-first egde search in the resulting graph (via networkx). Removing some nodes/edges before that is no problem.

Thirdly, performance.

Interesting! As opposed to @jklymak, I'm actually convinved that, here like in many other places, NumPy code can be made faster than C++. ... You know, I'll actually look into this, just for the fun of it. Thanks for the test code!

There are already 3 workarounds for this problem:

I don't understand how any of the three suggestions should work. tricontour for example has just the same problem.

[ianthomas23]

I will try to explain the workarounds further then.

The idea with all 3 workarounds is that you want your grid points along x < 0, y = 0 to have two z values, one value when viewed/approached from y > 0 and the other value when viewed/approached from y < 0. Because this isn't supported, instead you can have two sets of coincident grid points along x < 0, y = 0. One set of these points has "positive" values and is only connected to the grid in the y > 0 direction, the other set of these points has "negative" values and is only connected to the grid in the y < 0 direction. Essentially you are cutting the grid at x < 0, y = 0 so that the grid is discontinuous, but by duplicating the points the data is continuous in each side of the grid and when rendered it all meets up.

For the 3 workarounds in order:

1. The topology of a quad grid is fixed so you cannot cut half-way along it. But you can cut all the way along it at y = 0 and duplicate the y = 0 grid points in each grid for all x.
2. A (radius, angle) quad grid when transformed into (x, y) is wrapped around the origin, and if you choose your angles correctly the ends meet along x < 0, y = 0.
3. You can have whatever topology you want with a triangle grid, so the grid can start at x < 0, y = 0 and progress around the origin, in the positive y direction for example, to meet duplicates of the original points from the negative y direction.

Does this make it clearer?

[nschloe]

If I understand correctly, all of these workarounds require knowing where the branch cut is. Correct?

[ianthomas23]

Yes, correct.

[github-actions]

This issue has been marked "inactive" because it has been 365 days since the last comment. If this issue is still present in recent Matplotlib releases, or the feature request is still wanted, please leave a comment and this label will be removed. If there are no updates in another 30 days, this issue will be automatically closed, but you are free to re-open or create a new issue if needed. We value issue reports, and this procedure is meant to help us resurface and prioritize issues that have not been addressed yet, not make them disappear. Thanks for your help!