This is a continuation of #68. There we added containment checks (x ∈ S) for most set types. Some are still missing.

• HPolytope (#167)
• Intersection (#167)
• ConvexHull
  A sufficient test is to check containment in both generating sets, but for negative answers one has to continue, e.g., checking ray intersection from some vertex of one set in the other set.
• ExponentialProjectionMap
• MinkowskiSum
  One could subtract one set from the point and check for nonempty intersection with the other set. There is literature for the case that the sum consists of polytopes only, this could be interesting.

We have ρ(d::AbstractVector{N}, S::LazySet)::N where {N<:Real}.
This will crash for Ball2, even with rational direction. None of the below work:

ρ([1//1], Ball2([1//1], 1//1))
σ([1//1], Ball2([1//1], 1//1))
ρ([1//1], Ball2(1., 1.))
σ([1//1], Ball2(1., 1.))

The problem is that Ball2 (and Ballp) is only really usable with AbstractFloat types.
I guess there is no way around. Should we instead require that these types cannot even be created with numeric types that are above AbstractFloat? At least this would avoid confusion.

Following #7 , we add new @recipe functions for any lazy set, either exactly if possible or through a polyhedral approximation.

1. Rename the file to the singular form?
2. Outsource Line to a new file?
3. Should Line's field a::Vector{N} become a::Tuple{N, N} because the constructor requires 2D input?
   On the other hand, code-wise I do not see any problem with a Line having more than two dimensions. We only use it for 2D plots, but why should we restrict the general type? If we want to keep it 2D, I would rename it to Line2 or so.
4. Should we promote Line to LazySet?
5. (related) Should we have a LazySet type for a half-space? This is inspired by the doc for LinearConstraint.

Action points here:

• rename LinearConstraints.jl to Line.jl and promote Line to LazySet
• add LinearConstraint as alias to HalfSpace and remove the type
• outsource intersection function to a new file concrete_intersection.jl

We have a lot of subtypes, e.g., of LazySet.
We want all of them to provide a method for the function σ.
Hence we should check that those methods exist to prevent new types that may break things.

We can get all subtypes with the help of subtypes(LazySet) and use the method_exists function to check for existence:

for S in subtypes(LazySet)
    println(method_exists(σ, (Vector{Float64}, S)))
end

I propose that we add such tests for all abstract types, also in Reachability.

We could even create a function that automates this process if we have several common methods:
Add in LazySets.jl:

function require()::Vector{Function}
    return [
        S -> (σ, (Vector{Float64}, S)),
        S -> (my_func1, (S)),
        S -> (my_func2, (S, Real, Int))
    ]
end

Then in the unit tests we just have to go through this vector and apply the test to each of them.
We probably need a mechanism to allow several require methods for different abstract types, e.g., by changing the name.

(I have not tested the second one, the first one works.)

Julia still does not support deactivating evaluation of @assert macros.

I played around with the following example test.jl.

n = 10000000
@assert zeros(n) == zeros(n)

Performance:

@time include("test.jl")
0.146106 seconds (139 allocations: 152.596 MiB, 64.41% gc time)

Let us add the following lines on top that try to deactivate the check:

import Base.@assert
macro assert(arg) end

Performance:

@time include("test.jl")
0.002707 seconds (135 allocations: 9.344 KiB)

Now remove the whole @assert line.
Performance:

0.000995 seconds (112 allocations: 7.609 KiB)

You can see that writing @assert does not cost nothing (probably because it still calls the empty macro), but that it does not evaluate the code that is passed.
I find this good enough to be used.

For arbitrary 3D lazy sets, we lack an overapproximate in that dim that can be used for visualization. That would need a generalization of the Lotov-like algorithm, but in 3D.

In the meantime, we could:

• plot the result from ball-inf approximations (eps = Inf case)
• plot the cartesian decomposition of the given set, as valid (2+1 dim) overapproximation

Note: for the Singleton type this already works, e.g.

julia> using IterTools
julia> plot([Singleton(collect(vi)) for vi in IterTools.product([[1, -1] for i = 1:3]...)])

--

Edit

useful links:

• https://discourse.julialang.org/t/pyplot-3d-plotting/30838
• https://www.geeksforgeeks.org/how-to-draw-3d-cube-using-matplotlib-in-python/
• https://stackoverflow.com/questions/65482345/plotting-a-3d-surface-in-julia-using-either-plots-or-pyplot
• https://discourse.julialang.org/t/plotting-a-3d-surface/17143/16
• https://matplotlib.org/2.0.2/mpl_toolkits/mplot3d/api.html
• https://github.com/JuliaPy/PyPlot.jl#3d-plotting

Add a subset check ⊆(S1, S2), in particular for special types of S1.
In particular, for polytopes (and "subtypes"), since S2 is convex, it should suffice to check containment of all corners (see #68).

Also check if Julia has an alias for ⊆, otherwise add one such as isincluded or contains.

This happens if the lines are parallel or identical.

l1 = Line([1., 1.], 1.)
l2 = Line([1., 1.], 1.)
intersection(l1, l2) # Base.LinAlg.SingularException(2)
l3 = Line([1., 0.], 1.)
l4 = Line([1., 0.], 1.)
intersection(l3, l4) # Base.LinAlg.LAPACKException(2)

When operating with canonical basis vectors, it's a bit faster to use just the same single vector which is mutated, instead of e.g. [zeros(i-1); 1.; zeros(n-i)].

This can be applied to Approximations.ballinf_approximation

Creating the documentation warns about potentially missing docstrings.
However, they are defined and occur properly in the final documentation.

Most notably:

• LazySets.σ
• LazySets.apply_recipe

For the non-apply_recipe functions the problem is related to type parameters.

Currently the norm of BallInf is implemented as follows.

return maximum(map(x -> norm(x, p), vertices_list(B)))

We could instead find the corner that is furthest from the origin (by simply adding/subtracting radius from center[i], depending on the sign of center[i], for each dimension i) and then immediately return the norm of this one vertex:

return norm((@. c + unit_step(c) * r), p)

EDIT: I had to change sign to unit_step for the 0 case.

@mforets: Am I right?

A similar idea should be applicable for Hyperrectangle.

The symmetric interval hull is computed explicitly in Approximations.box_approximation_helper, through 2n support function evaluations.

An application of a new lazy (total or partial) type, symmetric interval hull, is to speedup the discretization phase; the explicit form can take up to 97% of the time for n>10000.

We should review that for d=0 the support vector returns a point in the set, or to consider the special cases that raise an error otherwise.

Just need to review where we actually use mutability and otherwise change it...

for support function computations, see: Reachability analysis of linear systems with linear continuous dynamics, page 19.

We have hrep -> vrep, but we are missing the converse.

See, e.g., BallInf.vertices_list(B::BallInf).

Have a look at the type Generator.

Types such as CartesianProductArray consist of several instances of LazySets.
Often all the elements have the same type, which we want to exploit by introducing a type parameter.
This helps in polymorphic function definitions where we specialize to subtypes and could improve general efficiency as the compiler can create a contiguous array.

This is probably intended, but not required.

We can then also simplify the radius function for the infinity norm: return (p == Inf) ? B.radius : ...

@mforets: Currently there is a failing doctest in master:

> File: .../LazySets/src/Ball1.jl

> Code block:                                                                                                       

\`\`\`jldoctest
julia> σ([0.,1], B)
1-element Array{Float64,1}:
 1.0
 0.0
\`\`\`

> Subexpression:                                                                                                    

    σ([0.,1], B)

> Output Diff:                                                                                                      

1-element 2-element Array{Float64,1}:
 1.0
 0.00.0
 1.0 

However, Travis does not complain.
See here how to potentially implement this.

• A default implementation could just pick any direction and ask for a support vector.
• The empy set should throw an exception.
• What would be a good name for this function? element? CC: @mforets

Review plot_polygon (-> plot ?) and test different backends: PyPlot, Gadfly, Plots.

Quite often we have the following code patterns:

d::Union{Vector{Float64}, SparseVector{Float64,Int64}}
M::Union{Matrix{Float64}, SparseMatrixCSC{Float64,Int64}}

We should use instead the common supertype:

d::AbstractVector{Float64}
M::AbstractMatrix{Float64}

This is easier to read and gives more flexibility (e.g., allow SubArrays if not already mentioned).
I do not expect negative performance impact; we should of course validate this.

Add a function is_intersection_empty(S1, S2; witness=false) that returns true if the sets S1 and S2 do not intersect.
If the flag witness is set, also produce a witness, i.e., a point (vector) x such that x ∈ S1 ∩ S2 if the intersection is nonempty.
The table below tracks the current combinations that we support:

Related: #139

let's create a new examples folder with some illustrative examples and pictures:

• polygon overapproximations => iterative refinement
• convex hulls
• interval hulls
• reachability algorithm for linear systems using zonotopes

cc: @asarvind

The complex zonotope is defined is a generalization of usual zonotopes, where the set of generators can contain complex-valued vectors.

To compute the support vector, does it make sense to use an analogue formula than that of the real-valued zonotope?

See also:

• Using complex zonotopes for stability verification. Arvind S. Adimoolam, Thao Dang.
• https://theses.hal.science/tel-01907688

In applications is sometimes useful to instantiate an Hyperrectangle as the pair of vectors low and high, instead of center and radius (both vectorial).

For example, consider:

julia> using LazySets
julia> B = BallInf([0.5, 0.5], 0.1)
julia> import LazySets.Approximations: overapproximate

The transformation to a polygon without refinement works as expected:

julia> vertices_list(overapproximate(B, 1.))
4-element Array{Array{Float64,1},1}:
 [0.6, 0.6]
 [0.4, 0.6]
 [0.4, 0.4]
 [0.6, 0.4]

However, extra vertices are created if we decrease the error tolerance:

julia> vertices_list(overapproximate(B, 1e-3))
5-element Array{Array{Float64,1},1}:
 [0.6, 0.6]
 [0.4, 0.6]
 [0.4, 0.4]
 [0.4, 0.4]
 [0.6, 0.4]

This is to be revised and possibly avoided, since many applications of
overapproximate would generate a lot of needless vertices.

in our use case, we need to handle small number of points (<~100).

see also: https://en.wikipedia.org/wiki/Convex_hull_algorithms

moreover, sorting the points is necessary for the supp vector implementation we plan to do, so a good candidate is Graham scan algorithm.

pointers to Julia libraries:

We can just check for radius 0 and only return the center in this case.
This is also an improvement if the radius is 0 only in some dimensions.

Types: Ball1, BallInf, Hyperrectangle.

In #13 we implemented a O(n) algorithm to compute the support vector of a VPolygon. Since the vertices are sorted in CCW fashion by construction, we can do better than that using a binary search, see references below.

• extreme points in convex polygons in the geomalgorithms site
• Joseph O'Rourke, Computational Geometry in C (2nd Edition), Sect 7.9 "Extreme Point of Convex Polygon" (1998).

Add a containment check ∈(x, S) for all sets.

For some we already have is_contained. I prefer in instead (the Julia alias for ∈).

for "v1.0" we need to finish polishing the docs. on a high-level:

• library of functions, include them all
• beginning example: convex hull, M-sum, lazy mexp
• finish "About" section
• add a wiki page describing the documentation style that we use
  (we use it for Reachability.jl as well; eventually this will be unified for all repositories)
• polyhedral approximations: resize figures and add text ssec 3
• decomposing an affine map tutorial (from paper)
• fast 2d LPs: resize fig and add Algorithm (from paper)
• small examples for a "getting started" section, use jupyter notebooks (?) << we used @examples blocks
• common set reps: review repeated docs of sigma << this is deferred to #72
• add Polyhedra to the docs
• clarify questions here

given a set of 2d vertices, the support vector of a VPolygon <:LazySet can be computed by:

• sorting vertices by polar angle (see function <=)
• performing binary search for the given direction

a convex hull preprocessing maybe needed, in general.

in a similar manner to HPolygonOpt, there can be a VPolygonOpt for warm-starting.

as suggested here, we can turn our examples into proper doctests.

with Documenter.jl, this boils down to using jldoctest in triple backticks for each example. i don't know if we need to have the line using LazySets in each code block.

Instead of

diameter(B::BallInf, p::Real=Inf)
diameter(H::Hyperrectangle, p::Real=Inf)
...

just define it once

diameter(S::LazySet, p::Real=Inf)

(In fact we have this, but it throws an error.)

@mforets: Any objection?

The code of iterative refinement should be included in the docs and explained (see iterative_refinement.md under docs/src/man).

In some cases we'd like to have eltype(X) where X is a LazySet, to return the underlying numeric type. Eg. eltype(BallInf(ones(2), 0.1)) would be Float64.

In principle this would just require

abstract type LazySet{N<:Number} end

But we should see if this makes sense for all LazySets, both set representations and operations.

The "VoidSet" is stated as the neutral element for "+" and the absorbing element for "*".

I don't understand the name (as a neutral element for "+"), it refers to a set that is empty. The empty set is the absorbing element for the Minkowski sum. Maybe "ZeroSet" is more relevant.

On the other hand, the empty set is the absorbing element for the Cartesian product.

I think it is a bit confusing, maybe we have to split both concepts?

It can be interesting to accept parametric number types so that we can work with:

• rationals (slower but exact)
• single precision (less precise but faster)
• complex numbers (eg. for complex zonotopes)

Let

Ω0 = LazySets.ConvexHull{LazySets.CartesianProductArray{LazySets.BallInf},LazySets.LinearMap{LazySets.CartesianProductArray{LazySets.BallInf}}}(LazySets.CartesianProductArray{LazySets.BallInf}(LazySets.BallInf[LazySets.BallInf([0.767292, 0.936613], 0.1), LazySets.BallInf([0.734104, 0.87296], 0.1)]), LazySets.LinearMap{LazySets.CartesianProductArray{LazySets.BallInf}}([1.92664 1.00674 1.0731 -0.995149; -2.05704 3.48059 0.0317863 1.83481; 0.990993 -1.97754 0.754192 -0.807085; -2.43723 0.782825 -3.99255 3.93324], LazySets.CartesianProductArray{LazySets.BallInf}(LazySets.BallInf[LazySets.BallInf([0.767292, 0.936613], 0.1), LazySets.BallInf([0.734104, 0.87296], 0.1)])))

then

julia> Approximations.decompose(Ω0).sfarray[1].constraints_list

4-element Array{LazySets.LinearConstraint,1}:
 LazySets.LinearConstraint([1.0, 0.0], 2.63908)  
 LazySets.LinearConstraint([0.0, 1.0], 4.04073)  
 LazySets.LinearConstraint([-1.0, 0.0], -2.04145)
 LazySets.LinearConstraint([0.0, -1.0], -2.5726) 

although

julia> Approximations.decompose(Ω0, Inf).sfarray[1].constraints_list

4-element Array{LazySets.LinearConstraint,1}:
 LazySets.LinearConstraint([1.0, 0.0], 2.84043)   
 LazySets.LinearConstraint([0.0, 1.0], 4.04709)   
 LazySets.LinearConstraint([-1.0, 0.0], -0.667292)
 LazySets.LinearConstraint([0.0, -1.0], -0.836613)

The only type union not affected by #19 exists in Polgyon.jl to handle HPolygons and HPolygonOpts simultaneously.

We should probably just make HPolygonOpt inherit from HPolygon.
This makes sense both from the code and from the conceptual side to me.
I am, however, not sure if Julia finds the correct function if one argument is an array of a subtype; I guess it does, but this has to be tested.

@mforets: What do you think?

Initial motivation: supertype for HPolygon and HPolygonOpt:

1. Avoid copy-pasting methods in both HPolygon and HPolygonOpt.
2. Allows to restrict a type to a polygon irrespective of the implementation.

The idea also translates to other types such as BoxSet(?) for Hyperrectangle/BallInf.

Another idea is to have different implementations of array-type sets (MinkowskiSumArray etc.) that store their dimension for faster access.

So the general message here is:
Introduce interface types and usually only work with them.
A crucial implication is that we should never access member/field variables directly (because different implementations may not have those members).
Instead there should be a function that is usually just a getter.