Merge pull request #2498 from ERGO-Code/no-col-dual

More informative logging for error checking and KKT errors in Highs::postsolve

Author: jajhall

check/TestPresolve.cpp

@@ -735,6 +735,32 @@ TEST_CASE("presolve-issue-2446", "[highs_test_presolve]") {
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kReduced);
 }
 
+TEST_CASE("presolve-solve-postsolve-no-col-dual", "[highs_test_presolve]") {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/afiro.mps";
+  highs.readModel(model_file);
+  highs.presolve();
+  HighsLp presolved_lp = highs.getPresolvedLp();
+  Highs highs1;
+  highs1.setOptionValue("output_flag", dev_run);
+  highs1.setOptionValue("presolve", kHighsOffString);
+  highs1.passModel(presolved_lp);
+  highs1.run();
+  HighsSolution solution = highs1.getSolution();
+
+  // Perform postsolve using the optimal solution and basis for the
+  // presolved model
+  REQUIRE(highs.postsolve(solution) == HighsStatus::kOk);
+
+  // If row duals are supplied, then column duals must also be suppplied
+  solution.col_dual.clear();
+  REQUIRE(highs.postsolve(solution) == HighsStatus::kError);
+
+  highs.resetGlobalScheduler(true);
+}
+
 TEST_CASE("presolve-egout-ac", "[highs_test_presolve]") {
   // Tests the case where, for this model when run_crossover is off,
   // sparsify is used to reduce the LP to empty. However, when
@@ -811,4 +837,4 @@ TEST_CASE("presolve-egout-ac", "[highs_test_presolve]") {
           lp_presolve_requires_basis_postsolve);
 
   h.resetGlobalScheduler(true);
-}
+}

highs/lp_data/Highs.cpp

@@ -4116,8 +4116,10 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
 
   // Must at least have a primal column solution of the right size
   if (HighsInt(solution.col_value.size()) != presolved_lp.num_col_) {
-    highsLogUser(options_.log_options, HighsLogType::kError,
-                 "Primal solution provided to postsolve is incorrect size\n");
+    highsLogUser(
+        options_.log_options, HighsLogType::kError,
+        "Primal solution provided to postsolve is of size %d rather than %d\n",
+        int(solution.col_value.size()), int(presolved_lp.num_col_));
     return HighsStatus::kError;
   }
   // Check any basis that is supplied
@@ -4179,7 +4181,7 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
           "Postsolve performed for MIP, but model status cannot be known\n");
     } else {
       highsLogUser(options_.log_options, HighsLogType::kError,
-                   "Postsolve return status is %d\n", (int)postsolve_status);
+                   "Postsolve return status is %d\n", int(postsolve_status));
       setHighsModelStatusAndClearSolutionAndBasis(
           HighsModelStatus::kPostsolveError);
     }
@@ -4194,10 +4196,20 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
         presolve_.data_.recovered_solution_.row_dual.size() > 0 ||
         presolve_.data_.recovered_solution_.dual_valid;
     if (dual_supplied) {
-      if (!isDualSolutionRightSize(presolved_lp,
-                                   presolve_.data_.recovered_solution_)) {
+      if (!isRowDualSolutionRightSize(presolved_lp,
+                                      presolve_.data_.recovered_solution_)) {
+        highsLogUser(options_.log_options, HighsLogType::kError,
+                     "Row dual solution provided to postsolve is of size %d "
+                     "rather than %d\n",
+                     int(solution.row_dual.size()), int(presolved_lp.num_row_));
+        return HighsStatus::kError;
+      }
+      if (!isColDualSolutionRightSize(presolved_lp,
+                                      presolve_.data_.recovered_solution_)) {
         highsLogUser(options_.log_options, HighsLogType::kError,
-                     "Dual solution provided to postsolve is incorrect size\n");
+                     "Column dual solution provided to postsolve is of size %d "
+                     "rather than %d\n",
+                     int(solution.col_dual.size()), int(presolved_lp.num_col_));
         return HighsStatus::kError;
       }
       presolve_.data_.recovered_solution_.dual_valid = true;
@@ -4288,6 +4300,9 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
         getKktFailures(this->options_, is_qp, this->model_.lp_,
                        this->model_.lp_.col_cost_, this->solution_, this->info_,
                        get_residuals);
+        highsLogUser(options_.log_options, HighsLogType::kInfo, "\n");
+        reportLpKktFailures(this->model_.lp_, this->options_, this->info_,
+                            "After postsolve");
         if (info_.num_primal_infeasibilities == 0 &&
             info_.num_dual_infeasibilities == 0) {
           model_status_ = HighsModelStatus::kOptimal;
@@ -4296,7 +4311,7 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
         }
         highsLogUser(
             options_.log_options, HighsLogType::kInfo,
-            "Pure postsolve yields primal %ssolution, but no basis: model "
+            "\nPure postsolve yields primal %ssolution, but no basis: model "
             "status is %s\n",
             solution_.dual_valid ? "and dual " : "",
             modelStatusToString(model_status_).c_str());

highs/lp_data/HighsInterface.cpp

@@ -2665,7 +2665,7 @@ HighsStatus Highs::lpKktCheck(const std::string& message) {
   //  highsLogUser(options.log_options, HighsLogType::kInfo,
   //               "Highs::lpKktCheck: %s\n", message.c_str());
   if (this->model_status_ == HighsModelStatus::kOptimal)
-    reportLpKktFailures(model_.lp_, options, info, "LP");
+    reportLpKktFailures(model_.lp_, options, info);
   // get_residuals is false when there is a valid basis, since
   // residual errors are assumed to be small, so
   // info.num_primal_residual_errors = -1, since they aren't

highs/lp_data/HighsSolution.cpp

@@ -1665,15 +1665,35 @@ bool isBasisConsistent(const HighsLp& lp, const HighsBasis& basis) {
   return num_basic_variables == lp.num_row_;
 }
 
+bool isColPrimalSolutionRightSize(const HighsLp& lp,
+                                  const HighsSolution& solution) {
+  return solution.col_value.size() == static_cast<size_t>(lp.num_col_);
+}
+
+bool isRowPrimalSolutionRightSize(const HighsLp& lp,
+                                  const HighsSolution& solution) {
+  return solution.row_value.size() == static_cast<size_t>(lp.num_row_);
+}
+
 bool isPrimalSolutionRightSize(const HighsLp& lp,
                                const HighsSolution& solution) {
-  return solution.col_value.size() == static_cast<size_t>(lp.num_col_) &&
-         solution.row_value.size() == static_cast<size_t>(lp.num_row_);
+  return isColPrimalSolutionRightSize(lp, solution) &&
+         isRowPrimalSolutionRightSize(lp, solution);
+}
+
+bool isColDualSolutionRightSize(const HighsLp& lp,
+                                const HighsSolution& solution) {
+  return solution.col_dual.size() == static_cast<size_t>(lp.num_col_);
+}
+
+bool isRowDualSolutionRightSize(const HighsLp& lp,
+                                const HighsSolution& solution) {
+  return solution.row_dual.size() == static_cast<size_t>(lp.num_row_);
 }
 
 bool isDualSolutionRightSize(const HighsLp& lp, const HighsSolution& solution) {
-  return solution.col_dual.size() == static_cast<size_t>(lp.num_col_) &&
-         solution.row_dual.size() == static_cast<size_t>(lp.num_row_);
+  return isColDualSolutionRightSize(lp, solution) &&
+         isRowDualSolutionRightSize(lp, solution);
 }
 
 bool isSolutionRightSize(const HighsLp& lp, const HighsSolution& solution) {
@@ -1687,7 +1707,7 @@ bool isBasisRightSize(const HighsLp& lp, const HighsBasis& basis) {
 }
 
 void reportLpKktFailures(const HighsLp& lp, const HighsOptions& options,
-                         const HighsInfo& info, const std::string& solver) {
+                         const HighsInfo& info, const std::string& message) {
   const HighsLogOptions& log_options = options.log_options;
   double primal_feasibility_tolerance = options.primal_feasibility_tolerance;
   double dual_feasibility_tolerance = options.dual_feasibility_tolerance;
@@ -1714,7 +1734,8 @@ void reportLpKktFailures(const HighsLp& lp, const HighsOptions& options,
   HighsLogType log_type =
       has_kkt_failures ? HighsLogType::kWarning : HighsLogType::kInfo;
 
-  highsLogUser(log_options, log_type, "LP solution KKT conditions\n");
+  highsLogUser(log_options, log_type, "LP solution KKT conditions%s%s\n",
+               message == "" ? "" : ": ", message == "" ? "" : message.c_str());
 
   highsLogUser(
       log_options, HighsLogType::kInfo,

highs/lp_data/HighsSolution.h

@@ -148,14 +148,24 @@ void resetModelStatusAndHighsInfo(HighsModelStatus& model_status,
                                   HighsInfo& highs_info);
 bool isBasisConsistent(const HighsLp& lp, const HighsBasis& basis);
 
+bool isColPrimalSolutionRightSize(const HighsLp& lp,
+                                  const HighsSolution& solution);
+bool isRowPrimalSolutionRightSize(const HighsLp& lp,
+                                  const HighsSolution& solution);
 bool isPrimalSolutionRightSize(const HighsLp& lp,
                                const HighsSolution& solution);
+
+bool isColDualSolutionRightSize(const HighsLp& lp,
+                                const HighsSolution& solution);
+bool isRowDualSolutionRightSize(const HighsLp& lp,
+                                const HighsSolution& solution);
 bool isDualSolutionRightSize(const HighsLp& lp, const HighsSolution& solution);
+
 bool isSolutionRightSize(const HighsLp& lp, const HighsSolution& solution);
 bool isBasisRightSize(const HighsLp& lp, const HighsBasis& basis);
 
 void reportLpKktFailures(const HighsLp& lp, const HighsOptions& options,
                          const HighsInfo& highs_info,
-                         const std::string& solver = "");
+                         const std::string& message = "");
 
 #endif  // LP_DATA_HIGHSSOLUTION_H_