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()

[moritz]

misleading comment deleted

[moritz]

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

[moritz]

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

[moritz]

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

[git]

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

​3e1ebc0

use the correct base ring

[moritz]

Now it works..

[moritz]

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.

[moritz]

By the way a better plotting result is obtained by

sage: p = P.plot(xmin=-3, xmax=8, ymin=-2, ymax=4) 
sage: p.show(dpi=500)

But giving an explicit bounding box circumvents the first problem..

[git]

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

​85bb29e

added doctests

[jipilab]

There seems to be a merge conflict. Could you rebase it?

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​3274eca	normalizing for rational input, when converting polyhedra to RDF
​ef00c9a	use the correct base ring
​54ff7b8	added doctests
​455e83e	indentation of docstring

[moritz]

This should now merge cleanly. (There was just an issue with some trailing whitespace that I had removed in a line that was changed in rc0.)

[jipilab]

Hi Moritz,

There are a couple of missing spaces between binary operator.

Let's see what the bot thinks.

[jipilab]

Hi Moritz,

The bot is happy and it looks good to me. If you correct the 3 spacing around operators in the parent.py, you can set it to positive review on my behalf.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​bc1dd36	normalizing for rational input, when converting polyhedra to RDF
​abdfa7e	use the correct base ring
​1ca25df	added doctests
​6c628ed	indentation of docstring
​736f8a3	rebase and whitespace

[moritz]

Thank for your review, JP!