Update localmip. Run in parallel

Author: Opt-Mucca

highs/lp_data/HighsOptions.h

@@ -480,6 +480,7 @@ struct HighsOptionsStruct {
   bool mip_heuristic_run_root_reduced_cost;
   bool mip_heuristic_run_zi_round;
   bool mip_heuristic_run_shifting;
+  bool mip_heuristic_run_local_mip;
   double mip_min_logging_interval;
   std::string mip_lp_solver;
   std::string mip_ipm_solver;
@@ -638,6 +639,7 @@ struct HighsOptionsStruct {
         mip_heuristic_run_root_reduced_cost(false),
         mip_heuristic_run_zi_round(false),
         mip_heuristic_run_shifting(false),
+        mip_heuristic_run_local_mip(false),
         mip_min_logging_interval(0.0),
         mip_lp_solver(""),
         mip_ipm_solver(""),
@@ -1217,6 +1219,11 @@ class HighsOptions : public HighsOptionsStruct {
                                        &mip_heuristic_run_shifting, false);
     records.push_back(record_bool);
 
+    record_bool = new OptionRecordBool("mip_heuristic_run_local_mip",
+                                       "Use the Local MIP heuristic", advanced,
+                                       &mip_heuristic_run_local_mip, true);
+    records.push_back(record_bool);
+
     record_bool = new OptionRecordBool(
         "mip_allow_cut_separation_at_nodes",
         "Whether cut separation at nodes is permitted", advanced,

highs/mip/HighsMipSolverData.cpp

@@ -343,6 +343,27 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
   return HighsModelStatus::kNotset;
 }
 
+void HighsMipSolverData::startLocalMipComputation(
+    const highs::parallel::TaskGroup& taskGroup, const HighsDomain& globaldom,
+    std::vector<double>& localMipSol) {
+  localMipSuccess = false;
+  taskGroup.spawn([&]() {
+    HighsRandom randgen(mipsolver.options_mip_->random_seed + numRestarts);
+    if (heuristics.localMip(globaldom, randgen, localMipSol)) {
+      localMipSuccess = true;
+    }
+  });
+}
+
+void HighsMipSolverData::finishLocalMipComputation(
+    const highs::parallel::TaskGroup& taskGroup,
+    const std::vector<double>& localMipSol) {
+  taskGroup.sync();
+  if (localMipSuccess) {
+    trySolution(localMipSol, kSolutionSourceLocalMip);
+  }
+}
+
 void HighsMipSolverData::startAnalyticCenterComputation(
     const highs::parallel::TaskGroup& taskGroup) {
   taskGroup.spawn([&]() {
@@ -1927,6 +1948,14 @@ void HighsMipSolverData::evaluateRootNode() {
     startSymmetryDetection(tg, symData);
     analysis.mipTimerStop(kMipClockStartSymmetryDetection);
   }
+  std::vector<double> localMipSol(mipsolver.numCol());
+  HighsDomain domainSnapshot = domain;
+  if (mipsolver.options_mip_->mip_heuristic_run_local_mip &&
+      !mipsolver.submip) {
+    analysis.mipTimerStart(kMipClockStartLocalMipComputation);
+    startLocalMipComputation(tg, domainSnapshot, localMipSol);
+    analysis.mipTimerStop(kMipClockStartLocalMipComputation);
+  }
   if (compute_analytic_centre && !analyticCenterComputed) {
     if (analysis.analyse_mip_time)
       highsLogUser(
@@ -2038,8 +2067,6 @@ void HighsMipSolverData::evaluateRootNode() {
   if (mipsolver.options_mip_->mip_heuristic_run_shifting)
     heuristics.shifting(firstlpsol);
 
-  heuristics.localMip();
-
   heuristics.flushStatistics();
 
   analysis.mipTimerStart(kMipClockEvaluateRootLp);
@@ -2414,11 +2441,6 @@ void HighsMipSolverData::evaluateRootNode() {
   if (lower_bound <= upper_limit) {
     if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart &&
         mipsolver.options_mip_->presolve != kHighsOffString) {
-      if (!analyticCenterComputed && compute_analytic_centre) {
-        analysis.mipTimerStart(kMipClockFinishAnalyticCentreComputation);
-        finishAnalyticCenterComputation(tg);
-        analysis.mipTimerStop(kMipClockFinishAnalyticCentreComputation);
-      }
       double fixingRate = percentageInactiveIntegers();
       if (fixingRate >= 2.5 + 7.5 * mipsolver.submip ||
           (!mipsolver.submip && fixingRate > 0 && numRestarts == 0)) {
@@ -2441,6 +2463,13 @@ void HighsMipSolverData::evaluateRootNode() {
       }
     }
 
+    if (!mipsolver.submip &&
+        mipsolver.options_mip_->mip_heuristic_run_local_mip) {
+      analysis.mipTimerStart(kMipClockFinishLocalMipComputation);
+      finishLocalMipComputation(tg, localMipSol);
+      analysis.mipTimerStop(kMipClockFinishLocalMipComputation);
+    }
+
     if (detectSymmetries) {
       finishSymmetryDetection(tg, symData);
       analysis.mipTimerStart(kMipClockEvaluateRootLp);

highs/mip/HighsMipSolverData.h

@@ -86,6 +86,7 @@ struct HighsMipSolverData {
   bool rowMatrixSet;
   bool analyticCenterComputed;
   HighsModelStatus analyticCenterStatus;
+  bool localMipSuccess;
   bool detectSymmetries;
   HighsInt numRestarts;
   HighsInt numRestartsRoot;
@@ -173,6 +174,7 @@ struct HighsMipSolverData {
         rowMatrixSet(false),
         analyticCenterComputed(false),
         analyticCenterStatus(HighsModelStatus::kNotset),
+        localMipSuccess(false),
         detectSymmetries(false),
         numRestarts(0),
         numRestartsRoot(0),
@@ -221,6 +223,12 @@ struct HighsMipSolverData {
   HighsModelStatus feasibilityJump();
   HighsModelStatus trivialHeuristics();
 
+  void startLocalMipComputation(const highs::parallel::TaskGroup& taskGroup,
+                                const HighsDomain& globaldom,
+                                std::vector<double>& localMipSol);
+  void finishLocalMipComputation(const highs::parallel::TaskGroup& taskGroup,
+                                 const std::vector<double>& localMipSol);
+
   void startAnalyticCenterComputation(
       const highs::parallel::TaskGroup& taskGroup);
   void finishAnalyticCenterComputation(

highs/mip/HighsPrimalHeuristics.cpp

@@ -1590,9 +1590,9 @@ void HighsPrimalHeuristics::centralRounding() {
                        kSolutionSourceCentralRounding);
 }
 
-bool HighsPrimalHeuristics::localMip() {
-  const HighsDomain& globaldom = mipsolver.mipdata_->domain;
-  std::vector<double> intsol(mipsolver.numCol());
+bool HighsPrimalHeuristics::localMip(const HighsDomain& globaldom,
+                                     HighsRandom& randgen,
+                                     std::vector<double>& intsol) {
   bool found_feas_before = false;
   double feastol = mipsolver.mipdata_->feastol;
   const HighsSparseMatrix& a_matrix = mipsolver.model_->a_matrix_;
@@ -1602,6 +1602,8 @@ bool HighsPrimalHeuristics::localMip() {
   // Initialise values and reserve vectors
   HighsInt num_violated_rows_sample = 250;
   HighsInt num_satisfied_rows_sample = 20;
+  HighsInt max_iters =
+      std::min(250000, std::max(5000, mipsolver.numNonzero() / 2));
   HighsInt iters = 0;
   HighsCDouble bestobj = kHighsInf;
   HighsCDouble obj = 0;
@@ -1785,6 +1787,8 @@ bool HighsPrimalHeuristics::localMip() {
       if (delta < feastol) return 0;
     }
     if (integral) delta = std::floor(delta + feastol);
+    assert(sol[c] + dir * delta > globaldom.col_lower_[c] - feastol);
+    assert(sol[c] + dir * delta < globaldom.col_upper_[c] + feastol);
     return delta * dir;
   };
 
@@ -1809,6 +1813,10 @@ bool HighsPrimalHeuristics::localMip() {
       }
     }
     if (best_col == -1) return false;
+    assert(sol[best_col] + one_opt_deltas[best_col] >
+           globaldom.col_lower_[best_col] - feastol);
+    assert(sol[best_col] + one_opt_deltas[best_col] <
+           globaldom.col_upper_[best_col] + feastol);
     apply_move(best_col, one_opt_deltas[best_col]);
     one_opt_deltas[best_col] = one_opt_calc_col_delta(best_col);
     HighsInt start = a_matrix.start_[best_col];
@@ -2091,7 +2099,7 @@ bool HighsPrimalHeuristics::localMip() {
     bestobj = obj;
   };
 
-  auto terminate = [&]() { return iters > 5000; };
+  auto terminate = [&]() { return iters > max_iters; };
 
   calc_activities(false);
   bool last_iter_feas = false;
@@ -2110,7 +2118,7 @@ bool HighsPrimalHeuristics::localMip() {
       momentum = 1;
       recalc_objective();
       // If solution is improving then store it
-      // TODO: Undo the transformation here??
+      // TODO: Undo the transformation here?? (Need parallel branch changes)
       if (obj + feastol < bestobj) {
         found_feas_before = true;
         save_sol();
@@ -2145,14 +2153,7 @@ bool HighsPrimalHeuristics::localMip() {
     //      viol.viol_index.size(), found_feas_before,
     //      static_cast<double>(bestobj));
   }
-  // Try and add incumbent to storage
-  if (found_feas_before) {
-    recalc_objective();
-    mipsolver.mipdata_->addIncumbent(intsol, static_cast<double>(obj),
-                                     kSolutionSourceLocalMip);
-    return true;
-  }
-  return false;
+  return found_feas_before;
 }
 
 #if 0

highs/mip/HighsPrimalHeuristics.h

@@ -57,7 +57,8 @@ class HighsPrimalHeuristics {
 
   void centralRounding();
 
-  bool localMip();
+  bool localMip(const HighsDomain& globaldom, HighsRandom& randgen,
+                std::vector<double>& intsol);
 
   void flushStatistics();
 

highs/mip/MipTimer.h

@@ -40,12 +40,14 @@ enum iClockMip {
   // Evaluate root node
   kMipClockStartSymmetryDetection,
   kMipClockStartAnalyticCentreComputation,
+  kMipClockStartLocalMipComputation,
   kMipClockEvaluateRootLp,
   kMipClockSeparateLpCuts,
   kMipClockRandomizedRounding,
   kMipClockPerformRestart,
   kMipClockRootSeparation,
   kMipClockFinishAnalyticCentreComputation,
+  kMipClockFinishLocalMipComputation,
   kMipClockRootCentralRounding,
   kMipClockRootSeparationRound0,
   kMipClockRootHeuristicsReducedCost,
@@ -176,6 +178,8 @@ class MipTimer {
         timer_pointer->clock_def("Start symmetry detection");
     clock[kMipClockStartAnalyticCentreComputation] =
         timer_pointer->clock_def("A-centre - start");
+    clock[kMipClockStartLocalMipComputation] =
+        timer_pointer->clock_def("Local Mip - start");
     clock[kMipClockEvaluateRootLp] =
         timer_pointer->clock_def("Evaluate root LP");
     clock[kMipClockSeparateLpCuts] =
@@ -188,6 +192,8 @@ class MipTimer {
         timer_pointer->clock_def("Root separation");
     clock[kMipClockFinishAnalyticCentreComputation] =
         timer_pointer->clock_def("A-centre - finish");
+    clock[kMipClockFinishLocalMipComputation] =
+        timer_pointer->clock_def("Local Mip - finish");
     clock[kMipClockRootCentralRounding] =
         timer_pointer->clock_def("Root central rounding");
     clock[kMipClockRootSeparationRound0] =
@@ -414,12 +420,14 @@ class MipTimer {
     const std::vector<HighsInt> mip_clock_list{
         kMipClockStartSymmetryDetection,
         kMipClockStartAnalyticCentreComputation,
+        kMipClockStartLocalMipComputation,
         kMipClockEvaluateRootLp,
         kMipClockSeparateLpCuts,
         kMipClockRandomizedRounding,
         kMipClockPerformRestart,
         kMipClockRootSeparation,
         kMipClockFinishAnalyticCentreComputation,
+        kMipClockFinishLocalMipComputation,
         kMipClockRootCentralRounding,
         kMipClockRootSeparationRound0,
         kMipClockRootHeuristicsReducedCost,
@@ -509,12 +517,14 @@ class MipTimer {
     const std::vector<HighsInt> mip_clock_list{
         kMipClockStartSymmetryDetection,
         kMipClockStartAnalyticCentreComputation,
+        kMipClockStartLocalMipComputation,
         kMipClockEvaluateRootLp,
         kMipClockSeparateLpCuts,
         kMipClockRandomizedRounding,
         kMipClockPerformRestart,
         kMipClockRootSeparation,
         kMipClockFinishAnalyticCentreComputation,
+        kMipClockFinishLocalMipComputation,
         kMipClockRootCentralRounding,
         kMipClockRootSeparationRound0,
         kMipClockRootHeuristicsReducedCost,