Merge pull request ERGO-Code#2034 from ERGO-Code/fix-2001

Added code to remove explicit slacks

Author: jajhall

check/TestPresolve.cpp

@@ -601,3 +601,41 @@ TEST_CASE("write-presolved-model", "[highs_test_presolve]") {
   REQUIRE(highs.getInfo().simplex_iteration_count == -1);
   std::remove(presolved_model_file.c_str());
 }
+
+TEST_CASE("presolve-slacks", "[highs_test_presolve]") {
+  // This LP reduces to empty, because the equation is a doubleton
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 1;
+  lp.col_cost_ = {1, 0};
+  lp.col_lower_ = {0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf};
+  lp.row_lower_ = {1};
+  lp.row_upper_ = {1};
+  lp.a_matrix_.start_ = {0, 1, 2};
+  lp.a_matrix_.index_ = {0, 0};
+  lp.a_matrix_.value_ = {1, 1};
+  Highs h;
+  //  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.presolve() == HighsStatus::kOk);
+  REQUIRE(h.getPresolvedLp().num_col_ == 0);
+  REQUIRE(h.getPresolvedLp().num_row_ == 0);
+
+  lp.num_col_ = 4;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {-10, -25, 0, 0};
+  lp.col_lower_ = {0, 0, 0, 0};
+  lp.col_upper_ = {kHighsInf, kHighsInf, kHighsInf, kHighsInf};
+  lp.row_lower_ = {80, 120};
+  lp.row_upper_ = {80, 120};
+  lp.a_matrix_.start_ = {0, 2, 4, 5, 6};
+  lp.a_matrix_.index_ = {0, 1, 0, 1, 0, 1};
+  lp.a_matrix_.value_ = {1, 1, 2, 4, 1, 1};
+  REQUIRE(h.setOptionValue("presolve_remove_slacks", true) == HighsStatus::kOk);
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.presolve() == HighsStatus::kOk);
+  REQUIRE(h.getPresolvedLp().num_col_ == 2);
+  REQUIRE(h.getPresolvedLp().num_row_ == 2);
+}

src/lp_data/Highs.cpp

@@ -1252,6 +1252,7 @@ HighsStatus Highs::run() {
     // Run solver.
     bool have_optimal_solution = false;
     // ToDo Put solution of presolved problem in a separate method
+
     switch (model_presolve_status_) {
       case HighsPresolveStatus::kNotPresolved: {
         ekk_instance_.lp_name_ = "Original LP";
@@ -1527,6 +1528,10 @@ HighsStatus Highs::run() {
             options_ = save_options;
             if (return_status == HighsStatus::kError)
               return returnFromRun(return_status, undo_mods);
+            if (postsolve_iteration_count > 0)
+              highsLogUser(options_.log_options, HighsLogType::kInfo,
+                           "Required %d simplex iterations after postsolve\n",
+                           int(postsolve_iteration_count));
           }
         } else {
           highsLogUser(log_options, HighsLogType::kError,

src/lp_data/HighsOptions.h

@@ -372,6 +372,7 @@ struct HighsOptionsStruct {
   HighsInt presolve_substitution_maxfillin;
   HighsInt presolve_rule_off;
   bool presolve_rule_logging;
+  bool presolve_remove_slacks;
   bool simplex_initial_condition_check;
   bool no_unnecessary_rebuild_refactor;
   double simplex_initial_condition_tolerance;
@@ -507,6 +508,7 @@ struct HighsOptionsStruct {
         presolve_substitution_maxfillin(0),
         presolve_rule_off(0),
         presolve_rule_logging(false),
+        presolve_remove_slacks(false),
         simplex_initial_condition_check(false),
         no_unnecessary_rebuild_refactor(false),
         simplex_initial_condition_tolerance(0.0),
@@ -1324,6 +1326,11 @@ class HighsOptions : public HighsOptionsStruct {
         advanced, &presolve_rule_logging, false);
     records.push_back(record_bool);
 
+    record_bool = new OptionRecordBool("presolve_remove_slacks",
+                                       "Remove slacks after presolve", advanced,
+                                       &presolve_remove_slacks, false);
+    records.push_back(record_bool);
+
     record_int = new OptionRecordInt(
         "presolve_substitution_maxfillin",
         "Maximal fillin allowed for substitutions in presolve", advanced,

src/presolve/HPresolve.cpp

@@ -4307,6 +4307,10 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
       break;
     }
 
+    // Now consider removing slacks
+    if (options->presolve_remove_slacks)
+      HPRESOLVE_CHECKED_CALL(removeSlacks(postsolve_stack));
+
     report();
   } else {
     highsLogUser(options->log_options, HighsLogType::kInfo,
@@ -4322,6 +4326,56 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::removeSlacks(
+    HighsPostsolveStack& postsolve_stack) {
+  // SingletonColumns data structure appears not to be retained
+  // throughout presolve
+  for (HighsInt iCol = 0; iCol != model->num_col_; ++iCol) {
+    if (colDeleted[iCol]) continue;
+    if (colsize[iCol] != 1) continue;
+    if (model->integrality_[iCol] == HighsVarType::kInteger) continue;
+    HighsInt coliter = colhead[iCol];
+    HighsInt iRow = Arow[coliter];
+    assert(Acol[coliter] == iCol);
+    assert(!rowDeleted[iRow]);
+    if (model->row_lower_[iRow] != model->row_upper_[iRow]) continue;
+    double lower = model->col_lower_[iCol];
+    double upper = model->col_upper_[iCol];
+    double cost = model->col_cost_[iCol];
+    double rhs = model->row_lower_[iRow];
+    double coeff = Avalue[coliter];
+    assert(coeff);
+    // Slack is s = (rhs - a^Tx)/coeff
+    //
+    // Constraint bounds become:
+    //
+    // For coeff > 0 [rhs - coeff * upper, rhs - coeff * lower]
+    //
+    // For coeff < 0 [rhs - coeff * lower, rhs - coeff * upper]
+    model->row_lower_[iRow] =
+        coeff > 0 ? rhs - coeff * upper : rhs - coeff * lower;
+    model->row_upper_[iRow] =
+        coeff > 0 ? rhs - coeff * lower : rhs - coeff * upper;
+    if (cost) {
+      // Cost is (cost * rhs / coeff) + (col_cost - (cost/coeff) row_values)^Tx
+      double multiplier = cost / coeff;
+      for (const HighsSliceNonzero& nonzero : getRowVector(iRow)) {
+        HighsInt local_iCol = nonzero.index();
+        double local_value = nonzero.value();
+        model->col_cost_[local_iCol] -= multiplier * local_value;
+      }
+      model->offset_ += multiplier * rhs;
+    }
+    //
+    postsolve_stack.slackColSubstitution(iRow, iCol, rhs, getRowVector(iRow));
+
+    markColDeleted(iCol);
+
+    unlink(coliter);
+  }
+  return Result::kOk;
+}
+
 HPresolve::Result HPresolve::checkTimeLimit() {
   assert(timer);
   if (options->time_limit < kHighsInf && timer->read() >= options->time_limit)

src/presolve/HPresolve.h

@@ -268,6 +268,8 @@ class HPresolve {
 
   Result presolve(HighsPostsolveStack& postsolve_stack);
 
+  Result removeSlacks(HighsPostsolveStack& postsolve_stack);
+
   Result checkTimeLimit();
 
   Result checkLimits(HighsPostsolveStack& postsolve_stack);

src/presolve/HighsPostsolveStack.cpp

@@ -13,6 +13,7 @@
 #include <numeric>
 
 #include "lp_data/HConst.h"
+#include "lp_data/HighsModelUtils.h"  // For debugging #2001
 #include "lp_data/HighsOptions.h"
 #include "util/HighsCDouble.h"
 #include "util/HighsUtils.h"
@@ -1351,4 +1352,72 @@ void HighsPostsolveStack::DuplicateColumn::transformToPresolvedSpace(
   primalSol[col] = primalSol[col] + colScale * primalSol[duplicateCol];
 }
 
+void HighsPostsolveStack::SlackColSubstitution::undo(
+    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
+    HighsSolution& solution, HighsBasis& basis) {
+  bool debug_print = false;
+  // May have to determine row dual and basis status unless doing
+  // primal-only transformation in MIP solver, in which case row may
+  // no longer exist if it corresponds to a removed cut, so have to
+  // avoid exceeding array bounds of solution.row_value
+  bool isModelRow = static_cast<size_t>(row) < solution.row_value.size();
+
+  // compute primal values
+  double colCoef = 0;
+  HighsCDouble rowValue = 0;
+  for (const auto& rowVal : rowValues) {
+    if (rowVal.index == col)
+      colCoef = rowVal.value;
+    else
+      rowValue += rowVal.value * solution.col_value[rowVal.index];
+  }
+
+  assert(colCoef != 0);
+  // Row values aren't fully postsolved, so why do this?
+  if (isModelRow)
+    solution.row_value[row] =
+        double(rowValue + colCoef * solution.col_value[col]);
+
+  solution.col_value[col] = double((rhs - rowValue) / colCoef);
+
+  // If no dual values requested, return here
+  if (!solution.dual_valid) return;
+
+  // Row retains its dual value, and column has this dual value scaled by coeff
+  if (isModelRow) solution.col_dual[col] = -solution.row_dual[row] / colCoef;
+
+  // Set basis status if necessary
+  if (!basis.valid) return;
+
+  // If row is basic, then slack is basic, otherwise row retains its status
+  if (isModelRow) {
+    HighsBasisStatus save_row_basis_status = basis.row_status[row];
+    if (basis.row_status[row] == HighsBasisStatus::kBasic) {
+      basis.col_status[col] = HighsBasisStatus::kBasic;
+      basis.row_status[row] =
+          computeRowStatus(solution.row_dual[row], RowType::kEq);
+    } else if (basis.row_status[row] == HighsBasisStatus::kLower) {
+      basis.col_status[col] =
+          colCoef > 0 ? HighsBasisStatus::kUpper : HighsBasisStatus::kLower;
+    } else {
+      basis.col_status[col] =
+          colCoef > 0 ? HighsBasisStatus::kLower : HighsBasisStatus::kUpper;
+    }
+    if (debug_print)
+      printf(
+          "HighsPostsolveStack::SlackColSubstitution::undo OgRowStatus = %s; "
+          "RowStatus = %s; ColStatus = %s\n",
+          utilBasisStatusToString(save_row_basis_status).c_str(),
+          utilBasisStatusToString(basis.row_status[row]).c_str(),
+          utilBasisStatusToString(basis.col_status[col]).c_str());
+    if (basis.col_status[col] == HighsBasisStatus::kLower) {
+      assert(solution.col_dual[col] > -options.dual_feasibility_tolerance);
+    } else if (basis.col_status[col] == HighsBasisStatus::kUpper) {
+      assert(solution.col_dual[col] < options.dual_feasibility_tolerance);
+    }
+  } else {
+    basis.col_status[col] = HighsBasisStatus::kNonbasic;
+  }
+}
+
 }  // namespace presolve

src/presolve/HighsPostsolveStack.h

@@ -219,6 +219,16 @@ class HighsPostsolveStack {
     void transformToPresolvedSpace(std::vector<double>& primalSol) const;
   };
 
+  struct SlackColSubstitution {
+    double rhs;
+    HighsInt row;
+    HighsInt col;
+
+    void undo(const HighsOptions& options,
+              const std::vector<Nonzero>& rowValues, HighsSolution& solution,
+              HighsBasis& basis);
+  };
+
   /// tags for reduction
   enum class ReductionType : uint8_t {
     kLinearTransform,
@@ -234,6 +244,7 @@ class HighsPostsolveStack {
     kForcingColumnRemovedRow,
     kDuplicateRow,
     kDuplicateColumn,
+    kSlackColSubstitution,
   };
 
   HighsDataStack reductionValues;
@@ -323,6 +334,19 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kFreeColSubstitution);
   }
 
+  template <typename RowStorageFormat>
+  void slackColSubstitution(HighsInt row, HighsInt col, double rhs,
+                            const HighsMatrixSlice<RowStorageFormat>& rowVec) {
+    rowValues.clear();
+    for (const HighsSliceNonzero& rowVal : rowVec)
+      rowValues.emplace_back(origColIndex[rowVal.index()], rowVal.value());
+
+    reductionValues.push(
+        SlackColSubstitution{rhs, origRowIndex[row], origColIndex[col]});
+    reductionValues.push(rowValues);
+    reductionAdded(ReductionType::kSlackColSubstitution);
+  }
+
   template <typename ColStorageFormat>
   void doubletonEquation(HighsInt row, HighsInt colSubst, HighsInt col,
                          double coefSubst, double coef, double rhs,
@@ -710,6 +734,13 @@ class HighsPostsolveStack {
           reduction.undo(options, solution, basis);
           break;
         }
+        case ReductionType::kSlackColSubstitution: {
+          SlackColSubstitution reduction;
+          reductionValues.pop(rowValues);
+          reductionValues.pop(reduction);
+          reduction.undo(options, rowValues, solution, basis);
+          break;
+        }
         default:
           printf("Reduction case %d not handled\n",
                  int(reductions[i - 1].first));
@@ -887,6 +918,13 @@ class HighsPostsolveStack {
           reductionValues.pop(reduction);
           reduction.undo(options, solution, basis);
         }
+        case ReductionType::kSlackColSubstitution: {
+          SlackColSubstitution reduction;
+          reductionValues.pop(rowValues);
+          reductionValues.pop(reduction);
+          reduction.undo(options, rowValues, solution, basis);
+          break;
+        }
       }
     }
 #ifdef DEBUG_EXTRA

src/qpsolver/dantzigpricing.hpp

@@ -55,7 +55,6 @@ class DantzigPricing : public Pricing {
       // clang-format off
       : runtime(rt), basis(bas), redcosts(rc) {};
   // clang-format on
-
   HighsInt price(const QpVector& x, const QpVector& gradient) {
     HighsInt minidx = chooseconstrainttodrop(redcosts.getReducedCosts());
     return minidx;