add examples for non-rational polyhedra

docs/src/PolyhedralGeometry/intro.md

@@ -27,16 +27,9 @@ OSCAR, which can be found [here](https://www.oscar-system.org/tutorials/Polyhedr
 
 The objects from polyhedral geometry operate on a given type, which (usually) resembles a field. This is indicated by the template parameter, e.g. the properties of a `Polyhedron{QQFieldElem}` are rational numbers of type `QQFieldElem`, if applicable.
 Supported scalar types are `FieldElem` and `Float64`, but some functionality might not work properly if the parent `Field` does not satisfy certain mathematic conditions, like being ordered.
-When constructing a polyhedral object from scratch, for the "simpler" types `QQFieldElem` and `Float64` it suffices to pass the `Type`, but more complex `FieldElem`s require a parent `Field` object. This can be set by either passing the desired `Field` instead of the type, or by inserting the type and have a matching `FieldElem` in your input data. If no type or field is given, the scalar type defaults to `QQFieldElem`.
+When constructing a polyhedral object from scratch, for the "simpler" types `QQFieldElem` and `Float64` it suffices to pass the `Type`, but more complex `FieldElem`s require a parent `Field` object. This can be set by either passing the desired `Field` instead of the type, or by inserting the type and have a matching `FieldElem` in your input data. If no type or field is given, the scalar type will be deduced from the input data and defaults to `QQFieldElem` for primitive types like `Int64`.
 
-The parent `Field` of the coefficients of an object `O` with coefficients of type `T` can be retrieved with the `coefficient_field` function, and it holds `elem_type(coefficient_field(O)) == T`.
-
-```@docs
-coefficient_field(x::PolyhedralObject)
-```
-
-!!! warning
-    Support for fields other than the rational numbers is currently in an experimental stage.
+The parent `Field` of the coefficients of an object `O` with coefficients of type `T` can be retrieved with the [`coefficient_field`](@ref) function, and it holds `elem_type(coefficient_field(O)) == T`.
 
 These three lines result in the same polytope over rational numbers. Besides the general support mentioned above, naming a `Field` explicitly is encouraged because it allows user control and increases efficiency.
 ```jldoctest
@@ -48,7 +41,49 @@ true
 
 julia> P == convex_hull([1 0 0; 0 0 1]) # `Field` defaults to `QQ`
 true
+```
 
+Working with polytopes over algebraic number fields requires an embedded number field. This can be constructed with [`embedded_number_field`](@ref) from a polynomial and a choice for the root.
+
+```jldoctest
+julia> Qx, x = QQ[:x];
+
+julia> F,a = embedded_number_field(x^2-3//4, 0.866)
+(Embedded field
+Number field of degree 2 over QQ
+at
+Real embedding with 0.87 of number field, a)
+
+julia> triangle = convex_hull(F, [0 0; 1 0; 1//2 a])
+Polyhedron in ambient dimension 2 with EmbeddedAbsSimpleNumFieldElem type coefficients
+
+julia> volume(triangle)
+1//2*a (0.43)
+```
+
+The algebraic closure of the rational numbers also comes with an embedding.
+
+```jldoctest
+julia> K = algebraic_closure(QQ);
+
+julia> h = sqrt(K(3//4))
+{a2: 0.866025}
+
+julia> mat = matrix(K, [0 0; 1 0; 1//2 h])
+[    {a1: 0}       {a1: 0}]
+[ {a1: 1.00}       {a1: 0}]
+[{a1: 0.500}   {a2: 0.866}]
+
+julia> triangle = convex_hull(mat)
+Polyhedron in ambient dimension 2 with QQBarFieldElem type coefficients
+
+julia> volume(triangle)
+{a2: 0.433013}
+```
+
+```@docs
+coefficient_field(x::PolyhedralObject)
+embedded_number_field
 ```
 
 ## Type compatibility

src/AlgebraicGeometry/ToricVarieties/NormalToricVarieties/properties.jl

@@ -120,7 +120,7 @@ false
 ```
 """
 @attr Bool has_torusfactor(v::NormalToricVarietyType) =
-  Polymake.common.rank(rays(v)) < ambient_dim(v)
+  Polymake.common.rank(matrix_for_polymake(rays(v))) < ambient_dim(v)
 
 @doc raw"""
     is_orbifold(v::NormalToricVarietyType)

src/PolyhedralGeometry/Cone/constructors.jl

@@ -154,11 +154,11 @@ function cone_from_inequalities(
   non_redundant::Bool=false,
 )
   parent_field, scalar_type = _determine_parent_and_scalar(f, I, E)
-  IM = -linear_matrix_for_polymake(I)
+  IM = -linear_matrix_for_polymake(parent_field, I)
   EM = if isnothing(E) || isempty(E)
     Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2))
   else
-    linear_matrix_for_polymake(E)
+    linear_matrix_for_polymake(parent_field, E)
   end
 
   if non_redundant
@@ -209,7 +209,7 @@ function cone_from_equations(
   non_redundant::Bool=false,
 )
   parent_field, scalar_type = _determine_parent_and_scalar(f, E)
-  EM = linear_matrix_for_polymake(E)
+  EM = linear_matrix_for_polymake(parent_field, E)
   IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2))
   return cone_from_inequalities(f, IM, EM; non_redundant=non_redundant)
 end

src/PolyhedralGeometry/Cone/properties.jl

@@ -13,7 +13,7 @@ _ray_cone(U::Type{RayVector{T}}, C::Cone{T}, i::Base.Integer) where {T<:scalar_t
   ray_vector(coefficient_field(C), view(pm_object(C).RAYS, i, :))::U
 
 _vector_matrix(::Val{_ray_cone}, C::Cone; homogenized=false) =
-  homogenized ? homogenize(pm_object(C).RAYS, 0) : pm_object(C).RAYS
+  homogenized ? homogenize(coefficient_field(C), pm_object(C).RAYS, 0) : pm_object(C).RAYS
 
 _matrix_for_polymake(::Val{_ray_cone}) = _vector_matrix
 
@@ -609,7 +609,11 @@ _lineality_cone(
   ray_vector(coefficient_field(C), view(pm_object(C).LINEALITY_SPACE, i, :))::U
 
 _generator_matrix(::Val{_lineality_cone}, C::Cone; homogenized=false) =
-  homogenized ? homogenize(pm_object(C).LINEALITY_SPACE, 0) : pm_object(C).LINEALITY_SPACE
+  if homogenized
+    homogenize(coefficient_field(C), pm_object(C).LINEALITY_SPACE, 0)
+  else
+    pm_object(C).LINEALITY_SPACE
+  end
 
 _matrix_for_polymake(::Val{_lineality_cone}) = _generator_matrix
 

src/PolyhedralGeometry/Optimization/linear_program.jl

@@ -30,18 +30,19 @@ function linear_program(
     throw(ArgumentError("convention must be set to :min or :max."))
   end
   ambDim = ambient_dim(P)
-  objective = coefficient_field(P).(objective)
+  cf = coefficient_field(P)
+  objective = cf.(objective)
   size(objective, 1) == ambDim || error("objective has wrong dimension.")
   lp = Polymake.polytope.LinearProgram{_scalar_type_to_polymake(T)}(;
-    LINEAR_OBJECTIVE=homogenize(objective, k)
+    LINEAR_OBJECTIVE=homogenize(cf, objective, k)
   )
   if convention == :max
     Polymake.attach(lp, "convention", "max")
   elseif convention == :min
     Polymake.attach(lp, "convention", "min")
   end
   Polymake.add(pm_object(P), "LP", lp)
-  LinearProgram{T}(P, lp, convention, coefficient_field(P))
+  LinearProgram{T}(P, lp, convention, cf)
 end
 
 linear_program(

src/PolyhedralGeometry/Optimization/mixed_integer_linear_program.jl

@@ -37,17 +37,19 @@ function mixed_integer_linear_program(
     integer_variables = 1:ambDim
   end
   size(objective, 1) == ambDim || error("objective has wrong dimension.")
-  objective = coefficient_field(P).(objective)
+  cf = coefficient_field(P)
+  objective = cf.(objective)
   milp = Polymake.polytope.MixedIntegerLinearProgram{_scalar_type_to_polymake(T)}(;
-    LINEAR_OBJECTIVE=homogenize(objective, k), INTEGER_VARIABLES=Vector(integer_variables)
+    LINEAR_OBJECTIVE=homogenize(cf, objective, k),
+    INTEGER_VARIABLES=Vector(integer_variables),
   )
   if convention == :max
     Polymake.attach(milp, "convention", "max")
   elseif convention == :min
     Polymake.attach(milp, "convention", "min")
   end
   Polymake.add(pm_object(P), "MILP", milp)
-  MixedIntegerLinearProgram{T}(P, milp, convention, coefficient_field(P))
+  MixedIntegerLinearProgram{T}(P, milp, convention, cf)
 end
 
 function mixed_integer_linear_program(

src/PolyhedralGeometry/PolyhedralComplex/constructors.jl

@@ -82,11 +82,11 @@ function polyhedral_complex(
   non_redundant::Bool=false,
 )
   parent_field, scalar_type = _determine_parent_and_scalar(f, vr, L)
-  points = homogenized_matrix(vr, 1)
+  points = homogenized_matrix(parent_field, vr, 1)
   LM = if isnothing(L) || isempty(L)
-    zero_matrix(QQ, 0, size(points, 2))
+    zero_matrix(parent_field, 0, size(points, 2))
   else
-    homogenized_matrix(L, 0)
+    homogenized_matrix(parent_field, L, 0)
   end
 
   # Rays and Points are homogenized and combined and
@@ -191,12 +191,12 @@ function polyhedral_complex(F::PolyhedralFan{T}) where {T<:scalar_types}
     if Polymake.exists(pmo_in, prop)
       inprop = Polymake.give(pmo_in, prop)
       outprop = if inprop isa Polymake.IncidenceMatrix
-        Polymake.IncidenceMatrix(nrows(inprop), ncols(inprop) + 1)
+        Polymake.IncidenceMatrix(Polymake.nrows(inprop), Polymake.ncols(inprop) + 1)
       else
-        Polymake.SparseMatrix(nrows(inprop), ncols(inprop) + 1)
+        Polymake.spzeros(Int, Polymake.nrows(inprop), Polymake.ncols(inprop) + 1)
       end
       outprop[1:end, 2:end] = inprop
-      for i in 1:nrows(inprop)
+      for i in 1:Polymake.nrows(inprop)
         outprop[i, 1] = 1
       end
       Polymake.take(pmo_out, prop, outprop)

src/PolyhedralGeometry/PolyhedralComplex/properties.jl

@@ -463,7 +463,7 @@ julia> n_maximal_polyhedra(PC)
 2
 ```
 """
-n_maximal_polyhedra(PC::PolyhedralComplex) = pm_object(PC).N_MAXIMAL_POLYTOPES
+n_maximal_polyhedra(PC::PolyhedralComplex) = pm_object(PC).N_MAXIMAL_POLYTOPES::Int
 
 @doc raw"""
     is_simplicial(PC::PolyhedralComplex)
@@ -559,8 +559,10 @@ julia> for p in P1s
 function polyhedra_of_dim(
   PC::PolyhedralComplex{T}, polyhedron_dim::Int
 ) where {T<:scalar_types}
+  t = Polyhedron{T}
   n = polyhedron_dim - lineality_dim(PC) + 1
-  n < 0 && return nothing
+  polyhedron_dim > dim(PC) && return _empty_subobjectiterator(t, PC)
+  n <= 0 && return _empty_subobjectiterator(t, PC)
   pfaces = Polymake.fan.cones_of_dim(pm_object(PC), n)
   nfaces = Polymake.nrows(pfaces)
   rfaces = Vector{Int64}()

src/PolyhedralGeometry/PolyhedralFan/properties.jl

@@ -76,7 +76,11 @@ _ray_fan(U::Type{RayVector{T}}, PF::_FanLikeType, i::Base.Integer) where {T<:sca
   ray_vector(coefficient_field(PF), view(pm_object(PF).RAYS, i, :))::U
 
 _vector_matrix(::Val{_ray_fan}, PF::_FanLikeType; homogenized=false) =
-  homogenized ? homogenize(pm_object(PF).RAYS, 0) : pm_object(PF).RAYS
+  if homogenized
+    homogenize(coefficient_field(PF), pm_object(PF).RAYS, 0)
+  else
+    pm_object(PF).RAYS
+  end
 
 _matrix_for_polymake(::Val{_ray_fan}) = _vector_matrix
 
@@ -213,7 +217,7 @@ julia> cones(PF, 2)
 function cones(PF::_FanLikeType, cone_dim::Int)
   l = cone_dim - length(lineality_space(PF))
   t = Cone{_get_scalar_type(PF)}
-  (l < 0 || dim(PF) == -1) && return _empty_subobjectiterator(t, PF)
+  (l < 0 || dim(PF) == -1 || cone_dim > dim(PF)) && return _empty_subobjectiterator(t, PF)
 
   if l == 0
     return SubObjectIterator{t}(
@@ -739,7 +743,11 @@ _lineality_fan(
   ray_vector(coefficient_field(PF), view(pm_object(PF).LINEALITY_SPACE, i, :))::U
 
 _generator_matrix(::Val{_lineality_fan}, PF::_FanLikeType; homogenized=false) =
-  homogenized ? homogenize(pm_object(PF).LINEALITY_SPACE, 0) : pm_object(PF).LINEALITY_SPACE
+  if homogenized
+    homogenize(coefficient_field(PF), pm_object(PF).LINEALITY_SPACE, 0)
+  else
+    pm_object(PF).LINEALITY_SPACE
+  end
 
 _matrix_for_polymake(::Val{_lineality_fan}) = _generator_matrix
 

src/PolyhedralGeometry/Polyhedron/constructors.jl

@@ -139,14 +139,14 @@ function polyhedron(
 )
   parent_field, scalar_type = _determine_parent_and_scalar(f, I, E)
   if isnothing(I) || _isempty_halfspace(I)
-    EM = affine_matrix_for_polymake(E)
+    EM = affine_matrix_for_polymake(parent_field, E)
     IM = Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(EM, 2))
   else
-    IM = -affine_matrix_for_polymake(I)
+    IM = -affine_matrix_for_polymake(parent_field, I)
     EM = if isnothing(E) || _isempty_halfspace(E)
       Polymake.Matrix{_scalar_type_to_polymake(scalar_type)}(undef, 0, size(IM, 2))
     else
-      affine_matrix_for_polymake(E)
+      affine_matrix_for_polymake(parent_field, E)
     end
   end
 
@@ -259,14 +259,30 @@ function convex_hull(
   L::Union{AbstractCollection[RayVector],Nothing}=nothing;
   non_redundant::Bool=false,
 )
+  if V isa AbstractVector
+    eltype(V) <: Union{AffineHyperplane,AffineHalfspace,LinearHalfspace,LinearHyperplane} &&
+      throw(ArgumentError("Cannot use `convex_hull` for halfspace of hyperplane input"))
+    eltype(V) <: RayVector && throw(ArgumentError("First argument must not contain rays"))
+  end
+  if R isa AbstractVector
+    eltype(R) <: Union{AffineHyperplane,AffineHalfspace,LinearHalfspace,LinearHyperplane} &&
+      throw(ArgumentError("Cannot use `convex_hull` for halfspace of hyperplane input"))
+    eltype(R) <: PointVector &&
+      throw(ArgumentError("Second argument must not contain points"))
+  end
+
   parent_field, scalar_type = _determine_parent_and_scalar(f, V, R, L)
   # Rays and Points are homogenized and combined and
   # Lineality is homogenized
-  points = stack(homogenized_matrix(V, 1), homogenized_matrix(R, 0))
+  points = stack(
+    parent_field,
+    homogenized_matrix(parent_field, V, 1),
+    homogenized_matrix(parent_field, R, 0),
+  )
   lineality = if isnothing(L) || isempty(L)
-    zero_matrix(QQ, 0, size(points, 2))
+    zero_matrix(parent_field, 0, size(points, 2))
   else
-    homogenized_matrix(L, 0)
+    homogenized_matrix(parent_field, L, 0)
   end
 
   # These matrices are in the right format for polymake.

src/PolyhedralGeometry/Polyhedron/standard_constructions.jl

@@ -181,12 +181,13 @@ julia> vertices(P)
 ```
 """
 function orbit_polytope(V::AbstractCollection[PointVector], G::PermGroup)
-  Vhom = stack(homogenized_matrix(V, 1), nothing)
+  parent_field, scalar_type = _determine_parent_and_scalar(_guess_fieldelem_type(V), V)
+  Vhom = homogenized_matrix(parent_field, V, 1)
   @req size(Vhom, 2) == degree(G) + 1 "Dimension of points and group degree need to be the same"
   generators = PermGroup_to_polymake_array(G)
   pmGroup = Polymake.group.PermutationAction(; GENERATORS=generators)
   pmPolytope = Polymake.polytope.orbit_polytope(Vhom, pmGroup)
-  return Polyhedron{QQFieldElem}(pmPolytope)
+  return Polyhedron{scalar_type}(pmPolytope)
 end
 
 @doc raw"""
@@ -2398,20 +2399,19 @@ vectors.
 The following creates the vertices of a parallelogram with the origin as its vertex barycenter.
 ```jldoctest
 julia> zonotope_vertices_fukuda_matrix([1 1 0; 1 1 1])
-pm::Matrix<pm::Rational>
--1 -1 -1/2
-0 0 -1/2
-0 0 1/2
-1 1 1/2
+[-1   -1   -1//2]
+[ 0    0   -1//2]
+[ 0    0    1//2]
+[ 1    1    1//2]
 ```
 """
 function zonotope_vertices_fukuda_matrix(M::Union{MatElem,AbstractMatrix})
-  A = M
-  if eltype(M) <: Union{Integer,ZZRingElem}
-    A = Polymake.Matrix{Polymake.Rational}(convert(Polymake.PolymakeType, M))
-  end
-  return Oscar.dehomogenize(
-    Polymake.polytope.zonotope_vertices_fukuda(Oscar.homogenized_matrix(A, 1))
+  parent_field, _ = _determine_parent_and_scalar(_guess_fieldelem_type(M), M)
+  return matrix(
+    parent_field,
+    dehomogenize(
+      Polymake.polytope.zonotope_vertices_fukuda(homogenized_matrix(parent_field, M, 1))
+    ),
   )
 end
 

src/PolyhedralGeometry/SubdivisionOfPoints/constructors.jl

@@ -62,9 +62,9 @@ function subdivision_of_points(
   f::scalar_type_or_field, points::AbstractCollection[PointVector], cells::IncidenceMatrix
 )
   @req size(points)[1] == ncols(cells) "Number of points must be the same as columns of IncidenceMatrix"
-  hpts = homogenized_matrix(points, 1)
-  @req allunique(eachrow(hpts)) "Points must be unique"
   parent_field, scalar_type = _determine_parent_and_scalar(f, points)
+  hpts = homogenized_matrix(parent_field, points, 1)
+  @req allunique(eachrow(hpts)) "Points must be unique"
   arr = Polymake.@convert_to Array{Set{Int}} Polymake.common.rows(cells)
   SubdivisionOfPoints{scalar_type}(
     Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}(;
@@ -105,9 +105,9 @@ function subdivision_of_points(
   f::scalar_type_or_field, points::AbstractCollection[PointVector], weights::AbstractVector
 )
   @req size(points)[1] == length(weights) "Number of points must equal number of weights"
-  hpts = homogenized_matrix(points, 1)
-  @req allunique(eachrow(hpts)) "Points must be unique"
   parent_field, scalar_type = _determine_parent_and_scalar(f, points, weights)
+  hpts = homogenized_matrix(parent_field, points, 1)
+  @req allunique(eachrow(hpts)) "Points must be unique"
   SubdivisionOfPoints{scalar_type}(
     Polymake.fan.SubdivisionOfPoints{_scalar_type_to_polymake(scalar_type)}(;
       POINTS=hpts, WEIGHTS=weights