Handle appended rows in presolve (parallel rows reduction)

check/TestPresolveRules.cpp

@@ -3,13 +3,22 @@
 #include "HCheckConfig.h"
 #include "Highs.h"
 #include "catch.hpp"
+#include "presolve/HPresolve.h"
+#include "presolve/HighsPostsolveStack.h"
 
 const bool dev_run = false;
 
+void solveAndCheck(const std::string& message, const HighsLp& lp, Highs& h,
+                   const std::string& solver, bool use_presolve,
+                   const HighsInt require_presolved_model_num_col = -1,
+                   const HighsInt require_presolved_model_num_row = -1,
+                   const HighsInt require_presolved_model_num_nz = -1);
+
 void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
-                   const HighsInt require_presolved_model_num_col = 0,
-                   const HighsInt require_presolved_model_num_row = 0,
-                   const HighsInt require_presolved_model_num_nz = 0);
+                   const std::vector<std::string>& solvers,
+                   const HighsInt require_presolved_model_num_col = -1,
+                   const HighsInt require_presolved_model_num_row = -1,
+                   const HighsInt require_presolved_model_num_nz = -1);
 
 TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   HighsLp lp;
@@ -59,89 +68,134 @@ TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   lp.a_matrix_.start_ = {0, lp.num_col_};
   lp.a_matrix_.index_.resize(lp.num_col_);
   std::iota(lp.a_matrix_.index_.begin(), lp.a_matrix_.index_.end(), 0);
+  const std::vector<std::string> solvers = {kSimplexString, kIpmString,
+                                            kHiPdlpString};
   if (lp1) {
     lp.col_cost_.assign(lp.num_col_, 1);
     lp.a_matrix_.value_.assign(lp.num_col_, 1);
-    presolveOffOn("Capturing neos-787933 issue", lp, h);
+    presolveOffOn("Capturing neos-787933 issue", lp, h, solvers);
   }
   if (lp1a) {
     lp.col_cost_ = {2, 1};
     lp.a_matrix_.value_.assign(lp.num_col_, 1);
-    presolveOffOn("Variant A neos-787933 issue", lp, h);
+    presolveOffOn("Variant A neos-787933 issue", lp, h, solvers);
   }
   if (lp1b) {
     lp.col_cost_ = {-2, -1};
     lp.a_matrix_.value_.assign(lp.num_col_, 1);
-    presolveOffOn("Variant B neos-787933 issue", lp, h);
+    presolveOffOn("Variant B neos-787933 issue", lp, h, solvers);
   }
   lp.clear();
 
   h.resetGlobalScheduler(true);
 }
 
-void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
+TEST_CASE("test-parallel-rows-cut-ordering", "[highs_test_presolve_rules]") {
+  // Rows 0 and 1 are parallel (both [1, 1]). Row 0 is marked as a
+  // cut. detectParallelRowsAndCols must remove the cut row (0) and
+  // keep the non-cut row (1), not the other way around.
+  // Row 2 involves only col 0, breaking column parallelism.
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 3;
+  lp.sense_ = ObjSense::kMinimize;
+  lp.col_cost_ = {1, 2};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {10, 10};
+  lp.row_lower_ = {-kHighsInf, -kHighsInf, -kHighsInf};
+  lp.row_upper_ = {5, 5, 3};
+  lp.a_matrix_.num_col_ = lp.num_col_;
+  lp.a_matrix_.num_row_ = lp.num_row_;
+  lp.a_matrix_.format_ = MatrixFormat::kColwise;
+  lp.a_matrix_.start_ = {0, 3, 5};
+  lp.a_matrix_.index_ = {0, 1, 2, 0, 1};
+  lp.a_matrix_.value_ = {1, 1, 1, 1, 1};
+
+  HighsOptions options;
+  options.presolve_rule_test = kPresolveRuleParallelRowsAndCols;
+
+  HighsTimer timer;
+  timer.start();
+
+  presolve::HighsPostsolveStack postsolve_stack;
+  postsolve_stack.initializeIndexMaps(lp.num_row_, lp.num_col_);
+  // Mark parallel row 0 as a cut
+  postsolve_stack.setRowType(0,
+                             presolve::HighsPostsolveStack::OrigRowType::kCut);
+
+  presolve::HPresolve presolve;
+  presolve.setInput(lp, options, -1, &timer);
+  REQUIRE(presolve.okSetupPresolveDataStructures());
+  HighsModelStatus status = presolve.run(postsolve_stack);
+  timer.stop();
+  REQUIRE(status == HighsModelStatus::kNotset);
+  // One row must have been removed
+  REQUIRE(lp.num_row_ == 1);
+  // The surviving row must be original row 1 (non-cut), not row 0 (cut)
+  REQUIRE(postsolve_stack.getOrigRowIndex(0) == 1);
+  REQUIRE(!postsolve_stack.isCutRow(0));
+}
+
+void solveAndCheck(const std::string& message, const HighsLp& lp, Highs& h,
+                   const std::string& solver, bool use_presolve,
                    const HighsInt require_presolved_model_num_col,
                    const HighsInt require_presolved_model_num_row,
                    const HighsInt require_presolved_model_num_nz) {
   const HighsRunData& run_data = h.getRunData();
-  bool presolve_on = false;
-  // If the model reduces to empty, then the output from different
-  // solvers cannot be tested
-  const bool reduce_to_empty = require_presolved_model_num_col == 0 &&
-                               require_presolved_model_num_row == 0;
-  const HighsInt to_k = reduce_to_empty ? 2 : 5;
-  for (int k = 0; k < to_k; k++) {
-    std::string solver = kSimplexString;
-    std::string run_crossover = kHighsOnString;
-    bool basis_postsolve = true;
-    if (k == 0) {
-      // Presolve off - to get the optimal solution to debug
-      // presolve
-      presolve_on = false;
-    } else {
-      presolve_on = true;
-      if (k == 1) {
-        solver = kSimplexString;
-      } else if (k == 2) {
-        solver = kIpmString;
-      } else if (k == 3) {
-        solver = kIpmString;
-        run_crossover = kHighsOffString;
-        basis_postsolve = false;
-      } else {
-        solver = kHiPdlpString;
-        basis_postsolve = false;
-      }
-    }
-    std::string presolve = presolve_on ? kHighsOnString : kHighsOffString;
-    h.setOptionValue(kPresolveString, presolve);
-    h.setOptionValue(kRunCrossoverString, run_crossover);
-    h.setOptionValue(kSolverString, solver);
-    if (dev_run)
-      printf(
-          "\n============\n%s: presolve = %s; solver = %s%s\n============\n\n",
-          message.c_str(), presolve.c_str(), solver.c_str(),
-          solver == kIpmString ? ("; run_crossover = " + run_crossover).c_str()
-                               : "");
-    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
-    h.run();
-    if (dev_run) h.writeSolution("", 1);
-    if (presolve_on) {
-      // Ensure that the model is reduced to empty
+  std::string run_crossover = kHighsOnString;
+  bool basis_postsolve = true;
+  if (solver == kIpmString) {
+    run_crossover = kHighsOffString;
+    basis_postsolve = false;
+  } else if (solver == kHiPdlpString) {
+    basis_postsolve = false;
+  }
+  std::string presolve = use_presolve ? kHighsOnString : kHighsOffString;
+  h.setOptionValue(kPresolveString, presolve);
+  h.setOptionValue(kRunCrossoverString, run_crossover);
+  h.setOptionValue(kSolverString, solver);
+  if (dev_run)
+    printf("\n============\n%s: presolve = %s; solver = %s%s\n============\n\n",
+           message.c_str(), presolve.c_str(), solver.c_str(),
+           solver == kIpmString ? ("; run_crossover = " + run_crossover).c_str()
+                                : "");
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  h.run();
+  if (dev_run) h.writeSolution("", 1);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().num_primal_infeasibilities == 0);
+  REQUIRE(h.getInfo().num_dual_infeasibilities == 0);
+  if (use_presolve) {
+    // Ensure that the model is reduced as expected
+    if (require_presolved_model_num_col >= 0)
       REQUIRE(run_data.presolved_model_num_col ==
               require_presolved_model_num_col);
+    if (require_presolved_model_num_row >= 0)
       REQUIRE(run_data.presolved_model_num_row ==
               require_presolved_model_num_row);
+    if (require_presolved_model_num_nz >= 0)
       REQUIRE(run_data.presolved_model_num_nz ==
               require_presolved_model_num_nz);
-      // Ensure that dual postsolve is correct
-      REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
-      REQUIRE(h.getInfo().num_primal_infeasibilities == 0);
-      REQUIRE(h.getInfo().num_dual_infeasibilities == 0);
+    if (require_presolved_model_num_col == 0 &&
+        require_presolved_model_num_row == 0)
       REQUIRE(h.getInfo().simplex_iteration_count == 0);
-      // Ensure that any basis postsolve is correct
-      if (basis_postsolve)
-        REQUIRE(run_data.num_simplex_iterations_after_postsolve == 0);
-    }
+    // Ensure that any basis postsolve is correct
+    if (basis_postsolve)
+      REQUIRE(run_data.num_simplex_iterations_after_postsolve == 0);
+  }
+}
+
+void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
+                   const std::vector<std::string>& solvers,
+                   const HighsInt require_presolved_model_num_col,
+                   const HighsInt require_presolved_model_num_row,
+                   const HighsInt require_presolved_model_num_nz) {
+  // Presolve off - to get the optimal solution to debug presolve
+  solveAndCheck(message, lp, h, kSimplexString, false);
+  // Presolve on with each solver
+  for (const std::string& solver : solvers) {
+    solveAndCheck(message, lp, h, solver, true, require_presolved_model_num_col,
+                  require_presolved_model_num_row,
+                  require_presolved_model_num_nz);
   }
 }

highs/mip/HighsCliqueTable.cpp

@@ -1290,15 +1290,10 @@ void HighsCliqueTable::extractCliques(HighsMipSolver& mipsolver,
 
   HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
 
-  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
+  for (HighsInt i : mipsolver.mipdata_->postSolveStack.getNonCutRows()) {
     HighsInt start = mipsolver.mipdata_->ARstart_[i];
     HighsInt end = mipsolver.mipdata_->ARstart_[i + 1];
 
-    if (mipsolver.mipdata_->postSolveStack.isCutRow(i)) {
-      if (!mipsolver.mipdata_->postSolveStack.hasAppendedRows()) break;
-      continue;
-    }
-
     // catch set packing and partitioning constraints that already have the form
     // of a clique without transformations and add those cliques with the rows
     // being recorded

highs/mip/HighsMipSolverData.cpp

@@ -1493,19 +1493,13 @@ void HighsMipSolverData::basisTransfer() {
     firstrootbasis.alien = true;
     firstrootbasis.useful = true;
 
-    for (HighsInt i = 0;
-         i < static_cast<HighsInt>(postSolveStack.getOrigRowIndex().size());
-         ++i) {
-      if (!postSolveStack.isOrigRow(i)) break;
+    for (HighsInt i : postSolveStack.getOrigRows()) {
       HighsBasisStatus status =
           mipsolver.rootbasis->row_status[postSolveStack.getOrigRowIndex(i)];
       firstrootbasis.row_status[i] = status;
     }
 
-    for (HighsInt i = 0;
-         i < static_cast<HighsInt>(postSolveStack.getOrigColIndex().size());
-         ++i) {
-      if (!postSolveStack.isOrigCol(i)) break;
+    for (HighsInt i : postSolveStack.getOrigCols()) {
       HighsBasisStatus status =
           mipsolver.rootbasis->col_status[postSolveStack.getOrigColIndex(i)];
       firstrootbasis.col_status[i] = status;

highs/mip/HighsPseudocost.cpp

@@ -46,8 +46,7 @@ HighsPseudocost::HighsPseudocost(const HighsMipSolver& mipsolver)
     conflict_avg_score =
         mipsolver.pscostinit->conflict_avg_score * mipsolver.numCol();
 
-    for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
-      if (!mipsolver.mipdata_->postSolveStack.isOrigCol(i)) continue;
+    for (HighsInt i : mipsolver.mipdata_->postSolveStack.getOrigCols()) {
       HighsInt origCol = mipsolver.mipdata_->postSolveStack.getOrigColIndex(i);
 
       pseudocostup[i] = mipsolver.pscostinit->pseudocostup[origCol];
@@ -115,8 +114,7 @@ HighsPseudocostInitialization::HighsPseudocostInitialization(
   HighsInt ncols = pscost.pseudocostup.size();
   conflict_avg_score /= ncols * pscost.conflict_weight;
 
-  for (HighsInt i = 0; i != ncols; ++i) {
-    if (!postsolveStack.isOrigCol(i)) continue;
+  for (HighsInt i : postsolveStack.getOrigCols()) {
     pseudocostup[postsolveStack.getOrigColIndex(i)] = pscost.pseudocostup[i];
     pseudocostdown[postsolveStack.getOrigColIndex(i)] =
         pscost.pseudocostdown[i];