Aggregate implications

check/TestPresolveRules.cpp

@@ -98,6 +98,200 @@ TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   h.resetGlobalScheduler(true);
 }
 
+TEST_CASE("test-implied-bound-aggregation", "[highs_test_presolve_rules]") {
+  // Three VLB rows on a continuous variable z with binary variables
+  // x0, x1, x2 in a clique (from assignment equation). Probing should
+  // discover the clique and aggregation should merge the three VLBs
+  // into one stronger row.
+  //
+  // min z
+  // x0 + x1 + x2 = 1
+  // z >= 10*x0  (stored as -z + 10*x0 <= 0)
+  // z >= 20*x1  (stored as -z + 20*x1 <= 0)
+  // z >= 30*x2  (stored as -z + 30*x2 <= 0)
+  // z >= 0, x0,x1,x2 binary
+  HighsLp lp;
+  lp.num_col_ = 4;
+  lp.num_row_ = 4;
+  lp.sense_ = ObjSense::kMinimize;
+  lp.col_cost_ = {1, 0, 0, 0};
+  lp.col_lower_ = {0, 0, 0, 0};
+  lp.col_upper_ = {kHighsInf, 1, 1, 1};
+  lp.integrality_ = {HighsVarType::kContinuous, HighsVarType::kInteger,
+                     HighsVarType::kInteger, HighsVarType::kInteger};
+  lp.row_lower_ = {1, -kHighsInf, -kHighsInf, -kHighsInf};
+  lp.row_upper_ = {1, 0, 0, 0};
+  lp.a_matrix_.format_ = MatrixFormat::kColwise;
+  lp.a_matrix_.num_col_ = 4;
+  lp.a_matrix_.num_row_ = 4;
+  lp.a_matrix_.start_ = {0, 3, 5, 7, 9};
+  lp.a_matrix_.index_ = {1, 2, 3, 0, 1, 0, 2, 0, 3};
+  lp.a_matrix_.value_ = {-1, -1, -1, 1, 10, 1, 20, 1, 30};
+
+  highs::parallel::initialize_scheduler(1);
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("presolve_rule_test", kPresolveRuleProbing);
+  highs.passModel(lp);
+
+  HighsCallback callback(&highs);
+  const HighsOptions& options = highs.getOptions();
+  HighsSolution solution;
+  HighsProfiling profiling;
+
+  HighsMipSolver mipsolver(callback, options, lp, solution);
+  mipsolver.timer_.start();
+  profiling.initialize(mipsolver.timer_, true, true);
+  mipsolver.setProfiling(&profiling);
+  mipsolver.mipdata_ =
+      std::unique_ptr<HighsMipSolverData>(new HighsMipSolverData(mipsolver));
+  mipsolver.mipdata_->init();
+  mipsolver.mipdata_->setupDomainPropagation();
+
+  // Populate VLBs: z >= 10*x0 (row 1), z >= 20*x1 (row 2), z >= 30*x2 (row 3)
+  HighsImplications& implications = mipsolver.mipdata_->implications;
+  implications.addVLB(0, 1, 10.0, 0.0, 1);
+  implications.addVLB(0, 2, 20.0, 0.0, 2);
+  implications.addVLB(0, 3, 30.0, 0.0, 3);
+
+  presolve::HighsPostsolveStack& postsolve_stack =
+      mipsolver.mipdata_->postSolveStack;
+
+  presolve::HPresolve presolve;
+  presolve.setInput(mipsolver, -1);
+  REQUIRE(presolve.okSetupPresolveDataStructures());
+  HighsModelStatus status = presolve.run(postsolve_stack);
+  REQUIRE(status == HighsModelStatus::kNotset);
+
+  // The three VLBs form a clique, so aggregation merges them into one
+  // row: z - 10*x0 - 20*x1 - 30*x2 >= 0
+  const HighsLp& presolved = *mipsolver.model_;
+  HighsInt agg_row = -1;
+  for (HighsInt i = 0; i < mipsolver.numRow(); i++) {
+    if (postsolve_stack.getOrigRowIndex(i) == 0) continue;
+    REQUIRE(agg_row == -1);
+    agg_row = i;
+  }
+  REQUIRE(agg_row >= 0);
+  REQUIRE(presolved.row_lower_[agg_row] == 0.0);
+  REQUIRE(presolved.row_upper_[agg_row] == kHighsInf);
+
+  std::vector<double> coeffs(4, 0.0);
+  for (HighsInt j = 0; j < presolved.num_col_; j++) {
+    for (HighsInt p = presolved.a_matrix_.start_[j];
+         p < presolved.a_matrix_.start_[j + 1]; p++) {
+      if (presolved.a_matrix_.index_[p] != agg_row) continue;
+      coeffs[postsolve_stack.getOrigColIndex(j)] =
+          presolved.a_matrix_.value_[p];
+    }
+  }
+  REQUIRE(coeffs[0] == 1.0);
+  REQUIRE(coeffs[1] == -10.0);
+  REQUIRE(coeffs[2] == -20.0);
+  REQUIRE(coeffs[3] == -30.0);
+
+  HighsTaskExecutor::shutdown(true);
+}
+
+TEST_CASE("test-implied-bound-aggregation-vub", "[highs_test_presolve_rules]") {
+  // Three VUB rows on a continuous variable z with binary variables
+  // x0, x1, x2 where x0+x1+x2 = 2 (so at most one xi can be 0).
+  // The complement literals (1-x0), (1-x1), (1-x2) form a clique.
+  // VUBs have positive coefficients so the activating literal is x=0
+  // (val=0), exercising the complement offset path in mergeCliques.
+  //
+  // max z
+  // x0 + x1 + x2 = 2
+  // z <= 10*x0 + 20  (stored as z - 10*x0 <= 20)
+  // z <= 20*x1 + 10  (stored as z - 20*x1 <= 10)
+  // z <= 30*x2        (stored as z - 30*x2 <= 0)
+  // 0 <= z <= 30, x0,x1,x2 binary
+  HighsLp lp;
+  lp.num_col_ = 4;
+  lp.num_row_ = 4;
+  lp.sense_ = ObjSense::kMaximize;
+  lp.col_cost_ = {1, 0, 0, 0};
+  lp.col_lower_ = {0, 0, 0, 0};
+  lp.col_upper_ = {30, 1, 1, 1};
+  lp.integrality_ = {HighsVarType::kContinuous, HighsVarType::kInteger,
+                     HighsVarType::kInteger, HighsVarType::kInteger};
+  lp.row_lower_ = {2, -kHighsInf, -kHighsInf, -kHighsInf};
+  lp.row_upper_ = {2, 20, 10, 0};
+  lp.a_matrix_.format_ = MatrixFormat::kColwise;
+  lp.a_matrix_.num_col_ = 4;
+  lp.a_matrix_.num_row_ = 4;
+  lp.a_matrix_.start_ = {0, 3, 5, 7, 9};
+  lp.a_matrix_.index_ = {1, 2, 3, 0, 1, 0, 2, 0, 3};
+  lp.a_matrix_.value_ = {1, 1, 1, 1, -10, 1, -20, 1, -30};
+
+  highs::parallel::initialize_scheduler(1);
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("presolve_rule_test", kPresolveRuleProbing);
+  highs.passModel(lp);
+
+  HighsCallback callback(&highs);
+  const HighsOptions& options = highs.getOptions();
+  HighsSolution solution;
+  HighsProfiling profiling;
+
+  HighsMipSolver mipsolver(callback, options, lp, solution);
+  mipsolver.timer_.start();
+  profiling.initialize(mipsolver.timer_, true, true);
+  mipsolver.setProfiling(&profiling);
+  mipsolver.mipdata_ =
+      std::unique_ptr<HighsMipSolverData>(new HighsMipSolverData(mipsolver));
+  mipsolver.mipdata_->init();
+  mipsolver.mipdata_->setupDomainPropagation();
+
+  // Populate VUBs with positive coef (val=0 complement path):
+  // z <= 10*x0 + 20 (row 1), z <= 20*x1 + 10 (row 2), z <= 30*x2 (row 3)
+  HighsImplications& implications = mipsolver.mipdata_->implications;
+  implications.addVUB(0, 1, 10.0, 20.0, 30.0, false, 1);
+  implications.addVUB(0, 2, 20.0, 10.0, 30.0, false, 2);
+  implications.addVUB(0, 3, 30.0, 0.0, 30.0, false, 3);
+
+  presolve::HighsPostsolveStack& postsolve_stack =
+      mipsolver.mipdata_->postSolveStack;
+
+  presolve::HPresolve presolve;
+  presolve.setInput(mipsolver, -1);
+  REQUIRE(presolve.okSetupPresolveDataStructures());
+  HighsModelStatus status = presolve.run(postsolve_stack);
+  REQUIRE(status == HighsModelStatus::kNotset);
+
+  // The three VUBs form a clique, so aggregation merges them into one
+  // row: z - 10*x0 - 20*x1 - 30*x2 <= -30
+  const HighsLp& presolved = *mipsolver.model_;
+  HighsInt agg_row = -1;
+  for (HighsInt i = 0; i < mipsolver.numRow(); i++) {
+    if (postsolve_stack.getOrigRowIndex(i) == 0) continue;
+    REQUIRE(agg_row == -1);
+    agg_row = i;
+  }
+  REQUIRE(agg_row >= 0);
+  REQUIRE(presolved.row_lower_[agg_row] == -kHighsInf);
+  REQUIRE(presolved.row_upper_[agg_row] == -30.0);
+
+  std::vector<double> coeffs(4, 0.0);
+  for (HighsInt j = 0; j < presolved.num_col_; j++) {
+    for (HighsInt p = presolved.a_matrix_.start_[j];
+         p < presolved.a_matrix_.start_[j + 1]; p++) {
+      if (presolved.a_matrix_.index_[p] != agg_row) continue;
+      coeffs[postsolve_stack.getOrigColIndex(j)] =
+          presolved.a_matrix_.value_[p];
+    }
+  }
+  REQUIRE(coeffs[0] == 1.0);
+  REQUIRE(coeffs[1] == -10.0);
+  REQUIRE(coeffs[2] == -20.0);
+  REQUIRE(coeffs[3] == -30.0);
+
+  HighsTaskExecutor::shutdown(true);
+}
+
 /*
 TEST_CASE("test-weakly-dominated-col-upper", "[highs_test_presolve_rules]") {
   Highs h;

highs/mip/HighsCliqueTable.cpp

@@ -985,6 +985,33 @@ void HighsCliqueTable::extractCliques(
   }
 }
 
+template <bool Sort, typename Comparator>
+HighsInt HighsCliqueTable::extendClique(
+    std::vector<HighsInt>& neighbourhoodInds, std::vector<CliqueVar>& clqVars,
+    HighsInt seedPos, HighsInt extensionEnd, const Comparator& candidateOrder) {
+  for (HighsInt i = seedPos; i < extensionEnd; ++i) {
+    if (Sort)
+      pdqsort_branchless(clqVars.begin() + i, clqVars.begin() + extensionEnd,
+                         candidateOrder);
+    CliqueVar v = clqVars[i];
+    HighsInt extensionStart = i + 1;
+    extensionEnd =
+        partitionNeighbourhood(neighbourhoodInds, numNeighbourhoodQueries, v,
+                               clqVars.data() + extensionStart,
+                               extensionEnd - extensionStart) +
+        extensionStart;
+  }
+  return extensionEnd;
+}
+
+HighsInt HighsCliqueTable::extendClique(
+    std::vector<HighsInt>& neighbourhoodInds, std::vector<CliqueVar>& clqVars,
+    HighsInt seedPos, HighsInt extensionEnd) {
+  auto noop = [](const CliqueVar&, const CliqueVar&) { return false; };
+  return extendClique<false>(neighbourhoodInds, clqVars, seedPos, extensionEnd,
+                             noop);
+}
+
 void HighsCliqueTable::cliquePartition(std::vector<CliqueVar>& clqVars,
                                        std::vector<HighsInt>& partitionStart) {
   randgen.shuffle(clqVars.data(), clqVars.size());
@@ -995,22 +1022,12 @@ void HighsCliqueTable::cliquePartition(std::vector<CliqueVar>& clqVars,
   HighsInt numClqVars = clqVars.size();
   partitionStart.clear();
   partitionStart.reserve(clqVars.size());
-  HighsInt extensionEnd = numClqVars;
-  partitionStart.push_back(0);
-  for (HighsInt i = 0; i < numClqVars; ++i) {
-    if (i == extensionEnd) {
-      partitionStart.push_back(i);
-      extensionEnd = numClqVars;
-    }
-    CliqueVar v = clqVars[i];
-    HighsInt extensionStart = i + 1;
-    extensionEnd =
-        partitionNeighbourhood(neighbourhoodInds, numNeighbourhoodQueries, v,
-                               clqVars.data() + extensionStart,
-                               extensionEnd - extensionStart) +
-        extensionStart;
-  }
 
+  HighsInt pos = 0;
+  while (pos < numClqVars) {
+    partitionStart.push_back(pos);
+    pos = extendClique(neighbourhoodInds, clqVars, pos, numClqVars);
+  }
   partitionStart.push_back(numClqVars);
 }
 
@@ -1062,6 +1079,60 @@ void HighsCliqueTable::cliquePartition(const std::vector<double>& objective,
   partitionStart.push_back(numClqVars);
 }
 
+void HighsCliqueTable::cliqueCover(std::vector<CliqueVar>& clqVars,
+                                   std::vector<std::vector<CliqueVar>>& cover) {
+  cover.clear();
+  if (clqVars.empty()) return;
+
+  auto candidateOrder = [&](const CliqueVar& a, const CliqueVar& b) {
+    return numcliquesvar[a.index()] > numcliquesvar[b.index()];
+  };
+
+  std::vector<HighsInt> neighbourhoodInds;
+  neighbourhoodInds.reserve(clqVars.size());
+
+  HighsInt numClqVars = clqVars.size();
+  HighsInt pos = 0;
+
+  while (pos < numClqVars) {
+    HighsInt cliqueEnd = extendClique(neighbourhoodInds, clqVars, pos,
+                                      numClqVars, candidateOrder);
+
+    std::vector<CliqueVar> clique(clqVars.begin() + pos,
+                                  clqVars.begin() + cliqueEnd);
+
+    for (HighsInt i = 0; i < pos; ++i) {
+      if (std::all_of(clique.begin(), clique.end(),
+                      [&](const CliqueVar& clqvar) {
+                        return haveCommonClique(clqVars[i], clqvar);
+                      }))
+        clique.push_back(clqVars[i]);
+    }
+
+    cover.push_back(std::move(clique));
+    pos = cliqueEnd;
+  }
+
+  // count variable occurrences
+  HighsHashTable<HighsInt, HighsInt> count;
+  for (const auto& clique : cover)
+    for (const auto& var : clique) count[var.index()]++;
+
+  // remove cliques where every member appears in at least one other clique
+  size_t kept = 0;
+  for (auto& clique : cover) {
+    if (std::all_of(clique.begin(), clique.end(), [&](const CliqueVar& var) {
+          return count[var.index()] >= 2;
+        })) {
+      for (const auto& var : clique) count[var.index()]--;
+    } else {
+      std::swap(cover[kept], clique);
+      kept++;
+    }
+  }
+  cover.resize(kept);
+}
+
 bool HighsCliqueTable::foundCover(HighsDomain& globaldom, CliqueVar v1,
                                   CliqueVar v2) {
   HighsInt commonclique = findCommonCliqueId(v1, v2);

highs/mip/HighsCliqueTable.h

@@ -160,6 +160,16 @@ class HighsCliqueTable {
                           int64_t& numNeighbourhoodqueries, CliqueVar v,
                           CliqueVar* q, HighsInt N) const;
 
+  template <bool Sort = true, typename Comparator>
+  HighsInt extendClique(std::vector<HighsInt>& neighbourhoodInds,
+                        std::vector<CliqueVar>& clqVars, HighsInt seedPos,
+                        HighsInt extensionEnd,
+                        const Comparator& candidateOrder);
+
+  HighsInt extendClique(std::vector<HighsInt>& neighbourhoodInds,
+                        std::vector<CliqueVar>& clqVars, HighsInt seedPos,
+                        HighsInt extensionEnd);
+
  public:
   int64_t numNeighbourhoodQueries;
 
@@ -263,6 +273,9 @@ class HighsCliqueTable {
                        std::vector<CliqueVar>& clqVars,
                        std::vector<HighsInt>& partitionStart);
 
+  void cliqueCover(std::vector<CliqueVar>& clqVars,
+                   std::vector<std::vector<CliqueVar>>& cover);
+
   bool foundCover(HighsDomain& globaldom, CliqueVar v1, CliqueVar v2);
 
   void extractCliques(HighsMipSolver& mipsolver, bool transformRows = true);