Merge pull request #2027 from ERGO-Code/fix-1537

Added trivial heuristics

check/TestMipSolver.cpp

@@ -677,8 +677,17 @@ TEST_CASE("IP-infeasible-unbounded", "[highs_test_mip_solver]") {
       // Solve
       highs.passModel(lp);
       highs.run();
-      REQUIRE(highs.getModelStatus() ==
-              HighsModelStatus::kUnboundedOrInfeasible);
+      HighsModelStatus required_model_status;
+      if (k == 0) {
+        if (l == 0) {
+          required_model_status = HighsModelStatus::kInfeasible;
+        } else {
+          required_model_status = HighsModelStatus::kUnbounded;
+        }
+      } else {
+        required_model_status = HighsModelStatus::kUnboundedOrInfeasible;
+      }
+      REQUIRE(highs.getModelStatus() == required_model_status);
     }
     highs.setOptionValue("presolve", kHighsOnString);
   }

src/mip/HighsMipSolver.cpp

@@ -175,6 +175,17 @@ void HighsMipSolver::run() {
       cleanupSolve();
       return;
     }
+    // Apply the trivial heuristics
+    analysis_.mipTimerStart(kMipClockTrivialHeuristics);
+    HighsModelStatus model_status = mipdata_->trivialHeuristics();
+    if (modelstatus_ == HighsModelStatus::kNotset &&
+        model_status == HighsModelStatus::kInfeasible) {
+      // trivialHeuristics can spot trivial infeasibility, so act on it
+      modelstatus_ = model_status;
+      cleanupSolve();
+      return;
+    }
+    analysis_.mipTimerStop(kMipClockTrivialHeuristics);
     if (analysis_.analyse_mip_time & !submip)
       highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                    "MIP-Timing: %11.2g - starting evaluate root node\n",

src/mip/HighsMipSolverData.cpp

@@ -60,6 +60,18 @@ std::string HighsMipSolverData::solutionSourceToString(
   } else if (solution_source == kSolutionSourceUnbounded) {
     if (code) return "U";
     return "Unbounded";
+  } else if (solution_source == kSolutionSourceTrivialZ) {
+    if (code) return "z";
+    return "Trivial zero";
+  } else if (solution_source == kSolutionSourceTrivialL) {
+    if (code) return "l";
+    return "Trivial lower";
+  } else if (solution_source == kSolutionSourceTrivialU) {
+    if (code) return "u";
+    return "Trivial upper";
+  } else if (solution_source == kSolutionSourceTrivialP) {
+    if (code) return "p";
+    return "Trivial point";
     //  } else if (solution_source == kSolutionSourceOpt1) {
     //    if (code) return "1";
     //    return "1-opt";
@@ -136,6 +148,155 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
   return addIncumbent(solution, double(obj), solution_source);
 }
 
+bool HighsMipSolverData::solutionRowFeasible(
+    const std::vector<double>& solution) const {
+  for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) {
+    HighsCDouble c_double_rowactivity = HighsCDouble(0.0);
+
+    HighsInt start = ARstart_[i];
+    HighsInt end = ARstart_[i + 1];
+
+    for (HighsInt j = start; j != end; ++j)
+      c_double_rowactivity += HighsCDouble(solution[ARindex_[j]] * ARvalue_[j]);
+
+    double rowactivity = double(c_double_rowactivity);
+    if (rowactivity > mipsolver.rowUpper(i) + feastol) return false;
+    if (rowactivity < mipsolver.rowLower(i) - feastol) return false;
+  }
+  return true;
+}
+
+HighsModelStatus HighsMipSolverData::trivialHeuristics() {
+  //  printf("\nHighsMipSolverData::trivialHeuristics() Number of continuous
+  //  columns is %d\n",
+  //	 int(continuous_cols.size()));
+  if (continuous_cols.size() > 0) return HighsModelStatus::kNotset;
+  const HighsInt num_try_heuristic = 4;
+  const std::vector<int> heuristic_source = {
+      kSolutionSourceTrivialZ, kSolutionSourceTrivialL, kSolutionSourceTrivialU,
+      kSolutionSourceTrivialP};
+
+  std::vector<double> col_lower = mipsolver.model_->col_lower_;
+  std::vector<double> col_upper = mipsolver.model_->col_upper_;
+  const std::vector<double>& row_lower = mipsolver.model_->row_lower_;
+  const std::vector<double>& row_upper = mipsolver.model_->row_upper_;
+  const HighsSparseMatrix& matrix = mipsolver.model_->a_matrix_;
+  // Determine the following properties, according to which some
+  // trivial heuristics are duplicated or fail immediately
+  bool all_integer_lower_non_positive = true;
+  bool all_integer_lower_zero = true;
+  bool all_integer_lower_finite = true;
+  bool all_integer_upper_finite = true;
+  for (HighsInt integer_col = 0; integer_col < numintegercols; integer_col++) {
+    HighsInt iCol = integer_cols[integer_col];
+    // Round bounds in to nearest integer
+    col_lower[iCol] = std::ceil(col_lower[iCol]);
+    col_upper[iCol] = std::floor(col_upper[iCol]);
+    // If bounds are inconsistent then MIP is infeasible
+    if (col_lower[iCol] > col_upper[iCol]) return HighsModelStatus::kInfeasible;
+
+    if (col_lower[iCol] > 0) all_integer_lower_non_positive = false;
+    if (col_lower[iCol]) all_integer_lower_zero = false;
+    if (col_lower[iCol] <= -kHighsInf) all_integer_lower_finite = false;
+    if (col_upper[iCol] >= kHighsInf) all_integer_upper_finite = false;
+    // Only continue if one of the properties still holds
+    if (!(all_integer_lower_non_positive || all_integer_lower_zero ||
+          all_integer_upper_finite))
+      break;
+  }
+  const bool all_integer_boxed =
+      all_integer_lower_finite && all_integer_upper_finite;
+  //  printf(
+  //      "Trying trivial heuristics\n"
+  //      "   all_integer_lower_non_positive = %d\n"
+  //      "   all_integer_lower_zero = %d\n"
+  //      "   all_integer_upper_finite = %d\n"
+  //      "   all_integer_boxed = %d\n",
+  //      all_integer_lower_non_positive, all_integer_lower_zero,
+  //      all_integer_upper_finite, all_integer_boxed);
+  const double feasibility_tolerance =
+      mipsolver.options_mip_->mip_feasibility_tolerance;
+  // Loop through the trivial heuristics
+  std::vector<double> solution(mipsolver.model_->num_col_);
+  for (HighsInt try_heuristic = 0; try_heuristic < num_try_heuristic;
+       try_heuristic++) {
+    if (try_heuristic == 0) {
+      // First heuristic is to see whether all-zero for integer
+      // variables is feasible
+      //
+      // If there is a positive lower bound then the heuristic fails
+      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++) {
+        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);
+    } else if (try_heuristic == 1) {
+      // Second heuristic is to see whether all-lower for integer
+      // variables (if distinct from all-zero) is feasible
+      if (all_integer_lower_zero) continue;
+      // Trivially feasible for columns
+      if (!solutionRowFeasible(col_lower)) continue;
+      solution = col_lower;
+    } else if (try_heuristic == 2) {
+      // Third heuristic is to see whether all-upper for integer
+      // variables is feasible
+      //
+      // If there is an infinite upper bound then the heuristic fails
+      if (!all_integer_upper_finite) continue;
+      // Trivially feasible for columns
+      if (!solutionRowFeasible(col_upper)) continue;
+      solution = col_upper;
+    } else if (try_heuristic == 3) {
+      // Fourth heuristic is to see whether the "lock point" is feasible
+      if (!all_integer_boxed) continue;
+      for (HighsInt integer_col = 0; integer_col < numintegercols;
+           integer_col++) {
+        HighsInt iCol = integer_cols[integer_col];
+        HighsInt num_positive_values = 0;
+        HighsInt num_negative_values = 0;
+        for (HighsInt iEl = matrix.start_[iCol]; iEl < matrix.start_[iCol + 1];
+             iEl++) {
+          if (matrix.value_[iEl] > 0)
+            num_positive_values++;
+          else
+            num_negative_values++;
+        }
+        solution[iCol] = num_positive_values > num_negative_values
+                             ? col_lower[iCol]
+                             : col_upper[iCol];
+      }
+      // Trivially feasible for columns
+      if (!solutionRowFeasible(solution)) continue;
+    }
+
+    HighsCDouble cdouble_obj = 0.0;
+    for (HighsInt iCol = 0; iCol < mipsolver.model_->num_col_; iCol++)
+      cdouble_obj += mipsolver.colCost(iCol) * solution[iCol];
+    double obj = double(cdouble_obj);
+    const double save_upper_bound = upper_bound;
+    const bool new_incumbent =
+        addIncumbent(solution, obj, heuristic_source[try_heuristic]);
+    const bool lc_report = true;
+    if (lc_report) {
+      printf("Trivial heuristic %d has succeeded: objective = %g",
+             int(try_heuristic), obj);
+      if (new_incumbent) {
+        printf("; upper bound from %g to %g\n", save_upper_bound, upper_bound);
+      } else {
+        printf("\n");
+      }
+    }
+  }
+  return HighsModelStatus::kNotset;
+}
+
 void HighsMipSolverData::startAnalyticCenterComputation(
     const highs::parallel::TaskGroup& taskGroup) {
   taskGroup.spawn([&]() {
@@ -1312,9 +1473,12 @@ void HighsMipSolverData::printSolutionSourceKey() {
   // not a solution source, but used to force the last logging line to
   // be printed
   const int last_enum = kSolutionSourceCount - 1;
-  int half_list = last_enum / 2;
+  int third_list = 5;
+  int two_third_list = 10;
+  std::vector<int> limits = {third_list, two_third_list, last_enum};
+
   ss.str(std::string());
-  for (int k = 0; k < half_list; k++) {
+  for (int k = 0; k < limits[0]; k++) {
     if (k == 0) {
       ss << "\nSrc: ";
     } else {
@@ -1323,21 +1487,22 @@ void HighsMipSolverData::printSolutionSourceKey() {
     ss << solutionSourceToString(k) << " => "
        << solutionSourceToString(k, false);
   }
-  highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, "%s\n",
-               ss.str().c_str());
-
-  ss.str(std::string());
-  for (int k = half_list; k < last_enum; k++) {
-    if (k == half_list) {
-      ss << "     ";
-    } else {
-      ss << "; ";
+  highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo,
+               "%s;\n", ss.str().c_str());
+  for (int line = 0; line < 2; line++) {
+    ss.str(std::string());
+    for (int k = limits[line]; k < limits[line + 1]; k++) {
+      if (k == limits[line]) {
+        ss << "     ";
+      } else {
+        ss << "; ";
+      }
+      ss << solutionSourceToString(k) << " => "
+         << solutionSourceToString(k, false);
     }
-    ss << solutionSourceToString(k) << " => "
-       << solutionSourceToString(k, false);
+    highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo,
+                 "%s%s\n", ss.str().c_str(), line == 0 ? ";" : "");
   }
-  highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, "%s\n",
-               ss.str().c_str());
 }
 
 void HighsMipSolverData::printDisplayLine(const int solution_source) {

src/mip/HighsMipSolverData.h

@@ -47,6 +47,10 @@ enum MipSolutionSource : int {
   kSolutionSourceSolveLp,
   kSolutionSourceEvaluateNode,
   kSolutionSourceUnbounded,
+  kSolutionSourceTrivialZ,
+  kSolutionSourceTrivialL,
+  kSolutionSourceTrivialU,
+  kSolutionSourceTrivialP,
   //  kSolutionSourceOpt1,
   //  kSolutionSourceOpt2,
   kSolutionSourceCleanup,
@@ -199,6 +203,9 @@ struct HighsMipSolverData {
     domain.addConflictPool(conflictPool);
   }
 
+  bool solutionRowFeasible(const std::vector<double>& solution) const;
+  HighsModelStatus trivialHeuristics();
+
   void startAnalyticCenterComputation(
       const highs::parallel::TaskGroup& taskGroup);
   void finishAnalyticCenterComputation(

src/mip/MipTimer.h

@@ -24,6 +24,7 @@ enum iClockMip {
   kMipClockInit,
   kMipClockRunPresolve,
   kMipClockRunSetup,
+  kMipClockTrivialHeuristics,
   kMipClockEvaluateRootNode,
   kMipClockPerformAging0,
   kMipClockSearch,
@@ -95,6 +96,8 @@ class MipTimer {
     clock[kMipClockInit] = timer_pointer->clock_def("Initialise");
     clock[kMipClockRunPresolve] = timer_pointer->clock_def("Run presolve");
     clock[kMipClockRunSetup] = timer_pointer->clock_def("Run setup");
+    clock[kMipClockTrivialHeuristics] =
+        timer_pointer->clock_def("Trivial heuristics");
     clock[kMipClockEvaluateRootNode] =
         timer_pointer->clock_def("Evaluate root node");
     clock[kMipClockPerformAging0] = timer_pointer->clock_def("Perform aging 0");
@@ -249,11 +252,14 @@ class MipTimer {
   };
 
   void reportMipLevel1Clock(const HighsTimerClock& mip_timer_clock) {
-    const std::vector<HighsInt> mip_clock_list{
-        kMipClockInit,          kMipClockRunPresolve,
-        kMipClockRunSetup,      kMipClockEvaluateRootNode,
-        kMipClockPerformAging0, kMipClockSearch,
-        kMipClockPostsolve};
+    const std::vector<HighsInt> mip_clock_list{kMipClockInit,
+                                               kMipClockRunPresolve,
+                                               kMipClockRunSetup,
+                                               kMipClockTrivialHeuristics,
+                                               kMipClockEvaluateRootNode,
+                                               kMipClockPerformAging0,
+                                               kMipClockSearch,
+                                               kMipClockPostsolve};
     reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock,
                        kMipClockTotal, tolerance_percent_report);
   };