Cap probing work and store reduced-cost bounds in a flat vector.

highs/lp_data/HighsSolution.cpp

@@ -66,8 +66,11 @@ void getLpKktFailures(const HighsOptions& options, const HighsLp& lp,
   getKktFailures(options, false, lp, lp.col_cost_, solution, highs_info,
                  get_residuals);
   getPrimalDualBasisErrors(options, lp, solution, basis, primal_dual_errors);
-  getPrimalDualGlpsolErrors(options, lp, lp.col_cost_, solution,
-                            primal_dual_errors);
+  // glpsol-style errors are diagnostic only (debug reporting); skip in the
+  // hot path of repeated LP solves
+  if (options.highs_debug_level > kHighsDebugLevelCheap)
+    getPrimalDualGlpsolErrors(options, lp, lp.col_cost_, solution,
+                              primal_dual_errors);
 }
 
 void getKktFailures(const HighsOptions& options, const bool is_qp,

highs/mip/HighsDomain.cpp

@@ -2390,8 +2390,38 @@ bool HighsDomain::propagate() {
     changedboundsize = std::max(
         changedboundsize, cutpoolprop.cutpool->getMatrix().nonzeroCapacity());
 
-  std::unique_ptr<HighsDomainChange[]> changedbounds(
-      new HighsDomainChange[changedboundsize]);
+  // 10xb: reuse a thread-local scratch buffer instead of allocating
+  // 2*nnz entries (often > mmap threshold => mmap/munmap + page faults)
+  // on every propagate() call. Falls back to a fresh allocation when
+  // re-entered.
+  struct PropBuf10 {
+    std::unique_ptr<HighsDomainChange[]> buf;
+    size_t cap = 0;
+    bool inUse = false;
+  };
+  static thread_local PropBuf10 propBuf10;
+  std::unique_ptr<HighsDomainChange[]> ownedBuf10;
+  HighsDomainChange* changedbounds;
+  bool usingShared10 = false;
+  if (!propBuf10.inUse) {
+    if (propBuf10.cap < changedboundsize) {
+      propBuf10.buf.reset(new HighsDomainChange[changedboundsize]);
+      propBuf10.cap = changedboundsize;
+    }
+    propBuf10.inUse = true;
+    usingShared10 = true;
+    changedbounds = propBuf10.buf.get();
+  } else {
+    ownedBuf10.reset(new HighsDomainChange[changedboundsize]);
+    changedbounds = ownedBuf10.get();
+  }
+  struct Release10 {
+    bool on;
+    PropBuf10& b;
+    ~Release10() {
+      if (on) b.inUse = false;
+    }
+  } release10{usingShared10, propBuf10};
 
   while (havePropagationRows()) {
     if (objProp_.isActive()) objProp_.propagate();

highs/mip/HighsImplications.cpp

@@ -1,3 +1,5 @@
+#include <cstdlib>
+#include <cstdio>
 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 /*                                                                       */
 /*    This file is part of the HiGHS linear optimization suite           */
@@ -67,6 +69,7 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) {
 
   HighsInt stackimplicend = domchgstack.size();
   numImplications += stackimplicend;
+  probeWork10 += numImplications;
   mipsolver.mipdata_->getPseudoCost().addInferenceObservation(
       col, numImplications, val);
 

highs/mip/HighsImplications.h

@@ -32,6 +32,11 @@ class HighsImplications {
   int64_t numVarBounds;
   int64_t maxVarBounds;
 
+ public:
+  int64_t probeWork10 = 0;  // 10xb: total propagated changes in probing
+
+ private:
+
   bool computeImplications(HighsInt col, bool val);
 
  public:

highs/mip/HighsMipSolver.cpp

@@ -681,12 +681,32 @@ void HighsMipSolver::run() {
     return mipdata_->workers[i].getGlobalDomain().infeasible() || pruned;
   };
 
+  // Adaptive node separation (cf. SCIP's separator frequencies and
+  // Achterberg 2007, Ch. 8; Wesselmann & Stuhl 2012 on cut-loop effort):
+  // node cuts are only worth their LP re-solves when they move the node bound.
+  // After a node round closes < 1% of the node's gap to the incumbent, skip
+  // separation on the next 2^k nodes (k = consecutive ineffective rounds,
+  // capped at 6); any effective round resets k.
+  std::vector<HighsInt> sepaFails(mipdata_->workers.size(), 0),
+      sepaSkip(mipdata_->workers.size(), 0);
   auto separateAndStoreBasis = [&](HighsInt i) -> bool {
     HighsMipWorker& worker = mipdata_->workers[i];
-    if (options_mip_->mip_allow_cut_separation_at_nodes) {
+    if (options_mip_->mip_allow_cut_separation_at_nodes && sepaSkip[i] > 0) {
+      --sepaSkip[i];
+      worker.getCutPool().performAging();
+    } else if (options_mip_->mip_allow_cut_separation_at_nodes) {
       if (!mipdata_->parallelLockActive())
         profiling_->start(kMipClockNodeSearchSeparation);
+      const double objBefore = worker.getLpRelaxation().getObjective();
       worker.sepa_ptr_->separate(worker.search_ptr_->getLocalDomain());
+      const double objAfter = worker.getLpRelaxation().getObjective();
+      const double gap = mipdata_->upper_limit - objBefore;
+      if (gap < kHighsInf && gap > 0 &&
+          !(objAfter - objBefore >= 0.01 * gap)) {
+        sepaFails[i] = std::min<HighsInt>(sepaFails[i] + 1, 6);
+        sepaSkip[i] = (HighsInt{1} << sepaFails[i]) - 1;
+      } else
+        sepaFails[i] = 0;
       if (!mipdata_->parallelLockActive())
         profiling_->stop(kMipClockNodeSearchSeparation);
     } else {

highs/mip/HighsMipSolverData.cpp

@@ -2166,7 +2166,7 @@ void HighsMipSolverData::evaluateRootNode(HighsMipWorker& worker) {
   sepa.setLpRelaxation(&getLp());
 
   while (getLp().scaledOptimal(status) &&
-         !getLp().getFractionalIntegers().empty() && stall < 3) {
+         !getLp().getFractionalIntegers().empty() && stall < 2) {
     printDisplayLine();
 
     if (checkLimits()) {

highs/mip/HighsPrimalHeuristics.cpp

@@ -1,3 +1,4 @@
+#include <cstdlib>
 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 /*                                                                       */
 /*    This file is part of the HiGHS linear optimization suite           */
@@ -79,6 +80,14 @@ bool HighsPrimalHeuristics::solveSubMip(
     double fixingRate, std::vector<double> colLower,
     std::vector<double> colUpper, HighsInt maxleaves, HighsInt maxnodes,
     HighsInt stallnodes) {
+  // Limit recursive sub-MIP nesting (SCIP disables LNS heuristics inside
+  // sub-SCIPs, cf. Berthold 2014 "Heuristic algorithms in global MINLP
+  // solvers", ch. on LNS): deep chains of RINS/RENS inside sub-MIPs cost
+  // root time without improving the outer incumbent much.
+  {
+    static const int kMaxLvl = 1;
+    if (mipsolver.submip_level >= kMaxLvl) return true;
+  }
   HighsOptions submipoptions = *mipsolver.options_mip_;
   HighsLp submip = lp;
 
@@ -312,6 +321,19 @@ class HeuristicNeighbourhood {
 };
 
 void HighsPrimalHeuristics::rootReducedCost(HighsMipWorker& worker) {
+  // Gate by model type: on models dominated by general (non-binary)
+  // integers with wide domains, fixing by lurking reduced-cost bounds only
+  // tightens bounds instead of fixing variables, so the sub-MIP neighbourhood
+  // stays large and expensive (cf. Berthold 2014: LNS pays off when the
+  // fixing rate yields a small sub-problem).
+  {
+    static const double kMaxGenInt = 0.5;
+    HighsInt ngen = 0;
+    const auto& dom = worker.getGlobalDomain();
+    for (HighsInt col : mipsolver.mipdata_->integral_cols)
+      if (dom.col_upper_[col] - dom.col_lower_[col] > 1.5) ++ngen;
+    if (ngen > kMaxGenInt * mipsolver.mipdata_->integral_cols.size()) return;
+  }
   std::vector<std::pair<double, HighsDomainChange>> lurkingBounds =
       mipsolver.mipdata_->redcostfixing.getLurkingBounds(
           mipsolver, worker.getGlobalDomain());
@@ -368,8 +390,7 @@ void HighsPrimalHeuristics::rootReducedCost(HighsMipWorker& worker) {
 
   solveSubMip(worker, *mipsolver.model_, mipsolver.mipdata_->firstrootbasis,
               fixingRate, localdom.col_lower_, localdom.col_upper_,
-              500,  // std::max(50, int(0.05 *
-                    // (mipsolver.mipdata_->num_leaves))),
+              500,
               200 + static_cast<HighsInt>(mipsolver.mipdata_->num_nodes / 20),
               12);
 }
@@ -1122,7 +1143,37 @@ void HighsPrimalHeuristics::randomizedRounding(
   HighsRandom& randgen =
       mipsolver.mipdata_->parallelLockActive() ? worker.randgen : this->randgen;
 
+  // Batch phase: integer columns that are already integral in the LP
+  // solution would be rounded to their LP value anyway; fix them all first
+  // and propagate once instead of once per column (the LP point stays
+  // feasible under these fixings). Then round the fractional ones one by one.
+  // Only on general-integer dominated models (where per-column propagation
+  // over wide domains is the bottleneck); binary models keep the original
+  // sequential scheme, whose propagation steers later roundings.
+  HighsInt ngenint = 0;
+  for (HighsInt i : intcols)
+    if (localdom.col_upper_[i] - localdom.col_lower_[i] > 1.5) ++ngenint;
+  if (2 * ngenint > (HighsInt)intcols.size()) {
+    const double ftol = mipsolver.mipdata_->feastol;
+    for (HighsInt i : intcols) {
+      double r = std::round(relaxationsol[i]);
+      if (std::abs(relaxationsol[i] - r) > ftol) continue;
+      if (localdom.col_lower_[i] == localdom.col_upper_[i]) continue;
+      if (r < localdom.col_lower_[i] || r > localdom.col_upper_[i]) continue;
+      localdom.fixCol(i, r, HighsDomain::Reason::branching());
+      if (localdom.infeasible()) break;
+    }
+    if (!localdom.infeasible()) localdom.propagate();
+    if (localdom.infeasible()) {
+      localdom.conflictAnalysis(worker.getConflictPool(),
+                                worker.getGlobalDomain(),
+                                worker.getPseudocost());
+      return;
+    }
+  }
+
   for (HighsInt i : intcols) {
+    if (localdom.col_lower_[i] == localdom.col_upper_[i]) continue;
     double intval;
     if (mipsolver.mipdata_->uplocks[i] == 0)
       intval = std::ceil(relaxationsol[i] - mipsolver.mipdata_->feastol);
@@ -1134,6 +1185,9 @@ void HighsPrimalHeuristics::randomizedRounding(
     intval = std::min(localdom.col_upper_[i], intval);
     intval = std::max(localdom.col_lower_[i], intval);
 
+    // already fixed by earlier propagation: fixCol would be a no-op and the
+    // propagation queue is unchanged, so skip both calls
+    if (localdom.col_lower_[i] == localdom.col_upper_[i]) continue;
     localdom.fixCol(i, intval, HighsDomain::Reason::branching());
     if (localdom.infeasible()) {
       localdom.conflictAnalysis(worker.getConflictPool(),

highs/mip/HighsRedcostFixing.cpp

@@ -1,3 +1,5 @@
+#include <cstdio>
+#include <cstdlib>
 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 /*                                                                       */
 /*    This file is part of the HiGHS linear optimization suite           */
@@ -45,9 +47,11 @@ void HighsRedcostFixing::propagateRootRedcost(const HighsMipSolver& mipsolver) {
         lurkingColUpper[col].begin(),
         lurkingColUpper[col].upper_bound(mipsolver.mipdata_->lower_bound));
 
+    // By the staircase invariant (see addRootRedcost) the first active entry
+    // carries the tightest bound; later active entries are implied by it.
     for (auto it =
              lurkingColLower[col].lower_bound(mipsolver.mipdata_->upper_limit);
-         it != lurkingColLower[col].end(); ++it) {
+         it != lurkingColLower[col].end(); it = lurkingColLower[col].end()) {
       if (it->second > mipsolver.mipdata_->getDomain().col_lower_[col]) {
         mipsolver.mipdata_->getDomain().changeBound(
             HighsBoundType::kLower, col, (double)it->second,
@@ -56,9 +60,11 @@ void HighsRedcostFixing::propagateRootRedcost(const HighsMipSolver& mipsolver) {
       }
     }
 
+    // By the staircase invariant (see addRootRedcost) the first active entry
+    // carries the tightest bound; later active entries are implied by it.
     for (auto it =
              lurkingColUpper[col].lower_bound(mipsolver.mipdata_->upper_limit);
-         it != lurkingColUpper[col].end(); ++it) {
+         it != lurkingColUpper[col].end(); it = lurkingColUpper[col].end()) {
       if (it->second < mipsolver.mipdata_->getDomain().col_upper_[col]) {
         mipsolver.mipdata_->getDomain().changeBound(
             HighsBoundType::kUpper, col, (double)it->second,
@@ -230,8 +236,8 @@ void HighsRedcostFixing::addRootRedcost(const HighsMipSolver& mipsolver,
       [&](HighsInt col, HighsInt direction, HighsInt bound, HighsInt otherBound,
           bool isOtherBoundFinite, double lpObjective, double redCost,
           HighsInt maxNumSteps, HighsInt maxNumStepsExp,
-          std::multimap<double, HighsInt>& lurkingBounds,
-          std::multimap<double, HighsInt>& otherLurkingBounds) {
+          HighsLurkingStair& lurkingBounds,
+          HighsLurkingStair& otherLurkingBounds) {
         if (direction * redCost == kHighsInf) {
           lurkingBounds.clear();
           otherLurkingBounds.clear();
@@ -252,7 +258,26 @@ void HighsRedcostFixing::addRootRedcost(const HighsMipSolver& mipsolver,
         double shift = direction * (1 - 10 * mipsolver.mipdata_->feastol);
         step *= direction;
 
-        for (HighsInt lurkingBound = bound;
+        // requiredCutoffBound grows linearly with the step index j, so jump
+        // straight to the first j whose cutoff can exceed the lower bound
+        // (the in-loop test still guards the boundary). Exact.
+        HighsInt startBound = bound;
+        {
+          double denom = double(step) * redCost;
+          double lbnd = mipsolver.mipdata_->lower_bound +
+                        mipsolver.mipdata_->feastol;
+          if (denom > 0 && lbnd > -kHighsInf) {
+            double j = std::floor((lbnd - lpObjective - shift * redCost) / denom) - 1;
+            if (j > 0) {
+              double maxj = double(direction * (lastBound - bound)) /
+                                double(direction * step) + 1;
+              if (j > maxj) j = maxj;
+              startBound = bound + HighsInt(j) * step;
+            }
+          }
+        }
+        size_t posIdx = 0;
+        for (HighsInt lurkingBound = startBound;
              direction * lurkingBound <= direction * lastBound;
              lurkingBound += step) {
           double fracBound = lurkingBound - bound + shift;
@@ -261,30 +286,48 @@ void HighsRedcostFixing::addRootRedcost(const HighsMipSolver& mipsolver,
               mipsolver.mipdata_->lower_bound + mipsolver.mipdata_->feastol)
             continue;
 
-          // check if we already have a better lurking bound stored
-          bool useful = true;
-          auto pos = lurkingBounds.lower_bound(requiredCutoffBound);
-          for (auto it = pos; it != lurkingBounds.end(); ++it) {
-            useful =
-                direction * it->second >= direction * (lurkingBound + step);
-            if (!useful) break;
-          }
-          if (!useful) continue;
-
-          // we have no better lurking bound stored; store this lurking bound
-          // and check if it dominates one that is already stored
-          auto it = lurkingBounds.emplace_hint(pos, requiredCutoffBound,
-                                               lurkingBound);
-
-          auto i = lurkingBounds.begin();
-          while (i != it) {
-            if (direction * i->second >= direction * lurkingBound) {
-              auto del = i++;
-              lurkingBounds.erase(del);
-            } else {
-              ++i;
+          // Invariant (kept by this function and by the prefix erases in
+          // propagateRootRedcost): along increasing key, direction*bound is
+          // strictly increasing. Hence the entry at pos carries the smallest
+          // bound among keys >= requiredCutoffBound, and the entries dominated
+          // by the new one form a contiguous block right before pos. This
+          // replaces the former O(n) scans per insertion (O(n^2) per column
+          // with up to 1024 steps) by O(log n + #erased); result is identical.
+          // forward-moving position: requiredCutoffBound is increasing in
+          // the loop, so lower_bound only ever advances.
+          auto& vec = lurkingBounds.v;
+          while (posIdx < vec.size() && vec[posIdx].first < requiredCutoffBound)
+            ++posIdx;
+          if (posIdx < vec.size() &&
+              !(direction * vec[posIdx].second >=
+                direction * (lurkingBound + step))) {
+            double inc = double(step) * redCost;
+            if (inc > 0) {
+              double k =
+                  std::floor((vec[posIdx].first - requiredCutoffBound) / inc);
+              if (k >= 2) {
+                double maxk = double(direction * (lastBound - lurkingBound)) /
+                              double(direction * step);
+                if (k - 1 > maxk) k = maxk + 1;
+                if (k >= 2) lurkingBound += HighsInt(k - 1) * step;
+              }
             }
+            continue;
+          }
+
+          // dominated entries form a contiguous block right before posIdx
+          size_t q = posIdx;
+          while (q > 0 && direction * vec[q - 1].second >=
+                              direction * lurkingBound)
+            --q;
+          if (q < posIdx) {
+            vec[q] = std::make_pair(requiredCutoffBound, lurkingBound);
+            vec.erase(vec.begin() + q + 1, vec.begin() + posIdx);
+          } else {
+            vec.insert(vec.begin() + q,
+                       std::make_pair(requiredCutoffBound, lurkingBound));
           }
+          posIdx = q + 1;
         }
       };