Improve triangulate and volume methods for polytopes over inexact rings

[gh-LaisRast]

The volume and triangulate methods for polytopes over inexact rings can fail with ZeroDivisionError.

For example, define the following polytope:

sage: pts = [(-0.1113445378, -0.55567226889999999, 1),
....:        (1.0063025210000001, -0.49684873950000003, 0),
....:        (1, 0, 1),
....:        (2, 0, 0),
....:        (-0.56836734690000001, -0.13673469390000001, 1),
....:        (-0.53979591839999996, 0.92040816329999997, 0),
....:        (0, 1, 1),
....:        (0, 2, 0)]
sage: P = Polyhedron(backend='cdd', base_ring=RDF, vertices=pts)

Trying to triangulate it or to compute its volume, we get a zero division error:

sage: P.volume()
Traceback (most recent call last)
...
ZeroDivisionError: input matrix must be nonsingular

sage: P.triangulate()
Traceback (most recent call last)
...
ZeroDivisionError: input matrix must be nonsingular

The problem happens inside the method placing_triangulation of PointConfiguration. Specifically, in the following line:

        simplices = [frozenset(P.contained_simplex(large=True))]

This line produces a list with different ordering in different sage sessions. So you may need to run the above code more than one time to reproduce the problem.

(The bug was reported by Günter Rote. It appeared when he was trying to compute the volume of the polytope P given above.)

Original discussion:

• ​sage-devel, 2020-10-02, Polytope triangulation fails sometimes on polytopes with non-precise vertices

To get the volume, a workaround is suggested by comment:3 below:

sage: P.change_ring(QQ).volume().n()
1.98224843509406

This ticket proposes to catch the ZeroDivisionError and fail with in informative error message and a suggestion to work over the rationals.

[gh-kliem]

Ok. Traced it back:

2044                 v = point.reduced_affine_vector() - origin.reduced_affine_vector()
2045                 if v*normal>0:
2046                     visible_facets.append(facet)

in geometry/triangulation/point_triangulation.py. I'm not sure, what is the best alternative. I replaced 0 by 0.000001 and it works just fine.

How is one supposed to check this? Is there a good way to check whether the scalar product of two vectors is positive, which considers their norm etc. in case of inexactness?

[dimpase]

this is a notoriously tricky problem. I'm inclined to close it as "won't fix, use exact arithmetic instead".

[dimpase]

here is a link to CGAL manual (CGAL was developed by computational geometry specialists, no reasons to expect we can beat them): ​https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation_3VariabilityDependingonthe

In a nutshell, you need exact predicates, even if you work with floating point data. Meanwhile:

sage: tmpRDF=RDF                                                                                                                                               
sage: RDF=QQ                                                                                                                                                   
sage: P = Polyhedron(backend='cdd', base_ring=RDF, vertices=[(-RDF(0.1113445378), -RDF(0.55567226889999999), RDF(1)), (RDF(1.0063025210000001), -RDF(0.49684873
....: 950000003), RDF(0)), (RDF(1), RDF(0), RDF(1)), (RDF(2), RDF(0), RDF(0)), (-RDF(0.56836734690000001), -RDF(0.13673469390000001), RDF(1)), (-RDF(0.53979591
....: 839999996), RDF(0.92040816329999997), RDF(0)), (RDF(0), RDF(1), RDF(1)), (RDF(0), RDF(2), RDF(0))]) 
....: sage: P.triangulate()                                                                                                                                    
(<0,1,2,3>, <0,1,3,4>, <0,2,3,6>, <0,3,4,6>, <1,3,4,7>, <1,4,5,7>, <3,4,6,7>, <4,5,6,7>)
sage: P                                                                                                                                                        
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 8 vertices
sage: P.volume()                                                                                                                                               
239866285104588490583706258416397108079995566644827095107690547050507/121007175921017411055114023964368753091725419133201274175573961383650
sage: P.volume().n()                                                                                                                                           
1.98224843509406
sage: RDF=tmpRDF
sage: Polyhedron(P.vertices()[0:4]).volume().n()                                                                                                               
1.62487829158265e-17      

so indeed some simplices there are very thin.

[gh-kliem]

Ok. I agree with you that if we really want to support triangulation for reals, then we should do via some other package. I also don't think this is a priority. (Note that I solved #29176 by using cdds incidence matrix, because our way of figuring when a scalar is zero, is just not as good as theirs.)

In the meantime: How about changing backend to field (inexpensive) for computing the volume as you suggested?

I really don't like the traceback and it is not documented that we do not expect this to work. So naturally people will file a bug report.

Replying to dimpase:

here is a link to CGAL manual (CGAL was developed by computational geometry specialists, no reasons to expect we can beat them): ​https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation_3VariabilityDependingonthe

In a nutshell, you need exact predicates, even if you work with floating point data. Meanwhile:

sage: tmpRDF=RDF                                                                                                                                               
sage: RDF=QQ                                                                                                                                                   
sage: P = Polyhedron(backend='cdd', base_ring=RDF, vertices=[(-RDF(0.1113445378), -RDF(0.55567226889999999), RDF(1)), (RDF(1.0063025210000001), -RDF(0.49684873
....: 950000003), RDF(0)), (RDF(1), RDF(0), RDF(1)), (RDF(2), RDF(0), RDF(0)), (-RDF(0.56836734690000001), -RDF(0.13673469390000001), RDF(1)), (-RDF(0.53979591
....: 839999996), RDF(0.92040816329999997), RDF(0)), (RDF(0), RDF(1), RDF(1)), (RDF(0), RDF(2), RDF(0))]) 
....: sage: P.triangulate()                                                                                                                                    
(<0,1,2,3>, <0,1,3,4>, <0,2,3,6>, <0,3,4,6>, <1,3,4,7>, <1,4,5,7>, <3,4,6,7>, <4,5,6,7>)
sage: P                                                                                                                                                        
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 8 vertices
sage: P.volume()                                                                                                                                               
239866285104588490583706258416397108079995566644827095107690547050507/121007175921017411055114023964368753091725419133201274175573961383650
sage: P.volume().n()                                                                                                                                           
1.98224843509406
sage: RDF=tmpRDF
sage: Polyhedron(P.vertices()[0:4]).volume().n()                                                                                                               
1.62487829158265e-17      

so indeed some simplices there are very thin.

[dimpase]

Documenting this, perhaps printing a warning, is a good idea. To me, the whole idea of doing inexact computations like this in a naive way is a bit problematic. IMHO it's a legacy of the times where such computations were too slow in exact arithmetic.

[slelievre]

More on the thin simplices in comment:3.

Start with the points from the initial report:

sage: pts = [(-0.1113445378, -0.55567226889999999, 1),
....:        (1.0063025210000001, -0.49684873950000003, 0),
....:        (1, 0, 1),
....:        (2, 0, 0),
....:        (-0.56836734690000001, -0.13673469390000001, 1),
....:        (-0.53979591839999996, 0.92040816329999997, 0),
....:        (0, 1, 1),
....:        (0, 2, 0)]

Construct the corresponding polytope over RDF:

sage: P = Polyhedron(backend='cdd', base_ring=RDF, vertices=pts)

Plotting it reveals a polyhedron with six quadrilateral faces:

sage: p = P.plot(threejs_flat_shading=True)
sage: p.show(frame=False)
Launched html viewer for Graphics3d Object

Changing the base ring to the rationals splits some of the faces:

sage: Q = P.change_ring(QQ)
sage: q = Q.plot(threejs_flat_shading=True)
sage: q.show(frame=False)

The simplex from comment:3 corresponds to one of those split faces:

sage: S = Polyhedron(Q.vertices()[0:4])
sage: S.volume().n()
1.62487829158265e-17
sage: s = S.plot(threejs_flat_shading=True)
sage: s.show(frame=False)

To see all the simplices, first triangulate:

sage: T = Q.triangulate()
sage: T
(<0,1,2,3>, <0,1,3,4>, <0,2,3,6>, <0,3,4,6>,
 <1,3,4,7>, <1,4,5,7>, <3,4,6,7>, <4,5,6,7>)

Plot all the simplices together:

sage: simplices = [Polyhedron([Q.Vrepresentation(i) for i in t]) for t in T]
sage: opts = {'threejs_flat_shading': True, 'alpha': 0.2}
sage: plot_simplex = lambda s: s.plot(fill=(random(), random(), random()), **opts)
sage: simplex_plots = [plot_simplex(s) for s in simplices]
sage: sum(simplex_plots).show(frame=False)

Or individually; two of them are very flat and correspond to faces of P:

sage: skeleton = Q.plot(fill=False)
sage: _ = any((skeleton + s).show() for s in simplex_plots)

[slelievre]

I propose to repurpose this ticket to:

• Improve triangulate and volume for polytopes over inexact rings

The improvement could consist in

• catching any ZeroDivisionError and fail with a more helpful error message
• for volume computations, suggest trying P.change_ring(QQ).volume().n()

[slelievre]

I changed the summary and milestone. Revert if you disagree.

[dimpase]

It makes no sense to talk about facets of polyhedron over an inexact ring.

sage: P = Polyhedron(backend='cdd', base_ring=RDF, vertices=pts)                                                                                                     
sage: P.facets()                                                                                                                                                     
(A 3-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices,
 A 3-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices,
 A 3-dimensional face of a Polyhedron in RDF^3 defined as the convex hull of 4 vertices)

3-dimensional facet of a 3-dimensional polytope, yeah, cool ?! (I guess Sage should throw an error here)

You can plot it, as per comment:10, and get an approximate picture of it, but the reality is more like it has 9 facets, only 4 of them 4-gons:

sage: P.change_ring(QQ).facets()                                                                                                                                     
(A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 4 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices,
 A 2-dimensional face of a Polyhedron in QQ^3 defined as the convex hull of 3 vertices)

[gh-kliem]

The combinatorics are fine. As mentioned earlier ‘cdd‘ computes the incidence matrix just fine (in many scenarios) and we now use this instead of recomputing. Faces should just be set up with the correct dimension from the combinatorics. That is a two line fix, which should be done.

Of course we cannot beat the backend. If the backend fails, we can't do anything (which is not the case here).

[dimpase]

Replying to gh-kliem:

The combinatorics are fine. As mentioned earlier ‘cdd‘ computes the incidence matrix just fine (in many scenarios) and we now use this instead of recomputing. Faces should just be set up with the correct dimension from the combinatorics. That is a two line fix, which should be done.

No, my point is that you can't really talk about combinatorics here - as you need a new definition of what an "approximate facet" is.

That the backend thinks that there 4 vectors with RDF coordinates are coplanar only makes sense if your affine hyperplanes are actually pairs of sufficiently close parallel hyperplanes. There are no real facets in the RDF world, only these "thick" facets. (Needless to say, vertices are also not points, they are balls of radius depending on the precision of your RDF).

[gh-kliem]

‘cdd‘ takes care of those things and figures out an incidence matrix for it.

I think we should just use it and e.g. set up faces with the dimension determined by the indices.

As for the volume, center or center of mass, it does not make a difference wether vertices are are coplanar or not.

Anyway, I don't exactly get your point. Do you disagree with the direction the ticket goes? Do you want to point out that no matter how much we work, inexact polyhedra will remain inperfect?

[dimpase]

Replying to gh-kliem:

‘cdd‘ takes care of those things and figures out an incidence matrix for it.

I doubt that cdd uses the needed extra precision etc to get a robust answer. One can potentially end up with an incidence matrix that just cannot correspond to a polytope.

I think we should just use it and e.g. set up faces with the dimension determined by the indices.

I am not sure. I'd rather be on a safe side and error out if the dimension is wrong.

As for the volume, center or center of mass, it does not make a difference wether vertices are are coplanar or not.

IMHO the volume is computed by triangulating first, and this step may fail.

[gh-kliem]

Replying to dimpase:

Replying to gh-kliem:

‘cdd‘ takes care of those things and figures out an incidence matrix for it.

I doubt that cdd uses the needed extra precision etc to get a robust answer. One can potentially end up with an incidence matrix that just cannot correspond to a polytope.

This is why we compute everything twice with cdd. If we start from the vertices, we check that feeding in the output (facets), we get somewhat the original vertices back. So we check whether the computation of cdd is at least consistent (this step already fails for many non-trivial polyhedra).

I think we should just use it and e.g. set up faces with the dimension determined by the indices.

I am not sure. I'd rather be on a safe side and error out if the dimension is wrong.

Of course the answer is wrong, if we require that points are exactly coplanar. If you don't like the answer that cdd is giving you, you just should not use it. You are free to use exact arithmetic anytime you want.

As for the volume, center or center of mass, it does not make a difference wether vertices are are coplanar or not.

IMHO the volume is computed by triangulating first, and this step may fail.

As explained above for the volume you can just change to QQ and get an answer which is approximately correct. That in QQ some of the faces break, won't matter for the volume.

[dimpase]

Does the consistency of cdd output guarantees that the answer mathematically makes sense, e.g. does not violate Euler formula?

My point is that Sage should not silently lead the user up the garden path into a place with facets of wrong dimension and all that.

[gh-kliem]

As mentioned earlier, this wrong dimension of facets is simply because we compute the dimension by the rank of a matrix corresponding to the vertices of a face. This computation assumes exactness, which cdd does not. cdd already determines the dimensions of the faces by the incidences matrix. We could use that and resolve the inconsistency.

No, the consistency does not guarantee Euler's formular to hold. It could easily end up with a combintorial type, that does not have a realization at all. This is the risk of using inexact data. There is two ways to deal with it:

1. Treat all coordinates as exact coordinates: You will likely break facets. You can do so, by changing the base ring to QQ.
2. Consider points to lie on a hyperplane, even if they are slightly off.

I think, we have already agreed on, that the first approach should be used for volume etc and that we will simply not support triangulation for the second approach (unless it accidentally works). I would even go so far, as to completely ban triangulations at least for inexact non-simplicial polytopes, because it will likely contradict with the combinatorics that cdd figured.

I would personally never work with inexact polyhedra, but I can see that they sometimes arise and sometimes you want a good guess, what this polyhedron actually is.

What is it, you are suggesting?

• Completely ban inexact polyhedra?
• Always raise a warning, no matter if we detected inconsistency or not? Something as: "inexact polyhedra may fail dramatically"
• Something else?

[mkoeppe]

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.