Colormaps and contours for complex_plot

[davidlowryduda]

This ticket rewrites complex_plot to accept an optional matplotlib-compatible colormap and optional arguments that change the domain coloring to include "contours". I've been using this code in https://github.com/davidlowryduda/phase_mag_plot, and a (slower, older) variant of this code were used in my paper https://arxiv.org/abs/2002.05234.

Typical usage would go like

sage: complex_plot(lambda x: x^2 - x + 2, (-10, 10), (-10, 10), aspect_ratio=1
                   cmap='viridis')

To include magnitude-based contours, one would write something like

sage: complex_plot(lambda x: x^2 - x + 2, (-10, 10), (-10, 10),
                   contoured=True)

This also includes the domain coloring tiling patterns from
Wegert and Semmler, Phase plots of complex func-
tions: a journey in illustration (Notices AMS 2010). This would be written like

sage: complex_plot(lambda x: x^2 - x + 2, (-10, 10), (-10, 10),
                   tiled=True)

Colormaps and contours can be combined or used separately. Without these arguments, the previous complex_plot functionality continues unchanged.

Component: graphics

Keywords: complex plotting

Author: David Lowry-Duda

Branch/Commit: e1b7ca9

Reviewer: Markus Wageringel, Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/33416

[davidlowryduda]

Author: David Lowry-Duda

[davidlowryduda]

Changed keywords from none to complex plotting

[davidlowryduda]

Branch: u/DavidLowry/colormaps_and_contours_for_complex_plot

[davidlowryduda]

Commit: 57414a3

[davidlowryduda]

New commits:

5ba67cb	Rename complex_to_rgb to direct_complex_to_rgb
803cc28	Add colormaps and contours to complex plots
2be3a30	Cythonize complex_plotting loops
2d2fe45	Format documentation correctly
57414a3	Merge branch 'colorcomplexplot' into t/33416/colormaps_and_contours_for_complex_plot

[tscrim]

comment:4

Thank you for the contribution. Some comments:

• I am not sure I like 'default' not being the default in Sage. Perhaps 'matplotlib'? If we can't agree on a better name, this will be okay, but I might want to put some effort towards this.

• Is there an easy way to extend the contours to have different user input? I am thinking two parameters: one for choosing linear versus log scaling, and then another for the intervals. I am thinking about if someone is looking at where a function changes more slowly (a more trivial example, but something like f(x) = z^2).

• Returns -> Return

• Revert this change

  -    OUTPUT:
  +    OUTPUT::

• In cyclic_mag_to_lightness(), use `\log(2)` (latex) instead of ``log(2)`` (code formatted). Also add an r before the first """.

• Similar changes as well:

  -    - ``r`` - a non-negative real number
  -    - ``arg`` - a real number
  +    - ``r`` -- a non-negative real number
  +    - ``arg`` -- a real number

• Add a one-line description to direct_complex_to_rgb() and cmap_complex_to_rgb().

• And similar changes

  -    - ``contoured`` -- boolean (default: False) - causes magnitude to be
  -      indicated through contour-like adjustments to lightness.
  +    - ``contoured`` -- boolean (default: ``False``); causes magnitude to be
  +      indicated through contour-like adjustments to lightness
   
  -    - ``tiled`` -- boolean (default: False) - causes magnitude and argument to
  -      be indicated through contour-like adjustments to lightness.
  +    - ``tiled`` -- boolean (default: ``False``); causes magnitude and argument to
  +      be indicated through contour-like adjustments to lightness

• In cmap_complex_to_rgb(), if you do r""", then you can change the \\ to \. Similar changes elsewhere too.

• And similar

  -for i from 0 <= i < imax:
  +for i in range(imax):

[mwageringel]

comment:5

Thank you for this addition. This will be nice to have in Sage. I have been playing around with this and have some remarks as well.

• The first thing I noticed is the aliasing effects. These are a bit unfortunate, but the way Sage interacts with Matplotlib makes it rather difficult to avoid them. At the time Sage creates the plot, the eventual size of the figure is not known, so upsampling of the image can not in general be avoided. I assume one might get fewer aliasing effects if it were possible to remap the colors of the upsampled image (so that the contours could be applied to the final image), but currently this seems very difficult to accomplish with Matplotlib. Therefore, I think the way you implemented it is the best we can do for now, even if it means that one needs a large number of plot points to obtain a smooth image.

• Renaming the function complex_to_rgb to direct_complex_to_rgb is not necessary. It is probably not used by anyone outside of Sage, but it is better to avoid breaking backward compatibility if we can. Also, functions that are meant for internal use should start with an underscore.

• Where possible we should use Numpy, unless there is a good reason not to. For example, the function clamp can be replaced by numpy.clip.

• I have tried to speed up the slow step in cmap_complex_to_rgb using Numpy. The following computation should be equivalent and faster, but please check this for correctness.

diff --git a/src/sage/plot/complex_plot.pyx b/src/sage/plot/complex_plot.pyx
index a3551584ce..d116887b83 100644
--- a/src/sage/plot/complex_plot.pyx
+++ b/src/sage/plot/complex_plot.pyx
@@ -387,14 +387,13 @@ def cmap_complex_to_rgb(z_values, cmap=None, contoured=False, tiled=False):
     normalized_colors = cmap((args + PI) / (2 * PI))
     normalized_colors = normalized_colors[:,:,:3]  # discard alpha channel
     lightdeltas = als[:,:,1]
-    lightdeltas = lightdeltas # normalize lightness adjustment
-    arg_d_s = np.dstack((normalized_colors, lightdeltas))

     if tiled or contoured:
+        arg_d_s = np.dstack((normalized_colors, lightdeltas))
         # add contours: convert to hls, adjust lightness, convert back
         rgbs = manual_contoured_cmap_complex_to_rgb(arg_d_s)
     else:
-        rgbs = manual_smooth_cmap_complex_to_rgb(arg_d_s)
+        rgbs = manual_smooth_cmap_complex_to_rgb(normalized_colors, lightdeltas)

     # Apply mask, making nan_indices white
     rgbs[nan_indices] = 1
@@ -451,7 +450,7 @@ cdef inline double clamp(double x):
     return x


-def manual_smooth_cmap_complex_to_rgb(rgb_d_s):
+def manual_smooth_cmap_complex_to_rgb(rgb, delta):
     """
     Returns an rgb array from given array of `(r, g, b, delta)`.

@@ -508,30 +507,10 @@ def manual_smooth_cmap_complex_to_rgb(rgb_d_s):
         array([[[0.75  , 0.8125, 0.875 ]]])
     """
     import numpy as np
-
-    cdef unsigned int i, j, imax, jmax
-    cdef double r, g, b, delta, h, l, s
-
-    imax = len(rgb_d_s)
-    jmax = len(rgb_d_s[0])
-
-    cdef cnumpy.ndarray[cnumpy.float_t, ndim=3, mode='c'] rgb = np.empty(dtype=float, shape=(imax, jmax, 3))
-
-    sig_on()
-    for i from 0 <= i < imax:
-        row = rgb_d_s[i]
-        for j from 0 <= j < jmax:
-            r, g, b, delta = row[j]
-            delta = 0.5 + delta/2.0
-            if delta < 0.5:
-                rgb[i][j][0] = clamp(2 * delta * r)
-                rgb[i][j][1] = clamp(2 * delta * g)
-                rgb[i][j][2] = clamp(2 * delta * b)
-            else:
-                rgb[i][j][0] = clamp(2*(delta - 1.0)*(1 - r) + 1)
-                rgb[i][j][1] = clamp(2*(delta - 1.0)*(1 - g) + 1)
-                rgb[i][j][2] = clamp(2*(delta - 1.0)*(1 - b) + 1)
-    sig_off()
+    delta = delta[:,:,np.newaxis]
+    delta_pos = delta > 0.0
+    rgb = (1.0 - np.abs(delta))*(rgb - delta_pos) + delta_pos
+    rgb = np.clip(rgb, 0.0, 1.0)
     return rgb

-• `rgb[i,j,k]` is faster than `rgb[i][j][k]` as the latter creates temporary arrays.

-• stacking the arrays was not necessary and cost time, so I kept `rgb` and `delta` separate.

Now, this is not much slower than the old implementation anymore:

sage: from sage.plot.complex_plot import *
....: arr = list(map(list, matrix.random(CDF, 800)))
....: %timeit res = direct_complex_to_rgb(arr)
....: %timeit res = cmap_complex_to_rgb(arr)
60.5 ms ± 1.34 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
144 ms ± 829 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

previously:

....: %timeit res = cmap_complex_to_rgb(arr)
1.13 s ± 9.25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

At this point, I think it makes sense to completely remove the old implementation and replace it by using the colormap hsv together with a sign change of the input data.

• For now I have focused on the colormap functionality. I have not looked at the contouring implementation yet, but noticed a problem with this example:

complex_plot(x, (-2,2), (-2,2), cmap='twilight', contoured=True)

The twilight colormap should be cyclic continuous, so the plot does not look right.

• For the documentation, I would also suggest to replace one of the examples that specified a colormap to use twilight instead because it is cyclic.

[sagetrac-git]

Changed commit from 57414a3 to d5b77df

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d5b77df

Incorporate suggestions from tscrim

[davidlowryduda]

comment:7

Replying to @tscrim:

Thank you for your suggestions and the documentation fixes. I've incorporated these changes.

• Is there an easy way to extend the contours to have different user input? I am thinking two parameters: one for choosing linear versus log scaling, and then another for the intervals. I am thinking about if someone is looking at where a function changes more slowly (a more trivial example, but something like f(x) = z^2).

It's not possible (as far as I can tell) to give the user the ability to choose the parameter for scaling. For example, I chose contours at 2^n for logarithmic scaling. Fortunately, logarithmic scaling is pretty forgiving since logarithms grow so slowly. For linear scaling, the difference between scaling at 1 and 10, say, is really dramatic. I didn't put in a default linear scaling option for this reason.

The reason it's nontrivial to include is that these are used heavily in the compiled portion of the cython. I think this is the same reason why there is no runtime option to change the exponent for mag_to_lightness in the current complex_plot.

Of course, I could be wrong.

• And similar

  -for i from 0 <= i < imax:
  +for i in range(imax):

Oh, thank you! I hadn't realized that cython changed to convert range-based loops into C-loops. (It seems this happened in 2018, and I'm behind the times). This is much nicer!

I'm currently examining the comments from @mwageringel, but haven't acted on those comments yet.

---

New commits:

d5b77df

Incorporate suggestions from tscrim

[tscrim]

comment:8

Replying to @davidlowryduda:

Replying to @tscrim:

Thank you for your suggestions and the documentation fixes. I've incorporated these changes.

• Is there an easy way to extend the contours to have different user input? I am thinking two parameters: one for choosing linear versus log scaling, and then another for the intervals. I am thinking about if someone is looking at where a function changes more slowly (a more trivial example, but something like f(x) = z^2).

It's not possible (as far as I can tell) to give the user the ability to choose the parameter for scaling. For example, I chose contours at 2^n for logarithmic scaling. Fortunately, logarithmic scaling is pretty forgiving since logarithms grow so slowly. For linear scaling, the difference between scaling at 1 and 10, say, is really dramatic. I didn't put in a default linear scaling option for this reason.

Someone might utilize that difference. With having an extra parameter that the user can control, this becomes a useful tool that I believe should be easy to implement.

The reason it's nontrivial to include is that these are used heavily in the compiled portion of the cython. I think this is the same reason why there is no runtime option to change the exponent for mag_to_lightness in the current complex_plot.

This doesn't make sense to me. You can just make your hardcoded values into input parameters, or have the lightness be a function (with the default being the hardcoded function).

[davidlowryduda]

comment:9

Replying to @tscrim:

This doesn't make sense to me. You can just make your hardcoded values into input parameters, or have the lightness be a function (with the default being the hardcoded function).

I don't know what I was thinking. Of course you're right. This is easy to implement. I'll do that too.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

b63fd38	Speed up cmap coloring with numpy
c768530	Use numpy whenever possible
06a921d	Parametrize contours

[sagetrac-git]

Changed commit from d5b77df to 06a921d

[davidlowryduda]

Example of aliasing artifacts for linear contours

[davidlowryduda]

comment:11

Attachment: result.png

I have now included several of the suggestions from @mwageringel. In particular, I have changed many of the implementations to now be in numpy, including writing numpy forms of HLS->RGB and RGB->HLS translations. I'm not particularly great at numpy, but nonetheless this is markedly faster than my previous version.

At this point, I think it makes sense to completely remove the old implementation and replace it by using the colormap hsv together with a sign change of the input data.

I almost did this completely. I changed the default complex_plot to use my new color implementation. There is a very small difference that can be seen by comparing the outputs of

sage: complex_plot(x^3, (-2, 2), (-2, 2))
sage: complex_plot(x^3, (-2, 2), (-2, 2), cmap=None) # old default

This seems close enough to me that I'd be happy to remove the old implementation entirely.

Is there an easy way to extend the contours to have different user input? I am thinking two parameters: one for choosing linear versus log scaling, and then another for the intervals. I am thinking about if someone is looking at where a function changes more slowly (a more trivial example, but something like f(x) = z^2).

These commits also include this implementation. It is now possible to choose between linear vs log scaling, to choose the intervals between contours, the number of phases in the tiling method, and the rate of change of magnitude when there are no contours. I note that linear scaling can often lead quickly to Moire patterns once the rate of change of the function gets to be large. For a concrete example, I've included the two outputs of the plots from

sage: complex_plot(x^3, (-5, 5), (-5, 5), cmap=None, contoured=True, contour_type='linear', contour_base=10, plot_points=800)
sage: complex_plot(x^3, (-5, 5), (-5, 5), cmap=None, contoured=True, contour_type='linear', contour_base=10, plot_points=800)

I'll add one such example and a warning to the documentation.

It is also now possible to adjust the exponent controlling how quickly the magnitude tends to white or black. This can be seen for example in

sage: complex_plot(x^3, (-5, 5), (-5, 5), cmap=None, plot_points=300)
sage: complex_plot(x^3, (-5, 5), (-5, 5), cmap=None, plot_points=300, dark_rate=1.0)

[davidlowryduda]

Attachment: twilight.png

twilight example

[davidlowryduda]

comment:12

Replying to @mwageringel:

• For now I have focused on the colormap functionality. I have not looked at the contouring implementation yet, but noticed a problem with this example:

complex_plot(x, (-2,2), (-2,2), cmap='twilight', contoured=True)

The twilight colormap should be cyclic continuous, so the plot does not look right.

The output looks about as I expect. For concreteness, I've included a version (with 300 plot points) here:

There are two interactions governing the colors. First, there are contours affecting magnitudes. Then there is the colormap. The phase is mapped to the color, so we should expect on any circle centered at the origin to go through a cyclic continuous colormap (but each circle might be adjusted in lightness by the contouring). In a commit that will immediately after I write this, it's also possible to choose contour visibility (by adjusting dark_rate). With twilight, I find that setting dark_rate = 0.2 works better. Setting dark_rate = 0.0 now gives a pure phase plot.

• For the documentation, I would also suggest to replace one of the examples that specified a colormap to use twilight instead because it is cyclic.

I think this is a good idea! I think I'm going to add many examples.

[sagetrac-git]

Changed commit from 06a921d to 104b339

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

104b339

Also allow changing contour visibility

[tscrim]

comment:14

Great, thank you. Looks very pretty. I just have a few trivial doc comments:

And similar:

     .. SEEALSO::
 
-        :func:`sage.plot.complex_plot.mag_to_lightness`
-        :func:`sage.plot.complex_plot.mag_and_arg_to_lightness`
-        :func:`sage.plot.complex_plot.cyclic_linear_mag_to_lightness`
+        - :func:`sage.plot.complex_plot.mag_to_lightness`
+        - :func:`sage.plot.complex_plot.mag_and_arg_to_lightness`
+        - :func:`sage.plot.complex_plot.cyclic_linear_mag_to_lightness`

For INPUT: blocks, please format as:

- ``arg1`` -- positive integer; a description
- ``arg2`` -- boolean (default: ``False``); some description

Instead of using ``... \times 3``, use `n \times 3`.

[mwageringel]

comment:15

Replying to @davidlowryduda:

I almost did this completely. I changed the default complex_plot to use my new color implementation. There is a very small difference that can be seen by comparing the outputs of

sage: complex_plot(x^3, (-2, 2), (-2, 2))
sage: complex_plot(x^3, (-2, 2), (-2, 2), cmap=None) # old default

This seems close enough to me that I'd be happy to remove the old implementation entirely.

The difference is very small, but I cannot fully explain where the difference comes from. On closer inspection, I also noticed that the matplotlib hsv color map is discretized (which I did not expect), so one can see a difference when the phase varies very little. For example:

graphics_array([complex_plot(x^(1/10)*e^(2*pi*I/10), (-2, 2), (-2, 2)), complex_plot(x^(1/10)*e^(2*pi*I/10), (-2, 2), (-2, 2), cmap=None)])

A way around this would be creating our own hsv color map with more than 256 sampling points, but maybe that goes too far for this ticket. For now, should we just keep the old implementation?

Replying to @davidlowryduda:

The output looks about as I expect. For concreteness, I've included a version (with 300 plot points) here:

Yes, I now obtain the same result. When I reported the problem, there was a discontinuity in the color spectrum at about 10 degrees, but it seems to be fixed now.

I still need to take a closer look at the implementation.

[tscrim]

comment:16

I wouldn’t worry too much about minor aliasing differences when there is slight variation (some monitors might not be able to resolve these differences either). That being said, I am not opposed to keeping the current implementation since it does provide slightly better performance and it is already there and working.

[sagetrac-git]

Changed commit from 104b339 to f205f94

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d864930

Remove unnecessary allocation

[mwageringel]

Changed commit from d864930 to 378755d

[mwageringel]

Changed branch from u/DavidLowry/colormaps_and_contours_for_complex_plot to u/gh-mwageringel/33416

[mwageringel]

comment:27

Thank you for all the changes. I have fixed one more typo in the documentation and have also run the tests locally.

---

New commits:

378755d

33416: fix another typo

[vbraun]

comment:28

File "src/sage/plot/complex_plot.pyx", line 256, in sage.plot.complex_plot.?
Failed example:
    complex_to_rgb([[0, 1 + 1j, -3 - 4j]], tiled=True, contour_base=5, nphases=15)
Expected:
    array([[[1.        , 0.        , 0.        ],
            [0.87741543, 0.65806157, 0.        ],
            [0.        , 0.11124213, 0.97156143]]])
Got:
    array([[[1.        , 0.        , 0.        ],
            [0.87741543, 0.65806157, 0.        ],
            [0.        , 0.08261755, 0.72156143]]])

[mwageringel]

comment:29

Thank you, Volker.

@davidlowry: To explain the problem a bit further, the problem with this test is that it is not numerically stable. I cannot replicate the problem directly, but adding a small perturbation to one of the values demonstrates it:

sage: from sage.plot.complex_plot import complex_to_rgb
sage: complex_to_rgb([[0, 1 + 1j, -3 - 4j + 1e-9]], tiled=True, contour_base=5, nphases=15)
[[[1.         0.         0.        ]
  [0.87741543 0.65806157 0.        ]
  [0.         0.08261755 0.72156143]]]

Probably, Volker observed the issue on a 32-bit machine, which can lead to numerical differences like this.

Is this a problem with the test or with the implementation, i.e. should the implementation of complex_to_rgb yield numerically stable results? I guess the tiling can never be completely stable, but seeing this problem with only three values is a bit unexpected. If it is just the test, maybe you could replace it by a different one.

[davidlowryduda]

comment:30

Oh, thank you Markus. I couldn't replicate this and didn't understand. But you're right!

Looking at this test case now, I see that this is a poor choice of parameter. I think I chose them to be "random small numbers", but the law of small numbers bites here.

The reason this test acting like this is because a rgb value is being associated to a point lying exactly on a contour boundary. Here, contours are associated on powers of contour_base=5, and the point has absolute value exactly 5. The actual evaluation of contour boundaries occurs after taking logs, and thus it's reasonable to expect a certain amount of instability.

Away from the contours, these tests should be fine and pretty stable. I'll replace this test with a different one. I'll also make sure that I didn't make this same error elsewhere.

[davidlowryduda]

Changed branch from u/gh-mwageringel/33416 to u/DavidLowry/33416

[davidlowryduda]

Changed commit from 378755d to 9338f8c

[davidlowryduda]

New commits:

9338f8c

Choose tests not directly on contours

[tscrim]

comment:33

I actually wonder how much we have to worry about numerical noise. Perhaps we can modify these doctests to not have explicit values? Maybe just the first two or three sigfigs with ...?

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

7654433

Add in small error tolerance to numerical tests

[sagetrac-git]

Changed commit from 9338f8c to 7654433

[davidlowryduda]

comment:35

@tscrim That makes sense to me. I've incorporated a pretty general absolute tolerance into the tests.

[tscrim]

comment:36

Thanks. I changed some numbers to verify that the tolerance works. Two last little changes: The doctests on lines 1104 and 1114 should be marked as # long time since they take >3s on my (new fast) desktop. Once that is done, you can (re)set a positive review on my behalf.

[slel]

comment:37

Minor suggestion:

 .. [LD2021] David Lowry-Duda.
             *Visualizing modular forms*.
-            Arithmetic Geometry, Number Theory, and Computation, Simons
-            Symposia, Springer, 2021.
+            Arithmetic Geometry, Number Theory, and Computation.
+            Simons Symposia, Springer, 2021.
             :arxiv:`2002.05234`

[sagetrac-git]

Changed commit from 7654433 to e1b7ca9

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

e1b7ca9

Mark long tests and adjust formatting of reference

[davidlowryduda]

comment:39

Thanks @tscrim and @slel. I've made these changes. And at @tscrim's suggestion, I'm also marking this positively reviewed.

[vbraun]

Changed branch from u/DavidLowry/33416 to e1b7ca9