#2739: added ZonotopeMD and support for Cartesian product in normal

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" => [

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)

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

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

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

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

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

src/Sets/Zonotope/ZonotopeMD/center.jl

@@ -0,0 +1,3 @@
+function center(Z::ZonotopeMD)
+    return Z.center
+end

src/Sets/Zonotope/ZonotopeMD/convert.jl

@@ -0,0 +1 @@
+Zonotope(Z::ZonotopeMD) = convert(Zonotope, Z)

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

src/Sets/Zonotope/ZonotopeMD/isoperationtype.jl

@@ -0,0 +1 @@
+isoperationtype(::Type{<:ZonotopeMD}) = false

src/Sets/Zonotope/ZonotopeMD/ngens.jl

@@ -0,0 +1 @@
+ngens(Z::ZonotopeMD) = size(Z.M, 2) + length(Z.d)

src/Sets/Zonotope/ZonotopeModule.jl

@@ -3,15 +3,16 @@ 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
 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

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

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