Plots for unbounded sets

[schillic]

We have the undounded set types HalfSpace, Hyperplane, Line, HPolyhedron, and Universe. ~I think we should print an error message when trying to plot such sets. Currently the plotting crashes when trying to overapproximate the sets.~ EDIT: We show an error now.

I also thought about adding a common interface/abstract supertype like UnboundedSet. However, I do not see any other application at the moment. And it would have the dual problem to what we had with ConvexSet, namely that set operations involving these set types would not be UnboundedSets.

[mforets]

it would be nice to plot such sets in 2D. but one needs to pass xmin/xmax values for the plots.

[schillic]

Good point!

[schillic]

This function plots a half-space by first obtaining the x- and y-limits and then computing the two y-coordinates at the left and right border. (One caveat: If the half-space has coefficient 0 for dimension y, the code does not work.) Finally it plots the given line segment and fills the area on the correct side of the line segment.

function plot_halfspace(H)
    pa = p.series_list[1].plotattributes
    x, y = pa[:x], pa[:y]
    y = [(H.a[1] * z - H.b) / (-H.a[2]) for z in x]
    fillrange = minimum(y)
    plot(x, y, fillrange=fillrange)
end

Note that if the axis ranges are extended later, the filling is not extended, so in the example you get a triangle.

[schillic]

See #1422 in Plots.jl and the referenced issues. They have hspan and vspan, which could be a start.

[ueliwechsler]

I think, something like this would work

function bounded_HPolytope_from_unbounded(P::HPolyhedron{N}, limit=1000) where {N}
    A, a = tosimplehrep(P)
    n = LazySets.dim(P)
    B = [Matrix{N}(I,n,n); -Matrix{N}(I,n,n)]
    b = limit*ones(N,2*n)
    return HPolytope([A;B], [a;b])
end

and add a condition here to create a HPolytope from the unbounded HPolyhedron https://github.com/JuliaReach/LazySets.jl/blob/c4e52db11b6c90543026a34d08571e7ae537bdf5/src/AbstractPolyhedron_functions.jl#L404

But where would I change the limits of the resulting plot? Can I do it directly in the plot_recipe or do I need to do this on a higher level?

[schillic]

But where would I change the limits of the resulting plot? Can I do it directly in the plot_recipe or do I need to do this on a higher level?

There are the lims/xlims/ylims keywords, and in https://github.com/JuliaReach/LazySets.jl/issues/576#issuecomment-477794073 I think I gave a way to obtain the current bounds.

One way to implement all this is probably to first obtain the plot bounds before plotting unbounded sets, then plotting the bounded part you computed, and then enforcing the previous bounds again. (This would not work if you plot the unbounded set as the first set).

[mforets]

i haven't checked, maybe Polyhedra already supports this feature?

[ueliwechsler]

As far as I know (Documentation and trials) it is not supported.

[ueliwechsler]

I have some first draft, but have not figured out how to use it with the plot_recipes. It is still hacky, but if think not to bad ;) The maximum of the existing shapes in the plot defines the region, in which the unbounded Polyhedra will be "bounded". If there is no plot, it starts with ylims=xlims=(-10,10). As a result, if the limits of the plot change after plotting, it has to be redrawn. The gif at the end shows, how it could look in practice.

using Plots
using LazySets
function bounded_approximation(P::HPolyhedron{N}, emax, emin) where {N}
    A, a = tosimplehrep(P)
    n = LazySets.dim(P)
    B = [Matrix{N}(I,n,n); -Matrix{N}(I,n,n)]
    # b = [x_max, y_max; -x_min, -y_min]
    b = [emax...; -emin...]
    return HPolytope([A;B], [a;b])
end

function plot_extrem_points()
    # gives you the extrem point and not predefined limits
    p = plot!()
    xaxis = p[1][:xaxis]
    yaxis = p[1][:yaxis]
    emax = [xaxis[:extrema].emax; yaxis[:extrema].emax]
    emin = [xaxis[:extrema].emin; yaxis[:extrema].emin]
    if Inf ∈ [emax...; emin...]
        # rand(2) decrease the possibility of hitting a singular point
        emax = 10*ones(2) + rand(2)
        emin = -10*ones(2) + + rand(2)
    end
    return emax, emin
end

function plot_polyhedron!(P::HPolyhedron{N}; kwargs...) where {N}
    emax, emin = plot_extrem_points()
    # TODO: it could be that bounded_approximation returns a emptyset!!
    P_bounded = bounded_approximation(P, emax, emin)
    while isempty(P_bounded)
        # recompute P_bounded wiht new emax, emin
        center = [(emax[i]+emin[i])/2 for i=1:2]
        radius_high = emax - center
        radius_low = emin - center
        emax, emin = emax + radius_high, emin + radius_low
        P_bounded = bounded_approximation(P, emax, emin)
    end
    plot!(P_bounded; kwargs...)
end

LazySets.plot_polyhedron!(HPolyhedron([0.5 0.5], -[9.0]))

[schillic]

Nice! Given that this is better than nothing, we can go for it.

Problems

1. You have to statically guess the plot area. [-10, 10] is probably a good range in most cases, but it is easy to think of a counterexample.
2. If I understand correctly, the user needs to actively redraw, which is annoying.
3. The previous versions of the unbounded sets are not removed, so with transparency enabled you still see those.

For problem 1 you can check whether the set is visible at all (i.e., intersects with the current plot-range box). I think this is the problem you mentioned above:

TODO: it could be that bounded_approximation returns a emptyset!!

Unfortunately, there is no easy way to find out in which direction the plot needs to be expanded. Maybe you can sample a point from the set and only expand in the corresponding x/y direction, but this is only a heuristics. Of course in general you not only want to see a part of the set but "the interesting part". This is more difficult (maybe even only heuristically achievable).

For problems 2 and 3 a user can put all the sets into a vector and plot them all at once. First I would plot the bounded sets to get an initial range. Then I would compute the "bounding" of the unbounded sets based on this range. Combining with the proposal to solve problem 1 above, one could then check whether one has to enlarge the bounds. If not, one can finally use the old plot recipe's code to plot the new vector of purely bounded sets.

Alternatively, for problems 2 and 3 there should be a way to hook into the update procedure of a plot and, if the range has changed, to replace the previous sets. The problem might be (I do not know) that this may not happen in Plots.jl but rather in the backends.

Side notes

1. I think radius_low is just -radius_high by construction, so you just need one of them.

2. I found this function, which preprocesses shape data and enables infinite horizontal and vertical shapes. But it relies on the fact that only one of x/y is infinite. So we cannot use that. I think the first two lines also show the better way to obtain the plot range.

[ueliwechsler]

As far as I know (Documentation and trials) it is not supported.

As an addition to this comment, I found here https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/master/docs/src/plot.md this sentence

A 3-dimensional polyhedron can be visualized with either MeshCat or Makie. Unbounded polyhedron are supported by truncating the polyhedron into a polytope and not triangularizing the faces in the directions of unbounded rays.

But for 2D-set, there is no information on plotting unbounded sets.

[schillic]

A 3-dimensional polyhedron can be visualized with either MeshCat or Makie. Unbounded polyhedron are supported by truncating the polyhedron into a polytope and not triangularizing the faces in the directions of unbounded rays.

But for 2D-set, there is no information on plotting unbounded sets.

I guess this is something that happens in Polyhedra, and that code should actually be independent of the dimension.

[ueliwechsler]

I think it is because for 2D sets Plots.jl is used and Makie.jl for 3D sets.

[schillic]

If Makie supports unbounded sets, that would be optimal, but it is not clear from that sentence.

Otherwise, the bounding has to happen on our end, and for that I still think my comment above is the most straightforward thing to do: assume the user has defined the proper plot ranges (which corresponds to a box in 2D and a cube in 3D, i.e., in general a hyperrectangle), read them out, bound the set by intersecting with the hyperrectangle, and then plot it.

For plotting a vector of sets at once, we should prioritize bounded sets such that we get reasonable plot ranges.

I see @ueliwechsler's approach here as a possible extension.

An alternative extension would be to have new functions uplot(X) and uplot!(X) that store all sets in a global variable. As long as the storage does not contain an unbounded set, uplot!(X) just calls plot!(X). As soon as there is an unbounded set in the storage, whenever uplot! is called again, it executes the following procedure: read out the current plot bounds, clear the plot, and replot everything. Pseudocode:

global plot_storage = PlotStorage()

mutable struct PlotStorage
    sets::Vector{LazySet}
    isbounded::Bool

    function PlotStorage() = new(LazySet[], true))  # default constructor
end

function uplot(X)
    # reset global variable
    plot_storage.sets = [X]
    plot_storage.isbounded = isbounded(X)
    # plot set
    plot(X)
end

function uplot!(X)
    push!(plot_storage.sets, X)

    if plot_storage.isbounded
        # just plot new set
        plot_storage.isbounded = isbounded(X)
        plot!(X)
    else
        # replot all sets
        *clean_plot()*  # clean plot but retain current bounds
        for Xi in plot_storage.sets
            plot!(Xi)
        end
    end
end

[ueliwechsler]

Otherwise, the bounding has to happen on our end, and for that I still think my comment above is the most straightforward thing to do:

Does your method plot_halfspace apply to polytopes with several halfspaces? Because I think a method, where we plot unbounded sets as bounded sets that have the limit well beyond the x and y limits of the plot pane would be a really simple solution.

I'd like to elaborate on my idea described in https://github.com/JuliaReach/LazySets.jl/issues/576#issuecomment-522757875 and show how the new implementation would handle the drawbacks mentioned in https://github.com/JuliaReach/LazySets.jl/issues/576#issuecomment-522865909 and what new drawbacks would appear.

Consider the following implementation:

const PLOT_LIMIT = 1000

function bounded_approximation(P::HPolyhedron{N}) where {N}
    A, a = tosimplehrep(P)
    n = LazySets.dim(P)
    B = [Matrix{N}(I,n,n); -Matrix{N}(I,n,n)]
    b = repeat([PLOT_LIMIT],2*n)
    return HPolytope([A;B], [a;b])
end

function plot_polyhedron!(P::HPolyhedron{N}; kwargs...) where {N}
    p = plot!()
    xlim = xlims(p)
    ylim = ylims(p)
    P_bounded = bounded_approximation(P)
    if isempty(P_bounded)
        error("resulting bounded polyhedron derived from the unbounded polyhedron "
              *"is an empty set. try increasing the plotting limit PLOT_LIMIT")
    end
    plot!(P_bounded; kwargs...)
    xlims!(xlim)
    ylims!(ylim)
end

We define a constant PLOT_LIMIT which defines the size of the bounded_approximation of the unbounded set. In plot_polyhedron we extract the current xlims and ylims and after plotting the bounded approximation of the set, we set the limits to the old limits.

As a result:

Problems

1. You have to statically guess the plot area. [-10, 10] is probably a good range in most cases, but it is easy to think of a counterexample.

We do not guess the limits, we just keep the existing ones. Since it is infeasible hard anyways to guess the range, we just keep the existing ones and let the user decide. If there is already a set plotted the range does not change. In a next step, we could add a heuristic for estimating the interesting range.

2. If I understand correctly, the user needs to actively redraw, which is annoying.

This would be avoided since the set is large enough.

3. The previous versions of the unbounded sets are not removed, so with transparency enabled you still see those.

This would also not matter anymore.

However, there is one general disadvantage to this approach.

• After plotting an unbounded set with this approach, the plot windows does not automatically adjust to the plotting windows (i.e. xlims and ylims). See the example below. I am not sure how to avoid this.

Another bug is apparent which is IMO independent of my implementation:

• If the plot window is zoomed in to strong, only the line but not the area is plotted. Also, see the example below.

A =[ -0.00642001   0.451886;
      0.238524     0.141031;
     -0.0829217   -0.178183]
b = [-0.325, -0.0978, 0.307]
P = HPolytope(A,b)
P1 = HPolyhedron([0.5 0.5], -[1.50])
P2 = HPolyhedron([-0.5 0.5], -[9.0])
P3 = HPolyhedron([-1 0.], -[9.0])

plot(P)
# if the plot view is to small, there is no area plottted
plot_polyhedron!(P1)
plot_polyhedron!(P2)
plot_polyhedron!(P3)
ylims!((-15,15)); xlims!((-15,15))
# the plot is not updated anymore! (we don't see the singleton)
plot!(Singleton([15,15]))

plot()
plot_polyhedron!(P1)
ylims!((-4,2)); xlims!((-4,2))
# if the limit is to small, there is no area plottted
ylims!((-3,2)); xlims!((-3,2))

[schillic]

In #1909 I sent a proposal. I do not say that it is perfect but it is very simple and it is integrated in the Plots environment, which allows you to "just call plot/plot! as usual". Let me know what you think.

[schillic]

Does your method plot_halfspace apply to polytopes with several halfspaces?

The code above did not, but the idea described further down does (see the PR, which works with any set for which we can compute an intersection (which is certainly the case for AbstractPolyhedra)).