fix 'linear_map' for Universe

Author: schillic

docs/src/lib/sets/Universe.md

@@ -84,6 +84,7 @@ Undocumented implementations:
 * [`reflect`](@ref reflect(::LazySet))
 * [`distance`](@ref distance(::AbstractVector, ::LazySet))
 * [`∈`](@ref ∈(::AbstractVector, ::LazySet))
+* [`linear_map`](@ref linear_map(::AbstractMatrix, ::LazySet))
 * [`permute`](@ref permute(::LazySet, ::AbstractVector{Int}))
 * [`project`](@ref project(::LazySet, ::AbstractVector{Int}))
 * [`rationalize`](@ref rationalize(::LazySet))
@@ -115,7 +116,6 @@ Inherited from [`LazySet`](@ref):
 * [`affine_map`](@ref affine_map(::AbstractMatrix, ::LazySet, ::AbstractVector))
 * [`exponential_map`](@ref exponential_map(::AbstractMatrix, ::LazySet))
 * [`is_interior_point`](@ref is_interior_point(::AbstractVector, ::LazySet))
-* [`linear_map`](@ref linear_map(::AbstractMatrix, ::LazySet))
 * [`sample`](@ref sample(::LazySet, ::Int=1))
 * [`≈`](@ref ≈(::LazySet, ::LazySet))
 * [`isdisjoint`](@ref isdisjoint(::LazySet, ::LazySet))

src/Sets/Universe/UniverseModule.jl

@@ -14,9 +14,9 @@ using ReachabilityBase.Require: require
                         constraints_list, diameter, dim, isbounded,
                         isboundedtype, isempty, isoperationtype, isuniversal,
                         norm, radius, rand, rectify, reflect, volume, ∈,
-                        permute, project, scale, scale!, ρ, σ, translate,
-                        translate!, cartesian_product, convex_hull, difference,
-                        distance, intersection, isdisjoint, ⊆,
+                        linear_map, permute, project, scale, scale!, ρ, σ,
+                        translate, translate!, cartesian_product, convex_hull,
+                        difference, distance, intersection, isdisjoint, ⊆,
                         linear_combination, minkowski_difference, minkowski_sum
 @reexport import ..LazySets: constrained_dimensions, linear_map_inverse,
                              rationalize, tosimplehrep, triangulate
@@ -49,6 +49,7 @@ include("rectify.jl")
 include("reflect.jl")
 include("volume.jl")
 include("in.jl")
+include("linear_map.jl")
 include("permute.jl")
 include("project.jl")
 include("scale.jl")

src/Sets/Universe/linear_map.jl

@@ -0,0 +1,41 @@
+"""
+# Extended help
+
+    linear_map(M::AbstractMatrix, U::Universe)
+
+### Output
+
+Either a `Universe`, or an `HPolyhedron` if the map contains zero rows.
+"""
+function linear_map(M::AbstractMatrix, U::Universe)
+    m, n = size(M)
+    @assert n == dim(U) "a linear map of size $(size(M)) cannot be " *
+                        "applied to a set of dimension $(dim(U))"
+
+    N = promote_type(eltype(M), eltype(U))
+    first_zero_row = 0
+    for i in 1:m
+        if iszero(M[i, :])
+            first_zero_row = i
+            break
+        end
+    end
+    if first_zero_row == 0
+        return Universe{N}(m)
+    end
+
+    require(@__MODULE__, :LazySets; fun_name="linear_map")
+    clist = Vector{HalfSpace{N,SingleEntryVector{N}}}()
+    _push_zero_constraint!(clist, N, first_zero_row, m)
+    for i in (first_zero_row + 1):m
+        if iszero(M[i, :])
+            _push_zero_constraint!(clist, N, i, m)
+        end
+    end
+    return HPolyhedron(clist)
+end
+
+function _push_zero_constraint!(clist, N, i, n)
+    push!(clist, HalfSpace(SingleEntryVector(i, n, one(N)), zero(N)))
+    push!(clist, HalfSpace(SingleEntryVector(i, n, N(-1)), zero(N)))
+end

test/Sets/Universe.jl

@@ -191,14 +191,12 @@ for N in [Float64, Float32, Rational{Int}]
     # affine_map
     @test_throws DimensionMismatch affine_map(ones(N, 2, 3), U, N[1, 1])
     @test_throws DimensionMismatch affine_map(ones(N, 2, 2), U, N[1])
-    @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        # TODO this should work, even without Polyhedra
-        @test_broken affine_map(ones(N, 2, 2), U, N[1, 1])
-        # U2 = affine_map(ones(N, 2, 2), U, N[1, 1])
-        # @test isidentical(U, U2)
-        # U2 = affine_map(ones(N, 3, 2), U, N[1, 1, 3])
-        # @test isidentical(U3, U2)
-    end
+    U2 = affine_map(ones(N, 2, 2), U, N[1, 1])
+    @test isidentical(U, U2)
+    X = affine_map(N[1 2; 0 0], U, N[1, 1])
+    @test X isa HPolyhedron{N} && isequivalent(X, Hyperplane(N[0, 1], N(1)))
+    U2 = affine_map(ones(N, 3, 2), U, N[1, 1, 3])
+    @test isidentical(U3, U2)
 
     # distance (between point and set)
     @test_throws DimensionMismatch distance(U, N[0])
@@ -228,16 +226,14 @@ for N in [Float64, Float32, Rational{Int}]
 
     # linear_map
     @test_throws DimensionMismatch linear_map(ones(N, 2, 3), U)
-    @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        # TODO this should work, even without Polyhedra
-        @test_broken linear_map(ones(N, 2, 2), U)
-        # U2 = linear_map(ones(N, 2, 2), U)
-        # @test_broken isidentical(U, U2)
-        @test_broken linear_map(ones(N, 3, 2), U)
-        # U2 = linear_map(ones(N, 3, 2), U)
-        # @test U2 isa HPolyhedron{N}  # TODO this should change
-        # @test_broken isidentical(U3, U2)
-    end
+    U2 = linear_map(ones(N, 2, 2), U)
+    @test isidentical(U, U2)
+    U2 = linear_map(ones(N, 3, 2), U)
+    @test isidentical(U3, U2)
+    P = linear_map(N[1 0; 0 0], U)
+    @test P isa HPolyhedron{N} && isequivalent(P, Hyperplane(N[0, 1], N(0)))
+    P = linear_map(N[0 0; 0 0; 0 0], U)
+    @test P isa HPolyhedron{N} && isequivalent(P, ZeroSet{N}(3))
 
     # linear_map_inverse
     U2 = LazySets.linear_map_inverse(ones(N, 2, 3), U)
@@ -461,5 +457,5 @@ for N in [Float64, Float32]
 
     # exponential_map
     U2 = exponential_map(ones(N, 2, 2), U)
-    @test_broken isidentical(U, U2)  # TODO this should change
+    @test isidentical(U, U2)
 end