Trac #20413: InteractiveLPBackend: Use standard-form transformation, objective_constant_term, change default base_ring to QQ

Follow-up on #20296, #20311, #20301.
Updating `InteractiveLPBackend` to:
 - use the standard-form transformation,
 - simplify its code slightly by making use of the new
`objective_constant_term` handling in `InteractiveLPProblem`
 - make the example of optimizing over the dodecahedron more natural and
remove use of backend methods.
 - change default `base_ring` to QQ -- a much saner default because it's
so much faster

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

Author: Release Manager

src/sage/numerical/backends/generic_backend.pyx

@@ -1169,7 +1169,7 @@ def default_mip_solver(solver = None):
 
         - ``InteractiveLPProblem`` (``solver="InteractiveLP"``).  A didactical
           implementation of the revised simplex method in Sage.  It works over
-          any exact ordered field, the default is ``AA``.
+          any exact ordered field, the default is ``QQ``.
 
         ``solver`` should then be equal to one of ``"GLPK"``,
         ``"Coin"``, ``"CPLEX"``,  ``"CVXOPT"``, ``"Gurobi"``, ``"PPL"`, or
@@ -1294,7 +1294,7 @@ cpdef GenericBackend get_solver(constraint_generation = False, solver = None, ba
 
         - ``InteractiveLPProblem`` (``solver="InteractiveLP"``).  A didactical
           implementation of the revised simplex method in Sage.  It works over
-          any exact ordered field, the default is ``AA``.
+          any exact ordered field, the default is ``QQ``.
 
         ``solver`` should then be equal to one of ``"GLPK"``, ``"Coin"``,
         ``"CPLEX"``, ``"CVXOPT"``,``"Gurobi"``, ``"PPL"``, ``"InteractiveLP"``,

src/sage/numerical/backends/interactivelp_backend.pxd

@@ -15,6 +15,7 @@ cdef class InteractiveLPBackend(GenericBackend):
     cdef object prob_name
 
     cdef object lp_std_form
+    cdef object std_form_transformation
     cdef object final_dictionary
     cdef int verbosity
 

src/sage/numerical/backends/interactivelp_backend.pyx

@@ -52,19 +52,22 @@ cdef class InteractiveLPBackend:
         This backend can work with irrational algebraic numbers::
 
             sage: poly = polytopes.dodecahedron(base_ring=AA)
-            sage: lp = poly.to_linear_program(solver='InteractiveLP')
-            sage: b = lp.get_backend()
-            sage: for k in range(3): b.variable_lower_bound(k, 0)
-            sage: b.set_objective([1, 1, 1])
+            sage: lp, x = poly.to_linear_program(solver='InteractiveLP', return_variable=True)
+            sage: lp.set_objective(x[0] + x[1] + x[2])
             sage: lp.solve()
             2.291796067500631?
-            sage: [b.get_variable_value(k) for k in range(3)]
+            sage: lp.get_values(x[0], x[1], x[2])
             [0.763932022500211?, 0.763932022500211?, 0.763932022500211?]
+            sage: lp.set_objective(x[0] - x[1] - x[2])
+            sage: lp.solve()
+            2.291796067500631?
+            sage: lp.get_values(x[0], x[1], x[2])
+            [0.763932022500211?, -0.763932022500211?, -0.763932022500211?]
         """
 
         if base_ring is None:
-            from sage.rings.all import AA
-            base_ring = AA
+            from sage.rings.all import QQ
+            base_ring = QQ
 
         self.lp = InteractiveLPProblem([], [], [], base_ring=base_ring)
         self.set_verbosity(0)
@@ -74,8 +77,6 @@ cdef class InteractiveLPBackend:
         else:
             self.set_sense(-1)
 
-        self.obj_constant_term = 0
-
         self.row_names = []
 
     cpdef base_ring(self):
@@ -91,7 +92,7 @@ cdef class InteractiveLPBackend:
             sage: from sage.numerical.backends.generic_backend import get_solver
             sage: p = get_solver(solver = "InteractiveLP")
             sage: p.base_ring()
-            Algebraic Real Field
+            Rational Field
         """
         return self.lp.base_ring()
 
@@ -196,7 +197,7 @@ cdef class InteractiveLPBackend:
             sage: p.objective_coefficient(1)
             1
         """
-        A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+        A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
         cdef int vtype = int(binary) + int(continuous) + int(integer)
         if  vtype == 0:
             continuous = True
@@ -219,7 +220,7 @@ cdef class InteractiveLPBackend:
         x = tuple(x) + (name,)
         self.lp = InteractiveLPProblem(A, b, c, x,
                                        constraint_types, variable_types,
-                                       problem_type, ring)
+                                       problem_type, ring, objective_constant_term=d)
         return self.ncols() - 1
 
     cpdef int add_variables(self, int n, lower_bound=0, upper_bound=None, binary=False, continuous=True, integer=False, obj=None, names=None) except -1:
@@ -298,23 +299,24 @@ cdef class InteractiveLPBackend:
         else:
             raise NotImplementedError()
 
-    def _AbcxCVPR(self):
+    def _AbcxCVPRd(self):
         """
         Retrieve all problem data from the LP.
 
         EXAMPLE::
 
             sage: from sage.numerical.backends.generic_backend import get_solver
             sage: p = get_solver(solver = "InteractiveLP")
-            sage: p._AbcxCVPR()
-            ([], (), (), (), (), (), 'max', Algebraic Real Field)
+            sage: p._AbcxCVPRd()
+            ([], (), (), (), (), (), 'max', Rational Field, 0)
         """
         A, b, c, x = self.lp.Abcx()
         constraint_types = self.lp.constraint_types()
         variable_types = self.lp.variable_types()
         problem_type = self.lp.problem_type()
         base_ring = self.lp.base_ring()
-        return A, b, c, x, constraint_types, variable_types, problem_type, base_ring
+        d = self.lp.objective_constant_term()
+        return A, b, c, x, constraint_types, variable_types, problem_type, base_ring, d
 
     cpdef set_sense(self, int sense):
         """
@@ -337,14 +339,14 @@ cdef class InteractiveLPBackend:
             sage: p.is_maximization()
             False
         """
-        A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+        A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
         if sense == +1:
             problem_type = "max"
         else:
             problem_type = "min"
         self.lp = InteractiveLPProblem(A, b, c, x,
                                        constraint_types, variable_types,
-                                       problem_type, ring)
+                                       problem_type, ring, objective_constant_term=d)
 
     cpdef objective_coefficient(self, int variable, coeff=None):
         """
@@ -372,12 +374,38 @@ cdef class InteractiveLPBackend:
         if coeff is None:
             return self.lp.objective_coefficients()[variable]
         else:
-            A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+            A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
             c = list(c)
             c[variable] = coeff
             self.lp = InteractiveLPProblem(A, b, c, x,
                                            constraint_types, variable_types,
-                                           problem_type, ring)
+                                           problem_type, ring, objective_constant_term=d)
+
+    cpdef objective_constant_term(self, d=None):
+        """
+        Set or get the constant term in the objective function
+
+        INPUT:
+
+        - ``d`` (double) -- its coefficient.  If `None` (default), return the current value.
+
+        EXAMPLE::
+
+            sage: from sage.numerical.backends.generic_backend import get_solver
+            sage: p = get_solver(solver = "InteractiveLP")
+            sage: p.objective_constant_term()
+            0
+            sage: p.objective_constant_term(42)
+            sage: p.objective_constant_term()
+            42
+        """
+        if d is None:
+            return self.lp.objective_constant_term()
+        else:
+            A, b, c, x, constraint_types, variable_types, problem_type, ring, _ = self._AbcxCVPRd()
+            self.lp = InteractiveLPProblem(A, b, c, x,
+                                           constraint_types, variable_types,
+                                           problem_type, ring, objective_constant_term=d)
 
     cpdef set_objective(self, list coeff, d = 0):
         """
@@ -423,12 +451,11 @@ cdef class InteractiveLPBackend:
             -47/5
 
         """
-        A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+        A, b, _, x, constraint_types, variable_types, problem_type, ring, _ = self._AbcxCVPRd()
         c = coeff
         self.lp = InteractiveLPProblem(A, b, c, x,
                                        constraint_types, variable_types,
-                                       problem_type, ring)
-        self.obj_constant_term = d
+                                       problem_type, ring, objective_constant_term=d)
 
     cpdef set_verbosity(self, int level):
         """
@@ -470,13 +497,13 @@ cdef class InteractiveLPBackend:
             sage: p.get_values([x,y])
             [0, 3]
         """
-        A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+        A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
         A = A.delete_rows((i,))
         b = list(b); del b[i]
         constraint_types=list(constraint_types); del constraint_types[i]
         self.lp = InteractiveLPProblem(A, b, c, x,
                                        constraint_types, variable_types,
-                                       problem_type, ring)
+                                       problem_type, ring, objective_constant_term=d)
 
     cpdef add_linear_constraint(self, coefficients, lower_bound, upper_bound, name=None):
         """
@@ -511,7 +538,7 @@ cdef class InteractiveLPBackend:
             sage: p.row_name(1)
             'foo'
         """
-        A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+        A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
         if lower_bound is None:
            if upper_bound is None:
                raise ValueError("At least one of lower_bound and upper_bound must be provided")
@@ -537,7 +564,7 @@ cdef class InteractiveLPBackend:
 
         self.lp = InteractiveLPProblem(A, b, c, x,
                                        constraint_types, variable_types,
-                                       problem_type, ring)
+                                       problem_type, ring, objective_constant_term=d)
 
 
     cpdef add_col(self, list indices, list coeffs):
@@ -633,7 +660,8 @@ cdef class InteractiveLPBackend:
         """
         ## FIXME: standard_form should allow to pass slack names (which we would take from row_names).
         ## FIXME: Perhaps also pass the problem name as objective name
-        lp_std_form = self.lp_std_form = self.lp.standard_form()
+        lp_std_form, transformation = self.lp.standard_form(transformation=True)
+        self.lp_std_form, self.std_form_transformation = lp_std_form, transformation
         output = lp_std_form.run_revised_simplex_method()
         ## FIXME: Display output as a side effect if verbosity is high enough
         d = self.final_dictionary = lp_std_form.final_revised_dictionary()
@@ -674,7 +702,7 @@ cdef class InteractiveLPBackend:
         v = d.objective_value()
         if self.lp_std_form.is_negative():
             v = - v
-        return self.obj_constant_term + v
+        return v
 
     cpdef get_variable_value(self, int variable):
         """
@@ -701,9 +729,8 @@ cdef class InteractiveLPBackend:
             sage: p.get_variable_value(1)
             3/2
         """
-        if str(self.lp.decision_variables()[variable]) != str(self.lp_std_form.decision_variables()[variable]):
-            raise NotImplementedError("Undoing the standard-form transformation is not implemented")
-        return self.final_dictionary.basic_solution()[variable]
+        solution = self.std_form_transformation(self.final_dictionary.basic_solution())
+        return solution[variable]
 
     cpdef int ncols(self):
         """
@@ -1027,12 +1054,12 @@ cdef class InteractiveLPBackend:
         else:
             if value != bounds[1]:
                 bounds = (bounds[0], value)
-                A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+                A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
                 variable_types = list(variable_types)
                 variable_types[index] = self._variable_type_from_bounds(*bounds)
                 self.lp = InteractiveLPProblem(A, b, c, x,
                                                constraint_types, variable_types,
-                                               problem_type, ring)
+                                               problem_type, ring, objective_constant_term=d)
 
     cpdef variable_lower_bound(self, int index, value = False):
         """
@@ -1071,12 +1098,12 @@ cdef class InteractiveLPBackend:
         else:
             if value != bounds[0]:
                 bounds = (value, bounds[1])
-                A, b, c, x, constraint_types, variable_types, problem_type, ring = self._AbcxCVPR()
+                A, b, c, x, constraint_types, variable_types, problem_type, ring, d = self._AbcxCVPRd()
                 variable_types = list(variable_types)
                 variable_types[index] = self._variable_type_from_bounds(*bounds)
                 self.lp = InteractiveLPProblem(A, b, c, x,
                                                constraint_types, variable_types,
-                                               problem_type, ring)
+                                               problem_type, ring, objective_constant_term=d)
 
     cpdef bint is_variable_basic(self, int index):
         """

src/sage/numerical/mip.pyx

@@ -598,7 +598,8 @@ cdef class MixedIntegerLinearProgram(SageObject):
             sage: p = MixedIntegerLinearProgram(solver='ppl')
             sage: p.base_ring()
             Rational Field
-            sage: p = MixedIntegerLinearProgram(solver='InteractiveLP')
+            sage: from sage.rings.all import AA
+            sage: p = MixedIntegerLinearProgram(solver='InteractiveLP', base_ring=AA)
             sage: p.base_ring()
             Algebraic Real Field
             sage: d = polytopes.dodecahedron()