Merge pull request #2695 from fwesselm/enumerateSolutions4

Solution enumeration

Author: jajhall

highs/lp_data/HConst.h

@@ -274,7 +274,8 @@ enum PresolveRuleType : int {
   kPresolveRuleParallelRowsAndCols,
   kPresolveRuleSparsify,
   kPresolveRuleProbing,
-  kPresolveRuleMax = kPresolveRuleProbing,
+  kPresolveRuleEnumeration,
+  kPresolveRuleMax = kPresolveRuleEnumeration,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };

highs/lp_data/HighsModelUtils.cpp

@@ -1455,6 +1455,8 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) {
     return "Sparsify";
   } else if (rule_type == kPresolveRuleProbing) {
     return "Probing";
+  } else if (rule_type == kPresolveRuleEnumeration) {
+    return "Enumeration";
   }
   assert(1 == 0);
   return "????";

highs/mip/HighsDomain.cpp

@@ -1548,10 +1548,7 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
       if (recordRedundantRows_ &&
           mip->row_lower_[mip->a_matrix_.index_[i]] != -kHighsInf &&
           mip->row_upper_[mip->a_matrix_.index_[i]] == kHighsInf)
-        updateRedundantRows(mip->a_matrix_.index_[i], HighsInt{1},
-                            activitymininf_[mip->a_matrix_.index_[i]],
-                            activitymin_[mip->a_matrix_.index_[i]],
-                            mip->row_lower_[mip->a_matrix_.index_[i]]);
+        updateRedundantRows(mip->a_matrix_.index_[i]);
 
       if (deltamin <= 0) {
         updateThresholdLbChange(col, newbound, mip->a_matrix_.value_[i],
@@ -1600,10 +1597,7 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound,
       if (recordRedundantRows_ &&
           mip->row_lower_[mip->a_matrix_.index_[i]] == -kHighsInf &&
           mip->row_upper_[mip->a_matrix_.index_[i]] != kHighsInf)
-        updateRedundantRows(mip->a_matrix_.index_[i], HighsInt{-1},
-                            activitymaxinf_[mip->a_matrix_.index_[i]],
-                            activitymax_[mip->a_matrix_.index_[i]],
-                            mip->row_upper_[mip->a_matrix_.index_[i]]);
+        updateRedundantRows(mip->a_matrix_.index_[i]);
 
       if (deltamax >= 0) {
         updateThresholdLbChange(col, newbound, mip->a_matrix_.value_[i],
@@ -1703,10 +1697,7 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
       if (recordRedundantRows_ &&
           mip->row_lower_[mip->a_matrix_.index_[i]] == -kHighsInf &&
           mip->row_upper_[mip->a_matrix_.index_[i]] != kHighsInf)
-        updateRedundantRows(mip->a_matrix_.index_[i], HighsInt{-1},
-                            activitymaxinf_[mip->a_matrix_.index_[i]],
-                            activitymax_[mip->a_matrix_.index_[i]],
-                            mip->row_upper_[mip->a_matrix_.index_[i]]);
+        updateRedundantRows(mip->a_matrix_.index_[i]);
 
       if (deltamax >= 0) {
         updateThresholdUbChange(col, newbound, mip->a_matrix_.value_[i],
@@ -1758,10 +1749,7 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound,
       if (recordRedundantRows_ &&
           mip->row_lower_[mip->a_matrix_.index_[i]] != -kHighsInf &&
           mip->row_upper_[mip->a_matrix_.index_[i]] == kHighsInf)
-        updateRedundantRows(mip->a_matrix_.index_[i], HighsInt{1},
-                            activitymininf_[mip->a_matrix_.index_[i]],
-                            activitymin_[mip->a_matrix_.index_[i]],
-                            mip->row_lower_[mip->a_matrix_.index_[i]]);
+        updateRedundantRows(mip->a_matrix_.index_[i]);
 
       if (deltamin <= 0) {
         updateThresholdUbChange(col, newbound, mip->a_matrix_.value_[i],
@@ -1845,11 +1833,8 @@ void HighsDomain::recomputeCapacityThreshold(HighsInt row) {
   }
 }
 
-void HighsDomain::updateRedundantRows(HighsInt row, HighsInt direction,
-                                      HighsInt numinf, HighsCDouble activity,
-                                      double bound) {
-  if (numinf != 0 || direction * activity <=
-                         direction * bound + mipsolver->mipdata_->feastol) {
+void HighsDomain::updateRedundantRows(HighsInt row) {
+  if (!isRedundantRow(row)) {
     // row that was found to be redundant should not be non-redundant
     assert(redundantRows_.find(row) == nullptr);
     return;
@@ -1870,6 +1855,12 @@ double HighsDomain::getRedundantRowValue(HighsInt row) const {
   }
 }
 
+bool HighsDomain::isRedundantRow(HighsInt row) const {
+  return (
+      getMinActivity(row) >= mipsolver->model_->row_lower_[row] - feastol() &&
+      getMaxActivity(row) <= mipsolver->model_->row_upper_[row] + feastol());
+}
+
 void HighsDomain::markPropagateCut(Reason reason) {
   switch (reason.type) {
     case Reason::kUnknown:
@@ -1957,7 +1948,7 @@ double HighsDomain::doChangeBound(const HighsDomainChange& boundchg) {
       if (!infeasible_)
         updateActivityLbChange(boundchg.column, oldbound, boundchg.boundval);
 
-      if (!changedcolsflags_[boundchg.column]) {
+      if (!isChangedCol(boundchg.column)) {
         changedcolsflags_[boundchg.column] = 1;
         changedcols_.push_back(boundchg.column);
       }
@@ -1969,7 +1960,7 @@ double HighsDomain::doChangeBound(const HighsDomainChange& boundchg) {
       if (!infeasible_)
         updateActivityUbChange(boundchg.column, oldbound, boundchg.boundval);
 
-      if (!changedcolsflags_[boundchg.column]) {
+      if (!isChangedCol(boundchg.column)) {
         changedcolsflags_[boundchg.column] = 1;
         changedcols_.push_back(boundchg.column);
       }

highs/mip/HighsDomain.h

@@ -324,8 +324,7 @@ class HighsDomain {
 
   void recomputeCapacityThreshold(HighsInt row);
 
-  void updateRedundantRows(HighsInt row, HighsInt direction, HighsInt numInf,
-                           HighsCDouble activity, double bound);
+  void updateRedundantRows(HighsInt row);
 
   double doChangeBound(const HighsDomainChange& boundchg);
 
@@ -444,18 +443,19 @@ class HighsDomain {
 
     changedcols_.erase(
         std::remove_if(changedcols_.begin(), changedcols_.end(),
-                       [&](HighsInt i) { return !changedcolsflags_[i]; }),
+                       [&](HighsInt i) { return !isChangedCol(i); }),
         changedcols_.end());
   }
 
-  void clearChangedCols(HighsInt start) {
-    HighsInt end = changedcols_.size();
-    for (HighsInt i = start; i != end; ++i)
+  void clearChangedCols(size_t start) {
+    for (size_t i = start; i != changedcols_.size(); ++i)
       changedcolsflags_[changedcols_[i]] = 0;
 
     changedcols_.resize(start);
   }
 
+  bool isChangedCol(HighsInt col) const { return changedcolsflags_[col] != 0; }
+
   void markPropagate(HighsInt row);
 
   bool isActive(const HighsDomainChange& domchg) const {
@@ -647,6 +647,8 @@ class HighsDomain {
   double getRedundantRowValue(HighsInt row) const;
 
   void setRecordRedundantRows(bool val) { recordRedundantRows_ = val; };
+
+  bool isRedundantRow(HighsInt row) const;
 };
 
 #endif

highs/mip/HighsImplications.cpp

@@ -300,9 +300,7 @@ bool HighsImplications::runProbing(HighsInt col, HighsInt& numReductions) {
   if (globaldomain.isBinary(col) && !implicationsCached(col, 1) &&
       !implicationsCached(col, 0) &&
       mipsolver.mipdata_->cliquetable.getSubstitution(col) == nullptr) {
-    bool infeasible;
-
-    infeasible = computeImplications(col, 1);
+    bool infeasible = computeImplications(col, 1);
     if (globaldomain.infeasible()) return true;
     if (infeasible) return true;
     if (mipsolver.mipdata_->cliquetable.getSubstitution(col) != nullptr)

highs/mip/HighsMipSolverData.cpp

@@ -98,14 +98,14 @@ std::string HighsMipSolverData::solutionSourceToString(
 
 bool HighsMipSolverData::checkSolution(
     const std::vector<double>& solution) const {
-  for (HighsInt i = 0; i != mipsolver.model_->num_col_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
     if (solution[i] < mipsolver.model_->col_lower_[i] - feastol) return false;
     if (solution[i] > mipsolver.model_->col_upper_[i] + feastol) return false;
     if (mipsolver.isColInteger(i) && fractionality(solution[i]) > feastol)
       return false;
   }
 
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     double rowactivity = 0.0;
 
     HighsInt start = ARstart_[i];
@@ -125,7 +125,7 @@ std::vector<std::tuple<HighsInt, HighsInt, double>>
 HighsMipSolverData::getInfeasibleRows(
     const std::vector<double>& solution) const {
   std::vector<std::tuple<HighsInt, HighsInt, double>> infeasibleRows;
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     HighsInt start = ARstart_[i];
     HighsInt end = ARstart_[i + 1];
 
@@ -149,11 +149,11 @@ HighsMipSolverData::getInfeasibleRows(
 
 bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
                                      const int solution_source) {
-  if (int(solution.size()) != mipsolver.model_->num_col_) return false;
+  if (int(solution.size()) != mipsolver.numCol()) return false;
 
   HighsCDouble obj = 0;
 
-  for (HighsInt i = 0; i != mipsolver.model_->num_col_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
     if (solution[i] < mipsolver.model_->col_lower_[i] - feastol) return false;
     if (solution[i] > mipsolver.model_->col_upper_[i] + feastol) return false;
     if (mipsolver.isColInteger(i) && fractionality(solution[i]) > feastol)
@@ -162,7 +162,7 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
     obj += mipsolver.colCost(i) * solution[i];
   }
 
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     double rowactivity = 0.0;
 
     HighsInt start = ARstart_[i];
@@ -180,7 +180,7 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
 
 bool HighsMipSolverData::solutionRowFeasible(
     const std::vector<double>& solution) const {
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     HighsCDouble c_double_rowactivity = HighsCDouble(0.0);
 
     HighsInt start = ARstart_[i];
@@ -260,7 +260,7 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
   const double feasibility_tolerance =
       mipsolver.options_mip_->mip_feasibility_tolerance;
   // Loop through the trivial heuristics
-  std::vector<double> solution(mipsolver.model_->num_col_);
+  std::vector<double> solution(mipsolver.numCol());
   for (HighsInt try_heuristic = 0; try_heuristic < num_try_heuristic;
        try_heuristic++) {
     if (try_heuristic == 0) {
@@ -271,15 +271,15 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
       if (!all_integer_lower_non_positive) continue;
       // Determine whether a zero row activity is feasible
       bool heuristic_failed = false;
-      for (HighsInt iRow = 0; iRow < mipsolver.model_->num_row_; iRow++) {
+      for (HighsInt iRow = 0; iRow < mipsolver.numRow(); iRow++) {
         if (row_lower[iRow] > feasibility_tolerance ||
             row_upper[iRow] < -feasibility_tolerance) {
           heuristic_failed = true;
           break;
         }
       }
       if (heuristic_failed) continue;
-      solution.assign(mipsolver.model_->num_col_, 0);
+      solution.assign(mipsolver.numCol(), 0);
     } else if (try_heuristic == 1) {
       // Second heuristic is to see whether all-lower for integer
       // variables (if distinct from all-zero) is feasible
@@ -320,7 +320,7 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
     }
 
     HighsCDouble cdouble_obj = 0.0;
-    for (HighsInt iCol = 0; iCol < mipsolver.model_->num_col_; iCol++)
+    for (HighsInt iCol = 0; iCol < mipsolver.numCol(); iCol++)
       cdouble_obj += mipsolver.colCost(iCol) * solution[iCol];
     double obj = double(cdouble_obj);
     const double save_upper_bound = upper_bound;
@@ -692,8 +692,7 @@ void HighsMipSolverData::removeFixedIndices() {
 }
 
 void HighsMipSolverData::init() {
-  postSolveStack.initializeIndexMaps(mipsolver.model_->num_row_,
-                                     mipsolver.model_->num_col_);
+  postSolveStack.initializeIndexMaps(mipsolver.numRow(), mipsolver.numCol());
   mipsolver.orig_model_ = mipsolver.model_;
   feastol = mipsolver.options_mip_->mip_feasibility_tolerance;
   epsilon = mipsolver.options_mip_->small_matrix_value;
@@ -887,11 +886,11 @@ void HighsMipSolverData::runSetup() {
     }
   }
 
-  rowintegral.resize(mipsolver.model_->num_row_);
+  rowintegral.resize(mipsolver.numRow());
 
   // compute the maximal absolute coefficients to filter propagation
-  maxAbsRowCoef.resize(mipsolver.model_->num_row_);
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  maxAbsRowCoef.resize(mipsolver.numRow());
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     double maxabsval = 0.0;
 
     HighsInt start = ARstart_[i];
@@ -1074,7 +1073,7 @@ void HighsMipSolverData::runSetup() {
         debugSolution.debugOrigSolution);
     debugSolution.debugSolObjective = 0;
     HighsCDouble debugsolobj = 0.0;
-    for (HighsInt i = 0; i != mipsolver.model_->num_col_; ++i)
+    for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
       debugsolobj +=
           mipsolver.colCost(i) * HighsCDouble(debugSolution.debugSolution[i]);
     debugSolution.debugSolObjective = static_cast<double>(debugsolobj);
@@ -1301,7 +1300,7 @@ void HighsMipSolverData::performRestart() {
     root_basis.valid = true;
     root_basis.useful = true;
 
-    for (HighsInt i = 0; i < mipsolver.model_->num_col_; ++i)
+    for (HighsInt i = 0; i < mipsolver.numCol(); ++i)
       root_basis.col_status[postSolveStack.getOrigColIndex(i)] =
           basis.col_status[i];
 
@@ -2570,8 +2569,8 @@ void HighsMipSolverData::setupDomainPropagation() {
   pseudocost = HighsPseudocost(mipsolver);
 
   // compute the maximal absolute coefficients to filter propagation
-  maxAbsRowCoef.resize(mipsolver.model_->num_row_);
-  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+  maxAbsRowCoef.resize(mipsolver.numRow());
+  for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
     double maxabsval = 0.0;
 
     HighsInt start = ARstart_[i];

highs/mip/HighsPrimalHeuristics.cpp

@@ -1475,8 +1475,8 @@ void HighsPrimalHeuristics::feasibilityPump() {
   std::vector<double> fracintcost;
   std::vector<HighsInt> fracintset;
 
-  std::vector<HighsInt> mask(mipsolver.model_->num_col_, 1);
-  std::vector<double> cost(mipsolver.model_->num_col_, 0.0);
+  std::vector<HighsInt> mask(mipsolver.numCol(), 1);
+  std::vector<double> cost(mipsolver.numCol(), 0.0);
 
   lprelax.getLpSolver().setOptionValue("simplex_strategy",
                                        kSimplexStrategyPrimal);

highs/mip/HighsSearch.cpp

@@ -412,7 +412,7 @@ HighsInt HighsSearch::selectBranchingCandidate(int64_t maxSbIters,
 
     auto analyzeSolution = [&](double objdelta,
                                const std::vector<double>& sol) {
-      HighsInt numChangedCols = localdom.getChangedCols().size();
+      size_t numChangedCols = localdom.getChangedCols().size();
       HighsInt domchgStackSize = localdom.getDomainChangeStack().size();
       const auto& domchgstack = localdom.getDomainChangeStack();
 
@@ -1517,7 +1517,7 @@ bool HighsSearch::backtrack(bool recoverBasis) {
         // repropagate the node, as it may have become infeasible due to
         // conflicts
         HighsInt oldNumDomchgs = localdom.getNumDomainChanges();
-        HighsInt oldNumChangedCols = localdom.getChangedCols().size();
+        size_t oldNumChangedCols = localdom.getChangedCols().size();
         localdom.propagate();
         if (!localdom.infeasible() &&
             oldNumDomchgs != localdom.getNumDomainChanges()) {
@@ -1564,7 +1564,7 @@ bool HighsSearch::backtrack(bool recoverBasis) {
     if (fallbackbranch)
       currnode.branching_point = currnode.branchingdecision.boundval;
 
-    HighsInt numChangedCols = localdom.getChangedCols().size();
+    size_t numChangedCols = localdom.getChangedCols().size();
     bool passStabilizerToChildNode =
         orbitsValidInChildNode(currnode.branchingdecision);
     localdom.changeBound(currnode.branchingdecision);
@@ -1694,7 +1694,7 @@ bool HighsSearch::backtrackPlunge(HighsNodeQueue& nodequeue) {
       currnode.branching_point = currnode.branchingdecision.boundval;
 
     HighsInt domchgPos = domchgstack.size();
-    HighsInt numChangedCols = localdom.getChangedCols().size();
+    size_t numChangedCols = localdom.getChangedCols().size();
     bool passStabilizerToChildNode =
         orbitsValidInChildNode(currnode.branchingdecision);
     localdom.changeBound(currnode.branchingdecision);

highs/mip/MipTimer.h

@@ -28,6 +28,7 @@ enum iClockMip {
   kMipClockSearch,
   // Search
   kMipClockProbingPresolve,
+  kMipClockEnumerationPresolve,
   kMipClockPerformAging1,
   kMipClockDive,
   kMipClockOpenNodesToQueue0,
@@ -237,6 +238,9 @@ class MipTimer {
     clock[kMipClockProbingPresolve] =
         timer_pointer->clock_def("Probing - presolve");
 
+    clock[kMipClockEnumerationPresolve] =
+        timer_pointer->clock_def("Enumeration - presolve");
+
     // Search - Should correspond to kMipClockSearch
     clock[kMipClockPerformAging1] = timer_pointer->clock_def("Perform aging 1");
     clock[kMipClockDive] = timer_pointer->clock_def("Dive");
@@ -391,7 +395,8 @@ class MipTimer {
   };
 
   void reportMipPresolveClock(const HighsTimerClock& mip_timer_clock) {
-    const std::vector<HighsInt> mip_clock_list{kMipClockProbingPresolve};
+    const std::vector<HighsInt> mip_clock_list{kMipClockProbingPresolve,
+                                               kMipClockEnumerationPresolve};
     reportMipClockList("MipPrslv", mip_clock_list, mip_timer_clock,
                        kMipClockRunPresolve, tolerance_percent_report);
   };