Now ensuring that IPM isn't used to try to get a basic solution when computing rays

check/TestRays.cpp

@@ -716,7 +716,6 @@ TEST_CASE("Rays-464a", "[highs_test_rays]") {
   REQUIRE(highs.setOptionValue("allow_unbounded_or_infeasible",
                                allow_unbounded_or_infeasible) ==
           HighsStatus::kOk);
-
   REQUIRE(highs.setOptionValue("presolve", kHighsOffString) ==
           HighsStatus::kOk);
   std::string presolve_status = "off";
@@ -962,6 +961,8 @@ TEST_CASE("Rays-2415", "[highs_test_rays]") {
           HighsStatus::kOk);
   REQUIRE(h.getInfo().primal_solution_status != kSolutionStatusFeasible);
   REQUIRE(h.getInfo().primal_solution_status == kSolutionStatusNone);
+
+  h.resetGlobalScheduler(true);
 }
 
 TEST_CASE("Rays-2415-primal", "[highs_test_rays]") {
@@ -1007,4 +1008,70 @@ TEST_CASE("Rays-2415-primal", "[highs_test_rays]") {
     lp.integrality_ = {HighsVarType::kInteger};
     if (dev_run) printf("\nFor unbounded MIP\n");
   }
+
+  h.resetGlobalScheduler(true);
+}
+
+TEST_CASE("issue-3413", "[highs_test_rays]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  // An infeasible LP
+  std::vector<HighsInt> index = {0};
+  std::vector<double> value = {1};
+  h.addCol(0.0, 0.0, kHighsInf, 0, nullptr, nullptr);
+  h.addRow(-kHighsInf, -1, 1, index.data(), value.data());
+  h.setOptionValue(kSolverString, kIpmString);
+  h.run();
+  bool has_dual_ray = false;
+  std::vector<double> dual_ray_value = {NAN};
+  h.getDualRay(has_dual_ray, dual_ray_value.data());
+  REQUIRE(dual_ray_value[0] == -1);
+  h.clear();
+
+  // An unbounded LP
+  h.setOptionValue("output_flag", dev_run);
+  h.addCol(-1.0, 0.0, kHighsInf, 0, nullptr, nullptr);
+  h.addRow(1, kHighsInf, 1, index.data(), value.data());
+  h.setOptionValue(kSolverString, kIpmString);
+  h.run();
+  bool has_primal_ray = false;
+  std::vector<double> primal_ray_value = {NAN};
+  h.getPrimalRay(has_primal_ray, primal_ray_value.data());
+  REQUIRE(primal_ray_value[0] == 1);
+  h.clear();
+
+  // An infeasible problem with a semi-integer variable
+  h.setOptionValue("output_flag", dev_run);
+  index = {0, 1};
+  value = {1, -1};
+  h.addVar(2.0, 3.0);
+  h.addVar(4.0, 4.0);
+  h.changeColIntegrality(0, HighsVarType::kSemiInteger);
+  h.addRow(0.0, kHighsInf, 2, index.data(), value.data());
+  h.run();
+
+  has_dual_ray = false;
+  dual_ray_value = {NAN, NAN};
+  h.getDualRay(has_dual_ray, dual_ray_value.data());
+  REQUIRE(dual_ray_value[0] == 1);
+  h.clear();
+
+  // An unbounded problem with a semi-integer variable
+  h.setOptionValue("output_flag", dev_run);
+  index = {1};
+  value = {1};
+  h.addCol(-1.0, 2.0, 3.0, 0, nullptr, nullptr);
+  h.addCol(-1.0, 0.0, kHighsInf, 0, nullptr, nullptr);
+  h.changeColIntegrality(0, HighsVarType::kSemiInteger);
+  h.addRow(1, kHighsInf, 1, index.data(), value.data());
+  h.run();
+
+  has_primal_ray = false;
+  primal_ray_value = {NAN, NAN};
+  h.writeModel("");
+  h.writeSolution("", 1);
+  h.getPrimalRay(has_primal_ray, primal_ray_value.data());
+  REQUIRE(primal_ray_value[0] == 0);
+  REQUIRE(primal_ray_value[1] == 1);
+  h.resetGlobalScheduler(true);
 }

highs/lp_data/HighsInterface.cpp

@@ -1590,6 +1590,8 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
   HighsHessian hessian;
   bool solve_relaxation;
   std::string presolve;
+  std::string solver;
+  bool allow_unbounded_or_infeasible;
   bool solve_feasibility_problem = false;
   const bool is_qp = model_.isQp();
 
@@ -1607,14 +1609,18 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
       }
       highsLogUser(options_.log_options, HighsLogType::kInfo,
                    "Solving LP to try to compute dual ray\n");
-      // Save the column costs, integrality, any Hessian and the
-      // presolve setting
+      // Save the column costs, any Hessian, the relaxation, presolve
+      // and solver settings, so that the feasibility problem is
+      // solved to yield a basic solution
       col_cost = lp.col_cost_;
       if (is_qp) hessian = model_.hessian_;
-      this->getOptionValue("presolve", presolve);
-      this->getOptionValue("solve_relaxation", solve_relaxation);
+      solve_relaxation = this->options_.solve_relaxation;
+      presolve = this->options_.presolve;
+      solver = this->options_.solver;
+      allow_unbounded_or_infeasible =
+          this->options_.allow_unbounded_or_infeasible;
       solve_feasibility_problem = true;
-      // Zero the costs, integrality and Hessian
+      // Zero the costs and Hessian
       std::vector<double> zero_costs;
       zero_costs.assign(lp.num_col_, 0);
       // Take a copy of the primal ray record, since this will be
@@ -1630,9 +1636,22 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
         HighsHessian zero_hessian;
         this->passHessian(zero_hessian);
       }
-      this->setOptionValue("presolve", kHighsOffString);
-      this->setOptionValue("solve_relaxation", true);
+      //
+      this->options_.solve_relaxation = true;
+      this->options_.presolve = kHighsOffString;
+      this->options_.solver = kSimplexString;
+      this->options_.allow_unbounded_or_infeasible = false;
       HighsStatus call_status = this->run();
+      // Restore the column costs, any Hessian, the relaxation,
+      // presolve and solver settings
+      this->options_.solve_relaxation = solve_relaxation;
+      this->options_.presolve = presolve;
+      this->options_.solver = solver;
+      this->options_.allow_unbounded_or_infeasible =
+          allow_unbounded_or_infeasible;
+
+      lp.col_cost_ = col_cost;
+      if (is_qp) model_.hessian_ = hessian;
       if (call_status != HighsStatus::kOk) return_status = call_status;
       has_dual_ray = ekk_instance_.dual_ray_record_.index != kNoRayIndex;
       has_invert = ekk_instance_.status_.has_invert;
@@ -1677,10 +1696,7 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
     // information has been lost. Reverting the objective function via
     // Highs calls clears info_, so better to just copy the data
     // directly and set the info_ entries that are no longer valid
-    lp.col_cost_ = col_cost;
-    if (is_qp) model_.hessian_ = hessian;
-    this->setOptionValue("presolve", presolve);
-    this->setOptionValue("solve_relaxation", solve_relaxation);
+    //
     // The relaxation for an infeasible MIP may be feasible - so no
     // ray is generated - so make sure (#2415) that the primal
     // solution status is reset
@@ -1690,7 +1706,8 @@ HighsStatus Highs::getDualRayInterface(bool& has_dual_ray,
     this->info_.objective_function_value = 0;
     this->info_.invalidateDualKkt();
     if (has_dual_ray) {
-      assert(this->info_.num_primal_infeasibilities > 0);
+      assert(this->info_.num_primal_infeasibilities > 0 ||
+             this->info_.num_semi_infeasibilities > 0);
       assert(this->model_status_ == HighsModelStatus::kInfeasible);
     } else {
       // If someone has tried to get a dual ray for a feasible problem
@@ -1720,8 +1737,9 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
   assert(!lp.is_moved_);
   has_primal_ray = ekk_instance_.primal_ray_record_.index != kNoRayIndex;
 
-  std::string presolve;
   bool solve_relaxation;
+  std::string presolve;
+  std::string solver;
   bool allow_unbounded_or_infeasible;
   bool solve_unboundedness_problem = false;
 
@@ -1740,15 +1758,24 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
       }
       highsLogUser(options_.log_options, HighsLogType::kInfo,
                    "Solving LP to try to compute primal ray\n");
-      this->getOptionValue("presolve", presolve);
-      this->getOptionValue("solve_relaxation", solve_relaxation);
-      this->getOptionValue("allow_unbounded_or_infeasible",
-                           allow_unbounded_or_infeasible);
+      // Save the relaxation, presolve and solver settings
+      solve_relaxation = this->options_.solve_relaxation;
+      presolve = this->options_.presolve;
+      solver = this->options_.solver;
+      allow_unbounded_or_infeasible =
+          this->options_.allow_unbounded_or_infeasible;
       solve_unboundedness_problem = true;
-      this->setOptionValue("presolve", kHighsOffString);
-      this->setOptionValue("solve_relaxation", true);
-      this->setOptionValue("allow_unbounded_or_infeasible", false);
+      this->options_.solve_relaxation = true;
+      this->options_.presolve = kHighsOffString;
+      this->options_.solver = kSimplexString;
+      this->options_.allow_unbounded_or_infeasible = false;
       HighsStatus call_status = this->run();
+      // Restore the relaxation, presolve and solver settings
+      this->options_.solve_relaxation = solve_relaxation;
+      this->options_.presolve = presolve;
+      this->options_.solver = solver;
+      this->options_.allow_unbounded_or_infeasible =
+          allow_unbounded_or_infeasible;
       if (call_status != HighsStatus::kOk) return_status = call_status;
       has_primal_ray = ekk_instance_.primal_ray_record_.index != kNoRayIndex;
       has_invert = ekk_instance_.status_.has_invert;
@@ -1813,11 +1840,6 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
       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(is_mip || this->info_.num_dual_infeasibilities > 0);
       assert(this->model_status_ == HighsModelStatus::kUnbounded);