Trac #20415: Polyhedron.to_linear_program should select solver by bas…

…e_ring

This is a follow-up on #20301.
#20296 provides a way to select a suitable solver by its `base_ring`.

`Polyhedron.to_linear_program` should do that.

Note, for exact polyhedra, this will change the behavior of
`to_linear_program` from inexact computation to exact computation, which
will be much slower (but note, speed is irrelevant for this method,
because when the `Polyhedron` was set up, already the double description
has been computed eagerly.)

To restore the old behavior, the user would need to pass an explicit
`base_ring=RDF` or `solver=...` argument to `to_linear_program`.
(The patch adds many examples that illustrate this.)

URL: http://trac.sagemath.org/20415
Reported by: mkoeppe
Ticket author(s): Matthias Koeppe
Reviewer(s): Dima Pasechnik

src/sage/geometry/polyhedron/base.py

@@ -1032,20 +1032,25 @@ def n_lines(self):
         """
         return len(self.lines())
 
-    def to_linear_program(self, solver=None, return_variable=False):
+    def to_linear_program(self, solver=None, return_variable=False, base_ring=None):
         r"""
-        Return the polyhedron as a :class:`MixedIntegerLinearProgram`.
+        Return a linear optimization problem over the polyhedron in the form of
+        a :class:`MixedIntegerLinearProgram`.
 
         INPUT:
 
-        - ``solver`` -- select a solver (data structure). See the documentation
+        - ``solver`` -- select a solver (MIP backend). See the documentation
           of for :class:`MixedIntegerLinearProgram`. Set to ``None`` by default.
 
         - ``return_variable`` -- (default: ``False``) If ``True``, return a tuple
           ``(p, x)``, where ``p`` is the :class:`MixedIntegerLinearProgram` object
           and ``x`` is the vector-valued MIP variable in this problem, indexed
           from 0.  If ``False``, only return ``p``.
 
+        - ``base_ring`` -- select a field over which the linear program should be
+          set up.  Use ``RDF`` to request a fast inexact (floating point) solver
+          even if ``self`` is exact.
+
         Note that the :class:`MixedIntegerLinearProgram` object will have the
         null function as an objective to be maximized.
 
@@ -1055,17 +1060,56 @@ def to_linear_program(self, solver=None, return_variable=False):
             polyhedron associated with a :class:`MixedIntegerLinearProgram`
             object.
 
-        EXAMPLE::
+        EXAMPLES:
+
+        Exact rational linear program::
 
             sage: p = polytopes.cube()
             sage: p.to_linear_program()
             Mixed Integer Program  ( maximization, 3 variables, 6 constraints )
             sage: lp, x = p.to_linear_program(return_variable=True)
             sage: lp.set_objective(2*x[0] + 1*x[1] + 39*x[2])
             sage: lp.solve()
-            42.0
+            42
             sage: lp.get_values(x[0], x[1], x[2])
-            [1.0, 1.0, 1.0]
+            [1, 1, 1]
+
+        Floating-point linear program::
+
+            sage: lp, x = p.to_linear_program(return_variable=True, base_ring=RDF)
+            sage: lp.set_objective(2*x[0] + 1*x[1] + 39*x[2])
+            sage: lp.solve()
+            42.0
+
+        Irrational algebraic linear program over an embedded number field::
+
+            sage: p=polytopes.icosahedron()
+            sage: lp, x = p.to_linear_program(return_variable=True)
+            sage: lp.set_objective(x[0] + x[1] + x[2])
+            sage: lp.solve()
+            1/4*sqrt5 + 3/4
+
+        Same example with floating point::
+
+            sage: lp, x = p.to_linear_program(return_variable=True, base_ring=RDF)
+            sage: lp.set_objective(x[0] + x[1] + x[2])
+            sage: lp.solve() # tol 1e-5
+            1.3090169943749475
+
+        Same example with a specific floating point solver::
+
+            sage: lp, x = p.to_linear_program(return_variable=True, solver='GLPK')
+            sage: lp.set_objective(x[0] + x[1] + x[2])
+            sage: lp.solve() # tol 1e-8
+            1.3090169943749475
+
+        Irrational algebraic linear program over `AA`::
+
+            sage: p=polytopes.icosahedron(base_ring=AA)
+            sage: lp, x = p.to_linear_program(return_variable=True)
+            sage: lp.set_objective(x[0] + x[1] + x[2])
+            sage: lp.solve()
+            1.309016994374948?
 
         TESTS::
 
@@ -1077,17 +1121,13 @@ def to_linear_program(self, solver=None, return_variable=False):
             sage: p.to_linear_program(solver='PPL')
             Traceback (most recent call last):
             ...
-            NotImplementedError: Cannot use PPL on exact irrational data.
+            TypeError: The PPL backend only supports rational data.
         """
-        from sage.rings.rational_field import QQ
-        R = self.base_ring()
-        if (solver is not None and
-            solver.lower() == 'ppl' and
-            R.is_exact() and (not R == QQ)):
-            raise NotImplementedError('Cannot use PPL on exact irrational data.')
-
+        if base_ring is None:
+            base_ring = self.base_ring()
+        base_ring = base_ring.fraction_field()
         from sage.numerical.mip import MixedIntegerLinearProgram
-        p = MixedIntegerLinearProgram(solver=solver)
+        p = MixedIntegerLinearProgram(solver=solver, base_ring=base_ring)
         x = p.new_variable(real=True, nonnegative=False)
 
         for ineqn in self.inequalities_list():

src/sage/numerical/backends/generic_backend.pyx

@@ -1301,7 +1301,13 @@ cpdef GenericBackend get_solver(constraint_generation = False, solver = None, ba
         or ``None``. If ``solver=None`` (default),
         the default solver is used (see ``default_mip_solver`` method).
 
-    - ``base_ring`` -- Request a solver that works over this field.
+    - ``base_ring`` -- If not ``None``, request a solver that works over this
+        (ordered) field.  If ``base_ring`` is not a field, its fraction field
+        is used.
+
+        For example, is ``base_ring=ZZ`` is provided, the solver will work over
+        the rational numbers.  This is unrelated to whether variables are
+        constrained to be integers or not.
 
     - ``constraint_generation`` -- Only used when ``solver=None``.
 
@@ -1325,6 +1331,8 @@ cpdef GenericBackend get_solver(constraint_generation = False, solver = None, ba
         Real Double Field
         sage: p = get_solver(base_ring=QQ); p
         <sage.numerical.backends.ppl_backend.PPLBackend object at ...>
+        sage: p = get_solver(base_ring=ZZ); p
+        <sage.numerical.backends.ppl_backend.PPLBackend object at ...>
         sage: p.base_ring()
         Rational Field
         sage: p = get_solver(base_ring=AA); p
@@ -1346,6 +1354,7 @@ cpdef GenericBackend get_solver(constraint_generation = False, solver = None, ba
         solver = default_mip_solver()
 
         if base_ring is not None:
+            base_ring = base_ring.fraction_field()
             from sage.rings.all import QQ, RDF
             if base_ring is QQ:
                 solver = "Ppl"
@@ -1384,7 +1393,7 @@ cpdef GenericBackend get_solver(constraint_generation = False, solver = None, ba
 
     elif solver == "Ppl":
         from sage.numerical.backends.ppl_backend import PPLBackend
-        return PPLBackend()
+        return PPLBackend(base_ring=base_ring)
 
     elif solver == "Interactivelp":
         from sage.numerical.backends.interactivelp_backend import InteractiveLPBackend

src/sage/numerical/backends/ppl_backend.pyx

@@ -41,15 +41,29 @@ cdef class PPLBackend(GenericBackend):
     # Common denominator for objective function in self.mip (not for the constant term)
     cdef Integer obj_denominator
 
-    def __cinit__(self, maximization = True):
+    def __cinit__(self, maximization = True, base_ring = None):
         """
         Constructor
 
         EXAMPLE::
 
             sage: p = MixedIntegerLinearProgram(solver = "PPL")
+
+        TESTS:
+
+        Raise an error if a ``base_ring`` is requested that is not supported::
+
+            sage: p = MixedIntegerLinearProgram(solver = "PPL", base_ring=AA)
+            Traceback (most recent call last):
+            ...
+            TypeError: The PPL backend only supports rational data.
         """
 
+        if base_ring is not None:
+            from sage.rings.all import QQ
+            if base_ring is not QQ:
+                raise TypeError('The PPL backend only supports rational data.')
+
         self.Matrix = []
         self.row_lower_bound = []
         self.row_upper_bound = []