diff --git a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md index cb24e34ae9..9f2e72e7ea 100644 --- a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md +++ b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md @@ -1,12 +1,16 @@ ```@meta -CurrentModule = LazySets +CurrentModule = LazySets.SimpleSparsePolynomialZonotopeModule ``` # [SimpleSparsePolynomialZonotope](@id def_SimpleSparsePolynomialZonotope) ```@docs SimpleSparsePolynomialZonotope -PolynomialZonotope +``` + +## Operations + +```@docs center(::SimpleSparsePolynomialZonotope) convex_hull(::SimpleSparsePolynomialZonotope) expmat(::SimpleSparsePolynomialZonotope) @@ -20,6 +24,16 @@ quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope) where {N, MT<:Abst quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} ``` +```@meta +CurrentModule = LazySets +``` + +Alias: + +```@docs +PolynomialZonotope +``` + Inherited from [`AbstractPolynomialZonotope`](@ref): * [`dim`](@ref dim(::AbstractPolynomialZonotope)) * [`order`](@ref dim(::AbstractPolynomialZonotope)) diff --git a/src/LazySets.jl b/src/LazySets.jl index d533fe4a3b..5fb43205ce 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -138,8 +138,21 @@ include("Sets/Singleton/SingletonModule.jl") include("Sets/LineSegment/LineSegmentModule.jl") @reexport using ..LineSegmentModule: LineSegment -include("Sets/SimpleSparsePolynomialZonotope.jl") include("Sets/SparsePolynomialZonotope.jl") + +include("Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl") +@reexport using ..SimpleSparsePolynomialZonotopeModule: SimpleSparsePolynomialZonotope, + SSPZ, + quadratic_map +using ..SimpleSparsePolynomialZonotopeModule: _remove_redundant_generators_polyzono + +""" + PolynomialZonotope = SimpleSparsePolynomialZonotope + +Alias for `SimpleSparsePolynomialZonotope`. +""" +const PolynomialZonotope = SimpleSparsePolynomialZonotope + include("Sets/VPolygon.jl") include("Sets/VPolytope.jl") include("Sets/Tetrahedron.jl") diff --git a/src/Sets/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope.jl deleted file mode 100644 index fa4dbdb04a..0000000000 --- a/src/Sets/SimpleSparsePolynomialZonotope.jl +++ /dev/null @@ -1,451 +0,0 @@ -export SimpleSparsePolynomialZonotope, PolynomialZonotope, expmat, nparams, - quadratic_map, remove_redundant_generators - -""" - SimpleSparsePolynomialZonotope{N, VN<:AbstractVector{N}, - MN<:AbstractMatrix{N}, - ME<:AbstractMatrix{Int}} - <: AbstractPolynomialZonotope{N} - -Type that represents a sparse polynomial zonotope that is *simple* in the sense -that there is no distinction between independent and dependent generators. - -A simple sparse polynomial zonotope ``\\mathcal{PZ} ⊂ ℝ^n`` is -represented by the set -```math -\\mathcal{PZ} = \\left\\{x ∈ ℝ^n : x = c + ∑_{i=1}^h \\left(∏_{k=1}^p α_k^{E_{k, i}} \\right) g_i,~~ α_k ∈ [-1, 1]~~ ∀ i = 1,…,p \\right\\}, -``` -where ``c ∈ ℝ^n`` is the offset vector (or center), -``G ∈ ℝ^{n × h}`` is the generator matrix with columns ``g_i`` -(each ``g_i`` is called a *generator*), and where ``E ∈ \\mathbb{N}^{p×h}_{≥0}`` -is the exponent matrix with matrix elements ``E_{k, i}``. - -### Fields - -- `c` -- offset vector -- `G` -- generator matrix -- `E` -- exponent matrix - -### Notes - -Sparse polynomial zonotopes were introduced in [1]. The *simple* variation -was defined in [2]. - -- [1] N. Kochdumper and M. Althoff. *Sparse Polynomial Zonotopes: A Novel Set -Representation for Reachability Analysis*. Transactions on Automatic Control, -2021. -- [2] N. Kochdumper. *Challenge Problem 5: Polynomial Zonotopes in Julia.* -JuliaReach and JuliaIntervals Days 3, 2021. -""" -struct SimpleSparsePolynomialZonotope{N,VN<:AbstractVector{N}, - MN<:AbstractMatrix{N}, - ME<:AbstractMatrix{Int}} <: - AbstractPolynomialZonotope{N} - c::VN - G::MN - E::ME - - function SimpleSparsePolynomialZonotope(c::VN, G::MN, - E::ME) where {N, - VN<:AbstractVector{N}, - MN<:AbstractMatrix{N}, - ME<:AbstractMatrix{Int}} - @assert length(c) == size(G, 1) throw(DimensionMismatch("c and G " * - "should have the same number of rows")) - @assert size(G, 2) == size(E, 2) throw(DimensionMismatch("G and E " * - "should have the same number of columns")) - @assert all(>=(0), E) throw(ArgumentError("E should contain " * - "non-negative integers")) - - return new{N,VN,MN,ME}(c, G, E) - end -end - -""" - PolynomialZonotope = SimpleSparsePolynomialZonotope - -Alias for `SimpleSparsePolynomialZonotope`. - -### Notes - -Another shorthand is `SSPZ`. -""" -const PolynomialZonotope = SimpleSparsePolynomialZonotope - -# short-hand -const SSPZ = SimpleSparsePolynomialZonotope - -function isoperationtype(P::Type{<:SimpleSparsePolynomialZonotope}) - return false -end - -""" - ngens(P::SimpleSparsePolynomialZonotope) - -Return the number of generators of a simple sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The number of generators of `P`. - -### Notes - -This number corresponds to the number of monomials in the polynomial -representation of `P`. -""" -ngens(P::SSPZ) = size(P.G, 2) - -""" - nparams(P::SimpleSparsePolynomialZonotope) - -Return the number of parameters in the polynomial representation of a simple -sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The number of parameters in the polynomial representation of P. - -### Notes - -This number corresponds to the number of rows in the exponent matrix ``E`` (`p` -in the mathematical set definition). - -### Examples - -```jldoctest -julia> S = SimpleSparsePolynomialZonotope([2.0, 0], [1 2;2 2.], [1 4;1 2]) -SimpleSparsePolynomialZonotope{Float64, Vector{Float64}, Matrix{Float64}, Matrix{Int64}}([2.0, 0.0], [1.0 2.0; 2.0 2.0], [1 4; 1 2]) - -julia> nparams(S) -2 -``` -""" -nparams(P::SSPZ) = size(P.E, 1) - -function polynomial_order(P::SSPZ) - return maximum(sum, eachcol(expmat(P))) -end - -""" - center(P::SimpleSparsePolynomialZonotope) - -Return the center of a simple sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The center of `P`. -""" -center(P::SSPZ) = P.c - -""" - genmat(P::SimpleSparsePolynomialZonotope) - -Return the matrix of generators of a simple sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The matrix of generators of `P`. -""" -genmat(P::SSPZ) = P.G - -""" - expmat(P::SimpleSparsePolynomialZonotope) - -Return the matrix of exponents of a simple sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The matrix of exponents, where each column is a multidegree. - -### Notes - -In the exponent matrix, each row corresponds to a parameter (``\alpha_k`` in the -mathematical set definition) and each column corresponds to a monomial. - -### Examples - -```jldoctest -julia> S = SimpleSparsePolynomialZonotope([2.0, 0], [1 2;2 2.], [1 4;1 2]) -SimpleSparsePolynomialZonotope{Float64, Vector{Float64}, Matrix{Float64}, Matrix{Int64}}([2.0, 0.0], [1.0 2.0; 2.0 2.0], [1 4; 1 2]) - -julia> expmat(S) -2×2 Matrix{Int64}: - 1 4 - 1 2 -``` -""" -expmat(P::SSPZ) = P.E - -""" - linear_map(M::AbstractMatrix, P::SimpleSparsePolynomialZonotope) - -Apply the linear map `M` to a simple sparse polynomial zonotope. - -### Input - -- `M` -- matrix -- `P` -- simple sparse polynomial zonotope - -### Output - -The set resulting from applying the linear map `M` to `P`. -""" -function linear_map(M::AbstractMatrix, P::SSPZ) - return SimpleSparsePolynomialZonotope(M * center(P), M * genmat(P), expmat(P)) -end - -""" - quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) - where {N, MT<:AbstractMatrix{N}} - -Return the quadratic map of a simple sparse polynomial zonotope. - -### Input - -- `Q` -- vector of square matrices -- `S` -- simple sparse polynomial zonotope - -### Output - -The quadratic map of `P` represented as a simple sparse polynomial zonotope. - -### Algorithm - -This method implements Proposition 12 in [1]. -See also Proposition 3.1.30 in [2]. - -[1] N. Kochdumper, M. Althoff. *Sparse polynomial zonotopes: A novel set -representation for reachability analysis*. 2021 -[2] N. Kochdumper. *Extensions of polynomial zonotopes and their application to -verification of cyber-physical systems*. 2021. -""" -function quadratic_map(Q::Vector{MT}, - S::SimpleSparsePolynomialZonotope) where {N,MT<:AbstractMatrix{N}} - m = length(Q) - c = center(S) - h = ngens(S) - G = genmat(S) - E = expmat(S) - - cnew = similar(c, m) - Gnew = similar(G, m, h^2 + h) - QiG = similar(Q) - @inbounds for (i, Qi) in enumerate(Q) - cnew[i] = dot(c, Qi, c) - Gnew[i, 1:h] = c' * (Qi + Qi') * G - QiG[i] = Qi * G - end - - Enew = repeat(E, 1, h + 1) - @inbounds for i in 1:h - idxstart = h * i + 1 - idxend = (i + 1) * h - Enew[:, idxstart:idxend] .+= E[:, i] - for j in eachindex(QiG) - Gnew[j, idxstart:idxend] = G[:, i]' * QiG[j] - end - end - Z = SimpleSparsePolynomialZonotope(cnew, Gnew, Enew) - return remove_redundant_generators(Z) -end - -""" - quadratic_map(Q::Vector{MT}, S1::SimpleSparsePolynomialZonotope, - S2::SimpleSparsePolynomialZonotope) - where {N, MT<:AbstractMatrix{N}} - -Return the quadratic map of two simple sparse polynomial zonotopes. -The quadratic map is the set -```math - \\{x \\mid xᵢ = s₁ᵀQᵢs₂, s₁ ∈ S₁, s₂ ∈ S₂, Qᵢ ∈ Q\\}. -``` - -### Input - -- `Q` -- vector of square matrices -- `S1` -- simple sparse polynomial zonotope -- `S2` -- simple sparse polynomial zonotope - -### Output - -The quadratic map of the given simple sparse polynomial zonotopes represented as -a simple sparse polynomial zonotope. - -### Algorithm - -This method implements Proposition 3.1.30 in [1]. - -[1] N. Kochdumper. *Extensions of polynomial zonotopes and their application to -verification of cyber-physical systems*. 2021. -""" -function quadratic_map(Q::Vector{MT}, S1::SimpleSparsePolynomialZonotope, - S2::SimpleSparsePolynomialZonotope) where {N,MT<:AbstractMatrix{N}} - @assert nparams(S1) == nparams(S2) - - c1 = center(S1) - c2 = center(S2) - G1 = genmat(S1) - G2 = genmat(S2) - E1 = expmat(S1) - E2 = expmat(S2) - - c = [dot(c1, Qi, c2) for Qi in Q] - - Ghat1 = reduce(vcat, c2' * Qi' * G1 for Qi in Q) - Ghat2 = reduce(vcat, c1' * Qi * G2 for Qi in Q) - - Gbar = reduce(hcat, reduce(vcat, gj' * Qi * G2 for Qi in Q) for gj in eachcol(G1)) - Ebar = reduce(hcat, E2 .+ e1j for e1j in eachcol(E1)) - - G = hcat(Ghat1, Ghat2, Gbar) - E = hcat(E1, E2, Ebar) - - return remove_redundant_generators(SimpleSparsePolynomialZonotope(c, G, E)) -end - -""" - remove_redundant_generators(S::SimpleSparsePolynomialZonotope) - -Remove redundant generators from a simple sparse polynomial zonotope. - -### Input - -- `S` -- simple sparse polynomial zonotope - -### Output - -A new simple sparse polynomial zonotope such that redundant generators have been -removed. - -## Notes - -The result uses dense arrays irrespective of the array type of `S`. - -### Algorithm - -Let `G` be the generator matrix and `E` the exponent matrix of `S`. The -following simplifications are performed: - -- Zero columns in `G` and the corresponding columns in `E` are removed. -- For zero columns in `E`, the corresponding column in `G` is summed to the - center. -- Repeated columns in `E` are grouped together by summing the corresponding - columns in `G`. -""" -function remove_redundant_generators(S::SimpleSparsePolynomialZonotope) - c, G, E = _remove_redundant_generators_polyzono(center(S), genmat(S), expmat(S)) - - return SimpleSparsePolynomialZonotope(c, G, E) -end - -function _remove_redundant_generators_polyzono(c, G, E) - Gnew = Matrix{eltype(G)}(undef, size(G, 1), 0) - Enew = Matrix{eltype(E)}(undef, size(E, 1), 0) - cnew = copy(c) - - visited_exps = Dict{Vector{Int},Int}() - @inbounds for (gi, ei) in zip(eachcol(G), eachcol(E)) - all(isapproxzero, gi) && continue - if iszero(ei) - cnew += gi - elseif haskey(visited_exps, ei) # repeated exponent - idx = visited_exps[ei] - Gnew[:, idx] += gi - else - Gnew = hcat(Gnew, gi) - Enew = hcat(Enew, ei) - visited_exps[ei] = size(Enew, 2) - end - end - - return cnew, Gnew, Enew -end - -""" - rand(::Type{SimpleSparsePolynomialZonotope}; - [N]::Type{<:Real}=Float64, [dim]::Int=2, [nparams]::Int=2, - [maxdeg]::Int=3, [num_generators]::Int=-1, - [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) - -Create a random simple sparse polynomial zonotope. - -### Input - -- `Zonotope` -- type for dispatch -- `N` -- (optional, default: `Float64`) numeric type -- `dim` -- (optional, default: 2) dimension -- `nparams` -- (optional, default: 2) number of parameters -- `maxdeg` -- (optional, default: 3) maximum degree for each parameter -- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator -- `seed` -- (optional, default: `nothing`) seed for reseeding -- `num_generators` -- (optional, default: `-1`) number of generators of the - zonotope (see comment below) - -### Output - -A random simple sparse polynomial zonotope. - -### Algorithm - -All numbers are normally distributed with mean 0 and standard deviation 1. - -The number of generators can be controlled with the argument `num_generators`. -For a negative value we choose a random number in the range `dim:2*dim` (except -if `dim == 1`, in which case we only create a single generator). Note that the -final number of generators may be lower if redundant monomials are generated. -""" -function rand(::Type{SimpleSparsePolynomialZonotope}; - N::Type{<:Real}=Float64, - dim::Int=2, - nparams::Int=2, - maxdeg::Int=3, - rng::AbstractRNG=GLOBAL_RNG, - seed::Union{Int,Nothing}=nothing, - num_generators::Int=-1) - rng = reseed!(rng, seed) - center = randn(rng, N, dim) - if num_generators < 0 - num_generators = (dim == 1) ? 1 : rand(rng, dim:(2 * dim)) - end - generators = randn(rng, N, dim, num_generators) - expmat = rand(rng, 0:maxdeg, nparams, num_generators) - SSPZ = SimpleSparsePolynomialZonotope(center, generators, expmat) - return remove_redundant_generators(SSPZ) -end - -""" - convex_hull(P::SimpleSparsePolynomialZonotope) - -Compute the convex hull of a simple sparse polynomial zonotope. - -### Input - -- `P` -- simple sparse polynomial zonotope - -### Output - -The tightest convex simple sparse polynomial zonotope containing `P`. -""" -function convex_hull(P::SimpleSparsePolynomialZonotope) - return linear_combination(P, P) -end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotope.jl b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotope.jl new file mode 100644 index 0000000000..0a23af8736 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotope.jl @@ -0,0 +1,59 @@ +""" + SimpleSparsePolynomialZonotope{N, VN<:AbstractVector{N}, + MN<:AbstractMatrix{N}, + ME<:AbstractMatrix{Int}} + <: AbstractPolynomialZonotope{N} + +Type that represents a sparse polynomial zonotope that is *simple* in the sense +that there is no distinction between independent and dependent generators. + +A simple sparse polynomial zonotope ``\\mathcal{PZ} ⊂ ℝ^n`` is +represented by the set +```math +\\mathcal{PZ} = \\left\\{x ∈ ℝ^n : x = c + ∑_{i=1}^h \\left(∏_{k=1}^p α_k^{E_{k, i}} \\right) g_i,~~ α_k ∈ [-1, 1]~~ ∀ i = 1,…,p \\right\\}, +``` +where ``c ∈ ℝ^n`` is the offset vector (or center), +``G ∈ ℝ^{n × h}`` is the generator matrix with columns ``g_i`` +(each ``g_i`` is called a *generator*), and where ``E ∈ \\mathbb{N}^{p×h}_{≥0}`` +is the exponent matrix with matrix elements ``E_{k, i}``. + +### Fields + +- `c` -- offset vector +- `G` -- generator matrix +- `E` -- exponent matrix + +### Notes + +Sparse polynomial zonotopes were introduced in [1]. The *simple* variation +was defined in [2]. + +- [1] N. Kochdumper and M. Althoff. *Sparse Polynomial Zonotopes: A Novel Set +Representation for Reachability Analysis*. Transactions on Automatic Control, +2021. +- [2] N. Kochdumper. *Challenge Problem 5: Polynomial Zonotopes in Julia.* +JuliaReach and JuliaIntervals Days 3, 2021. +""" +struct SimpleSparsePolynomialZonotope{N,VN<:AbstractVector{N}, + MN<:AbstractMatrix{N}, + ME<:AbstractMatrix{Int}} <: + AbstractPolynomialZonotope{N} + c::VN + G::MN + E::ME + + function SimpleSparsePolynomialZonotope(c::VN, G::MN, + E::ME) where {N, + VN<:AbstractVector{N}, + MN<:AbstractMatrix{N}, + ME<:AbstractMatrix{Int}} + @assert length(c) == size(G, 1) throw(DimensionMismatch("c and G " * + "should have the same number of rows")) + @assert size(G, 2) == size(E, 2) throw(DimensionMismatch("G and E " * + "should have the same number of columns")) + @assert all(>=(0), E) throw(ArgumentError("E should contain " * + "non-negative integers")) + + return new{N,VN,MN,ME}(c, G, E) + end +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl new file mode 100644 index 0000000000..e24110b0f9 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl @@ -0,0 +1,38 @@ +module SimpleSparsePolynomialZonotopeModule + +using Reexport + +using ..LazySets: AbstractPolynomialZonotope +using Random: AbstractRNG, GLOBAL_RNG +using ReachabilityBase.Distribution: reseed! +using ReachabilityBase.Comparison: isapproxzero +using LinearAlgebra: dot + +@reexport import ..API: convex_hull, center, isoperationtype, rand, linear_map +@reexport import ..LazySets: expmat, genmat, ngens, nparams, polynomial_order, + remove_redundant_generators +@reexport using ..API + +export SimpleSparsePolynomialZonotope, + SSPZ, + quadratic_map + +include("SimpleSparsePolynomialZonotope.jl") + +# short-hand +const SSPZ = SimpleSparsePolynomialZonotope + +include("center.jl") +include("convex_hull.jl") +include("genmat.jl") +include("isoperationtype.jl") +include("ngens.jl") +include("polynomial_order.jl") +include("remove_redundant_generators.jl") +include("expmat.jl") +include("nparams.jl") +include("quadratic_map.jl") +include("rand.jl") +include("linear_map.jl") + +end # module diff --git a/src/Sets/SimpleSparsePolynomialZonotope/center.jl b/src/Sets/SimpleSparsePolynomialZonotope/center.jl new file mode 100644 index 0000000000..9110debec3 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/center.jl @@ -0,0 +1,14 @@ +""" + center(P::SimpleSparsePolynomialZonotope) + +Return the center of a simple sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The center of `P`. +""" +center(P::SSPZ) = P.c diff --git a/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl b/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl new file mode 100644 index 0000000000..3914715e42 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl @@ -0,0 +1,16 @@ +""" + convex_hull(P::SimpleSparsePolynomialZonotope) + +Compute the convex hull of a simple sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The tightest convex simple sparse polynomial zonotope containing `P`. +""" +function convex_hull(P::SimpleSparsePolynomialZonotope) + return linear_combination(P, P) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/expmat.jl b/src/Sets/SimpleSparsePolynomialZonotope/expmat.jl new file mode 100644 index 0000000000..dfa72ec655 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/expmat.jl @@ -0,0 +1,31 @@ +""" + expmat(P::SimpleSparsePolynomialZonotope) + +Return the matrix of exponents of a simple sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The matrix of exponents, where each column is a multidegree. + +### Notes + +In the exponent matrix, each row corresponds to a parameter (``\alpha_k`` in the +mathematical set definition) and each column corresponds to a monomial. + +### Examples + +```jldoctest +julia> S = SimpleSparsePolynomialZonotope([2.0, 0], [1 2;2 2.], [1 4;1 2]) +SimpleSparsePolynomialZonotope{Float64, Vector{Float64}, Matrix{Float64}, Matrix{Int64}}([2.0, 0.0], [1.0 2.0; 2.0 2.0], [1 4; 1 2]) + +julia> expmat(S) +2×2 Matrix{Int64}: + 1 4 + 1 2 +``` +""" +expmat(P::SSPZ) = P.E diff --git a/src/Sets/SimpleSparsePolynomialZonotope/genmat.jl b/src/Sets/SimpleSparsePolynomialZonotope/genmat.jl new file mode 100644 index 0000000000..d61cb7ed45 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/genmat.jl @@ -0,0 +1,14 @@ +""" + genmat(P::SimpleSparsePolynomialZonotope) + +Return the matrix of generators of a simple sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The matrix of generators of `P`. +""" +genmat(P::SSPZ) = P.G diff --git a/src/Sets/SimpleSparsePolynomialZonotope/isoperationtype.jl b/src/Sets/SimpleSparsePolynomialZonotope/isoperationtype.jl new file mode 100644 index 0000000000..56205dd39a --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/isoperationtype.jl @@ -0,0 +1,3 @@ +function isoperationtype(::Type{<:SimpleSparsePolynomialZonotope}) + return false +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/linear_map.jl b/src/Sets/SimpleSparsePolynomialZonotope/linear_map.jl new file mode 100644 index 0000000000..f7b6a6e3dc --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/linear_map.jl @@ -0,0 +1,17 @@ +""" + linear_map(M::AbstractMatrix, P::SimpleSparsePolynomialZonotope) + +Apply the linear map `M` to a simple sparse polynomial zonotope. + +### Input + +- `M` -- matrix +- `P` -- simple sparse polynomial zonotope + +### Output + +The set resulting from applying the linear map `M` to `P`. +""" +function linear_map(M::AbstractMatrix, P::SSPZ) + return SimpleSparsePolynomialZonotope(M * center(P), M * genmat(P), expmat(P)) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/ngens.jl b/src/Sets/SimpleSparsePolynomialZonotope/ngens.jl new file mode 100644 index 0000000000..33a848285f --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/ngens.jl @@ -0,0 +1,19 @@ +""" + ngens(P::SimpleSparsePolynomialZonotope) + +Return the number of generators of a simple sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The number of generators of `P`. + +### Notes + +This number corresponds to the number of monomials in the polynomial +representation of `P`. +""" +ngens(P::SSPZ) = size(P.G, 2) diff --git a/src/Sets/SimpleSparsePolynomialZonotope/nparams.jl b/src/Sets/SimpleSparsePolynomialZonotope/nparams.jl new file mode 100644 index 0000000000..ffa6d6ad83 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/nparams.jl @@ -0,0 +1,30 @@ +""" + nparams(P::SimpleSparsePolynomialZonotope) + +Return the number of parameters in the polynomial representation of a simple +sparse polynomial zonotope. + +### Input + +- `P` -- simple sparse polynomial zonotope + +### Output + +The number of parameters in the polynomial representation of P. + +### Notes + +This number corresponds to the number of rows in the exponent matrix ``E`` (`p` +in the mathematical set definition). + +### Examples + +```jldoctest +julia> S = SimpleSparsePolynomialZonotope([2.0, 0], [1 2;2 2.], [1 4;1 2]) +SimpleSparsePolynomialZonotope{Float64, Vector{Float64}, Matrix{Float64}, Matrix{Int64}}([2.0, 0.0], [1.0 2.0; 2.0 2.0], [1 4; 1 2]) + +julia> nparams(S) +2 +``` +""" +nparams(P::SSPZ) = size(P.E, 1) diff --git a/src/Sets/SimpleSparsePolynomialZonotope/polynomial_order.jl b/src/Sets/SimpleSparsePolynomialZonotope/polynomial_order.jl new file mode 100644 index 0000000000..04ea03787c --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/polynomial_order.jl @@ -0,0 +1,3 @@ +function polynomial_order(P::SSPZ) + return maximum(sum, eachcol(expmat(P))) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/quadratic_map.jl b/src/Sets/SimpleSparsePolynomialZonotope/quadratic_map.jl new file mode 100644 index 0000000000..908fec81c0 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/quadratic_map.jl @@ -0,0 +1,108 @@ +""" + quadratic_map(Q::Vector{MT}, S::SimpleSparsePolynomialZonotope) + where {N, MT<:AbstractMatrix{N}} + +Return the quadratic map of a simple sparse polynomial zonotope. + +### Input + +- `Q` -- vector of square matrices +- `S` -- simple sparse polynomial zonotope + +### Output + +The quadratic map of `P` represented as a simple sparse polynomial zonotope. + +### Algorithm + +This method implements Proposition 12 in [1]. +See also Proposition 3.1.30 in [2]. + +[1] N. Kochdumper, M. Althoff. *Sparse polynomial zonotopes: A novel set +representation for reachability analysis*. 2021 +[2] N. Kochdumper. *Extensions of polynomial zonotopes and their application to +verification of cyber-physical systems*. 2021. +""" +function quadratic_map(Q::Vector{MT}, + S::SimpleSparsePolynomialZonotope) where {N,MT<:AbstractMatrix{N}} + m = length(Q) + c = center(S) + h = ngens(S) + G = genmat(S) + E = expmat(S) + + cnew = similar(c, m) + Gnew = similar(G, m, h^2 + h) + QiG = similar(Q) + @inbounds for (i, Qi) in enumerate(Q) + cnew[i] = dot(c, Qi, c) + Gnew[i, 1:h] = c' * (Qi + Qi') * G + QiG[i] = Qi * G + end + + Enew = repeat(E, 1, h + 1) + @inbounds for i in 1:h + idxstart = h * i + 1 + idxend = (i + 1) * h + Enew[:, idxstart:idxend] .+= E[:, i] + for j in eachindex(QiG) + Gnew[j, idxstart:idxend] = G[:, i]' * QiG[j] + end + end + Z = SimpleSparsePolynomialZonotope(cnew, Gnew, Enew) + return remove_redundant_generators(Z) +end + +""" + quadratic_map(Q::Vector{MT}, S1::SimpleSparsePolynomialZonotope, + S2::SimpleSparsePolynomialZonotope) + where {N, MT<:AbstractMatrix{N}} + +Return the quadratic map of two simple sparse polynomial zonotopes. +The quadratic map is the set +```math + \\{x \\mid xᵢ = s₁ᵀQᵢs₂, s₁ ∈ S₁, s₂ ∈ S₂, Qᵢ ∈ Q\\}. +``` + +### Input + +- `Q` -- vector of square matrices +- `S1` -- simple sparse polynomial zonotope +- `S2` -- simple sparse polynomial zonotope + +### Output + +The quadratic map of the given simple sparse polynomial zonotopes represented as +a simple sparse polynomial zonotope. + +### Algorithm + +This method implements Proposition 3.1.30 in [1]. + +[1] N. Kochdumper. *Extensions of polynomial zonotopes and their application to +verification of cyber-physical systems*. 2021. +""" +function quadratic_map(Q::Vector{MT}, S1::SimpleSparsePolynomialZonotope, + S2::SimpleSparsePolynomialZonotope) where {N,MT<:AbstractMatrix{N}} + @assert nparams(S1) == nparams(S2) + + c1 = center(S1) + c2 = center(S2) + G1 = genmat(S1) + G2 = genmat(S2) + E1 = expmat(S1) + E2 = expmat(S2) + + c = [dot(c1, Qi, c2) for Qi in Q] + + Ghat1 = reduce(vcat, c2' * Qi' * G1 for Qi in Q) + Ghat2 = reduce(vcat, c1' * Qi * G2 for Qi in Q) + + Gbar = reduce(hcat, reduce(vcat, gj' * Qi * G2 for Qi in Q) for gj in eachcol(G1)) + Ebar = reduce(hcat, E2 .+ e1j for e1j in eachcol(E1)) + + G = hcat(Ghat1, Ghat2, Gbar) + E = hcat(E1, E2, Ebar) + + return remove_redundant_generators(SimpleSparsePolynomialZonotope(c, G, E)) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/rand.jl b/src/Sets/SimpleSparsePolynomialZonotope/rand.jl new file mode 100644 index 0000000000..d2b5e0704d --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/rand.jl @@ -0,0 +1,51 @@ +""" + rand(::Type{SimpleSparsePolynomialZonotope}; + [N]::Type{<:Real}=Float64, [dim]::Int=2, [nparams]::Int=2, + [maxdeg]::Int=3, [num_generators]::Int=-1, + [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) + +Create a random simple sparse polynomial zonotope. + +### Input + +- `Zonotope` -- type for dispatch +- `N` -- (optional, default: `Float64`) numeric type +- `dim` -- (optional, default: 2) dimension +- `nparams` -- (optional, default: 2) number of parameters +- `maxdeg` -- (optional, default: 3) maximum degree for each parameter +- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator +- `seed` -- (optional, default: `nothing`) seed for reseeding +- `num_generators` -- (optional, default: `-1`) number of generators of the + zonotope (see comment below) + +### Output + +A random simple sparse polynomial zonotope. + +### Algorithm + +All numbers are normally distributed with mean 0 and standard deviation 1. + +The number of generators can be controlled with the argument `num_generators`. +For a negative value we choose a random number in the range `dim:2*dim` (except +if `dim == 1`, in which case we only create a single generator). Note that the +final number of generators may be lower if redundant monomials are generated. +""" +function rand(::Type{SimpleSparsePolynomialZonotope}; + N::Type{<:Real}=Float64, + dim::Int=2, + nparams::Int=2, + maxdeg::Int=3, + rng::AbstractRNG=GLOBAL_RNG, + seed::Union{Int,Nothing}=nothing, + num_generators::Int=-1) + rng = reseed!(rng, seed) + center = randn(rng, N, dim) + if num_generators < 0 + num_generators = (dim == 1) ? 1 : rand(rng, dim:(2 * dim)) + end + generators = randn(rng, N, dim, num_generators) + expmat = rand(rng, 0:maxdeg, nparams, num_generators) + SSPZ = SimpleSparsePolynomialZonotope(center, generators, expmat) + return remove_redundant_generators(SSPZ) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/remove_redundant_generators.jl b/src/Sets/SimpleSparsePolynomialZonotope/remove_redundant_generators.jl new file mode 100644 index 0000000000..ebe230281b --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/remove_redundant_generators.jl @@ -0,0 +1,57 @@ +""" + remove_redundant_generators(S::SimpleSparsePolynomialZonotope) + +Remove redundant generators from a simple sparse polynomial zonotope. + +### Input + +- `S` -- simple sparse polynomial zonotope + +### Output + +A new simple sparse polynomial zonotope such that redundant generators have been +removed. + +## Notes + +The result uses dense arrays irrespective of the array type of `S`. + +### Algorithm + +Let `G` be the generator matrix and `E` the exponent matrix of `S`. The +following simplifications are performed: + +- Zero columns in `G` and the corresponding columns in `E` are removed. +- For zero columns in `E`, the corresponding column in `G` is summed to the + center. +- Repeated columns in `E` are grouped together by summing the corresponding + columns in `G`. +""" +function remove_redundant_generators(S::SimpleSparsePolynomialZonotope) + c, G, E = _remove_redundant_generators_polyzono(center(S), genmat(S), expmat(S)) + + return SimpleSparsePolynomialZonotope(c, G, E) +end + +function _remove_redundant_generators_polyzono(c, G, E) + Gnew = Matrix{eltype(G)}(undef, size(G, 1), 0) + Enew = Matrix{eltype(E)}(undef, size(E, 1), 0) + cnew = copy(c) + + visited_exps = Dict{Vector{Int},Int}() + @inbounds for (gi, ei) in zip(eachcol(G), eachcol(E)) + all(isapproxzero, gi) && continue + if iszero(ei) + cnew += gi + elseif haskey(visited_exps, ei) # repeated exponent + idx = visited_exps[ei] + Gnew[:, idx] += gi + else + Gnew = hcat(Gnew, gi) + Enew = hcat(Enew, ei) + visited_exps[ei] = size(Enew, 2) + end + end + + return cnew, Gnew, Enew +end diff --git a/test/Sets/SimpleSparsePolynomialZonotope.jl b/test/Sets/SimpleSparsePolynomialZonotope.jl index a31ee85622..e6ff34ab5b 100644 --- a/test/Sets/SimpleSparsePolynomialZonotope.jl +++ b/test/Sets/SimpleSparsePolynomialZonotope.jl @@ -14,6 +14,7 @@ for N in [Float64, Float32, Rational{Int}] @test ngens(S) == 2 @test nparams(S) == 2 @test order(S) == 1 // 1 + @test polynomial_order(S) == 6 Z = Zonotope(N[3.0, 1], N[1 1; 2 1.0]) @test overapproximate(S, Zonotope) == Z @@ -96,6 +97,9 @@ for N in [Float64, Float32, Rational{Int}] @test expmat(Sred) == [4 1 3 7 5 2 4 8 6 3 5 9] + + # isoperationtype + @test !isoperationtype(SimpleSparsePolynomialZonotope) end for Z in [rand(Zonotope), rand(Hyperrectangle)]