diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index b0a4dbc95b..7b676dc011 100644 --- a/check/TestPresolve.cpp +++ b/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); } diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index c1f1f48b72..87d6309d6f 100644 --- a/src/lp_data/Highs.cpp +++ b/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 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 diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 6e2be6cc78..8103f4db7b 100644 --- a/src/lp_data/HighsOptions.h +++ b/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, diff --git a/src/presolve/HighsPostsolveStack.cpp b/src/presolve/HighsPostsolveStack.cpp index bec3e4945f..79d582a7f4 100644 --- a/src/presolve/HighsPostsolveStack.cpp +++ b/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& rowValues, + const std::vector& colValues, HighsSolution& solution, + HighsBasis& basis) { + // a (removed) cut may have been used in this reduction. + bool isModelRow = static_cast(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(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 diff --git a/src/presolve/HighsPostsolveStack.h b/src/presolve/HighsPostsolveStack.h index 2200977b5b..79a658f977 100644 --- a/src/presolve/HighsPostsolveStack.h +++ b/src/presolve/HighsPostsolveStack.h @@ -219,6 +219,19 @@ class HighsPostsolveStack { void transformToPresolvedSpace(std::vector& primalSol) const; }; + struct SlackColSubstitution { + double rhs; + double colCost; + HighsInt row; + HighsInt col; + RowType rowType; + + void undo(const HighsOptions& options, + const std::vector& rowValues, + const std::vector& 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 + void slackColSubstitution(HighsInt row, HighsInt col, double rhs, + double colCost, RowType rowType, + const HighsMatrixSlice& rowVec, + const HighsMatrixSlice& 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 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