From f26aa1b07d68b9aa9d3204ff835152a1a7dd6432 Mon Sep 17 00:00:00 2001 From: mforets Date: Sat, 13 Oct 2018 00:10:07 -0300 Subject: [PATCH 1/6] remove visible upper_bounds --- src/Approximations/Approximations.jl | 2 +- src/Approximations/box_approximations.jl | 46 ++-- src/Approximations/overapproximate.jl | 51 ++--- src/Intersection.jl | 272 +++++++++++------------ src/is_intersection_empty.jl | 6 +- 5 files changed, 161 insertions(+), 216 deletions(-) diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 85f3fdaddd..536f2ec3b7 100644 --- a/src/Approximations/Approximations.jl +++ b/src/Approximations/Approximations.jl @@ -37,6 +37,6 @@ include("box_approximations.jl") include("template_directions.jl") include("overapproximate.jl") include("decompositions.jl") -include("upper_bounds.jl") +#include("upper_bounds.jl") end # module diff --git a/src/Approximations/box_approximations.jl b/src/Approximations/box_approximations.jl index bf9f2bd3e3..920872c6f7 100644 --- a/src/Approximations/box_approximations.jl +++ b/src/Approximations/box_approximations.jl @@ -3,15 +3,13 @@ # =================================== """ - box_approximation(S::LazySet; upper_bound::Bool=false)::Hyperrectangle + box_approximation(S::LazySet)::Hyperrectangle Overapproximate a convex set by a tight hyperrectangle. ### Input - `S` -- convex set -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? ### Output @@ -24,20 +22,19 @@ of the given set in the canonical directions, and the lengths of the sides can be recovered from the distance among support functions in the same directions. """ function box_approximation(S::LazySet{N}; - upper_bound::Bool=false )::Union{Hyperrectangle{N}, EmptySet{N}} where N<:Real - (c, r) = box_approximation_helper(S; upper_bound=upper_bound) - if upper_bound && r[1] < 0 + (c, r) = box_approximation_helper(S) + if r[1] < 0 return EmptySet{N}() end return Hyperrectangle(c, r) end # special case: Hyperrectangle -box_approximation(S::Hyperrectangle; upper_bound::Bool=false) = S +box_approximation(S::Hyperrectangle) = S # special case: other rectangle -box_approximation(S::AbstractHyperrectangle; upper_bound::Bool=false) = +box_approximation(S::AbstractHyperrectangle) = Hyperrectangle(center(S), radius_hyperrectangle(S)) """ @@ -48,7 +45,7 @@ Alias for `box_approximation`. interval_hull = box_approximation """ - box_approximation_symmetric(S::LazySet{N}; upper_bound::Bool=false + box_approximation_symmetric(S::LazySet{N} )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real} @@ -57,8 +54,6 @@ Overapproximate a convex set by a tight hyperrectangle centered in the origin. ### Input - `S` -- convex set -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? ### Output @@ -70,11 +65,10 @@ The center of the box is the origin, and the radius is obtained by computing the maximum value of the support function evaluated at the canonical directions. """ function box_approximation_symmetric(S::LazySet{N}; - upper_bound::Bool=false )::Union{Hyperrectangle{N}, EmptySet{N}} where {N<:Real} - (c, r) = box_approximation_helper(S; upper_bound=upper_bound) - if upper_bound && r[1] < 0 + (c, r) = box_approximation_helper(S) + if r[1] < 0 return EmptySet{N}() end return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r) @@ -88,7 +82,7 @@ Alias for `box_approximation_symmetric`. symmetric_interval_hull = box_approximation_symmetric """ - box_approximation_helper(S::LazySet{N}; upper_bound::Bool=false + box_approximation_helper(S::LazySet{N}; ) where {N<:Real} Common code of `box_approximation` and `box_approximation_symmetric`. @@ -96,8 +90,6 @@ Common code of `box_approximation` and `box_approximation_symmetric`. ### Input - `S` -- convex set -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? ### Output @@ -115,24 +107,22 @@ The lengths of the sides can be recovered from the distance among support functions in the same directions. """ @inline function box_approximation_helper(S::LazySet{N}; - upper_bound::Bool=false ) where {N<:Real} zero_N = zero(N) one_N = one(N) - ρ_rec = upper_bound ? ρ_upper_bound : ρ n = dim(S) c = Vector{N}(undef, n) r = Vector{N}(undef, n) d = zeros(N, n) @inbounds for i in 1:n d[i] = one_N - htop = ρ_rec(d, S) + htop = ρ(d, S) d[i] = -one_N - hbottom = -ρ_rec(d, S) + hbottom = -ρ(d, S) d[i] = zero_N c[i] = (htop + hbottom) / 2 r[i] = (htop - hbottom) / 2 - if upper_bound && r[i] < 0 + if r[i] < 0 # contradicting bounds => set is empty # terminate with first radius entry being negative r[1] = r[i] @@ -143,7 +133,7 @@ functions in the same directions. end """ - ballinf_approximation(S::LazySet{N}; upper_bound::Bool=false + ballinf_approximation(S::LazySet{N}; )::BallInf{N} where {N<:Real} Overapproximate a convex set by a tight ball in the infinity norm. @@ -151,8 +141,6 @@ Overapproximate a convex set by a tight ball in the infinity norm. ### Input - `S` -- convex set -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? ### Output @@ -164,26 +152,24 @@ The center and radius of the box are obtained by evaluating the support function of the given convex set along the canonical directions. """ function ballinf_approximation(S::LazySet{N}; - upper_bound::Bool=false )::Union{BallInf{N}, EmptySet{N}} where {N<:Real} zero_N = zero(N) one_N = one(N) - ρ_rec = upper_bound ? ρ_upper_bound : ρ n = dim(S) c = Vector{N}(undef, n) r = zero_N d = zeros(N, n) @inbounds for i in 1:n d[i] = one_N - htop = ρ_rec(d, S) + htop = ρ(d, S) d[i] = -one_N - hbottom = -ρ_rec(d, S) + hbottom = -ρ(d, S) d[i] = zero_N c[i] = (htop + hbottom) / 2 rcur = (htop - hbottom) / 2 if (rcur > r) r = rcur - elseif upper_bound && rcur < 0 + elseif rcur < 0 # contradicting bounds => set is empty return EmptySet{N}() end diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 83680ed173..1e5a17c12b 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -30,7 +30,6 @@ Otherwise the result is an ε-close approximation as a polygon. - `S` -- convex set, assumed to be two-dimensional - `HPolygon` -- type for dispatch - `ε` -- (optional, default: `Inf`) error bound -- `upper_bound` -- (optional, default: `false`) currently ignored ### Output @@ -38,8 +37,7 @@ A polygon in constraint representation. """ function overapproximate(S::LazySet{N}, ::Type{<:HPolygon}, - ε::Real=Inf; - upper_bound::Bool=false + ε::Real=Inf )::HPolygon where {N<:Real} @assert dim(S) == 2 if ε == Inf @@ -59,15 +57,14 @@ function overapproximate(S::LazySet{N}, end """ - overapproximate(S::LazySet, ε::Real; [upper_bound]::Bool=false)::HPolygon + overapproximate(S::LazySet, ε::Real)::HPolygon -Alias for `overapproximate(S, HPolygon, ε, upper_bound=upper_bound)`. +Alias for `overapproximate(S, HPolygon, ε)`. """ function overapproximate(S::LazySet, ε::Real; - upper_bound::Bool=false )::HPolygon - return overapproximate(S, HPolygon, ε; upper_bound=upper_bound) + return overapproximate(S, HPolygon, ε) end """ @@ -80,8 +77,6 @@ Return an approximation of a given set as a hyperrectangle. - `S` -- set - `Hyperrectangle` -- type for dispatch -- `upper_bound` -- (optional, default: `false`) use overapproximation in - support function computation? ### Output @@ -89,9 +84,8 @@ A hyperrectangle. """ function overapproximate(S::LazySet, ::Type{<:Hyperrectangle}; - upper_bound::Bool=false )::Union{Hyperrectangle, EmptySet} - box_approximation(S; upper_bound=upper_bound) + box_approximation(S) end """ @@ -100,14 +94,12 @@ end Alias for `overapproximate(S, Hyperrectangle)`. """ overapproximate(S::LazySet; - upper_bound::Bool=false )::Union{Hyperrectangle, EmptySet} = - overapproximate(S, Hyperrectangle; upper_bound=upper_bound) + overapproximate(S, Hyperrectangle) """ overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}}, - ::Type{<:Zonotope}; - [upper_bound]::Bool=false)::Zonotope where {N<:Real} + ::Type{<:Zonotope})::Zonotope where {N<:Real} Overapproximate the convex hull of two zonotopes. @@ -115,7 +107,6 @@ Overapproximate the convex hull of two zonotopes. - `S` -- convex hull of two zonotopes of the same order - `Zonotope` -- type for dispatch -- `upper_bound` -- (optional, default: `false`) currently ignored ### Algorithm @@ -134,8 +125,7 @@ zonotope, which is in general expensive in high dimensions. This is further inve in: *Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on Discrete Algorithms, pp. 803-812*. """ function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}}, - ::Type{<:Zonotope}; - upper_bound::Bool=false)::Zonotope where {N<:Real} + ::Type{<:Zonotope})::Zonotope where {N<:Real} Z1, Z2 = S.X, S.Y @assert order(Z1) == order(Z2) center = (Z1.center+Z2.center)/2 @@ -145,8 +135,7 @@ end """ overapproximate(X::LazySet, - dir::AbstractDirections; - [upper_bound]::Bool=false + dir::AbstractDirections )::HPolytope Overapproximating a set with template directions. @@ -155,8 +144,6 @@ Overapproximating a set with template directions. - `X` -- set - `dir` -- direction representation -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? ### Output @@ -164,22 +151,19 @@ A `HPolytope` overapproximating the set `X` with the directions from `dir`. """ function overapproximate(X::LazySet{N}, dir::AbstractDirections{N}; - upper_bound::Bool=false )::HPolytope{N} where N - ρ_rec = upper_bound ? ρ_upper_bound : ρ halfspaces = Vector{LinearConstraint{N}}() sizehint!(halfspaces, length(dir)) H = HPolytope(halfspaces) for d in dir - addconstraint!(H, LinearConstraint(d, ρ_rec(d, X))) + addconstraint!(H, LinearConstraint(d, ρ_(d, X))) end return H end """ overapproximate(S::LazySet{N}, - ::Type{Interval}; - [upper_bound]::Bool=false + ::Type{Interval} ) where {N<:Real} Return the overapproximation of a real unidimensional set with an interval. @@ -188,15 +172,13 @@ Return the overapproximation of a real unidimensional set with an interval. - `S` -- one-dimensional set - `Interval` -- type for dispatch -- `upper_bound` -- (optional, default: `false`) currently ignored ### Output An interval. """ function overapproximate(S::LazySet{N}, - ::Type{Interval}; - upper_bound::Bool=false + ::Type{Interval} ) where {N<:Real} @assert dim(S) == 1 lo = σ([-one(N)], S)[1] @@ -207,7 +189,6 @@ end """ overapproximate(cap::Intersection{N, <:LazySet, S}, dir::AbstractDirections{N}; - [upper_bound]::Bool=false, kwargs... ) where {N<:Real, S<:AbstractPolytope{N}} @@ -218,8 +199,6 @@ polytope given a set of template directions. - `cap` -- intersection of a compact set and a polytope - `dir` -- template directions -- `upper_bound` -- (optional, default: `false`) use overapproximation in support - function computation? - `kwargs` -- additional arguments that are passed to the support function algorithm @@ -253,7 +232,6 @@ intersection is empty. """ function overapproximate(cap::Intersection{N, <:LazySet, S}, dir::AbstractDirections{N}; - upper_bound::Bool=false, kwargs... ) where {N<:Real, S<:AbstractPolytope{N}} @@ -263,12 +241,11 @@ function overapproximate(cap::Intersection{N, <:LazySet, S}, Hi = constraints_list(P) m = length(Hi) Q = HPolytope{N}() - ρ_rec = upper_bound ? ρ_upper_bound : ρ for di in dir - ρ_X_Hi_min = ρ_rec(di, X ∩ Hi[1], kwargs...) + ρ_X_Hi_min = ρ(di, X ∩ Hi[1], kwargs...) for i in 2:m - ρ_X_Hi = ρ_rec(di, X ∩ Hi[i], kwargs...) + ρ_X_Hi = ρ(di, X ∩ Hi[i], kwargs...) if ρ_X_Hi < ρ_X_Hi_min ρ_X_Hi_min = ρ_X_Hi end diff --git a/src/Intersection.jl b/src/Intersection.jl index a29c5acc0a..fb1ebb9148 100644 --- a/src/Intersection.jl +++ b/src/Intersection.jl @@ -110,40 +110,141 @@ direction. The support vector in the given direction. """ function σ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real} - # TODO document behavior if the direction has norm zero - # TODO error message if the intersection is empty! - # TODO implement - error("the exact support vector of an intersection is not implemented yet") + error("the exact support vector of an intersection is not implemented") end """ - ρ(d::AbstractVector{N}, cap::Intersection{N}; - upper_bound=false, kwargs...) where {N<:Real} + ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real} -Return the support function of the intersection of two convex sets in a given -direction. +Return an upper bound on the support function of the intersection of two convex +sets in a given direction. ### Input -- `d` -- direction -- `cap` -- intersection of two convex sets -- `upper_bound` -- (optional, default: `false`) if `false`, compute the support - function exactly; otherwise use an overapproximative algorithm +- `d` -- direction +- `cap` -- intersection of two convex sets ### Output -The support function in the given direction. +An uper bound on the support function in the given direction. + +### Algorithm + +The support function of an intersection of ``X`` and ``Y`` is upper bounded by +the minimum of the support functions of ``X`` and ``Y``. +""" +function ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real} + return min(ρ(d, cap.X), ρ(d, cap.Y)) +end + +""" + ρ(d::AbstractVector{N}, + cap::Intersection{N, + <:LazySet{N}, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; + [algorithm]::String="line_search", + [check_intersection]::Bool=true, + [kwargs...]) where N<:Real + +Return the support function of the intersection of a compact set and a +half-space/hyperplane/line in a given direction. + +### Input + +- `d` -- direction +- `cap` -- lazy intersection of a compact set and a half-space/hyperplane/ + line +- `algorithm` -- (optional, default: `"line_search"`): the algorithm to + calculate the support function; valid options are: + + * `"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 + support function + +### Output + +The scalar value of the support function of the set `cap` in the given +direction. + +### Notes + +It is assumed that the set `cap.X` is compact. + +The `check_intersection` flag can be useful if it is known in advance that the +intersection is non-empty. + +Any additional number of arguments to the algorithm backend can be passed as +keyword arguments. + +### Algorithm + +The algorithms are based on solving the associated optimization problem + +```math +\\min_\\{ λ ∈ D_h \\} ρ(ℓ - λa, X) + λb. +``` +where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or +``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane. + +For additional information we refer to: + +- [G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with + Support Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809). +- [C. Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous + Dynamics, PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2). +- [T. Rockafellar, R. Wets. + Variational Analysis](https://www.springer.com/us/book/9783540627722). """ -function ρ(d::AbstractVector{N}, cap::Intersection{N}; - upper_bound=false, kwargs...) where {N<:Real} - if upper_bound - return LazySets.Approximations.ρ_upper_bound(d, cap; kwargs...) +function ρ(d::AbstractVector{N}, + cap::Intersection{N, + <:LazySet{N}, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; + algorithm::String="line_search", + check_intersection::Bool=true, + kwargs...) where N<:Real + + X = cap.X # compact set + H = cap.Y # halfspace or hyperplane + + # if the intersection is empty => stop + if check_intersection && is_intersection_empty(X, H, false) + error("the intersection is empty") + end + + if algorithm == "line_search" + @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " * + "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("the exact support function of an intersection is not implemented; " * - "try using `upper_bound=true`") + error("algorithm $(algorithm) unknown") end + return s end +# Symmetric case +ρ(ℓ::AbstractVector{N}, + cap::Intersection{N, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, + <:LazySet{N}}; + algorithm::String="line_search", check_intersection::Bool=true, + kwargs...) where N<:Real = + ρ(ℓ, cap.Y ∩ cap.X; algorithm=algorithm, + check_intersection=check_intersection, kwargs...) + """ ∈(x::AbstractVector{N}, cap::Intersection{N})::Bool where {N<:Real} @@ -316,115 +417,9 @@ function ∈(x::AbstractVector{N}, ia::IntersectionArray{N})::Bool where {N<:Rea return true end -""" - ρ(d::AbstractVector{N}, - cap::Intersection{N, - <:LazySet{N}, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; - [algorithm]::String="line_search", - [check_intersection]::Bool=true, - [kwargs...]) where N<:Real - -Return the support function of the intersection of a compact set and a -half-space/hyperplane/line in a given direction. - -### Input - -- `d` -- direction -- `cap` -- lazy intersection of a compact set and a half-space/hyperplane/ - line -- `algorithm` -- (optional, default: `"line_search"`): the algorithm to - calculate the support function; valid options are: - - * `"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 - support function - -### Output - -The scalar value of the support function of the set `cap` in the given -direction. - -### Notes - -It is assumed that the set `cap.X` is compact. - -The `check_intersection` flag can be useful if it is known in advance that the -intersection is non-empty. - -Any additional number of arguments to the algorithm backend can be passed as -keyword arguments. - -### Algorithm - -The algorithms are based on solving the associated optimization problem - -```math -\\min_\\{ λ ∈ D_h \\} ρ(ℓ - λa, X) + λb. -``` -where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or -``D_h = \\{ λ : λ ∈ \\mathbb{R} \\}`` if ``H`` is a hyperplane. - -For additional information we refer to: - -- [G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with - Support Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809). -- [C. Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous - Dynamics, PhD thesis](https://tel.archives-ouvertes.fr/tel-00422569v2). -- [T. Rockafellar, R. Wets. - Variational Analysis](https://www.springer.com/us/book/9783540627722). -""" -function ρ(d::AbstractVector{N}, - cap::Intersection{N, - <:LazySet{N}, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; - algorithm::String="line_search", - check_intersection::Bool=true, - upper_bound::Bool=false, - kwargs...) where N<:Real - - X = cap.X # compact set - H = cap.Y # halfspace or hyperplane - - # if the intersection is empty => stop - if check_intersection && - is_intersection_empty(X, H, false; upper_bound=upper_bound) - error("the intersection is empty") - end - - if algorithm == "line_search" - @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " * - "the package 'Optim' to be loaded" - (s, _) = _line_search(d, X, H; upper_bound=upper_bound, 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; upper_bound=upper_bound, kwargs...) - else - error("algorithm $(algorithm) unknown") - end - return s -end - -# Symmetric case -ρ(ℓ::AbstractVector{N}, - cap::Intersection{N, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, - <:LazySet{N}}; - algorithm::String="line_search", check_intersection::Bool=true, - kwargs...) where N<:Real = - ρ(ℓ, cap.Y ∩ cap.X; algorithm=algorithm, - check_intersection=check_intersection, kwargs...) +# ================================== +# Algorithms for lazy intersection +# ================================== function load_optim_intersection() return quote @@ -432,7 +427,7 @@ return quote import Optim """ - _line_search(ℓ, X, H; [upper_bound]=false, [kwargs...]) + _line_search(ℓ, X, H; [kwargs...]) Given a compact and convex set ``X`` and a halfspace ``H = \\{x: a^T x ≤ b \\}`` or a hyperplane ``H = \\{x: a^T x = b \\}``, calculate: @@ -448,9 +443,6 @@ where ``D_h = \\{ λ : λ ≥ 0 \\}`` if ``H`` is a half-space or - `ℓ` -- direction - `X` -- set - `H` -- halfspace or hyperplane -- `upper_bound` -- (optional, default: `false`) if `false`, compute the support - function exactly; otherwise use an overapproximative - algorithm ### Output @@ -495,14 +487,12 @@ 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, Line}; - upper_bound::Bool=false, kwargs...) +function _line_search(ℓ, X, H::Union{HalfSpace, Hyperplane, Line}; kwargs...) options = Dict(kwargs) # Initialization a, b = H.a, H.b - ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ - f(λ) = ρ_rec(ℓ - λ[1] * a, X) + λ[1] * b + f(λ) = ρ(ℓ - λ[1] * a, X) + λ[1] * b if haskey(options, :lower) lower = pop!(options, :lower) @@ -541,7 +531,6 @@ end # load_optim [lazy_linear_map]=false, [lazy_2d_intersection]=true, [algorithm_2d_intersection]=nothing, - [upper_bound]::Bool=false, [kwargs...]) where {N} Given a compact and convex set ``X`` and a hyperplane ``H = \\{x: n ⋅ x = γ \\}``, @@ -561,9 +550,6 @@ rank-2 projection ``Π_{nℓ} X`` and the line ``Lγ = \\{(x, y): x = γ \\}``. - `algorithm_2d_intersection` -- (optional, default: `nothing`) if given, fixes the support function algorithm used for the intersection in 2D; otherwise the default is implied -- `upper_bound` -- (optional, default: `false`) if `false`, compute the - support function exactly; otherwise use an - overapproximative algorithm ### Output @@ -597,7 +583,6 @@ function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}}; lazy_linear_map=false, lazy_2d_intersection=true, algorithm_2d_intersection=nothing, - upper_bound::Bool=false, kwargs...) where {N} n = H.a # normal vector to the hyperplane @@ -612,10 +597,9 @@ function _projection(ℓ, X, H::Union{Hyperplane{N}, Line{N}}; Xnℓ = lazy_linear_map ? LinearMap(Πnℓ, X) : linear_map(Πnℓ, X) Xnℓ⋂Lγ = lazy_2d_intersection ? Intersection(Xnℓ, Lγ) : intersection(Xnℓ, Lγ) - ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ if algorithm_2d_intersection == nothing - return ρ_rec(y_dir, Xnℓ⋂Lγ; kwargs...) + return ρ(y_dir, Xnℓ⋂Lγ; kwargs...) else - return ρ_rec(y_dir, Xnℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...) + return ρ(y_dir, Xnℓ⋂Lγ, algorithm=algorithm_2d_intersection; kwargs...) end end diff --git a/src/is_intersection_empty.jl b/src/is_intersection_empty.jl index fb059518f4..8436f9d4de 100644 --- a/src/is_intersection_empty.jl +++ b/src/is_intersection_empty.jl @@ -721,17 +721,15 @@ the half-space: ``σ(-a) ∈ hs``. The support vector is thus also a witness. Optional keyword arguments can be passed to the `ρ` function. In particular, if -`X` is a lazy intersection, the option `upper_bound=true` may be useful. +`X` is a lazy intersection, options can be passed to the line search algorithm. """ function is_intersection_empty(X::LazySet{N}, hs::HalfSpace{N}, witness::Bool=false; - upper_bound::Bool=false, kwargs... )::Union{Bool, Tuple{Bool,Vector{N}}} where N<:Real if !witness - ρ_rec = upper_bound ? LazySets.Approximations.ρ_upper_bound : ρ - return -ρ_rec(-hs.a, X; kwargs...) > hs.b + return -ρ(-hs.a, X; kwargs...) > hs.b end # for witness production, we compute the support vector instead From c3bb22695a3a435d2bb649f7012cd33bf373262b Mon Sep 17 00:00:00 2001 From: mforets Date: Sat, 13 Oct 2018 00:23:07 -0300 Subject: [PATCH 2/6] add lazy intersection w/polytopes --- docs/src/lib/approximations.md | 6 - src/Approximations/Approximations.jl | 1 - src/Approximations/upper_bounds.jl | 241 --------------------------- src/Intersection.jl | 61 +++++++ test/runtests.jl | 1 - 5 files changed, 61 insertions(+), 249 deletions(-) diff --git a/docs/src/lib/approximations.md b/docs/src/lib/approximations.md index 6273600771..85093af9bb 100644 --- a/docs/src/lib/approximations.md +++ b/docs/src/lib/approximations.md @@ -66,9 +66,3 @@ BoxDiagDirections ``` See also `overapproximate(X::LazySet, dir::AbstractDirections)::HPolytope`. - -## Upper bounds - -```@docs -ρ_upper_bound -``` diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 536f2ec3b7..fa0f70c51a 100644 --- a/src/Approximations/Approximations.jl +++ b/src/Approximations/Approximations.jl @@ -37,6 +37,5 @@ include("box_approximations.jl") include("template_directions.jl") include("overapproximate.jl") include("decompositions.jl") -#include("upper_bounds.jl") end # module diff --git a/src/Approximations/upper_bounds.jl b/src/Approximations/upper_bounds.jl index 9d42bbb3af..139597f9cb 100644 --- a/src/Approximations/upper_bounds.jl +++ b/src/Approximations/upper_bounds.jl @@ -1,243 +1,2 @@ -export ρ_upper_bound -""" - ρ_upper_bound(d::AbstractVector{N}, - X::LazySet{N}; - kwargs...) where {N<:Real} -Return an upper bound of the support function of a given set. - -### Input - -- `d` -- direction -- `X` -- set - -### Output - -An upper bound of the support function of the given set. - -### Algorithm - -The default implementation of `ρ_upper_bound` is the exact `ρ(d, X)`. -""" -function ρ_upper_bound(d::AbstractVector{N}, - X::LazySet{N}; - kwargs...) where {N<:Real} - return ρ(d, X) -end - -""" - ρ_upper_bound(d::AbstractVector{N}, - X::LazySet{N}; - kwargs...) where {N<:Real} - -Return an upper bound of the support function of a given set. - -### Input - -- `d` -- direction -- `X` -- set - -### Output - -An upper bound of the support function of the given set. - -### Algorithm - -The default implementation of `ρ_upper_bound` for lazy operations is to pass -all keyword arguments to `ρ(d, X; kwargs...)` if that method exists. -""" -function ρ_upper_bound(d::AbstractVector{N}, - X::Union{LinearMap{N}}; - kwargs...) where {N<:Real} - return ρ(d, X, upper_bound=true, kwargs...) -end - -""" - ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N}; kwargs...) where {N<:Real} - -Return an upper bound of the support function of the intersection of two sets. - -### Input - -- `d` -- direction -- `cap` -- intersection - -### Output - -An upper bound of the support function of the given intersection. - -### Algorithm - -The support function of an intersection of ``X`` and ``Y`` is upper bounded by -the minimum of the support functions of ``X`` and ``Y``. -""" -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N}; kwargs...) where {N<:Real} - return min(ρ_upper_bound(d, cap.X; kwargs...), ρ_upper_bound(d, cap.Y; kwargs...)) -end - -""" - ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}}; - kwargs...) where {N<:Real} - -Return an upper bound of the intersection between a compact set and a -polytope along a given direction. - -### Input - -- `d` -- direction -- `cap` -- intersection of a compact set and a polytope -- `kwargs` -- additional arguments that are passed to the support function algorithm - -### Output - -An upper bound of the support function of the given intersection. - -### Algorithm - -The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for each -half-space in the set `P` and then take the minimum. This gives an overapproximation -of the exact support function. - -This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection -for Reachability Computations with Support -Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809). - -### Notes - -This method relies on having available the `constraints_list` of the polytope -`P`. - -This method of overapproximation can return a non-empty set even if the original -intersection is empty. -""" -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}}; - kwargs...) where {N<:Real} - - X = cap.X # compact set - P = cap.Y # polytope - return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) -end - -# disambiguation -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, <:AbstractPolytope{N}, <:AbstractPolytope{N}}; - kwargs...) where {N<:Real} - X = cap.X # compact set - P = cap.Y # polytope - return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) -end - -# symmetric function -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, <:AbstractPolytope{N}, <:LazySet{N}}; - kwargs...) where {N<:Real} - return ρ_upper_bound(d, cap.Y ∩ cap.X; kwargs...) -end - -""" - ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, <:LazySet{N}, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; - [algorithm]::String="line_search", - [check_intersection]::Bool=true, - [kwargs...]) where N<:Real - -Return an upper bound of the support function of the intersection of a compact -set and a half-space or a hyperplane or a line in a given direction. - -### Input - -- `d` -- direction -- `cap` -- lazy intersection of a compact set and a half-space/hyperplane/ - line -- `algorithm` -- (optional, default: `"line_search"`): the algorithm to - calculate the support function; valid options are: - - * `"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 - support function - -### Output - -The scalar value of the support function of the set `cap` in the given -direction. - -### Notes - -It is assumed that the set `cap.X` is compact. - -The `check_intersection` flag can be useful if it is known in advance that the -intersection is non-empty. - -Any additional number of arguments to the algorithm backend can be passed as -keyword arguments. - -### Algorithm - -We fall back to the implementation of `ρ` with the same arguments. -""" -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, - <:LazySet{N}, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; - algorithm::String="line_search", - check_intersection::Bool=true, - kwargs... - ) where N<:Real - return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection, - upper_bound=true, kwargs...) -end - -# symmetric case -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, - <:LazySet{N}}; - algorithm::String="line_search", - check_intersection::Bool=true, - kwargs... - ) where N<:Real - return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection, - upper_bound=true, kwargs...) -end - -# disambiguation -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, - <:AbstractPolytope{N}}; - algorithm::String="line_search", - check_intersection::Bool=true, - kwargs... - ) where N<:Real - return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection, - upper_bound=true, kwargs...) -end - -# disambiguation symmetric case -function ρ_upper_bound(d::AbstractVector{N}, - cap::Intersection{N, - <:AbstractPolytope{N}, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; - algorithm::String="line_search", - check_intersection::Bool=true, - kwargs... - ) where N<:Real - return ρ(d, cap; algorithm=algorithm, check_intersection=check_intersection, - upper_bound=true, kwargs...) -end diff --git a/src/Intersection.jl b/src/Intersection.jl index fb1ebb9148..27a80d8547 100644 --- a/src/Intersection.jl +++ b/src/Intersection.jl @@ -245,6 +245,67 @@ end ρ(ℓ, cap.Y ∩ cap.X; algorithm=algorithm, check_intersection=check_intersection, kwargs...) +""" + ρ(d::AbstractVector{N}, + cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}}; + kwargs...) where {N<:Real} + +Return an upper bound of the intersection between a compact set and a +polytope along a given direction. + +### Input + +- `d` -- direction +- `cap` -- intersection of a compact set and a polytope +- `kwargs` -- additional arguments that are passed to the support function algorithm + +### Output + +An upper bound of the support function of the given intersection. + +### Algorithm + +The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for each +half-space in the set `P` and then take the minimum. This gives an overapproximation +of the exact support function. + +This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection +for Reachability Computations with Support +Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809). + +### Notes + +This method relies on having available the `constraints_list` of the polytope +`P`. + +This method of overapproximation can return a non-empty set even if the original +intersection is empty. +""" +function ρ(d::AbstractVector{N}, + cap::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}}; + kwargs...) where {N<:Real} + + X = cap.X # compact set + P = cap.Y # polytope + return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) +end + +# symmetric function +function ρ(d::AbstractVector{N}, + cap::Intersection{N, <:AbstractPolytope{N}, <:LazySet{N}}; + kwargs...) where {N<:Real} + return ρ(d, cap.Y ∩ cap.X; kwargs...) +end + +# disambiguation +function ρ(d::AbstractVector{N}, + cap::Intersection{N, <:AbstractPolytope{N}, <:AbstractPolytope{N}}; + kwargs...) where {N<:Real} + X = cap.X # compact set + P = cap.Y # polytope + return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) +end + """ ∈(x::AbstractVector{N}, cap::Intersection{N})::Bool where {N<:Real} diff --git a/test/runtests.jl b/test/runtests.jl index becffaf7c1..9e51174df5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -96,7 +96,6 @@ if test_suite_basic @time @testset "LazySets.Approximations.ballinf_approximation" begin include("unit_ballinf_approximation.jl") end @time @testset "LazySets.Approximations.radiusdiameter" begin include("unit_radiusdiameter.jl") end @time @testset "LazySets.Approximations.decompose" begin include("unit_decompose.jl") end - @time @testset "LazySets.Approximations.upper_bounds" begin include("unit_upper_bounds.jl") end # ======================== # Testing method ambiguity From 0bc04abbdc72950de9f589c78e96920017aabe4b Mon Sep 17 00:00:00 2001 From: mforets Date: Sat, 13 Oct 2018 00:49:40 -0300 Subject: [PATCH 3/6] fix _ in function name --- src/Approximations/overapproximate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 1e5a17c12b..02049a5051 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -156,7 +156,7 @@ function overapproximate(X::LazySet{N}, sizehint!(halfspaces, length(dir)) H = HPolytope(halfspaces) for d in dir - addconstraint!(H, LinearConstraint(d, ρ_(d, X))) + addconstraint!(H, LinearConstraint(d, ρ(d, X))) end return H end From 3c86327177edbb0f50edce0cf7f60258213a338b Mon Sep 17 00:00:00 2001 From: mforets Date: Sat, 13 Oct 2018 16:25:33 -0300 Subject: [PATCH 4/6] update test of its --- test/unit_plot.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/unit_plot.jl b/test/unit_plot.jl index bc1f848865..634685c931 100644 --- a/test/unit_plot.jl +++ b/test/unit_plot.jl @@ -47,8 +47,11 @@ for N in [Float64, Rational{Int}, Float32] # binary set operations ch = ConvexHull(b1, bi) cha = ConvexHullArray([b1, bi]) + + using Optim its = Intersection(b1, bi) itsa = IntersectionArray([b1, bi]) + ms = MinkowskiSum(b1, bi) msa = MinkowskiSumArray([b1, bi]) cms = CacheMinkowskiSum([b1, bi]) @@ -79,7 +82,7 @@ for N in [Float64, Rational{Int}, Float32] @test_throws ErrorException plot(l) # TODO see #576 plot(ch) plot(cha) - @test_throws ErrorException plot(its) # TODO not implemented yet + plot(its) @test_throws ErrorException plot(itsa) # TODO not implemented yet plot(sih) plot(lm) From 0cd55be356eb7f9ee0b23eedcd4e1da12db1efab Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 13 Oct 2018 22:58:29 +0200 Subject: [PATCH 5/6] move tests (no Rational support) --- test/unit_plot.jl | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/test/unit_plot.jl b/test/unit_plot.jl index 634685c931..59c109462a 100644 --- a/test/unit_plot.jl +++ b/test/unit_plot.jl @@ -1,4 +1,4 @@ -using Plots +using Plots, Optim for N in [Float64, Rational{Int}, Float32] p0 = zero(N) @@ -48,10 +48,6 @@ for N in [Float64, Rational{Int}, Float32] ch = ConvexHull(b1, bi) cha = ConvexHullArray([b1, bi]) - using Optim - its = Intersection(b1, bi) - itsa = IntersectionArray([b1, bi]) - ms = MinkowskiSum(b1, bi) msa = MinkowskiSumArray([b1, bi]) cms = CacheMinkowskiSum([b1, bi]) @@ -82,8 +78,6 @@ for N in [Float64, Rational{Int}, Float32] @test_throws ErrorException plot(l) # TODO see #576 plot(ch) plot(cha) - plot(its) - @test_throws ErrorException plot(itsa) # TODO not implemented yet plot(sih) plot(lm) plot(ms) @@ -112,8 +106,6 @@ for N in [Float64, Rational{Int}, Float32] @test_throws ErrorException plot(hs, ε) # TODO see #576/#578 @test_throws ErrorException plot(hp, ε) # TODO see #576/#578 @test_throws ErrorException plot(l, ε) # TODO see #576/#578 - @test_throws ErrorException plot(its, ε) # TODO not implemented yet - @test_throws ErrorException plot(itsa, ε) # TODO not implemented yet plot(ch, ε) plot(cha, ε) plot(sih, ε) @@ -172,3 +164,28 @@ for N in [Float64, Float32] @test_throws ErrorException plot(b2, ε) # TODO see #578 end end + +# set types that only work with Float64 +for N in [Float64] + v0 = zeros(N, 2) + p1 = one(N) + + b1 = Ball1(v0, p1) + bi = BallInf(v0, p1) + + # binary set operations + its = Intersection(b1, bi) + itsa = IntersectionArray([b1, bi]) + + # ----- + # plots + # ----- + + plot(its) + @test_throws ErrorException plot(itsa) # TODO not implemented yet + + # ε-close + ε = N(1e-2) + plot(its) + @test_throws ErrorException plot(itsa, ε) # TODO not implemented yet +end From 67327e7e3892bb06cdffaa8e735fbfec47305895 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 13 Oct 2018 23:30:21 +0200 Subject: [PATCH 6/6] disambiguation --- src/Intersection.jl | 92 ++++++++++++++++++++++++++++++--------------- 1 file changed, 61 insertions(+), 31 deletions(-) diff --git a/src/Intersection.jl b/src/Intersection.jl index 27a80d8547..3bb0d39620 100644 --- a/src/Intersection.jl +++ b/src/Intersection.jl @@ -137,6 +137,33 @@ function ρ(d::AbstractVector{N}, cap::Intersection{N}) where {N<:Real} return min(ρ(d, cap.X), ρ(d, cap.Y)) end +function ρ_helper(d::AbstractVector{N}, + cap::Intersection{N}, + algorithm::String, + check_intersection::Bool; + kwargs...) where N<:Real + X = cap.X # compact set + H = cap.Y # halfspace or hyperplane or line + + # if the intersection is empty => stop + if check_intersection && is_intersection_empty(X, H, false) + error("the intersection is empty") + end + + if algorithm == "line_search" + @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " * + "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 + """ ρ(d::AbstractVector{N}, cap::Intersection{N, @@ -212,38 +239,19 @@ function ρ(d::AbstractVector{N}, algorithm::String="line_search", check_intersection::Bool=true, kwargs...) where N<:Real - - X = cap.X # compact set - H = cap.Y # halfspace or hyperplane - - # if the intersection is empty => stop - if check_intersection && is_intersection_empty(X, H, false) - error("the intersection is empty") - end - - if algorithm == "line_search" - @assert isdefined(Main, :Optim) "the algorithm $algorithm needs " * - "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 + return ρ_helper(d, cap, algorithm, check_intersection, kwargs...) end -# Symmetric case -ρ(ℓ::AbstractVector{N}, - cap::Intersection{N, - <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, - <:LazySet{N}}; - algorithm::String="line_search", check_intersection::Bool=true, - kwargs...) where N<:Real = - ρ(ℓ, cap.Y ∩ cap.X; algorithm=algorithm, - check_intersection=check_intersection, kwargs...) +# symmetric method +function ρ(ℓ::AbstractVector{N}, + cap::Intersection{N, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, + <:LazySet{N}}; + algorithm::String="line_search", + check_intersection::Bool=true, + kwargs...) where N<:Real + return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection, kwargs...) +end """ ρ(d::AbstractVector{N}, @@ -290,7 +298,7 @@ function ρ(d::AbstractVector{N}, return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) end -# symmetric function +# symmetric method function ρ(d::AbstractVector{N}, cap::Intersection{N, <:AbstractPolytope{N}, <:LazySet{N}}; kwargs...) where {N<:Real} @@ -306,6 +314,28 @@ function ρ(d::AbstractVector{N}, return minimum([ρ(d, X ∩ Hi; kwargs...) for Hi in constraints_list(P)]) end +# disambiguation +function ρ(d::AbstractVector{N}, + cap::Intersection{N, + <:AbstractPolytope{N}, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}; + algorithm::String="line_search", + check_intersection::Bool=true, + kwargs...) where N<:Real + return ρ_helper(d, cap, algorithm, check_intersection; kwargs...) +end + +# symmetric method +function ρ(d::AbstractVector{N}, + cap::Intersection{N, + <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}, + <:AbstractPolytope{N}}; + algorithm::String="line_search", + check_intersection::Bool=true, + kwargs...) where N<:Real + return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection; kwargs...) +end + """ ∈(x::AbstractVector{N}, cap::Intersection{N})::Bool where {N<:Real}