invertible linear_map for H-polyhedra

docs/src/lib/representations.md

@@ -439,6 +439,7 @@ tosimplehrep(::HPoly{N}) where {N<:Real}
 tohrep(::HPoly{N}) where {N<:Real}
 isempty(::HPoly{N}) where {N<:Real}
 cartesian_product(::HPoly{N}, ::HPoly{N}) where {N<:Real}
+linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}}
 tovrep(::HPoly{N}) where {N<:Real}
 polyhedron(::HPoly{N}) where {N<:Real}
 remove_redundant_constraints
@@ -450,9 +451,6 @@ Inherited from [`LazySet`](@ref):
 * [`radius`](@ref radius(::LazySet, ::Real))
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 
-Inherited from [`AbstractPolytope`](@ref):
-* [`linear_map`](@ref linear_map(::AbstractMatrix{N}, ::AbstractPolytope{N}) where {N<:Real})
-
 #### Polytopes in constraint representation
 
 The following methods are specific for `HPolytope`.

docs/src/man/set_operations.md

@@ -82,8 +82,8 @@ The table entries have the following meaning.
 | `EmptySet`                   | x   | i | x | x          | x | x       | x         |            | x    | x      | x        |
 | `HalfSpace`                  | x   | x | x | x          | x | x       | x         |            |      |        | i        |
 | `HPolygon`/`HPolygonOpt`     | i   | i | x | i          | i | i       | i         | i          |      |        | i        |
-| `HPolyhedron`                | x   | x | x | i          | x | x       | x         |            |      |        | i        |
-| `HPolytope`                  | x   | x | x | i          | x | x       | i         | i          |      |        | i        |
+| `HPolyhedron`                | x   | x | x | i          | x | x       | x         | x          |      |        | i        |
+| `HPolytope`                  | x   | x | x | i          | x | x       | i         | x          |      |        | i        |
 | `Hyperplane`                 | x   | x | x | x          | x | x       | x         |            |      |        | i        |
 | `Hyperrectangle`             | i   | i | i | i          | i | i       | i         | i          | i    | i      | i        |
 | `Interval`                   | x   | i | x | x          | x | i       | i         | i          | i    | i      | i        |

src/HPolyhedron.jl

@@ -16,6 +16,7 @@ export HPolyhedron,
        vertices_list,
        singleton_list,
        isempty,
+       linear_map,
        remove_redundant_constraints,
        remove_redundant_constraints!,
        constrained_dimensions
@@ -499,6 +500,60 @@ function remove_redundant_constraints!(P::PT;
     return P
 end
 
+"""
+    linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}}
+
+Concrete linear map of a polyhedron in constraint representation.
+
+### Input
+
+- `M` -- matrix
+- `P` -- polyhedron in constraint representation
+
+### Output
+
+A polyhedron of the same type as the input (`PT`).
+
+### Algorithm
+
+If the matrix ``M`` is invertible (which we check with a sufficient condition),
+then ``y = M x`` implies ``x = \\text{inv}(M) y`` and we transform the
+constraint system ``A x ≤ b`` to ``A \\text{inv}(M) y ≤ b``.
+"""
+function linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}}
+    if !isinvertible_sufficient(M)
+        if P isa HPolyhedron
+            error("linear maps for polyhedra need to be invertible")
+        end
+        # use the implementation for general polytopes
+        return invoke(linear_map, Tuple{typeof(M), AbstractPolytope{N}}, M, P)
+    end
+    # matrix is invertible
+    invM = inv(M)
+    constraints = Vector{LinearConstraint{N}}(undef, length(constraints_list(P)))
+    for c in constraints_list(P)
+        push!(constraints, LinearConstraint(vec(c.a' * invM), c.b))
+    end
+    return PT(constraints)
+end
+
+"""
+    copy(P::PT) where {N, PT<:HPoly{N}}
+
+Create a copy of a polyhedron.
+
+### Input
+
+- `P` -- polyhedron
+
+### Output
+
+The polyhedron obtained by copying the constraints in `P` using `Base.copy`.
+"""
+function copy(P::PT) where {N, PT<:HPoly{N}}
+    return PT(copy(P.constraints))
+end
+
 # ========================================================
 # External methods that require Polyhedra.jl to be loaded
 # ========================================================

test/unit_Polyhedron.jl

@@ -54,6 +54,9 @@ for N in [Float64, Rational{Int}, Float32]
     @test constrained_dimensions(
         HPolyhedron{N}([LinearConstraint(N[1, 0], N(1))])) == [1]
 
+    # concrete linear map with invertible matrix
+    linear_map(N[2 3; 1 2], p)
+
     if test_suite_polyhedra
         # conversion to and from Polyhedra's VRep data structure
         cl = constraints_list(HPolyhedron(polyhedron(p)))
@@ -69,6 +72,9 @@ for N in [Float64, Rational{Int}, Float32]
         @test !isempty(P)
         addconstraint!(P, LinearConstraint(N[-1, 0], N(-1)))  # x >= 1
         @test isempty(P)
+
+        # concrete linear map with noninvertible matrix throws an error
+        @test_throws ErrorException linear_map(N[2 3; 0 0], P)
     end
 end
 

test/unit_Polytope.jl

@@ -86,11 +86,11 @@ for N in [Float64, Rational{Int}, Float32]
     end
 
     # remove redundant constraints
-    P = HPolytope([HalfSpace([1.0, 0.0], 1.0),
-                   HalfSpace([0.0, 1.0], 1.0),
-                   HalfSpace([-1.0, -0.0], 1.0),
-                   HalfSpace([-0.0, -1.0], 1.0),
-                   HalfSpace([2.0, 0.0], 2.0)]) # redundant
+    P = HPolytope([HalfSpace(N[1, 0], N(1)),
+                   HalfSpace(N[0, 1], N(1)),
+                   HalfSpace(N[-1, -0], N(1)),
+                   HalfSpace(N[-0, -1], N(1)),
+                   HalfSpace(N[2, 0], N(2))]) # redundant
 
     Pred = remove_redundant_constraints(P)
     @test length(Pred.constraints) == 4
@@ -100,6 +100,12 @@ for N in [Float64, Rational{Int}, Float32]
     remove_redundant_constraints!(P)
     @test length(P.constraints) == 4
 
+    # concrete linear map
+    linear_map(N[2 3; 1 2], P)  # invertible matrix
+    if test_suite_polyhedra
+        linear_map(N[2 3; 0 0], P)  # noninvertible matrix
+    end
+
     # -----
     # V-rep
     # -----