cut generation: additional options, query instance features

docs/bibliography.rst

@@ -1,9 +1,13 @@
 .. rubric::
 
+.. [AGY05] K. Andersen, G. Cornuéjols and Y. Li, Reduce-and-Split Cuts: Improving the Performance of Mixed-Integer Gomory Cuts. *Management Science*, 51, 1593-1732, 2005.
+
 .. [Ande73] L.G. Anderson: A Simulation Study of some Dynamic Channel Assignment Algorithms in a High Capacity Mobile Telecommunications System. *IEEE Transactions on Communications*, 21:1294--1301. 1973.
 
 .. [Appleg06] D.L. Applegate, R.E. Bixby, V. Chvatal and W.J. Cook. *The traveling salesman problem: a computational study*. Princeton university press. 2006.
 
+.. [AtRa02] A Atamtürk, D Rajan. On splittable and unsplittable flow capacitated network design arc–set polyhedra. *Mathematical Programming*, 92(2), 315-333, 2002.
+
 .. [Bala75] E. Balas. Facets of the knapsack polytope. *Mathematical Programming*, 8, 146-164, 1975.
 
 .. [BCC93] E.Balas, S. Ceria and G. Cornuéjols. A lift-and-project cutting plane algorithm for mixed 0–1 programs. *Mathematical Programming*, 58, 295-324, 1993.
@@ -22,6 +26,8 @@
 
 .. [Fisch14] M. Fischetti, M. Monaci. Exploiting erraticism in search. *Operations Research*, 62:114--122, 2014.
 
+.. [FiSa11] M. Fischetti, and D. Salvagnin. A relax-and-cut framework for Gomory mixed-integer cuts. *Mathematical Programming Computation*, 3(2), 79-102, 2011.
+
 .. [Gomo69] R.E. Gomory. Some polyhedra related to combinatorial problems. *Linear Algebra and its Applications*, 2(4):451-558, 1969.
 
 .. [Junge09] M. Jünger, T.M. Liebling, D. Naddef, G. Nemhauser, W.R. Pulleyblank, G. Reinelt, G. Rinaldi and Wolsey, L.A. *50 Years of integer programming 1958-2008: From the early years to the state-of-the-art*. Springer Science & Business Media. 2009.

docs/classes.rst

@@ -112,5 +112,7 @@ Useful functions
 .. autofunction:: mip.minimize
 .. autofunction:: mip.maximize
 .. autofunction:: mip.xsum
+.. autofunction:: mip.compute_features
+.. autofunction:: mip.features
 
 

examples/extract_features_mip.py

@@ -0,0 +1,25 @@
+"""This example reads a MIP (in .lp or .mps), solves its linear programming
+relaxation and then tests the impact of adding different types of cutting
+planes."""
+
+import sys
+from mip import Model, features, compute_features
+import mip
+
+lp_path = ""
+
+#  using test data, replace with your instance
+lp_path = mip.__file__.replace("mip/__init__.py", "test/data/1443_0-9.lp").replace(
+    "mip\\__init__.py", "test\\data\\1443_0-9.lp"
+)
+
+m = Model()
+if m.solver_name.upper() in ["GRB", "GUROBI"]:
+    print("This feature is only supported in CBC.")
+else:
+    m.read(lp_path)
+
+    print("instance features:")
+    X = compute_features(m)
+    for i, fn in enumerate(features()):
+        print("%s: %g" % (fn, X[i]))

examples/gen_cuts.py

@@ -0,0 +1,46 @@
+"""Example of use of general cutting planes in a small MIP."""
+
+import sys
+from mip import Model, INTEGER, maximize, CutType, OptimizationStatus
+
+larger_diff = -1
+best_cut = None
+
+for ct in CutType:
+    print("Trying cut type: {}".format(ct.name))
+
+    m = Model()
+    if m.solver_name.upper() in ["GRB", "GUROBI"]:
+        print("This feature is currently not supported in Gurobi.")
+    else:
+        m.verbose = 0
+
+        x = m.add_var_tensor(shape=(2, 1), name="x", var_type=INTEGER)
+
+        m.objective = maximize(2 * x[0] + x[1])
+
+        m += 7 * x[0] + x[1] <= 28
+        m += -x[0] + 3 * x[1] <= 7
+        m += -8 * x[0] - 9 * x[1] <= -32
+
+        m.optimize(relax=True)
+        olr = m.objective_value
+
+        cp = m.generate_cuts([ct])
+        if cp and cp.cuts:
+            print("{} cuts generated:".format(len(cp.cuts)))
+            for c in cp.cuts:
+                print("  " + str(c))
+
+        if cp.cuts:
+            m += cp
+            m.optimize(relax=True)
+
+            print("Dual bound now: {}".format(m.objective_value))
+            assert m.status == OptimizationStatus.OPTIMAL
+            diff = m.objective_value - olr
+            if diff > larger_diff:
+                larger_diff = diff
+                best_cut = ct
+
+print("Best cut: {}".format(best_cut))

examples/gen_cuts_mip.py

@@ -0,0 +1,64 @@
+"""This example reads a MIP (in .lp or .mps), solves its linear programming
+relaxation and then tests the impact of adding different types of cutting
+planes. In the end, the it informs which cut generator produced the best bound
+improvement."""
+
+from textwrap import shorten
+import sys
+from mip import Model, CutType, OptimizationStatus
+import mip
+
+lp_path = ""
+
+#  using test data
+lp_path = mip.__file__.replace("mip/__init__.py", "test/data/1443_0-9.lp").replace(
+    "mip\\__init__.py", "test\\data\\1443_0-9.lp"
+)
+
+m = Model()
+if m.solver_name.upper() in ["GRB", "GUROBI"]:
+    print("This feature is currently not supported in Gurobi.")
+else:
+    m.read(lp_path)
+
+    m.verbose = 0
+    m.optimize(relax=True)
+    print("Original LP bound: {}".format(m.objective_value))
+
+    best_impr = -1
+    best_cut = ""
+
+    for ct in CutType:
+        print()
+        m2 = m.copy()
+        m2.verbose = 0
+        m2.optimize(relax=True)
+        assert (
+            m2.status == OptimizationStatus.OPTIMAL
+            and abs(m2.objective_value - m.objective_value) <= 1e-4
+        )
+        print("Searching for violated {} inequalities ...".format(ct.name))
+        cp = m2.generate_cuts([ct])
+        if cp and cp.cuts:
+            print("{} cuts found:".format(len(cp.cuts)))
+            for c in cp.cuts[0 : min(10, len(cp.cuts))]:
+                print("  {}".format(shorten(str(c), width=90, placeholder="...")))
+            m2 += cp
+            m2.optimize(relax=True)
+            perc_impr = (
+                abs(m2.objective_value - m.objective_value)
+                / max(abs(m2.objective_value), abs(m.objective_value))
+            ) * 100
+
+            if perc_impr > best_impr:
+                best_impr = perc_impr
+                best_cut = ct
+
+            print(
+                "Linear programming relaxation bound now: %g, improvement of %.2f"
+                % (m2.objective_value, perc_impr)
+            )
+        else:
+            continue
+
+    print("Best cut:  {}   improved dual bound: {:.2f} ".format(best_cut, best_impr))

mip/cbc.py

@@ -361,16 +361,27 @@
 
     int Cbc_solveLinearProgram(Cbc_Model *model);
 
+    /*! Type of cutting plane */
     enum CutType {
-      CT_Gomory         = 0,  /*! Gomory cuts obtained from the tableau */
-      CT_MIR            = 1,  /*! Mixed integer rounding cuts */
-      CT_ZeroHalf       = 2,  /*! Zero-half cuts */
-      CT_Clique         = 3,  /*! Clique cuts */
-      CT_KnapsackCover  = 4,  /*! Knapsack cover cuts */
-      CT_LiftAndProject = 5   /*! Lift and project cuts */
+      CT_Probing          =  0,  /*! Cuts generated evaluating the impact of fixing bounds for integer variables */
+      CT_Gomory           =  1,  /*! Gomory cuts obtained from the tableau, implemented by John Forrest  */
+      CT_GMI              =  2,  /*! Gomory cuts obtained from the tableau, implementation from Giacomo Nannicini focusing on safer cuts */
+      CT_LaGomory         =  3,  /*! Additional gomory cuts, simplification of 'A Relax-and-Cut Framework for Gomory's Mixed-Integer Cuts' by Matteo Fischetti & Domenico Salvagnin */
+      CT_RedSplit         =  4,  /*! Reduce and split cuts, implemented by Francois Margot */
+      CT_RedSplitG        =  5,  /*! Reduce and split cuts, implemented by Giacomo Nannicini */
+      CT_FlowCover        =  6,  /*! Flow cover cuts */
+      CT_MIR              =  7,  /*! Mixed-integer rounding cuts */
+      CT_TwoMIR           =  8,  /*! Two-phase Mixed-integer rounding cuts */
+      CT_LaTwoMIR         =  9,  /*! Lagrangean relaxation for two-phase Mixed-integer rounding cuts, as in CT_LaGomory */
+      CT_LiftAndProject   = 10,  /*! Lift and project cuts */
+      CT_ResidualCapacity = 11,  /*! Residual capacity cuts */
+      CT_ZeroHalf         = 12,  /*! Zero-half cuts */
+      CT_Clique           = 13,  /*! Clique cuts */
+      CT_OddWheel         = 14,  /*! Lifted odd-hole inequalities */
+      CT_KnapsackCover    = 15,  /*! Knapsack cover cuts */
     };
 
-    void Cgl_generateCuts( void *osiClpSolver, enum CutType ct, void *osiCuts, int strength );
+    void Cbc_generateCuts( Cbc_Model *cbcModel, enum CutType ct, void *oc, int depth, int pass );
 
     void Cbc_strengthenPacking(Cbc_Model *model);
 
@@ -534,6 +545,12 @@
     } CGNeighbors;
 
     CGNeighbors CG_conflictingNodes(Cbc_Model *model, void *cgraph, size_t node);
+
+    void Cbc_computeFeatures(Cbc_Model *model, double *features);
+
+    int Cbc_nFeatures();
+
+    const char *Cbc_featureName(int i);
     """
     )
 
@@ -584,9 +601,13 @@
 Cbc_setDblParam = cbclib.Cbc_setDblParam
 Cbc_getSolverPtr = cbclib.Cbc_getSolverPtr
 
-Cgl_generateCuts = cbclib.Cgl_generateCuts
+Cbc_generateCuts = cbclib.Cbc_generateCuts
 Cbc_solveLinearProgram = cbclib.Cbc_solveLinearProgram
 
+Cbc_computeFeatures = cbclib.Cbc_computeFeatures
+Cbc_nFeatures = cbclib.Cbc_nFeatures
+Cbc_featureName = cbclib.Cbc_featureName
+
 OsiCuts_new = cbclib.OsiCuts_new
 OsiCuts_addRowCut = cbclib.OsiCuts_addRowCut
 OsiCuts_addGlobalRowCut = cbclib.OsiCuts_addGlobalRowCut
@@ -882,23 +903,45 @@ def var_set_obj(self, var: "Var", value: numbers.Real):
     def generate_cuts(
         self,
         cut_types: Optional[List[CutType]] = None,
+        depth: int = 0,
+        npass: int = 0,
         max_cuts: int = maxsize,
         min_viol: numbers.Real = 1e-4,
     ) -> CutPool:
         cp = CutPool()
         cbc_model = self._model
         osi_solver = Cbc_getSolverPtr(cbc_model)
         osi_cuts = OsiCuts_new()
+        nbin = 0
+        for v in self.model.vars:
+            if v.var_type == BINARY:
+                nbin += 1
 
         if not cut_types:
-            cut_types = [e for e in CutType]
+            cut_types = [
+                CutType.GOMORY,
+                CutType.MIR,
+                CutType.TWO_MIR,
+                CutType.KNAPSACK_COVER,
+                CutType.CLIQUE,
+                CutType.ZERO_HALF,
+                CutType.ODD_WHEEL,
+            ]
         for cut_type in cut_types:
             if self.__verbose >= 1:
                 logger.info(
                     "searching for violated " "{} cuts ... ".format(cut_type.name)
                 )
+
+            if nbin == 0:
+                if cut_type in [CutType.CLIQUE, CutType.ODD_WHEEL]:
+                    continue
+
             nc1 = OsiCuts_sizeRowCuts(osi_cuts)
-            Cgl_generateCuts(osi_solver, int(cut_type.value), osi_cuts, int(1))
+
+            Cbc_generateCuts(
+                self._model, int(cut_type.value), osi_cuts, int(depth), int(npass),
+            )
             nc2 = OsiCuts_sizeRowCuts(osi_cuts)
             if self.__verbose >= 1:
                 logger.info("{} found.\n".format(nc2 - nc1))
@@ -1551,6 +1594,17 @@ def constr_get_pi(self, constr: Constr) -> Optional[numbers.Real]:
     def constr_get_slack(self, constr: Constr) -> Optional[numbers.Real]:
         return self.__slack[constr.idx]
 
+    def feature_values(self) -> List[float]:
+        n = int(Cbc_nFeatures())
+        fv = ffi.new("double[%d]" % n)
+        Cbc_computeFeatures(self._model, fv)
+        return [float(fv[i]) for i in range(n)]
+
+
+def feature_names() -> List[str]:
+    n = int(Cbc_nFeatures())
+    return [ffi.string(Cbc_featureName(i)).decode("utf-8") for i in range(n)]
+
 
 class ModelOsi(Model):
     def __init__(self, osi_ptr):

mip/constants.py

@@ -40,25 +40,62 @@
 
 # cutting planes types
 class CutType(Enum):
-    """ Types of cuts that can be generated"""
 
-    GOMORY = 0
-    """Gomory Mixed Integer cuts [Gomo69]_ ."""
+    """Types of cuts that can be generated. Each cut type is an implementation
+    in the `COIN-OR Cut Generation Library <https://github.com/coin-or/Cgl>`_.
+    For some cut types multiple implementations are available. Sometimes these
+    implementations were designed with different objectives: for the generation
+    of Gomory cutting planes, for example, the GMI cuts are focused on numerical
+    stability, while Forrest's implementation (GOMORY) is more integrated into
+    the CBC code."""
 
-    MIR = 1
+    PROBING = 0
+    """Cuts generated evaluating the impact of fixing bounds for integer
+    variables"""
+
+    GOMORY = 1
+    """Gomory Mixed Integer cuts [Gomo69]_, as implemented by John Forrest."""
+
+    GMI = 2
+    """Gomory Mixed Integer cuts [Gomo69]_, as implemented by Giacomo
+    Nannicini, focusing on numerically safer cuts."""
+
+    RED_SPLIT = 3
+    """Reduce and split cuts [AGY05]_, implemented by Francois Margot."""
+
+    RED_SPLIT_G = 4
+    """Reduce and split cuts [AGY05]_, implemented by Giacomo Nannicini."""
+
+    FLOW_COVER = 5
+    """Lifted Simple Generalized Flow Cover Cut Generator."""
+
+    MIR = 6
     """Mixed-Integer Rounding cuts [Marc01]_."""
 
-    ZERO_HALF = 2
+    TWO_MIR = 7
+    """Two-phase Mixed-integer rounding cuts."""
+
+    LATWO_MIR = 8
+    """Lagrangean relaxation for two-phase Mixed-integer rounding cuts, as in
+    LAGomory"""
+
+    LIFT_AND_PROJECT = 9
+    """Lift-and-project cuts [BCC93]_, implemented by Pierre Bonami."""
+
+    RESIDUAL_CAPACITY = 10
+    """Residual capacity cuts [AtRa02]_, implemented by Francisco Barahona."""
+
+    ZERO_HALF = 11
     """Zero/Half cuts [Capr96]_."""
 
-    CLIQUE = 3
+    CLIQUE = 12
     """Clique cuts [Padb73]_."""
 
-    KNAPSACK_COVER = 4
-    """Knapsack cover cuts [Bala75]_."""
+    ODD_WHEEL = 13
+    """Lifted odd-hole inequalities."""
 
-    LIFT_AND_PROJECT = 5
-    """Lift-and-project cuts [BCC93]_."""
+    KNAPSACK_COVER = 14
+    """Knapsack cover cuts [Bala75]_."""
 
 
 # optimization status

mip/entities.py

@@ -190,7 +190,7 @@ def __str__(self) -> str:
                 result.append(str(abs(coeff)) if abs(coeff) != 1 else "")
                 result.append("{var} ".format(**locals()))
 
-        if hasattr(self, "__sense"):
+        if self.__sense:
             if self.__sense == mip.EQUAL:
                 result.append(" = ")
             if self.__sense == mip.LESS_OR_EQUAL: