Added skeleton SlackColSubstitution as presolve/postsolve action

check/TestPresolve.cpp

@@ -616,7 +616,7 @@ TEST_CASE("presolve-slacks", "[highs_test_presolve]") {
   lp.a_matrix_.index_ = {0, 0};
   lp.a_matrix_.value_ = {1, 1};
   Highs h;
-  h.setOptionValue("output_flag", dev_run);
+  //  h.setOptionValue("output_flag", dev_run);
   REQUIRE(h.passModel(lp) == HighsStatus::kOk);
   REQUIRE(h.presolve() == HighsStatus::kOk);
   REQUIRE(h.getPresolvedLp().num_col_ == 0);
@@ -636,6 +636,6 @@ TEST_CASE("presolve-slacks", "[highs_test_presolve]") {
   REQUIRE(h.passModel(lp) == HighsStatus::kOk);
   REQUIRE(h.run() == HighsStatus::kOk);
   REQUIRE(h.presolve() == HighsStatus::kOk);
-  //  REQUIRE(h.getPresolvedLp().num_col_ == 0);
-  //  REQUIRE(h.getPresolvedLp().num_row_ == 0);
+  REQUIRE(h.getPresolvedLp().num_col_ == 2);
+  REQUIRE(h.getPresolvedLp().num_row_ == 2);
 }

src/lp_data/Highs.cpp

@@ -1251,6 +1251,40 @@ HighsStatus Highs::run() {
     // Run solver.
     bool have_optimal_solution = false;
     // ToDo Put solution of presolved problem in a separate method
+
+    if (!this->options_.presolve_remove_slacks) {
+      HighsLp& reduced_lp = presolve_.getReducedProblem();
+      HighsInt num_double_slack = 0;
+      HighsInt num_slack = 0;
+      HighsInt num_zero_cost_slack = 0;
+      HighsInt num_unit_coeff_slack = 0;
+      double min_slack_coeff = kHighsInf;
+      double max_slack_coeff = -kHighsInf;
+      std::vector<bool> found_slack;
+      found_slack.assign(reduced_lp.num_row_, false);
+      for (HighsInt iCol = 0; iCol < reduced_lp.num_col_; iCol++) {
+	HighsInt nnz = reduced_lp.a_matrix_.start_[iCol+1] - reduced_lp.a_matrix_.start_[iCol];
+	if (nnz != 1) continue;
+	HighsInt iRow = reduced_lp.a_matrix_.index_[reduced_lp.a_matrix_.start_[iCol]];
+	if (found_slack[iRow]) {
+	  num_double_slack++;
+	  continue;
+	}
+	if (reduced_lp.row_lower_[iRow] != reduced_lp.row_upper_[iRow]) continue;
+	double coeff = std::fabs(reduced_lp.a_matrix_.value_[reduced_lp.a_matrix_.start_[iCol]]);
+	if (coeff == 1.0) num_unit_coeff_slack++;
+	min_slack_coeff = std::min(coeff,min_slack_coeff);
+	max_slack_coeff = std::max(coeff,max_slack_coeff);
+	found_slack[iRow] = true;
+	num_slack++;
+	if (reduced_lp.col_cost_[iCol] == 0) num_zero_cost_slack++;
+      }
+      printf("grepSlack,model,col,slack,unit coeff,zero_cost,double,min coeff, max_coeff\n");
+      printf("grepSlack,%s,%d, %d, %d, %d, %d, %g, %g\n", this->model_.lp_.model_name_.c_str(),
+	     int(reduced_lp.num_col_), int(num_slack), int(num_unit_coeff_slack), int(num_zero_cost_slack), int(num_double_slack),
+	     min_slack_coeff, max_slack_coeff);
+    }
+
     switch (model_presolve_status_) {
       case HighsPresolveStatus::kNotPresolved: {
         ekk_instance_.lp_name_ = "Original LP";
@@ -1275,7 +1309,7 @@ HighsStatus Highs::run() {
         break;
       }
       case HighsPresolveStatus::kReduced: {
-        HighsLp& reduced_lp = presolve_.getReducedProblem();
+	HighsLp& reduced_lp = presolve_.getReducedProblem();
         reduced_lp.setMatrixDimensions();
         if (kAllowDeveloperAssert) {
           // Validate the reduced LP

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, true);//false);
+    records.push_back(record_bool);
+
     record_int = new OptionRecordInt(
         "presolve_substitution_maxfillin",
         "Maximal fillin allowed for substitutions in presolve", advanced,

src/presolve/HighsPostsolveStack.cpp

@@ -1351,4 +1351,52 @@ void HighsPostsolveStack::DuplicateColumn::transformToPresolvedSpace(
   primalSol[col] = primalSol[col] + colScale * primalSol[duplicateCol];
 }
 
+void HighsPostsolveStack::SlackColSubstitution::undo(
+    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
+    const std::vector<Nonzero>& colValues, HighsSolution& solution,
+    HighsBasis& basis) {
+  // a (removed) cut may have been used in this reduction.
+  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;
+
+  // compute the row dual value such that reduced cost of basic column is 0
+  if (isModelRow) {
+    solution.row_dual[row] = 0;
+    HighsCDouble dualval = colCost;
+    for (const auto& colVal : colValues) {
+      if (static_cast<size_t>(colVal.index) < solution.row_dual.size())
+        dualval -= colVal.value * solution.row_dual[colVal.index];
+    }
+    solution.row_dual[row] = double(dualval / colCoef);
+  }
+
+  solution.col_dual[col] = 0;
+
+  // set basis status if necessary
+  if (!basis.valid) return;
+
+  basis.col_status[col] = HighsBasisStatus::kBasic;
+  if (isModelRow)
+    basis.row_status[row] = computeRowStatus(solution.row_dual[row], rowType);
+}
+
 }  // namespace presolve

src/presolve/HighsPostsolveStack.h

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