.is_inscribed() makes a bad assumption in Polyhedron

[jipilab]

There is an error in the code for is_inscribed, where one assume that any d neighbors of a vertex necessarily span the space. This is wrong.

The following (most likely not minimal example) reproduces this:

sage: P = Polyhedron(vertices=[(-130658298093891402635075/416049251842505144482473,
....: 177469511761879509172000/1248147755527515433447419,
....: 485550543257132133136169/2496295511055030866894838,
....: 2010744967797898733758669/2496295511055030866894838),
....: (-146945725603929909850/706333405676769433081,
....: -84939725782618445000/706333405676769433081,
....: 560600045283000988081/1412666811353538866162,
....: 969778382942371268081/1412666811353538866162),
....: (-46275018824497300/140422338198040641,
....: -5747688262110000/46807446066013547, 1939357556329/7033601552658,
....: 1939357556329/7033601552658), (-17300/59929, -10000/59929, 39929/119858,
....: 39929/119858), (-4700/32209, -10000/32209, 12209/64418, 12209/64418),
....: (QQ(0), QQ(0), QQ(0), QQ(1)), (QQ(0), QQ(0), 1/2, 1/2), (300/10027,
....: -10000/30081, 10081/60162, 10081/60162), (112393975400/1900567733649,
....: 117311600000/633522577883, 43678681/95197362, 43678681/95197362),
....: (6109749955400/133380598418321, 37106807920000/133380598418321,
....: 2677964249/6680888498, 2677964249/6680888498),
....: (29197890764005600/402876806828660641,
....: -2150510776960000/402876806828660641,
....: 398575785274740641/805753613657321282,
....: 398575785274740641/805753613657321282),
....: (5576946899441759759983005325/110078073300232813237456943251,
....: -29071211718677797926570478000/110078073300232813237456943251,
....: 59439312069347378584317232001/220156146600465626474913886502,
....: 181346577228466312205473034501/220156146600465626474913886502),
....: (150040732779124914266530235300/6774574358246204311268446913881,
....: -2813827375989039189507000218000/6774574358246204311268446913881,
....: 1260217414021285074925933133881/13549148716492408622536893827762,
....: 3232518047094242684574253773881/13549148716492408622536893827762),
....: (3816349407976279597850158016285000/88842127448735433741180809504357161,
....: 27965821247423216557301387453968000/88842127448735433741180809504357161,
....: 68546256000224819256028677086357161/177684254897470867482361619008714322,
....: 86062257922545755787315412690197161/177684254897470867482361619008714322)])
sage: P.is_inscribed()
Traceback (most recent call last):
...
ZeroDivisionError: rational division by zero

The matrix a has zero determinant and causes this error. One should pick the neighbors more carefully to fix this.

[gh-LaisRast]

With out getting deeply into the algorithm, a naive solution to this problem is to iterate over all possible subsets (of size dimension) of vertex_neighbors, and test if the resulting polytope (convex hull of vertex and that subset of neighbors) is a simplex. But this however can be costly.

-            vertex_neighbors = vertex.neighbors()
+            vertex_neighbors = list(vertex.neighbors())

         # The following simplex is full-dimensional because `self` is assumed
-        # to be: every vertex has at least `dimension` neighbors and they form
-        # a full simplex with `vertex`.
-        simplex_vertices = [vertex] + [next(vertex_neighbors) for i in range(dimension)]
+        # to be: every `vertex` has a set of `dimension` neighbors, such that
+        # the convex hull of that `vertex` and that set of neighbors is
+        # a full-dimensional simplex.
+        for C in Combinations(range(1,dimension+1)):
+            simplex_vertices = [vertex] + [vertex_neighbors[i] for i in C]
+            if Polyhedron(simplex_vertices).is_full_dimensional():
+                break

Just because I was testing it, this is what I believe the smallest example where the current implementation of the method is_inscribed() fails. This is a pyramid over a pyramid over a square (with a specific ordering of the vertices)

sage: V =[
....:  [0, -1, 0, 0],
....:  [0, 0, -1, 0],
....:  [0, 0, 0, -1],
....:  [0, 0, +1, 0],
....:  [0, 0, 0, +1],
....:  [+1, 0, 0, 0]
....:  ]
....: P = Polyhedron(V)
....: P.is_inscribed()
Traceback (most recent call last):
...
ZeroDivisionError: rational division by zero

[tmonteil]

Replying to gh-LaisRast:

This is indeed overkill: a set of vectors form a matroid, so you can select a maximal independent set of vectors greedily, leading to a linear-time algorithm.

[tmonteil]

New commits:

​12a15b8

#28464 : fix the construction of a simplex of vertices in Polyhedron_base.is_inscribed

[jipilab]

Salut,

Looks good to me! I also like the second example, which is simpler. I was thinking about trying to reproduce the error with a smaller one, but then it all depends on the internal order given to vertices... But it is good in any case to have that one too.

I set this to positive review. Error fixed!

Merci pour la réparation rapide!

[chapoton]

moving milestone to 9.0 (after release of 8.9)