From b7e2446123f65f5dde93e3d915ab590c25e2b8ed Mon Sep 17 00:00:00 2001 From: schillic Date: Wed, 9 Jan 2019 19:35:45 +0100 Subject: [PATCH 1/8] sort docs --- docs/src/lib/utils.md | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index 9c23be1553..c1ed403ca2 100644 --- a/docs/src/lib/utils.md +++ b/docs/src/lib/utils.md @@ -7,32 +7,29 @@ end # Utility functions -```@docs -sign_cadlag -ispermutation -@neutral -@absorbing -@declare_array_version -``` - ## Helpers for internal use only ### Functions and Macros ```@docs -@neutral_absorbing -@array_neutral -@array_absorbing -cross_product(::AbstractMatrix{N}) where {N<:Real} -get_radius! an_element_helper -σ_helper binary_search_constraints +cross_product(::AbstractMatrix{N}) where {N<:Real} +get_radius! +ispermutation nonzero_indices samedir +sign_cadlag _random_zero_sum_vector remove_duplicates_sorted! reseed +σ_helper +@neutral +@absorbing +@neutral_absorbing +@declare_array_version +@array_neutral +@array_absorbing ``` ### Types From 6b9e5bb86ef5938bfe58867ded00396bd4d7a4b4 Mon Sep 17 00:00:00 2001 From: schillic Date: Wed, 9 Jan 2019 19:37:48 +0100 Subject: [PATCH 2/8] add isinvertible_sufficient function --- docs/src/lib/utils.md | 1 + src/compat.jl | 2 +- src/helper_functions.jl | 25 +++++++++++++++++++++++++ test/unit_util.jl | 4 ++++ 4 files changed, 31 insertions(+), 1 deletion(-) diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index c1ed403ca2..bdd75979eb 100644 --- a/docs/src/lib/utils.md +++ b/docs/src/lib/utils.md @@ -16,6 +16,7 @@ an_element_helper binary_search_constraints cross_product(::AbstractMatrix{N}) where {N<:Real} get_radius! +isinvertible_sufficient ispermutation nonzero_indices samedir diff --git a/src/compat.jl b/src/compat.jl index 150787b8fc..f1ee86fdd7 100644 --- a/src/compat.jl +++ b/src/compat.jl @@ -8,7 +8,7 @@ using Compat: copyto!, axes, argmax, @warn import Compat.String using Compat.LinearAlgebra import Compat.LinearAlgebra: norm, checksquare, LAPACKException, - SingularException, × + SingularException, ×, cond import Compat.InteractiveUtils.subtypes export _At_mul_B diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 908dc140ac..13f5c9abc6 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -86,6 +86,31 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T return true end +""" + isinvertible_sufficient(M::AbstractMatrix; [cond_tol]::Number=1e6) + +A sufficient check of a matrix being invertible (or nonsingular). + +### Input + +- `M` -- matrix +- `cond_tol` -- (optional, default: `1e6`) tolerance of matrix condition + +### Output + +If the result is `true`, `M` is invertible. +If the result is `false`, this function could not conclude. + +### Algorithm + +We check whether the +[matrix condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices) +`cond(M)` is below some prescribed tolerance. +""" +function isinvertible_sufficient(M::AbstractMatrix; cond_tol::Number=1e6) + return cond(M) < cond_tol +end + """ remove_duplicates_sorted!(v::AbstractVector) diff --git a/test/unit_util.jl b/test/unit_util.jl index 63356d1db8..e9fe12c2f6 100644 --- a/test/unit_util.jl +++ b/test/unit_util.jl @@ -14,4 +14,8 @@ for _dummy_ in 1:1 # avoid global variable warnings 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]] + + # invertible matrix + @test LazySets.isinvertible_sufficient([2 3; 1 2]) + @test !LazySets.isinvertible_sufficient([2 3; 0 0]) end From 09deb3f872e8ed548cf5aa19908edd679ce09779 Mon Sep 17 00:00:00 2001 From: schillic Date: Wed, 9 Jan 2019 19:51:20 +0100 Subject: [PATCH 3/8] faster boundedness check for invertible LinearMap --- src/LinearMap.jl | 6 ++++++ test/unit_LinearMap.jl | 7 ++++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/LinearMap.jl b/src/LinearMap.jl index 468dae5233..ec3b2ec347 100644 --- a/src/LinearMap.jl +++ b/src/LinearMap.jl @@ -216,12 +216,18 @@ Determine whether a linear map is bounded. ### Algorithm We first check if the matrix is zero or the wrapped set is bounded. +If not, we perform a sufficient check whether the matrix is invertible. +If the matrix is invertible, then the map being bounded is equivalent to the +wrapped set being bounded, and hence the map is unbounded. Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref). """ function isbounded(lm::LinearMap)::Bool if iszero(lm.M) || isbounded(lm.X) return true end + if isinvertible_sufficient(lm.M) + return false + end return isbounded_unit_dimensions(lm) end diff --git a/test/unit_LinearMap.jl b/test/unit_LinearMap.jl index 78572e70cf..481cce021d 100644 --- a/test/unit_LinearMap.jl +++ b/test/unit_LinearMap.jl @@ -53,9 +53,10 @@ for N in [Float64, Rational{Int}, Float32] @test lm1_copy.X == lm1.X # boundedness - @test isbounded(lm) - @test isbounded(zeros(N, 2, 2) * HalfSpace(N[1, 1], N(1))) - @test !isbounded(ones(N, 2, 2) * HalfSpace(N[1, 1], N(1))) + @test isbounded(ones(N, 2, 2) * Singleton(N[1, 2])) # bounded set + @test isbounded(zeros(N, 2, 2) * HalfSpace(N[1, 1], N(1))) # zero map + @test !isbounded(N[2 3; 1 2] * HalfSpace(N[1, 1], N(1))) # invertible matrix + @test !isbounded(N[2 3; 0 0] * HalfSpace(N[1, 1], N(1))) # singular matrix (expensive check) # isempty @test !isempty(lm) From ac7d775061f22f07a918b4ce77f2885cb72dc941 Mon Sep 17 00:00:00 2001 From: schillic Date: Wed, 9 Jan 2019 20:50:15 +0100 Subject: [PATCH 4/8] invertible linear_map for H-polyhedra --- docs/src/lib/representations.md | 4 +-- docs/src/man/set_operations.md | 4 +-- src/HPolyhedron.jl | 55 +++++++++++++++++++++++++++++++++ test/unit_Polyhedron.jl | 6 ++++ test/unit_Polytope.jl | 6 ++++ 5 files changed, 70 insertions(+), 5 deletions(-) diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md index a77f594ae4..c8295505ce 100644 --- a/docs/src/lib/representations.md +++ b/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`. diff --git a/docs/src/man/set_operations.md b/docs/src/man/set_operations.md index ff9ca48e2c..498ad3c604 100644 --- a/docs/src/man/set_operations.md +++ b/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 | diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 45c75b4135..44c4c45a3a 100644 --- a/src/HPolyhedron.jl +++ b/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 # ======================================================== diff --git a/test/unit_Polyhedron.jl b/test/unit_Polyhedron.jl index efe5f84e0f..e48e84dbec 100644 --- a/test/unit_Polyhedron.jl +++ b/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 diff --git a/test/unit_Polytope.jl b/test/unit_Polytope.jl index e8845b94cf..1afe9c310c 100644 --- a/test/unit_Polytope.jl +++ b/test/unit_Polytope.jl @@ -106,6 +106,12 @@ for N in [Float64, Rational{Int}, Float32] @test BallInf(N[0, 0], N(1)) ⊆ P @test !(BallInf(N[0, 0], N(1.01)) ⊆ P) + # 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 # ----- From f1712700d6efe01af699e1872b450cd0f92adaac Mon Sep 17 00:00:00 2001 From: schillic Date: Thu, 17 Jan 2019 12:56:19 +0100 Subject: [PATCH 5/8] remove copy method again --- src/HPolyhedron.jl | 26 +++++--------------------- 1 file changed, 5 insertions(+), 21 deletions(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 44c4c45a3a..e4d5b42625 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -2,8 +2,7 @@ using MathProgBase, GLPKMathProgInterface import Base: isempty, rand, - convert, - copy + convert export HPolyhedron, dim, σ, ∈, @@ -438,8 +437,10 @@ A new polyhedron obtained by removing the redundant constraints in `P`. See `remove_redundant_constraints!`. """ function remove_redundant_constraints(P::PT; - backend=GLPKSolverLP()) where {N, PT<:HPoly{N}} - return remove_redundant_constraints!(copy(P), backend=backend) + backend=GLPKSolverLP() + ) where {N, PT<:HPoly{N}} + return remove_redundant_constraints!(PT(copy(constraints_list(P))), + backend=backend) end """ @@ -537,23 +538,6 @@ function linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}} 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 # ======================================================== From d32e77a28a27e7e684b17bf00715d272163fe542 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 19 Jan 2019 09:49:31 +0100 Subject: [PATCH 6/8] rename isinvertible_sufficient --- docs/src/lib/utils.md | 2 +- src/HPolyhedron.jl | 2 +- src/LinearMap.jl | 2 +- src/helper_functions.jl | 4 ++-- test/unit_util.jl | 4 ++-- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index bdd75979eb..cb655153ae 100644 --- a/docs/src/lib/utils.md +++ b/docs/src/lib/utils.md @@ -16,7 +16,7 @@ an_element_helper binary_search_constraints cross_product(::AbstractMatrix{N}) where {N<:Real} get_radius! -isinvertible_sufficient +isinvertible ispermutation nonzero_indices samedir diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index e4d5b42625..0b85e43168 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -522,7 +522,7 @@ 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 !isinvertible(M) if P isa HPolyhedron error("linear maps for polyhedra need to be invertible") end diff --git a/src/LinearMap.jl b/src/LinearMap.jl index ec3b2ec347..29bfe15d1c 100644 --- a/src/LinearMap.jl +++ b/src/LinearMap.jl @@ -225,7 +225,7 @@ function isbounded(lm::LinearMap)::Bool if iszero(lm.M) || isbounded(lm.X) return true end - if isinvertible_sufficient(lm.M) + if isinvertible(lm.M) return false end return isbounded_unit_dimensions(lm) diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 13f5c9abc6..aeb285f4fc 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -87,7 +87,7 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T end """ - isinvertible_sufficient(M::AbstractMatrix; [cond_tol]::Number=1e6) + isinvertible(M::AbstractMatrix; [cond_tol]::Number=1e6) A sufficient check of a matrix being invertible (or nonsingular). @@ -107,7 +107,7 @@ We check whether the [matrix condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices) `cond(M)` is below some prescribed tolerance. """ -function isinvertible_sufficient(M::AbstractMatrix; cond_tol::Number=1e6) +function isinvertible(M::AbstractMatrix; cond_tol::Number=1e6) return cond(M) < cond_tol end diff --git a/test/unit_util.jl b/test/unit_util.jl index e9fe12c2f6..1a06f150d4 100644 --- a/test/unit_util.jl +++ b/test/unit_util.jl @@ -16,6 +16,6 @@ for _dummy_ in 1:1 # avoid global variable warnings @test vectors == [[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]] # invertible matrix - @test LazySets.isinvertible_sufficient([2 3; 1 2]) - @test !LazySets.isinvertible_sufficient([2 3; 0 0]) + @test LazySets.isinvertible([2 3; 1 2]) + @test !LazySets.isinvertible([2 3; 0 0]) end From c6049aeccf09f26861a3116611ba8c692c1d935e Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 19 Jan 2019 09:59:59 +0100 Subject: [PATCH 7/8] add optional argument to isinvertible's callers --- src/HPolyhedron.jl | 16 +++++++++++----- src/LinearMap.jl | 10 ++++++---- src/helper_functions.jl | 10 +++++++--- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 0b85e43168..846e322321 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -502,14 +502,17 @@ function remove_redundant_constraints!(P::PT; end """ - linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}} + linear_map(M::AbstractMatrix{N}, P::PT; [cond_tol=DEFAULT_COND_TOL]::Number) + where {N<:Real, PT<:HPoly{N}} Concrete linear map of a polyhedron in constraint representation. ### Input -- `M` -- matrix -- `P` -- polyhedron in constraint representation +- `M` -- matrix +- `P` -- polyhedron in constraint representation +- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether + the matrix is invertible) ### Output @@ -521,8 +524,11 @@ 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(M) +function linear_map(M::AbstractMatrix{N}, + P::PT; + cond_tol::Number=DEFAULT_COND_TOL + ) where {N<:Real, PT<:HPoly{N}} + if !isinvertible(M; cond_tol=cond_tol) if P isa HPolyhedron error("linear maps for polyhedra need to be invertible") end diff --git a/src/LinearMap.jl b/src/LinearMap.jl index 29bfe15d1c..0b4d156371 100644 --- a/src/LinearMap.jl +++ b/src/LinearMap.jl @@ -201,13 +201,15 @@ function ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real} end """ - isbounded(lm::LinearMap)::Bool + isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool Determine whether a linear map is bounded. ### Input -- `lm` -- linear map +- `lm` -- linear map +- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether + the matrix is invertible) ### Output @@ -221,11 +223,11 @@ If the matrix is invertible, then the map being bounded is equivalent to the wrapped set being bounded, and hence the map is unbounded. Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref). """ -function isbounded(lm::LinearMap)::Bool +function isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool if iszero(lm.M) || isbounded(lm.X) return true end - if isinvertible(lm.M) + if isinvertible(lm.M; cond_tol=cond_tol) return false end return isbounded_unit_dimensions(lm) diff --git a/src/helper_functions.jl b/src/helper_functions.jl index aeb285f4fc..e76066ac30 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -1,3 +1,6 @@ +# default tolerance for matrix condition number (see 'isinvertible') +DEFAULT_COND_TOL = 1e6 + """ sign_cadlag(x::N)::N where {N<:Real} @@ -87,14 +90,15 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T end """ - isinvertible(M::AbstractMatrix; [cond_tol]::Number=1e6) + isinvertible(M::AbstractMatrix; [cond_tol]::Number=DEFAULT_COND_TOL) A sufficient check of a matrix being invertible (or nonsingular). ### Input - `M` -- matrix -- `cond_tol` -- (optional, default: `1e6`) tolerance of matrix condition +- `cond_tol` -- (optional, default: `DEFAULT_COND_TOL`) tolerance of matrix + condition ### Output @@ -107,7 +111,7 @@ We check whether the [matrix condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices) `cond(M)` is below some prescribed tolerance. """ -function isinvertible(M::AbstractMatrix; cond_tol::Number=1e6) +function isinvertible(M::AbstractMatrix; cond_tol::Number=DEFAULT_COND_TOL) return cond(M) < cond_tol end From a1564367b295290e6d48902b65f9197272898ac9 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sat, 19 Jan 2019 13:00:15 +0100 Subject: [PATCH 8/8] Update src/helper_functions.jl Co-Authored-By: schillic --- src/helper_functions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/helper_functions.jl b/src/helper_functions.jl index e76066ac30..eb0ee89e02 100644 --- a/src/helper_functions.jl +++ b/src/helper_functions.jl @@ -1,5 +1,5 @@ # default tolerance for matrix condition number (see 'isinvertible') -DEFAULT_COND_TOL = 1e6 +const DEFAULT_COND_TOL = 1e6 """ sign_cadlag(x::N)::N where {N<:Real}