Merge pull request #2417 from ERGO-Code/fix-2415

Fix 2415

Author: jajhall

check/TestRays.cpp

@@ -4,7 +4,7 @@
 #include "catch.hpp"
 #include "lp_data/HConst.h"
 
-const bool dev_run = false;
+const bool dev_run = true;  // false;
 const double zero_ray_value_tolerance = 1e-14;
 
 void reportRay(std::string message, HighsInt dim, double* computed,
@@ -963,3 +963,48 @@ TEST_CASE("Rays-2415", "[highs_test_rays]") {
   REQUIRE(h.getInfo().primal_solution_status != kSolutionStatusFeasible);
   REQUIRE(h.getInfo().primal_solution_status == kSolutionStatusNone);
 }
+
+TEST_CASE("Rays-2415-primal", "[highs_test_rays]") {
+  HighsLp lp;
+  lp.num_col_ = 1;
+  lp.num_row_ = 1;
+  lp.col_cost_ = {-1};
+  lp.col_lower_ = {0};
+  lp.col_upper_ = {kHighsInf};
+  lp.row_lower_ = {0};
+  lp.row_upper_ = {kHighsInf};
+  lp.a_matrix_.start_ = {0, 1};
+  lp.a_matrix_.index_ = {0};
+  lp.a_matrix_.value_ = {1};
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  if (dev_run) printf("For unbounded LP\n");
+  for (HighsInt k = 0; k < 2; k++) {
+    bool is_mip = k == 1;
+    HighsInt require_dual_solution_status = kSolutionStatusInfeasible;
+    if (is_mip) require_dual_solution_status = kSolutionStatusNone;
+    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+    REQUIRE(h.run() == HighsStatus::kOk);
+    if (dev_run)
+      printf("Solution values are col_value[0] = %g; row_value[0] = %g\n",
+             h.getSolution().col_value[0], h.getSolution().row_value[0]);
+    if (dev_run)
+      printf("Solution duals are col_dual[0] = %g; row_dual[0] = %g\n",
+             h.getSolution().col_dual[0], h.getSolution().row_dual[0]);
+    if (dev_run)
+      printf("Dual Solution status is %d\n",
+             int(h.getInfo().dual_solution_status));
+    REQUIRE(h.getInfo().dual_solution_status != kSolutionStatusFeasible);
+    REQUIRE(h.getInfo().dual_solution_status == require_dual_solution_status);
+
+    bool has_primal_ray;
+    std::vector<double> primal_ray_value(lp.num_col_);
+    REQUIRE(h.getPrimalRay(has_primal_ray, primal_ray_value.data()) ==
+            HighsStatus::kOk);
+    REQUIRE(h.getInfo().dual_solution_status != kSolutionStatusFeasible);
+    REQUIRE(h.getInfo().dual_solution_status == require_dual_solution_status);
+    if (k == 1) break;
+    lp.integrality_ = {HighsVarType::kInteger};
+    if (dev_run) printf("\nFor unbounded MIP\n");
+  }
+}

highs/lp_data/HighsInterface.cpp

@@ -1774,7 +1774,6 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
       this->getOptionValue("allow_unbounded_or_infeasible",
                            allow_unbounded_or_infeasible);
       solve_unboundedness_problem = true;
-      lp.integrality_.clear();
       this->setOptionValue("presolve", kHighsOffString);
       this->setOptionValue("solve_relaxation", true);
       this->setOptionValue("allow_unbounded_or_infeasible", false);
@@ -1835,14 +1834,21 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
       return_status = HighsStatus::kOk;
     }
   }
+  const bool is_mip = this->model_.isMip();
   if (solve_unboundedness_problem) {
+    if (is_mip) {
+      // Unboundedness LP has been solved, but that will give dual
+      // solution status kInfeasible which, for a MIP is not correct
+      this->info_.dual_solution_status = SolutionStatus::kSolutionStatusNone;
+      this->info_.invalidateDualKkt();
+    }
     // Restore the option values
     this->setOptionValue("presolve", presolve);
     this->setOptionValue("solve_relaxation", solve_relaxation);
     this->setOptionValue("allow_unbounded_or_infeasible",
                          allow_unbounded_or_infeasible);
     if (has_primal_ray) {
-      assert(this->info_.num_dual_infeasibilities > 0);
+      assert(is_mip || this->info_.num_dual_infeasibilities > 0);
       assert(this->model_status_ == HighsModelStatus::kUnbounded);
     }
   }
@@ -2602,7 +2608,8 @@ HighsStatus Highs::lpKktCheck(const std::string& message) {
                    primal_dual_errors, get_residuals);
   //  highsLogUser(options.log_options, HighsLogType::kInfo,
   //               "Highs::lpKktCheck: %s\n", message.c_str());
-  reportLpKktFailures(model_.lp_, options, info, "LP");
+  if (this->model_status_ == HighsModelStatus::kOptimal)
+    reportLpKktFailures(model_.lp_, options, info, "LP");
   // 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