Merge pull request #3896 from JuliaReach/schillic/linear_map_universe

Fix and improve `linear_map` for `Universe` with invertible map

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))
@@ -116,7 +117,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))
 * [`isapprox`](@ref isapprox(::LazySet, ::LazySet))
 * [`isdisjoint`](@ref isdisjoint(::LazySet, ::LazySet))

src/Sets/Universe/UniverseModule.jl

@@ -3,9 +3,10 @@ module UniverseModule
 using Reexport, Requires
 
 using ..LazySets: LazySet, AbstractPolyhedron, default_polyhedra_backend,
-                  _witness_result_empty, @validate, @validate_commutative
+                  _linear_map_polyhedron, _witness_result_empty, @validate,
+                  @validate_commutative
 using Random: AbstractRNG, GLOBAL_RNG
-using ReachabilityBase.Arrays: SingleEntryVector
+using ReachabilityBase.Arrays: SingleEntryVector, isinvertible
 using ReachabilityBase.Distribution: reseed!
 using ReachabilityBase.Iteration: EmptyIterator
 using ReachabilityBase.Require: require
@@ -14,9 +15,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 +50,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,18 @@
+"""
+# Extended help
+
+    linear_map(M::AbstractMatrix, U::Universe)
+
+### Output
+
+Either a `Universe` if the map is invertible, or the result of the default implementation, which
+typically yields an `HPolyhedron`.
+"""
+@validate function linear_map(M::AbstractMatrix, U::Universe)
+    if isinvertible(M)
+        N = eltype(U)
+        m = size(M, 1)
+        return Universe{N}(m)
+    end
+    return _linear_map_polyhedron(M, U)
+end

test/Sets/Universe.jl

@@ -209,6 +209,8 @@ for N in @tN([Float64, Float32, Rational{Int}])
     # affine_map (part 1)
     @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])
+    U2 = affine_map(N[1 0; 0 1], U, N[1, 1])
+    @test isidentical(U2, U)
 
     # distance (between point and set)
     @test_throws DimensionMismatch distance(U, N[0])
@@ -241,6 +243,8 @@ for N in @tN([Float64, Float32, Rational{Int}])
 
     # linear_map (part 1)
     @test_throws DimensionMismatch linear_map(ones(N, 2, 1), U)
+    U2 = linear_map(N[1 0; 0 1], U)
+    @test isidentical(U2, U)
 
     # linear_map_inverse
     U2 = LazySets.linear_map_inverse(ones(N, 2, 3), U)
@@ -459,34 +463,25 @@ for N in @tN([Float64, Float32])
     @test_throws MethodError rationalize(U2)
 
     # affine_map (part 2)
-    @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        @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 LazySet{N}
-            @test isequivalent(X, Hyperplane(N[0, 0, 1], N(3)))
-        end
+    @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
 
     # exponential_map
     U2 = exponential_map(ones(N, 2, 2), U)
-    @test_broken isidentical(U2, U)  # TODO this should change
+    @test isidentical(U2, U)
 
     # linear_map (part 2)
-    @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
-        U2 = linear_map(N[1 0; 0 1], U)
-        @test U2 isa HPolyhedron{N}  # TODO this should change (HPolyhedron has no dimension)
-        @test_broken isequivalent(U2, U)
-        @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}
-            @test isequivalent(X, Hyperplane(N[0, 0, 1], N(0)))
-        end
+    @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
 end
 
@@ -496,32 +491,34 @@ for N in [Float64]
     # affine_map (part 3)
     @static if isdefined(@__MODULE__, :Polyhedra) && isdefined(@__MODULE__, :CDDLib)
         @static if VERSION < v"1.12"
-            # TODO this should work, even without Polyhedra
+            # 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}
-            @test isequivalent(X, Line2D(N[1, -1], N(2)))
+            @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}
-            @test isequivalent(X, Singleton(N[2, 0]))
+            @test X isa HPolyhedron{N} && isequivalent(X, Singleton(N[2, 0]))
         end
     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}
-            @test isequivalent(X, Line2D(N[1, -1], N(0)))
+            @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}
-            @test isequivalent(X, ZeroSet{N}(2))
+            @test X isa HPolyhedron{N} && isequivalent(X, ZeroSet{N}(2))
         end
     end
 end