[PolyhedralGeometry] Check containment of points in hyperplanes and h…

…alfspaces (#4392) * [PolyhedralGeometry] Check containment of points in hyperplanes and halfspaces * [PolyhedralGeometry] Streamline docs a bit * docs: fix signature for facets ref --------- Co-authored-by: Benjamin Lorenz <[email protected]>
oscar-system · Jan 3, 2025 · e5a71df · e5a71df
1 parent cf2e23f
commit e5a71df
Show file tree

Hide file tree

Showing 9 changed files with 44 additions and 24 deletions.
diff --git a/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md b/docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md
@@ -7,7 +7,7 @@ CurrentModule = Oscar
 ## Geometric data
 
 ```@docs
-facets(P::Polyhedron)
+facets(as::Type{T}, P::Polyhedron{S}) where {S<:scalar_types,T<:Union{AffineHalfspace{S},AffineHyperplane{S},Pair{R,S} where R,Polyhedron{S}}}
 vertices(as::Type{PointVector{T}}, P::Polyhedron{T}) where {T<:scalar_types}
 rays(as::Type{RayVector{T}}, P::Polyhedron{T}) where {T<:scalar_types}
 rays_modulo_lineality(P::Polyhedron{T}) where T<:scalar_types

diff --git a/docs/src/PolyhedralGeometry/Polyhedra/constructions.md b/docs/src/PolyhedralGeometry/Polyhedra/constructions.md
@@ -19,7 +19,7 @@ polyhedron(::Oscar.scalar_type_or_field, A::AnyVecOrMat, b::AbstractVector)
 polyhedron(::Oscar.scalar_type_or_field, I::Union{Nothing, AbstractCollection[AffineHalfspace]}, E::Union{Nothing, AbstractCollection[AffineHyperplane]} = nothing)
 ```
 
-The complete $H$-representation can be retrieved using [`facets`](@ref facets(as::Type{T}, P::Polyhedron{S}) where {S<:scalar_types,T<:Union{AffineHalfspace{S},Pair{R,S} where R,Polyhedron{S}}})
+The complete $H$-representation can be retrieved using [`facets`](@ref facets(as::Type{T}, P::Polyhedron{S}) where {S<:scalar_types,T<:Union{AffineHalfspace{S},AffineHyperplane{S},Pair{R,S} where R,Polyhedron{S}}})
 and [`affine_hull`](@ref affine_hull(P::Polyhedron{T}) where {T<:scalar_types}):
 ```jldoctest
 julia> P = polyhedron(([-1 0; 1 0], [0,1]), ([0 1], [0]))

diff --git a/docs/src/PolyhedralGeometry/cones.md b/docs/src/PolyhedralGeometry/cones.md
@@ -37,6 +37,7 @@ ambient_dim(C::Cone)
 Base.in(v::AbstractVector, C::Cone)
 Base.issubset(C0::Cone{T}, C1::Cone{T}) where T<:scalar_types
 facet_degrees(C::Cone)
+facets(as::Type{<:Union{LinearHalfspace{T},LinearHyperplane{T},Cone{T}}}, C::Cone{T}) where {T<:scalar_types}
 f_vector(C::Cone)
 hilbert_basis(C::Cone{QQFieldElem})
 codim(C::Cone)

diff --git a/src/PolyhedralGeometry/Cone/properties.jl b/src/PolyhedralGeometry/Cone/properties.jl
@@ -535,6 +535,7 @@ Return the facets of `C` in the format defined by `as`.
 
 The allowed values for `as` are
 * `Halfspace` (or its subtype `LinearHalfspace`),
+* `Hyperplane` (or its subtype `LinearHyperplane1),
 * `Cone`.
 
 # Examples
@@ -549,13 +550,19 @@ julia> f = facets(Halfspace, c)
 -x_2 + x_3 <= 0
 ```
 """
-facets(as::Type{<:Union{LinearHalfspace{T},Cone{T}}}, C::Cone{T}) where {T<:scalar_types} =
+facets(
+  as::Type{<:Union{LinearHalfspace{T},LinearHyperplane{T},Cone{T}}}, C::Cone{T}
+) where {T<:scalar_types} =
   SubObjectIterator{as}(C, _facet_cone, n_facets(C))
 
 _facet_cone(
   U::Type{LinearHalfspace{T}}, C::Cone{T}, i::Base.Integer
 ) where {T<:scalar_types} =
   linear_halfspace(coefficient_field(C), -pm_object(C).FACETS[[i], :])::U
+_facet_cone(
+  U::Type{LinearHyperplane{T}}, C::Cone{T}, i::Base.Integer
+) where {T<:scalar_types} =
+  linear_hyperplane(coefficient_field(C), -pm_object(C).FACETS[[i], :])::U
 
 _facet_cone(::Type{Cone{T}}, C::Cone{T}, i::Base.Integer) where {T<:scalar_types} =
   Cone{T}(Polymake.polytope.facet(pm_object(C), i - 1), coefficient_field(C))
@@ -572,6 +579,8 @@ facets(C::Cone{T}) where {T<:scalar_types} = facets(LinearHalfspace{T}, C)
 
 facets(::Type{<:Halfspace}, C::Cone{T}) where {T<:scalar_types} =
   facets(LinearHalfspace{T}, C)
+facets(::Type{<:Hyperplane}, C::Cone{T}) where {T<:scalar_types} =
+  facets(LinearHyperplane{T}, C)
 
 facets(::Type{Cone}, C::Cone{T}) where {T<:scalar_types} = facets(Cone{T}, C)
 

diff --git a/src/PolyhedralGeometry/Polyhedron/properties.jl b/src/PolyhedralGeometry/Polyhedron/properties.jl
@@ -492,6 +492,7 @@ Return the facets of `P` in the format defined by `as`.
 
 The allowed values for `as` are
 * `Halfspace` (or its subtype `AffineHalfspace`),
+* `Hyperplane` (or its subtype `AffineHyperplane`),
 * `Polyhedron`,
 * `Pair`.
 
@@ -521,7 +522,10 @@ x_3 <= 1
 """
 facets(
   as::Type{T}, P::Polyhedron{S}
-) where {S<:scalar_types,T<:Union{AffineHalfspace{S},Pair{R,S} where R,Polyhedron{S}}} =
+) where {
+  S<:scalar_types,
+  T<:Union{AffineHalfspace{S},AffineHyperplane{S},Pair{R,S} where R,Polyhedron{S}},
+} =
   SubObjectIterator{as}(P, _facet_polyhedron, n_facets(P))
 
 function _facet_polyhedron(
@@ -545,6 +549,12 @@ function _facet_polyhedron(
     coefficient_field(P),
   )
 end
+function _facet_polyhedron(
+  U::Type{AffineHyperplane{S}}, P::Polyhedron{S}, i::Base.Integer
+) where {S<:scalar_types}
+  h = decompose_hdata(view(pm_object(P).FACETS, [_facet_index(pm_object(P), i)], :))
+  return affine_hyperplane(coefficient_field(P), h[1], h[2][])::U
+end
 
 _affine_inequality_matrix(::Val{_facet_polyhedron}, P::Polyhedron) =
   -_remove_facet_at_infinity(pm_object(P))
@@ -582,30 +592,12 @@ facets(::Type{<:Pair}, P::Polyhedron{T}) where {T<:scalar_types} =
 facets(::Type{Polyhedron}, P::Polyhedron{T}) where {T<:scalar_types} =
   facets(Polyhedron{T}, P)
 
-@doc raw"""
-    facets(P::Polyhedron)
-
-Return the facets of `P` as halfspaces.
-
-# Examples
-We can retrieve the six facets of the 3-dimensional cube this way:
-```jldoctest
-julia> C = cube(3);
-
-julia> facets(C)
-6-element SubObjectIterator{AffineHalfspace{QQFieldElem}} over the halfspaces of R^3 described by:
--x_1 <= 1
-x_1 <= 1
--x_2 <= 1
-x_2 <= 1
--x_3 <= 1
-x_3 <= 1
-```
-"""
 facets(P::Polyhedron{T}) where {T<:scalar_types} = facets(AffineHalfspace{T}, P)
 
 facets(::Type{<:Halfspace}, P::Polyhedron{T}) where {T<:scalar_types} =
   facets(AffineHalfspace{T}, P)
+facets(::Type{<:Hyperplane}, P::Polyhedron{T}) where {T<:scalar_types} =
+  facets(AffineHyperplane{T}, P)
 
 function _facet_index(P::Polymake.BigObject, i::Base.Integer)
   i < _facet_at_infinity(P) && return i

diff --git a/src/PolyhedralGeometry/iterators.jl b/src/PolyhedralGeometry/iterators.jl
@@ -222,6 +222,15 @@ function Base.:(==)(x::Hyperplane, y::Hyperplane)
   return (r .* ax == ay) && (r * negbias(x) == negbias(y))
 end
 
+Base.in(x::AbstractVector, y::Hyperplane) = (dot(x, normal_vector(y)) == negbias(y))
+Base.in(x::AbstractVector, y::Halfspace) = (dot(x, normal_vector(y)) <= negbias(y))
+# A ray vector needs a base point for containment in an affine space, so we
+# just error when this combination is tested.
+Base.in(x::RayVector, y::T) where {T<:Union{AffineHalfspace,
+    AffineHyperplane}} =
+  throw(ArgumentError("Containment of RayVector in affine spaces is not
+                      well-defined."))
+
 Base.hash(x::T, h::UInt) where {T<:Union{AffineHalfspace,AffineHyperplane}} =
   hash((x.a, x.b), hash(T, h))
 

diff --git a/test/PolyhedralGeometry/cone.jl b/test/PolyhedralGeometry/cone.jl
@@ -116,6 +116,8 @@
 
     @test n_facets(Cone5) == 4
     @test relative_interior_point(Cone1) == f.([1//2, 1//2])
+    @test length(findall(f->[1,0,0] in f, facets(Hyperplane, Cone5))) == 2
+    @test length(findall(f->[1,0,0] in f, facets(Halfspace, Cone5))) == 4
   end
 
   @testset "constructors" begin

diff --git a/test/PolyhedralGeometry/iterators.jl b/test/PolyhedralGeometry/iterators.jl
@@ -41,4 +41,9 @@
     soi2 = SubObjectIterator{PointVector{QQFieldElem}}(c, _testSOI, 4)
     @test_throws ArgumentError point_matrix(soi2)
   end
+  @testset "RayVector" begin
+    rv = ray_vector([1,0,0])
+    F = facets(cube(3))
+    @test_throws ArgumentError rv in F[1]
+  end
 end
diff --git a/test/PolyhedralGeometry/polyhedron.jl b/test/PolyhedralGeometry/polyhedron.jl
@@ -205,6 +205,8 @@
 
     @test dim(full) == ambient_dim(full)
     @test lineality_dim(full) == 3
+    @test length(findall(f-> [1,0] in f, facets(Hyperplane, Q0))) == 2
+    @test length(findall(f-> [1,0] in f, facets(Halfspace, Q0))) == 3
   end
 
   @testset "volume" begin