PolyhedralGeometry: Compute triangulations of polytope (#848)

* PolyhedralGeometry: Compute triangulations 

- all or only regular triangulations
- Input can be polytope or point set
- Backend used is TOPCOM via polymake

Co-authored-by: Benjamin Lorenz <benlorenz@users.noreply.github.com>

docs/src/PolyhedralGeometry/Polyhedra/auxiliary.md

@@ -44,6 +44,7 @@ h_vector(P::Polyhedron)
 ## Other
 
 ```@docs
+all_triangulations
 boundary_lattice_points(P::Polyhedron)
 contains(P::Polyhedron, v::AbstractVector)
 ehrhart_polynomial(P::Polyhedron)
@@ -62,6 +63,8 @@ polarize(P::Polyhedron)
 project_full(P::Polyhedron)
 print_constraints(A::AnyVecOrMat, b::AbstractVector; trivial::Bool = false)
 print_constraints(P::Polyhedron; trivial::Bool = false)
+regular_triangulations
+secondary_polytope
 solve_ineq(A::fmpz_mat, b::fmpz_mat)
 solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat, d::fmpz_mat)
 solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)

src/PolyhedralGeometry/main.jl

@@ -20,6 +20,7 @@ export AffineHalfspace,
     affine_equation_matrix,
     affine_hull,
     affine_inequality_matrix,
+    all_triangulations,
     archimedean_solid,
     ambient_dim,
     bipyramid,
@@ -111,13 +112,15 @@ export AffineHalfspace,
     project_full,
     pyramid,
     recession_cone,
+    regular_triangulations,
     relative_interior_point,
     save_cone,
     save_linearprogram,
     save_polyhedralfan,
     save_polyhedron,
     save_subdivisionofpoints,
     secondary_cone,
+    secondary_polytope,
     simplex,
     solve_lp,
     starsubdivision,
@@ -159,6 +162,7 @@ include("Groups.jl")
 include("Serialization.jl")
 include("Visualization.jl")
 include("solving_integrally.jl")
+include("triangulations.jl")
 
 # Some temporary aliases to avoid breaking all current PRs
 pm_cone(C::Cone) = pm_object(C)

src/PolyhedralGeometry/triangulations.jl

@@ -0,0 +1,158 @@
+@doc Markdown.doc"""
+    all_triangulations(pts::AnyVecOrMat)
+
+Compute all triangulations on the points given as the rows of `pts`.
+
+The return type is a `Vector{Vector{Vector{Int}}}` where each
+`Vector{Vector{Int}}` encodes a triangulation, in which a `Vector{Int}` encodes
+a simplex as the set of indices of the vertices of the simplex. I.e. the
+`Vector{Int}` `[1,2,4]` corresponds to the simplex that is the convex hull of
+the first, second, and fourth input point.
+
+# Examples
+```jldoctest
+julia> c = cube(2,0,1)
+A polyhedron in ambient dimension 2
+
+julia> V = vertices(c)
+4-element SubObjectIterator{PointVector{Polymake.Rational}}:
+ [0, 0]
+ [1, 0]
+ [0, 1]
+ [1, 1]
+
+julia> all_triangulations(V)
+2-element Vector{Vector{Vector{Int64}}}:
+ [[1, 2, 3], [2, 3, 4]]
+ [[1, 2, 4], [1, 3, 4]]
+```
+"""
+function all_triangulations(pts::Union{SubObjectIterator{<:PointVector}, AbstractMatrix, Oscar.MatElem})
+    input = homogenized_matrix(pts, 1)
+    PC = Polymake.polytope.PointConfiguration(POINTS=input)
+    triangs = Polymake.polytope.topcom_all_triangulations(PC)
+    result = [[[e for e in simplex] for simplex in triang] for triang in triangs]
+    return Polymake.to_one_based_indexing(result)
+end
+
+
+@doc Markdown.doc"""
+    all_triangulations(P::Polyhedron)
+
+Compute all triangulations that can be formed using the vertices of the given
+polytope `P`.
+
+The return type is a `Vector{Vector{Vector{Int}}}` where each
+`Vector{Vector{Int}}` encodes a triangulation, in which a `Vector{Int}` encodes
+a simplex as the set of indices of the vertices of the simplex. I.e. the
+`Vector{Int}` `[1,2,4]` corresponds to the simplex that is the convex hull of
+the first, second, and fourth input point.
+
+# Examples
+```jldoctest
+julia> c = cube(2,0,1)
+A polyhedron in ambient dimension 2
+
+julia> all_triangulations(c)
+2-element Vector{Vector{Vector{Int64}}}:
+ [[1, 2, 3], [2, 3, 4]]
+ [[1, 2, 4], [1, 3, 4]]
+```
+"""
+all_triangulations(P::Polyhedron) = all_triangulations(vertices(P))
+
+
+@doc Markdown.doc"""
+    regular_triangulations(pts::AnyVecOrMat)
+
+Compute all regular triangulations on the points given as the rows of `pts`.
+
+A triangulation is regular if it can be induced by weights, i.e. attach a
+weight to every point, take the convex hull of these new vectors and then take
+the subdivision corresponding to the facets visible from below (lower
+envelope).
+
+The return type is a `Vector{Vector{Vector{Int}}}` where each
+`Vector{Vector{Int}}` encodes a triangulation, in which a `Vector{Int}` encodes
+a simplex as the set of indices of the vertices of the simplex. I.e. the
+`Vector{Int}` `[1,2,4]` corresponds to the simplex that is the convex hull of
+the first, second, and fourth input point.
+
+# Examples
+```jldoctest
+julia> c = cube(2,0,1)
+A polyhedron in ambient dimension 2
+
+julia> V = vertices(c)
+4-element SubObjectIterator{PointVector{Polymake.Rational}}:
+ [0, 0]
+ [1, 0]
+ [0, 1]
+ [1, 1]
+
+julia> regular_triangulations(V)
+2-element Vector{Vector{Vector{Int64}}}:
+ [[1, 2, 3], [2, 3, 4]]
+ [[1, 3, 4], [1, 2, 4]]
+```
+"""
+function regular_triangulations(pts::Union{SubObjectIterator{<:PointVector}, AbstractMatrix, Oscar.MatElem})
+    input = homogenized_matrix(pts, 1)
+    PC = Polymake.polytope.PointConfiguration(POINTS=input)
+    triangs = Polymake.polytope.topcom_regular_triangulations(PC)
+    result = [[[e for e in simplex] for simplex in triang] for triang in triangs]
+    return Polymake.to_one_based_indexing(result)
+end
+
+
+@doc Markdown.doc"""
+    regular_triangulations(P::Polyhedron)
+
+Compute all regular triangulations that can be formed using the vertices of the
+given polytope `P`.
+
+A triangulation is regular if it can be induced by weights, i.e. attach a
+weight to every point, take the convex hull of these new vectors and then take
+the subdivision corresponding to the facets visible from below (lower
+envelope).
+
+The return type is a `Vector{Vector{Vector{Int}}}` where each
+`Vector{Vector{Int}}` encodes a triangulation, in which a `Vector{Int}` encodes
+a simplex as the set of indices of the vertices of the simplex. I.e. the
+`Vector{Int}` `[1,2,4]` corresponds to the simplex that is the convex hull of
+the first, second, and fourth input point.
+
+# Examples
+```jldoctest
+julia> c = cube(2,0,1)
+A polyhedron in ambient dimension 2
+
+julia> regular_triangulations(c)
+2-element Vector{Vector{Vector{Int64}}}:
+ [[1, 2, 3], [2, 3, 4]]
+ [[1, 3, 4], [1, 2, 4]]
+```
+"""
+regular_triangulations(P::Polyhedron) = regular_triangulations(vertices(P))
+
+
+@doc Markdown.doc"""
+    secondary_polytope(P::Polyhedron)
+
+Compute the secondary polytope of a polyhedron, i.e. the convex hull of all the
+gkz vectors of all its (regular) triangulations. A triangulation here means
+only using the vertices of `P`.
+
+# Examples
+Compute the secondary polytope of the cube.
+```jldoctest
+julia> c = cube(3)
+A polyhedron in ambient dimension 3
+
+julia> sc = secondary_polytope(c)
+A polyhedron in ambient dimension 8
+```
+"""
+function secondary_polytope(P::Polyhedron)
+    return Polyhedron(Polymake.polytope.secondary_polytope(pm_object(P)))
+end