Merge pull request #1228 from JuliaReach/mforets/1226

#1226 - Convert cartesian product of zonotopes to a new zonotope
JuliaReach · Mar 12, 2019 · 4173669 · 4173669
2 parents 1225724 + 89d9d1a
commit 4173669
Show file tree

Hide file tree

Showing 4 changed files with 38 additions and 0 deletions.
diff --git a/docs/src/lib/conversion.md b/docs/src/lib/conversion.md
@@ -39,4 +39,5 @@ convert(::Type{VPolytope}, ::HPolytope)
 convert(::Type{Zonotope}, ::AbstractHyperrectangle)
 convert(::Type{IntervalArithmetic.IntervalBox}, ::AbstractHyperrectangle)
 convert(::Type{Hyperrectangle}, ::IntervalArithmetic.IntervalBox)
+convert(::Type{Zonotope}, ::CartesianProduct{N, Zonotope{N}, Zonotope{N}}) where {N<:Real}
 ```
diff --git a/src/compat.jl b/src/compat.jl
@@ -19,6 +19,7 @@ export ×
     @inline _At_mul_B(A, B) = At_mul_B(A, B)
     @inline _At_ldiv_B(A, B) = At_ldiv_B(A, B)
     expmat = expm
+    blockdiag = Base.SparseArrays.blkdiag
 else
     using SparseArrays
     using Random

diff --git a/src/convert.jl b/src/convert.jl
@@ -455,6 +455,36 @@ function convert(::Type{Hyperrectangle},
      return Hyperrectangle(low=l, high=h)
 end
 
+"""
+    convert(::Type{Zonotope}, cp::CartesianProduct{N, Zonotope{N}, Zonotope{N}}) where {N<:Real}
+
+Converts the cartesian product of two zonotopes to a new zonotope.
+
+### Input
+
+- `Zonotope` -- type used for dispatch
+- `S`        -- cartesian product of two zonotopes
+
+### Output
+
+A zonotope.
+
+### Algorithm
+
+The cartesian product is obtained by:
+
+- Concatenating the centers of each input zonotope.
+- Arranging the generators in block-diagional fashion, and filled with zeros
+  in the off-diagonal; for this reason, the generator matrix of the returned
+  zonotope is built as a sparse matrix.
+"""
+function convert(::Type{Zonotope}, cp::CartesianProduct{N, Zonotope{N}, Zonotope{N}}) where {N<:Real}
+    Z1, Z2 = cp.X, cp.Y
+    c = vcat(Z1.center, Z2.center)
+    G = blockdiag(sparse(Z1.generators), sparse(Z2.generators))
+    return Zonotope(c, G)
+end
+
 """
     convert(::Type{HPolytope}, H::AbstractHyperrectangle)
 

diff --git a/test/unit_Zonotope.jl b/test/unit_Zonotope.jl
@@ -106,6 +106,12 @@ for N in [Float64, Rational{Int}, Float32]
     Z1, Z2 = split(Z, 1) # in this case the splitting is exact
     @test Z1 ⊆ Z && Z2 ⊆ Z
 
+    # converts the cartesian product of two zonotopes to a new zonotope
+    Z1 = Zonotope(N[0], hcat(N[1]))
+    Z2 = Zonotope(N[1/2], hcat(N[1/2]))
+    Z = convert(Zonotope, Z1×Z2)
+    @test Z isa Zonotope && Z.center == N[0, 1/2] && Matrix(Z.generators) == N[1 0; 0 1/2]
+
     # list of constraints
     Z = Zonotope(zeros(N, 3), Matrix(N(1)*I, 3, 3))
     B = BallInf(zeros(N, 3), N(1)) # equivalent to Z