Merge pull request #2490 from ERGO-Code/debugsol

Expand MIP debug solution capabilities

Author: jajhall

highs/mip/HighsDebugSol.cpp

@@ -35,30 +35,44 @@ void HighsDebugSol::activate() {
     if (file) {
       std::string varname;
       double varval;
+      std::string line;
+      bool incolsection;
       std::map<std::string, int> nametoidx;
 
-      for (HighsInt i = 0; i != mipsolver->model_->num_col_; ++i)
-        nametoidx["c" + std::to_string(i)] = i;
+      for (HighsInt i = 0; i != mipsolver->orig_model_->num_col_; ++i)
+        nametoidx[mipsolver->orig_model_->col_names_[i]] = i;
 
-      debugSolution.resize(mipsolver->model_->num_col_, 0.0);
-      while (!file.eof()) {
-        file >> varname;
+      debugOrigSolution.resize(mipsolver->orig_model_->num_col_, 0.0);
+      while (std::getline(file, line)) {
+        // Check for start and stop markers
+        if (line.find("# Columns") != std::string::npos) {
+          incolsection = true;
+          continue;
+        }
+        if (line.find("# Rows") != std::string::npos) {
+          break;
+        }
+        if (!incolsection) continue;
+
+        std::istringstream linestream(line);
+        linestream >> varname;
         auto it = nametoidx.find(varname);
         if (it != nametoidx.end()) {
-          file >> varval;
+          linestream >> varval;
+          debugOrigSolution[it->second] = varval;
           highsLogDev(mipsolver->options_mip_->log_options, HighsLogType::kInfo,
                       "%s = %g\n", varname.c_str(), varval);
-          debugSolution[it->second] = varval;
+          debugOrigSolution[it->second] = varval;
         }
-
-        file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
       }
+      debugSolution = debugOrigSolution;
 
       HighsCDouble debugsolobj = 0.0;
-      for (HighsInt i = 0; i != mipsolver->model_->num_col_; ++i)
-        debugsolobj += mipsolver->model_->col_cost_[i] * debugSolution[i];
+      for (HighsInt i = 0; i != mipsolver->orig_model_->num_col_; ++i)
+        debugsolobj += mipsolver->orig_model_->col_cost_[i] *
+                       HighsCDouble(debugOrigSolution[i]);
 
-      debugSolObjective = double(debugsolobj);
+      debugSolObjective = double(debugsolobj + mipsolver->orig_model_->offset_);
       debugSolActive = true;
       printf("debug sol active\n");
       registerDomain(mipsolver->mipdata_->domain);

highs/mip/HighsDebugSol.h

@@ -29,6 +29,7 @@ class HighsLp;
 struct HighsDebugSol {
   const HighsMipSolver* mipsolver;
   double debugSolObjective;
+  std::vector<double> debugOrigSolution;
   std::vector<double> debugSolution;
   bool debugSolActive;
   std::unordered_map<const HighsDomain*, std::multiset<HighsDomainChange>>

highs/mip/HighsMipSolver.cpp

@@ -87,6 +87,11 @@ void HighsMipSolver::run() {
   analysis_.mipTimerStart(kMipClockInit);
   mipdata_->init();
   analysis_.mipTimerStop(kMipClockInit);
+#ifdef HIGHS_DEBUGSOL
+  mipdata_->debugSolution.activate();
+  bool debugSolActive = false;
+  std::swap(mipdata_->debugSolution.debugSolActive, debugSolActive);
+#endif
   analysis_.mipTimerStart(kMipClockRunPresolve);
   mipdata_->runPresolve(options_mip_->presolve_reduction_limit);
   analysis_.mipTimerStop(kMipClockRunPresolve);
@@ -119,6 +124,9 @@ void HighsMipSolver::run() {
     highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                  "MIP-Timing: %11.2g - starting  setup\n", timer_.read());
   analysis_.mipTimerStart(kMipClockRunSetup);
+#ifdef HIGHS_DEBUGSOL
+  mipdata_->debugSolution.debugSolActive = debugSolActive;
+#endif
   mipdata_->runSetup();
   analysis_.mipTimerStop(kMipClockRunSetup);
   if (analysis_.analyse_mip_time && !submip)

highs/mip/HighsMipSolverData.cpp

@@ -694,11 +694,6 @@ void HighsMipSolverData::init() {
 }
 
 void HighsMipSolverData::runPresolve(const HighsInt presolve_reduction_limit) {
-#ifdef HIGHS_DEBUGSOL
-  bool debugSolActive = false;
-  std::swap(debugSolution.debugSolActive, debugSolActive);
-#endif
-
   mipsolver.timer_.start(mipsolver.timer_.presolve_clock);
   presolve::HPresolve presolve;
   if (!presolve.okSetInput(mipsolver, presolve_reduction_limit)) {
@@ -709,13 +704,6 @@ void HighsMipSolverData::runPresolve(const HighsInt presolve_reduction_limit) {
     presolve_status = presolve.getPresolveStatus();
   }
   mipsolver.timer_.stop(mipsolver.timer_.presolve_clock);
-
-#ifdef HIGHS_DEBUGSOL
-  debugSolution.debugSolActive = debugSolActive;
-  if (debugSolution.debugSolActive) debugSolution.registerDomain(domain);
-  assert(!debugSolution.debugSolActive ||
-         checkSolution(debugSolution.debugSolution));
-#endif
 }
 
 void HighsMipSolverData::runSetup() {
@@ -1022,10 +1010,18 @@ void HighsMipSolverData::runSetup() {
   heuristics.setupIntCols();
 
 #ifdef HIGHS_DEBUGSOL
-  if (numRestarts == 0) {
-    debugSolution.activate();
-    assert(!debugSolution.debugSolActive ||
-           checkSolution(debugSolution.debugSolution));
+  if (debugSolution.debugSolActive) {
+    debugSolution.debugSolution.clear();
+    debugSolution.debugSolution = postSolveStack.getReducedPrimalSolution(
+        debugSolution.debugOrigSolution);
+    debugSolution.debugSolObjective = 0;
+    HighsCDouble debugsolobj = 0.0;
+    for (HighsInt i = 0; i != mipsolver.model_->num_col_; ++i)
+      debugsolobj +=
+          mipsolver.colCost(i) * HighsCDouble(debugSolution.debugSolution[i]);
+    debugSolution.debugSolObjective = static_cast<double>(debugsolobj);
+    debugSolution.registerDomain(domain);
+    assert(checkSolution(debugSolution.debugSolution));
   }
 #endif
 
@@ -1218,6 +1214,10 @@ void HighsMipSolverData::performRestart() {
   presolvedModel = lp.getLp();
   presolvedModel.offset_ = offset;
   presolvedModel.integrality_ = std::move(integrality);
+#ifdef HIGHS_DEBUGSOL
+  bool debugSolActive = false;
+  std::swap(debugSolution.debugSolActive, debugSolActive);
+#endif
 
   const HighsBasis& basis = firstrootbasis;
   if (basis.valid) {
@@ -1321,6 +1321,9 @@ void HighsMipSolverData::performRestart() {
   }
   // Bounds are currently in the original space since presolve will have
   // changed offset_
+#ifdef HIGHS_DEBUGSOL
+  debugSolution.debugSolActive = debugSolActive;
+#endif
   runSetup();
 
   postSolveStack.removeCutsFromModel(numCuts);