Polyhedral: rework input type handling, add support for QQBarField ob…

…jects (#3308)

* require Polymake 0.11.12
* fix augment for some type issues
* helpers for QQBarField
* rework type detection for polyhedral objects, dont fallback to QQ for proper fields
* add type detection to constructors
* adjust visualize for more types
* tests, adjust autodetected type, add qqbar example

Project.toml

@@ -37,7 +37,7 @@ JSON3 = "1.13.2"
 LazyArtifacts = "1.6"
 Nemo = "0.41.3"
 Pkg = "1.6"
-Polymake = "0.11.8"
+Polymake = "0.11.12"
 Preferences = "1"
 Random = "1.6"
 RandomExtensions = "0.4.3"

docs/src/PolyhedralGeometry/Polyhedra/constructions.md

@@ -61,7 +61,7 @@ This is a standard triangle, defined via a (redundant) $V$-representation  and
 its unique minimal $H$-representation:
 
 ```jldoctest
-julia> T = convex_hull([ 0 0 ; 1 0 ; 0 1; 0 1/2 ])
+julia> T = convex_hull([ 0 0 ; 1 0 ; 0 1; 0 1//2 ])
 Polyhedron in ambient dimension 2
 
 julia> halfspace_matrix_pair(facets(T))

src/PolyhedralGeometry/Cone/constructors.jl

@@ -77,8 +77,8 @@ function positive_hull(f::scalar_type_or_field, R::AbstractCollection[RayVector]
 end
 # Redirect everything to the above constructor, use QQFieldElem as default for the
 # scalar type T.
-positive_hull(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(QQFieldElem, R, L; non_redundant=non_redundant)
-cone(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(QQFieldElem, R, L; non_redundant=non_redundant)
+positive_hull(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant)
+cone(R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(_guess_fieldelem_type(R, L), R, L; non_redundant=non_redundant)
 cone(f::scalar_type_or_field, R::AbstractCollection[RayVector], L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = positive_hull(f, R, L; non_redundant=non_redundant)
 cone(f::scalar_type_or_field, x...) = positive_hull(f, x...)
 
@@ -166,7 +166,7 @@ end
 
 cone_from_inequalities(x...) = cone_from_inequalities(QQFieldElem, x...)
 
-cone_from_equations(E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) = cone_from_equations(QQFieldElem, E; non_redundant = non_redundant)
+cone_from_equations(E::AbstractCollection[LinearHyperplane]; non_redundant::Bool = false) = cone_from_equations(_guess_fieldelem_type(E), E; non_redundant = non_redundant)
 
 """
     pm_object(C::Cone)

src/PolyhedralGeometry/PolyhedralComplex/constructors.jl

@@ -106,13 +106,13 @@ function polyhedral_complex(f::scalar_type_or_field,
     ), parent_field)
   end
 end
-# default scalar type: `QQFieldElem`
+# default scalar type: guess from input, fallback QQ
 polyhedral_complex(polyhedra::IncidenceMatrix, 
                 vr::AbstractCollection[PointVector], 
                 far_vertices::Union{Vector{Int}, Nothing} = nothing, 
                 L::Union{AbstractCollection[RayVector], Nothing} = nothing;
                 non_redundant::Bool = false) =
-  polyhedral_complex(QQFieldElem, polyhedra, vr, far_vertices, L; non_redundant=non_redundant)
+  polyhedral_complex(_guess_fieldelem_type(vr, L), polyhedra, vr, far_vertices, L; non_redundant=non_redundant)
 
 function polyhedral_complex(f::scalar_type_or_field, v::AbstractCollection[PointVector], vi::IncidenceMatrix, r::AbstractCollection[RayVector], ri::IncidenceMatrix, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false)
   vr = [unhomogenized_matrix(v); unhomogenized_matrix(r)]

src/PolyhedralGeometry/PolyhedralFan/constructors.jl

@@ -91,8 +91,8 @@ function polyhedral_fan(f::scalar_type_or_field,
   end
 end
 polyhedral_fan(f::scalar_type_or_field, Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool = false) = polyhedral_fan(f, Incidence, Rays, nothing; non_redundant)
-polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector], LS::Union{AbstractCollection[RayVector], Nothing}; non_redundant::Bool = false) = polyhedral_fan(QQFieldElem, Incidence, Rays, LS; non_redundant)
-polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool = false) = polyhedral_fan(QQFieldElem, Incidence, Rays; non_redundant)
+polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector], LS::Union{AbstractCollection[RayVector], Nothing}; non_redundant::Bool = false) = polyhedral_fan(_guess_fieldelem_type(Rays, LS), Incidence, Rays, LS; non_redundant)
+polyhedral_fan(Incidence::IncidenceMatrix, Rays::AbstractCollection[RayVector]; non_redundant::Bool = false) = polyhedral_fan(_guess_fieldelem_type(Rays), Incidence, Rays; non_redundant)
 
 """
     pm_object(PF::PolyhedralFan)

src/PolyhedralGeometry/Polyhedron/constructors.jl

@@ -20,9 +20,9 @@ struct Polyhedron{T<:scalar_types} <: PolyhedralObject{T} #a real polymake polyh
     Polyhedron{QQFieldElem}(p::Polymake.BigObject) = new{QQFieldElem}(p, QQ)
 end
 
-# default scalar type: `QQFieldElem`
-polyhedron(A, b) = polyhedron(QQFieldElem, A, b)
-polyhedron(A) = polyhedron(QQFieldElem, A)
+# default scalar type: guess from input and fall back to `QQFieldElem`
+polyhedron(A, b) = polyhedron(_guess_fieldelem_type(A,b), A, b)
+polyhedron(A) = polyhedron(_guess_fieldelem_type(A), A)
 
 @doc raw"""
     polyhedron(P::Polymake.BigObject)
@@ -226,7 +226,7 @@ function convex_hull(f::scalar_type_or_field, V::AbstractCollection[PointVector]
   end
 end
 
-convex_hull(V::AbstractCollection[PointVector], R::Union{AbstractCollection[RayVector], Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = convex_hull(QQFieldElem, V, R, L; non_redundant=non_redundant)
+convex_hull(V::AbstractCollection[PointVector], R::Union{AbstractCollection[RayVector], Nothing} = nothing, L::Union{AbstractCollection[RayVector], Nothing} = nothing; non_redundant::Bool = false) = convex_hull(_guess_fieldelem_type(V,R,L), V, R, L; non_redundant=non_redundant)
 
 ###############################################################################
 ###############################################################################

src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl

@@ -13,11 +13,11 @@ struct SubdivisionOfPoints{T} <: PolyhedralObject{T}
 end
 
 
-# default scalar type: `QQFieldElem`
+# default scalar type: guess from input, fallback to QQ
 subdivision_of_points(points::AbstractCollection[PointVector], cells) =
-  subdivision_of_points(QQFieldElem, points, cells)
+  subdivision_of_points(_guess_fieldelem_type(points), points, cells)
 subdivision_of_points(points::AbstractCollection[PointVector], weights::AbstractVector) =
-  subdivision_of_points(QQFieldElem, points, weights)
+  subdivision_of_points(_guess_fieldelem_type(points), points, weights)
 
 # Automatic detection of corresponding OSCAR scalar type;
 # Avoid, if possible, to increase type stability

src/PolyhedralGeometry/helpers.jl

@@ -233,8 +233,9 @@ end
 
 function augment(vec::AbstractVector, val)
     s = size(vec)
+    @req s[1] > 0 "cannot homogenize empty vector"
     res = similar(vec, (s[1] + 1,))
-    res[1] = val
+    res[1] = val + zero(first(vec))
     res[2:end] = vec
     return assure_vector_polymake(res)
 end
@@ -407,47 +408,76 @@ function _embedded_quadratic_field(r::ZZRingElem)
     end
 end
 
-function _find_parent_field(::Type{T}, x, y...) where T <: scalar_types
-    f = _find_parent_field(T, x)
-    elem_type(f) == T && return f
-    return _find_parent_field(T, y...)
-end
-function _find_parent_field(::Type{T}, x::AbstractArray{<:FieldElem}) where T <: scalar_types
-    for el in x
-        el isa T && return parent(el)
+function _check_field_polyhedral(::Type{T}) where T
+  @req !(T <: NumFieldElem) "Number fields must be embedded, e.g. via `embedded_number_field`."
+  @req hasmethod(isless, (T, T)) "Field must be ordered and have `isless` method."
+end
+
+_find_elem_type(x::Nothing) = Any
+_find_elem_type(x::Any) = typeof(x)
+_find_elem_type(x::Type) = x
+_find_elem_type(x::Polymake.Rational) = QQFieldElem
+_find_elem_type(x::Polymake.Integer) = ZZRingElem
+_find_elem_type(x::AbstractArray) = reshape(_find_elem_type.(x), :)
+_find_elem_type(x::Tuple) = vcat(_find_elem_type.(x)...)
+_find_elem_type(x::AbstractArray{<:AbstractArray}) = vcat(_find_elem_type.(x)...)
+_find_elem_type(x::MatElem) = elem_type(base_ring(x))
+
+function _guess_fieldelem_type(x...)
+  types = filter(!=(Any), vcat(_find_elem_type.(x)...))
+  T = QQFieldElem
+  for t in types
+    if t == Float64
+      return Float64
+    elseif promote_type(t, T) != T
+      T = t
     end
-    return QQ
-end
-function _find_parent_field(::Type{T}, x::MatElem{<:FieldElem}) where T <: scalar_types
-    f = base_ring(x)
-    elem_type(f) == T && return f
-    return QQ
-end
-_find_parent_field(::Type{T}, x::Tuple{<:AnyVecOrMat, <:Any}) where T <: scalar_types = _find_parent_field(T, x...)
-_find_parent_field(::Type{T}, x::FieldElem) where T <: scalar_types = x isa T ? parent(x) : QQ
-_find_parent_field(::Type{T}, x::Number) where T <: scalar_types = QQ
-_find_parent_field(::Type{T}) where T <: scalar_types = QQ
-# _find_parent_field() = QQ
-_find_parent_field(::Type{T}, x::AbstractArray{<:AbstractArray}) where T <: scalar_types = _find_parent_field(T, x...)
-function _find_parent_field(::Type{T}, x::AbstractArray) where T <: scalar_types
-    for el in x
-        el isa T && return parent(el)
+  end
+  return T
+end
+
+_parent_or_coefficient_field(::Type{Float64}, x...) = AbstractAlgebra.Floats{Float64}()
+_parent_or_coefficient_field(::Type{ZZRingElem}, x...) = ZZ
+
+_parent_or_coefficient_field(r::Base.RefValue{<:Union{FieldElem, ZZRingElem}}, x...) = parent(r.x)
+_parent_or_coefficient_field(v::AbstractArray{T}) where T<:Union{FieldElem, ZZRingElem} = _parent_or_coefficient_field(T, v)
+
+# QQ is done as a special case here to avoid ambiguities
+function _parent_or_coefficient_field(::Type{T}, e::Any) where T<:FieldElem
+  hasmethod(parent, (typeof(e),)) && elem_type(parent(e)) <: T ?
+    parent(e) :
+    missing
+end
+
+function _parent_or_coefficient_field(::Type{T}, c::MatElem) where T<:FieldElem
+  elem_type(base_ring(c)) <: T ? base_ring(c) : missing
+end
+
+function _parent_or_coefficient_field(::Type{T}, c::AbstractArray) where T<:FieldElem
+  first([collect(skipmissing(_parent_or_coefficient_field.(Ref(T), c))); missing])
+end
+
+function _parent_or_coefficient_field(::Type{T}, x, y...) where T <: scalar_types
+  for c in (x, y...)
+    p = _parent_or_coefficient_field(T, c)
+    if p !== missing && elem_type(p) <: T
+      return p
     end
-    return QQ
+  end
+  missing
+end
+
+function _determine_parent_and_scalar(f::Union{Field, ZZRing}, x...)
+  _check_field_polyhedral(elem_type(f))
+  return (f, elem_type(f))
 end
 
-_determine_parent_and_scalar(f::Union{Field, ZZRing}, x...) = (f, elem_type(f))
-# isempty(x) => standard/trivial field?
 function _determine_parent_and_scalar(::Type{T}, x...) where T <: scalar_types
-    if T == QQFieldElem
-        f = QQ
-    elseif T == Float64
-        f = AbstractAlgebra.Floats{Float64}()
-    else
-        pf = _find_parent_field(T, x...)
-        f = pf == QQ ? throw(ArgumentError("Scalars of type $T require specification of a parent field. Please pass the desired Field instead of the type or have a $T contained in your input data.")) : pf
-    end
-    return (f, T)
+  T == QQFieldElem && return (QQ, QQFieldElem)
+  p = _parent_or_coefficient_field(T, x...)
+  @req p !== missing "Scalars of type $T require specification of a parent field. Please pass the desired Field instead of the type or have a $T contained in your input data."
+  _check_field_polyhedral(elem_type(p))
+  return (p, elem_type(p))
 end
 
 function _detect_default_field(::Type{Hecke.EmbeddedNumFieldElem{AbsSimpleNumFieldElem}}, p::Polymake.BigObject)
@@ -539,10 +569,6 @@ function _promote_scalar_field(a::AbstractArray{<:FieldElem})
     return _promote_scalar_field(parent.(a)...)
 end
 
-_parent_or_coefficient_field(r::Base.RefValue{<:Union{FieldElem, ZZRingElem}}) = parent(r.x)
-
-_parent_or_coefficient_field(v::AbstractVector{T}) where T<:Union{FieldElem, ZZRingElem} = _determine_parent_and_scalar(T, v)[1]
-
 function _promoted_bigobject(::Type{T}, obj::PolyhedralObject{U}) where {T <: scalar_types, U <: scalar_types}
   T == U ? pm_object(obj) : Polymake.common.convert_to{_scalar_type_to_polymake(T)}(pm_object(obj))
 end
@@ -561,6 +587,15 @@ function Polymake._fieldelem_to_float(e::EmbeddedNumFieldElem)
    return Float64(real(embedding(parent(e))(data(e), 32)))
 end
 
+function Polymake._fieldelem_to_float(e::QQBarFieldElem)
+   return Float64(ArbField(64)(e))
+end
+
+function Polymake._fieldelem_from_rational(::QQBarField, r::Rational{BigInt})
+  return QQBarFieldElem(QQFieldElem(r))
+end
+
+
 # convert a Polymake.BigObject's scalar from QuadraticExtension to OscarNumber (Polytope only)
 
 function _polyhedron_qe_to_on(x::Polymake.BigObject, f::Field)

src/PolyhedralGeometry/iterators.jl

@@ -48,7 +48,8 @@ for (T, _t) in ((:PointVector, :point_vector), (:RayVector, :ray_vector))
 
     Base.BroadcastStyle(::Type{<:$T}) = Broadcast.ArrayStyle{$T}()
 
-    _parent_or_coefficient_field(po::$T) = coefficient_field(po)
+    _parent_or_coefficient_field(::Type{TT}, po::$T{<:TT}) where TT <: FieldElem = coefficient_field(po)
+    _find_elem_type(po::$T) = elem_type(coefficient_field(po))
 
     function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$T}}, ::Type{ElType}) where ElType<:scalar_types_extended
       e = bc.f(first.(bc.args)...)
@@ -159,6 +160,9 @@ Return the `$($Hlin)` `H(a)`, which is given by a vector `a` such that
 
     coefficient_field(h::$Habs) = base_ring(h.a)
 
+    _find_elem_type(h::$Habs) = elem_type(coefficient_field(h))
+    _parent_or_coefficient_field(::Type{T}, h::$Habs{<:T}) where T <: FieldElem = coefficient_field(h)
+
   end
 
 end

src/PolyhedralGeometry/visualization.jl

@@ -1,13 +1,13 @@
 @doc raw"""
-    visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}, SubdivisionOfPoints{T}}) where T<:Union{QQFieldElem, Float64, EmbeddedNumFieldElem}
+    visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}, SubdivisionOfPoints{T}}) where T<:Union{FieldElem, Float64}
 
 Visualize a polyhedral object of dimension at most four (in 3-space).
 In dimensions up to 3 a usual embedding is shown.
 Four-dimensional polytopes are visualized as a Schlegel diagram, which is a projection onto one of the facets; e.g., see Chapter 5 of [Zie95](@cite).
 
 In higher dimensions there is no standard method; use projections to lower dimensions or try ideas from [GJRW10](@cite).
 """
-function visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}}) where T<:Union{QQFieldElem, Float64, EmbeddedNumFieldElem}
+function visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, PolyhedralComplex{T}}) where T<:Union{Float64, FieldElem}
   d = ambient_dim(P)
   b = P isa Polyhedron
   if b && d == 4
@@ -25,7 +25,7 @@ function visualize(P::Union{Polyhedron{T}, Cone{T}, PolyhedralFan{T}, Polyhedral
   Polymake.visual(pmo)
 end
 
-function visualize(P::SubdivisionOfPoints{T}) where T<:Union{QQFieldElem, Float64, EmbeddedNumFieldElem}
+function visualize(P::SubdivisionOfPoints{T}) where T<:Union{FieldElem, Float64}
   d = ambient_dim(P)
   @req d <= 3 "Can not visualize $(typeof(P)) of ambient dimension $d. Supported range: 1 <= d <= 3"
   # polymake will by default use 0:n-1 as labels so we assign labels

test/PolyhedralGeometry/extended.jl

@@ -157,9 +157,9 @@
         @test polyhedron(matrix(QQ, [-1 0 0; 0 -1 0; 0 0 -1]), [0, 0, 0]) == Pos_poly
 
         # testing different input types
-        @test polyhedron([-1 0 0; 0 -1 0; 0 0 -1], Float64[0, 0, 0]) == Pos_poly
-        @test polyhedron(Float64[-1 0 0; 0 -1 0; 0 0 -1], [0, 0, 0]) == Pos_poly
-        
+        @test polyhedron([-1 0 0; 0 -1 0; 0 0 -1], Float64[0, 0, 0]) isa Polyhedron{Float64}
+        @test polyhedron(Float64[-1 0 0; 0 -1 0; 0 0 -1], [0, 0, 0]) isa Polyhedron{Float64}
+
         let y = convex_hull([0, 0, 0], [1, 0, 0], [[0, 1, 0], [0, 0, 1]])
             @test polyhedron([-1 0 0], [0]) == y
             @test polyhedron([-1 0 0], 0) == y

test/PolyhedralGeometry/scalar_types.jl

@@ -22,7 +22,11 @@
     pq2 = E(roots(z^2 - (3 - r^2))[1])
     p = (pq1 + pq2)//2
     r = E(r)
-
+
+    e = one(QQBarFieldElem)
+    s, c = sinpi(2*e/5), cospi(2*e/5)
+    pts = [ [e, e], [s,c], [c,s] ]
+
     v = [0 1 p; 0 -1 p; r 0 q; -r 0 q; 0 r -q; 0 -r -q; 1 0 -p; -1 0 -p]
     V = matrix(E, v)
     sd = convex_hull(E, v; non_redundant = true)
@@ -48,6 +52,12 @@
     nfc =  polyhedral_fan(E, maximal_cones(IncidenceMatrix,nf), rays(nf))
     @test is_regular(nfc)
 
+    @test convex_hull(pts) isa Polyhedron{QQBarFieldElem}
+    qp = convex_hull(pts)
+
+    @test number_of_vertices(qp) == 3
+    @test number_of_facets(qp) == 3
+
     @testset "Scalar detection" begin
         let j = johnson_solid(12)
             @test j isa Polyhedron{QQFieldElem}