Create sepa stats

Author: Opt-Mucca

highs/mip/HighsCutGeneration.cpp

@@ -8,7 +8,7 @@
 #include "mip/HighsCutGeneration.h"
 
 #include "../extern/pdqsort/pdqsort.h"
-#include "HighsDomain.h"
+#include "mip/HighsDomain.h"
 #include "mip/HighsMipSolverData.h"
 #include "mip/HighsTransformedLp.h"
 #include "util/HighsIntegers.h"

highs/mip/HighsCutPool.cpp

@@ -171,8 +171,8 @@ void HighsCutPool::performAging() {
 
   for (HighsInt i = 0; i != cutIndexEnd; ++i) {
     // Catch buffered changes (should only occur in parallel case)
-    // TODO: This misses the case where a cut is added then deleted before aging
-    // TODO: has been called once. We'd miss resetting the age in this case.
+    // TODO MT: Misses the case where a cut is added then deleted before aging
+    // TODO MT: has been called once. We'd miss resetting the age in this case.
     if (numLps_[i] > 0 && ages_[i] >= 0) {
       // Cut has been added to the LP, but age changes haven't been made
       --ageDistribution[ages_[i]];
@@ -245,7 +245,9 @@ void HighsCutPool::separate(const std::vector<double>& sol, HighsDomain& domain,
 
   for (HighsInt i = 0; i < nrows; ++i) {
     // cuts with an age of -1 are already in the LP and are therefore skipped
-    // TODO: Parallel case here loops over cuts potentially added in current LP
+    // TODO MT: Parallel case tries to add cuts already in current LP.
+    // TODO MT: Inefficient. Not sure what happens if added twice.
+    // TODO MT: The cut shouldn't have enough violation to be added though.
     if (ages_[i] < 0) continue;
 
     HighsInt start = matrix_.getRowStart(i);
@@ -402,8 +404,9 @@ void HighsCutPool::separate(const std::vector<double>& sol, HighsDomain& domain,
           break;
         }
       } else {
-        // TODO MT: Is this safe for the future? If we copy an LP with a cut
-        // from a local pool then this is not thread safe
+        // TODO MT: This assumes the cuts in the pool are not changing during,
+        // TODO MT: this query, i.e., the worker's pool and the global pool.
+        // TODO MT: Currently safe, but doesn't generalise to all designs.
         if (getParallelism(p.second, cutset.cutindices[i],
                            cutpools[cutset.cutpools[i]]) > maxpar) {
           discard = true;
@@ -498,8 +501,7 @@ void HighsCutPool::separateLpCutsAfterRestart(HighsCutSet& cutset) {
 HighsInt HighsCutPool::addCut(const HighsMipSolver& mipsolver, HighsInt* Rindex,
                               double* Rvalue, HighsInt Rlen, double rhs,
                               bool integral, bool propagate,
-                              bool extractCliques, bool isConflict,
-                              HighsCutPool* globalpool) {
+                              bool extractCliques, bool isConflict) {
   mipsolver.mipdata_->debugSolution.checkCut(Rindex, Rvalue, Rlen, rhs);
 
   sortBuffer.resize(Rlen);
@@ -527,10 +529,6 @@ HighsInt HighsCutPool::addCut(const HighsMipSolver& mipsolver, HighsInt* Rindex,
 
   if (isDuplicate(h, normalization, Rindex, Rvalue, Rlen, rhs)) return -1;
 
-  if (globalpool != nullptr &&
-      globalpool->isDuplicate(h, normalization, Rindex, Rvalue, Rlen, rhs))
-    return -1;
-
   // if (Rlen > 0.15 * matrix_.numCols())
   //   printf("cut with len %d not propagated\n", Rlen);
   if (propagate) {
@@ -645,6 +643,8 @@ void HighsCutPool::syncCutPool(const HighsMipSolver& mipsolver,
       std::vector<double> vals(Rvalue, Rvalue + Rlen);
       syncpool.addCut(mipsolver, idxs.data(), vals.data(), Rlen, rhs_[i],
                       rowintegral[i]);
+      // TODO MT: Should I check whether the cut is accepted before changing
+      // hasSynced?
       hasSynced_[i] = true;
     }
   }

highs/mip/HighsCutPool.h

@@ -174,8 +174,7 @@ class HighsCutPool {
   HighsInt addCut(const HighsMipSolver& mipsolver, HighsInt* Rindex,
                   double* Rvalue, HighsInt Rlen, double rhs,
                   bool integral = false, bool propagate = true,
-                  bool extractCliques = true, bool isConflict = false,
-                  HighsCutPool* globalpool = nullptr);
+                  bool extractCliques = true, bool isConflict = false);
 
   HighsInt getRowLength(HighsInt row) const {
     return matrix_.getRowEnd(row) - matrix_.getRowStart(row);

highs/mip/HighsMipSolver.cpp

@@ -734,15 +734,22 @@ void HighsMipSolver::run() {
     applyTask(doSeparate, tg, true, search_indices);
     analysis_.mipTimerStop(kMipClockNodeSearchSeparation);
 
-    for (const HighsInt i : search_indices) {
-      // Sync numNeighbourhoodQueries
+    auto syncSepaStats = [&](HighsMipWorker& worker) {
       mipdata_->cliquetable.getNumNeighbourhoodQueries() +=
-          mipdata_->workers[i].numNeighbourhoodQueries;
-      mipdata_->workers[i].numNeighbourhoodQueries = 0;
-      if (mipdata_->workers[i].getGlobalDomain().infeasible()) {
-        mipdata_->workers[i].search_ptr_->cutoffNode();
+          worker.sepa_stats.numNeighbourhoodQueries;
+      worker.sepa_stats.numNeighbourhoodQueries = 0;
+      mipdata_->sepa_lp_iterations += worker.sepa_stats.sepa_lp_iterations;
+      mipdata_->total_lp_iterations += worker.sepa_stats.sepa_lp_iterations;
+      worker.sepa_stats.sepa_lp_iterations = 0;
+    };
+
+    for (const HighsInt i : search_indices) {
+      HighsMipWorker& worker = mipdata_->workers[i];
+      syncSepaStats(worker);
+      if (worker.getGlobalDomain().infeasible()) {
+        worker.search_ptr_->cutoffNode();
         analysis_.mipTimerStart(kMipClockOpenNodesToQueue1);
-        mipdata_->workers[i].search_ptr_->openNodesToQueue(mipdata_->nodequeue);
+        worker.search_ptr_->openNodesToQueue(mipdata_->nodequeue);
         analysis_.mipTimerStop(kMipClockOpenNodesToQueue1);
         mipdata_->nodequeue.clear();
         mipdata_->pruned_treeweight = 1.0;

highs/mip/HighsMipWorker.cpp

@@ -28,7 +28,6 @@ HighsMipWorker::HighsMipWorker(const HighsMipSolver& mipsolver,
   search_ptr_ =
       std::unique_ptr<HighsSearch>(new HighsSearch(*this, getPseudocost()));
   sepa_ptr_ = std::unique_ptr<HighsSeparation>(new HighsSeparation(*this));
-  numNeighbourhoodQueries = 0;
   search_ptr_->setLpRelaxation(lp_);
   sepa_ptr_->setLpRelaxation(lp_);
 }

highs/mip/HighsMipWorker.h

@@ -22,6 +22,12 @@ class HighsSearch;
 
 class HighsMipWorker {
  public:
+  struct SepaStatistics {
+    SepaStatistics() : numNeighbourhoodQueries(0), sepa_lp_iterations(0) {}
+
+    int64_t numNeighbourhoodQueries;
+    int64_t sepa_lp_iterations;
+  };
   const HighsMipSolver& mipsolver_;
   const HighsMipSolverData& mipdata_;
 
@@ -43,11 +49,10 @@ class HighsMipWorker {
   std::vector<std::tuple<std::vector<double>, double, int>> solutions_;
 
   HighsPrimalHeuristics::Statistics heur_stats;
+  SepaStatistics sepa_stats;
 
   HighsRandom randgen;
 
-  int64_t numNeighbourhoodQueries;
-
   HighsMipWorker(const HighsMipSolver& mipsolver, HighsLpRelaxation* lp,
                  HighsDomain* domain, HighsCutPool* cutpool,
                  HighsConflictPool* conflictpool, HighsPseudocost* pseudocost);
@@ -67,7 +72,6 @@ class HighsMipWorker {
 
   HighsPseudocost& getPseudocost() const { return *pseudocost_; };
 
-
   bool addIncumbent(const std::vector<double>& sol, double solobj,
                     int solution_source);
 

highs/mip/HighsSeparation.cpp

@@ -59,8 +59,6 @@ HighsInt HighsSeparation::separationRound(HighsDomain& propdomain,
     if (&propdomain == &mipdata.domain)
       mipdata.cliquetable.cleanupFixed(mipdata.domain);
 
-    // TODO: Currently adding a check for both. Should only need to check
-    // mipworker
     if (mipworker_.getGlobalDomain().infeasible()) {
       status = HighsLpRelaxation::Status::kInfeasible;
       propdomain.clearChangedCols();
@@ -85,9 +83,7 @@ HighsInt HighsSeparation::separationRound(HighsDomain& propdomain,
     return numBoundChgs;
   };
 
-  // TODO MT: Look into delta implications (probing for global info locally and
-  // buffer it)
-  // TODO MT: Disabled timers because they fail for parallel mode
+  // TODO MT: Look into delta implications
   // lp->getMipSolver().analysis_.mipTimerStart(implBoundClock);
   mipdata.implications.separateImpliedBounds(
       *lp, lp->getSolution().col_value, *mipworker_.cutpool_, mipdata.feastol,
@@ -101,15 +97,13 @@ HighsInt HighsSeparation::separationRound(HighsDomain& propdomain,
   else
     ncuts += numboundchgs;
 
-  // TODO: This can be enabled if randgen and cliquesubsumption are disabled for
-  // parallel case
   // lp->getMipSolver().analysis_.mipTimerStart(cliqueClock);
   mipdata.cliquetable.separateCliques(
       lp->getMipSolver(), sol.col_value, *mipworker_.cutpool_, mipdata.feastol,
       mipdata.parallelLockActive() ? mipworker_.randgen
                                    : mipdata.cliquetable.getRandgen(),
       mipdata.parallelLockActive()
-          ? mipworker_.numNeighbourhoodQueries
+          ? mipworker_.sepa_stats.numNeighbourhoodQueries
           : mipdata.cliquetable.getNumNeighbourhoodQueries());
   // lp->getMipSolver().analysis_.mipTimerStop(cliqueClock);
 
@@ -183,16 +177,16 @@ void HighsSeparation::separate(HighsDomain& propdomain) {
     while (lp->getObjective() < mipsolver.mipdata_->optimality_limit) {
       double lastobj = lp->getObjective();
 
-      size_t nlpiters = -lp->getNumLpIterations();
+      int64_t nlpiters = -lp->getNumLpIterations();
       HighsInt ncuts = separationRound(propdomain, status);
       nlpiters += lp->getNumLpIterations();
 
-      // replace with mipworker iterations field
-      // mipsolver.mipdata_->sepa_lp_iterations += nlpiters;
-      // mipsolver.mipdata_->total_lp_iterations += nlpiters;
-
-      // todo:ig  more stats for separation iterations?
-      mipworker_.heur_stats.lp_iterations += nlpiters;
+      if (mipsolver.mipdata_->parallelLockActive()) {
+        mipworker_.sepa_stats.sepa_lp_iterations += nlpiters;
+      } else {
+        mipsolver.mipdata_->sepa_lp_iterations += nlpiters;
+        mipsolver.mipdata_->total_lp_iterations += nlpiters;
+      }
 
       // printf("separated %" HIGHSINT_FORMAT " cuts\n", ncuts);
 
@@ -217,11 +211,7 @@ void HighsSeparation::separate(HighsDomain& propdomain) {
     //        (HighsInt)status);
     lp->performAging(true);
 
-    // mipsolver.mipdata_->cutpool.performAging();
-    // ig: using worker cutpool
-    // TODO MT: Is this thread safe? Depends if LP is only copied at the start.
-    if (!mipsolver.mipdata_->parallelLockActive()) {
-      mipworker_.cutpool_->performAging();
-    }
+    // TODO MT: If LP is only copied at start this should be thread safe.
+    mipworker_.cutpool_->performAging();
   }
 }