Merge pull request #3135 from ERGO-Code/fix-3118-2

Fix 3118

Author: jajhall

check/TestCheckSolution.cpp

@@ -78,7 +78,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
       std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   Highs highs;
   const HighsInfo& info = highs.getInfo();
-  if (dev_run) printf("\n********************\nSolving from scratch\n");
+  if (dev_run)
+    printf(
+        "\n********************\nSolving from scratch\n********************\n");
   highs.setOptionValue("output_flag", dev_run);
   highs.readModel(model_file);
   HighsLp lp = highs.getLp();
@@ -92,7 +94,7 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
 
   std::string solution_file = test_name + model + ".sol";
   if (dev_run) REQUIRE(highs.writeSolution("") == HighsStatus::kOk);
-  ;
+
   REQUIRE(highs.writeSolution(solution_file) == HighsStatus::kOk);
 
   highs.clear();
@@ -102,7 +104,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   bool valid, integral, feasible;
   if (test0) {
     if (dev_run)
-      printf("\n***************************\nSolving from saved solution\n");
+      printf(
+          "\n***************************\nSolving from saved "
+          "solution\n***************************\n");
     highs.setOptionValue("output_flag", dev_run);
     highs.readModel(model_file);
 
@@ -124,7 +128,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   const bool test1 = other_tests;
   if (test1) {
     if (dev_run)
-      printf("\n***************************\nSolving from solution file\n");
+      printf(
+          "\n**************************\nSolving from solution "
+          "file\n**************************\n");
     highs.setOptionValue("output_flag", dev_run);
     highs.readModel(model_file);
 
@@ -147,8 +153,8 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   if (test2) {
     if (dev_run)
       printf(
-          "\n***************************\nSolving from saved integer "
-          "solution\n");
+          "\n***********************************\nSolving from saved integer "
+          "solution\n***********************************\n");
     highs.setOptionValue("output_flag", dev_run);
     highs.readModel(model_file);
 
@@ -180,7 +186,8 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   if (test3) {
     if (dev_run)
       printf(
-          "\n***************************\nSolving from column solution file\n");
+          "\n*********************************\nSolving from column solution "
+          "file\n*********************************\n");
     std::string column_solution_file =
         std::string(HIGHS_DIR) + "/check/instances/flugpl_integer.sol";
 
@@ -206,8 +213,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   if (test4) {
     if (dev_run)
       printf(
-          "\n***************************\nSolving from illegal column solution "
-          "file\n");
+          "\n*****************************************\nSolving from illegal "
+          "column solution "
+          "file\n*****************************************\n");
     std::string column_solution_file =
         std::string(HIGHS_DIR) + "/check/instances/flugpl_illegal_integer.sol";
 
@@ -225,8 +233,9 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
     HighsSolution starting_solution = optimal_solution;
     if (dev_run)
       printf(
-          "\n***************************\nSolving from partial integer "
-          "solution\n");
+          "\n*************************************\nSolving from partial "
+          "integer "
+          "solution\n*************************************\n");
     highs.setOptionValue("output_flag", dev_run);
     highs.readModel(model_file);
 
@@ -256,8 +265,8 @@ TEST_CASE("check-set-mip-solution", "[highs_check_solution]") {
   if (test6) {
     if (dev_run)
       printf(
-          "\n***************************\nSolving from sparse integer "
-          "solution\n");
+          "\n************************************\nSolving from sparse integer "
+          "solution\n************************************\n");
     HighsInt num_integer_variable = 0;
     for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
       if (lp.integrality_[iCol] == HighsVarType::kInteger)

check/TestMipSolver.cpp

@@ -1483,3 +1483,125 @@ TEST_CASE("issue-2975", "[highs_test_mip_solver]") {
 
   highs.resetGlobalScheduler(true);
 }
+
+TEST_CASE("issue-3118", "[highs_test_mip_solver]") {
+  const double M = 1e10;
+  //   min    x +   y
+  //   s.t.   x + M*y = 1
+  //        M*x +   y = 1
+  //          x, y binary
+  //   Initial "solution" x = y = 1/M
+  //
+  // This problem is found infeasible in presolve, but the point x = y
+  // = 1/M is integer feasible to within the tolerances, so is
+  // acccepted as an optimal solution to the problem
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+
+  HighsInt x = 0;
+  HighsInt y = 1;
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_lower_ = {0., 0.};
+  lp.col_upper_ = {1., 1.};
+  lp.col_cost_ = {1., 1.};
+  lp.integrality_ = {HighsVarType::kInteger, HighsVarType::kInteger};
+  lp.row_lower_ = {1., 1.};
+  lp.row_upper_ = {1., 1.};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {x, y, x, y};
+  lp.a_matrix_.value_ = {1., M, M, 1};
+  highs.passModel(lp);
+
+  std::vector<double> solution_values(lp.num_col_, 1 / M);
+  highs.setSolution(2, nullptr, solution_values.data());
+
+  highs.run();
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(std::abs(2 / M - highs.getInfo().objective_function_value) < 1e-9);
+  if (dev_run) highs.writeSolution("", 1);
+
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("issue-3118a", "[highs_test_mip_solver]") {
+  const double M = 1e10;
+  //   max f = x
+  //   s.t.    x - M*y  = -1
+  //           x       <=  b
+  //           x, y integer in [0, 2]x[0, 1]
+  //
+  // The MIP is not feasible over the integers for any b
+  //
+  // With initial "solution" x = 0; y = 1/M
+  //
+  // For b = 0: (0, 1/M) with f = 0 is the only feasible solution of
+  // the relaxation, and it's also feasible for the MIP, so claiming
+  // it is optimal when presolve identifies infeasibility is clearly
+  // justified
+  //
+  // For b = 1: (1, 2/M) with f = 1 is the optimal solution of the
+  // relaxation, and it's also feasible for the MIP. However, claiming
+  // that the initial "solution" (0, 1/M) with f = 0 is optimal when
+  // presolve identifies infeasibility is still justified, because
+  // it's a point that is feasible for a MIP that is deemed infeasible
+  // by presolve
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+
+  HighsInt x = 0;
+  HighsInt y = 1;
+  HighsLp lp;
+  lp.sense_ = ObjSense::kMaximize;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {2, 1};
+  lp.col_cost_ = {1, 0};
+  lp.integrality_ = {HighsVarType::kInteger, HighsVarType::kInteger};
+  lp.row_lower_ = {-1, -kHighsInf};
+  lp.row_upper_ = {-1, 0};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.start_ = {0, 2, 3};
+  lp.a_matrix_.index_ = {x, y, x};
+  lp.a_matrix_.value_ = {1., -M, 1};
+
+  for (HighsInt k = 0; k < 2; k++) {
+    HighsInt b = k == 0 ? 0 : 1;
+    lp.row_upper_[1] = b;
+
+    highs.passModel(lp);
+
+    // Solve as MIP
+    if (dev_run)
+      printf("================\nCase b = %d (MIP)\n================\n", int(b));
+    highs.setOptionValue("solve_relaxation", false);
+    std::vector<double> solution_values = {0, 1 / M};
+    highs.setSolution(2, nullptr, solution_values.data());
+
+    highs.run();
+    if (dev_run) highs.writeSolution("", 1);
+    REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+
+    // Solve as LP
+    if (dev_run)
+      printf("===============\nCase b = %d (LP)\n===============\n", int(b));
+    highs.setOptionValue("solve_relaxation", true);
+    highs.clearSolver();
+    highs.run();
+    highs.writeSolution("", 1);
+    REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+    double lp_objective_value = highs.getInfo().objective_function_value;
+
+    solution_values = highs.getSolution().col_value;
+    highs.setSolution(2, nullptr, solution_values.data());
+
+    bool valid, integral, feasible;
+    REQUIRE(highs.assessPrimalSolution(valid, integral, feasible) ==
+            HighsStatus::kOk);
+  }
+
+  highs.resetGlobalScheduler(true);
+}

highs/Highs.h

@@ -266,8 +266,11 @@ class Highs {
 
   /**
    * @brief Assess the validity, integrality and feasibility of the
-   * current primal solution. Of value after calling
-   * Highs::readSolution
+   * current primal solution. Row values are computed and checked
+   * against what's in Highs::solution_.row_value and, if the
+   * differences exceed a tolerance, valid returns false.  If any of
+   * valid, integral or feasible is false, then assessPrimalSolution
+   * returns HighsStatus::kWarning.
    */
   HighsStatus assessPrimalSolution(bool& valid, bool& integral,
                                    bool& feasible) const;
@@ -1245,7 +1248,15 @@ class Highs {
   HighsStatus setSolution(const HighsSolution& solution);
 
   /**
-   * @brief Pass a sparse primal solution
+   * @brief Pass a primal solution. If index is not a null pointer,
+   * then it is assumed that value contains num_entries of packed
+   * values, with index defing the corresponding primal solution
+   * components. If index is a null pointer, then value is assumed to
+   * be a full primal solution.
+   *
+   * It allows a full primal solution to be passed (for MIPs) without
+   * requiring either a HighsSolution that contains (empty vectors of)
+   * spurious dual information, or a full list of indices
    */
   HighsStatus setSolution(const HighsInt num_entries, const HighsInt* index,
                           const double* value);

highs/lp_data/Highs.cpp

@@ -2575,15 +2575,37 @@ HighsStatus Highs::getColOrRowName(const HighsLp& lp, const bool is_col,
 
 HighsStatus Highs::setSolution(const HighsInt num_entries,
                                const HighsInt* index, const double* value) {
+  const bool packed = index != nullptr;
+  if (packed) {
+    const bool ok_num_entries = num_entries >= 0;
+    if (!ok_num_entries) {
+      highsLogUser(options_.log_options, HighsLogType::kError,
+                   "setSolution: Packed user solution has %d < 0 entries\n",
+                   int(num_entries));
+      assert(ok_num_entries);
+      return HighsStatus::kError;
+    }
+  } else {
+    const bool ok_num_entries = num_entries == model_.lp_.num_col_;
+    if (!ok_num_entries) {
+      highsLogUser(options_.log_options, HighsLogType::kError,
+                   "setSolution: Packed user solution has %d != %d entries\n",
+                   int(num_entries), int(model_.lp_.num_col_));
+      assert(ok_num_entries);
+      return HighsStatus::kError;
+    }
+  }
   HighsStatus return_status = HighsStatus::kOk;
   if (model_.lp_.num_col_ == 0) return return_status;
   // Warn about duplicates in index
   HighsInt num_duplicates = 0;
   std::vector<bool> is_set;
   is_set.assign(model_.lp_.num_col_, false);
-  for (HighsInt iX = 0; iX < num_entries; iX++) {
-    HighsInt iCol = index[iX];
+  const HighsInt to_ix = packed ? num_entries : model_.lp_.num_col_;
+  for (HighsInt iX = 0; iX < to_ix; iX++) {
+    HighsInt iCol = packed ? index[iX] : iX;
     if (iCol < 0 || iCol >= model_.lp_.num_col_) {
+      assert(packed);
       highsLogUser(options_.log_options, HighsLogType::kError,
                    "setSolution: User solution index %d has value %d out of "
                    "range [0, %d)\n",
@@ -2605,6 +2627,7 @@ HighsStatus Highs::setSolution(const HighsInt num_entries,
     is_set[iCol] = true;
   }
   if (num_duplicates > 0) {
+    assert(packed);
     highsLogUser(options_.log_options, HighsLogType::kWarning,
                  "setSolution: User set of indices has %d duplicate%s: last "
                  "value used\n",
@@ -2616,8 +2639,8 @@ HighsStatus Highs::setSolution(const HighsInt num_entries,
   // user, and insert the values determined by the user
   HighsSolution new_solution;
   new_solution.col_value.assign(model_.lp_.num_col_, kHighsUndefined);
-  for (HighsInt iX = 0; iX < num_entries; iX++) {
-    HighsInt iCol = index[iX];
+  for (HighsInt iX = 0; iX < to_ix; iX++) {
+    HighsInt iCol = packed ? index[iX] : iX;
     new_solution.col_value[iCol] = value[iX];
   }
   return interpretCallStatus(options_.log_options, setSolution(new_solution),
@@ -3990,6 +4013,7 @@ HighsStatus Highs::completeSolutionFromDiscreteAssignment() {
     HighsStatus return_status = assessLpPrimalSolution(
         "", options_, lp, solution_, valid, integral, feasible);
     assert(return_status != HighsStatus::kError);
+    if (return_status == HighsStatus::kError) return HighsStatus::kError;
     assert(valid);
     // If the current solution is integer feasible, then it can be
     // used by MIP solver to get a primal bound
@@ -4050,7 +4074,8 @@ HighsStatus Highs::completeSolutionFromDiscreteAssignment() {
       // There are no continuous variables, so no feasible solution can be
       // deduced
       highsLogUser(options_.log_options, HighsLogType::kInfo,
-                   "User-supplied values of discrete variables cannot yield "
+                   "No continuous variables, so user-supplied values of "
+                   "discrete variables cannot yield "
                    "feasible solution\n");
       call_run = false;
     } else {