diff --git a/Project.toml b/Project.toml index 6d22c0290c..b8888fc29d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LazySets" uuid = "b4f0291d-fe17-52bc-9479-3d1a343d9043" -version = "1.12.0" +version = "1.13.0" [deps] Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -26,6 +26,7 @@ Optim = "≥ 0.14.1" Polyhedra = "≥ 0.3.4" RecipesBase = "≥ 0.3.0" Requires = "≥ 0.4.0" +TaylorModels = "≥ 0.1.1" [extras] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" @@ -40,6 +41,9 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" [targets] -test = ["Expokit", "Pkg", "IntervalArithmetic", "GLPKMathProgInterface", "RecipesBase", "MathProgBase", "Requires", "Optim", "Polyhedra", "Documenter", "Plots", "GR"] +test = ["Documenter", "Expokit", "GLPKMathProgInterface", "GR", + "IntervalArithmetic", "MathProgBase", "Optim", "Pkg", "Plots", + "Polyhedra", "RecipesBase", "Requires", "TaylorModels"] diff --git a/README.md b/README.md index 01b38a58b6..2188b80857 100644 --- a/README.md +++ b/README.md @@ -12,25 +12,36 @@ ## Resources - [Manual](http://juliareach.github.io/LazySets.jl/latest/) -- [Contributing](https://juliareach.github.io/LazySets.jl/latest/about.html#Contributing-1) +- [Contributing](https://juliareach.github.io/LazySets.jl/latest/about/#Contributing-1) - [Release notes of tagged versions](https://github.com/JuliaReach/LazySets.jl/releases) - [Release notes of the development version](https://github.com/JuliaReach/LazySets.jl/wiki/Release-log-tracker) - [Publications](https://juliareach.github.io/Reachability.jl/latest/publications.html) +- [Developers](https://juliareach.github.io/LazySets.jl/latest/about/#Credits-1) ## Installing -This package requires Julia v0.6 or later. Refer to the [official documentation](https://julialang.org/downloads) -on how to install and run Julia in your system. +This package requires Julia v1.0 or later. +Refer to the [official documentation](https://julialang.org/downloads) on how to +install and run Julia on your system. -To install the latest release of this package, use the following command inside Julia's REPL: +Depending on your needs, choose an appropriate command from the following list +and enter it in Julia's REPL. +To activate the `pkg` mode, type `]` (and to leave it, type `exit()`). + +#### [Install the latest release version](https://julialang.github.io/Pkg.jl/v1/managing-packages/#Adding-registered-packages-1) + +```julia +pkg> add LazySets +``` + +#### Install the latest development version ```julia -using Pkg -Pkg.add("LazySets") +pkg> add LazySets#master ``` -If you want to install the latest development version, do: +#### [Clone the package for development](https://julialang.github.io/Pkg.jl/v1/managing-packages/#Developing-packages-1) ```julia -Pkg.clone("https://github.com/JuliaReach/LazySets.jl.git") +pkg> dev LazySets ``` diff --git a/docs/Project.toml b/docs/Project.toml index 76eaed3858..de86e12157 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,6 +10,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" [compat] BenchmarkTools = "≥ 0.4.2" @@ -23,3 +24,4 @@ Plots = "≥ 0.25.1" Polyhedra = "≥ 0.3.4" RecipesBase = "≥ 0.3.0" StaticArrays = "≥ 0.10.3" +TaylorModels = "≥ 0.1.1" diff --git a/docs/src/about.md b/docs/src/about.md index 34ae0810e2..194fdee83f 100644 --- a/docs/src/about.md +++ b/docs/src/about.md @@ -111,9 +111,9 @@ Here we list the names of the maintainers of the `LazySets.jl` library, as well - Kostiantyn Potomkin - Frédéric Viry -### Colaborators +### Acknowledgements -We are also grateful to the following persons for enlightening discussions: +We are grateful to the following persons for enlightening discussions: - [Sergiy Bogomolov](https://www.sergiybogomolov.com/) - [Goran Frehse](https://sites.google.com/site/frehseg/) diff --git a/docs/src/lib/binary_functions.md b/docs/src/lib/binary_functions.md index 05e3c5bd51..06c3f972e5 100644 --- a/docs/src/lib/binary_functions.md +++ b/docs/src/lib/binary_functions.md @@ -37,6 +37,8 @@ is_intersection_empty(::Universe{N}, ::LazySet{N}, ::Bool=false) where {N<:Real} is_intersection_empty(::Complement{N}, ::LazySet{N}, ::Bool=false) where {N<:Real} is_intersection_empty(::Zonotope{N}, ::Zonotope{N}, ::Bool=false) where {N<:Real} is_intersection_empty(::Interval{N}, ::Interval{N}, ::Bool=false) where {N<:Real} +is_intersection_empty(X::CartesianProductArray{N}, Y::HPolyhedron{N}) where {N<:Real} +is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N} ``` ## Convex hull @@ -64,6 +66,7 @@ intersection(::UnionSet{N}, ::LazySet{N}) where {N<:Real} intersection(::UnionSetArray{N}, ::LazySet{N}) where {N<:Real} intersection(::Universe{N}, ::LazySet{N}) where {N<:Real} intersection(::AbstractPolyhedron{N}, ::ResetMap{N}) where {N<:Real} +intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N<:Real} ``` ## Subset check diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md index afa14eae19..64155c7b0d 100644 --- a/docs/src/lib/interfaces.md +++ b/docs/src/lib/interfaces.md @@ -76,6 +76,7 @@ RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractVector{VN}, ::N=N(1e-3), ```@docs ==(::LazySet, ::LazySet) +≈(::LazySet, ::LazySet) copy(::LazySet) ``` diff --git a/docs/src/lib/operations.md b/docs/src/lib/operations.md index 8498512ed9..17a57e8ae2 100644 --- a/docs/src/lib/operations.md +++ b/docs/src/lib/operations.md @@ -52,6 +52,10 @@ isempty(::CartesianProductArray) constraints_list(::CartesianProductArray{N}) where {N<:Real} vertices_list(::CartesianProductArray{N}) where {N<:Real} array(::CartesianProductArray{N, S}) where {N<:Real, S<:LazySet{N}} +block_structure(cpa::CartesianProductArray{N}) where {N} +block_to_dimension_indices(cpa::CartesianProductArray{N}, vars::Vector{Int}) where {N} +substitute_blocks(low_dim_cpa::CartesianProductArray{N}, orig_cpa::CartesianProductArray{N}, +blocks::Vector{Tuple{Int,Int}}) where {N} ``` Inherited from [`LazySet`](@ref): * [`norm`](@ref norm(::LazySet, ::Real)) diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index 04f07b9e24..1a12ecf0a8 100644 --- a/docs/src/lib/utils.md +++ b/docs/src/lib/utils.md @@ -35,6 +35,7 @@ same_block_structure substitute substitute! σ_helper +get_constrained_lowdimset @neutral @absorbing @neutral_absorbing diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 76bd81697c..023279fb2d 100644 --- a/src/Approximations/Approximations.jl +++ b/src/Approximations/Approximations.jl @@ -34,5 +34,6 @@ include("box_approximations.jl") include("template_directions.jl") include("overapproximate.jl") include("decompositions.jl") +include("init.jl") end # module diff --git a/src/Approximations/init.jl b/src/Approximations/init.jl new file mode 100644 index 0000000000..335002246b --- /dev/null +++ b/src/Approximations/init.jl @@ -0,0 +1,7 @@ +function __init__() + @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" load_taylormodels() +end + +function load_taylormodels() + eval(load_taylormodels_overapproximation()) +end diff --git a/src/Approximations/iterative_refinement.jl b/src/Approximations/iterative_refinement.jl index 22e7487d5e..e466c29c85 100644 --- a/src/Approximations/iterative_refinement.jl +++ b/src/Approximations/iterative_refinement.jl @@ -67,7 +67,7 @@ end """ new_approx(S::LazySet, p1::Vector{N}, d1::Vector{N}, p2::Vector{N}, - d2::Vector{N}) where {N<:Real} + d2::Vector{N}) where {N<:AbstractFloat} Create a `LocalApproximation` instance for the given excerpt of a polygonal approximation. @@ -85,7 +85,7 @@ approximation. A local approximation of `S` in the given directions. """ function new_approx(S::LazySet, p1::Vector{N}, d1::Vector{N}, p2::Vector{N}, - d2::Vector{N}) where {N<:Real} + d2::Vector{N}) where {N<:AbstractFloat} if norm(p1-p2, 2) <= TOL(N) # this approximation cannot be refined and we set q = p1 by convention ap = LocalApproximation{N}(p1, d1, p2, d2, p1, false, zero(N)) @@ -189,7 +189,7 @@ end """ approximate(S::LazySet{N}, - ε::N)::PolygonalOverapproximation{N} where {N<:Real} + ε::N)::PolygonalOverapproximation{N} where {N<:AbstractFloat} Return an ε-close approximation of the given 2D convex set (in terms of Hausdorff distance) as an inner and an outer approximation composed by sorted @@ -204,9 +204,8 @@ local `Approximation2D`. An ε-close approximation of the given 2D convex set. """ -function approximate(S::LazySet{N}, - ε::N)::PolygonalOverapproximation{N} where {N<:Real} - +function approximate(S::LazySet{N}, ε::N + )::PolygonalOverapproximation{N} where {N<:AbstractFloat} # initialize box directions pe = σ(DIR_EAST(N), S) pn = σ(DIR_NORTH(N), S) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 9f54aa77f3..de10114c98 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -1,3 +1,4 @@ +using LazySets: block_to_dimension_indices, substitute_blocks, get_constrained_lowdimset """ overapproximate(X::S, ::Type{S}) where {S<:LazySet} @@ -27,9 +28,9 @@ Otherwise the result is an ε-close approximation as a polygon. ### Input -- `S` -- convex set, assumed to be two-dimensional -- `HPolygon` -- type for dispatch -- `ε` -- (optional, default: `Inf`) error bound +- `S` -- convex set, assumed to be two-dimensional +- `HPolygon` -- type for dispatch +- `ε` -- (optional, default: `Inf`) error bound ### Output @@ -159,7 +160,7 @@ A hyperrectangle. If `c` and `r` denote the center and vector radius of a hyperrectangle `H`, a tight hyperrectangular overapproximation of `M * H` is obtained by transforming `c ↦ M*c` and `r ↦ abs.(M) * c`, where `abs.(⋅)` denotes the element-wise absolute -value operator. +value operator. """ function overapproximate(lm::LinearMap{N, <:AbstractHyperrectangle{N}}, ::Type{Hyperrectangle}) where {N<:Real} @@ -220,7 +221,7 @@ further investigated in [2]. function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}}, ::Type{<:Zonotope})::Zonotope where {N<:Real} Z1, Z2 = S.X, S.Y - + # reduce to the same order if possible if order(Z1) != order(Z2) min_order = min(order(Z1), order(Z2)) @@ -241,7 +242,7 @@ function _overapproximate_convex_hull_zonotope(Z1::Zonotope{N}, Z2::Zonotope{N}) c = (Z1.center + Z2.center)/N(2) # the case of equal order is treated separately to avoid a slicing (this creates a copy) - if order(Z1) == order(Z2) + if order(Z1) == order(Z2) G = hcat(Z1.generators .+ Z2.generators, Z1.center - Z2.center, Z1.generators .- Z2.generators)/N(2) @@ -334,7 +335,7 @@ function overapproximate(X::LazySet{N}, end """ - overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real} + overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real} Return the overapproximation of a real unidimensional set with an interval. @@ -352,11 +353,10 @@ An interval. The method relies on the exact conversion to `Interval`. Two support function evaluations are needed in general. """ -function overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real} +function overapproximate(S::LazySet{N}, ::Type{<:Interval}) where {N<:Real} @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`" return convert(Interval, S) end - function overapproximate_cap_helper(X::LazySet{N}, # compact set P::AbstractPolyhedron{N}, # polyhedron dir::AbstractDirections{N}; @@ -516,3 +516,456 @@ function overapproximate(cap::Intersection{N, return overapproximate(swap(cap), dir; kwargs...) end +# ========================================== +# Functionality that requires TaylorModels +# ========================================== + +# function to be loaded by Requires +function load_taylormodels_overapproximation() + +return quote + +using .TaylorModels: Taylor1, TaylorN, TaylorModelN, TaylorModel1, + polynomial, remainder, domain, + normalize_taylor, linear_polynomial, + constant_term, evaluate, mid + +# helper functions +@inline get_linear_coeffs(p::Taylor1) = linear_polynomial(p).coeffs[2:end] +@inline get_linear_coeffs(p::TaylorN) = linear_polynomial(p).coeffs[2].coeffs + +""" + overapproximate(vTM::Vector{TaylorModel1{T, S}}, + ::Type{Zonotope}) where {T, S} + +Overapproximate a taylor model in one variable with a zonotope. + +### Input + +- `vTM` -- `TaylorModel1` +- `Zonotope` -- type for dispatch + +### Output + +A zonotope that overapproximates the range of the given taylor model. + +### Examples + +If the polynomials are linear, this functions exactly transforms to a zonotope. +However, the nonlinear case necessarily introduces overapproximation error. +Consider the linear case first: + +```julia +julia> using LazySets, TaylorModels + +julia> const IA = IntervalArithmetic; + +julia> I = IA.Interval(-0.5, 0.5) # interval remainder +[-0.5, 0.5] + +julia> x₀ = IA.Interval(0.0) # expansion point +[0, 0] + +julia> D = IA.Interval(-3.0, 1.0) +[-3, 1] + +julia> p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial + 2.0 + 1.0 t + 𝒪(t³) + +julia> p2 = Taylor1([0.9, 3.0], 2) # define another linear polynomial + 0.9 + 3.0 t + 𝒪(t³) + +julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2]] +2-element Array{TaylorModel1{Float64,Float64},1}: + 2.0 + 1.0 t + [-0.5, 0.5] + 0.9 + 3.0 t + [-0.5, 0.5] +``` + +Here, `vTM` is a taylor model vector, since each component is a taylor model in +one variable (`TaylorModel1`). Using `overapproximate(vTM, Zonotope)` we can +compute its associated zonotope in generator representation: + +```julia +julia> using LazySets.Approximations + +julia> Z = overapproximate(vTM, Zonotope); + +julia> center(Z) +2-element Array{Float64,1}: + 1.0 + -2.1 + +julia> Matrix(generators(Z)) +2×3 Array{Float64,2}: + 2.0 0.5 0.0 + 6.0 0.0 0.5 +``` + +Note how the generators of this zonotope mainly consist of two pieces: one comes +from the linear part of the polynomials, and another one that corresponds to the +interval remainder. This conversion gives the same upper and lower bounds as the +range evaluation using interval arithmetic: + +```julia +julia> X = box_approximation(Z) +Hyperrectangle{Float64}([1.0, -2.1], [2.5, 6.5]) + +julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom) +[-1.5, 3.5] × [-8.60001, 4.40001] + +julia> H = convert(Hyperrectangle, Y) # this IntevalBox is the same as X +Hyperrectangle{Float64}([1.0, -2.1], [2.5, 6.5]) +``` +However, the zonotope returns better results if we want to approximate the `TM`, +since it is not axis-aligned: + +```julia +julia> d = [-0.35, 0.93]; + +julia> ρ(d, Z) < ρ(d, X) +true +``` + +This function also works if the polynomials are non-linear; for example suppose +that we add a third polynomial with a quadratic term: + +```julia +julia> p3 = Taylor1([0.9, 3.0, 1.0], 3); + +julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2, p3]] +3-element Array{TaylorModel1{Float64,Float64},1}: + 2.0 + 1.0 t + [-0.5, 0.5] + 0.9 + 3.0 t + [-0.5, 0.5] + 0.9 + 3.0 t + 1.0 t² + [-0.5, 0.5] + +julia> Z = overapproximate(vTM, Zonotope); + +julia> center(Z) +3-element Array{Float64,1}: + 1.0 + -2.1 + 0.8999999999999999 + +julia> Matrix(generators(Z)) +3×4 Array{Float64,2}: + 2.0 0.5 0.0 0.0 + 6.0 0.0 0.5 0.0 + 6.0 0.0 0.0 6.5 +``` +The fourth and last generator corresponds to the addition between the interval +remainder and the box overapproximation of the nonlinear part of `p3` over the +domain. + +### Algorithm + +Let ``\\text{vTM} = (p, I)`` be a vector of ``m`` taylor models, where ``I`` +is the interval remainder in ``\\mathbb{R}^m``. Let ``p_{lin}`` +(resp. ``p_{nonlin}``) correspond to the linear (resp. nonlinear) part of each +scalar polynomial. + +The range of ``\\text{vTM}`` can be enclosed by a zonotope with center ``c`` +and matrix of generators ``G``, ``Z = ⟨c, G⟩``, by performing a conservative +linearization of ``\\text{vTM}``: + +```math + vTM' = (p', I') := (p_{lin} − p_{nonlin} , I + \\text{Int}(p_{nonlin})). +``` + +This algorithm proceeds in two steps: + +1- Conservatively linearize ``\\text{vTM}`` as above and compute a box + overapproximation of the nonlinear part. +2- Transform the linear taylor model to a zonotope exactly through variable + normalization onto the symmetric intervals ``[-1, 1]``. +""" +function overapproximate(vTM::Vector{TaylorModel1{T, S}}, + ::Type{Zonotope}) where {T, S} + m = length(vTM) + + # preallocations + c = Vector{T}(undef, m) # center of the zonotope + gen_lin = Matrix{T}(undef, m, 1) # generator of the linear part + gen_rem = Vector{T}(undef, m) # generators of the remainder + + # compute overapproximation + return _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem) +end + +""" + overapproximate(vTM::Vector{TaylorModelN{N, T, S}}, + ::Type{Zonotope}) where {N,T, S} + + +Overapproximate a multivariate taylor model with a zonotope. + +### Input + +- `vTM` -- `TaylorModelN` +- `Zonotope` -- type for dispatch + +### Output + +A zonotope that overapproximates the range of the given taylor model. + +### Examples + +Consider a vector of two 2-dimensional taylor models of order 2 and 4 +respectively. + +```julia +julia> using LazySets, LazySets.Approximations, TaylorModels + +julia> const IA = IntervalArithmetic; + +julia> x₁, x₂ = set_variables(Float64, ["x₁", "x₂"], order=8) +2-element Array{TaylorN{Float64},1}: + 1.0 x₁ + 𝒪(‖x‖⁹) + 1.0 x₂ + 𝒪(‖x‖⁹) + +julia> x₀ = IntervalBox(0..0, 2) # expansion point +[0, 0] × [0, 0] + +julia> Dx₁ = IA.Interval(0.0, 3.0) # domain for x₁ +[0, 3] + +julia> Dx₂ = IA.Interval(-1.0, 1.0) # domain for x₂ +[-1, 1] + +julia> D = Dx₁ × Dx₂ # take the cartesian product of the domain on each variable +[0, 3] × [-1, 1] + +julia> r = IA.Interval(-0.5, 0.5) # interval remainder +[-0.5, 0.5] + +julia> p1 = 1 + x₁^2 - x₂ + 1.0 - 1.0 x₂ + 1.0 x₁² + 𝒪(‖x‖⁹) + +julia> p2 = x₂^3 + 3x₁^4 + x₁ + 1 + 1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + 𝒪(‖x‖⁹) + +julia> vTM = [TaylorModelN(pi, r, x₀, D) for pi in [p1, p2]] +2-element Array{TaylorModelN{2,Float64,Float64},1}: + 1.0 - 1.0 x₂ + 1.0 x₁² + [-0.5, 0.5] + 1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + [-0.5, 0.5] + +julia> Z = overapproximate(vTM, Zonotope); + +julia> center(Z) +2-element Array{Float64,1}: + 5.5 + 124.0 + +julia> Matrix(generators(Z)) +2×4 Array{Float64,2}: + 0.0 -1.0 5.0 0.0 + 1.5 0.0 0.0 123.0 +``` + +### Algorithm + +We refer to the algorithm description for the univariate case. +""" +function overapproximate(vTM::Vector{TaylorModelN{N, T, S}}, + ::Type{Zonotope}) where {N,T, S} + m = length(vTM) + n = N # number of variables is get_numvars() in TaylorSeries + + # preallocations + c = Vector{T}(undef, m) # center of the zonotope + gen_lin = Matrix{T}(undef, m, n) # generator of the linear part + gen_rem = Vector{T}(undef, m) # generators for the remainder + + # compute overapproximation + return _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem) +end + +function _overapproximate_vTM_zonotope!(vTM, c, gen_lin, gen_rem) + @inbounds for (i, x) in enumerate(vTM) + xpol, xdom = polynomial(x), domain(x) + + # linearize the TM + pol_lin = constant_term(xpol) + linear_polynomial(xpol) + rem_nonlin = remainder(x) + + # build an overapproximation of the nonlinear terms + pol_nonlin = xpol - pol_lin + rem_nonlin += evaluate(pol_nonlin, xdom) + + # normalize the linear polynomial to the symmetric interval [-1, 1] + Q = normalize_taylor(pol_lin, xdom, true) + + # build the generators + α = mid(rem_nonlin) + c[i] = constant_term(Q) + α # constant terms + gen_lin[i, :] = get_linear_coeffs(Q) # linear terms + gen_rem[i] = abs(rem_nonlin.hi - α) + end + return Zonotope(c, hcat(gen_lin, Diagonal(gen_rem))) +end + +end # quote +end # load_taylormodels_overapproximation + +# ========================================== +# Lazy linear maps of cartesian products +# ========================================== + +""" + overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{CartesianProductArray{N, S}} + ) where {N, S<:LazySet{N}} + +Decompose a lazy linear map of a cartesian product array while keeping the +original block structure. + +### Input + +- `lm` -- lazy linear map of cartesian product array +- `CartesianProductArray` -- type for dispatch + +### Output + +A `CartesianProductArray` representing the decomposed linear map. +""" +function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{CartesianProductArray{N, S}} + ) where {N, S<:LazySet{N}} + cpa = array(lm.X) + arr = Vector{S}(undef, length(cpa)) + return _overapproximate_lm_cpa!(arr, lm.M, cpa, S) +end + +""" + overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{<:CartesianProductArray}, + dir::Type{<:AbstractDirections}) where {N} + +Decompose a lazy linear map of a cartesian product array with template +directions while keeping the original block structure. + +### Input + +- `lm` -- lazy linear map of a cartesian product array +- `CartesianProductArray` -- type for dispatch +- `dir` -- template directions for overapproximation + +### Output + +A `CartesianProductArray` representing the decomposed linear map. +""" +function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{<:CartesianProductArray}, + dir::Type{<:AbstractDirections}) where {N} + cpa = array(lm.X) + arr = Vector{HPolytope{N}}(undef, length(cpa)) + return _overapproximate_lm_cpa!(arr, lm.M, cpa, dir) +end + +""" + overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{<:CartesianProductArray}, + set_type::Type{<:LazySet}) where {N} + +Decompose a lazy linear map of a cartesian product array with a given set type +while keeping the original block structure. + +### Input + +- `lm` -- lazy linear map of a cartesian product array +- `CartesianProductArray` -- type for dispatch +- `set_type` -- set type for overapproximation + +### Output + +A `CartesianProductArray` representing the decomposed linear map. +""" +function overapproximate(lm::LinearMap{N, <:CartesianProductArray{N}}, + ::Type{<:CartesianProductArray}, + set_type::Type{<:LazySet}) where {N} + cpa = array(lm.X) + arr = Vector{set_type{N}}(undef, length(cpa)) + return _overapproximate_lm_cpa!(arr, lm.M, cpa, set_type) +end + +function _overapproximate_lm_cpa!(arr, M, cpa, overapprox_option) + # construct Minkowski sum for block row + function _block_row(cpa::Vector{S}, M::AbstractMatrix{N}, + row_range::UnitRange{Int}) where {N, S<:LazySet{N}} + arr_inner = Vector{LinearMap{N, <:S}}(undef, length(cpa)) + col_start_ind, col_end_ind = 1, 0 + @inbounds for (j, bj) in enumerate(cpa) + col_end_ind += dim(bj) + arr_inner[j] = + LinearMap(M[row_range, col_start_ind:col_end_ind], bj) + col_start_ind = col_end_ind + 1 + end + return MinkowskiSumArray(arr_inner) + end + + row_start_ind, row_end_ind = 1, 0 + @inbounds for (i, bi) in enumerate(cpa) + row_end_ind += dim(bi) + ms = _block_row(cpa, M, row_start_ind:row_end_ind) + arr[i] = overapproximate(ms, overapprox_option) + row_start_ind = row_end_ind + 1 + end + + return CartesianProductArray(arr) +end + +""" + overapproximate(cap::Intersection{N, + <:CartesianProductArray{N}, + <:AbstractPolyhedron{N}}, + ::Type{CartesianProductArray}, oa) where {N} + +Return the intersection of the cartesian product of a finite number of convex sets and a polyhedron. + +### Input + + - `cap` -- Lazy intersection of cartesian product array and polyhedron + - `CartesianProductArray` -- type for dispatch + - `oa` -- template directions + +### Output + +The intersection between `cpa` and `P` with overapproximation in each constrained block. + +### Algorithm + +The intersection is only needed to be taken in the elements of the cartesian product array (subsets of +variables, or "blocks") which are constrained in `P`. +Hence we first search for constrained blocks and then take the intersection of +a lower-dimensional Cartesian product of these blocks with the projection of `Y` +onto the variables of these blocks. (This projection is syntactic and exact.) +The result is a `CartesianProductArray` with the same block structure as in `X`. +However, when we decompose back our set we overapproximate during projection operation, +therefore we need to specify overapproximation strategy (it is Hyperrectangle by default) +""" +function overapproximate(cap::Intersection{N, + <:CartesianProductArray{N}, + <:AbstractPolyhedron{N}}, + ::Type{CartesianProductArray}, oa) where {N} + + cpa, P = cap.X, cap.Y + + cpa_low_dim, vars, block_structure, blocks = get_constrained_lowdimset(cpa, P) + + low_intersection = intersection(cpa_low_dim, project(P, vars)) + + if isempty(low_intersection) + return EmptySet{N}() + end + + decomposed_low_set = decompose(low_intersection, block_structure, oa) + + return substitute_blocks(decomposed_low_set, cpa, blocks) +end + +# symmetric method +function overapproximate(cap::Intersection{N, + <:AbstractPolyhedron{N}, + <:CartesianProductArray{N}}, + ::Type{CartesianProductArray}, oa) where {N} + overapproximate(Intersection(cap.Y, cap.X), oa) +end diff --git a/src/CartesianProduct.jl b/src/CartesianProduct.jl index 40fb385ef8..6910e8c48b 100644 --- a/src/CartesianProduct.jl +++ b/src/CartesianProduct.jl @@ -3,7 +3,8 @@ import Base: *, ∈, isempty export CartesianProduct, CartesianProductArray, CartesianProduct!, - array + array, + same_block_structure """ CartesianProduct{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N} @@ -562,3 +563,143 @@ function same_block_structure(x::AbstractVector{S1}, y::AbstractVector{S2} end return true end + +""" + block_structure(cpa::CartesianProductArray{N}) where {N} + +Returns an array containing the dimension ranges of each block in a `CartesianProductArray`. + +### Input + +- `cpa` -- Cartesian product array + +### Output + +A vector of ranges + +### Example + +```jldoctest +julia> using LazySets: block_structure + +julia> cpa = CartesianProductArray([BallInf(zeros(n), 1.0) for n in [3, 1, 2]]); + +julia> block_structure(cpa) +3-element Array{UnitRange{Int64},1}: + 1:3 + 4:4 + 5:6 +``` +""" +function block_structure(cpa::CartesianProductArray{N}) where {N} + result = Vector{UnitRange{Int}}(undef, length(array(cpa))) + start_index = 1 + @inbounds for (i, bi) in enumerate(array(cpa)) + end_index = start_index + dim(bi) - 1 + result[i] = start_index:end_index + start_index = end_index + 1 + end + return result +end + + +""" + block_to_dimension_indices(cpa::CartesianProductArray{N}, vars::Vector{Int}) where {N} + +Returns a vector mapping block index `i` +to tuple `(f, l)` such that either `f = l = -1` or `f` is the first dimension index and `l` is the last dimension +index of the `i`-th block, depending on whether one of the block's dimension indices is specified in `vars`. + +### Input + +- `cpa` -- Cartesian product array +- `vars` -- list containing the variables of interest, sorted in ascending order + +### Output + +A vector of tuples, where values in tuple relate to range of dimensions in the i-th block. + +### Example + +```jldoctest +julia> using LazySets: block_to_dimension_indices + +julia> cpa = CartesianProductArray([BallInf(zeros(n), 1.0) for n in [1, 3, 2, 3]]); + +julia> block_to_dimension_indices(cpa, [2, 4, 8]) +(Tuple{Int64,Int64}[(-1, -1), (2, 4), (-1, -1), (7, 9)], 2) +``` +This vector represents the mapping "second block from dimension 2 to dimension 4, +fourth block from dimension 7 to dimension 9." +These blocks contain the dimensions specified in `[2, 4, 8]`. +""" +function block_to_dimension_indices(cpa::CartesianProductArray{N}, vars::Vector{Int}) where {N} + result = fill((-1, -1), length(array(cpa))) + non_empty_length = 0 + start_index, end_index = 1, 0 + v_i = 1 + for i in 1:length(cpa.array) + end_index += dim(cpa.array[i]) + if v_i <= length(vars) && vars[v_i] <= end_index + result[i] = (start_index, end_index) + non_empty_length += 1 + while v_i <= length(vars) && vars[v_i] <= end_index + v_i += 1 + end + end + if v_i > length(vars) + break + end + start_index = end_index + 1 + end + return result, non_empty_length +end + +#same but for all variables +function block_to_dimension_indices(cpa::CartesianProductArray{N}) where {N} + result = Vector{Tuple{Int, Int}}(undef, length(cpa.array)) + + start_index, end_index = 1, 0 + for i in 1:length(cpa.array) + end_index += dim(cpa.array[i]) + result[i] = (start_index, end_index) + start_index = end_index + 1 + end + return result, length(cpa.array) +end + +""" + substitute_blocks(low_dim_cpa::CartesianProductArray{N}, + orig_cpa::CartesianProductArray{N}, + blocks::Vector{Tuple{Int,Int}}) where {N} + +Return merged Cartesian Product Array between original CPA and some low-dimensional CPA, +which represents updated subset of variables in specified blocks. + +### Input + +- `low_dim_cpa` -- low-dimensional cartesian product array +- `orig_cpa` -- original high-dimensional Cartesian product array +- `blocks` -- index of the first variable in each block of `orig_cpa` + +### Output + +Merged cartesian product array +""" +function substitute_blocks(low_dim_cpa::CartesianProductArray{N}, + orig_cpa::CartesianProductArray{N}, + blocks::Vector{Tuple{Int,Int}}) where {N} + + array = Vector{LazySet{N}}(undef, length(orig_cpa.array)) + index = 1 + for bi in 1:length(orig_cpa.array) + start_ind, end_index = blocks[bi] + if start_ind == -1 + array[bi] = orig_cpa.array[bi] + else + array[bi] = low_dim_cpa.array[index] + index += 1 + end + end + return CartesianProductArray(array) +end diff --git a/src/LazySet.jl b/src/LazySet.jl index 1438fb3b68..4e9d46c59c 100644 --- a/src/LazySet.jl +++ b/src/LazySet.jl @@ -1,4 +1,4 @@ -import Base: ==, copy +import Base: ==, ≈, copy import Random.rand export LazySet, @@ -293,12 +293,10 @@ function an_element(S::LazySet{N}) where {N<:Real} return σ(sparsevec([1], [one(N)], dim(S)), S) end - """ ==(X::LazySet, Y::LazySet) -Return whether two LazySets of the same type are exactly equal by recursively -comparing their fields until a mismatch is found. +Return whether two LazySets of the same type are exactly equal. ### Input @@ -312,10 +310,16 @@ comparing their fields until a mismatch is found. ### Notes The check is purely syntactic and the sets need to have the same base type. -I.e. `X::VPolytope == Y::HPolytope` returns `false` even if `X` and `Y` represent the -same polytope. However `X::HPolytope{Int64} == Y::HPolytope{Float64}` is a valid comparison. +For instance, `X::VPolytope == Y::HPolytope` returns `false` even if `X` and `Y` +represent the same polytope. +However `X::HPolytope{Int64} == Y::HPolytope{Float64}` is a valid comparison. + +### Algorithm + +We recursively compare the fields of `X` and `Y` until a mismatch is found. ### Examples + ```jldoctest julia> HalfSpace([1], 1) == HalfSpace([1], 1) true @@ -341,6 +345,61 @@ function ==(X::LazySet, Y::LazySet) return true end +""" + ≈(X::LazySet, Y::LazySet) + +Return whether two LazySets of the same type are approximately equal. + +### Input + +- `X` -- any `LazySet` +- `Y` -- another `LazySet` of the same type as `X` + +### Output + +- `true` iff `X` is equal to `Y`. + +### Notes + +The check is purely syntactic and the sets need to have the same base type. +For instance, `X::VPolytope ≈ Y::HPolytope` returns `false` even if `X` and `Y` +represent the same polytope. +However `X::HPolytope{Int64} ≈ Y::HPolytope{Float64}` is a valid comparison. + +### Algorithm + +We recursively compare the fields of `X` and `Y` until a mismatch is found. + +### Examples + +```jldoctest +julia> HalfSpace([1], 1) ≈ HalfSpace([1], 1) +true + +julia> HalfSpace([1], 1) ≈ HalfSpace([1.00000001], 0.99999999) +true + +julia> HalfSpace([1], 1) ≈ HalfSpace([1.0], 1.0) +true + +julia> Ball1([0.], 1.) ≈ Ball2([0.], 1.) +false +``` +""" +function ≈(X::LazySet, Y::LazySet) + # if the common supertype of X and Y is abstract, they cannot be compared + if isabstracttype(promote_type(typeof(X), typeof(Y))) + return false + end + + for f in fieldnames(typeof(X)) + if !(getfield(X, f) ≈ getfield(Y, f)) + return false + end + end + + return true +end # hook into random API function rand(rng::AbstractRNG, ::SamplerType{T}) where T<:LazySet diff --git a/src/LazySets.jl b/src/LazySets.jl index 6f83324ff2..6d01e160d7 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -11,6 +11,7 @@ using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle import InteractiveUtils: subtypes export Approximations +export × # =================== # Auxiliary functions @@ -88,8 +89,14 @@ include("Rectification.jl") # ============================= # Conversions between set types # ============================= +include("intersection_helper.jl") include("convert.jl") +# ===================== +# Approximations module +# ===================== +include("Approximations/Approximations.jl") + # =========================== # Concrete operations on sets # =========================== @@ -103,11 +110,6 @@ include("difference.jl") # ======= include("aliases.jl") -# ===================== -# Approximations module -# ===================== -include("Approximations/Approximations.jl") - # ========================== # Parallel algorithms module # ========================== diff --git a/src/LinearMap.jl b/src/LinearMap.jl index 82bcb4e0f7..32d0eba809 100644 --- a/src/LinearMap.jl +++ b/src/LinearMap.jl @@ -316,13 +316,14 @@ function isempty(lm::LinearMap)::Bool end """ - vertices_list(lm::LinearMap{N})::Vector{Vector{N}} where {N<:Real} + vertices_list(lm::LinearMap{N}; prune::Bool=true)::Vector{Vector{N}} where {N<:Real} Return the list of vertices of a (polyhedral) linear map. ### Input - `lm` -- linear map +- `prune` -- (optional, default: `true`) if true removes redundant vertices ### Output @@ -333,7 +334,7 @@ A list of vertices. We assume that the underlying set `X` is polyhedral. Then the result is just the linear map applied to the vertices of `X`. """ -function vertices_list(lm::LinearMap{N})::Vector{Vector{N}} where {N<:Real} +function vertices_list(lm::LinearMap{N}; prune::Bool=true)::Vector{Vector{N}} where {N<:Real} # for a zero map, the result is just the list containing the origin if iszero(lm.M) return [zeros(N, dim(lm))] @@ -349,7 +350,7 @@ function vertices_list(lm::LinearMap{N})::Vector{Vector{N}} where {N<:Real} push!(vlist, lm.M * v) end - return vlist + return prune ? convex_hull(vlist) : vlist end """ diff --git a/src/concrete_convex_hull.jl b/src/concrete_convex_hull.jl index f022fe08f5..3ac54df7d0 100644 --- a/src/concrete_convex_hull.jl +++ b/src/concrete_convex_hull.jl @@ -28,8 +28,8 @@ The convex hull as a list of vectors with the coordinates of the points. ### Algorithm -A pre-processing step treats the cases with `0`, `1` and `2` points in any dimension. -For more than `3` points, the algorithm used depends on the dimension. +A pre-processing step treats the cases with `0`, `1`, `2` and `3` points in any dimension. +For `4` or more points, the algorithm used depends on the dimension. For the one-dimensional case we return the minimum and maximum points, in that order. @@ -99,41 +99,91 @@ function convex_hull!(points::Vector{VN}; n = length(first(points)) - # two points - if m == 2 - if n == 1 - p1, p2 = points[1], points[2] - if p1 == p2 # check for redundancy - pop!(points) - elseif p1[1] > p2[1] - points[1], points[2] = p2, p1 - end - elseif n == 2 - # special case, see #876 - p1, p2 = points[1], points[2] - if p1 == p2 # check for redundancy - pop!(points) - elseif p1 <= p2 - nothing - else - points[1], points[2] = p2, p1 - end + if n == 1 # dimensional check + if m == 2 + # two points case in 1d + return _two_points_1d!(points) + else + # general case in 1d + return _convex_hull_1d!(points) end - return points + elseif n == 2 + if m == 2 + # two points case in 2d + return _two_points_2d!(points) + elseif m == 3 + # three points case in 2d + return _three_points_2d!(points) + else + # general case in 2d + return _convex_hull_2d!(points, algorithm=algorithm) + end + else + # general case in nd + return _convex_hull_nd!(points, backend=backend) + end +end + +function _two_points_1d!(points) + p1, p2 = points[1], points[2] + if p1 == p2 + # check for redundancy + pop!(points) + elseif p1[1] > p2[1] + points[1], points[2] = p2, p1 end + return points +end - # =========================== - # Dispatch for general cases - # =========================== +function _two_points_2d!(points) - # _convex_hull_x can assume that there are at least three points - if n == 1 - return _convex_hull_1d!(points) - elseif n == 2 - return _convex_hull_2d!(points, algorithm=algorithm) + # special case, see #876 + p1, p2 = points[1], points[2] + if p1 == p2 + # check for redundancy + pop!(points) + elseif p1 <= p2 + nothing else - return _convex_hull_nd!(points, backend=backend) + points[1], points[2] = p2, p1 + end + return points +end + +function _three_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real} + # Algorithm: the function takes three points and uses the formula + # from here: https://stackoverflow.com/questions/2122305/convex-hull-of-4-points/2122620#2122620 + # to decide if the points are ordered in a counter-clockwise fashion or not, the result is saved + # in the 'turn' boolean, then returns the points in ordered ccw acting according to 'turn'. For the + # cases where the points are collinear we pass the points with the minimum and maximum first + # component to the function for two points in 2d(_two_points_2d), if those are equal, we do the same + # but with the second component. + A, B, C = points[1], points[2], points[3] + turn = (A[2] - B[2]) * C[1] + (B[1] - A[1]) * C[2] + (A[1] * B[2] - B[1] * A[2]) + + if isapproxzero(turn) + # ABC are collinear + if isapprox(A[1], B[1]) && isapprox(B[1], C[1]) && isapprox(C[1], A[1]) + # points are approximately equal in their first component + if isapprox(A[2], B[2]) && isapprox(B[2], C[2]) && isapprox(C[2], A[2]) + # all points are approximately equal + pop!(points) + pop!(points) + return points + else + points[1], points[2] = points[argmin([A[2], B[2], C[2]])], points[argmax([A[2], B[2], C[2]])] + end + else + points[1], points[2] = points[argmin([A[1], B[1], C[1]])], points[argmax([A[1], B[1], C[1]])] + end + pop!(points) + _two_points_2d!(points) + elseif turn < zero(N) + # ABC is CW + points[1], points[2], points[3] = C, B, A end + # else ABC is CCW => nothing to do + return points end function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:AbstractVector{N}} diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl index fe43ccf754..6e4037a75f 100644 --- a/src/concrete_intersection.jl +++ b/src/concrete_intersection.jl @@ -538,8 +538,12 @@ function intersection(P1::AbstractPolyhedron{N}, # remove the redundancies removehredundancy!(Qph) - # convert back to HPOLY - return convert(HPOLY, Qph) + if isempty(Qph) + return EmptySet{N}() + else + # convert back to HPOLY + return convert(HPOLY, Qph) + end end end @@ -674,7 +678,7 @@ Return the intersection of a lazy linear map and a convex set. - `L` -- linear map - `S` -- convex set - + ### Output The polytope obtained by the intersection of `l.M * L.X` and `S`. @@ -778,3 +782,36 @@ function intersection(rm::ResetMap{N, <:AbstractPolytope}, P::AbstractPolyhedron{N}) where {N<:Real} return intersection(P, rm) end + +function intersection(U::Universe{N}, X::CartesianProductArray{N}) where {N<:Real} + return X +end + +# symmetric method +function intersection(X::CartesianProductArray{N}, U::Universe{N}) where {N<:Real} + return intersection(U, X) +end + +""" + intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) + +Return the intersection between cartesian products of a finite number of convex sets. + +### Input + + - `X` -- cartesian product of a finite number of convex sets + - `Y` -- cartesian product of a finite number of convex sets + +### Output + +The decomposed set which represents concrete intersection between `X` and `Y` + +### Algorithm + +This algorithm intersect corresponding blocks between sets. +""" +function intersection(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N<:Real} + @assert same_block_structure(array(X), array(Y)) "block structure has to be the same" + + return CartesianProductArray([intersection(array(X)[i], array(Y)[i]) for i in eachindex(X.array)]) +end diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 7982262cec..5f6b932f47 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -542,9 +542,9 @@ julia> subtypes(AbstractPolytope, true) function subtypes(interface, concrete::Bool) subtypes_to_test = subtypes(interface) - + # do not seek the concrete subtypes further - if !concrete + if !concrete return sort(subtypes_to_test, by=string) end diff --git a/src/intersection_helper.jl b/src/intersection_helper.jl new file mode 100644 index 0000000000..ca83261de9 --- /dev/null +++ b/src/intersection_helper.jl @@ -0,0 +1,49 @@ +""" + get_constrained_lowdimset(cpa::CartesianProductArray{N, S}, P::AbstractPolyhedron{N}) where {N<:Real, S<:LazySet{N}} + +Preprocess step for intersection between cartesian product array and H-polyhedron. +Returns low-dimensional HPolytope in constrained dimensions of original cpa, +constrained variables and variables in corresponding blocks, original block structure +of low-dimensional set and list of constrained blocks. + +### Input + +- `cpa` -- Cartesian product array of convex sets +- `P` -- polyhedron + +### Output + +A tuple of low-dimensional set, list of constrained dimensions, original block +structure of low-dimensional set and corresponding blocks indices. +""" +function get_constrained_lowdimset(cpa::CartesianProductArray{N, S}, P::AbstractPolyhedron{N}) where {N<:Real, S<:LazySet{N}} + + if isbounded(P) + blocks, non_empty_length = block_to_dimension_indices(cpa) + else + constrained_vars = constrained_dimensions(P) + blocks, non_empty_length = block_to_dimension_indices(cpa, constrained_vars) + end + + array = Vector{S}() + sizehint!(array, non_empty_length) + vars = Vector{Int}() + block_structure = Vector{UnitRange{Int}}() + sizehint!(block_structure, non_empty_length) + + last_var = 1 + for i in 1:length(blocks) + start_index, end_index = blocks[i] + block_end = last_var + end_index - start_index + if start_index != -1 + push!(array, cpa.array[i]) + append!(vars, start_index : end_index) + push!(block_structure, last_var : block_end) + last_var = block_end + 1 + end + end + + cpa_low_dim = HPolytope(constraints_list(CartesianProductArray(array))); + + return cpa_low_dim, vars, block_structure, blocks +end diff --git a/src/is_intersection_empty.jl b/src/is_intersection_empty.jl index 5ea84ddd62..c657afa30d 100644 --- a/src/is_intersection_empty.jl +++ b/src/is_intersection_empty.jl @@ -1,3 +1,4 @@ +using LazySets.Approximations: project export is_intersection_empty, isdisjoint @@ -1365,3 +1366,57 @@ function is_intersection_empty(X::LazySet{N}, {N<:Real} return is_intersection_empty(C, X, witness) end + +""" + is_intersection_empty(X::CartesianProductArray{N, S}, + P::HPolyhedron{N}) where {N<:Real, S<:LazySet{N}} + +Check whether a Cartesian product array intersects with a H-polyhedron. + +### Input + +- `X` -- Cartesian product array of convex sets +- `P` -- H-polyhedron + +### Output + +`true` iff ``X ∩ Y = ∅ ``. +""" +function is_intersection_empty(X::CartesianProductArray{N}, + P::HPolyhedron{N}) where {N<:Real} + cpa_low_dim, vars, block_structure = get_constrained_lowdimset(X, P) + + return isdisjoint(cpa_low_dim, project(P, vars)) +end + +# symmetric method +function is_intersection_empty(P::HPolyhedron{N}, + X::CartesianProductArray{N, S}) where {N<:Real, S<:LazySet{N}} + is_intersection_empty(X, P) +end + +""" + is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N} + +Check whether two Cartesian products of a finite number of convex sets do not intersect. + +### Input + +- `X` -- Cartesian product array of convex sets +- `Y` -- Cartesian product array of convex sets + +### Output + +`true` iff ``X ∩ Y = ∅ ``. +""" +function is_intersection_empty(X::CartesianProductArray{N}, Y::CartesianProductArray{N}) where {N} + @assert same_block_structure(array(X), array(Y)) "block structure has to be the same" + + for i in 1:length(X.array) + if isdisjoint(X.array[i], Y.array[i]) + return true + end + end + + return false +end diff --git a/src/plot_recipes.jl b/src/plot_recipes.jl index 6caa210c27..14e460c9b2 100644 --- a/src/plot_recipes.jl +++ b/src/plot_recipes.jl @@ -8,7 +8,7 @@ DEFAULT_COLOR = :auto DEFAULT_ALPHA = 0.5 DEFAULT_LABEL = "" DEFAULT_GRID = true -DEFAULT_ASPECT_RATIO = 1.0 +DEFAULT_ASPECT_RATIO = :none PLOT_PRECISION = 1e-3 PLOT_POLAR_DIRECTIONS = 40 @@ -68,7 +68,9 @@ julia> plot(Bs, 1e-2) # faster but less accurate than the previous call if fast label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end seriesalpha --> DEFAULT_ALPHA seriescolor --> DEFAULT_COLOR seriestype --> :shape @@ -153,7 +155,9 @@ julia> plot(B, 1e-2) # faster but less accurate than the previous call else label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end seriesalpha --> DEFAULT_ALPHA seriescolor --> DEFAULT_COLOR @@ -192,7 +196,9 @@ julia> plot(Singleton([0.5, 1.0])) ) where {N<:Real} label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end seriesalpha --> DEFAULT_ALPHA seriescolor --> DEFAULT_COLOR seriestype := :scatter @@ -246,7 +252,9 @@ julia> plot(L, marker=0) ε::N=zero(N)) where {N<:Real} label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end seriesalpha --> DEFAULT_ALPHA linecolor --> DEFAULT_COLOR markercolor --> DEFAULT_COLOR @@ -269,7 +277,9 @@ Plot an empty set. @recipe function plot_emptyset(∅::EmptySet{N}, ε::N=zero(N)) where {N<:Real} label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end plot_recipe(∅) end @@ -325,7 +335,9 @@ julia> plot(X, -1., 100) # equivalent to the above line ) where {N<:Real} label --> DEFAULT_LABEL grid --> DEFAULT_GRID - aspect_ratio --> DEFAULT_ASPECT_RATIO + if DEFAULT_ASPECT_RATIO != :none + aspect_ratio --> DEFAULT_ASPECT_RATIO + end seriesalpha --> DEFAULT_ALPHA seriescolor --> DEFAULT_COLOR seriestype := :shape diff --git a/test/runtests.jl b/test/runtests.jl index f6147d986e..b5b6b68b83 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,12 @@ using LazySets, LazySets.Approximations, Test, LinearAlgebra, SparseArrays + import IntervalArithmetic, Expokit, Optim +const IA = IntervalArithmetic using IntervalArithmetic: IntervalBox +import TaylorModels +using TaylorModels: set_variables, TaylorModelN + # conversion between numeric types include("to_N.jl") @@ -88,6 +93,7 @@ if test_suite_basic @time @testset "LazySets.CartesianProduct" begin include("unit_CartesianProduct.jl") end @time @testset "LazySets.ResetMap" begin include("unit_ResetMap.jl") end @time @testset "LazySets.SymmetricIntervalHull" begin include("unit_SymmetricIntervalHull.jl") end + @time @testset "LazySets.concrete_convex_hull" begin include("unit_convex_hull.jl") end # ====================== # Testing set interfaces diff --git a/test/unit_CartesianProduct.jl b/test/unit_CartesianProduct.jl index cf6812f3b2..1133c199d7 100644 --- a/test/unit_CartesianProduct.jl +++ b/test/unit_CartesianProduct.jl @@ -1,3 +1,5 @@ +using LazySets.Approximations: overapproximate + for N in [Float64, Float32, Rational{Int}] # Cartesian Product of a centered 1D BallInf and a centered 2D BallInf # Here a 3D BallInf @@ -261,7 +263,7 @@ for N in [Float64, Float32, Rational{Int}] Hcp = convert(CartesianProductArray{N, Interval{N}}, H) @test Hcp isa CartesianProductArray && all([array(Hcp)[i] == Interval(N(-1), N(2*i-1)) for i in 1:3]) - + # ================ # common functions # ================ @@ -273,3 +275,21 @@ for N in [Float64, Float32, Rational{Int}] EmptySet @test b × e == e × b == cpa × e == e × cpa == e × e == e end + +for N in [Float64, Float32] + # is_intersection_empty + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + h1 = Hyperrectangle(low=N[3, 4], high=N[5, 7]) + h2 = Hyperrectangle(low=N[5, 5], high=N[6, 8]) + cpa1 = CartesianProductArray([i1, i2, h1]) + cpa2 = CartesianProductArray([i1, i2, h2]) + G = HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(1))]) + G_empty = HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(-1))]) + cpa1_box = overapproximate(cpa1) + cpa2_box = overapproximate(cpa2) + + @test is_intersection_empty(cpa1, cpa2) == is_intersection_empty(cpa1_box, cpa2_box) == false + @test is_intersection_empty(cpa1, G_empty) == is_intersection_empty(cpa1_box, G_empty) == true + @test is_intersection_empty(cpa1, G) == is_intersection_empty(Approximations.overapproximate(cpa1), G) == false +end diff --git a/test/unit_LinearMap.jl b/test/unit_LinearMap.jl index 9b49dcb501..026785998a 100644 --- a/test/unit_LinearMap.jl +++ b/test/unit_LinearMap.jl @@ -84,6 +84,11 @@ for N in [Float64, Rational{Int}, Float32] vlist = vertices_list(LinearMap(M, b)) @test vlist == [zeros(N, 2)] + # test that redundant vertices are removed by default (see #1355) + X = Interval(N(0), N(1)) × ZeroSet{N}(1) + M = N[0 1; 0 2] + @test vertices_list(M * X) == [N[0, 0]] + if test_suite_polyhedra # constraints_list b = BallInf(N[0, 0], N(1)) diff --git a/test/unit_Polytope.jl b/test/unit_Polytope.jl index efac9572a5..b023133202 100644 --- a/test/unit_Polytope.jl +++ b/test/unit_Polytope.jl @@ -222,26 +222,24 @@ if test_suite_polyhedra LinearConstraint(N[-1], N(-1))]) @test_throws ErrorException σ(N[1], p_infeasible) - # intersection - A = [N(0) N(1); N(1) N(0); N(2) N(2)] + # empty intersection + A = [N(0) N(-1); N(-1) N(0); N(1) N(1)] b = N[0, 0, 1] p1 = HPolytope(A, b) - A = [N(0) N(-1); N(-1) N(0); N(1) N(1)] - b = N[-0.25, -0.25, -0] + A = [N(0) N(1); N(1) N(0); N(-1) N(-1)] + b = N[2, 2, -2] p2 = HPolytope(A, b) - cap = intersection(p1, p2) - # currently broken, see #565 + @test intersection(p1, p2) isa EmptySet{N} + @test intersection(p1, p2; backend=Polyhedra.default_library(2, N)) isa EmptySet{N} + # @test intersection(p1, p2; backend=CDDLib.Library()) isa EmptySet{N} + # commented because we do not load CDDLib at the moment # intersection with half-space - hs = HalfSpace(N[2, 2], N(-1)) + hs = HalfSpace(N[2, -2], N(-1)) c1 = constraints_list(intersection(p1, hs)) c2 = constraints_list(intersection(hs, p1)) @test length(c1) == 3 && ispermutation(c1, c2) - # convex hull - ch = convex_hull(p1, p2) - # currently broken, see #566 - # Cartesian product A = [N(1) N(-1)]' b = N[1, 0] diff --git a/test/unit_convex_hull.jl b/test/unit_convex_hull.jl index e9be3b3530..3916e30d62 100644 --- a/test/unit_convex_hull.jl +++ b/test/unit_convex_hull.jl @@ -33,16 +33,33 @@ for N in [Float64, Rational{Int}] convex_hull!(points_2D) # check in-place version @test points_2D == [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]] + # three vertex case in 2 dimensions + function iscounterclockwise(result, correct_expr) + # this function checks if the result equals any of the correct answer cyclic permutation + result == correct_expr || result == circshift(correct_expr, 1) || result == circshift(correct_expr, 2) + end + ccw_points = [N[1, 1], N[-1, 1], N[-1, 0]] + ccw_p = convex_hull!(ccw_points) + ccw_expr = [N[1, 1], N[-1, 1], N[-1, 0]] + @test iscounterclockwise(ccw_p, ccw_expr) # points sorted in a counter-clockwise fashion + cw_points = [N[-1, 1], N[1, 1], N[-1, 0]] + cw_p = convex_hull!(cw_points) + cw_expr = [N[1, 1], N[-1, 1], N[-1, 0]] + @test iscounterclockwise(cw_p, cw_expr) # points sorted in clockwise fashion + @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3]]), [N[1, 1], N[3, 3]]) # points aligned + @test ispermutation(convex_hull!([N[0, 1], N[0, 2], N[0, 3]]), [N[0, 1], N[0, 3]]) # three points on a vertical line + @test convex_hull!([N[0, 1], N[0, 1], N[0, 1]]) == [N[0, 1]] # three equal points + # higher dimension if test_suite_polyhedra && N != Float32 # no backend supporting Float32 points_3D = [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(0), N(0), N(0)], [N(1), N(2), N(3)]] - @test convex_hull(points_3D) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(convex_hull(points_3D), [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) convex_hull!(points_3D) # check in-place version - @test points_3D == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(points_3D, [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) end @@ -53,7 +70,7 @@ for N in [Float64, Rational{Int}] V1 = VPolygon([[N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) V2 = VPolygon([[N(-1), N(-1)], [N(1/2), N(1/2)]]) ch = convex_hull(V1, V2) - @test vertices_list(ch) == ispermutation([[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) + @test ispermutation(vertices_list(ch), [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]) if test_suite_polyhedra && N != Float32 # no backend supporting Float32 V1 = VPolytope([[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)]]) @@ -61,7 +78,7 @@ for N in [Float64, Rational{Int}] [N(1), N(2), N(3)]]) ch = convex_hull(V1, V2) @test ch isa VPolytope - @test vertices_list(ch) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)], + @test ispermutation(vertices_list(ch), [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]) end diff --git a/test/unit_overapproximate.jl b/test/unit_overapproximate.jl index e318ea8273..fd6c00647b 100644 --- a/test/unit_overapproximate.jl +++ b/test/unit_overapproximate.jl @@ -1,3 +1,5 @@ +using LazySets.Approximations: project + for N in [Float64, Rational{Int}, Float32] # overapproximating a set of type T1 with an unsupported type T2 is the # identity if T1 = T2 @@ -47,12 +49,48 @@ for N in [Float64, Rational{Int}, Float32] h2 = Hyperrectangle(N[2.5, 4.5], N[1/2, 1/2]) H = overapproximate(h1 × h2, Hyperrectangle) # defaults to convert method @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5] - + # overapproximation of the lazy linear map of a hyperrectangular set H = Hyperrectangle(N[0, 0], N[1/2, 1]) M = Diagonal(N[2, 2]) OA = overapproximate(M*H, Hyperrectangle) @test OA isa Hyperrectangle && OA.center == N[0, 0] && OA.radius == N[1, 2] + + #overapproximation of Minkowski sum of linear maps for each block in the row block + i1 = Interval(N[0, 1]) + h = Hyperrectangle(low=N[3, 4], high=N[5, 7]) + M = N[1 2 3; 4 5 6; 7 8 9] + cpa = CartesianProductArray([i1, h]) + lm = M * cpa + + oa = overapproximate(lm, Hyperrectangle) + oa_box = overapproximate(lm, Approximations.BoxDirections) + d_oa_d_hp = overapproximate(lm, CartesianProductArray{N, Hyperrectangle{N}}) + d_oa_d_box = overapproximate(lm, CartesianProductArray, Approximations.BoxDirections) + oa_d_hp = overapproximate(d_oa_d_hp) + oa_d_box = overapproximate(d_oa_d_box, Approximations.BoxDirections) + + @test oa == oa_d_hp + @test oa_box == oa_d_box + + for (oax, set_type) in [(d_oa_d_hp, Hyperrectangle), (d_oa_d_box, HPolytope)] + @test oax isa CartesianProductArray + arr = oax.array + @test length(arr) == 2 && dim(arr[1]) == 1 && dim(arr[2]) == 2 + @test all(X -> X isa set_type, arr) + end + + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + i3 = Interval(N[1, 4]) + cpa = CartesianProductArray([i1, i2, i3]) + M = N[1 2 0; 0 1 0; 0 1 1] + lm = M * cpa + d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) + oa = overapproximate(lm) + @test overapproximate(d_oa) == oa + @test typeof(d_oa) == CartesianProductArray{N, Interval{N}} + end # tests that do not work with Rational{Int} @@ -129,4 +167,76 @@ for N in [Float64, Float32] Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope @test Y_polygon ⊆ Y_zonotope @test !(Y_zonotope ⊆ Y_polygon) + + # decomposed intersection between cartesian product array and abstract polyhedron + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + h1 = Hyperrectangle(low=N[3, 4], high=N[5, 7]) + h2 = Hyperrectangle(low=N[5, 5], high=N[6, 8]) + cpa1 = CartesianProductArray([i1, i2, h1]) + o_cpa1 = overapproximate(cpa1) + cpa2 = CartesianProductArray([i1, i2, h2]) + o_cpa2 = overapproximate(cpa2) + G = HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(1))]) + G_3 = HPolyhedron([HalfSpace(N[0, 0, 1, 0], N(1))]) + G_3_neg = HPolyhedron([HalfSpace(N[0, 0, -1, 0], N(0))]) + + d_int_g = Intersection(cpa1, G) + d_int_g_3 = Intersection(cpa1, G_3) + d_int_g_3_neg = Intersection(cpa1, G_3_neg) + d_int_cpa = intersection(cpa1, cpa2) + l_int_cpa = intersection(o_cpa1, o_cpa2) + o_d_int_g = overapproximate(d_int_g, CartesianProductArray, Hyperrectangle) + o_d_int_g_3 = overapproximate(d_int_g_3, CartesianProductArray, Hyperrectangle) + o_d_int_g_3_neg = overapproximate(d_int_g_3_neg, CartesianProductArray, Hyperrectangle) + + @test overapproximate(o_d_int_g) == overapproximate(d_int_g) + @test overapproximate(o_d_int_g_3) == overapproximate(d_int_g_3) == EmptySet{N}() + @test overapproximate(o_d_int_g_3_neg) == overapproximate(d_int_g_3_neg) + @test overapproximate(d_int_cpa, Hyperrectangle) == l_int_cpa + @test all([X isa CartesianProductArray for X in [d_int_cpa, o_d_int_g, o_d_int_g_3_neg]]) +end + +for N in [Float64] # due to sparse vectors: a = sparse(Float32[1 -1; 1 1];); a \ Float32[4, 10] + # decomposed linear map approximation + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + M = N[1 2; 0 1] + cpa = CartesianProductArray([i1, i2]) + lm = M * cpa + d_oa = overapproximate(lm, CartesianProductArray{N, Interval{N}}) + oa = overapproximate(lm, OctDirections) + @test oa ⊆ d_oa + + # ======================================= + # Zonotope overapprox. of a Taylor model + # ======================================= + x₁, x₂, x₃ = set_variables(N, ["x₁", "x₂", "x₃"], order=5) + Dx₁ = IA.Interval(N(1.0), N(3.0)) + Dx₂ = IA.Interval(N(-1.0), N(1.0)) + Dx₃ = IA.Interval(N(-1.0), N(0.0)) + D = Dx₁ × Dx₂ × Dx₃ # domain + x0 = IntervalBox(IA.mid.(D)...) + I = IA.Interval(N(0.0), N(0.0)) # interval remainder + p₁ = 1 + x₁ - x₂ + p₂ = x₃ - x₁ + vTM = [TaylorModels.TaylorModelN(pi, I, x0, D) for pi in [p₁, p₂]] + Z1 = overapproximate(vTM, Zonotope) + @test center(Z1) == N[3, -2.5] + @test Matrix(generators(Z1)) == N[1 -1 0; -1 0 0.5] + + # decomposed intersection between cartesian product array and abstract polyhedron + i1 = Interval(N[0, 1]) + i2 = Interval(N[2, 3]) + h1 = Hyperrectangle(low=N[3, 4], high=N[5, 7]) + cpa = CartesianProductArray([i1, i2, h1]) + G_comb = HPolyhedron([HalfSpace(N[1, 1, 0, 0], N(2.5))]) + int_g_comb = Intersection(cpa, G_comb) + o_d_int_g_comb = overapproximate(int_g_comb, CartesianProductArray, Hyperrectangle) + o_bd_int_g_comb = overapproximate(int_g_comb, BoxDirections) + @test overapproximate(o_d_int_g_comb) ≈ overapproximate(o_bd_int_g_comb) + projection = project(o_bd_int_g_comb, [1, 2], BoxDirections) + @test vertices_list(projection) ≈ [[0.5, 2.5], [0., 2.5], [0., 2], [0.5, 2.]] + projection = project(overapproximate(int_g_comb, OctDirections), [1, 2], OctDirections) + @test vertices_list(projection) ≈ [[0., 2.5], [0., 2.], [0.5, 2.]] end