Opened 6 months ago
Last modified 2 months ago
#30699 new enhancement
Improve triangulate and volume methods for polytopes over inexact rings
Reported by: | gh-LaisRast | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-9.4 |
Component: | geometry | Keywords: | PointConfiguration, triangulation, volume |
Cc: | gh-kliem, jipilab, mkoeppe, dimpase, slelievre | Merged in: | |
Authors: | Reviewers: | ||
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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:
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.
Change History (21)
comment:1 Changed 6 months ago by
comment:2 Changed 6 months ago by
this is a notoriously tricky problem. I'm inclined to close it as "won't fix, use exact arithmetic instead".
comment:3 follow-up: ↓ 5 Changed 6 months ago by
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.
comment:4 Changed 6 months ago by
- Milestone changed from sage-9.2 to sage-wishlist
- Type changed from defect to enhancement
comment:5 in reply to: ↑ 3 Changed 6 months ago by
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 cdd
s 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-17so indeed some simplices there are very thin.
comment:6 Changed 6 months ago by
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.
comment:7 Changed 6 months ago by
- Description modified (diff)
comment:8 Changed 6 months ago by
- Cc slelievre added
- Description modified (diff)
comment:9 Changed 6 months ago by
- Description modified (diff)
comment:10 Changed 6 months ago by
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)
comment:11 Changed 6 months ago by
I propose to repurpose this ticket to:
- Improve
triangulate
andvolume
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()
comment:12 Changed 6 months ago by
- Description modified (diff)
- Milestone changed from sage-wishlist to sage-9.3
- Summary changed from bug in triangulation of PointConfiguration to Improve triangulate and volume methods for polytopes over inexact rings
I changed the summary and milestone. Revert if you disagree.
comment:13 Changed 6 months ago by
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)
comment:14 follow-up: ↓ 15 Changed 6 months ago by
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).
comment:15 in reply to: ↑ 14 Changed 6 months ago by
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).
comment:16 follow-up: ↓ 17 Changed 6 months ago by
‘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?
comment:17 in reply to: ↑ 16 ; follow-up: ↓ 18 Changed 6 months ago by
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.
comment:18 in reply to: ↑ 17 Changed 6 months ago by
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.
comment:19 Changed 6 months ago by
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.
comment:20 Changed 6 months ago by
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:
- Treat all coordinates as exact coordinates: You will likely break facets. You can do so, by changing the base ring to
QQ
. - 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?
comment:21 Changed 2 months ago by
- Milestone changed from sage-9.3 to sage-9.4
Setting new milestone based on a cursory review of ticket status, priority, and last modification date.
Ok. Traced it back:
in
geometry/triangulation/point_triangulation.py
. I'm not sure, what is the best alternative. I replaced0
by0.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?