#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"`  -- evaluate the support function by reducing the problem to
+                         the 2D intersection of a rank 2 linear transformation of
+                         the given compat set in the plane generated by the given
+                         direction `d` and the hyperplane's normal vector `n`;
+                         only valid for intersection with a hyperplane   
 
 - `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,38 @@ 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}}; kwargs...) where {N}
+
+Given a compact and convex set ``X`` and a hyperplane ``H = \\{x: n ⋅ x = γ \\}``,
+calculate ...
+
+### Input
+
+- `ℓ`      -- direction
+- `X`      -- set
+- `H`      -- hyperplane
+
+### Output
+
+The support function of ``X ∩ H`` along direction ``ℓ``.
+"""
+function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}};
+                     lazy_linear_map=false,
+                     lazy_2d_intersection=true,
+                     kwargs...) where {N}
+
+    options = Dict(kwargs)
+    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)
+    cap = lazy_2d_intersection ? Intersection(Snℓ, Lγ) : intersection(Snℓ, Lγ)
+    return ρ(y_dir, cap; kwargs...)
+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)