Now timing knapsack est/solve and (also) looking for knapsack after presolve; formatted

Author: jajhall

highs/lp_data/HighsLp.cpp

@@ -85,18 +85,57 @@ bool HighsLp::equalButForNames(const HighsLp& lp) const {
 }
 
 bool HighsLp::equalButForScalingAndNames(const HighsLp& lp) const {
+  bool equal;
   bool equal_vectors = true;
-  equal_vectors = this->num_col_ == lp.num_col_ && equal_vectors;
-  equal_vectors = this->num_row_ == lp.num_row_ && equal_vectors;
-  equal_vectors = this->sense_ == lp.sense_ && equal_vectors;
-  equal_vectors = this->offset_ == lp.offset_ && equal_vectors;
-  equal_vectors = this->model_name_ == lp.model_name_ && equal_vectors;
-  equal_vectors = this->col_cost_ == lp.col_cost_ && equal_vectors;
-  equal_vectors = this->col_upper_ == lp.col_upper_ && equal_vectors;
-  equal_vectors = this->col_lower_ == lp.col_lower_ && equal_vectors;
-  equal_vectors = this->row_upper_ == lp.row_upper_ && equal_vectors;
-  equal_vectors = this->row_lower_ == lp.row_lower_ && equal_vectors;
+  equal = this->num_col_ == lp.num_col_;
+  equal_vectors = equal && equal_vectors;
 #ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal col dim\n");
+#endif
+  equal = this->num_row_ == lp.num_row_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal row dim\n");
+#endif
+  equal = this->sense_ == lp.sense_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal sense\n");
+#endif
+  equal = this->offset_ == lp.offset_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal offset\n");
+#endif
+  equal = this->model_name_ == lp.model_name_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal model name\n");
+#endif
+  equal = this->col_cost_ == lp.col_cost_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal col cost\n");
+#endif
+  equal = this->col_lower_ == lp.col_lower_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal col lower\n");
+#endif
+  equal = this->col_upper_ == lp.col_upper_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal col upper\n");
+#endif
+  equal = this->row_lower_ == lp.row_lower_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal row lower\n");
+#endif
+  equal = this->row_upper_ == lp.row_upper_;
+  equal_vectors = equal && equal_vectors;
+#ifndef NDEBUG
+  if (!equal) printf("HighsLp::equalButForNames: unequal row upper\n");
   if (!equal_vectors) printf("HighsLp::equalButForNames: Unequal vectors\n");
 #endif
   const bool equal_matrix = this->a_matrix_ == lp.a_matrix_;

highs/mip/HighsMipSolver.cpp

@@ -88,22 +88,31 @@ void HighsMipSolver::run() {
   mipdata_->init();
   analysis_.mipTimerStop(kMipClockInit);
 
-  // Determine whether this is a knapsack problem and, if so, at least
-  // update the data on knapsack sub-MIPs
-  HighsInt capacity = 0;
+  // If this is a knapsack problem, don't do MIP logging in
+  // cleanupSolve
+  const bool mip_logging_for_knapsack = false;
+
+  // Look for knapsack before presolve
+  analysis_.mipTimerStart(kMipClockKnapsack);
   if (mipdata_->mipIsKnapsack()) {
-    mipdata_->knapsack_data_.num_problem++;
-    mipdata_->knapsack_data_.sum_variables += orig_model_->num_col_;
-    mipdata_->knapsack_data_.sum_capacity += mipdata_->knapsack_capacity_;
     // Solve as a knapsack
     HighsStatus call_status = mipdata_->heuristics.solveMipKnapsack();
     assert(call_status == HighsStatus::kOk);
-    const bool mip_logging = false;
-    cleanupSolve(mip_logging);
+    cleanupSolve(mip_logging_for_knapsack);
+    analysis_.mipTimerStop(kMipClockKnapsack);
     return;
   }
+  analysis_.mipTimerStop(kMipClockKnapsack);
 
+  const HighsInt num_col_before_presolve = model_->num_col_;
+  const HighsInt num_row_before_presolve = model_->num_row_;
   analysis_.mipTimerStart(kMipClockRunPresolve);
+  // Run presolve
+  //
+  // Note that, even with presolve=off and presolve_reduction_limit=0,
+  // reductions can occur. This is because presolve scales the MIP -
+  // and must do - using HighsPresolve::transformColumn. Note that
+  // presolve is where a maximization is converted to a minimization.
   mipdata_->runPresolve(options_mip_->presolve_reduction_limit);
   analysis_.mipTimerStop(kMipClockRunPresolve);
   analysis_.mipTimerStop(kMipClockPresolve);
@@ -128,6 +137,21 @@ void HighsMipSolver::run() {
     cleanupSolve();
     return;
   }
+  // Possibly look for knapsack after presolve
+  const bool reduced = num_col_before_presolve > model_->num_col_ ||
+                       num_row_before_presolve > model_->num_row_;
+  if (reduced) {
+    analysis_.mipTimerStart(kMipClockKnapsack);
+    if (mipdata_->mipIsKnapsack()) {
+      // Solve as a knapsack
+      HighsStatus call_status = mipdata_->heuristics.solveMipKnapsack();
+      assert(call_status == HighsStatus::kOk);
+      cleanupSolve(mip_logging_for_knapsack);
+      analysis_.mipTimerStop(kMipClockKnapsack);
+      return;
+    }
+    analysis_.mipTimerStop(kMipClockKnapsack);
+  }
 
   analysis_.mipTimerStart(kMipClockSolve);
 

highs/mip/HighsMipSolverData.cpp

@@ -2653,7 +2653,9 @@ void HighsMipSolverData::callbackUserSolution(
   }
 }
 
-bool HighsMipSolverData::mipIsKnapsack(const bool logging) {
+bool HighsMipSolverData::mipIsKnapsack(const bool silent) {
+  // Silent is to prevent duplicate logging and data collection when
+  // using assert(mipIsKnapsack(true));
   const HighsLp& lp = *(mipsolver.model_);
   // Has to have one constraint
   if (lp.num_row_ != 1) return false;
@@ -2671,7 +2673,7 @@ bool HighsMipSolverData::mipIsKnapsack(const bool logging) {
       HighsIntegers::integralScale(lp.a_matrix_.value_, 1e-6, 1e-6);
   if (this->knapsack_integral_scale_ == 0) return false;
   if (this->knapsack_integral_scale_ > 1000000) return false;
-  if (logging)
+  if (!silent)
     highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo,
                  "MIP is a knapsack problem with with scale %d\n",
                  int(this->knapsack_integral_scale_));
@@ -2683,6 +2685,11 @@ bool HighsMipSolverData::mipIsKnapsack(const bool logging) {
   const double capacity_margin = 1e-6;
   this->knapsack_capacity_ = std::floor(double_capacity + capacity_margin);
   // Problem is knapsack!
+  if (!silent) {
+    this->knapsack_data_.num_problem++;
+    this->knapsack_data_.sum_variables += lp.num_col_;
+    this->knapsack_data_.sum_capacity += this->knapsack_capacity_;
+  }
   return true;
 }
 

highs/mip/HighsMipSolverData.h

@@ -313,7 +313,7 @@ struct HighsMipSolverData {
   void callbackUserSolution(
       const double mipsolver_objective_value,
       const userMipSolutionCallbackOrigin user_solution_callback_origin);
-  bool mipIsKnapsack(const bool logging = true);
+  bool mipIsKnapsack(const bool silent = false);
 };
 
 #endif

highs/mip/HighsPrimalHeuristics.cpp

@@ -1786,7 +1786,8 @@ HighsStatus HighsPrimalHeuristics::solveMipKnapsackReturn(
 HighsStatus HighsPrimalHeuristics::solveMipKnapsack() {
   const HighsLp& lp = *(mipsolver.model_);
   const HighsLogOptions& log_options = mipsolver.options_mip_->log_options;
-  assert(mipsolver.mipdata_->mipIsKnapsack(false));
+  // Check silently that this really is a knapsack
+  assert(mipsolver.mipdata_->mipIsKnapsack(true));
 
   const bool upper = lp.row_upper_[0] < kHighsInf;
   const HighsInt constraint_sign = upper ? 1 : -1;

highs/mip/MipTimer.h

@@ -19,6 +19,7 @@ enum iClockMip {
   kMipClockPostsolve,
   // Level 1
   kMipClockInit,
+  kMipClockKnapsack,
   kMipClockRunPresolve,
   kMipClockRunSetup,
   kMipClockFeasibilityJump,
@@ -134,6 +135,7 @@ class MipTimer {
 
     // Level 1 - Should correspond to kMipClockTotal
     clock[kMipClockInit] = timer_pointer->clock_def("Initialise");
+    clock[kMipClockKnapsack] = timer_pointer->clock_def("Knapsack test+solve");
     clock[kMipClockRunPresolve] = timer_pointer->clock_def("Run presolve");
     clock[kMipClockRunSetup] = timer_pointer->clock_def("Run setup");
     clock[kMipClockFeasibilityJump] =
@@ -321,6 +323,7 @@ class MipTimer {
 
   void reportMipLevel1Clock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockInit,
+                                               kMipClockKnapsack,
                                                kMipClockRunPresolve,
                                                kMipClockRunSetup,
                                                kMipClockFeasibilityJump,
@@ -330,7 +333,7 @@ class MipTimer {
                                                kMipClockSearch,
                                                kMipClockPostsolve};
     reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal, tolerance_percent_report);
+                       kMipClockTotal);  //, tolerance_percent_report);
   };
 
   void reportMipSolveLpClock(const HighsTimerClock& mip_timer_clock) {
@@ -441,8 +444,9 @@ class MipTimer {
                    const HighsTimerClock& mip_timer_clock, const bool header,
                    const bool end_line) {
     const std::vector<HighsInt> mip_clock_list{
-        kMipClockRunPresolve, kMipClockEvaluateRootNode,
-        kMipClockDivePrimalHeuristics, kMipClockTheDive, kMipClockNodeSearch};
+        kMipClockKnapsack,         kMipClockRunPresolve,
+        kMipClockEvaluateRootNode, kMipClockDivePrimalHeuristics,
+        kMipClockTheDive,          kMipClockNodeSearch};
     csvMipClockList("csvMIP", model_name, mip_clock_list, mip_timer_clock,
                     kMipClockTotal, header, end_line);
   };

highs/presolve/HPresolve.cpp

@@ -4339,7 +4339,7 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
   analysis_.setup(this->model, this->options, this->numDeletedRows,
                   this->numDeletedCols);
 
-  if (options->presolve != kHighsOffString) {
+  if (options->presolve != kHighsOffString && reductionLimit > 0) {
     if (mipsolver) mipsolver->mipdata_->cliquetable.setPresolveFlag(true);
     if (!mipsolver || mipsolver->mipdata_->numRestarts == 0)
       highsLogUser(options->log_options, HighsLogType::kInfo,
@@ -4529,10 +4529,18 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
     report();
   } else {
     highsLogUser(options->log_options, HighsLogType::kInfo,
-                 "\nPresolve is switched off\n");
+                 "\nPresolve: Model reduction is switched off\n");
   }
 
+  this->num_true_reductions_ = postsolve_stack.numReductions();
   if (mipsolver != nullptr) scaleMIP(postsolve_stack);
+  const size_t num_scaling_reductions =
+      postsolve_stack.numReductions() - this->num_true_reductions_;
+  if (num_scaling_reductions)
+    highsLogUser(options->log_options, HighsLogType::kInfo,
+                 "\nPresolve: MIP scaling requires %" PRId64
+                 " reduction operation%s\n",
+                 num_scaling_reductions, num_scaling_reductions > 1 ? "s" : "");
 
   // analysePresolveRuleLog() should return true - no errors
   assert(analysis_.analysePresolveRuleLog());
@@ -4680,11 +4688,12 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) {
     if (options->presolve != kHighsOffString &&
         reductionLimit < kHighsSize_tInf) {
       highsLogUser(options->log_options, HighsLogType::kInfo,
-                   "Presolve performed %" PRId64 " of %" PRId64
-                   " permitted reductions\n",
-                   postsolve_stack.numReductions(), reductionLimit);
+                   "\nPresolve: Performed %" PRId64 " of %" PRId64
+                   " permitted true reductions\n",
+                   this->num_true_reductions_, reductionLimit);
     }
   };
+  this->num_true_reductions_ = 0;
   switch (presolve(postsolve_stack)) {
     case Result::kStopped:
     case Result::kOk:

highs/presolve/HPresolve.h

@@ -132,6 +132,7 @@ class HPresolve {
   HighsInt oldNumCol;
   HighsInt oldNumRow;
   bool probingEarlyAbort;
+  int64_t num_true_reductions_;
 
   enum class Result {
     kOk,