Merge pull request #2126 from ERGO-Code/fix-2125

Fix 2125

FEATURES.md

@@ -17,3 +17,5 @@ Fixed bug in presolve when pointers stored in HighsMatrixSlice get invalidated w
 Primal and dual residual tolerances - applied following IPM or PDLP solution - now documented as options
 
 Highs::getCols (Highs::getRows) now runs in linear time if the internal constraint matrix is stored column-wise (row-wise). Added ensureColwise/Rowwise to the Highs class, the C API and highspy so that users can set the internal constraint matrix storage orientation
+
+When columns and rows are deleted from the incumbent LP after a basic solution has been found, HiGHS no longer invalidates the basis. Now it maintains the basic and nonbasic status of the remaining variables and constraints. When the model is re-solved, this information is used to construct a starting basis.

check/TestLpModification.cpp

@@ -2274,3 +2274,204 @@ TEST_CASE("row-wise-get-row-avgas", "[highs_data]") {
   testAvgasGetRow(h);
   testAvgasGetCol(h);
 }
+
+TEST_CASE("hot-start-after-delete", "[highs_data]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  const HighsLp& lp = h.getLp();
+  const HighsInfo& info = h.getInfo();
+  const HighsBasis& basis = h.getBasis();
+  const HighsSolution& solution = h.getSolution();
+  std::string model = "avgas";
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  h.readModel(model_file);
+  h.run();
+  HighsInt ieration_count0 = info.simplex_iteration_count;
+  if (dev_run)
+    printf("Initial solve takes %d iterations and yields objective = %g\n",
+           int(info.simplex_iteration_count), info.objective_function_value);
+
+  HighsInt max_dim = std::max(lp.num_col_, lp.num_row_);
+  std::vector<double> cost(1);
+  std::vector<double> lower(1);
+  std::vector<double> upper(1);
+  HighsInt nnz;
+  std::vector<HighsInt> start(1);
+  std::vector<HighsInt> index(max_dim);
+  std::vector<double> value(max_dim);
+  HighsInt get_num;
+  HighsInt use_col, use_row;
+  for (HighsInt k = 0; k < 2; k++) {
+    if (dev_run) {
+      for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+        printf("Col %2d is %s\n", int(iCol),
+               basis.col_status[iCol] == HighsBasisStatus::kBasic ? "basic"
+                                                                  : "nonbasic");
+      printf("\n");
+    }
+    if (k == 0) {
+      use_col = 1;  // Nonbasic
+    } else {
+      use_col = 4;  // Basic
+    }
+    if (dev_run)
+      printf(
+          "\nDeleting and adding column %1d with status \"%s\" and value %g\n",
+          int(use_col),
+          h.basisStatusToString(basis.col_status[use_col]).c_str(),
+          solution.col_value[use_col]);
+
+    h.getCols(use_col, use_col, get_num, cost.data(), lower.data(),
+              upper.data(), nnz, start.data(), index.data(), value.data());
+
+    h.deleteCols(use_col, use_col);
+    if (dev_run) basis.printScalars();
+
+    h.addCol(cost[0], lower[0], upper[0], nnz, index.data(), value.data());
+
+    h.run();
+    if (dev_run)
+      printf(
+          "After deleting and adding column %1d, solve takes %d iterations and "
+          "yields objective = %g\n",
+          int(use_col), int(info.simplex_iteration_count),
+          info.objective_function_value);
+    REQUIRE(info.simplex_iteration_count < ieration_count0);
+  }
+
+  for (HighsInt k = 0; k < 2; k++) {
+    if (dev_run) {
+      for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++)
+        printf("Row %2d is %s\n", int(iRow),
+               basis.row_status[iRow] == HighsBasisStatus::kBasic ? "basic"
+                                                                  : "nonbasic");
+    }
+    if (k == 0) {
+      use_row = 1;  // Nonbasic
+    } else {
+      use_row = 8;  // Basic
+    }
+    if (dev_run)
+      printf("\nDeleting and adding row %1d with status \"%s\" and value %g\n",
+             int(use_row),
+             h.basisStatusToString(basis.row_status[use_row]).c_str(),
+             solution.row_value[use_row]);
+
+    h.getRows(use_row, use_row, get_num, lower.data(), upper.data(), nnz,
+              start.data(), index.data(), value.data());
+
+    h.deleteRows(use_row, use_row);
+    if (dev_run) basis.printScalars();
+
+    h.addRow(lower[0], upper[0], nnz, index.data(), value.data());
+
+    h.run();
+    if (dev_run)
+      printf(
+          "After deleting and adding row %1d, solve takes %d iterations and "
+          "yields objective = %g\n",
+          int(use_row), int(info.simplex_iteration_count),
+          info.objective_function_value);
+    REQUIRE(info.simplex_iteration_count < ieration_count0);
+  }
+  std::vector<HighsInt> set = {1, 3, 4};
+  HighsInt num_set_en = set.size();
+  cost.resize(num_set_en);
+  lower.resize(num_set_en);
+  upper.resize(num_set_en);
+  start.resize(num_set_en);
+  index.resize(num_set_en * max_dim);
+  value.resize(num_set_en * max_dim);
+
+  h.getCols(num_set_en, set.data(), get_num, cost.data(), lower.data(),
+            upper.data(), nnz, start.data(), index.data(), value.data());
+
+  h.deleteCols(num_set_en, set.data());
+  if (dev_run) basis.printScalars();
+
+  h.addCols(get_num, cost.data(), lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d columns in set, solve takes %d "
+        "iterations and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+
+  h.getRows(num_set_en, set.data(), get_num, lower.data(), upper.data(), nnz,
+            start.data(), index.data(), value.data());
+
+  h.deleteRows(num_set_en, set.data());
+  if (dev_run) basis.printScalars();
+
+  h.addRows(get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d rows in set, solve takes %d iterations "
+        "and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+  std::vector<HighsInt> mask;
+  mask.assign(max_dim, 0);
+  mask[1] = 1;
+  mask[4] = 1;
+  mask[5] = 1;
+
+  h.getCols(mask.data(), get_num, cost.data(), lower.data(), upper.data(), nnz,
+            start.data(), index.data(), value.data());
+
+  h.deleteCols(mask.data());
+  if (dev_run) basis.printScalars();
+
+  h.addCols(get_num, cost.data(), lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d columns in mask, solve takes %d "
+        "iterations and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+
+  mask.assign(max_dim, 0);
+  mask[1] = 1;
+  mask[4] = 1;
+  mask[5] = 1;
+  mask[8] = 1;
+  mask[9] = 1;
+  HighsInt num_mask_en = mask.size();
+  cost.resize(num_mask_en);
+  lower.resize(num_mask_en);
+  upper.resize(num_mask_en);
+  start.resize(num_mask_en);
+  index.resize(num_mask_en * max_dim);
+  value.resize(num_mask_en * max_dim);
+
+  h.getRows(mask.data(), get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.deleteRows(mask.data());
+  if (dev_run) basis.printScalars();
+
+  h.addRows(get_num, lower.data(), upper.data(), nnz, start.data(),
+            index.data(), value.data());
+
+  h.run();
+  if (dev_run)
+    printf(
+        "After deleting and adding %d rows in mask, solve takes %d iterations "
+        "and yields objective = %g\n",
+        int(get_num), int(info.simplex_iteration_count),
+        info.objective_function_value);
+  //  REQUIRE(info.simplex_iteration_count < ieration_count0);
+}

src/lp_data/HStruct.h

@@ -58,14 +58,27 @@ struct HotStart {
 };
 
 struct HighsBasis {
+  // Logical flags for a HiGHS basis:
+  //
+  // valid: has been factored by HiGHS
+  //
+  // alien: a basis that's been set externally, so cannot be assumed
+  // to even have the right number of basic and nonbasic variables
+  //
+  // useful: a basis that may be useful
+  //
+  // Need useful since, by default, a basis is alien but not useful
   bool valid = false;
   bool alien = true;
+  bool useful = false;
   bool was_alien = true;
   HighsInt debug_id = -1;
   HighsInt debug_update_count = -1;
   std::string debug_origin_name = "None";
   std::vector<HighsBasisStatus> col_status;
   std::vector<HighsBasisStatus> row_status;
+  void print(std::string message = "") const;
+  void printScalars(std::string message = "") const;
   void invalidate();
   void clear();
 };

src/lp_data/Highs.cpp

@@ -716,6 +716,7 @@ HighsStatus Highs::readBasis(const std::string& filename) {
   // Update the HiGHS basis and invalidate any simplex basis for the model
   basis_ = read_basis;
   basis_.valid = true;
+  basis_.useful = true;
   // Follow implications of a new HiGHS basis
   newHighsBasis();
   // Can't use returnFromHighs since...
@@ -1265,16 +1266,24 @@ HighsStatus Highs::solve() {
   // is available, simplex should surely be chosen.
   const bool solver_will_use_basis = options_.solver == kSimplexString ||
                                      options_.solver == kHighsChooseString;
-  if ((basis_.valid || options_.presolve == kHighsOffString ||
-       unconstrained_lp) &&
+  const bool has_basis = basis_.useful;
+  if (has_basis) {
+    assert(basis_.col_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_col_));
+    assert(basis_.row_status.size() ==
+           static_cast<size_t>(incumbent_lp.num_row_));
+  }
+  if (basis_.valid) assert(basis_.useful);
+
+  if ((has_basis || options_.presolve == kHighsOffString || unconstrained_lp) &&
       solver_will_use_basis) {
     // There is a valid basis for the problem, presolve is off, or LP
     // has no constraint matrix, and the solver will use the basis
     ekk_instance_.lp_name_ =
         "LP without presolve, or with basis, or unconstrained";
     // If there is a valid HiGHS basis, refine any status values that
     // are simply HighsBasisStatus::kNonbasic
-    if (basis_.valid) refineBasis(incumbent_lp, solution_, basis_);
+    if (basis_.useful) refineBasis(incumbent_lp, solution_, basis_);
     solveLp(incumbent_lp,
             "Solving LP without presolve, or with basis, or unconstrained",
             this_solve_original_lp_time);
@@ -1416,6 +1425,7 @@ HighsStatus Highs::solve() {
         basis_.debug_origin_name = "Presolve to empty";
         basis_.valid = true;
         basis_.alien = false;
+        basis_.useful = true;
         basis_.was_alien = false;
         solution_.value_valid = true;
         solution_.dual_valid = true;
@@ -1551,6 +1561,7 @@ HighsStatus Highs::solve() {
           solution_.dual_valid = true;
           // Set basis and its status
           basis_.valid = true;
+          basis_.useful = true;
           basis_.col_status = presolve_.data_.recovered_basis_.col_status;
           basis_.row_status = presolve_.data_.recovered_basis_.row_status;
           basis_.debug_origin_name += ": after postsolve";
@@ -2339,6 +2350,7 @@ HighsStatus Highs::setBasis(const HighsBasis& basis,
     basis_ = basis;
   }
   basis_.valid = true;
+  basis_.useful = true;
   if (origin != "") basis_.debug_origin_name = origin;
   assert(basis_.debug_origin_name != "");
   assert(!basis_.alien);
@@ -4157,6 +4169,7 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
       // Set basis and its status
       //
       // basis_.valid = true;
+      // basis_.useful = true;
       // basis_.col_status = presolve_.data_.recovered_basis_.col_status;
       // basis_.row_status = presolve_.data_.recovered_basis_.row_status;
       basis_.debug_origin_name += ": after postsolve";
@@ -4270,10 +4283,12 @@ void Highs::forceHighsSolutionBasisSize() {
   if (basis_.col_status.size() != static_cast<size_t>(model_.lp_.num_col_)) {
     basis_.col_status.resize(model_.lp_.num_col_);
     basis_.valid = false;
+    basis_.useful = false;
   }
   if (basis_.row_status.size() != static_cast<size_t>(model_.lp_.num_row_)) {
     basis_.row_status.resize(model_.lp_.num_row_);
     basis_.valid = false;
+    basis_.useful = false;
   }
 }