diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index 86c5dfc469..e4263c010e 100644 --- a/check/TestPresolve.cpp +++ b/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); +} diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index 823bff0047..c1c32546a6 100644 --- a/src/lp_data/Highs.cpp +++ b/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, diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index b6513c683a..cfb395a3fc 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, 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/HPresolve.cpp b/src/presolve/HPresolve.cpp index d588210a76..5360c9e092 100644 --- a/src/presolve/HPresolve.cpp +++ b/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) diff --git a/src/presolve/HPresolve.h b/src/presolve/HPresolve.h index fc859aa2a3..2be1b24218 100644 --- a/src/presolve/HPresolve.h +++ b/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); diff --git a/src/presolve/HighsPostsolveStack.cpp b/src/presolve/HighsPostsolveStack.cpp index bec3e4945f..2e78d9a1a7 100644 --- a/src/presolve/HighsPostsolveStack.cpp +++ b/src/presolve/HighsPostsolveStack.cpp @@ -13,6 +13,7 @@ #include #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& 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(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 diff --git a/src/presolve/HighsPostsolveStack.h b/src/presolve/HighsPostsolveStack.h index 2200977b5b..f3409d074c 100644 --- a/src/presolve/HighsPostsolveStack.h +++ b/src/presolve/HighsPostsolveStack.h @@ -219,6 +219,16 @@ class HighsPostsolveStack { void transformToPresolvedSpace(std::vector& primalSol) const; }; + struct SlackColSubstitution { + double rhs; + HighsInt row; + HighsInt col; + + void undo(const HighsOptions& options, + const std::vector& 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 + void slackColSubstitution(HighsInt row, HighsInt col, double rhs, + const HighsMatrixSlice& 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 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 diff --git a/src/qpsolver/dantzigpricing.hpp b/src/qpsolver/dantzigpricing.hpp index 05b0d9c00c..3c305529b2 100644 --- a/src/qpsolver/dantzigpricing.hpp +++ b/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;