#676 - Add _projection method for set - hyperplane lazy intersection

docs/src/lib/operations.md

@@ -102,7 +102,7 @@ dim(::Intersection)
 σ(::AbstractVector{Real}, ::Intersection{Real})
 ∈(::AbstractVector{Real}, ::Intersection{Real})
 isempty(::Intersection)
-ρ(::AbstractVector{N}, ::Intersection{N, <:LazySet, S}) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane}}
+ρ(::AbstractVector{N}, ::Intersection{N, <:LazySet, S}) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane, Line}}
 ```
 
 Inherited from [`LazySet`](@ref):

src/Intersection.jl

@@ -293,7 +293,7 @@ end
       cap::Intersection{N, <:LazySet, S};
       algorithm::String="line_search",
       check_intersection::Bool=true,
-      kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane}}
+      kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane, Line}}
 
 Return the support function of the intersection of a compact set and a half-space
 or a hyperplane in a given direction.
@@ -308,6 +308,11 @@ or a hyperplane in a given direction.
     * `"line_search"` -- solve the associated univariate optimization problem
                          using a line search method (either Brent or the
                          Golden Section method) 
+    * `"projection"`  -- only valid for intersection with a hyperplane;
+                         evaluates the support function by reducing the problem to
+                         the 2D intersection of a rank 2 linear transformation of
+                         the given compact set in the plane generated by the given
+                         direction `d` and the hyperplane's normal vector `n`
 
 - `check_intersection` -- (optional, default: `true`) if `true`, check if the 
                           intersection is empty before actually calculating the
@@ -350,7 +355,7 @@ function ρ(d::AbstractVector{N},
            cap::Intersection{N, <:LazySet, S};
            algorithm::String="line_search",
            check_intersection::Bool=true,
-           kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane}}
+           kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane, Line}}
 
     X = cap.X    # compact set
     H = cap.Y    # halfspace or hyperplane
@@ -362,14 +367,24 @@ function ρ(d::AbstractVector{N},
 
     if algorithm == "line_search"
         @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " *
-            "the package 'Optim' to be loaded"
+                                        "the package 'Optim' to be loaded"
         (s, _) = _line_search(d, X, H; kwargs...)
+    elseif algorithm == "projection"
+        @assert H isa Hyperplane "the algorithm $algorithm cannot be used with a
+                                  $(typeof(H)); it only works with hyperplanes"
+        s = _projection(d, X, H; kwargs...)
     else
         error("algorithm $(algorithm) unknown")
     end
     return s
 end
 
+# Symmetric case
+ρ(ℓ::AbstractVector{N}, cap::Intersection{N, S, <:LazySet};
+  algorithm::String="line_search", check_intersection::Bool=true,
+  kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane, Line}} = ρ(ℓ, cap.Y ∩ cap.X;
+  algorithm=algorithm, check_intersection=check_intersection, kwargs...)
+
 function load_optim_intersection()
 return quote
 
@@ -436,7 +451,7 @@ julia> _line_search([1.0, 0.0], X, H, upper=1e3, method=GoldenSection())
 (1.0, 381.9660112501051)
 ```
 """
-function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane}; kwargs...)
+function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane, Line}; kwargs...)
     options = Dict(kwargs)
 
     # Initialization
@@ -448,7 +463,7 @@ function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane}; kwargs...)
     else
         if H isa HalfSpace
             lower = 0.0
-        elseif H isa Hyperplane
+        elseif (H isa Hyperplane) || (H isa Line)
             lower = -1e6 # "big": TODO relate with f(λ)
         end
     end
@@ -475,8 +490,80 @@ end # _line_search
 end # quote
 end # load_optim
 
-# Symmetric case
-ρ(ℓ::AbstractVector{N}, cap::Intersection{N, S, <:LazySet};
-  algorithm::String="line_search", check_intersection::Bool=true,
-  kwargs...) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane}} = ρ(ℓ, cap.Y ∩ cap.X;
-  algorithm=algorithm, check_intersection=check_intersection, kwargs...)
+"""
+    _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
+                lazy_linear_map=false,
+                lazy_2d_intersection=true,
+                algorithm_2d_intersection=nothing,
+                kwargs...) where {N}
+
+Given a compact and convex set ``X`` and a hyperplane ``H = \\{x: n ⋅ x = γ \\}``,
+calculate the support function of the intersection between the
+rank-2 projection ``Π_{nℓ} X`` and the line ``Lγ = \\{(x, y): x = γ \\}``.
+
+### Input
+
+- `ℓ`                    -- direction
+- `X`                    -- set
+- `H`                    -- hyperplane
+- `lazy_linear_map`      -- (optional, default: `false`) to perform the projection
+                            lazily or concretely
+- `lazy_2d_intersection` -- (optional, default: `true`) to perform the 2D
+                            intersection between the projected set and the line
+                            lazily or concretely
+- `algorithm_2d_intersection` -- (optional, default: `nothing`) if given, fixes the
+                                 support function algorithm used for the intersection
+                                 in 2D; otherwise the default is implied
+
+### Output
+
+The support function of ``X ∩ H`` along direction ``ℓ``.
+
+### Algorithm
+
+This projection method is based on Prop. 8.2, page 103, [C. Le Guernic.
+Reachability Analysis of Hybrid Systems with Linear Continuous Dynamics,
+PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2).
+
+In the original algorithm, Section 8.2 of Le Guernic's thesis, the linear map
+is performed concretely and the intersection is performed lazily (these are the
+default options in this algorithm, but here the four combinations are available).
+If the set ``X`` is a zonotope, its concrete projection is again a zonotope
+(sometimes called "zonogon"). The intersection between this zonogon and the line
+can be taken efficiently in a lazy way (see Section 8.2.2 of Le Guernic's thesis),
+if one uses dispatch on `ρ(y_dir, Sℓ⋂Lγ; kwargs...)` given that `Sℓ` is itself
+a zonotope.
+
+### Notes
+
+This function depends itself on the calculation of the support function of another
+set in two dimensions. Obviously one doesn't want to use again `algorithm="projection"`
+for this second calculation. The option `algorithm_2d_intersection` is such that,
+if it is not given, the default support function algorithm is used (e.g. `"line_search"`).
+You can still pass additional arguments to the `"line_search"` backend through the
+`kwargs`.
+"""
+function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
+                     lazy_linear_map=false,
+                     lazy_2d_intersection=true,
+                     algorithm_2d_intersection=nothing,
+                     kwargs...) where {N}
+
+    n = H.a                  # normal vector to the hyperplane
+    γ = H.b                  # displacement of the hyperplane
+    Πnℓ = vcat(n', ℓ')       # projection map
+
+    x_dir = [one(N), zero(N)]
+    y_dir = [zero(N), one(N)]
+
+    Lγ = Line(x_dir, γ)
+
+    Snℓ = lazy_linear_map ? LinearMap(Πnℓ, X) : linear_map(Πnℓ, X)
+    Sℓ⋂Lγ = lazy_2d_intersection ? Intersection(Snℓ, Lγ) : intersection(Snℓ, Lγ)
+
+    if algorithm_2d_intersection == nothing
+        return ρ(y_dir, Sℓ⋂Lγ; kwargs...)
+    else
+        return ρ(y_dir, Sℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...)
+    end
+end

test/unit_Intersection.jl

@@ -67,6 +67,19 @@ d = normalize([1.0, 0.0])
 using Optim
 @test ρ(d, X ∩ H) == ρ(d, H ∩ X) == 1.0
 
-# Hyperplane - Ball1 intersection
+# Ball1 vs HalfSpace intersection using default algorithm
 H = HalfSpace([1.0, 0.0], 0.0); # x = 0
 @test ρ(d, X ∩ H) < 1e-6 && ρ(d, H ∩ X) < 1e-6
+# specify  line search algorithm
+@test ρ(d, X ∩ H, algorithm="line_search") < 1e-6 &&
+      ρ(d, H ∩ X, algorithm="line_search") < 1e-6
+
+# Ball1 vs Hyperplane intersection
+H = Hyperplane([1.0, 0.0], 0.5); # x = 0.5
+@test isapprox(ρ(d, X ∩ H, algorithm="line_search"), 0.5, atol=1e-6)
+# For the projection algorithm, if the linear map is taken lazily we can use Ball1
+@test isapprox(ρ(d, X ∩ H, algorithm="projection", lazy_linear_map=true), 0.5, atol=1e-6)
+# But the default is to take the linear map concretely; in this case, we *may*
+# need Polyhedra (in the general case), for the concrete linear map. As a valid workaround
+# if we don't want to load Polyhedra here is to convert the given set to a polygon in V-representation
+@test isapprox(ρ(d, convert(VPolygon, X) ∩ H, algorithm="projection", lazy_linear_map=false), 0.5, atol=1e-6)