MIP Worker Refactor

check/TestCheckSolution.cpp

@@ -87,6 +87,7 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   HighsSolution optimal_solution = highs.getSolution();
 
   HighsInt scratch_num_nodes = info.mip_node_count;
+  HighsInt scratch_num_simplex = info.simplex_iteration_count;
   if (dev_run) printf("Num nodes = %d\n", int(scratch_num_nodes));
 
   std::string solution_file = test_name + model + ".sol";
@@ -113,7 +114,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
 
     highs.run();
     if (dev_run) printf("Num nodes = %d\n", int(info.mip_node_count));
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
 
@@ -132,7 +135,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
 
     highs.run();
     if (dev_run) printf("Num nodes = %d\n", int(info.mip_node_count));
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
 
@@ -162,7 +167,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
 
     highs.run();
     if (dev_run) printf("Num nodes = %d\n", int(info.mip_node_count));
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
 
@@ -185,7 +192,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
 
     highs.run();
     if (dev_run) printf("Num nodes = %d\n", int(info.mip_node_count));
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
 
@@ -232,7 +241,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
     return_status = highs.setSolution(starting_solution);
     REQUIRE(return_status == HighsStatus::kOk);
     highs.run();
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
 
@@ -283,7 +294,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
     return_status = highs.setSolution(num_entries, index.data(), value.data());
     REQUIRE(return_status == HighsStatus::kOk);
     highs.run();
-    REQUIRE(info.mip_node_count != scratch_num_nodes);
+    const bool different_search = info.mip_node_count != scratch_num_nodes ||
+            info.simplex_iteration_count != scratch_num_simplex;
+    REQUIRE(different_search);
     highs.clear();
   }
   assert(other_tests);

check/TestMipSolver.cpp

@@ -1404,3 +1404,25 @@ TEST_CASE("issue-2173", "[highs_test_mip_solver]") {
   const double optimal_objective = -26770.8075489;
   solve(highs, kHighsOnString, require_model_status, optimal_objective);
 }
+
+TEST_CASE("parallel-mip-determinism", "[highs_test_mip_solver]") {
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/bell5.mps";
+  HighsInt num_runs = 6;
+  std::vector<HighsInt> lp_iters(num_runs);
+  for (HighsInt i = 0; i < num_runs; i++) {
+    Highs highs;
+    highs.setOptionValue("output_flag", dev_run);
+    highs.setOptionValue("mip_rel_gap", 0);
+    highs.setOptionValue("threads", 2);
+    highs.setOptionValue("mip_search_concurrency", 2);
+    if (i % 2 == 0) highs.setOptionValue("mip_search_simulate_concurrency", 1);
+    highs.readModel(filename);
+    const HighsModelStatus require_model_status = HighsModelStatus::kOptimal;
+    const double optimal_objective = 8966406.491519;
+    solve(highs, kHighsOffString, require_model_status, optimal_objective);
+    lp_iters[i] = highs.getInfo().simplex_iteration_count;
+    if (i > 0) {
+      REQUIRE(lp_iters[i] == lp_iters[0]);
+    }
+  }
+}

cmake/sources-python.cmake

@@ -218,6 +218,7 @@ set(highs_sources_python
     highs/mip/HighsMipAnalysis.cpp
     highs/mip/HighsMipSolver.cpp
     highs/mip/HighsMipSolverData.cpp
+    highs/mip/HighsMipWorker.cpp
     highs/mip/HighsModkSeparator.cpp
     highs/mip/HighsNodeQueue.cpp
     highs/mip/HighsObjectiveFunction.cpp
@@ -342,6 +343,7 @@ set(highs_headers_python
     highs/mip/HighsMipAnalysis.h
     highs/mip/HighsMipSolver.h
     highs/mip/HighsMipSolverData.h
+    highs/mip/HighsMipWorker.h
     highs/mip/HighsModkSeparator.h
     highs/mip/HighsNodeQueue.h
     highs/mip/HighsObjectiveFunction.h

cmake/sources.cmake

@@ -372,6 +372,7 @@ set(highs_sources
     mip/HighsMipAnalysis.cpp
     mip/HighsMipSolver.cpp
     mip/HighsMipSolverData.cpp
+    mip/HighsMipWorker.cpp
     mip/HighsModkSeparator.cpp
     mip/HighsNodeQueue.cpp
     mip/HighsObjectiveFunction.cpp
@@ -499,6 +500,7 @@ set(highs_headers
     mip/HighsMipAnalysis.h
     mip/HighsMipSolver.h
     mip/HighsMipSolverData.h
+    mip/HighsMipWorker.h
     mip/HighsModkSeparator.h
     mip/HighsNodeQueue.h
     mip/HighsObjectiveFunction.h

highs/lp_data/HConst.h

@@ -37,6 +37,7 @@ const std::string kHighsChooseString = "choose";
 const std::string kHighsOnString = "on";
 const HighsInt kHighsMaxStringLength = 512;
 const HighsInt kSimplexConcurrencyLimit = 8;
+const HighsInt kMipSearchConcurrencyLimit = 8;
 const double kRunningAverageMultiplier = 0.05;
 
 const double kExcessivelySmallObjectiveCoefficient = 1e-4;

highs/lp_data/HighsOptions.h

@@ -504,6 +504,8 @@ struct HighsOptionsStruct {
   std::string mip_improving_solution_file;
   bool mip_root_presolve_only;
   HighsInt mip_lifting_for_probing;
+  HighsInt mip_search_concurrency;
+  bool mip_search_simulate_concurrency;
   bool mip_allow_cut_separation_at_nodes;
 
   // Logging callback identifiers
@@ -667,6 +669,8 @@ struct HighsOptionsStruct {
         mip_improving_solution_file(""),
         mip_root_presolve_only(false),
         mip_lifting_for_probing(-1),
+        mip_search_concurrency(0),
+        mip_search_simulate_concurrency(false),
         // clang-format off
         mip_allow_cut_separation_at_nodes(true) {};
   // clang-format on
@@ -1284,6 +1288,18 @@ class HighsOptions : public HighsOptionsStruct {
         kHighsInf);
     records.push_back(record_double);
 
+    record_int = new OptionRecordInt(
+        "mip_search_concurrency",
+        "Number of workers to create per thread for concurrent MIP search",
+        advanced, &mip_search_concurrency, 0, 2, kMipSearchConcurrencyLimit);
+    records.push_back(record_int);
+
+    record_bool = new OptionRecordBool(
+        "mip_search_simulate_concurrency",
+        "Simulate MIP search concurrency on a single thread", advanced,
+        &mip_search_simulate_concurrency, false);
+    records.push_back(record_bool);
+
     record_int = new OptionRecordInt(
         "ipm_iteration_limit", "Iteration limit for IPM solver", advanced,
         &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf);

highs/meson.build

@@ -235,6 +235,7 @@ _srcs = [
     'mip/HighsMipAnalysis.cpp',
     'mip/HighsMipSolver.cpp',
     'mip/HighsMipSolverData.cpp',
+    'mip/HighsMipWorker.cpp',
     'mip/HighsModkSeparator.cpp',
     'mip/HighsNodeQueue.cpp',
     'mip/HighsObjectiveFunction.cpp',

highs/mip/HighsCliqueTable.cpp

@@ -647,7 +647,7 @@ bool HighsCliqueTable::processNewEdge(HighsDomain& globaldom, CliqueVar v1,
 void HighsCliqueTable::addClique(const HighsMipSolver& mipsolver,
                                  CliqueVar* cliquevars, HighsInt numcliquevars,
                                  bool equality, HighsInt origin) {
-  HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
   mipsolver.mipdata_->debugSolution.checkClique(cliquevars, numcliquevars);
   const HighsInt maxNumCliqueVars = 100;
 
@@ -841,7 +841,7 @@ void HighsCliqueTable::extractCliques(
     HighsInt nbin, std::vector<HighsInt>& perm, std::vector<CliqueVar>& clique,
     double feastol) {
   HighsImplications& implics = mipsolver.mipdata_->implications;
-  HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
 
   perm.resize(inds.size());
   std::iota(perm.begin(), perm.end(), 0);
@@ -1081,7 +1081,7 @@ void HighsCliqueTable::extractCliquesFromCut(const HighsMipSolver& mipsolver,
   if (isFull()) return;
 
   HighsImplications& implics = mipsolver.mipdata_->implications;
-  HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
 
   const double feastol = mipsolver.mipdata_->feastol;
 
@@ -1276,7 +1276,7 @@ void HighsCliqueTable::extractCliques(HighsMipSolver& mipsolver,
 
   double rhs;
 
-  HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
 
   for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     HighsInt start = mipsolver.mipdata_->ARstart_[i];
@@ -1380,7 +1380,7 @@ void HighsCliqueTable::extractObjCliques(HighsMipSolver& mipsolver) {
   HighsInt nbin =
       mipsolver.mipdata_->objectiveFunction.getNumBinariesInObjective();
   if (nbin <= 1) return;
-  HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
   if (globaldom.getObjectiveLowerBound() == -kHighsInf) return;
 
   const double* vals;
@@ -1603,15 +1603,17 @@ void HighsCliqueTable::vertexInfeasible(HighsDomain& globaldom, HighsInt col,
 
 void HighsCliqueTable::separateCliques(const HighsMipSolver& mipsolver,
                                        const std::vector<double>& sol,
-                                       HighsCutPool& cutpool, double feastol) {
+                                       HighsCutPool& cutpool, double feastol,
+                                       HighsRandom& randgen,
+                                       int64_t& localNumNeighbourhoodQueries) {
   BronKerboschData data(sol);
   data.feastol = feastol;
   data.maxNeighbourhoodQueries = 1000000 +
                                  int64_t{100} * mipsolver.numNonzero() +
                                  mipsolver.mipdata_->total_lp_iterations * 1000;
   if (numNeighbourhoodQueries > data.maxNeighbourhoodQueries) return;
   data.maxNeighbourhoodQueries -= numNeighbourhoodQueries;
-  const HighsDomain& globaldom = mipsolver.mipdata_->domain;
+  const HighsDomain& globaldom = mipsolver.mipdata_->getDomain();
 
   for (HighsInt i : mipsolver.mipdata_->integral_cols) {
     if (colsubstituted[i] || colDeleted[i]) continue;
@@ -1698,9 +1700,9 @@ void HighsCliqueTable::separateCliques(const HighsMipSolver& mipsolver,
                    false, false);
   }
 
-  numNeighbourhoodQueries += data.numNeighbourhoodQueries;
+  localNumNeighbourhoodQueries += data.numNeighbourhoodQueries;
 
-  if (runcliquesubsumption) {
+  if (runcliquesubsumption && &randgen == &this->randgen) {
     for (std::vector<CliqueVar>& clique : data.cliques) {
       HighsInt nremoved = runCliqueSubsumption(globaldom, clique);
 
@@ -1841,7 +1843,7 @@ void HighsCliqueTable::cleanupFixed(HighsDomain& globaldom) {
   if (nfixings != oldnfixings) propagateAndCleanup(globaldom);
 }
 
-HighsInt HighsCliqueTable::getNumImplications(HighsInt col) {
+HighsInt HighsCliqueTable::getNumImplications(HighsInt col) const {
   // first count all cliques as one implication, so that cliques of size two
   // are accounted for already
   HighsInt i0 = CliqueVar(col, 0).index();
@@ -1860,7 +1862,7 @@ HighsInt HighsCliqueTable::getNumImplications(HighsInt col) {
   return numimplics;
 }
 
-HighsInt HighsCliqueTable::getNumImplications(HighsInt col, bool val) {
+HighsInt HighsCliqueTable::getNumImplications(HighsInt col, bool val) const {
   HighsInt iVal = CliqueVar(col, val).index();
 
   // each size two clique is one implication

highs/mip/HighsCliqueTable.h

@@ -182,6 +182,10 @@ class HighsCliqueTable {
 
   HighsInt getNumEntries() const { return numEntries; }
 
+  HighsRandom& getRandgen() { return randgen; }
+
+  int64_t& getNumNeighbourhoodQueries() { return numNeighbourhoodQueries; }
+
   HighsInt partitionNeighbourhood(std::vector<HighsInt>& neighbourhoodInds,
                                   int64_t& numNeighbourhoodqueries, CliqueVar v,
                                   CliqueVar* q, HighsInt N) const;
@@ -287,7 +291,8 @@ class HighsCliqueTable {
 
   void separateCliques(const HighsMipSolver& mipsolver,
                        const std::vector<double>& sol, HighsCutPool& cutpool,
-                       double feastol);
+                       double feastol, HighsRandom& randgen,
+                       int64_t& localNumNeighbourhoodQueries);
 
   std::vector<std::vector<CliqueVar>> separateCliques(
       const std::vector<double>& sol, const HighsDomain& globaldom,
@@ -300,9 +305,9 @@ class HighsCliqueTable {
 
   void addImplications(HighsDomain& domain, HighsInt col, HighsInt val);
 
-  HighsInt getNumImplications(HighsInt col);
+  HighsInt getNumImplications(HighsInt col) const;
 
-  HighsInt getNumImplications(HighsInt col, bool val);
+  HighsInt getNumImplications(HighsInt col, bool val) const;
 
   void runCliqueMerging(HighsDomain& globaldomain);