diff --git a/SYMBOLS_MANIFEST.txt b/SYMBOLS_MANIFEST.txt
index cc534ee1c..9a2f9c552 100644
--- a/SYMBOLS_MANIFEST.txt
+++ b/SYMBOLS_MANIFEST.txt
@@ -265,6 +265,7 @@ System`Context
System`Contexts
System`Continue
System`ContinuedFraction
+System`ContourPlot
System`Convert`B64Dump`B64Decode
System`Convert`B64Dump`B64Encode
System`ConvertersDump`$ExtensionMappings
diff --git a/mathics/builtin/box/graphics3d.py b/mathics/builtin/box/graphics3d.py
index dc112128a..e6e3d8eea 100644
--- a/mathics/builtin/box/graphics3d.py
+++ b/mathics/builtin/box/graphics3d.py
@@ -250,6 +250,17 @@ def _prepare_elements(self, elements, options, max_width=None):
self.background_color = elements.background_color
def calc_dimensions(final_pass=True):
+ # TODO: the code below is broken in any other case but Automatic
+ # because it calls elements.translate which is not implemented.
+ # Plots may pass specific plot ranges, triggering this deficiency
+ # and causing tests to fail The following line avoids this,
+ # and it should not change the behavior of any case which did
+ # previously fail with an exception.
+ #
+ # This code should be DRYed (together with the very similar code
+ # for the 2d case), and the missing .translate method added.
+ plot_range = ["System`Automatic"] * 3
+
if "System`Automatic" in plot_range:
xmin, xmax, ymin, ymax, zmin, zmax = elements.extent()
else:
@@ -291,7 +302,7 @@ def calc_dimensions(final_pass=True):
elif zmin == zmax:
zmin -= 1
zmax += 1
- elif isinstance(plot_range[1], list) and len(plot_range[1]) == 2:
+ elif isinstance(plot_range[1], list) and len(plot_range[2]) == 2:
zmin, zmax = list(map(float, plot_range[2]))
zmin = elements.translate((0, 0, zmin))[2]
zmax = elements.translate((0, 0, zmax))[2]
diff --git a/mathics/builtin/drawing/plot.py b/mathics/builtin/drawing/plot.py
index 1cebf9c4c..cf59c329d 100644
--- a/mathics/builtin/drawing/plot.py
+++ b/mathics/builtin/drawing/plot.py
@@ -37,6 +37,7 @@
SymbolLine,
SymbolLog10,
SymbolNone,
+ SymbolPlotRange,
SymbolRGBColor,
SymbolSequence,
SymbolStyle,
@@ -67,6 +68,7 @@
from mathics.eval.drawing.plot3d_vectorized import (
eval_ComplexPlot,
eval_ComplexPlot3D,
+ eval_ContourPlot,
eval_DensityPlot,
eval_Plot3D,
)
@@ -74,6 +76,7 @@
from mathics.eval.drawing.plot3d import (
eval_ComplexPlot,
eval_ComplexPlot3D,
+ eval_ContourPlot,
eval_DensityPlot,
eval_Plot3D,
)
@@ -467,7 +470,7 @@ def error(self, what, *args, **kwargs):
def __init__(self, expr, range_exprs, options, evaluation):
self.evaluation = evaluation
- # plot ranges
+ # plot ranges of the form {x,xmin,xmax} etc.
self.ranges = []
for range_expr in range_exprs:
if not range_expr.has_form("List", 3):
@@ -488,6 +491,19 @@ def __init__(self, expr, range_exprs, options, evaluation):
self.error(expr, "invrange", range_expr)
self.ranges.append(range)
+ # Contours option
+ contours = expr.get_option(options, "Contours", evaluation)
+ if contours is not None:
+ c = contours.to_python()
+ if not (
+ c == "System`Automatic"
+ or isinstance(c, int)
+ or isinstance(c, tuple)
+ and all(isinstance(cc, (int, float)) for cc in c)
+ ):
+ self.error(expr, "invcontour", contours)
+ self.contours = c
+
# Mesh option
mesh = expr.get_option(options, "Mesh", evaluation)
if mesh not in (SymbolNone, SymbolFull, SymbolAll):
@@ -579,6 +595,9 @@ class _Plot3D(Builtin):
"Plot range `1` must be of the form {variable, min, max}, "
"where max > min."
),
+ "invcontour": (
+ "Contours option must be Automatic, an integer, or a list of numbers."
+ ),
}
# Plot3D, ComplexPlot3D
@@ -633,6 +652,39 @@ def eval(
graphics = self.eval_function(plot_options, evaluation)
if not graphics:
return
+
+ # Expand PlotRange option using the {x,xmin,xmax} etc. range specifications
+ # Pythonize it, so Symbol becomes str, numeric becomes int or float
+ plot_range = self.get_option(options, str(SymbolPlotRange), evaluation)
+ plot_range = plot_range.to_python()
+ dim = 3 if self.graphics_class is Graphics3D else 2
+ if isinstance(plot_range, str):
+ # PlotRange -> Automatic becomes PlotRange -> {Automatic, ...}
+ plot_range = [str(SymbolAutomatic)] * dim
+ if isinstance(plot_range, (int, float)):
+ # PlotRange -> s becomes PlotRange -> {Automatic,...,{-s,s}}
+ pr = plot_range
+ plot_range = [str(SymbolAutomatic)] * dim
+ plot_range[-1] = [-pr, pr]
+ elif isinstance(plot_range, (list, tuple)) and isinstance(
+ plot_range[0], (int, float)
+ ):
+ # PlotRange -> {s0,s1} becomes PlotRange -> {Automatic,...,{s0,s1}}
+ pr = plot_range
+ plot_range = [str(SymbolAutomatic)] * dim
+ plot_range[-1] = pr
+
+ # now we have a list of length dim
+ # handle Automatic ~ {xmin,xmax} etc.
+ for i, (pr, r) in enumerate(zip(plot_range, plot_options.ranges)):
+ # TODO: this treats Automatic and Full as the same, which isn't quite right
+ if isinstance(pr, str) and not isinstance(r[1], complex):
+ plot_range[i] = r[1:] # extract {xmin,xmax} from {x,xmin,xmax}
+
+ # unpythonize and update PlotRange option
+ options[str(SymbolPlotRange)] = to_mathics_list(*plot_range)
+
+ # generate the Graphics[3D] result
graphics_expr = graphics.generate(
options_to_rules(options, self.graphics_class.options)
)
@@ -809,6 +861,32 @@ class ComplexPlot(_Plot3D):
graphics_class = Graphics
+class ContourPlot(_Plot3D):
+ """
+ :WMA link: https://reference.wolfram.com/language/ref/ContourPlot.html
+
+ - 'Contour'[$f$, {$x$, $x_{min}$, $x_{max}$}, {$y$, $y_{min}$, $y_{max}$}]
+
- creates a two-dimensional contour plot ofh $f$ over the region
+ $x$ ranging from $x_{min}$ to $x_{max}$ and $y$ ranging from $y_{min}$ to $y_{max}$.
+
+ See :Drawing Option and Option Values:
+ /doc/reference-of-built-in-symbols/graphics-and-drawing/drawing-options-and-option-values
+ for a list of Plot options.
+
+
+ """
+
+ requires = ["skimage"]
+ summary_text = "creates a contour plot"
+ expected_args = 3
+ options = _Plot3D.options2d | {"Contours": "Automatic"}
+ # TODO: other options?
+
+ many_functions = True
+ eval_function = staticmethod(eval_ContourPlot)
+ graphics_class = Graphics
+
+
class DensityPlot(_Plot3D):
"""
:WMA link: https://reference.wolfram.com/language/ref/DensityPlot.html
diff --git a/mathics/core/builtin.py b/mathics/core/builtin.py
index f91e4ea76..f1784391e 100644
--- a/mathics/core/builtin.py
+++ b/mathics/core/builtin.py
@@ -852,12 +852,14 @@ class UnavailableFunction:
def __init__(self, builtin):
self.name = builtin.get_name()
+ self.requires = builtin.requires
def __call__(self, **kwargs):
kwargs["evaluation"].message(
"General",
"pyimport", # see messages.py for error message definition
strip_context(self.name),
+ ", ".join(self.requires),
)
diff --git a/mathics/core/systemsymbols.py b/mathics/core/systemsymbols.py
index ae554c96d..69ace5548 100644
--- a/mathics/core/systemsymbols.py
+++ b/mathics/core/systemsymbols.py
@@ -228,6 +228,7 @@
SymbolPiecewise = Symbol("System`Piecewise")
SymbolPlot = Symbol("System`Plot")
SymbolPlotLabel = Symbol("System`PlotLabel")
+SymbolPlotRange = Symbol("System`PlotRange")
SymbolPlotRangeClipping = Symbol("System`PlotRangeClipping")
SymbolPlotRegion = Symbol("System`PlotRegion")
SymbolPlus = Symbol("System`Plus")
diff --git a/mathics/eval/drawing/plot3d.py b/mathics/eval/drawing/plot3d.py
index d297247c1..9edc31acb 100644
--- a/mathics/eval/drawing/plot3d.py
+++ b/mathics/eval/drawing/plot3d.py
@@ -506,3 +506,10 @@ def eval_ComplexPlot(
evaluation: Evaluation,
):
return None
+
+
+def eval_ContourPlot(
+ plot_options,
+ evaluation: Evaluation,
+):
+ return None
diff --git a/mathics/eval/drawing/plot3d_vectorized.py b/mathics/eval/drawing/plot3d_vectorized.py
index d08c50f3d..1d7e94674 100644
--- a/mathics/eval/drawing/plot3d_vectorized.py
+++ b/mathics/eval/drawing/plot3d_vectorized.py
@@ -9,8 +9,14 @@
from mathics.builtin.colors.color_internals import convert_color
from mathics.core.evaluation import Evaluation
+from mathics.core.expression import Expression
from mathics.core.symbols import strip_context
-from mathics.core.systemsymbols import SymbolNone, SymbolRGBColor
+from mathics.core.systemsymbols import (
+ SymbolEqual,
+ SymbolNone,
+ SymbolRGBColor,
+ SymbolSubtract,
+)
from mathics.timing import Timer
from .plot_compile import plot_compile
@@ -124,28 +130,64 @@ def compute_over_grid(nx, ny):
return graphics
+# color-blind friendly palette from https://davidmathlogic.com/colorblind
+# palette3 is for 3d plots, i.e. surfaces
+palette3 = [
+ (255, 176, 0), # orange
+ (100, 143, 255), # blue
+ (220, 38, 127), # red
+ (50, 150, 140), # green
+ (120, 94, 240), # purple
+ (254, 97, 0), # dark orange
+ (0, 114, 178), # dark blue
+]
+
+
+# palette 2 is for 2d plots, i.e. lines
+# same colors as palette3 but in a little different order
+palette2 = [
+ (100, 143, 255), # blue
+ (255, 176, 0), # orange
+ (50, 150, 140), # green
+ (220, 38, 127), # red
+ (120, 94, 240), # purple
+ (254, 97, 0), # dark orange
+ (0, 114, 178), # dark blue
+]
+
+
+def palette_color_directive(palette, i):
+ """returns a directive in a form suitable for graphics.add_directives"""
+ """ for setting the color of an entire entity such as a line or surface """
+ rgb = palette[i % len(palette)]
+ rgb = [c / 255.0 for c in rgb]
+ return [SymbolRGBColor, *rgb]
+
+
+@Timer("density_colors")
+def density_colors(zs):
+ """default color palette for DensityPlot and ContourPlot (f(x) form)"""
+ z_min, z_max = min(zs), max(zs)
+ zs = zs[:, np.newaxis] # allow broadcasting
+ # c_min, c_max = [0.3, 0.00, 0.3], [1.0, 0.95, 0.8]
+ c_min, c_max = [0.5, 0, 0.1], [1.0, 0.9, 0.5]
+ c_min, c_max = (
+ np.full((len(zs), 3), c_min),
+ np.full((len(zs), 3), c_max),
+ )
+ colors = ((zs - z_min) * c_max + (z_max - zs) * c_min) / (z_max - z_min)
+ return colors
+
+
@Timer("eval_Plot3D")
def eval_Plot3D(
plot_options,
evaluation: Evaluation,
):
def emit(graphics, i, xyzs, quads):
- # color-blind friendly palette from https://davidmathlogic.com/colorblind
- palette = [
- (255, 176, 0), # orange
- (100, 143, 255), # blue
- (220, 38, 127), # red
- (50, 150, 140), # green
- (120, 94, 240), # purple
- (254, 97, 0), # dark orange
- (0, 114, 178), # dark blue
- ]
-
# choose a color
- rgb = palette[i % len(palette)]
- rgb = [c / 255.0 for c in rgb]
- # graphics.add_color(SymbolRGBColor, rgb)
- graphics.add_directives([SymbolRGBColor, *rgb])
+ color_directive = palette_color_directive(palette3, i)
+ graphics.add_directives(color_directive)
# add a GraphicsComplex displaying a surface for this function
graphics.add_complex(xyzs, lines=None, polys=quads)
@@ -161,16 +203,7 @@ def eval_DensityPlot(
def emit(graphics, i, xyzs, quads):
# Fixed palette for now
# TODO: accept color options
- with Timer("compute colors"):
- zs = xyzs[:, 2]
- z_min, z_max = min(zs), max(zs)
- zs = zs[:, np.newaxis] # allow broadcasting
- c_min, c_max = [0.5, 0, 0.1], [1.0, 0.9, 0.5]
- c_min, c_max = (
- np.full((len(zs), 3), c_min),
- np.full((len(zs), 3), c_max),
- )
- colors = ((zs - z_min) * c_max + (z_max - zs) * c_min) / (z_max - z_min)
+ colors = density_colors(xyzs[:, 2])
# flatten the points and add the quads
graphics.add_complex(xyzs[:, 0:2], lines=None, polys=quads, colors=colors)
@@ -178,6 +211,85 @@ def emit(graphics, i, xyzs, quads):
return make_plot(plot_options, evaluation, dim=2, is_complex=False, emit=emit)
+@Timer("eval_ContourPlot")
+def eval_ContourPlot(
+ plot_options,
+ evaluation: Evaluation,
+):
+ import skimage.measure
+
+ # whether to show a background similar to density plot, except quantized
+ background = len(plot_options.functions) == 1
+ contour_levels = plot_options.contours
+
+ # convert fs of the form a==b to a-b, inplicit contour level 0
+ plot_options.functions = list(plot_options.functions) # so we can modify it
+ for i, f in enumerate(plot_options.functions):
+ if f.head == SymbolEqual:
+ f = Expression(SymbolSubtract, *f.elements[0:2])
+ plot_options.functions[i] = f
+ contour_levels = [0]
+ background = False
+
+ def emit(graphics, i, xyzs, quads):
+ # set line color
+ if background:
+ # showing a background, so just black lines
+ color_directive = [SymbolRGBColor, 0, 0, 0]
+ else:
+ # no background - each curve gets its own color
+ color_directive = palette_color_directive(palette2, i)
+ graphics.add_directives(color_directive)
+
+ # get data
+ nx, ny = plot_options.plotpoints
+ _, xmin, xmax = plot_options.ranges[0]
+ _, ymin, ymax = plot_options.ranges[1]
+ zs = xyzs[:, 2] # this is a linear list matching with quads
+
+ # process contour_levels
+ levels = contour_levels
+ zmin, zmax = np.min(zs), np.max(zs)
+ if isinstance(levels, str):
+ # TODO: need to pick "nice" number so levels have few digits
+ # an odd number ensures there is a contour at 0 if range is balanced
+ levels = 9
+ if isinstance(levels, int):
+ # computed contour levels have equal distance between them,
+ # and half that between first/last contours and zmin/zmax
+ dz = (zmax - zmin) / levels
+ levels = zmin + np.arange(levels) * dz + dz / 2
+
+ # one contour line per contour level
+ for level in levels:
+ # find contours and add lines
+ with Timer("contours"):
+ zgrid = zs.reshape((nx, ny)) # find_contours needs it as an array
+ contours = skimage.measure.find_contours(zgrid, level)
+
+ # add lines
+ for segment in contours:
+ segment[:, 0] = segment[:, 0] * ((xmax - xmin) / nx) + xmin
+ segment[:, 1] = segment[:, 1] * ((ymax - ymin) / ny) + ymin
+ graphics.add_linexyzs(segment)
+
+ # background is solid colors between contour lines
+ if background:
+ with Timer("contour background"):
+ # use masks and fancy indexing to assign (lo+hi)/2 to all zs between lo and hi
+ zs[zs <= levels[0]] = zmin
+ for lo, hi in zip(levels[:-1], levels[1:]):
+ zs[(lo < zs) & (zs <= hi)] = (lo + hi) / 2
+ zs[levels[-1] < zs] = zmax
+ colors = density_colors(zs) # same colors as density plot
+ graphics.add_complex(
+ xyzs[:, 0:2], lines=None, polys=quads, colors=colors
+ )
+
+ # plot_options.plotpoints = [n * 10 for n in plot_options.plotpoints]
+ return make_plot(plot_options, evaluation, dim=2, is_complex=False, emit=emit)
+
+
@Timer("complex colors")
def complex_colors(zs, s=None):
# hue depends on phase