Merge pull request #3953 from JuliaReach/schillic/hrep

Fix `linear_map` elimination algorithm with empty constraints

Author: schillic

src/Interfaces/AbstractPolyhedron.jl

@@ -73,8 +73,14 @@ end
 # generic function for the AbstractPolyhedron interface => returns an HPolyhedron
 function _linear_map_hrep_helper(M::AbstractMatrix, P::LazySet,
                                  algo::AbstractLinearMapAlgorithm)
-    constraints = _linear_map_hrep(M, P, algo)
-    return HPolyhedron(constraints)
+    clist = _linear_map_hrep(M, P, algo)
+    if isempty(clist)
+        N = promote_type(eltype(M), eltype(P))
+        n = size(M, 1)
+        return Universe{N}(n)
+    else
+        return HPolyhedron(clist)
+    end
 end
 
 # internal functions; defined here due to optional dependencies and submodules

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -508,7 +508,12 @@ function _linear_map_hrep(M::AbstractMatrix, P::AbstractPolyhedron, algo::Linear
 
     # append zeros to the existing constraints, in the last m-n coordinates
     # TODO: cast to common vector type instead of Vector(c.a), see #1942, #1952
-    cext = [HalfSpace(vcat(Vector(c.a), zeros(N, m - n)), c.b) for c in constraints_list(P)]
+    clist = constraints_list(P)
+    if isempty(clist)
+        cext = HalfSpace{N,Vector{N}}[]
+    else
+        cext = [HalfSpace(vcat(Vector(c.a), zeros(N, m - n)), c.b) for c in clist]
+    end
 
     # now fix the last m-n coordinates to zero
     id_out = Matrix(one(N) * I, m - n, m - n)
@@ -529,16 +534,20 @@ end
 function _linear_map_hrep(M::AbstractMatrix, P::AbstractPolyhedron, algo::LinearMapElimination)
     m, n = size(M)
     N = promote_type(eltype(M), eltype(P))
-    ₋Id_m = Matrix(-one(N) * I, m, m)
+    Id_neg = Matrix(-one(N) * I, m, m)
     backend = algo.backend
     method = algo.method
 
     # extend the polytope storing the y variables first
     # append zeros to the existing constraints, in the last m-n coordinates
     # TODO: cast to common vector type instead of hard-coding Vector(c.a), see #1942 and #1952
-    Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(N, m), Vector(c.a)), c.b)
-                for c in constraints_list(P)]
-    y_eq_Mx = [Polyhedra.HyperPlane(vcat(₋Id_m[i, :], Vector(M[i, :])), zero(N)) for i in 1:m]
+    clist = constraints_list(P)
+    if isempty(clist)
+        Ax_leq_b = Polyhedra.HalfSpace{N,Vector{N}}[]
+    else
+        Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(N, m), Vector(c.a)), c.b) for c in clist]
+    end
+    y_eq_Mx = [Polyhedra.HyperPlane(vcat(Id_neg[i, :], Vector(M[i, :])), zero(N)) for i in 1:m]
 
     Phrep = Polyhedra.hrep(y_eq_Mx, Ax_leq_b)
     Phrep = polyhedron(Phrep, backend) # define concrete subtype

src/Interfaces/LazySet_functions.jl

@@ -572,7 +572,6 @@ end
 The default implementation applies the functions `exp` and `linear_map`.
 """
 @validate function exponential_map(M::AbstractMatrix, X::LazySet)
-    n = dim(X)
     return linear_map(exp(M), X)
 end
 

test/Sets/Universe.jl

@@ -463,62 +463,38 @@ for N in @tN([Float64, Float32])
     @test_throws MethodError rationalize(U2)
 
     # affine_map (part 2)
-    @static if VERSION < v"1.12"
-        # TODO this should work with older versions, see below
-        @test_broken affine_map(N[1 0; 0 1; 0 0], U, N[1, 1, 3])
-    else
-        X = affine_map(N[1 0; 0 1; 0 0], U, N[1, 1, 3])
-        @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 0, 1], N(3)))
-    end
+    X = affine_map(N[1 0; 0 1; 0 0], U, N[1, 1, 3])
+    @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 0, 1], N(3)))
 
     # exponential_map
     U2 = exponential_map(ones(N, 2, 2), U)
     @test isidentical(U2, U)
 
     # linear_map (part 2)
-    @static if VERSION < v"1.12"
-        # TODO this should work with older versions, see below
-        @test_broken linear_map(N[1 0; 0 1; 0 0], U)
-    else
-        X = linear_map(N[1 0; 0 1; 0 0], U)
-        @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 0, 1], N(0)))
-    end
+    X = linear_map(N[1 0; 0 1; 0 0], U)
+    @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 0, 1], N(0)))
 end
 
 for N in [Float64]
     U = Universe{N}(2)
 
     # affine_map (part 3)
     @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        @static if VERSION < v"1.12"
-            # TODO these should work with older versions, see below
-            @test_broken affine_map(N[1 2; 0 0], U, N[1, 1])
-            @test_broken affine_map(ones(N, 2, 2), U, N[2, 0])
-            @test_broken affine_map(zeros(N, 2, 2), U, N[2, 0])
-        else
-            X = affine_map(N[1 2; 0 0], U, N[1, 1])  # projection to axis
-            @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 1], N(1)))
-            X = affine_map(ones(N, 2, 2), U, N[2, 0])  # projection to line
-            @test X isa HPolyhedron{N} && isequivalent(X, Line2D(N[1, -1], N(2)))
-            X = affine_map(zeros(N, 2, 2), U, N[2, 0])  # zero map
-            @test X isa HPolyhedron{N} && isequivalent(X, Singleton(N[2, 0]))
-        end
+        X = affine_map(N[1 2; 0 0], U, N[1, 1])  # projection to axis
+        @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 1], N(1)))
+        X = affine_map(ones(N, 2, 2), U, N[2, 0])  # projection to line
+        @test X isa HPolyhedron{N} && isequivalent(X, Line2D(N[1, -1], N(2)))
+        X = affine_map(zeros(N, 2, 2), U, N[2, 0])  # zero map
+        @test X isa HPolyhedron{N} && isequivalent(X, Singleton(N[2, 0]))
     end
 
     # linear_map (part 3)
     @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        @static if VERSION < v"1.12"
-            # TODO these should work with older versions, see below
-            @test_broken linear_map(N[1 2; 0 0], U)
-            @test_broken linear_map(ones(N, 2, 2), U)
-            @test_broken linear_map(zeros(N, 2, 2), U)
-        else
-            X = linear_map(N[1 2; 0 0], U)  # projection to axis
-            @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 1], N(0)))
-            X = linear_map(ones(N, 2, 2), U)  # projection to line
-            @test X isa HPolyhedron{N} && isequivalent(X, Line2D(N[1, -1], N(0)))
-            X = linear_map(zeros(N, 2, 2), U)  # zero map
-            @test X isa HPolyhedron{N} && isequivalent(X, ZeroSet{N}(2))
-        end
+        X = linear_map(N[1 2; 0 0], U)  # projection to axis
+        @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 1], N(0)))
+        X = linear_map(ones(N, 2, 2), U)  # projection to line
+        @test X isa HPolyhedron{N} && isequivalent(X, Line2D(N[1, -1], N(0)))
+        X = linear_map(zeros(N, 2, 2), U)  # zero map
+        @test X isa HPolyhedron{N} && isequivalent(X, ZeroSet{N}(2))
     end
 end