#416 - Conversion from Zonotope to HPolytope

docs/src/lib/representations.md

@@ -583,6 +583,7 @@ Zonotope
 rand(::Type{Zonotope})
 vertices_list(::Zonotope{N}) where {N<:Real}
 constraints_list(::Zonotope{N}) where {N<:Real}
+constraints_list(::Zonotope{N}) where {N<:AbstractFloat}
 center(::Zonotope{N}) where {N<:Real}
 order(::Zonotope)
 minkowski_sum(::Zonotope{N}, ::Zonotope{N}) where {N<:Real}

docs/src/lib/utils.md

@@ -23,6 +23,7 @@ ispermutation
 @neutral_absorbing
 @array_neutral
 @array_absorbing
+cross_product(::AbstractMatrix{N}) where {N<:Real}
 get_radius!
 an_element_helper
 σ_helper
@@ -39,4 +40,5 @@ reseed
 ```@docs
 CachedPair
 Approximations.UnitVector
+StrictlyIncreasingIndices
 ```

src/Zonotope.jl

@@ -439,7 +439,8 @@ function reduce_order(Z::Zonotope{N}, r)::Zonotope{N} where {N<:Real}
 end
 
 """
-    constraints_list(P::Zonotope{N})::Vector{LinearConstraint{N}} where {N<:Real}
+    constraints_list(P::Zonotope{N}
+                    )::Vector{LinearConstraint{N}} where {N<:Real}
 
 Return the list of constraints defining a zonotope.
 
@@ -449,14 +450,90 @@ Return the list of constraints defining a zonotope.
 
 ### Output
 
-The list of constraints of the polyhedron.
+The list of constraints of the zonotope.
 
 ### Algorithm
 
-This is a naive implementation that calculates all vertices and transforms
-to the H-representation of the zonotope. The transformation to the dual
-representation requires the concrete polyhedra package `Polyhedra`.
+This is the (inefficient) fallback implementation for rational numbers.
+It first computes the vertices and then converts the corresponding polytope
+to constraint representation.
 """
-function constraints_list(Z::Zonotope{N})::Vector{LinearConstraint{N}} where {N<:Real}
-    return constraints_list(convert(HPolytope, Z))
+function constraints_list(Z::Zonotope{N}
+                         )::Vector{LinearConstraint{N}} where {N<:Real}
+    return constraints_list(VPolytope(vertices_list(Z)))
+end
+
+"""
+    constraints_list(P::Zonotope{N}
+                    )::Vector{LinearConstraint{N}} where {N<:AbstractFloat}
+
+Return the list of constraints defining a zonotope.
+
+### Input
+
+- `Z` -- zonotope
+
+### Output
+
+The list of constraints of the zonotope.
+
+### Notes
+
+The algorithm assumes that no generator is redundant.
+The result has ``2 \\binom{p}{n-1}`` (with ``p`` being the number of generators
+and ``n`` being the ambient dimension) constraints, which is optimal under this
+assumption.
+
+### Algorithm
+
+We follow the algorithm presented in *Althoff, Stursberg, Buss: Computing
+Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes.
+2009.*
+
+The one-dimensional case is not covered by that algorithm; we manually handle
+this case, assuming that there is only one generator.
+"""
+function constraints_list(Z::Zonotope{N}
+                         )::Vector{LinearConstraint{N}} where {N<:AbstractFloat}
+    p = ngens(Z)
+    n = dim(Z)
+    if p < n
+        error("can only convert a zonotope of order >= 1")
+    end
+
+    G = Z.generators
+    m = binomial(p, n - 1)
+    constraints = Vector{LinearConstraint{N}}(undef, 2 * m)
+
+    # special handling of 1D case
+    if n == 1
+        if p > 1
+            error("1D-zonotope conversion only supports a single generator")
+        end
+
+        c = Z.center[1]
+        g = G[:, 1][1]
+        constraints[1] = LinearConstraint([N(1)], c + g)
+        constraints[2] = LinearConstraint([N(-1)], g - c)
+        return constraints
+    end
+
+    i = 0
+    c = Z.center
+    for columns in StrictlyIncreasingIndices(p, n-1)
+        i += 1
+        c⁺ = cross_product(view(G, :, columns))
+        normalize!(c⁺, 2)
+
+        Δd = sum(abs.(transpose(G) * c⁺))
+
+        d⁺ = dot(c⁺, c) + Δd
+        c⁻ = -c⁺
+        d⁻ = -d⁺ + 2 * Δd  # identical to dot(c⁻, c) + Δd
+
+        constraints[i] = LinearConstraint(c⁺, d⁺)
+        constraints[i + p] = LinearConstraint(c⁻, d⁻)
+    end
+    @assert i == m "expected 2*$m constraints, but only created 2*$i"
+    return constraints
 end

src/helper_functions.jl

@@ -225,3 +225,172 @@ function reseed(rng::AbstractRNG, seed::Union{Int, Nothing})::AbstractRNG
     end
     return rng
 end
+
+"""
+    cross_product(M::AbstractMatrix{N}) where {N<:Real}
+
+Compute the high-dimensional cross product of ``n-1`` ``n``-dimensional vectors.
+
+### Input
+
+- `M` -- ``n × n - 1``-dimensional matrix
+
+### Output
+
+A vector.
+
+### Algorithm
+
+The cross product is defined as follows:
+
+```math
+\\left[ \\dots, (-1)^{n+1} \\det(M^{[i]}), \\dots \\right]^T
+```
+where ``M^{[i]}`` is defined as ``M`` with the ``i``-th row removed.
+See *Althoff, Stursberg, Buss: Computing Reachable Sets of Hybrid Systems Using
+a Combination of Zonotopes and Polytopes. 2009.*
+"""
+function cross_product(M::AbstractMatrix{N})::Vector{N} where {N<:Real}
+    n = size(M, 1)
+    @assert 1 < n == size(M, 2) + 1 "the matrix must be n x (n-1) dimensional"
+
+    v = Vector{N}(undef, n)
+    minus = false
+    for i in 1:n
+        Mi = view(M, 1:n .!= i, :)  # remove i-th row
+        d = det(Mi)
+        if minus
+            v[i] = -d
+            minus = false
+        else
+            v[i] = d
+            minus = true
+        end
+    end
+    return v
+end
+
+"""
+    StrictlyIncreasingIndices
+
+Iterator over the vectors of `m` strictly increasing indices from 1 to `n`.
+
+### Fields
+
+- `n` -- size of the index domain
+- `m` -- number of indices to choose (resp. length of the vectors)
+
+### Notes
+
+The vectors are modified in-place.
+
+The iterator ranges over ``\\binom{n}{m}`` (`n` choose `m`) possible vectors.
+
+This implementation results in a lexicographic order with the last index growing
+first.
+
+### Examples
+
+```jldoctest
+julia> for v in LazySets.StrictlyIncreasingIndices(4, 2)
+           println(v)
+       end
+[1, 2]
+[1, 3]
+[1, 4]
+[2, 3]
+[2, 4]
+[3, 4]
+```
+"""
+struct StrictlyIncreasingIndices
+    n::Int
+    m::Int
+
+    function StrictlyIncreasingIndices(n::Int, m::Int)
+        @assert n >= m > 0 "require n >= m > 0"
+        new(n, m)
+    end
+end
+
+Base.eltype(::Type{StrictlyIncreasingIndices}) = Vector{Int}
+Base.length(sii::StrictlyIncreasingIndices) = binomial(sii.n, sii.m)
+
+@static if VERSION < v"0.7-"
+@eval begin
+
+# returns [1, 2, ..., m-2, m-1, m-1]
+function Base.start(sii::StrictlyIncreasingIndices)
+    v = [1:sii.m;]
+    if sii.n > sii.m
+        v[end] -= 1
+    end
+    return v
+end
+
+function Base.next(sii::StrictlyIncreasingIndices, state)
+    v = state
+    i = sii.m
+    diff = sii.n
+    if i == diff
+        return (v, nothing)
+    end
+    while v[i] == diff
+        i -= 1
+        diff -= 1
+    end
+    # update vector
+    v[i] += 1
+    for j in i+1:sii.m
+        v[j] = v[j-1] + 1
+    end
+    # detect termination: first index has maximum value
+    if i == 1 && v[1] == (sii.n - sii.m + 1)
+        return (v, nothing)
+    end
+    return (v, v)
+end
+
+Base.done(sii::StrictlyIncreasingIndices, state) = state == nothing
+
+end # @eval
+else
+@eval begin
+
+# initialization
+function Base.iterate(sii::StrictlyIncreasingIndices)
+    v = [1:sii.m;]
+    return (v, v)
+end
+
+# normal iteration
+function Base.iterate(sii::StrictlyIncreasingIndices, state::AbstractVector{Int})
+    v = state
+    i = sii.m
+    diff = sii.n
+    if i == diff
+        return nothing
+    end
+    while v[i] == diff
+        i -= 1
+        diff -= 1
+    end
+    # update vector
+    v[i] += 1
+    for j in i+1:sii.m
+        v[j] = v[j-1] + 1
+    end
+    # detect termination: first index has maximum value
+    if i == 1 && v[1] == (sii.n - sii.m + 1)
+        return (v, nothing)
+    end
+    return (v, v)
+end
+
+# termination
+function Base.iterate(sii::StrictlyIncreasingIndices, state::Nothing)
+    return nothing
+end
+
+end # @eval
+end # if

test/unit_Zonotope.jl

@@ -97,3 +97,19 @@ for N in [Float64, Rational{Int}, Float32]
     @test Z.center == N[2, 3] && diag(Z.generators) == N[4, 5]
     convert(Zonotope, BallInf(N[5, 3], N(2)))
 end
+
+for N in [Float64, Rational{Int}]
+    # conversion to HPolytope
+    # 1D
+    Z = Zonotope(N[0], Matrix{N}(I, 1, 1))
+    P = HPolytope(constraints_list(Z))
+    for d in [N[1], N[-1]]
+        @test ρ(d, P) == ρ(d, Z)
+    end
+    # 2D
+    Z = Zonotope(N[0, 0], Matrix{N}(I, 2, 2))
+    P = HPolytope(constraints_list(Z))
+    for d in LazySets.Approximations.BoxDiagDirections{N}(2)
+        @test ρ(d, P) == ρ(d, Z)
+    end
+end

test/unit_util.jl

@@ -7,4 +7,11 @@ for _dummy_ in 1:1 # avoid global variable warnings
     LazySets.reseed(rng, seed)
     n2 = rand(Int)
     @test n1 == n2
+
+    # StrictlyIncreasingIndices
+    vectors = Vector{AbstractVector{Int}}()
+    for v in LazySets.StrictlyIncreasingIndices(5, 4)
+        push!(vectors, copy(v))
+    end
+    @test vectors == [[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]
 end