Polyhedron RDF face lattice bug / intersection of polyhedra

[mkoeppe]

As noted first on #20570, current Sage is not able to plot some LPs with irrational data. This appears to be due to a bug in the Sage polyhedron code for RDF data (which InteractiveLP uses when the data are not rational).

            sage: poly = polytopes.regular_polygon(7)
            sage: lp, x = poly.to_linear_program(solver='InteractiveLP', return_variable=True)
            sage: lp.set_objective(x[0] + x[1])
            sage: b = lp.get_backend()
            sage: P = b.interactive_lp_problem()
            sage: p = P.plot() ### ERROR
            sage: p.show()

CC: @fchapoton @mo271 @jplab

Component: geometry

Keywords: polyhedron, days84

Author: Moritz Firsching

Branch/Commit: 736f8a3

Reviewer: Jean-Philippe Labbé

Issue created by migration from https://trac.sagemath.org/ticket/21270

[mo271]

comment:4

misleading comment deleted

[mo271]

Changed keywords from none to polyhedron, days84

[mo271]

comment:6

Apperently something goes wrong when taking the intersection of a polyhedron in RDF and one in QQ:
Q plays the role of box in the function plot_feasible_set of interactive_simplex_method.py.

sage: L=[[1349811595/3385908, -488029882/7345133],
....:  [1349811595/3385908, 1349811595/5078862],
....:  [-732044823/7345133, 1349811595/5078862],
....:  [-732044823/7345133, -488029882/7345133]]
....: 
sage: 
sage: Q=Polyhedron(L)
sage: P=Polyhedron(polytopes.regular_polygon(7).vertices_list(), base_ring=RDF)
sage: P.intersection(Q)
A -1-dimensional polyhedron in RDF^2 defined as the convex hull of 7 vertices
sage: Q.intersection(P)
The empty polyhedron in RDF^2

[mo271]

comment:7

Here is a smaller example of the failing intersection between QQ and RDF polyhedra:

sage: Q=Polyhedron(ieqs=[[-499999,1000000],[1499999,-1000000]])
sage: P=Polyhedron(ieqs=[[0,1.0],[1.0,-1.0]], base_ring=RDF)
sage: Q.intersection(P)
The empty polyhedron in RDF^1
sage: P.intersection(Q)
A 1-dimensional polyhedron in RDF^1 defined as the convex hull of 2 vertices
sage: Q=Polyhedron(ieqs=[[-4999999,10000000],[14999999,-10000000]])
sage: P.intersection(Q)

The last command then raises
ValueError: *Error: Numerical inconsistency is found. Use the GMP exact arithmetic., which is better than a wrong result, I guess. In any case one should perhaps do the conversions of the base ring more explicitly.

The relevant function in base.py reads:

        new_ieqs = self.inequalities() + other.inequalities()
        new_eqns = self.equations() + other.equations()
        parent = self.parent()
        try:
            return parent.element_class(parent, None, [new_ieqs, new_eqns])
        except TypeError as msg:
            if self.base_ring() is ZZ:
                parent = parent.base_extend(QQ)
                return parent.element_class(parent, None, [new_ieqs, new_eqns])
            else:
                raise TypeError(msg)

So in the example above we have in the new_ieqs variable:

(An inequality (1.0) x + 0.0 >= 0,
 An inequality (-1.0) x + 1.0 >= 0,
 An inequality (-10000000) x + 14999999 >= 0,
 An inequality (10000000) x - 4999999 >= 0)

Which leads to the unexpected behaviour

[mo271]

Branch: u/moritz/21270

[mo271]

comment:9

This commit fixes the instances of the problem, that I added.

The original thing is still not working, but this is now another bug (more later).

---

New commits:

83d5e8f

normalizing for rational input, when converting polyhedra to RDF

[mo271]

Commit: 83d5e8f

[sagetrac-git]

Changed commit from 83d5e8f to 3e1ebc0

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

3e1ebc0

use the correct base ring

[mo271]

comment:11

Now it works..

[mo271]

Description changed:

--- 
+++ 
@@ -6,7 +6,7 @@
             sage: lp.set_objective(x[0] + x[1])
             sage: b = lp.get_backend()
             sage: P = b.interactive_lp_problem()
-            sage: p = P.plot()    ### ERROR
-            sage: p.show()
+            sage: p = P.plot(xmin=-3, xmax=8, ymin=-2, ymax=4) ### ERROR
+            sage: p.show(dpi=400)
 ```
 

[mo271]

comment:13

I submitted a fix to this bug.

Explanation:
What went wrong was was happens if a polyhedron in QQ is intersection with a polyhedron in RDF; basically it takes the union of the inequalities. The inequalities are typically stored as linear inequality with coefficients in ZZ. Those coefficients can get very large and cause numerical trouble when converted to RDF. Thus it is advisable to convert the inequalities to RDF only after some renormalization. We choose to divide all terms of the inequality by the absolute values of the coefficient with the largest absolute value. (In case all the coefficients are zero, we leave it as is.)

This helped with the examples above. In the original question there was another problem: in the method _solve of the class InteractiveLPProblem( in the interactive_simplex_method.py file, taking the base ring to be something other than the base ring of the feasible region (which might have been coerced to something else) caused trouble. I changed this by simply setting

R = F.base_ring()

Then I was able to produce the plot.