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]

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.

[mo271]

comment:15

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..

[sagetrac-git]

Changed commit from 3e1ebc0 to 85bb29e

[sagetrac-git]

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

85bb29e

added doctests

[mo271]

Author: Moritz Firsching

[jplab]

comment:19

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

[sagetrac-git]

Changed commit from 85bb29e to 455e83e

[sagetrac-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

[mo271]

comment:22

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

[jplab]

comment:23

Hi Moritz,

There are a couple of missing spaces between binary operator.

Let's see what the bot thinks.

[jplab]

comment:24

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.

[jplab]

Reviewer: Jean-Philippe Labbé

[sagetrac-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

[sagetrac-git]

Changed commit from 455e83e to 736f8a3

[mo271]

comment:27

Thank for your review, JP!

[vbraun]

Changed branch from u/moritz/21270 to 736f8a3