diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index 8af6618f340..90ae416fa2b 100644 --- a/check/TestPresolve.cpp +++ b/check/TestPresolve.cpp @@ -910,3 +910,20 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") { REQUIRE(highs.presolve() == HighsStatus::kOk); REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible); } + +TEST_CASE("presolve-light", "[highs_test_presolve]") { + std::string model_file = + std::string(HIGHS_DIR) + "/check/instances/afiro.mps"; + Highs highs; + highs.setOptionValue("output_flag", dev_run); + highs.readModel(model_file); + for (HighsInt k = 0; k < 2; k++) { + REQUIRE(highs.presolve() == HighsStatus::kOk); + REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kReduced); + const HighsLp& presolved_lp = highs.getPresolvedLp(); + if (dev_run) printf("%s presolved LP has %d columns; %d rows and %d nonzeros\n", + k == 0 ? "Fully" : "Lightly", int(presolved_lp.num_col_), + int(presolved_lp.num_row_), int(presolved_lp.numNz())); + highs.setOptionValue("presolve_light", true); + } +} diff --git a/cmake/sources.cmake b/cmake/sources.cmake index a1683a8315a..f5c2468e749 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -558,6 +558,7 @@ set(highs_headers presolve/ICrashUtil.h presolve/ICrashX.h presolve/PresolveComponent.h + presolve/PresolveTimer.h qpsolver/a_asm.hpp qpsolver/a_quass.hpp qpsolver/basis.hpp diff --git a/highs/Highs.h b/highs/Highs.h index f8b046a99cc..a38593ff912 100644 --- a/highs/Highs.h +++ b/highs/Highs.h @@ -717,7 +717,7 @@ class Highs { * @brief Get the number of (constraint matrix) nonzeros in the incumbent * model */ - HighsInt getNumNz() const { return model_.lp_.a_matrix_.numNz(); } + HighsInt getNumNz() const { return model_.lp_.numNz(); } /** * @brief Get the number of Hessian matrix nonzeros in the incumbent model diff --git a/highs/interfaces/highs_c_api.cpp b/highs/interfaces/highs_c_api.cpp index bc5cd039916..b590c33ae2d 100644 --- a/highs/interfaces/highs_c_api.cpp +++ b/highs/interfaces/highs_c_api.cpp @@ -1184,7 +1184,7 @@ HighsInt Highs_getPresolvedNumRow(const void* highs) { } HighsInt Highs_getPresolvedNumNz(const void* highs) { - return ((Highs*)highs)->getPresolvedLp().a_matrix_.numNz(); + return ((Highs*)highs)->getPresolvedLp().numNz(); } // Gets pointers to all the public data members of HighsLp: avoids @@ -1231,7 +1231,7 @@ static HighsInt Highs_getHighsLpData(const HighsLp& lp, const HighsInt a_format, (desired_a_format == MatrixFormat::kRowwise && lp.a_matrix_.isRowwise())) { // Incumbent format is OK - *num_nz = lp.a_matrix_.numNz(); + *num_nz = lp.numNz(); if (a_start) memcpy(a_start, lp.a_matrix_.start_.data(), num_start_entries * sizeof(HighsInt)); diff --git a/highs/io/HighsIO.cpp b/highs/io/HighsIO.cpp index fc4d1e17398..8fb1cc33f7f 100644 --- a/highs/io/HighsIO.cpp +++ b/highs/io/HighsIO.cpp @@ -310,6 +310,11 @@ const std::string highsBoolToString(const bool b, const HighsInt field_width) { return b ? " true" : "false"; } +const std::string highsIntToPlural(const HighsInt i, const bool y) { + if (y) return i == 1 ? "y" : "ies"; + return i == 1 ? "" : "s"; +} + const std::string highsTimeToString(const double time) { return #ifndef NDEBUG diff --git a/highs/io/HighsIO.h b/highs/io/HighsIO.h index c35ad89dd06..91a21b8275c 100644 --- a/highs/io/HighsIO.h +++ b/highs/io/HighsIO.h @@ -111,6 +111,7 @@ std::string highsFormatToString(const char* format, ...); const std::string highsBoolToString(const bool b, const HighsInt field_width = 2); +const std::string highsIntToPlural(const HighsInt i, const bool y = false); const std::string highsInsertMdEscapes(const std::string& from_string); const std::string highsInsertMdId(const std::string& from_string); const std::string highsTimeToString(const double time); diff --git a/highs/lp_data/HConst.h b/highs/lp_data/HConst.h index ed5372ae54f..3776c43f8f5 100644 --- a/highs/lp_data/HConst.h +++ b/highs/lp_data/HConst.h @@ -88,12 +88,14 @@ enum HighsAnalysisLevel { kHighsAnalysisLevelNlaTime = 32, kHighsAnalysisLevelMipData = 64, kHighsAnalysisLevelMipTime = 128, + kHighsAnalysisLevelPresolveTime = 256, kHighsAnalysisLevelMin = kHighsAnalysisLevelNone, kHighsAnalysisLevelMax = kHighsAnalysisLevelModelData + kHighsAnalysisLevelSolverSummaryData + kHighsAnalysisLevelSolverRuntimeData + kHighsAnalysisLevelSolverTime + kHighsAnalysisLevelNlaData + kHighsAnalysisLevelNlaTime + - kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime + kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime + + kHighsAnalysisLevelPresolveTime }; enum class HighsVarType : uint8_t { @@ -275,7 +277,9 @@ enum PresolveRuleType : int { kPresolveRuleSparsify, kPresolveRuleProbing, kPresolveRuleEnumeration, - kPresolveRuleMax = kPresolveRuleEnumeration, + kPresolveRuleDualFixing, + kPresolveRuleColStuffing, + kPresolveRuleMax = kPresolveRuleColStuffing, kPresolveRuleLastAllowOff = kPresolveRuleMax, kPresolveRuleCount }; diff --git a/highs/lp_data/Highs.cpp b/highs/lp_data/Highs.cpp index e5ff245b658..750ecbadb14 100644 --- a/highs/lp_data/Highs.cpp +++ b/highs/lp_data/Highs.cpp @@ -1375,7 +1375,7 @@ HighsStatus Highs::calledOptimizeModel() { time += timer_.read(timer_.solve_clock); }; - const bool unconstrained_lp = incumbent_lp.a_matrix_.numNz() == 0; + const bool unconstrained_lp = incumbent_lp.numNz() == 0; assert(incumbent_lp.num_row_ || unconstrained_lp); const bool has_basis = basis_.useful; if (has_basis) { @@ -3676,8 +3676,8 @@ HighsPresolveStatus Highs::runPresolve(const bool force_lp_presolve, original_lp.num_col_ - reduced_lp.num_col_; presolve_.info_.n_rows_removed = original_lp.num_row_ - reduced_lp.num_row_; - presolve_.info_.n_nnz_removed = (HighsInt)original_lp.a_matrix_.numNz() - - (HighsInt)reduced_lp.a_matrix_.numNz(); + presolve_.info_.n_nnz_removed = + (HighsInt)original_lp.numNz() - (HighsInt)reduced_lp.numNz(); // Clear any scaling information inherited by the reduced LP reduced_lp.clearScale(); assert(lpDimensionsOk("RunPresolve: reduced_lp", reduced_lp, @@ -3687,7 +3687,7 @@ HighsPresolveStatus Highs::runPresolve(const bool force_lp_presolve, case HighsPresolveStatus::kReducedToEmpty: { presolve_.info_.n_cols_removed = original_lp.num_col_; presolve_.info_.n_rows_removed = original_lp.num_row_; - presolve_.info_.n_nnz_removed = (HighsInt)original_lp.a_matrix_.numNz(); + presolve_.info_.n_nnz_removed = (HighsInt)original_lp.numNz(); break; } default: diff --git a/highs/lp_data/HighsLp.h b/highs/lp_data/HighsLp.h index e3dc241593e..77236f406f2 100644 --- a/highs/lp_data/HighsLp.h +++ b/highs/lp_data/HighsLp.h @@ -64,6 +64,7 @@ class HighsLp { bool equalNames(const HighsLp& lp) const; bool equalScaling(const HighsLp& lp) const; bool isMip() const; + HighsInt numNz() const { return this->a_matrix_.numNz(); } bool hasSemiVariables() const; bool hasInfiniteCost(const double infinite_cost) const; bool hasMods() const; diff --git a/highs/lp_data/HighsLpUtils.cpp b/highs/lp_data/HighsLpUtils.cpp index d7cd3da831e..cfa9f34fff9 100644 --- a/highs/lp_data/HighsLpUtils.cpp +++ b/highs/lp_data/HighsLpUtils.cpp @@ -76,7 +76,7 @@ HighsStatus assessLp(HighsLp& lp, const HighsOptions& options) { // If the LP has no columns the matrix must be empty and there is // nothing left to test if (lp.num_col_ == 0) { - assert(!lp.a_matrix_.numNz()); + assert(!lp.numNz()); return HighsStatus::kOk; } // From here, any LP has lp.num_col_ > 0 and lp.a_matrix_.start_[lp.num_col_] @@ -92,7 +92,7 @@ HighsStatus assessLp(HighsLp& lp, const HighsOptions& options) { if (return_status == HighsStatus::kError) return return_status; // If entries have been removed from the matrix, resize the index // and value vectors to prevent bug in presolve - HighsInt lp_num_nz = lp.a_matrix_.numNz(); + HighsInt lp_num_nz = lp.numNz(); if ((HighsInt)lp.a_matrix_.index_.size() > lp_num_nz) lp.a_matrix_.index_.resize(lp_num_nz); if ((HighsInt)lp.a_matrix_.value_.size() > lp_num_nz) @@ -3102,7 +3102,7 @@ void reportPresolveReductions(const HighsLogOptions& log_options, const HighsLp& lp, const HighsLp& presolved_lp) { const HighsInt num_col_from = lp.num_col_; const HighsInt num_row_from = lp.num_row_; - const HighsInt num_nz_from = lp.a_matrix_.numNz(); + const HighsInt num_nz_from = lp.numNz(); HighsInt num_col_to = 0; HighsInt num_row_to = 0; HighsInt num_nz_to = 0; @@ -3124,7 +3124,7 @@ void reportPresolveReductions(const HighsLogOptions& log_options, case HighsPresolveStatus::kTimeout: { num_col_to = presolved_lp.num_col_; num_row_to = presolved_lp.num_row_; - num_nz_to = presolved_lp.a_matrix_.numNz(); + num_nz_to = presolved_lp.numNz(); message = presolve_status == HighsPresolveStatus::kTimeout ? "- Timeout" : ""; break; diff --git a/highs/lp_data/HighsModelUtils.cpp b/highs/lp_data/HighsModelUtils.cpp index a9310a4ea20..f20bb429481 100644 --- a/highs/lp_data/HighsModelUtils.cpp +++ b/highs/lp_data/HighsModelUtils.cpp @@ -598,10 +598,10 @@ void writeGlpsolSolution(FILE* file, const HighsOptions& options, assert(lp.row_names_.size() == static_cast(lp.num_row_)); // Determine number of nonzeros including the objective function // and, hence, determine whether there is an objective function - HighsInt num_nz = lp.a_matrix_.numNz(); + HighsInt num_nz = lp.numNz(); for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) if (lp.col_cost_[iCol]) num_nz++; - const bool empty_cost_row = num_nz == lp.a_matrix_.numNz(); + const bool empty_cost_row = num_nz == lp.numNz(); const bool has_objective = !empty_cost_row || model.hessian_.dim_; // Writes the solution using the GLPK raw style (defined in // api/wrsol.c) or pretty style (defined in api/prsol.c) @@ -702,7 +702,7 @@ void writeGlpsolSolution(FILE* file, const HighsOptions& options, const HighsInt glpsol_num_row = num_row + delta_num_row; // If the cost row isn't reported, then the number of nonzeros is // just the number in the constraint matrix - if (cost_row_location <= 0) num_nz = lp.a_matrix_.numNz(); + if (cost_row_location <= 0) num_nz = lp.numNz(); // Record the discrete nature of the model HighsInt num_integer = 0; HighsInt num_binary = 0; @@ -1515,6 +1515,10 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) { return "Probing"; } else if (rule_type == kPresolveRuleEnumeration) { return "Enumeration"; + } else if (rule_type == kPresolveRuleDualFixing) { + return "Dual fixing"; + } else if (rule_type == kPresolveRuleColStuffing) { + return "Col stuffing"; } assert(1 == 0); return "????"; diff --git a/highs/lp_data/HighsOptions.h b/highs/lp_data/HighsOptions.h index 5562deee868..df02d236895 100644 --- a/highs/lp_data/HighsOptions.h +++ b/highs/lp_data/HighsOptions.h @@ -421,6 +421,7 @@ struct HighsOptionsStruct { bool lp_presolve_requires_basis_postsolve; bool mps_parser_type_free; bool use_warm_start; + bool presolve_light; HighsInt keep_n_rows; HighsInt cost_scale_factor; HighsInt allowed_matrix_scale_factor; @@ -588,6 +589,7 @@ struct HighsOptionsStruct { lp_presolve_requires_basis_postsolve(false), mps_parser_type_free(false), use_warm_start(true), + presolve_light(false), keep_n_rows(0), cost_scale_factor(0), allowed_matrix_scale_factor(0), @@ -1452,6 +1454,11 @@ class HighsOptions : public HighsOptionsStruct { advanced, &use_warm_start, true); records.push_back(record_bool); + record_bool = new OptionRecordBool("presolve_light", + "Use only low-cost presolve rules", + advanced, &presolve_light, false); + records.push_back(record_bool); + record_int = new OptionRecordInt("keep_n_rows", "For multiple N-rows in MPS files: delete rows / " diff --git a/highs/lp_data/HighsSolve.cpp b/highs/lp_data/HighsSolve.cpp index 050e62f29f6..981662d7e4d 100644 --- a/highs/lp_data/HighsSolve.cpp +++ b/highs/lp_data/HighsSolve.cpp @@ -59,7 +59,7 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { } return return_status; }; - if (!solver_object.lp_.num_row_ || solver_object.lp_.a_matrix_.numNz() == 0) { + if (!solver_object.lp_.num_row_ || solver_object.lp_.numNz() == 0) { // LP is unconstrained due to having no rows or a zero constraint // matrix, so solve directly call_status = solveUnconstrainedLp(solver_object); @@ -208,11 +208,11 @@ HighsStatus solveUnconstrainedLp(const HighsOptions& options, const HighsLp& lp, resetModelStatusAndHighsInfo(model_status, highs_info); // Check that the LP really is unconstrained! - assert(lp.num_row_ == 0 || lp.a_matrix_.numNz() == 0); + assert(lp.num_row_ == 0 || lp.numNz() == 0); if (lp.num_row_ > 0) { // LP has rows, but should only be here if the constraint matrix // is zero - if (lp.a_matrix_.numNz() > 0) return HighsStatus::kError; + if (lp.numNz() > 0) return HighsStatus::kError; } highsLogUser(options.log_options, HighsLogType::kInfo, @@ -465,7 +465,7 @@ void assessExcessiveObjectiveBoundScaling(const HighsLogOptions log_options, double min_matrix_value = kHighsInf; double max_matrix_value = -kHighsInf; - const HighsInt num_matrix_nz = lp.a_matrix_.numNz(); + const HighsInt num_matrix_nz = lp.numNz(); for (HighsInt iEl = 0; iEl < num_matrix_nz; iEl++) assessFiniteNonzero(lp.a_matrix_.value_[iEl], min_matrix_value, max_matrix_value); diff --git a/highs/mip/MipTimer.h b/highs/mip/MipTimer.h index f1044de0658..3d55feec73c 100644 --- a/highs/mip/MipTimer.h +++ b/highs/mip/MipTimer.h @@ -11,7 +11,7 @@ #ifndef MIP_MIPTIMER_H_ #define MIP_MIPTIMER_H_ -// Clocks for profiling the MIP dual mip solver +// Clocks for profiling the MIP solver enum iClockMip { kMipClockTotal = 0, kMipClockPresolve, @@ -114,7 +114,7 @@ enum iClockMip { kNumMipClock //!< Number of MIP clocks }; -const double tolerance_percent_report = 0.1; +static const double kMipClockTolerancePercentReport = 0.1; class MipTimer { public: @@ -371,7 +371,7 @@ class MipTimer { kMipClockSearch, kMipClockPostsolve}; reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock, - kMipClockTotal, tolerance_percent_report); + kMipClockTotal, kMipClockTolerancePercentReport); }; void reportMipSolveLpClock(const HighsTimerClock& mip_timer_clock) { @@ -385,20 +385,20 @@ class MipTimer { kMipClockHipoSolveLp, kMipClockIpxSolveLp}; reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock, - kMipClockTotal); //, tolerance_percent_report); + kMipClockTotal); //, kMipClockTolerancePercentReport); }; void reportMipSubMipSolveClock(const HighsTimerClock& mip_timer_clock) { const std::vector mip_clock_list{kMipClockSubMipSolve}; reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock, - kMipClockTotal); //, tolerance_percent_report); + kMipClockTotal); //, kMipClockTolerancePercentReport); }; void reportMipPresolveClock(const HighsTimerClock& mip_timer_clock) { const std::vector mip_clock_list{kMipClockProbingPresolve, kMipClockEnumerationPresolve}; reportMipClockList("MipPrslv", mip_clock_list, mip_timer_clock, - kMipClockRunPresolve, tolerance_percent_report); + kMipClockRunPresolve, kMipClockTolerancePercentReport); }; void reportAltEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) { @@ -407,7 +407,7 @@ class MipTimer { kMipClockEvaluateRootNode2}; reportMipClockList( "AltEvaluateRootNode", mip_clock_list, mip_timer_clock, - kMipClockEvaluateRootNode); //, tolerance_percent_report); + kMipClockEvaluateRootNode); //, kMipClockTolerancePercentReport); }; void reportMipEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) { @@ -433,7 +433,7 @@ class MipTimer { }; reportMipClockList( "MipEvaluateRootNode", mip_clock_list, mip_timer_clock, - kMipClockEvaluateRootNode); //, tolerance_percent_report); + kMipClockEvaluateRootNode); //, kMipClockTolerancePercentReport); }; void reportMipRootSeparationClock(const HighsTimerClock& mip_timer_clock) { @@ -442,8 +442,9 @@ class MipTimer { kMipClockRootSeparationFinishAnalyticCentreComputation, kMipClockRootSeparationCentralRounding, kMipClockRootSeparationEvaluateRootLp}; - reportMipClockList("MipRootSeparation", mip_clock_list, mip_timer_clock, - kMipClockRootSeparation); //, tolerance_percent_report); + reportMipClockList( + "MipRootSeparation", mip_clock_list, mip_timer_clock, + kMipClockRootSeparation); //, kMipClockTolerancePercentReport); }; void reportMipSearchClock(const HighsTimerClock& mip_timer_clock) { @@ -455,7 +456,7 @@ class MipTimer { // kMipClock@ }; reportMipClockList("MipSerch", mip_clock_list, mip_timer_clock, - kMipClockSearch, tolerance_percent_report); + kMipClockSearch, kMipClockTolerancePercentReport); }; void reportMipDiveClock(const HighsTimerClock& mip_timer_clock) { @@ -463,7 +464,7 @@ class MipTimer { kMipClockDiveEvaluateNode, kMipClockDivePrimalHeuristics, kMipClockTheDive, kMipClockBacktrackPlunge, kMipClockPerformAging2}; reportMipClockList("MipDive_", mip_clock_list, mip_timer_clock, - kMipClockDive, tolerance_percent_report); + kMipClockDive, kMipClockTolerancePercentReport); }; void reportMipDivePrimalHeuristicsClock( @@ -472,7 +473,7 @@ class MipTimer { kMipClockDiveRandomizedRounding, kMipClockDiveRens, kMipClockDiveRins}; reportMipClockList("MipDivePrimalHeuristics", mip_clock_list, mip_timer_clock, kMipClockDivePrimalHeuristics, - tolerance_percent_report); + kMipClockTolerancePercentReport); }; void reportMipNodeSearchClock(const HighsTimerClock& mip_timer_clock) { @@ -481,8 +482,9 @@ class MipTimer { // kMipClockSearchBacktrack, kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1, kMipClockNodeSearchSeparation}; //, kMipClockStoreBasis}; - reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock, - kMipClockNodeSearch); //, tolerance_percent_report); + reportMipClockList( + "MipNodeSearch", mip_clock_list, mip_timer_clock, + kMipClockNodeSearch); //, kMipClockTolerancePercentReport); }; void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) { @@ -490,7 +492,7 @@ class MipTimer { kMipClockImplboundSepa, kMipClockCliqueSepa, kMipClockTableauSepa, kMipClockPathAggrSepa, kMipClockModKSepa}; reportMipClockList("MipSeparation", mip_clock_list, mip_timer_clock, - kMipClockTotal); //, tolerance_percent_report); + kMipClockTotal); //, kMipClockTolerancePercentReport); }; void csvMipClock(const std::string model_name, diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 4bfcb146f85..3293dadba6a 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -25,6 +25,7 @@ #include "mip/HighsObjectiveFunction.h" #include "mip/MipTimer.h" #include "presolve/HighsPostsolveStack.h" +#include "presolve/PresolveTimer.h" #include "test_kkt/DevKkt.h" #include "util/HFactor.h" #include "util/HighsCDouble.h" @@ -70,10 +71,16 @@ void HPresolve::debugPrintRow(HighsPostsolveStack& postsolve_stack, bool HPresolve::okSetInput(HighsLp& model_, const HighsOptions& options_, const HighsInt presolve_reduction_limit, HighsTimer* timer) { - model = &model_; - options = &options_; + this->model = &model_; + this->options = &options_; this->timer = timer; + // Set up the logic to allow presolve rules, and logging for their + // effectiveness + analysis_.setup(this->model, this->options, this->numDeletedRows, + this->numDeletedCols, silentLog(), this->timer); + analysis_.presolveTimerStart(kPresolveClockPresolve); + analysis_.presolveTimerStart(kPresolveClockSetupResize); if (!okResize(colLowerSource, model->num_col_, HighsInt{-1})) return false; if (!okResize(colUpperSource, model->num_col_, HighsInt{-1})) return false; if (!okResize(implColLower, model->num_col_, -kHighsInf)) return false; @@ -100,6 +107,9 @@ bool HPresolve::okSetInput(HighsLp& model_, const HighsOptions& options_, } else primal_feastol = options->mip_feasibility_tolerance; + analysis_.presolveTimerStop(kPresolveClockSetupResize); + + analysis_.presolveTimerStart(kPresolveClockSetupToCsc); if (model_.a_matrix_.isRowwise()) { // Does this even happen? assert(model_.a_matrix_.isColwise()); @@ -111,7 +121,9 @@ bool HPresolve::okSetInput(HighsLp& model_, const HighsOptions& options_, model->a_matrix_.start_)) return false; } + analysis_.presolveTimerStop(kPresolveClockSetupToCsc); + analysis_.presolveTimerStart(kPresolveClockSetupResize); // initialize everything as changed, but do not add all indices // since the first thing presolve will do is a scan for easy reductions // of each row and column and set the flag of processed columns to false @@ -126,6 +138,8 @@ bool HPresolve::okSetInput(HighsLp& model_, const HighsOptions& options_, if (!okResize(singleEquationChecked, model->num_row_)) return false; numDeletedCols = 0; numDeletedRows = 0; + analysis_.presolveTimerStop(kPresolveClockSetupResize); + analysis_.presolveTimerStart(kPresolveClockSetupInitialSubstitution); // initialize substitution opportunities for (HighsInt row = 0; row != model->num_row_; ++row) { if (!isDualImpliedFree(row)) continue; @@ -134,6 +148,7 @@ bool HPresolve::okSetInput(HighsLp& model_, const HighsOptions& options_, substitutionOpportunities.emplace_back(row, nonzero.index()); } } + analysis_.presolveTimerStop(kPresolveClockSetupInitialSubstitution); // Take value passed in as reduction limit, allowing different // values to be used for initial presolve, and after restart reductionLimit = @@ -1745,12 +1760,13 @@ HPresolve::Result HPresolve::runProbing(HighsPostsolveStack& postsolve_stack) { // Check for timeout tt = this->timer->read(); if (tt > options->time_limit) { - highsLogUser(options->log_options, HighsLogType::kInfo, - "Time limit reached in probing: " - "consider not using probing by setting option " - "presolve_rule_off to 2^%-d = %d\n", - int(kPresolveRuleProbing), - int(std::pow(int(2), int(kPresolveRuleProbing)))); + highsLogUser( + options->log_options, HighsLogType::kInfo, + "Time limit reached in probing: " + "consider not using probing by setting option " + "presolve_rule_off to 2^%-d = %d\n", + int(kPresolveRuleProbing), + int(std::pow(int(2), static_cast(kPresolveRuleProbing)))); return Result::kStopped; } @@ -3342,7 +3358,7 @@ HPresolve::Result HPresolve::singletonRow(HighsPostsolveStack& postsolve_stack, } HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack, - HighsInt col) { + HighsInt col, const bool timing) { assert(colsize[col] == 1); assert(!colDeleted[col]); HighsInt nzPos = colhead[col]; @@ -3350,17 +3366,26 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack, double colCoef = Avalue[nzPos]; if (rowsize[row] == 1) { + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColSingletonRow); HPRESOLVE_CHECKED_CALL(singletonRow(postsolve_stack, row);); if (!colDeleted[col]) { assert(colsize[col] == 0); - return emptyCol(postsolve_stack, col); + HPresolve::Result result = emptyCol(postsolve_stack, col); + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColSingletonRow); + return result; } + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColSingletonRow); return Result::kOk; } // detect strong / weak domination + if (timing) analysis_.presolveTimerStart(kPresolveClockSingletonColDominated); HPRESOLVE_CHECKED_CALL(detectDominatedCol(postsolve_stack, col, false)); + if (timing) analysis_.presolveTimerStop(kPresolveClockSingletonColDominated); if (colDeleted[col]) return Result::kOk; // check if variable is implied integer @@ -3369,28 +3394,64 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack, static_cast(convertImpliedInteger(col, row))); // dual fixing - HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col)); - if (colDeleted[col]) return Result::kOk; + if (analysis_.allow_rule_[kPresolveRuleDualFixing]) { + const bool logging_on = analysis_.logging_on_; + if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleDualFixing); + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColDualFixing); + HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col)); + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColDualFixing); + analysis_.logging_on_ = logging_on; + if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleDualFixing); + if (colDeleted[col]) return Result::kOk; + } // singleton column stuffing - HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col)); - if (colDeleted[col]) return Result::kOk; + if (analysis_.allow_rule_[kPresolveRuleColStuffing]) { + const bool logging_on = analysis_.logging_on_; + if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleColStuffing); + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColStuffing); + HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col)); + if (timing) analysis_.presolveTimerStop(kPresolveClockSingletonColStuffing); + analysis_.logging_on_ = logging_on; + if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleColStuffing); + if (colDeleted[col]) return Result::kOk; + }; // update column implied bounds + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColImpliedBounds); HPRESOLVE_CHECKED_CALL(updateColImpliedBounds(row, col, colCoef)); + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColImpliedBounds); // update row dual implied bounds - if (model->integrality_[col] != HighsVarType::kInteger) + if (model->integrality_[col] != HighsVarType::kInteger) { + if (timing) + analysis_.presolveTimerStart( + kPresolveClockSingletonColRowDualImpliedBounds); updateRowDualImpliedBounds(row, col, colCoef); - + if (timing) + analysis_.presolveTimerStop( + kPresolveClockSingletonColRowDualImpliedBounds); + } // now check if column is implied free within an equation and substitute the // column if that is the case + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColDualImpliedFree); if (isDualImpliedFree(row) && isImpliedFree(col) && analysis_.allow_rule_[kPresolveRuleFreeColSubstitution]) { if (model->integrality_[col] == HighsVarType::kInteger) { StatusResult impliedIntegral = isImpliedIntegral(col); HPRESOLVE_CHECKED_CALL(static_cast(impliedIntegral)); - if (!impliedIntegral) return Result::kOk; + if (!impliedIntegral) { + if (timing) + analysis_.presolveTimerStop( + kPresolveClockSingletonColDualImpliedFree); + return Result::kOk; + } } const bool logging_on = analysis_.logging_on_; @@ -3406,8 +3467,12 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack, analysis_.logging_on_ = logging_on; if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleFreeColSubstitution); + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColDualImpliedFree); return checkLimits(postsolve_stack); } + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColDualImpliedFree); // todo: check for zero cost singleton and remove return Result::kOk; @@ -4403,7 +4468,7 @@ HPresolve::Result HPresolve::emptyCol(HighsPostsolveStack& postsolve_stack, } HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack, - HighsInt col) { + HighsInt col, const bool timing) { assert(!colDeleted[col]); const bool logging_on = analysis_.logging_on_; @@ -4412,24 +4477,36 @@ HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack, HPRESOLVE_CHECKED_CALL(checkColBounds(col, &isFixed)); if (isFixed) { // remove fixed column + if (timing) analysis_.presolveTimerStart(kPresolveClockInitialColIsFixed); postsolve_stack.removedFixedCol(col, model->col_lower_[col], model->col_cost_[col], getColumnVector(col)); removeFixedCol(col); + if (timing) analysis_.presolveTimerStop(kPresolveClockInitialColIsFixed); return checkLimits(postsolve_stack); } - + HPresolve::Result result; switch (colsize[col]) { case 0: - return emptyCol(postsolve_stack, col); + if (timing) analysis_.presolveTimerStart(kPresolveClockInitialColIsEmpty); + result = emptyCol(postsolve_stack, col); + if (timing) analysis_.presolveTimerStop(kPresolveClockInitialColIsEmpty); + return result; case 1: - return singletonCol(postsolve_stack, col); + if (timing) + analysis_.presolveTimerStart(kPresolveClockInitialColIsSingleton); + result = singletonCol(postsolve_stack, col, timing); + if (timing) + analysis_.presolveTimerStop(kPresolveClockInitialColIsSingleton); + return result; default: break; } // detect strong / weak domination + if (timing) analysis_.presolveTimerStart(kPresolveClockInitialColDominated); HPRESOLVE_CHECKED_CALL(detectDominatedCol(postsolve_stack, col)); + if (timing) analysis_.presolveTimerStop(kPresolveClockInitialColDominated); if (colDeleted[col]) return Result::kOk; // column is not (weakly) dominated @@ -4470,7 +4547,11 @@ HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack, impliedDualRowBounds.getNumInfSumLowerOrig(col)); // check if variable is implied integer + if (timing) + analysis_.presolveTimerStart(kPresolveClockInitialColImpliedInteger); HPRESOLVE_CHECKED_CALL(static_cast(convertImpliedInteger(col))); + if (timing) + analysis_.presolveTimerStop(kPresolveClockInitialColImpliedInteger); // shift integral variables to have a lower bound of zero if (model->integrality_[col] != HighsVarType::kContinuous && @@ -4493,12 +4574,32 @@ HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack, } // dual fixing - HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col)); - if (colDeleted[col]) return Result::kOk; + if (analysis_.allow_rule_[kPresolveRuleDualFixing]) { + const bool logging_on = analysis_.logging_on_; + if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleDualFixing); + if (timing) + analysis_.presolveTimerStart(kPresolveClockSingletonColDualFixing); + HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col)); + if (timing) + analysis_.presolveTimerStop(kPresolveClockSingletonColDualFixing); + analysis_.logging_on_ = logging_on; + if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleDualFixing); + if (colDeleted[col]) return Result::kOk; + } // singleton column stuffing - HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col)); - if (colDeleted[col]) return Result::kOk; + if (analysis_.allow_rule_[kPresolveRuleColStuffing]) { + const bool logging_on = analysis_.logging_on_; + if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleColStuffing); + if (timing) + analysis_.presolveTimerStart(kPresolveClockInitialColSingletonStuffing); + HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col)); + if (timing) + analysis_.presolveTimerStop(kPresolveClockInitialColSingletonStuffing); + analysis_.logging_on_ = logging_on; + if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleColStuffing); + if (colDeleted[col]) return Result::kOk; + } // update dual implied bounds of all rows in given column if (model->integrality_[col] != HighsVarType::kInteger) @@ -5723,22 +5824,27 @@ HPresolve::Result HPresolve::initialRowAndColPresolve( // arrays are not initialized, also unset changedRowFlag so that the row will // be added to the changed row vector when it is changed after it was // processed + analysis_.presolveTimerStart(kPresolveClockInitialRow); for (HighsInt row = 0; row != model->num_row_; ++row) { if (rowDeleted[row]) continue; HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, row)); changedRowFlag[row] = false; } + analysis_.presolveTimerStop(kPresolveClockInitialRow); // same for the columns + analysis_.presolveTimerStart(kPresolveClockInitialCol); + const bool timing = analysis_.analyse_presolve_time_; for (HighsInt col = 0; col != model->num_col_; ++col) { if (colDeleted[col]) continue; // round and update bounds if (model->integrality_[col] != HighsVarType::kContinuous) HPRESOLVE_CHECKED_CALL( changeColBounds(col, model->col_lower_[col], model->col_upper_[col])); - HPRESOLVE_CHECKED_CALL(colPresolve(postsolve_stack, col)); + HPRESOLVE_CHECKED_CALL(colPresolve(postsolve_stack, col, timing)); changedColFlag[col] = false; } + analysis_.presolveTimerStop(kPresolveClockInitialCol); return checkLimits(postsolve_stack); } @@ -5748,15 +5854,25 @@ HPresolve::Result HPresolve::fastPresolveLoop( do { storeCurrentProblemSize(); + analysis_.presolveTimerStart(kPresolveClockFastLoopRowSingletons); HPRESOLVE_CHECKED_CALL(removeRowSingletons(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoopRowSingletons); + analysis_.presolveTimerStart(kPresolveClockFastLoopColSingletons); HPRESOLVE_CHECKED_CALL(presolveChangedRows(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoopColSingletons); + analysis_.presolveTimerStart(kPresolveClockFastLoopDoubletonEquations); HPRESOLVE_CHECKED_CALL(removeDoubletonEquations(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoopDoubletonEquations); + analysis_.presolveTimerStart(kPresolveClockFastLoopChangedRows); HPRESOLVE_CHECKED_CALL(presolveColSingletons(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoopChangedRows); + analysis_.presolveTimerStart(kPresolveClockFastLoopChangedCols); HPRESOLVE_CHECKED_CALL(presolveChangedCols(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoopChangedCols); } while (problemSizeReduction() > 0.01); @@ -5801,15 +5917,11 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { const bool silent = silentLog(); if (options->presolve != kHighsOffString) { - if (!silent) + if (!silent) { highsLogUser(options->log_options, HighsLogType::kInfo, "Presolving model\n"); + } } - // Set up the logic to allow presolve rules, and logging for their - // effectiveness - analysis_.setup(this->model, this->options, this->numDeletedRows, - this->numDeletedCols, silent); - if (options->presolve != kHighsOffString) { if (mipsolver) mipsolver->mipdata_->cliquetable.setPresolveFlag(true); @@ -5836,7 +5948,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { assert(this->timer); assert(this->timer->running()); + analysis_.presolveTimerStart(kPresolveClockInitial); HPRESOLVE_CHECKED_CALL(initialRowAndColPresolve(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockInitial); HighsInt numParallelRowColCalls = 0; // ReductionType::kEqualityRowAddition(s) has no basis postsolve, @@ -5870,7 +5984,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { report(); } + analysis_.presolveTimerStart(kPresolveClockFastLoop); HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoop); storeCurrentProblemSize(); @@ -5883,14 +5999,19 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { applyConflictGraphSubstitutions(postsolve_stack, numDelCol)); } - if (analysis_.allow_rule_[kPresolveRuleAggregator]) + if (analysis_.allow_rule_[kPresolveRuleAggregator]) { + analysis_.presolveTimerStart(kPresolveClockAggregator); HPRESOLVE_CHECKED_CALL(aggregator(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockAggregator); + } if (problemSizeReduction() > 0.05) continue; if (trySparsify && analysis_.allow_rule_[kPresolveRuleSparsify]) { HighsInt numNz = numNonzeros(); + analysis_.presolveTimerStart(kPresolveClockSparsify); HPRESOLVE_CHECKED_CALL(sparsify(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockSparsify); double nzReduction = 100.0 * (1.0 - (numNonzeros() / static_cast(numNz))); @@ -5898,12 +6019,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { highsLogDev(options->log_options, HighsLogType::kInfo, "Sparsify removed %.1f%% of nonzeros\n", nzReduction); - // #1710 exposes that this should not be - // - // fastPresolveLoop(postsolve_stack); - // - // but + analysis_.presolveTimerStart(kPresolveClockFastLoop); HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoop); } trySparsify = false; } @@ -5912,7 +6030,11 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { numParallelRowColCalls < 5) { if (shrinkProblemEnabled && (numDeletedCols >= model->num_col_ / 2 || numDeletedRows >= model->num_row_ / 2)) { + // analysis_.presolveTimerStart(kPresolveClock@); + // analysis_.presolveTimerStop(kPresolveClock@); + analysis_.presolveTimerStart(kPresolveClockShrinkProblem); shrinkProblem(postsolve_stack); + analysis_.presolveTimerStop(kPresolveClockShrinkProblem); toCSC(model->a_matrix_.value_, model->a_matrix_.index_, model->a_matrix_.start_); @@ -5920,12 +6042,16 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { model->a_matrix_.start_); } storeCurrentProblemSize(); + analysis_.presolveTimerStart(kPresolveClockParallelRowsAndCols); HPRESOLVE_CHECKED_CALL(detectParallelRowsAndCols(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockParallelRowsAndCols); ++numParallelRowColCalls; if (problemSizeReduction() > 0.05) continue; } + analysis_.presolveTimerStart(kPresolveClockFastLoop); HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoop); if (mipsolver != nullptr) { HighsInt num_strengthened = -1; @@ -5938,7 +6064,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { num_strengthened); } + analysis_.presolveTimerStart(kPresolveClockFastLoop); HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoop); if (mipsolver != nullptr && numCliquesBeforeProbing == -1) { numCliquesBeforeProbing = mipsolver->mipdata_->cliquetable.numCliques(); @@ -5971,7 +6099,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { if (!dependentEquationsCalled) { if (shrinkProblemEnabled && (numDeletedCols >= model->num_col_ / 2 || numDeletedRows >= model->num_row_ / 2)) { + analysis_.presolveTimerStart(kPresolveClockShrinkProblem); shrinkProblem(postsolve_stack); + analysis_.presolveTimerStop(kPresolveClockShrinkProblem); toCSC(model->a_matrix_.value_, model->a_matrix_.index_, model->a_matrix_.start_); @@ -5980,11 +6110,16 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { } storeCurrentProblemSize(); if (analysis_.allow_rule_[kPresolveRuleDependentEquations]) { + analysis_.presolveTimerStart(kPresolveClockDependentEquations); HPRESOLVE_CHECKED_CALL(removeDependentEquations(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockDependentEquations); dependentEquationsCalled = true; } - if (analysis_.allow_rule_[kPresolveRuleDependentFreeCols]) + if (analysis_.allow_rule_[kPresolveRuleDependentFreeCols]) { + analysis_.presolveTimerStart(kPresolveClockDependentFreeCol); HPRESOLVE_CHECKED_CALL(removeDependentFreeCols(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockDependentFreeCol); + } if (problemSizeReduction() > 0.05) continue; } @@ -6315,7 +6450,7 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) { postsolve_stack.debug_prev_row_upper = 0; // Presolve should only be called with a model that has a non-empty // constraint matrix unless it has no rows - assert(model->a_matrix_.numNz() || model->num_row_ == 0); + assert(model->numNz() || model->num_row_ == 0); auto reportReductions = [&]() { if (options->presolve != kHighsOffString && reductionLimit < kHighsSize_tInf) { @@ -6325,6 +6460,12 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) { postsolve_stack.numReductions(), reductionLimit); } }; + auto reportProfiling = [&]() { + // Presolve profiling not currently enabled for MIP + this->analysis_.presolveTimerStop(kPresolveClockPresolve); + this->analysis_.reportPresolveTimer(); + }; + switch (presolve(postsolve_stack)) { case Result::kStopped: case Result::kOk: @@ -6332,14 +6473,15 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) { case Result::kPrimalInfeasible: presolve_status_ = HighsPresolveStatus::kInfeasible; reportReductions(); + reportProfiling(); return HighsModelStatus::kInfeasible; case Result::kDualInfeasible: presolve_status_ = HighsPresolveStatus::kUnboundedOrInfeasible; reportReductions(); + reportProfiling(); return HighsModelStatus::kUnboundedOrInfeasible; } reportReductions(); - shrinkProblem(postsolve_stack); if (mipsolver != nullptr) { @@ -6392,6 +6534,8 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) { toCSC(model->a_matrix_.value_, model->a_matrix_.index_, model->a_matrix_.start_); + reportProfiling(); + if (model->num_col_ == 0) { // Reduced to empty if (mipsolver) { @@ -6456,19 +6600,30 @@ HPresolve::Result HPresolve::removeDependentEquations( const bool logging_on = analysis_.logging_on_; if (equations.empty()) return Result::kOk; + auto returnOk = [&]() { + analysis_.logging_on_ = logging_on; + if (logging_on) + analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations); + return Result::kOk; + }; + if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleDependentEquations); HighsSparseMatrix matrix; - matrix.num_col_ = equations.size(); + HighsInt num_equations = equations.size(); + matrix.num_col_ = num_equations; matrix.num_row_ = model->num_col_ + 1; - matrix.start_.resize(matrix.num_col_ + 1); + matrix.start_.resize(num_equations + 1); matrix.start_[0] = 0; - const HighsInt maxCapacity = numNonzeros() + matrix.num_col_; + const HighsInt maxCapacity = numNonzeros() + num_equations; matrix.value_.reserve(maxCapacity); matrix.index_.reserve(maxCapacity); - std::vector eqSet(matrix.num_col_); + std::vector eqSet(num_equations); + std::vector row_count; + row_count.assign(model->num_col_, 0); + HighsInt i = 0; for (const std::pair& p : equations) { HighsInt eq = p.second; @@ -6476,8 +6631,10 @@ HPresolve::Result HPresolve::removeDependentEquations( // add entries of equation for (const HighsSliceNonzero& nonz : getRowVector(eq)) { + HighsInt iCol = nonz.index(); + row_count[iCol]++; matrix.value_.push_back(nonz.value()); - matrix.index_.push_back(nonz.index()); + matrix.index_.push_back(iCol); } // add entry for artificial rhs column @@ -6488,7 +6645,22 @@ HPresolve::Result HPresolve::removeDependentEquations( matrix.start_[i] = matrix.value_.size(); } - std::vector colSet(matrix.num_col_); + // Find the number of (true) variables in the system of equations as + // the number of columns with entries in at least one equation + HighsInt num_variables = 0; + for (HighsInt iCol = 0; iCol < model->num_col_; iCol++) + if (row_count[iCol]) num_variables++; + HighsInt num_nz = matrix.numNz(); + const bool silent = silentLog(); + if (!silent) + highsLogUser(options->log_options, HighsLogType::kInfo, + "Considering dependency of %d equation%s in %d variable%s " + "with %d nonzero%s\n", + int(num_equations), highsIntToPlural(num_equations).c_str(), + int(num_variables), highsIntToPlural(num_variables).c_str(), + int(num_nz), highsIntToPlural(num_nz).c_str()); + // Identify any dependent equations + std::vector colSet(num_equations); std::iota(colSet.begin(), colSet.end(), 0); HFactor factor; factor.setup(matrix, colSet); @@ -6501,30 +6673,40 @@ HPresolve::Result HPresolve::removeDependentEquations( // // ToDo: This is strictly non-deterministic, but so conservative // that it'll only reap the cases when factor.build never finishes - const double time_limit = - std::max(1.0, std::min(0.01 * options->time_limit, 1000.0)); + const double kMaxDependentEquationsTime = 100; + const double time_limit = std::max( + 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); - const bool silent = silentLog(); // Determine rank deficiency of the equations if (!silent) highsLogUser(options->log_options, HighsLogType::kInfo, - "Dependent equations search running on %d equations with time " + "Dependent equations search running with time " "limit of %.2fs\n", - static_cast(matrix.num_col_), time_limit); + time_limit); double time_taken = -this->timer->read(); HighsInt build_return = factor.build(); time_taken += this->timer->read(); + // Analyse what's been removed + HighsInt num_removed_row = 0; + HighsInt num_removed_nz = 0; + HighsInt num_fictitious_rows_skipped = 0; if (build_return == kBuildKernelReturnTimeout) { // HFactor::build has timed out, so just return - if (!silent) + if (!silent) { + if (options->log_dev_level > 0) + highsLogUser( + options->log_options, HighsLogType::kInfo, + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g,Terminated\n", + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), static_cast(model->num_col_), + static_cast(num_nz), static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); highsLogUser(options->log_options, HighsLogType::kInfo, "Dependent equations search terminated after %.3gs due to " "expected time exceeding limit\n", time_taken); - analysis_.logging_on_ = logging_on; - if (logging_on) - analysis_.stopPresolveRuleLog(kPresolveRuleDependentFreeCols); - return Result::kOk; + } + return returnOk(); } else { double pct_off_timeout = 1e2 * std::fabs(time_taken - time_limit) / time_limit; @@ -6538,10 +6720,6 @@ HPresolve::Result HPresolve::removeDependentEquations( // build_return as rank_deficiency must be valid assert(build_return >= 0); const HighsInt rank_deficiency = build_return; - // Analyse what's been removed - HighsInt num_removed_row = 0; - HighsInt num_removed_nz = 0; - HighsInt num_fictitious_rows_skipped = 0; for (HighsInt k = 0; k < rank_deficiency; k++) { if (factor.var_with_no_pivot[k] >= 0) { HighsInt redundant_row = eqSet[factor.var_with_no_pivot[k]]; @@ -6553,22 +6731,29 @@ HPresolve::Result HPresolve::removeDependentEquations( num_fictitious_rows_skipped++; } } - if (!silent) - highsLogUser(options->log_options, HighsLogType::kInfo, - "Dependent equations search removed %d rows and %d nonzeros " - "in %.2fs (limit = %.2fs)\n", - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken, time_limit); - if (num_fictitious_rows_skipped) - highsLogDev(options->log_options, HighsLogType::kInfo, - ", avoiding %d fictitious rows", - static_cast(num_fictitious_rows_skipped)); - highsLogDev(options->log_options, HighsLogType::kInfo, "\n"); - - analysis_.logging_on_ = logging_on; - if (logging_on) - analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations); - return Result::kOk; + if (!silent) { + highsLogUser( + options->log_options, HighsLogType::kInfo, + "Search of %d equation%s with %d / %d variable%s and %d nonzero%s " + "removed %d dependent equation%s and %d nonzero%s " + "in %.2fs with bounds in (%.2fs, %.2fs) and limit = %.2fs", + // clang-format off + static_cast(num_equations), highsIntToPlural(num_equations).c_str(), + static_cast(num_variables), + static_cast(model->num_col_), highsIntToPlural(num_variables).c_str(), + static_cast(num_nz), highsIntToPlural(num_nz).c_str(), + static_cast(num_removed_row), highsIntToPlural(num_removed_row).c_str(), + static_cast(num_removed_nz), highsIntToPlural(num_removed_nz).c_str(), + // clang-format on + time_taken, factor.min_time_bound_, factor.max_time_bound_, time_limit); + if (num_fictitious_rows_skipped) + highsLogDev(options->log_options, HighsLogType::kInfo, + ", avoiding %d fictitious row%s", + static_cast(num_fictitious_rows_skipped), + highsIntToPlural(num_fictitious_rows_skipped).c_str()); + highsLogUser(options->log_options, HighsLogType::kInfo, "\n"); + } + return returnOk(); } HPresolve::Result HPresolve::removeDependentFreeCols( diff --git a/highs/presolve/HPresolve.h b/highs/presolve/HPresolve.h index 2bec6fa912b..185aaf31b1e 100644 --- a/highs/presolve/HPresolve.h +++ b/highs/presolve/HPresolve.h @@ -400,14 +400,16 @@ class HPresolve { Result emptyCol(HighsPostsolveStack& postsolve_stack, HighsInt col); - Result singletonCol(HighsPostsolveStack& postsolve_stack, HighsInt col); + Result singletonCol(HighsPostsolveStack& postsolve_stack, HighsInt col, + const bool timing = false); void substituteFreeCol(HighsPostsolveStack& postsolve_stack, HighsInt row, HighsInt col, bool relaxRowDualBounds = false); Result rowPresolve(HighsPostsolveStack& postsolve_stack, HighsInt row); - Result colPresolve(HighsPostsolveStack& postsolve_stack, HighsInt col); + Result colPresolve(HighsPostsolveStack& postsolve_stack, HighsInt col, + const bool timing = false); Result detectDominatedCol(HighsPostsolveStack& postsolve_stack, HighsInt col, bool handleSingletonRows = true); diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index f4724ebd056..6ed21b43192 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -7,58 +7,97 @@ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ #include "lp_data/HighsModelUtils.h" #include "presolve/HPresolve.h" +#include "presolve/PresolveTimer.h" void HPresolveAnalysis::setup(const HighsLp* model_, const HighsOptions* options_, const HighsInt& numDeletedRows_, const HighsInt& numDeletedCols_, - const bool silent) { + const bool silent, HighsTimer* timer) { model = model_; options = options_; numDeletedRows = &numDeletedRows_; numDeletedCols = &numDeletedCols_; + timer_ = timer; + analyse_presolve_time_ = + kHighsAnalysisLevelPresolveTime & options->highs_analysis_level && + !model_->isMip(); + if (analyse_presolve_time_) { + HighsTimerClock clock; + clock.timer_pointer_ = timer_; + PresolveTimer presolve_timer; + presolve_timer.initialisePresolveClocks(clock); + presolve_clocks_ = clock; + } + this->allow_rule_.assign(kPresolveRuleCount, true); + std::vector presolve_light_rule_off(kPresolveRuleCount, false); + if (options->presolve_light) { + // Define the rules not used in presolve_light mode + presolve_light_rule_off[kPresolveRuleDependentEquations] = true; + presolve_light_rule_off[kPresolveRuleDependentFreeCols] = true; + presolve_light_rule_off[kPresolveRuleAggregator] = true; + presolve_light_rule_off[kPresolveRuleParallelRowsAndCols] = true; + presolve_light_rule_off[kPresolveRuleSparsify] = true; + presolve_light_rule_off[kPresolveRuleProbing] = true; + presolve_light_rule_off[kPresolveRuleEnumeration] = true; + presolve_light_rule_off[kPresolveRuleDualFixing] = true; + presolve_light_rule_off[kPresolveRuleColStuffing] = true; + } - if (options->presolve_rule_off || options_->log_dev_level) { - // Some presolve rules are off + if (!silent && options_->log_dev_level) { + // State which rules can be off, and what bit to set + highsLogUser(options->log_options, HighsLogType::kInfo, + "Permitted suppression of presolve rules via " + "presolve_rule_off option:\n"); + HighsInt bit = + std::pow(int(2), static_cast(kPresolveRuleFirstAllowOff)); + for (HighsInt rule_type = kPresolveRuleFirstAllowOff; + rule_type < kPresolveRuleCount; rule_type++) { + // This is a rule that can be switched off + highsLogUser(options->log_options, HighsLogType::kInfo, + " Rule %2d (set bit %2d = %6d): %s\n", int(rule_type), + int(rule_type), int(bit), + utilPresolveRuleTypeToString(rule_type).c_str()); + bit *= 2; + } + } + if (options->presolve_rule_off || options->presolve_light) { + // Some presolve rules are off or presolve_light mode is being used // // Transform options->presolve_rule_off into logical settings in // allow_rule_[*], commenting on the rules switched off - if (!silent) { - if (options->presolve_rule_off) { - highsLogUser(options->log_options, HighsLogType::kInfo, - "Presolve rules not allowed:\n"); - } else { - highsLogUser(options->log_options, HighsLogType::kInfo, - "Permitted suppression of presolve rules via " - "presolve_rule_off option:\n"); - } - } + if (!silent) + highsLogUser(options->log_options, HighsLogType::kInfo, + "Presolve rules not allowed:\n"); HighsInt bit = 1; for (HighsInt rule_type = kPresolveRuleMin; rule_type < kPresolveRuleCount; rule_type++) { // Identify whether this rule is allowed - const bool allow = !(options->presolve_rule_off & bit); + const bool rule_off = (options->presolve_rule_off & bit) || + presolve_light_rule_off[rule_type]; if (rule_type >= kPresolveRuleFirstAllowOff) { // This is a rule that can be switched off, so comment // positively if it is off - allow_rule_[rule_type] = allow; - if (!silent) - if (!allow || - (!options->presolve_rule_off && options_->log_dev_level)) - highsLogUser(options->log_options, HighsLogType::kInfo, - " Rule %2d (set bit %2d = %5d): %s\n", - int(rule_type), int(rule_type), int(bit), - utilPresolveRuleTypeToString(rule_type).c_str()); - } else if (!allow && !silent) { + allow_rule_[rule_type] = !rule_off; + if (rule_off && !silent) + highsLogUser(options->log_options, HighsLogType::kInfo, + " Rule %2d (set bit %2d = %6d): %s\n", int(rule_type), + int(rule_type), int(bit), + utilPresolveRuleTypeToString(rule_type).c_str()); + } else if (rule_off) { // This is a rule that cannot be switched off so, if an - // attempt is made, don't allow it to be off and comment - // negatively - highsLogUser(options->log_options, HighsLogType::kWarning, - "Cannot disallow rule %2d (bit %2d = %5d): %s\n", - int(rule_type), int(rule_type), int(bit), - utilPresolveRuleTypeToString(rule_type).c_str()); + // attempt is made, don't allow it to be off and possibly + // comment negatively + if (!silent) + highsLogUser(options->log_options, HighsLogType::kWarning, + "Cannot disallow rule %2d (bit %2d = %5d): %s\n", + int(rule_type), int(rule_type), int(bit), + utilPresolveRuleTypeToString(rule_type).c_str()); + // Check that we're not here because presolve_light mode is + // being used + assert(!presolve_light_rule_off[rule_type]); } bit *= 2; } @@ -237,3 +276,26 @@ bool HPresolveAnalysis::analysePresolveRuleLog(const bool report) { } return true; } + +void HPresolveAnalysis::presolveTimerStart( + const HighsInt presolve_clock) const { + if (!analyse_presolve_time_) return; + HighsInt highs_timer_clock = presolve_clocks_.clock_[presolve_clock]; + presolve_clocks_.timer_pointer_->start(highs_timer_clock); +} + +void HPresolveAnalysis::presolveTimerStop(const HighsInt presolve_clock) const { + if (!analyse_presolve_time_) return; + HighsInt highs_timer_clock = presolve_clocks_.clock_[presolve_clock]; + presolve_clocks_.timer_pointer_->stop(highs_timer_clock); +} + +void HPresolveAnalysis::reportPresolveTimer() { + if (!analyse_presolve_time_) return; + PresolveTimer presolve_timer; + presolve_timer.reportPresolveCoreClock(model->model_name_, presolve_clocks_); + presolve_timer.reportPresolveInitialColPresolveClock(model->model_name_, + presolve_clocks_); + presolve_timer.reportPresolveSingletonColPresolveClock(model->model_name_, + presolve_clocks_); +} diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h index 9fa48f32d26..e61eec96ca9 100644 --- a/highs/presolve/HPresolveAnalysis.h +++ b/highs/presolve/HPresolveAnalysis.h @@ -8,10 +8,16 @@ /**@file presolve/HPresolveAnalysis.h * @brief */ -#ifndef PRESOLVE_HIGHS_PRESOLVE_ANALYSIS_H_ -#define PRESOLVE_HIGHS_PRESOLVE_ANALYSIS_H_ +#ifndef PRESOLVE_HPRESOLVEANALYSIS_H_ +#define PRESOLVE_HPRESOLVEANALYSIS_H_ + +#include "util/HighsTimer.h" class HPresolveAnalysis { + public: + HPresolveAnalysis() : timer_(nullptr), analyse_presolve_time_(false) {} + + HighsTimer* timer_; const HighsLp* model; const HighsOptions* options; const bool* allow_rule; @@ -33,20 +39,27 @@ class HPresolveAnalysis { HighsInt num_deleted_cols0_; HighsPresolveLog presolve_log_; + HighsTimerClock presolve_clocks_; + bool analyse_presolve_time_; + // for LP presolve // // Transform options->presolve_rule_off into logical settings in // allow_rule_[*], commenting on the rules switched off void setup(const HighsLp* model_, const HighsOptions* options_, const HighsInt& numDeletedRows_, const HighsInt& numDeletedCols_, - const bool silent); + const bool silent, HighsTimer* timer); + void setupPresolveTime(const HighsOptions& options); void resetNumDeleted(); std::string presolveReductionTypeToString(const HighsInt reduction_type); void startPresolveRuleLog(const HighsInt rule_type); void stopPresolveRuleLog(const HighsInt rule_type); bool analysePresolveRuleLog(const bool report = false); + void presolveTimerStart(const HighsInt presolve_clock = 0) const; + void presolveTimerStop(const HighsInt presolve_clock = 0) const; + void reportPresolveTimer(); friend class HPresolve; }; -#endif +#endif /* PRESOLVE_HPRESOLVEANALYSIS_H_ */ diff --git a/highs/presolve/PresolveComponent.cpp b/highs/presolve/PresolveComponent.cpp index dcaa16d7364..515b3b92ab9 100644 --- a/highs/presolve/PresolveComponent.cpp +++ b/highs/presolve/PresolveComponent.cpp @@ -37,6 +37,7 @@ HighsPresolveStatus PresolveComponent::run() { } presolve.run(data_.postSolveStack); + data_.presolve_log_ = presolve.getPresolveLog(); presolve_status_ = presolve.getPresolveStatus(); return presolve_status_; diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h new file mode 100644 index 00000000000..9bb8dc04930 --- /dev/null +++ b/highs/presolve/PresolveTimer.h @@ -0,0 +1,246 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* */ +/* This file is part of the HiGHS linear optimization suite */ +/* */ +/* Available as open-source under the MIT License */ +/* */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/**@file presolve/PresolveTimer.h + * @brief Indices of presolve iClocks + */ +#ifndef PRESOLVE_PRESOLVETIMER_H_ +#define PRESOLVE_PRESOLVETIMER_H_ + +// Clocks for profiling presolve +enum iClockPresolve { + kPresolveClockTotal = 0, + kPresolveClockPresolve, + kPresolveClockSetupResize, + kPresolveClockSetupToCsc, + kPresolveClockSetupInitialSubstitution, + kPresolveClockInitial, + kPresolveClockInitialRow, + kPresolveClockInitialCol, + kPresolveClockInitialColIsFixed, + kPresolveClockInitialColIsEmpty, + kPresolveClockInitialColIsSingleton, + kPresolveClockInitialColDominated, + kPresolveClockInitialColImpliedInteger, + kPresolveClockInitialColDualFixing, + kPresolveClockInitialColSingletonStuffing, + kPresolveClockSingletonColSingletonRow, + kPresolveClockSingletonColDominated, + kPresolveClockSingletonColDualFixing, + kPresolveClockSingletonColStuffing, + kPresolveClockSingletonColImpliedBounds, + kPresolveClockSingletonColRowDualImpliedBounds, + kPresolveClockSingletonColDualImpliedFree, + kPresolveClockFastLoop, + kPresolveClockFastLoopRowSingletons, + kPresolveClockFastLoopColSingletons, + kPresolveClockFastLoopDoubletonEquations, + kPresolveClockFastLoopChangedRows, + kPresolveClockFastLoopChangedCols, + kPresolveClockAggregator, + kPresolveClockSparsify, + kPresolveClockParallelRowsAndCols, + kPresolveClockDependentEquations, + kPresolveClockDependentFreeCol, + kPresolveClockShrinkProblem, + // kPresolveClock@, + kNumPresolveClock //!< Number of PRESOLVE clocks +}; + +static const double kPresolveClockTolerancePercentReport = 0.1; + +class PresolveTimer { + public: + void initialisePresolveClocks(HighsTimerClock& presolve_timer_clock) { + HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_; + std::vector& clock = presolve_timer_clock.clock_; + + clock.resize(kNumPresolveClock); + clock[kPresolveClockTotal] = 0; + clock[kPresolveClockPresolve] = timer_pointer->clock_def("Presolve"); + clock[kPresolveClockSetupResize] = + timer_pointer->clock_def("Setup: resize"); + clock[kPresolveClockSetupToCsc] = timer_pointer->clock_def("Setup: to CSC"); + clock[kPresolveClockSetupInitialSubstitution] = + timer_pointer->clock_def("Setup: initial substitution"); + clock[kPresolveClockInitial] = timer_pointer->clock_def("Initial"); + clock[kPresolveClockInitialRow] = timer_pointer->clock_def("Initial row"); + clock[kPresolveClockInitialCol] = timer_pointer->clock_def("Initial col"); + clock[kPresolveClockInitialColIsFixed] = + timer_pointer->clock_def("I-col: is fixed"); + clock[kPresolveClockInitialColIsEmpty] = + timer_pointer->clock_def("I-col: is empty"); + clock[kPresolveClockInitialColIsSingleton] = + timer_pointer->clock_def("I-col: is singleton"); + clock[kPresolveClockInitialColDominated] = + timer_pointer->clock_def("I-col: dominated"); + clock[kPresolveClockInitialColImpliedInteger] = + timer_pointer->clock_def("I-col: implied integer"); + clock[kPresolveClockInitialColDualFixing] = + timer_pointer->clock_def("I-col: dual fixing"); + clock[kPresolveClockInitialColSingletonStuffing] = + timer_pointer->clock_def("I-col: singleton stuffing"); + clock[kPresolveClockSingletonColSingletonRow] = + timer_pointer->clock_def("S-col: singleton row"); + clock[kPresolveClockSingletonColDominated] = + timer_pointer->clock_def("S-col: dominated"); + clock[kPresolveClockSingletonColDualFixing] = + timer_pointer->clock_def("S-col: dual fixing"); + clock[kPresolveClockSingletonColStuffing] = + timer_pointer->clock_def("S-col: singleton stuffing"); + clock[kPresolveClockSingletonColImpliedBounds] = + timer_pointer->clock_def("S-col: impl bounds"); + clock[kPresolveClockSingletonColRowDualImpliedBounds] = + timer_pointer->clock_def("S-col: row dual impl bounds"); + clock[kPresolveClockSingletonColDualImpliedFree] = + timer_pointer->clock_def("S-col: dual impl free"); + clock[kPresolveClockFastLoop] = timer_pointer->clock_def("Fast loop"); + clock[kPresolveClockFastLoopRowSingletons] = + timer_pointer->clock_def("Fast loop: row singletons"); + clock[kPresolveClockFastLoopColSingletons] = + timer_pointer->clock_def("Fast loop: col singletons"); + clock[kPresolveClockFastLoopDoubletonEquations] = + timer_pointer->clock_def("Fast loop: doubleton equations"); + clock[kPresolveClockFastLoopChangedRows] = + timer_pointer->clock_def("Fast loop: changed rows"); + clock[kPresolveClockFastLoopChangedCols] = + timer_pointer->clock_def("Fast loop: changed cols"); + clock[kPresolveClockAggregator] = timer_pointer->clock_def("Aggregator"); + clock[kPresolveClockSparsify] = timer_pointer->clock_def("Sparsify"); + clock[kPresolveClockParallelRowsAndCols] = + timer_pointer->clock_def("Parallel rows and cols"); + clock[kPresolveClockDependentEquations] = + timer_pointer->clock_def("Dependent equations"); + clock[kPresolveClockDependentFreeCol] = + timer_pointer->clock_def("Dependent free columns"); + clock[kPresolveClockShrinkProblem] = + timer_pointer->clock_def("Shrink problem"); + // clock[kPresolveClock@] = timer_pointer->clock_def("@"); + }; + + bool reportPresolveClockList( + const char* grepStamp, const std::vector presolve_clock_list, + const HighsTimerClock& presolve_timer_clock, + const HighsInt kPresolveClockIdeal = kPresolveClockPresolve, + const double tolerance_percent_report_ = -1) { + HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_; + if (!timer_pointer->printf_flag) return false; + const std::vector& clock = presolve_timer_clock.clock_; + HighsInt presolve_clock_list_size = presolve_clock_list.size(); + std::vector clockList; + clockList.resize(presolve_clock_list_size); + for (HighsInt en = 0; en < presolve_clock_list_size; en++) { + clockList[en] = clock[presolve_clock_list[en]]; + } + const double ideal_sum_time = + timer_pointer->clock_time[clock[kPresolveClockIdeal]]; + const double tolerance_percent_report = + tolerance_percent_report_ >= 0 ? tolerance_percent_report_ : 1e-8; + return timer_pointer->reportOnTolerance( + grepStamp, clockList, ideal_sum_time, tolerance_percent_report); + }; + + void csvPresolveClockList(const std::string& grep_query, + const std::string& model_name, + const std::vector presolve_clock_list, + const HighsTimerClock& presolve_timer_clock, + const HighsInt kPresolveClockIdeal, + const bool header, const bool end_line) { + HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_; + if (!timer_pointer->printf_flag) return; + const std::vector& clock = presolve_timer_clock.clock_; + const double ideal_sum_time = + timer_pointer->clock_time[clock[kPresolveClockIdeal]]; + if (ideal_sum_time < 1e-2) return; + const HighsInt num_clock = presolve_clock_list.size(); + if (header) { + printf("grep_%s,model,ideal", grep_query.c_str()); + for (HighsInt iX = 0; iX < num_clock; iX++) { + HighsInt iclock = clock[presolve_clock_list[iX]]; + printf(",%s", timer_pointer->clock_names[iclock].c_str()); + } + printf(",Unaccounted"); + if (end_line) printf("\n"); + return; + } + double sum_time = 0; + printf("grep_%s,%s,%11.4g", grep_query.c_str(), model_name.c_str(), + ideal_sum_time); + for (HighsInt iX = 0; iX < num_clock; iX++) { + HighsInt iclock = clock[presolve_clock_list[iX]]; + double time = timer_pointer->read(iclock); + sum_time += time; + printf(",%11.4g", time); + } + printf(",%11.4g", ideal_sum_time - sum_time); + if (end_line) printf("\n"); + } + + void reportPresolveCoreClock(const std::string& model_name, + const HighsTimerClock& presolve_timer_clock) { + const std::vector presolve_clock_list{ + kPresolveClockSetupResize, kPresolveClockSetupToCsc, + kPresolveClockSetupInitialSubstitution, + // kPresolveClockInitial, + kPresolveClockInitialRow, kPresolveClockInitialCol, + // kPresolveClockFastLoop, + kPresolveClockFastLoopRowSingletons, + kPresolveClockFastLoopColSingletons, + kPresolveClockFastLoopDoubletonEquations, + kPresolveClockFastLoopChangedRows, kPresolveClockFastLoopChangedCols, + kPresolveClockAggregator, kPresolveClockSparsify, + kPresolveClockParallelRowsAndCols, kPresolveClockDependentEquations, + kPresolveClockDependentFreeCol, kPresolveClockShrinkProblem + // kPresolveClock@ + }; + reportPresolveClockList("PresolveCore_", presolve_clock_list, + presolve_timer_clock, kPresolveClockPresolve, 0.1); + const bool csv_output = false; + if (csv_output) { + csvPresolveClockList("GrepPresolveCore_", model_name, presolve_clock_list, + presolve_timer_clock, kPresolveClockPresolve, true, + true); + csvPresolveClockList("GrepPresolveCore_", model_name, presolve_clock_list, + presolve_timer_clock, kPresolveClockPresolve, false, + true); + } + }; + + void reportPresolveInitialColPresolveClock( + const std::string& model_name, + const HighsTimerClock& presolve_timer_clock) { + const std::vector presolve_clock_list{ + kPresolveClockInitialColIsFixed, + kPresolveClockInitialColIsEmpty, + kPresolveClockInitialColIsSingleton, + kPresolveClockInitialColDominated, + kPresolveClockInitialColImpliedInteger, + kPresolveClockInitialColDualFixing, + kPresolveClockInitialColSingletonStuffing}; + reportPresolveClockList("PresolveInitialCol_", presolve_clock_list, + presolve_timer_clock, kPresolveClockInitialCol, + 0.1); + }; + + void reportPresolveSingletonColPresolveClock( + const std::string& model_name, + const HighsTimerClock& presolve_timer_clock) { + const std::vector presolve_clock_list{ + kPresolveClockSingletonColSingletonRow, + kPresolveClockSingletonColDominated, + kPresolveClockSingletonColDualFixing, + kPresolveClockSingletonColStuffing, + kPresolveClockSingletonColImpliedBounds, + kPresolveClockSingletonColRowDualImpliedBounds, + kPresolveClockSingletonColDualImpliedFree}; + reportPresolveClockList("PresolveSingletonCol_", presolve_clock_list, + presolve_timer_clock, + kPresolveClockInitialColIsSingleton, 0.1); + }; +}; + +#endif /* PRESOLVE_PRESOLVETIMER_H_ */ diff --git a/highs/simplex/HEkk.cpp b/highs/simplex/HEkk.cpp index 9993d20ec97..073db47d9b0 100644 --- a/highs/simplex/HEkk.cpp +++ b/highs/simplex/HEkk.cpp @@ -453,7 +453,7 @@ HighsStatus HEkk::dualize() { assert(lp_.a_matrix_.isColwise()); original_num_col_ = lp_.num_col_; original_num_row_ = lp_.num_row_; - original_num_nz_ = lp_.a_matrix_.numNz(); + original_num_nz_ = lp_.numNz(); original_offset_ = lp_.offset_; original_col_cost_ = lp_.col_cost_; original_col_lower_ = lp_.col_lower_; @@ -959,7 +959,7 @@ HighsStatus HEkk::undualize() { // Some sanity checks assert(lp_.num_col_ == original_num_col_); assert(lp_.num_row_ == original_num_row_); - assert(lp_.a_matrix_.numNz() == original_num_nz_); + assert(lp_.numNz() == original_num_nz_); HighsInt num_basic_variables = primal_basic_index.size(); bool num_basic_variables_ok = num_basic_variables == original_num_row_; if (!num_basic_variables_ok) diff --git a/highs/simplex/HSimplexNla.cpp b/highs/simplex/HSimplexNla.cpp index 4c9d6d6961b..bafdf8056bb 100644 --- a/highs/simplex/HSimplexNla.cpp +++ b/highs/simplex/HSimplexNla.cpp @@ -478,7 +478,7 @@ HighsDebugStatus HSimplexNla::debugCheckData(const std::string message) const { assert(!error_found); return HighsDebugStatus::kLogicalError; } - HighsInt nnz = check_lp.a_matrix_.numNz(); + HighsInt nnz = check_lp.numNz(); HighsInt error_el = -1; for (HighsInt iEl = 0; iEl < nnz; iEl++) { if (check_lp.a_matrix_.index_[iEl] != factor_Aindex[iEl]) { diff --git a/highs/simplex/HighsSimplexAnalysis.cpp b/highs/simplex/HighsSimplexAnalysis.cpp index 1d24b7b2fba..1e9f1aba809 100644 --- a/highs/simplex/HighsSimplexAnalysis.cpp +++ b/highs/simplex/HighsSimplexAnalysis.cpp @@ -1227,10 +1227,10 @@ void HighsSimplexAnalysis::updateInvertFormData(const HFactor& factor) { HighsInt kernel_invert_num_el = factor.invert_num_el - - (factor.basis_matrix_num_el - factor.kernel_num_el); + (factor.basis_matrix_num_el - factor.kernel_num_el) - factor.kernel_dim; assert(factor.kernel_num_el); - double kernel_fill_factor = - (1.0 * kernel_invert_num_el) / factor.kernel_num_el; + double kernel_fill_factor = (1.0 * kernel_invert_num_el) / + (factor.kernel_num_el + factor.kernel_dim); sum_kernel_fill_factor += kernel_fill_factor; running_average_kernel_fill_factor = 0.95 * running_average_kernel_fill_factor + 0.05 * kernel_fill_factor; diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index 31a16af2540..38b02ccc8d7 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -677,7 +677,11 @@ void HFactor::buildSimple() { b_start[iCol + 1] = BcountX; b_var[iCol] = iMat; } - // Record the number of elements in the basis matrix + // Record the number of elements in the basis matrix, remebering + // that BcountX is the number of entries in the nwork structural + // columns, so have to add num_row - nwork to get the entries in + // logical columns. In particular, if the basis matrix is an + // identity, BcountX = nwork = 0 basis_matrix_num_el = num_row - nwork + BcountX; // count1 = 0; @@ -819,16 +823,20 @@ void HFactor::buildSimple() { row_link_first.assign(num_basic + 1, -1); mr_count.assign(num_row, 0); HighsInt mr_countX = 0; - // Determine the number of entries in the kernel - kernel_num_el = 0; + // Determine the initial number of active nonzeros - to be updated + // as values are eliminated and fill-in or cancellation occur + num_active_nz_ = 0; + HighsInt check_nwork = 0; for (HighsInt iRow = 0; iRow < num_row; iRow++) { HighsInt count = mr_count_before[iRow]; if (count > 0) { + // In the active part of the kernel mr_start[iRow] = mr_countX; mr_space[iRow] = count * 2; mr_countX += count * 2; rlinkAdd(iRow, count); - kernel_num_el += count + 1; + num_active_nz_ += count; + check_nwork++; } } mr_index.resize(mr_countX); @@ -853,9 +861,11 @@ void HFactor::buildSimple() { const HighsInt iRow = b_index[k]; const double value = b_value[k]; if (mr_count_before[iRow] > 0) { + // In the active part of the kernel colInsert(iCol, iRow, value); rowInsert(iCol, iRow); } else { + // Above the active part of the kernel colStoreN(iCol, iRow, value); } } @@ -863,8 +873,11 @@ void HFactor::buildSimple() { clinkAdd(iCol, mc_count_a[iCol]); } build_synthetic_tick += (num_row + nwork + MCcountX) * 40 + mr_countX * 20; - // Record the kernel dimension + // Record the dimension and number of entries in the kernel. kernel_dim = nwork; + kernel_num_el = num_active_nz_; + min_time_bound_ = kHighsInf; + max_time_bound_ = 0; assert((HighsInt)this->refactor_info_.pivot_row.size() == num_basic - nwork); } @@ -875,40 +888,124 @@ HighsInt HFactor::buildKernel() { double fake_fill = 0; double fake_eliminate = 0; + HighsInt search_k = 0; const bool progress_report = false; // num_basic != num_row; const HighsInt progress_frequency = 10000; + + // Normally this->time_limit_ is kHighsInf, but if there is a finite + // time limit, it implies that buildKernel can bail out. This is + // (currently) used only with the dependent equations rule in + // presolve + // + // ToDo This bail-out is non-deterministic, but the pseudo-clock + // model for INVERT can be used to make it deterministic + const bool check_for_timeout = this->time_limit_ < kHighsInf; + // To know when to bail out of buildKernel a model of how long + // buildKernel will take is needed + // + // If a pivot is deferred, nwork is increased to force another loop, + // so have to count the number of deferred pivots to know how many + // iterations might be needed + HighsInt num_defer_pivot = 0; // Initial timer frequency: may be reduced if iterations get slow + const HighsInt max_timer_frequency = 1000; HighsInt timer_frequency = 100; - double previous_iteration_time = 0; + // Need to maintain an averate iteration time + double previous_iteration_time = build_timer_->read(); double average_iteration_time = 0; - const bool check_for_timeout = this->time_limit_ < kHighsInf; - HighsInt search_k = 0; - + // Need to maintain an average rate of change of the number of + // active nonzeros - in order to predict when it will go to zero + HighsInt previous_num_active_nz = num_active_nz_; + double average_num_active_nz_change_rate = 0; + // Very occasionally, vectors of indices and values have to be moved + // in order to be resized, in which caes the cost of an iteration is + // abnormally high and needs to be excluded from the average + // calculation, so record when it happens + bool resize_data_shift = false; + // Parameters for the running average claculation + double mu0, mu1; + + // Lambda functions that resize HighsInt/double vectors, with a + // check for them being moved as a result + auto resizeHighsInt = [&](std::vector& i_vector, + const HighsInt to_size) { + HighsInt* from_p = i_vector.data(); + i_vector.resize(to_size); + if (i_vector.data() != from_p) resize_data_shift = true; + }; + + auto resizeDouble = [&](std::vector& d_vector, + const HighsInt to_size) { + double* from_p = d_vector.data(); + d_vector.resize(to_size); + if (d_vector.data() != from_p) resize_data_shift = true; + }; + + // Work out the parameters for the running average claculation + auto runningAverageMu = [&]() { + mu0 = (1.0 * timer_frequency) / (10.0 * max_timer_frequency); + assert(mu0 <= 0.2); + mu1 = 1.0 - mu0; + }; + + runningAverageMu(); const HighsInt check_nwork = -11; while (nwork-- > 0) { // printf("\nnwork = %d\n", (int)nwork); if (nwork == check_nwork) { reportAsm(); } - // Determine whether to return due to exceeding the time limit - if (check_for_timeout && search_k % timer_frequency == 0) { + // Determine whether to return due to (expecting to) exceed the + // time limit + if (check_for_timeout && search_k > 0 && search_k % timer_frequency == 0) { + // Get the current iteration time, and update the running + // average (if there has been no data shift on resize) double current_time = build_timer_->read(); double time_difference = current_time - previous_iteration_time; previous_iteration_time = current_time; double iteration_time = time_difference / (1.0 * timer_frequency); - average_iteration_time = - 0.9 * average_iteration_time + 0.1 * iteration_time; - - if (time_difference > this->time_limit_ / 1e3) - timer_frequency = std::max(HighsInt(1), timer_frequency / 10); - HighsInt iterations_left = kernel_dim - search_k + 1; + if (!resize_data_shift) { + average_iteration_time = + mu1 * average_iteration_time + mu0 * iteration_time; + // Determine whether to reduce the timer frequency + if (timer_frequency > 1 && time_difference > this->time_limit_ / 1e1) { + timer_frequency = std::max(HighsInt(1), timer_frequency / 10); + runningAverageMu(); + } + } + // Determine the current rate of change of the number of active + // nonzeros + double num_active_nz_change_rate = + static_cast(num_active_nz_ - previous_num_active_nz) / + static_cast(timer_frequency); + previous_num_active_nz = num_active_nz_; + average_num_active_nz_change_rate = + mu1 * average_num_active_nz_change_rate + + mu0 * num_active_nz_change_rate; + + // Get an estimate of the number of iterations left, based on + // the bound and the average rate of change of the number of + // active nonzeros (if negative) + HighsInt iterations_left = kernel_dim - search_k + num_defer_pivot + 1; + if (average_num_active_nz_change_rate < -1) { + HighsInt active_nz_iterations_left = + -num_active_nz_ / average_num_active_nz_change_rate; + iterations_left = std::min(active_nz_iterations_left, iterations_left); + } double remaining_time_bound = average_iteration_time * iterations_left; double total_time_bound = current_time + remaining_time_bound; + // Update the record of bounds on total time - for logging in + // presolve + min_time_bound_ = std::min(total_time_bound, min_time_bound_); + max_time_bound_ = std::max(total_time_bound, max_time_bound_); + // Bail out if current or expected time exceeds the limit if (current_time > this->time_limit_ || total_time_bound > this->time_limit_) return kBuildKernelReturnTimeout; + // Clear the record of vectors of indices and values having to + // be moved in order to be resized + resize_data_shift = false; } - /** * 1. Search for the pivot */ @@ -1049,7 +1146,6 @@ HighsInt HFactor::buildKernel() { (int)rank_deficiency); return rank_deficiency; } - /** * 2. Elimination other elements by the pivot */ @@ -1061,6 +1157,9 @@ HighsInt HFactor::buildKernel() { // // Remove the pivot row index from the pivotal column of the // col-wise matrix. Also decreases the column count + // + // One active nonzero is lost + num_active_nz_--; double pivot_multiplier = colDelete(jColPivot, iRowPivot); // Remove the pivot column index from the pivotal row of the // row-wise matrix. Also decreases the row count @@ -1078,6 +1177,7 @@ HighsInt HFactor::buildKernel() { "Defer singular pivot = %11.4g\n", pivot_multiplier); // Matrix is singular, but defer return since other valid pivots // may exist. + num_defer_pivot++; assert(mr_count[iRowPivot] == original_pivotal_row_count - 1); if (mr_count[iRowPivot] == 0) { // The pivot corresponds to a singleton row. Entry is zeroed, @@ -1123,6 +1223,8 @@ HighsInt HFactor::buildKernel() { mr_count_before[iRow] = mr_count[iRow]; rowDelete(jColPivot, (int)iRow); } + // One active entry is lost for each entry in the pivotal column + num_active_nz_ -= (end_A - start_A); l_start.push_back(l_index.size()); fake_fill += 2 * mc_count_a[jColPivot]; @@ -1149,6 +1251,8 @@ HighsInt HFactor::buildKernel() { const HighsInt my_end = my_start + my_count - 1; double my_pivot = colDelete(iCol, iRowPivot); colStoreN(iCol, iRowPivot, my_pivot); + // One active entry is lost for each entry in the pivotal row + num_active_nz_--; // 2.4.2. Elimination on the overlapping part HighsInt nFillin = mwz_column_count; @@ -1167,6 +1271,8 @@ HighsInt HFactor::buildKernel() { mc_value[my_k] = value; } } + // One active entry is lost for each instance of cancellation + num_active_nz_ -= nCancel; fake_eliminate += mwz_column_count; fake_eliminate += nFillin * 2; @@ -1196,8 +1302,8 @@ HighsInt HFactor::buildKernel() { mc_space[iCol] += max(mc_space[iCol], nFillin); HighsInt p5 = mc_start[iCol] = mc_index.size(); HighsInt p7 = p5 + mc_space[iCol] - mc_count_n[iCol]; - mc_index.resize(p5 + mc_space[iCol]); - mc_value.resize(p5 + mc_space[iCol]); + resizeHighsInt(mc_index, p5 + mc_space[iCol]); + resizeDouble(mc_value, p5 + mc_space[iCol]); copy(&mc_index[p1], &mc_index[p2], &mc_index[p5]); copy(&mc_value[p1], &mc_value[p2], &mc_value[p5]); copy(&mc_index[p3], &mc_index[p4], &mc_index[p7]); @@ -1207,8 +1313,11 @@ HighsInt HFactor::buildKernel() { // 2.4.4.2 Fill into column copy for (HighsInt i = 0; i < mwz_column_count; i++) { HighsInt iRow = mwz_column_index[i]; - if (mwz_column_mark[iRow]) + if (mwz_column_mark[iRow]) { colInsert(iCol, iRow, -my_pivot * mwz_column_array[iRow]); + // One active entry is gained for each instance of fill-in + num_active_nz_++; + } } // 2.4.4.3 Fill into the row copy @@ -1221,7 +1330,7 @@ HighsInt HFactor::buildKernel() { HighsInt p2 = p1 + mr_count[iRow]; HighsInt p3 = mr_start[iRow] = mr_index.size(); mr_space[iRow] *= 2; - mr_index.resize(p3 + mr_space[iRow]); + resizeHighsInt(mr_index, p3 + mr_space[iRow]); copy(&mr_index[p1], &mr_index[p2], &mr_index[p3]); } rowInsert(iCol, iRow); @@ -1253,6 +1362,8 @@ HighsInt HFactor::buildKernel() { rlinkAdd(iRow, mr_count[iRow]); } } + // End of loop while(nwork-- > 0) + // Final execution has nwork = 0 } build_synthetic_tick += fake_search * 20 + fake_fill * 160 + fake_eliminate * 80; diff --git a/highs/util/HFactor.h b/highs/util/HFactor.h index 8e452846f89..3dfae4cfabf 100644 --- a/highs/util/HFactor.h +++ b/highs/util/HFactor.h @@ -336,6 +336,9 @@ class HFactor { HighsInt invert_num_el; HighsInt kernel_dim; HighsInt kernel_num_el; + HighsInt num_active_nz_; + double min_time_bound_; + double max_time_bound_; /** * Data of the factor diff --git a/highs/util/HighsTimer.h b/highs/util/HighsTimer.h index 85bb50b7181..4b64273056b 100644 --- a/highs/util/HighsTimer.h +++ b/highs/util/HighsTimer.h @@ -232,6 +232,21 @@ class HighsTimer { clock_time[i_clock] += time; } + /* + void logRunTime(const char* message) const { + if (!printf_flag) return; + double time = this->read(); + std::string time_string = + // std::to_string(time); +#ifndef NDEBUG + std::to_string(time); +#else + std::to_string(static_cast(time)); +#endif + printf("%-30s: %s\n", message, time_string.c_str()); + } + */ + /** * @brief Report timing information for the clock indices in the list */