Merge pull request #3167 from ERGO-Code/fix-3040

Introduces presolve "light"

Author: jajhall

check/TestPresolve.cpp

@@ -909,6 +909,158 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") {
   highs.readModel(model_file);
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible);
+
+  highs.resetGlobalScheduler(true);
+}
+
+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);
+  HighsInt presolved_num_col;
+  HighsInt presolved_num_row;
+  HighsInt presolved_num_nz;
+  for (HighsInt k = 0; k < 2; k++) {
+    REQUIRE(highs.presolve() == HighsStatus::kOk);
+    REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kReduced);
+    const HighsLp& presolved_lp = highs.getPresolvedLp();
+    if (k == 1) {
+      REQUIRE(presolved_lp.num_col_ > presolved_num_col);
+      REQUIRE(presolved_lp.num_row_ > presolved_num_row);
+      REQUIRE(presolved_lp.numNz() > presolved_num_nz);
+    }
+    presolved_num_col = presolved_lp.num_col_;
+    presolved_num_row = presolved_lp.num_row_;
+    presolved_num_nz = presolved_lp.numNz();
+
+    if (dev_run)
+      printf("%s presolved LP has %d columns; %d rows and %d nonzeros\n",
+             k == 0 ? "Fully" : "Lightly", int(presolved_num_col),
+             int(presolved_num_row), int(presolved_num_nz));
+    highs.setOptionValue("presolve_light", kHighsOnString);
+  }
+
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("presolve-initial-sweep-postsolve", "[highs_test_presolve]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  if (dev_run) {
+    highs.setOptionValue("log_dev_level", 1);
+    highs.setOptionValue("presolve_rule_logging", true);
+  }
+  HighsLp lp;
+  lp.num_col_ = 3;
+  lp.num_row_ = 1;
+  lp.col_cost_ = {0, 1, 2};
+  lp.col_lower_ = {1, 2, 0};
+  lp.col_upper_ = {1, kHighsInf, kHighsInf};
+  lp.row_lower_ = {-kHighsInf};
+  lp.row_upper_ = {5};
+  lp.a_matrix_.start_ = {0, 1, 2, 3};
+  lp.a_matrix_.index_ = {0, 0, 0};
+  lp.a_matrix_.value_ = {3, 2, -1};
+  highs.passModel(lp);
+  // Presolved to empty, so no simplex iterations if postsolve is correct
+  highs.run();
+  REQUIRE(highs.getInfo().simplex_iteration_count == 0);
+  if (dev_run) highs.writeSolution("", 1);
+
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("presolve-initial-sweep", "[highs_test_presolve]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  if (dev_run) {
+    highs.setOptionValue("log_dev_level", 1);
+    highs.setOptionValue("presolve_rule_logging", true);
+  }
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 1;
+  lp.col_cost_ = {-1, 1};
+  lp.col_lower_ = {1, 0};
+  lp.col_upper_ = {1, kHighsInf};
+  lp.row_lower_ = {-kHighsInf};
+  lp.row_upper_ = {5};
+  lp.a_matrix_.start_ = {0, 1, 1};
+  lp.a_matrix_.index_ = {0};
+  lp.a_matrix_.value_ = {3};
+  HighsStatus pass_model_status = HighsStatus::kOk;
+  for (HighsInt k = 0; k < 4; k++) {
+    if (dev_run) printf("\nPass k = %d\n==========\n", int(k));
+    REQUIRE(highs.passModel(lp) == pass_model_status);
+    highs.run();
+    if (k == 0) {
+      // Remove the empty column
+      lp.num_col_ = 1;
+      lp.col_cost_ = {-1};
+      lp.col_lower_ = {1};
+      lp.col_upper_ = {1};
+      lp.a_matrix_.start_ = {0, 1};
+    } else if (k == 1) {
+      REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+      lp.col_upper_ = {0};
+      pass_model_status = HighsStatus::kWarning;
+      // Can infeasible column bounds even reach presolve?
+    } else if (k == 2) {
+      REQUIRE(highs.getModelStatus() == HighsModelStatus::kInfeasible);
+      lp.col_lower_ = {1e+20};
+      lp.col_upper_ = {kHighsInf};
+      lp.row_upper_ = {kHighsInf};
+      pass_model_status = HighsStatus::kError;
+      // Why is this model accepted when error returns?
+    } else {
+      REQUIRE(highs.getModelStatus() == HighsModelStatus::kUnbounded);
+    }
+    if (dev_run) highs.writeSolution("", 1);
+  }
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("presolve-initial-sweep-all", "[highs_test_presolve]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("presolve_light", kHighsOnString);
+  if (dev_run) {
+    highs.setOptionValue("log_dev_level", 1);
+    highs.setOptionValue("presolve_rule_logging", true);
+  }
+  HighsLp lp;
+  lp.num_col_ = 7;
+  lp.num_row_ = 4;
+  lp.col_cost_ = {1, 1, 1, 1, 1, 1, 1};
+  lp.col_lower_ = {0, 1, 0, 1, 1, -kHighsInf, 0};
+  lp.col_upper_ = {1, 1, kHighsInf, 3, 1, 1, 1};
+  lp.row_lower_ = {2, 8, 10, 13};
+  lp.row_upper_ = {4, 9, 16, 26};
+  lp.a_matrix_.start_ = {0, 1, 5, 6, 6, 10, 12, 13};
+  lp.a_matrix_.index_ = {2, 0, 1, 2, 3, 0, 0, 1, 2, 3, 2, 3, 2};
+  lp.a_matrix_.value_ = {6, 1, 4, 7, 11, 2, 3, 5, 8, 12, 9, 13, 10};
+  // Cols 1 and 4 fixed at 1; col 3 empty (fixed at LB = 1) then
+  //
+  // Rows 0 and 3 singletons; row 1 empty
+  REQUIRE(highs.passModel(lp) == HighsStatus::kOk);
+
+  highs.run();
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+  if (dev_run) highs.writeSolution("", 1);
+
+  // Add a redundant row
+  std::vector<HighsInt> index = {0, 5, 6};
+  std::vector<double> value = {1, 1, 1};
+  highs.addRow(-kHighsInf, 4, 3, index.data(), value.data());
+
+  highs.setOptionValue("use_warm_start", false);
+  highs.run();
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+  if (dev_run) highs.writeSolution("", 1);
+
+  highs.resetGlobalScheduler(true);
 }
 
 TEST_CASE("bound_implied", "[highs_test_presolve]") {

cmake/sources.cmake

@@ -570,6 +570,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

@@ -781,7 +781,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

@@ -1189,7 +1189,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
@@ -1236,7 +1236,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/FilereaderLp.cpp

@@ -251,8 +251,8 @@ FilereaderRetcode FilereaderLp::readModelFromFile(const HighsOptions& options,
                            "coefficient%s in row %d (name \"%s\")\n",
                            int(iCol), lp.col_names_[iCol].c_str(),
                            int(zero_count[iRow]),
-                           zero_count[iRow] > 1 ? "s" : "", int(iRow),
-                           lp.row_names_[iRow].c_str());
+                           highsIntToPlural(zero_count[iRow]).c_str(),
+                           int(iRow), lp.row_names_[iRow].c_str());
             num_report++;
           }
           sum_num_duplicate += (num_ocurrence - 1);

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/Highs.cpp

@@ -1587,7 +1587,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) {
@@ -3968,8 +3968,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,
@@ -3979,7 +3979,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:
@@ -4000,9 +4000,10 @@ HighsPostsolveStatus Highs::runPostsolve() {
     return HighsPostsolveStatus::kNoPrimalSolutionError;
   const bool have_dual_solution =
       presolve_.data_.recovered_solution_.dual_valid;
-  presolve_.data_.postSolveStack.undo(options_,
-                                      presolve_.data_.recovered_solution_,
-                                      presolve_.data_.recovered_basis_);
+  const HighsInt report_3040_col = -21792;
+  presolve_.data_.postSolveStack.undo(
+      options_, presolve_.data_.recovered_solution_,
+      presolve_.data_.recovered_basis_, report_3040_col);
   // Compute the row activities
   assert(model_.lp_.a_matrix_.isColwise());
   calculateRowValuesQuad(model_.lp_, presolve_.data_.recovered_solution_);

highs/lp_data/HighsInterface.cpp

@@ -66,23 +66,20 @@ void Highs::reportModelStats() const {
   const HighsInt a_num_nz = lp.a_matrix_.numNz();
   const HighsInt q_num_nz = hessian.dim_ > 0 ? hessian.numNz() : 0;
   if (*log_options.log_dev_level) {
-    highsLogDev(log_options, HighsLogType::kInfo, "%4s      : %s\n",
+    highsLogDev(log_options, HighsLogType::kInfo, "%-4s      : %s\n",
                 problem_type.c_str(), lp.model_name_.c_str());
     highsLogDev(log_options, HighsLogType::kInfo,
-                "Row%s      : %" HIGHSINT_FORMAT "\n",
-                lp.num_row_ == 1 ? "" : "s", lp.num_row_);
+                "Rows      : %" HIGHSINT_FORMAT "\n", lp.num_row_);
     highsLogDev(log_options, HighsLogType::kInfo,
-                "Col%s      : %" HIGHSINT_FORMAT "\n",
-                lp.num_col_ == 1 ? "" : "s", lp.num_col_);
+                "Cols      : %" HIGHSINT_FORMAT "\n", lp.num_col_);
     if (q_num_nz) {
       highsLogDev(log_options, HighsLogType::kInfo,
                   "Matrix Nz : %" HIGHSINT_FORMAT "\n", a_num_nz);
       highsLogDev(log_options, HighsLogType::kInfo,
                   "Hessian Nz: %" HIGHSINT_FORMAT "\n", q_num_nz);
     } else {
       highsLogDev(log_options, HighsLogType::kInfo,
-                  "Nonzero%s  : %" HIGHSINT_FORMAT "\n",
-                  a_num_nz == 1 ? "" : "s", a_num_nz);
+                  "Nonzeros  : %" HIGHSINT_FORMAT "\n", a_num_nz);
     }
     if (num_integer)
       highsLogDev(log_options, HighsLogType::kInfo,
@@ -100,26 +97,28 @@ void Highs::reportModelStats() const {
     stats_line << problem_type;
     if (lp.model_name_.length()) stats_line << " " << lp.model_name_;
     stats_line << " has " << lp.num_row_ << " row"
-               << (lp.num_row_ == 1 ? "" : "s") << "; " << lp.num_col_ << " col"
-               << (lp.num_col_ == 1 ? "" : "s");
+               << highsIntToPlural(lp.num_row_) << "; " << lp.num_col_ << " col"
+               << highsIntToPlural(lp.num_col_);
     if (q_num_nz) {
       stats_line << "; " << a_num_nz << " matrix nonzero"
-                 << (a_num_nz == 1 ? "" : "s");
+                 << highsIntToPlural(a_num_nz);
       stats_line << "; " << q_num_nz << " Hessian nonzero"
-                 << (q_num_nz == 1 ? "" : "s");
+                 << highsIntToPlural(q_num_nz);
     } else {
       stats_line << "; " << a_num_nz << " nonzero"
-                 << (a_num_nz == 1 ? "" : "s");
+                 << highsIntToPlural(a_num_nz);
     }
     if (hessian.isOracle()) stats_line << "; Hessian as oracle";
     if (num_integer)
       stats_line << "; " << num_integer << " integer variable"
-                 << (a_num_nz == 1 ? "" : "s") << " (" << num_binary
+                 << highsIntToPlural(num_integer) << " (" << num_binary
                  << " binary)";
     if (num_semi_continuous)
-      stats_line << "; " << num_semi_continuous << " semi-continuous variables";
+      stats_line << "; " << num_semi_continuous << " semi-continuous variable"
+                 << highsIntToPlural(num_semi_continuous);
     if (num_semi_integer)
-      stats_line << "; " << num_semi_integer << " semi-integer variables";
+      stats_line << "; " << num_semi_integer << " semi-integer variable"
+                 << highsIntToPlural(num_semi_integer);
     highsLogUser(log_options, HighsLogType::kInfo, "%s\n",
                  stats_line.str().c_str());
   }

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;