Plotting sets using natural mathematical set expressions

[tapyu]

In convex optimization, it is quite common to come across a n-dimensional set formally written like that, for example,

C = \left\{ (x_1, x_2, \cdots, x_n, t) \in \mathbb{R}^{n+1} \mid \mathbf{x} \in \mathbb{R}^{n}, \lVert\mathbf{x}\rVert_{p} \leq t \right\} \subseteq \mathbb{R}^{n+1}

It would be amazing to have a Julia package that we give the set as it is written mathematically and get plot of the set (in a 3D space, at most, of course…). The output woud be something like that (this is not the resulting plot of the previous set, btw)

Is there a package that fill this gap? If the answer is no, what is the most convenient way to get the plot of such set in Julia?

[jd-foster]

You may be interested in the combination of Polyhedra.jl

• Optimization · Polyhedra
• Plot · Polyhedra

with Makie or MeshCat.jl.

[John_Gibson]

That’s a nice idea. Maybe aim for software that would give some kind of visual representation for some subset of mathematical definitions. It would be so valuable for students. It must be doable at some level (and probably errors) with current AI.

[tapyu]

Thanks to @jd-foster , I could gather required packages to plot sets. It is not that hard, but the user must be acquaint with different Julia’s packages, which makes it more confusing for beginners… I started off with a simpler set:

\left\{\mathbf{x} \in \mathbb{R}^3 \mid \mathbf{0} \preceq \mathbf{x} \preceq \mathbf{1} \right\}

My Julia code is

import JuMP, CDDLib, Polyhedra, GLMakie, Makie

n = 3
m = JuMP.Model()
JuMP.@variable(m, 0 ≤ x[1:n] ≤ 1) # put your mathematical set here!
poly = Polyhedra.polyhedron(m, CDDLib.Library(:exact))
m = Polyhedra.Mesh(poly)
Makie.mesh(m, color=:blue)

It yields this beautiful (but simple) 3D interactive plot

I am not sure how these packages relate each other, but as far as I noticed:

• GLMakie is the render backend of Makie
• Polyhedra is responsible for creating polyhedra on Julia
• JuMP seems to be a optimization package, which is weird as we are not optimizing anything here… Furthermore, it seems to have an overlap with the Convex.jl, which I personally prefer…

I have no idea what the purpose of CDDLib is…

[tapyu]

However, for a minimally more complicated set, e.g.,

\left\{ (x,y) \in \mathbb{R}_{++}^{2} \mid x/y \leq 1 \right\}

Julia code fails in plotting it trivially…

import JuMP, CDDLib, Polyhedra, Plots
m = JuMP.Model()
JuMP.@variable(m, 0 ≤ x[1:2]) # it is not possible to define opened intervals
JuMP.@constraint(m, x[1] ≤ x[2])
poly = Polyhedra.polyhedron(m, CDDLib.Library(:exact))
Plots.plot(poly)

It prompts the following warning

ERROR: Rays not supported yet in the 2D plotting recipe.

If there is a way to plot nonbounded set, it is not so simple as I predicted. We need to truncate the set in order to plot it…

import JuMP, CDDLib, Polyhedra, Plots
m = JuMP.Model()
JuMP.@variable(m, 0 ≤ x[1:2] ≤ 3) # it is not possible to define opened intervals, the value 3 was arbitrary
JuMP.@constraint(m, x[1] ≤ x[2])
poly = Polyhedra.polyhedron(m, CDDLib.Library(:exact))
Plots.plot(poly)

[tapyu]

Ok, this is harder when the constraint happens to be nonlinear… For example, I wrote this set

\left\{ (x,y) \in \mathbb{R}_{+}^{2} \mid xy \leq 1 \right\}

in Julia:

import JuMP, CDDLib, Polyhedra, Plots
m = JuMP.Model()
JuMP.@variable(m, 0.001 ≤ x[1:2] ≤ 3) # it is not possible to define opened intervals, the value 3 was arbitrary
JuMP.@NLconstraint(m, x[1]*x[2] ≤ 1)
poly = Polyhedra.polyhedron(m, CDDLib.Library(:exact))
Plots.plot(poly)

But it produced this nonsense set

I do not know where is the error. However, I suppose that, since this set is not a polyhedron anymore, the command Polyhedra.polyhedron(m, CDDLib.Library(:exact)) produces a wrong set. how to solve it?

[odow]

The name of Polyhedra.jl gives the game away; it does not support nonlinear functions.

As you’ve found out, the current setup is quite limited. It supports only finite volume polyhedra in 2 or 3 dimensions. If you want something closer to your set syntax, you’ll likely need to write this yourself. I don’t know if there is anything off-the-shelf that can do what you want.

To answer some of your other questions:

• JuMP is used as the modeling layer. It doesn’t always have to be used for optimization.
• CDDlib is used to generate the list of vertices of the polyhedron from the constraint representation in JuMP

[rafael.guerra]

Just in case, to plot this one:

\left\{ (x,y) \in \mathbb{R}^{2} \mid xy \leq 1 \right\}

You could use the ImplicitEquations.jl package:

using ImplicitEquations, Plots
f(x,y) = x * y
plot(f ≦ 1, fc=:blues, widen=false)         # \leqq[tab]

[rafael.guerra]

And this one:

\left\{ (x,y) \in \mathbb{R}_{++}^{2} \mid x/y \leq 1 \right\}

It could possibly be plotted as follows:

using ImplicitEquations, Plots
f(x,y) = x / y
g(x,y) = x
h(x,y) = y
plot((f ≦ 1) & (g ≫ 0) & (h ≫ 0), fc=:reds, widen=false)     # \leqq[tab] ,  \gg[tab]

[tapyu]

@rafael.guerra ImplicitEquations seems much more powerful than Polyhedra since the former allows more nonlinear sets (the necessity of requiring \gg instead > is kind of bewildering, but I will study more this pkg). Anyway, thank you for the sharing.

[rafael.guerra]

@tapyu

the necessity of requiring ≫ instead > is kind of bewildering

According to the author of the package, this is to avoid type pyracy.