Presolve profile

check/TestPresolve.cpp

@@ -910,3 +910,20 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") {
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible);
 }
+
+TEST_CASE("presolve-light", "[highs_test_presolve]") {
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/afiro.mps";
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.readModel(model_file);
+  for (HighsInt k = 0; k < 2; k++) {
+    REQUIRE(highs.presolve() == HighsStatus::kOk);
+    REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kReduced);
+    const HighsLp& presolved_lp = highs.getPresolvedLp();
+    if (dev_run) printf("%s presolved LP has %d columns; %d rows and %d nonzeros\n",
+			k == 0 ? "Fully" : "Lightly", int(presolved_lp.num_col_),
+			int(presolved_lp.num_row_), int(presolved_lp.numNz()));
+    highs.setOptionValue("presolve_light", true);
+  }
+}

cmake/sources.cmake

@@ -558,6 +558,7 @@ set(highs_headers
     presolve/ICrashUtil.h
     presolve/ICrashX.h
     presolve/PresolveComponent.h
+    presolve/PresolveTimer.h
     qpsolver/a_asm.hpp
     qpsolver/a_quass.hpp
     qpsolver/basis.hpp

highs/Highs.h

@@ -717,7 +717,7 @@ class Highs {
    * @brief Get the number of (constraint matrix) nonzeros in the incumbent
    * model
    */
-  HighsInt getNumNz() const { return model_.lp_.a_matrix_.numNz(); }
+  HighsInt getNumNz() const { return model_.lp_.numNz(); }
 
   /**
    * @brief Get the number of Hessian matrix nonzeros in the incumbent model

highs/interfaces/highs_c_api.cpp

@@ -1184,7 +1184,7 @@ HighsInt Highs_getPresolvedNumRow(const void* highs) {
 }
 
 HighsInt Highs_getPresolvedNumNz(const void* highs) {
-  return ((Highs*)highs)->getPresolvedLp().a_matrix_.numNz();
+  return ((Highs*)highs)->getPresolvedLp().numNz();
 }
 
 // Gets pointers to all the public data members of HighsLp: avoids
@@ -1231,7 +1231,7 @@ static HighsInt Highs_getHighsLpData(const HighsLp& lp, const HighsInt a_format,
         (desired_a_format == MatrixFormat::kRowwise &&
          lp.a_matrix_.isRowwise())) {
       // Incumbent format is OK
-      *num_nz = lp.a_matrix_.numNz();
+      *num_nz = lp.numNz();
       if (a_start)
         memcpy(a_start, lp.a_matrix_.start_.data(),
                num_start_entries * sizeof(HighsInt));

highs/io/HighsIO.cpp

@@ -310,6 +310,11 @@ const std::string highsBoolToString(const bool b, const HighsInt field_width) {
   return b ? " true" : "false";
 }
 
+const std::string highsIntToPlural(const HighsInt i, const bool y) {
+  if (y) return i == 1 ? "y" : "ies";
+  return i == 1 ? "" : "s";
+}
+
 const std::string highsTimeToString(const double time) {
   return
 #ifndef NDEBUG

highs/io/HighsIO.h

@@ -111,6 +111,7 @@ std::string highsFormatToString(const char* format, ...);
 const std::string highsBoolToString(const bool b,
                                     const HighsInt field_width = 2);
 
+const std::string highsIntToPlural(const HighsInt i, const bool y = false);
 const std::string highsInsertMdEscapes(const std::string& from_string);
 const std::string highsInsertMdId(const std::string& from_string);
 const std::string highsTimeToString(const double time);

highs/lp_data/HConst.h

@@ -88,12 +88,14 @@ enum HighsAnalysisLevel {
   kHighsAnalysisLevelNlaTime = 32,
   kHighsAnalysisLevelMipData = 64,
   kHighsAnalysisLevelMipTime = 128,
+  kHighsAnalysisLevelPresolveTime = 256,
   kHighsAnalysisLevelMin = kHighsAnalysisLevelNone,
   kHighsAnalysisLevelMax =
       kHighsAnalysisLevelModelData + kHighsAnalysisLevelSolverSummaryData +
       kHighsAnalysisLevelSolverRuntimeData + kHighsAnalysisLevelSolverTime +
       kHighsAnalysisLevelNlaData + kHighsAnalysisLevelNlaTime +
-      kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime
+      kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime +
+      kHighsAnalysisLevelPresolveTime
 };
 
 enum class HighsVarType : uint8_t {
@@ -275,7 +277,9 @@ enum PresolveRuleType : int {
   kPresolveRuleSparsify,
   kPresolveRuleProbing,
   kPresolveRuleEnumeration,
-  kPresolveRuleMax = kPresolveRuleEnumeration,
+  kPresolveRuleDualFixing,
+  kPresolveRuleColStuffing,
+  kPresolveRuleMax = kPresolveRuleColStuffing,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };

highs/lp_data/Highs.cpp

@@ -1375,7 +1375,7 @@ HighsStatus Highs::calledOptimizeModel() {
     time += timer_.read(timer_.solve_clock);
   };
 
-  const bool unconstrained_lp = incumbent_lp.a_matrix_.numNz() == 0;
+  const bool unconstrained_lp = incumbent_lp.numNz() == 0;
   assert(incumbent_lp.num_row_ || unconstrained_lp);
   const bool has_basis = basis_.useful;
   if (has_basis) {
@@ -3676,8 +3676,8 @@ HighsPresolveStatus Highs::runPresolve(const bool force_lp_presolve,
           original_lp.num_col_ - reduced_lp.num_col_;
       presolve_.info_.n_rows_removed =
           original_lp.num_row_ - reduced_lp.num_row_;
-      presolve_.info_.n_nnz_removed = (HighsInt)original_lp.a_matrix_.numNz() -
-                                      (HighsInt)reduced_lp.a_matrix_.numNz();
+      presolve_.info_.n_nnz_removed =
+          (HighsInt)original_lp.numNz() - (HighsInt)reduced_lp.numNz();
       // Clear any scaling information inherited by the reduced LP
       reduced_lp.clearScale();
       assert(lpDimensionsOk("RunPresolve: reduced_lp", reduced_lp,
@@ -3687,7 +3687,7 @@ HighsPresolveStatus Highs::runPresolve(const bool force_lp_presolve,
     case HighsPresolveStatus::kReducedToEmpty: {
       presolve_.info_.n_cols_removed = original_lp.num_col_;
       presolve_.info_.n_rows_removed = original_lp.num_row_;
-      presolve_.info_.n_nnz_removed = (HighsInt)original_lp.a_matrix_.numNz();
+      presolve_.info_.n_nnz_removed = (HighsInt)original_lp.numNz();
       break;
     }
     default:

highs/lp_data/HighsLp.h

@@ -64,6 +64,7 @@ class HighsLp {
   bool equalNames(const HighsLp& lp) const;
   bool equalScaling(const HighsLp& lp) const;
   bool isMip() const;
+  HighsInt numNz() const { return this->a_matrix_.numNz(); }
   bool hasSemiVariables() const;
   bool hasInfiniteCost(const double infinite_cost) const;
   bool hasMods() const;

highs/lp_data/HighsLpUtils.cpp

@@ -76,7 +76,7 @@ HighsStatus assessLp(HighsLp& lp, const HighsOptions& options) {
   // If the LP has no columns the matrix must be empty and there is
   // nothing left to test
   if (lp.num_col_ == 0) {
-    assert(!lp.a_matrix_.numNz());
+    assert(!lp.numNz());
     return HighsStatus::kOk;
   }
   // From here, any LP has lp.num_col_ > 0 and lp.a_matrix_.start_[lp.num_col_]
@@ -92,7 +92,7 @@ HighsStatus assessLp(HighsLp& lp, const HighsOptions& options) {
   if (return_status == HighsStatus::kError) return return_status;
   // If entries have been removed from the matrix, resize the index
   // and value vectors to prevent bug in presolve
-  HighsInt lp_num_nz = lp.a_matrix_.numNz();
+  HighsInt lp_num_nz = lp.numNz();
   if ((HighsInt)lp.a_matrix_.index_.size() > lp_num_nz)
     lp.a_matrix_.index_.resize(lp_num_nz);
   if ((HighsInt)lp.a_matrix_.value_.size() > lp_num_nz)
@@ -3102,7 +3102,7 @@ void reportPresolveReductions(const HighsLogOptions& log_options,
                               const HighsLp& lp, const HighsLp& presolved_lp) {
   const HighsInt num_col_from = lp.num_col_;
   const HighsInt num_row_from = lp.num_row_;
-  const HighsInt num_nz_from = lp.a_matrix_.numNz();
+  const HighsInt num_nz_from = lp.numNz();
   HighsInt num_col_to = 0;
   HighsInt num_row_to = 0;
   HighsInt num_nz_to = 0;
@@ -3124,7 +3124,7 @@ void reportPresolveReductions(const HighsLogOptions& log_options,
     case HighsPresolveStatus::kTimeout: {
       num_col_to = presolved_lp.num_col_;
       num_row_to = presolved_lp.num_row_;
-      num_nz_to = presolved_lp.a_matrix_.numNz();
+      num_nz_to = presolved_lp.numNz();
       message =
           presolve_status == HighsPresolveStatus::kTimeout ? "- Timeout" : "";
       break;

highs/lp_data/HighsModelUtils.cpp

@@ -598,10 +598,10 @@ void writeGlpsolSolution(FILE* file, const HighsOptions& options,
   assert(lp.row_names_.size() == static_cast<size_t>(lp.num_row_));
   // Determine number of nonzeros including the objective function
   // and, hence, determine whether there is an objective function
-  HighsInt num_nz = lp.a_matrix_.numNz();
+  HighsInt num_nz = lp.numNz();
   for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
     if (lp.col_cost_[iCol]) num_nz++;
-  const bool empty_cost_row = num_nz == lp.a_matrix_.numNz();
+  const bool empty_cost_row = num_nz == lp.numNz();
   const bool has_objective = !empty_cost_row || model.hessian_.dim_;
   // Writes the solution using the GLPK raw style (defined in
   // api/wrsol.c) or pretty style (defined in api/prsol.c)
@@ -702,7 +702,7 @@ void writeGlpsolSolution(FILE* file, const HighsOptions& options,
   const HighsInt glpsol_num_row = num_row + delta_num_row;
   // If the cost row isn't reported, then the number of nonzeros is
   // just the number in the constraint matrix
-  if (cost_row_location <= 0) num_nz = lp.a_matrix_.numNz();
+  if (cost_row_location <= 0) num_nz = lp.numNz();
   // Record the discrete nature of the model
   HighsInt num_integer = 0;
   HighsInt num_binary = 0;
@@ -1515,6 +1515,10 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) {
     return "Probing";
   } else if (rule_type == kPresolveRuleEnumeration) {
     return "Enumeration";
+  } else if (rule_type == kPresolveRuleDualFixing) {
+    return "Dual fixing";
+  } else if (rule_type == kPresolveRuleColStuffing) {
+    return "Col stuffing";
   }
   assert(1 == 0);
   return "????";

highs/lp_data/HighsOptions.h

@@ -421,6 +421,7 @@ struct HighsOptionsStruct {
   bool lp_presolve_requires_basis_postsolve;
   bool mps_parser_type_free;
   bool use_warm_start;
+  bool presolve_light;
   HighsInt keep_n_rows;
   HighsInt cost_scale_factor;
   HighsInt allowed_matrix_scale_factor;
@@ -588,6 +589,7 @@ struct HighsOptionsStruct {
         lp_presolve_requires_basis_postsolve(false),
         mps_parser_type_free(false),
         use_warm_start(true),
+        presolve_light(false),
         keep_n_rows(0),
         cost_scale_factor(0),
         allowed_matrix_scale_factor(0),
@@ -1452,6 +1454,11 @@ class HighsOptions : public HighsOptionsStruct {
                                        advanced, &use_warm_start, true);
     records.push_back(record_bool);
 
+    record_bool = new OptionRecordBool("presolve_light",
+                                       "Use only low-cost presolve rules",
+                                       advanced, &presolve_light, false);
+    records.push_back(record_bool);
+
     record_int =
         new OptionRecordInt("keep_n_rows",
                             "For multiple N-rows in MPS files: delete rows / "

highs/lp_data/HighsSolve.cpp

@@ -59,7 +59,7 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) {
     }
     return return_status;
   };
-  if (!solver_object.lp_.num_row_ || solver_object.lp_.a_matrix_.numNz() == 0) {
+  if (!solver_object.lp_.num_row_ || solver_object.lp_.numNz() == 0) {
     // LP is unconstrained due to having no rows or a zero constraint
     // matrix, so solve directly
     call_status = solveUnconstrainedLp(solver_object);
@@ -208,11 +208,11 @@ HighsStatus solveUnconstrainedLp(const HighsOptions& options, const HighsLp& lp,
   resetModelStatusAndHighsInfo(model_status, highs_info);
 
   // Check that the LP really is unconstrained!
-  assert(lp.num_row_ == 0 || lp.a_matrix_.numNz() == 0);
+  assert(lp.num_row_ == 0 || lp.numNz() == 0);
   if (lp.num_row_ > 0) {
     // LP has rows, but should only be here if the constraint matrix
     // is zero
-    if (lp.a_matrix_.numNz() > 0) return HighsStatus::kError;
+    if (lp.numNz() > 0) return HighsStatus::kError;
   }
 
   highsLogUser(options.log_options, HighsLogType::kInfo,
@@ -465,7 +465,7 @@ void assessExcessiveObjectiveBoundScaling(const HighsLogOptions log_options,
 
   double min_matrix_value = kHighsInf;
   double max_matrix_value = -kHighsInf;
-  const HighsInt num_matrix_nz = lp.a_matrix_.numNz();
+  const HighsInt num_matrix_nz = lp.numNz();
   for (HighsInt iEl = 0; iEl < num_matrix_nz; iEl++)
     assessFiniteNonzero(lp.a_matrix_.value_[iEl], min_matrix_value,
                         max_matrix_value);

highs/mip/MipTimer.h

@@ -11,7 +11,7 @@
 #ifndef MIP_MIPTIMER_H_
 #define MIP_MIPTIMER_H_
 
-// Clocks for profiling the MIP dual mip solver
+// Clocks for profiling the MIP solver
 enum iClockMip {
   kMipClockTotal = 0,
   kMipClockPresolve,
@@ -114,7 +114,7 @@ enum iClockMip {
   kNumMipClock  //!< Number of MIP clocks
 };
 
-const double tolerance_percent_report = 0.1;
+static const double kMipClockTolerancePercentReport = 0.1;
 
 class MipTimer {
  public:
@@ -371,7 +371,7 @@ class MipTimer {
                                                kMipClockSearch,
                                                kMipClockPostsolve};
     reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal, tolerance_percent_report);
+                       kMipClockTotal, kMipClockTolerancePercentReport);
   };
 
   void reportMipSolveLpClock(const HighsTimerClock& mip_timer_clock) {
@@ -385,20 +385,20 @@ class MipTimer {
         kMipClockHipoSolveLp,
         kMipClockIpxSolveLp};
     reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSubMipSolveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockSubMipSolve};
     reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipPresolveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockProbingPresolve,
                                                kMipClockEnumerationPresolve};
     reportMipClockList("MipPrslv", mip_clock_list, mip_timer_clock,
-                       kMipClockRunPresolve, tolerance_percent_report);
+                       kMipClockRunPresolve, kMipClockTolerancePercentReport);
   };
 
   void reportAltEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) {
@@ -407,7 +407,7 @@ class MipTimer {
                                                kMipClockEvaluateRootNode2};
     reportMipClockList(
         "AltEvaluateRootNode", mip_clock_list, mip_timer_clock,
-        kMipClockEvaluateRootNode);  //, tolerance_percent_report);
+        kMipClockEvaluateRootNode);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) {
@@ -433,7 +433,7 @@ class MipTimer {
     };
     reportMipClockList(
         "MipEvaluateRootNode", mip_clock_list, mip_timer_clock,
-        kMipClockEvaluateRootNode);  //, tolerance_percent_report);
+        kMipClockEvaluateRootNode);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipRootSeparationClock(const HighsTimerClock& mip_timer_clock) {
@@ -442,8 +442,9 @@ class MipTimer {
         kMipClockRootSeparationFinishAnalyticCentreComputation,
         kMipClockRootSeparationCentralRounding,
         kMipClockRootSeparationEvaluateRootLp};
-    reportMipClockList("MipRootSeparation", mip_clock_list, mip_timer_clock,
-                       kMipClockRootSeparation);  //, tolerance_percent_report);
+    reportMipClockList(
+        "MipRootSeparation", mip_clock_list, mip_timer_clock,
+        kMipClockRootSeparation);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSearchClock(const HighsTimerClock& mip_timer_clock) {
@@ -455,15 +456,15 @@ class MipTimer {
         //	kMipClock@
     };
     reportMipClockList("MipSerch", mip_clock_list, mip_timer_clock,
-                       kMipClockSearch, tolerance_percent_report);
+                       kMipClockSearch, kMipClockTolerancePercentReport);
   };
 
   void reportMipDiveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{
         kMipClockDiveEvaluateNode, kMipClockDivePrimalHeuristics,
         kMipClockTheDive, kMipClockBacktrackPlunge, kMipClockPerformAging2};
     reportMipClockList("MipDive_", mip_clock_list, mip_timer_clock,
-                       kMipClockDive, tolerance_percent_report);
+                       kMipClockDive, kMipClockTolerancePercentReport);
   };
 
   void reportMipDivePrimalHeuristicsClock(
@@ -472,7 +473,7 @@ class MipTimer {
         kMipClockDiveRandomizedRounding, kMipClockDiveRens, kMipClockDiveRins};
     reportMipClockList("MipDivePrimalHeuristics", mip_clock_list,
                        mip_timer_clock, kMipClockDivePrimalHeuristics,
-                       tolerance_percent_report);
+                       kMipClockTolerancePercentReport);
   };
 
   void reportMipNodeSearchClock(const HighsTimerClock& mip_timer_clock) {
@@ -481,16 +482,17 @@ class MipTimer {
         //      kMipClockSearchBacktrack,
         kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1,
         kMipClockNodeSearchSeparation};  //, kMipClockStoreBasis};
-    reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock,
-                       kMipClockNodeSearch);  //, tolerance_percent_report);
+    reportMipClockList(
+        "MipNodeSearch", mip_clock_list, mip_timer_clock,
+        kMipClockNodeSearch);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{
         kMipClockImplboundSepa, kMipClockCliqueSepa, kMipClockTableauSepa,
         kMipClockPathAggrSepa, kMipClockModKSepa};
     reportMipClockList("MipSeparation", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void csvMipClock(const std::string model_name,