diff --git a/docs/make.jl b/docs/make.jl index 611292db90..c20e2a00d5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -88,7 +88,8 @@ makedocs(; sitename="LazySets.jl", "VPolygon" => "lib/sets/VPolygon.md", "VPolytope" => "lib/sets/VPolytope.md", "ZeroSet" => "lib/sets/ZeroSet.md", - "Zonotope" => "lib/sets/Zonotope.md" + "Zonotope" => "lib/sets/Zonotope.md", + "ZonotopeMD" => "lib/sets/ZonotopeMD.md" # ], "Lazy Operations" => [ diff --git a/docs/src/lib/sets/ZonotopeMD.md b/docs/src/lib/sets/ZonotopeMD.md new file mode 100644 index 0000000000..2325caf9d7 --- /dev/null +++ b/docs/src/lib/sets/ZonotopeMD.md @@ -0,0 +1,27 @@ +```@meta +CurrentModule = LazySets.ZonotopeModule +``` + +# [ZonotopeMD](@id def_ZonotopeMD) + +```@docs +ZonotopeMD +``` + +## Conversion + +```julia +Zonotope(Z::ZonotopeMD) +``` + +## Operations + +```@docs +genmat(::ZonotopeMD) +cartesian_product(::ZonotopeMD, ::ZonotopeMD) +``` + +Undocumented implementations: + +* [`center`](@ref center) +* [`ngens`](@ref ngens) diff --git a/src/Interfaces/AbstractZonotope.jl b/src/Interfaces/AbstractZonotope.jl index 2027472e7e..16cac305bd 100644 --- a/src/Interfaces/AbstractZonotope.jl +++ b/src/Interfaces/AbstractZonotope.jl @@ -42,11 +42,12 @@ The subtypes of `AbstractZonotope` (including abstract interfaces): ```jldoctest; setup = :(using LazySets: subtypes) julia> subtypes(AbstractZonotope) -4-element Vector{Any}: +5-element Vector{Any}: AbstractHyperrectangle HParallelotope LineSegment Zonotope + ZonotopeMD ``` """ abstract type AbstractZonotope{N} <: AbstractCentrallySymmetricPolytope{N} end diff --git a/src/Interfaces/LazySet.jl b/src/Interfaces/LazySet.jl index 03b3b58533..9d1b958a74 100644 --- a/src/Interfaces/LazySet.jl +++ b/src/Interfaces/LazySet.jl @@ -64,7 +64,7 @@ If we only consider *concrete* subtypes, then: julia> concrete_subtypes = subtypes(LazySet, true); julia> length(concrete_subtypes) -53 +54 julia> println.(concrete_subtypes); AffineMap @@ -120,7 +120,7 @@ VPolygon VPolytope ZeroSet Zonotope - +ZonotopeMD ``` """ abstract type LazySet{N} end diff --git a/src/LazySets.jl b/src/LazySets.jl index dd6043cda5..e8f0cc09d5 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -218,7 +218,8 @@ include("Sets/Zonotope/ZonotopeModule.jl") @reexport using ..ZonotopeModule: Zonotope, remove_zero_generators, linear_map!, - split! + split!, + ZonotopeMD using ..ZonotopeModule: _split include("LazyOperations/UnionSet.jl") # must come before IntervalModule diff --git a/src/Sets/Zonotope/ZonotopeMD/ZonotopeMD.jl b/src/Sets/Zonotope/ZonotopeMD/ZonotopeMD.jl new file mode 100644 index 0000000000..1c4e22be2b --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/ZonotopeMD.jl @@ -0,0 +1,102 @@ +""" + struct ZonotopeMD{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}, DN<:AbstractVector{N}} <: AbstractZonotope{N} + +Type that represents a zonotope of order `k` in normal form. + +### Fields + +- `center::VN` — the center of the zonotope +- `M::MN` — matrix of general (non-axis-aligned) generators +- `d::DN` — vector of axis-aligned (diagonal) generators + +### Notes +A zonotope is of order `k` if it has `n * k` generators in `ℝⁿ`, where `n` is the ambient dimension. + +A zonotope of order `k` in *normal form* is defined as the set + +```math +Z = \\left\\{ x ∈ ℝ^n : x = c + Mξ + d ⊙ η, ~~ ξ ∈ [-1, 1]^m, ~~ η ∈ [-1, 1]^n \\right\\}, +``` + +where `M ∈ ℝ^{n×m}` is a matrix of general generators with `m = n*(k -1)` and `d ∈ ℝⁿ` is a vector of axis-aligned generators. +Equivalently, this can be seen as a zonotope with generator matrix `[M D]`, where `D` is the diagonal matrix +formed from the vector `d`. +ZonotopeMD can be constructed in two ways: by passing the full generator matrix `[M D]` in normal form +or by passing `M` and a vector `d` separately. + +### Examples + +Constructing a zonotope in normal form from a center, general generator matrix `M`, and diagonal vector `d`: + +```jldoctest zonotopemd_label +julia> c = [0.0, 0.0]; + +julia> M = [1.0 2.0; 3.0 1.0]; + +julia> d = [0.1, 0.2]; + +julia> Z = ZonotopeMD(c, M, d) +ZonotopeMD{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([0.0, 0.0], [1.0 2.0; 3.0 1.0], [0.1, 0.2]) + +julia> center(Z) +2-element Vector{Float64}: + 0.0 + 0.0 + +julia> genmat(Z) +2×4 SparseArrays.SparseMatrixCSC{Float64, Int64} with 6 stored entries: + 1.0 2.0 0.1 ⋅ + 3.0 1.0 ⋅ 0.2 +``` + +The generator matrix returned by `genmat` is the concatenation `[M D]`, where `D` is the diagonal matrix formed from `d`. +THe resulting matrix is stored as a sparse matrix. +Constructing the same zonotope by passing the full generator matrix `[M D]` directly: + +```jldoctest zonotopemd_label +julia> G = [1.0 2.0 0.1 0.0; + 3.0 1.0 0.0 0.2]; + +julia> Z2 = ZonotopeMD([0.0, 0.0], G) +ZonotopeMD{Float64, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([0.0, 0.0], [1.0 2.0; 3.0 1.0], [0.1, 0.2]) + +julia> genmat(Z2) == G +true +``` +You can also convert back to a standard `Zonotope` if needed: + +```jldoctest zonotopemd_label +julia> Zstd = Zonotope(Z) +Zonotope{Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}([0.0, 0.0], sparse([1, 2, 1, 2, 1, 2], [1, 1, 2, 2, 3, 4], [1.0, 3.0, 2.0, 1.0, 0.1, 0.2], 2, 4)) +``` + +""" +struct ZonotopeMD{N,VN<:AbstractVector{N},MN<:AbstractMatrix{N},DN<:AbstractVector{N}} <: + AbstractZonotope{N} + center::VN + M::MN + d::DN + + function ZonotopeMD(center::VN, M::MN, + d::DN) where {N, + VN<:AbstractVector{N}, + MN<:AbstractMatrix{N}, + DN<:AbstractVector{N}} + @assert length(center) == size(M, 1) == length(d) "Dimensions must match" + return new{N,VN,MN,DN}(center, M, d) + end +end + +# constructor from generator matrix +function ZonotopeMD(center::VN, G::AbstractMatrix{N}) where {N,VN<:AbstractVector{N}} + n, p = size(G) + @assert p % n == 0 "The generator matrix must contain a multiple of n columns" + @assert p >= 2n "Expected at least order 2 zonotope" + + M = G[:, 1:(p - n)] + D = G[:, (end - n + 1):end] + + @assert isdiag(D) "The last n columns of the generator matrix must be diagonal" + d = diag(D) + return ZonotopeMD(center, M, d) +end diff --git a/src/Sets/Zonotope/ZonotopeMD/cartesian_product.jl b/src/Sets/Zonotope/ZonotopeMD/cartesian_product.jl new file mode 100644 index 0000000000..ed0a78e8e9 --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/cartesian_product.jl @@ -0,0 +1,20 @@ + +""" + cartesian_product(Z1::ZonotopeMD, Z2::ZonotopeMD) + +Return the Cartesian product of two zonotopes in normal form (`ZonotopeMD`). + +### Input + +- `Z1`, `Z2` -- zonotopes in normal form (`ZonotopeMD`) + +### Output + +A new `ZonotopeMD` representing the Cartesian product `Z1 × Z2`. +""" +function cartesian_product(Z1::ZonotopeMD, Z2::ZonotopeMD) + c = vcat(Z1.center, Z2.center) + d = vcat(Z1.d, Z2.d) + M = blockdiag(sparse(Z1.M), sparse(Z2.M)) + return ZonotopeMD(c, M, d) +end diff --git a/src/Sets/Zonotope/ZonotopeMD/center.jl b/src/Sets/Zonotope/ZonotopeMD/center.jl new file mode 100644 index 0000000000..f3f476d73e --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/center.jl @@ -0,0 +1,3 @@ +function center(Z::ZonotopeMD) + return Z.center +end diff --git a/src/Sets/Zonotope/ZonotopeMD/convert.jl b/src/Sets/Zonotope/ZonotopeMD/convert.jl new file mode 100644 index 0000000000..41e7c84c8e --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/convert.jl @@ -0,0 +1 @@ +Zonotope(Z::ZonotopeMD) = convert(Zonotope, Z) diff --git a/src/Sets/Zonotope/ZonotopeMD/genmat.jl b/src/Sets/Zonotope/ZonotopeMD/genmat.jl new file mode 100644 index 0000000000..b966d29bc7 --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/genmat.jl @@ -0,0 +1,17 @@ +""" + genmat(Z::ZonotopeMD) + +Return the generator matrix of a ZonotopeMD. + +### Input + +- `Z` -- zonotope in normal form + +### Output + +A matrix where each column represents one generator of the zonotope `Z`. +""" +function genmat(Z::ZonotopeMD) + D = spdiagm(0 => Z.d) + return hcat(Z.M, D) +end diff --git a/src/Sets/Zonotope/ZonotopeMD/isoperationtype.jl b/src/Sets/Zonotope/ZonotopeMD/isoperationtype.jl new file mode 100644 index 0000000000..5e7bd5ee64 --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/isoperationtype.jl @@ -0,0 +1 @@ +isoperationtype(::Type{<:ZonotopeMD}) = false \ No newline at end of file diff --git a/src/Sets/Zonotope/ZonotopeMD/ngens.jl b/src/Sets/Zonotope/ZonotopeMD/ngens.jl new file mode 100644 index 0000000000..a6d74cc7fe --- /dev/null +++ b/src/Sets/Zonotope/ZonotopeMD/ngens.jl @@ -0,0 +1 @@ +ngens(Z::ZonotopeMD) = size(Z.M, 2) + length(Z.d) diff --git a/src/Sets/Zonotope/ZonotopeModule.jl b/src/Sets/Zonotope/ZonotopeModule.jl index 0bd767a04f..75a2ceb6a5 100644 --- a/src/Sets/Zonotope/ZonotopeModule.jl +++ b/src/Sets/Zonotope/ZonotopeModule.jl @@ -3,7 +3,8 @@ module ZonotopeModule using Reexport, Requires using ..LazySets: AbstractZonotope, generators_fallback, _scale_copy_inplace -using LinearAlgebra: mul! +using LinearAlgebra: mul!, Diagonal, isdiag, diag +using SparseArrays: blockdiag, sparse, spdiagm using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Arrays: ismultiple, remove_zero_columns, to_matrix, vector_type @@ -11,7 +12,7 @@ using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Require: require @reexport import ..API: center, high, isoperationtype, low, rand, - permute, scale, scale!, translate! + permute, scale, scale!, translate!, cartesian_product @reexport import ..LazySets: generators, genmat, ngens, reduce_order, remove_redundant_generators, togrep import Base: convert @@ -20,7 +21,8 @@ import Base: convert export Zonotope, remove_zero_generators, linear_map!, - split! + split!, + ZonotopeMD include("Zonotope.jl") @@ -47,4 +49,13 @@ include("convert.jl") include("init.jl") +#ZonotopeMD +include("ZonotopeMD/ZonotopeMD.jl") +include("ZonotopeMD/convert.jl") +include("ZonotopeMD/center.jl") +include("ZonotopeMD/genmat.jl") +include("ZonotopeMD/cartesian_product.jl") +include("ZonotopeMD/ngens.jl") +include("ZonotopeMD/isoperationtype.jl") + end # module diff --git a/test/Sets/ZonotopeMD.jl b/test/Sets/ZonotopeMD.jl new file mode 100644 index 0000000000..f07893d393 --- /dev/null +++ b/test/Sets/ZonotopeMD.jl @@ -0,0 +1,96 @@ +for N in [Float64, Float32, Rational{Int}] + # For order 2: + # center ∈ ℝ², + # M is 2×2 and d ∈ ℝ² + + # Direct construction via (center, M, d) + c = [N(1), N(2)] + M = [N(1) N(0); N(0) N(1)] + d = [N(1//10), N(2)] + Zmd = ZonotopeMD(c, M, d) + @test Zmd isa ZonotopeMD{N} + @test center(Zmd) == c + @test genmat(Zmd) == hcat(M, Diagonal(d)) + @test length(Zmd.center) == size(Zmd.M, 1) == length(Zmd.d) + + @test length(center(Zmd)) == 2 + @test genmat(Zmd) !== nothing + + # Construction from generator matrix + G = hcat(M, Diagonal(d)) + Zmd2 = ZonotopeMD(c, G) + @test Zmd2 == Zmd + + # Generator matrix with wrong shape + G_wrong = hcat(M, ones(N, 2, 1)) + @test_throws AssertionError ZonotopeMD(c, G_wrong) + + # Convert to standard Zonotope + Zstd = convert(Zonotope, Zmd) + @test center(Zstd) == c + @test genmat(Zstd) == G + + # Cartesian product + c2 = [N(3), N(4)] + M2 = [N(2) N(0); N(0) N(2)] + d2 = [N(3), N(4)] + Zmd_2 = ZonotopeMD(c2, M2, d2) + cp_md = cartesian_product(Zmd, Zmd_2) + cp_std = cartesian_product(Zstd, Zonotope(Zmd_2)) + @test isequivalent(cp_md, cp_std) + + # Apply a linear map + A = [N(1) N(1); N(0) N(1)] + Zlm = linear_map(A, Zmd) + expected_c = A * c + expected_genmat = A * hcat(M, Diagonal(d)) + @test center(Zlm) == expected_c + @test genmat(Zlm) == expected_genmat + + # For order 3: + # center ∈ ℝ², + # M is 2×4 and d ∈ ℝ². + + # Direct construction via (center, M, d) + c = [N(1), N(2)] + M1 = [N(1) N(0) N(2) N(0); N(0) N(1) N(3) N(0)] + d = [N(1//10), N(2)] + Zmd = ZonotopeMD(c, M1, d) + @test Zmd isa ZonotopeMD{N} + @test center(Zmd) == c + @test genmat(Zmd) == hcat(M1, Diagonal(d)) + @test length(Zmd.center) == size(Zmd.M, 1) == length(Zmd.d) + + @test length(center(Zmd)) == 2 + @test genmat(Zmd) !== nothing + + # Construction from generator matrix + G = hcat(M1, Diagonal(d)) + Zmd2 = ZonotopeMD(c, G) + @test Zmd2 == Zmd + + G_wrong = hcat(M1, ones(N, 2, 1)) # here p = 3 instead of 4 (2n) + @test_throws AssertionError ZonotopeMD(c, G_wrong) + + # Convert to standard Zonotope + Zstd = convert(Zonotope, Zmd) + @test center(Zstd) == c + @test genmat(Zstd) == G + + # Cartesian product + c2 = [N(3), N(4)] + M2 = [N(2) N(0) N(1) N(0); N(0) N(2) N(3) N(0)] + d2 = [N(3), N(4)] + Zmd_2 = ZonotopeMD(c2, M2, d2) + cp_md = cartesian_product(Zmd, Zmd_2) + cp_std = cartesian_product(Zstd, Zonotope(Zmd_2)) + @test isequivalent(cp_md, cp_std) + + # Apply a linear map + A = [N(1) N(1); N(0) N(1)] + Zlm = linear_map(A, Zmd) + expected_c = A * c + expected_genmat = A * hcat(M1, Diagonal(d)) + @test center(Zlm) == expected_c + @test genmat(Zlm) == expected_genmat +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6e06d05449..a862542270 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -150,6 +150,9 @@ if test_suite_basic @testset "LazySets.Zonotope" begin include("Sets/Zonotope.jl") end + @testset "LazySets.ZonotopeMD" begin + include("Sets/ZonotopeMD.jl") + end @testset "LazySets.ZeroSet" begin include("Sets/ZeroSet.jl") end