Merge pull request #2026 from ERGO-Code/mip-analysis

MIP analysis

FEATURES.md

@@ -9,6 +9,10 @@ The method `Highs::getDualObjectiveValue` now exitsts to compute the dual object
 The method `Highs::getStandardFormLp` now exists to return the incumbent LP in standard form - overlooking any integrality or Hessian. To determine the sizes of the vectors of data, the method is called without specifying pointers to the data arrays.
 
 
+Added documentation on the use of presolve when solving an incumbent model, and clarifying the use of the method `Highs::presolve`.
+
+
+
 
 
 

check/TestLpSolvers.cpp

@@ -635,7 +635,9 @@ TEST_CASE("standard-form-lp", "[highs_lp_solver]") {
   REQUIRE(highs.addRow(1.0, 1.0, 4, index.data(), value.data()) ==
           HighsStatus::kOk);
   REQUIRE(highs.changeObjectiveSense(ObjSense::kMaximize) == HighsStatus::kOk);
-  printf(
-      "\nNow test by adding a fixed column and a fixed row, and maximizing\n");
+  if (dev_run)
+    printf(
+        "\nNow test by adding a fixed column and a fixed row, and "
+        "maximizing\n");
   testStandardForm(highs.getLp());
 }

cmake/sources-python.cmake

@@ -208,6 +208,7 @@ set(highs_sources_python
     src/mip/HighsImplications.cpp
     src/mip/HighsLpAggregator.cpp
     src/mip/HighsLpRelaxation.cpp
+    src/mip/HighsMipAnalysis.cpp
     src/mip/HighsMipSolver.cpp
     src/mip/HighsMipSolverData.cpp
     src/mip/HighsModkSeparator.cpp
@@ -325,6 +326,7 @@ set(highs_headers_python
     src/mip/HighsImplications.h
     src/mip/HighsLpAggregator.h
     src/mip/HighsLpRelaxation.h
+    src/mip/HighsMipAnalysis.h
     src/mip/HighsMipSolver.h
     src/mip/HighsMipSolverData.h
     src/mip/HighsModkSeparator.h
@@ -339,6 +341,7 @@ set(highs_headers_python
     src/mip/HighsSeparator.h
     src/mip/HighsTableauSeparator.h
     src/mip/HighsTransformedLp.h
+    src/mip/MipTimer.h
     src/model/HighsHessian.h
     src/model/HighsHessianUtils.h
     src/model/HighsModel.h

cmake/sources.cmake

@@ -209,6 +209,7 @@ set(highs_sources
     mip/HighsImplications.cpp
     mip/HighsLpAggregator.cpp
     mip/HighsLpRelaxation.cpp
+    mip/HighsMipAnalysis.cpp
     mip/HighsMipSolver.cpp
     mip/HighsMipSolverData.cpp
     mip/HighsModkSeparator.cpp
@@ -329,6 +330,7 @@ set(highs_headers
     mip/HighsImplications.h
     mip/HighsLpAggregator.h
     mip/HighsLpRelaxation.h
+    mip/HighsMipAnalysis.h
     mip/HighsMipSolver.h
     mip/HighsMipSolverData.h
     mip/HighsModkSeparator.h
@@ -343,6 +345,7 @@ set(highs_headers
     mip/HighsSeparator.h
     mip/HighsTableauSeparator.h
     mip/HighsTransformedLp.h
+    mip/MipTimer.h
     model/HighsHessian.h
     model/HighsHessianUtils.h
     model/HighsModel.h

docs/src/guide/further.md

@@ -47,15 +47,15 @@ is [changeRowsBounds](@ref Modify-model-data). An individual matrix coefficient
 is changed by passing its row index, column index and new value to
 [changeCoeff](@ref Modify-model-data).
 
-### [Hot start](@id hot-start)
+## [Hot start](@id hot-start)
 
 It may be possible for HiGHS to start solving a model using data
 obtained by solving a related model, or supplied by a user. Whether
 this is possible depends on the the class of model being solved, the
 solver to be used, and the modifications (if any) that have been to
 the incumbent model since it was last solved.
 
-#### LP
+### LP
 
 To run HiGHS from a user-defined solution or basis, this is passed to HiGHS
 using the methods [setSolution](@ref Set-solution) or [setBasis](@ref Set-basis). The basis passed to HiGHS need not be complete
@@ -71,7 +71,7 @@ using the methods [setSolution](@ref Set-solution) or [setBasis](@ref Set-basis)
 * For nonbasic variables, it is unnecessary to specify whether they
   are at their lower or upper bound unless they are "boxed" variables.
 
-#### MIP
+### MIP
 
 If a (partial) feasible assignment of the integer variables is known,
 this can be passed to HiGHS via [setSolution](@ref Set-solution). If
@@ -81,4 +81,26 @@ integer variables is not complete). If a feasible solution is
 obtained, it will be used to provide the MIP solver with an initial
 primal bound when it run to solve for all integer variables.
 
-
+## [Presolve](@id guide-presolve)
+
+HiGHS has a sophisticated presolve procedure for LPs and MIPs that
+aims to reduce the dimension of the model that must be solved. In most
+cases, the time saved by solving the reduced model is very much
+greater than the time taken to perform presolve. Once he presolved
+model is solved, a postsolve procedure (of minimal computational cost)
+deduces the optimal solution to the original model. Hence presolve is
+performed by default. The only exception occus when there is a valid
+basis for an LP and the simplex solver is used. In this case the
+original LP is solved, starting from this basis. In cases where the
+use of presolve is found not to be advantageous, its use can be
+switched off by setting the [presolve](@ref option-presolve) option to
+"off".
+
+HiGHS has a method [presolve](@ref Presolve/postsolve) that performs presolve on
+the incumbent model, allowing the presolved model to be extracted or
+written to a file. This is intended for users who have their own
+solution technique that they wish to test using models presolved by
+HiGHS. Note that this does not affect how the incumbent model is
+solved. There are two corresponding [postsolve](@ref Presolve/postsolve)
+methods, according to whether there are just solution values, or also
+a basis.

docs/src/interfaces/python/example-py.md

@@ -222,3 +222,20 @@ print('Basis validity = ', h.basisValidityToString(info.basis_validity))
 ## Set basis
 
  * `setBasis`
+
+## Presolve/postsolve
+
+ * `presolve`
+ * `getPresolvedLp`
+ * `getPresolvedModel`
+ * `getPresolveLog`
+ * `getPresolveOrigColsIndex`
+ * `getPresolveOrigRowsIndex`
+ * `getModelPresolveStatus`
+ * `writePresolvedModel`
+ * `presolveStatusToString`
+ * `presolveRuleTypeToString`
+ * `postsolve`
+
+
+

docs/src/options/definitions.md

@@ -1,6 +1,6 @@
 # [List of options](@id option-definitions)
 
-## presolve
+## [presolve](@id option-presolve)
 - Presolve option: "off", "choose" or "on"
 - Type: string
 - Default: "choose"
@@ -236,7 +236,7 @@
 - Default: ""
 
 ## mip\_max\_leaves
-- MIP solver max number of leave nodes
+- MIP solver max number of leaf nodes
 - Type: integer
 - Range: {0, 2147483647}
 - Default: 2147483647

examples/plot_highs_log.py

@@ -101,7 +101,7 @@ def plot_highs_log(
     # Add fourth y-axis for Gap % values (scaled)
     ax4 = ax1.twinx()
     ax4.spines["right"].set_position(("outward", 90))
-    ax4.plot(time_values, gap_values, label="Gap %", color=gap_colour, linestyle="--", linewidth=0.5)
+    ax4.plot(time_values, gap_values, label="Gap %", color=gap_colour, linestyle="--", linewidth=1)
 #    ax4.set_ylabel("Gap.%", color=gap_colour, loc="top", labelpad=22)
     ax4.yaxis.label.set_rotation(0)
     ax4.tick_params(axis="y", labelcolor=gap_colour)

src/Highs.h

@@ -176,7 +176,9 @@ class Highs {
   HighsStatus readBasis(const std::string& filename);
 
   /**
-   * @brief Presolve the incumbent model
+   * @brief Presolve the incumbent model, allowing the presolved model
+   * to be extracted. Subsequent solution of the incumbent model will
+   * only use presolve if there is no valid basis
    */
   HighsStatus presolve();
 

src/ipm/ipx/lp_solver.cc

@@ -354,7 +354,11 @@ void LpSolver::ClearSolution() {
 }
 
 void LpSolver::InteriorPointSolve() {
-  control_.hLog("Interior Point Solve\n");
+    if (control_.runCentring()) {
+      control_.hLog("Interior point solve for analytic centre\n");
+    } else {
+      control_.hLog("Interior point solve\n");
+    }
 
     // Allocate new iterate and set tolerances for IPM termination test.
     iterate_.reset(new Iterate(model_));

src/lp_data/HConst.h

@@ -86,11 +86,14 @@ enum HighsAnalysisLevel {
   kHighsAnalysisLevelSolverTime = 8,
   kHighsAnalysisLevelNlaData = 16,
   kHighsAnalysisLevelNlaTime = 32,
+  kHighsAnalysisLevelMipData = 64,
+  kHighsAnalysisLevelMipTime = 128,
   kHighsAnalysisLevelMin = kHighsAnalysisLevelNone,
   kHighsAnalysisLevelMax =
       kHighsAnalysisLevelModelData + kHighsAnalysisLevelSolverSummaryData +
       kHighsAnalysisLevelSolverRuntimeData + kHighsAnalysisLevelSolverTime +
-      kHighsAnalysisLevelNlaData + kHighsAnalysisLevelNlaTime
+      kHighsAnalysisLevelNlaData + kHighsAnalysisLevelNlaTime +
+      kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime
 };
 
 enum class HighsVarType : uint8_t {

src/lp_data/Highs.cpp

@@ -3399,6 +3399,10 @@ HighsPresolveStatus Highs::runPresolve(const bool force_lp_presolve,
     // Presolved model is extracted now since it's part of solver,
     // which is lost on return
     HighsMipSolver solver(callback_, options_, original_lp, solution_);
+    // Start the MIP solver's total clock so that timeout in presolve
+    // can be identified
+    solver.timer_.start(timer_.total_clock);
+    // Only place that HighsMipSolver::runPresolve is called
     solver.runPresolve(options_.presolve_reduction_limit);
     presolve_return_status = solver.getPresolveStatus();
     // Assign values to data members of presolve_

src/lp_data/HighsModelUtils.cpp

@@ -104,7 +104,11 @@ std::string statusToString(const HighsBasisStatus status, const double lower,
       return "BS";
       break;
     case HighsBasisStatus::kUpper:
-      return "UB";
+      if (lower == upper) {
+        return "FX";
+      } else {
+        return "UB";
+      }
       break;
     case HighsBasisStatus::kZero:
       return "FR";

src/meson.build

@@ -232,6 +232,7 @@ _srcs = [
     'mip/HighsImplications.cpp',
     'mip/HighsLpAggregator.cpp',
     'mip/HighsLpRelaxation.cpp',
+    'mip/HighsMipAnalysis.cpp',
     'mip/HighsMipSolver.cpp',
     'mip/HighsMipSolverData.cpp',
     'mip/HighsModkSeparator.cpp',

src/mip/HighsImplications.cpp

@@ -13,6 +13,7 @@
 #include "../extern/pdqsort/pdqsort.h"
 #include "mip/HighsCliqueTable.h"
 #include "mip/HighsMipSolverData.h"
+#include "mip/MipTimer.h"
 
 bool HighsImplications::computeImplications(HighsInt col, bool val) {
   HighsDomain& globaldomain = mipsolver.mipdata_->domain;
@@ -539,7 +540,10 @@ void HighsImplications::separateImpliedBounds(
           (implicationsCached(col, 0) && implicationsCached(col, 1)))
         continue;
 
-      if (runProbing(col, numboundchgs)) {
+      mipsolver.analysis_.mipTimerStart(kMipClockProbingImplications);
+      const bool probing_result = runProbing(col, numboundchgs);
+      mipsolver.analysis_.mipTimerStop(kMipClockProbingImplications);
+      if (probing_result) {
         if (globaldomain.infeasible()) return;
       }
 

src/mip/HighsLpRelaxation.cpp

@@ -17,6 +17,7 @@
 #include "mip/HighsMipSolver.h"
 #include "mip/HighsMipSolverData.h"
 #include "mip/HighsPseudocost.h"
+#include "mip/MipTimer.h"
 #include "util/HighsCDouble.h"
 #include "util/HighsHash.h"
 
@@ -547,7 +548,9 @@ void HighsLpRelaxation::removeCuts(HighsInt ndelcuts,
     assert(lpsolver.getLp().num_row_ == (HighsInt)lprows.size());
     basis.debug_origin_name = "HighsLpRelaxation::removeCuts";
     lpsolver.setBasis(basis);
+    mipsolver.analysis_.mipTimerStart(kMipClockSimplexBasisSolveLp);
     lpsolver.run();
+    mipsolver.analysis_.mipTimerStop(kMipClockSimplexBasisSolveLp);
   }
 }
 
@@ -1045,7 +1048,22 @@ void HighsLpRelaxation::setObjectiveLimit(double objlim) {
 HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
   lpsolver.setOptionValue(
       "time_limit", lpsolver.getRunTime() + mipsolver.options_mip_->time_limit -
-                        mipsolver.timer_.read(mipsolver.timer_.solve_clock));
+                        mipsolver.timer_.read(mipsolver.timer_.total_clock));
+  // lpsolver.setOptionValue("output_flag", true);
+  const bool valid_basis = lpsolver.getBasis().valid;
+  const HighsInt simplex_solve_clock = valid_basis
+                                           ? kMipClockSimplexBasisSolveLp
+                                           : kMipClockSimplexNoBasisSolveLp;
+  const bool dev_report = false;
+  if (dev_report && !mipsolver.submip) {
+    if (valid_basis) {
+      printf("Solving LP (%7d, %7d) with    a valid basis\n",
+             int(lpsolver.getNumCol()), int(lpsolver.getNumRow()));
+    } else {
+      printf("Solving LP (%7d, %7d) without a valid basis\n",
+             int(lpsolver.getNumCol()), int(lpsolver.getNumRow()));
+    }
+  }
   const bool solver_logging = false;
   const bool detailed_simplex_logging = false;
   if (solver_logging) lpsolver.setOptionValue("output_flag", true);
@@ -1055,7 +1073,10 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
     lpsolver.setOptionValue("highs_analysis_level",
                             kHighsAnalysisLevelSolverRuntimeData);
   }
+
+  mipsolver.analysis_.mipTimerStart(simplex_solve_clock);
   HighsStatus callstatus = lpsolver.run();
+  mipsolver.analysis_.mipTimerStop(simplex_solve_clock);
 
   const HighsInfo& info = lpsolver.getInfo();
   HighsInt itercount = std::max(HighsInt{0}, info.simplex_iteration_count);
@@ -1153,7 +1174,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
                      "HighsLpRelaxation::run LP is unbounded with no basis, "
                      "but not returning Status::kError\n");
       if (info.primal_solution_status == kSolutionStatusFeasible)
-        mipsolver.mipdata_->trySolution(lpsolver.getSolution().col_value, 'T');
+        mipsolver.mipdata_->trySolution(lpsolver.getSolution().col_value,
+                                        kSolutionSourceUnbounded);
 
       return Status::kUnbounded;
     case HighsModelStatus::kUnknown:
@@ -1193,7 +1215,9 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
         // istanbul-no-cutoff
         ipm.setOptionValue("simplex_iteration_limit",
                            info.simplex_iteration_count);
+        mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp);
         ipm.run();
+        mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp);
         lpsolver.setBasis(ipm.getBasis(), "HighsLpRelaxation::run IPM basis");
         return run(false);
       }
@@ -1371,7 +1395,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
             objsum += roundsol[i] * mipsolver.colCost(i);
 
-          mipsolver.mipdata_->addIncumbent(roundsol, double(objsum), 'S');
+          mipsolver.mipdata_->addIncumbent(roundsol, double(objsum),
+                                           kSolutionSourceSolveLp);
           objsum = 0;
         }