From fdc149a9b435dbd3fab5ba25d9af361fd1831a0d Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 15:51:54 +0000 Subject: [PATCH 01/17] Make analytic centre calculation an option - allowing true calculation - and count failures --- src/ipm/ipx/conjugate_residuals.cc | 11 +++++++++++ src/ipm/ipx/lp_solver.cc | 20 ++++++++++---------- src/mip/HighsMipSolverData.cpp | 29 +++++++++++++++++++---------- src/mip/HighsMipSolverData.h | 2 ++ 4 files changed, 42 insertions(+), 20 deletions(-) diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index 9fe770e3b9..eb25be7118 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -26,6 +26,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (maxiter < 0) maxiter = m+100; + Int putative_large_iter = 100;//std::max(1000, maxiter/100); + maxiter = putative_large_iter; + Int large_iter = 0;//std::max(1000, maxiter/100); + Int report_frequency = 1;//100; + Int report_iter = 0; + // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -59,6 +65,11 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = IPX_ERROR_cr_matrix_not_posdef; break; } + if (iter_ > large_iter && iter_ > report_iter + report_frequency) { + printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d\n", int(putative_large_iter), int(maxiter), int(iter_)); + report_iter = iter_; + report_frequency *= 2; + } // Update lhs, residual and Cresidual. const double denom = Dot(Cstep,Cstep); diff --git a/src/ipm/ipx/lp_solver.cc b/src/ipm/ipx/lp_solver.cc index 145930c75b..6fb771587d 100644 --- a/src/ipm/ipx/lp_solver.cc +++ b/src/ipm/ipx/lp_solver.cc @@ -65,16 +65,16 @@ Int LpSolver::Solve() { // if ((info_.status_ipm == IPX_STATUS_optimal || // info_.status_ipm == IPX_STATUS_imprecise) && run_crossover_on) { if (run_crossover) { - if (run_crossover_on) { - control_.hLog("Running crossover as requested\n"); - } else if (run_crossover_choose) { - assert(info_.status_ipm == IPX_STATUS_imprecise); - control_.hLog("Running crossover since IPX is imprecise\n"); - } else { - assert(run_crossover_on || run_crossover_choose); - } - BuildCrossoverStartingPoint(); - RunCrossover(); + if (run_crossover_on) { + control_.hLog("Running crossover as requested\n"); + } else if (run_crossover_choose) { + assert(info_.status_ipm == IPX_STATUS_imprecise); + control_.hLog("Running crossover since IPX is imprecise\n"); + } else { + assert(run_crossover_on || run_crossover_choose); + } + BuildCrossoverStartingPoint(); + RunCrossover(); } if (basis_) { info_.ftran_sparse = basis_->frac_ftran_sparse(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index e84f9b5396..3dd88092c9 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -318,7 +318,8 @@ void HighsMipSolverData::startAnalyticCenterComputation( // "_ipm.mps"; printf("Calling ipm.writeModel(%s)\n", // file_name.c_str()); fflush(stdout); ipm.writeModel(file_name); // } - + printf("HighsMipSolverData::startAnalyticCenterComputation submip = %d analyticCenterFailed = %d; num_Col = %d; num_row = %d\n", + mipsolver.submip, analyticCenterFailed, int(ipm.getNumCol()), int(ipm.getNumRow())); mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp); ipm.run(); mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp); @@ -345,6 +346,7 @@ void HighsMipSolverData::finishAnalyticCenterComputation( fflush(stdout); } analyticCenterComputed = true; + analyticCenterFailed = analyticCenterStatus == HighsModelStatus::kSolveError; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -611,6 +613,7 @@ void HighsMipSolverData::init() { firstlpsolobj = -kHighsInf; rootlpsolobj = -kHighsInf; analyticCenterComputed = false; + analyticCenterFailed = false; analyticCenterStatus = HighsModelStatus::kNotset; maxTreeSizeLog2 = 0; numRestarts = 0; @@ -963,7 +966,10 @@ void HighsMipSolverData::runSetup() { } #endif - if (upper_limit == kHighsInf) analyticCenterComputed = false; + if (upper_limit == kHighsInf) { + analyticCenterComputed = false; + analyticCenterFailed = false; + } analyticCenterStatus = HighsModelStatus::kNotset; analyticCenter.clear(); @@ -1865,10 +1871,13 @@ void HighsMipSolverData::evaluateRootNode() { startSymmetryDetection(tg, symData); mipsolver.analysis_.mipTimerStop(kMipClockStartSymmetryDetection); } - if (compute_analytic_centre && !analyticCenterComputed) { - highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, - "MIP-Timing: %11.2g - starting analytic centre calculation\n", - mipsolver.timer_.read()); + if (compute_analytic_centre && !analyticCenterComputed && !analyticCenterFailed) { + if (mipsolver.analysis_.analyse_mip_time) { + highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, + "MIP-Timing: %11.2g - starting analytic centre calculation\n", + mipsolver.timer_.read()); + fflush(stdout); + } mipsolver.analysis_.mipTimerStart(kMipClockStartAnalyticCentreComputation); startAnalyticCenterComputation(tg); mipsolver.analysis_.mipTimerStop(kMipClockStartAnalyticCentreComputation); @@ -2045,8 +2054,8 @@ void HighsMipSolverData::evaluateRootNode() { mipsolver.analysis_.mipTimerStop(kMipClockSeparation); return; } - if (nseparounds >= 5 && !mipsolver.submip && !analyticCenterComputed && - compute_analytic_centre) { + if (nseparounds >= 5 && !mipsolver.submip && + compute_analytic_centre && !analyticCenterComputed) { if (checkLimits()) { mipsolver.analysis_.mipTimerStop(kMipClockSeparation); return; @@ -2141,7 +2150,7 @@ void HighsMipSolverData::evaluateRootNode() { rootlpsolobj = lp.getObjective(); lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters))); - if (!analyticCenterComputed && compute_analytic_centre) { + if (compute_analytic_centre && !analyticCenterComputed) { if (checkLimits()) return; mipsolver.analysis_.mipTimerStart(kMipClockFinishAnalyticCentreComputation); @@ -2263,7 +2272,7 @@ void HighsMipSolverData::evaluateRootNode() { if (lower_bound <= upper_limit) { if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart && mipsolver.options_mip_->presolve != kHighsOffString) { - if (!analyticCenterComputed && compute_analytic_centre) { + if (compute_analytic_centre && !analyticCenterComputed) { mipsolver.analysis_.mipTimerStart( kMipClockFinishAnalyticCentreComputation); finishAnalyticCenterComputation(tg); diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 38693bf361..8917b8dee2 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -84,6 +84,7 @@ struct HighsMipSolverData { bool cliquesExtracted; bool rowMatrixSet; bool analyticCenterComputed; + bool analyticCenterFailed; HighsModelStatus analyticCenterStatus; bool detectSymmetries; HighsInt numRestarts; @@ -170,6 +171,7 @@ struct HighsMipSolverData { cliquesExtracted(false), rowMatrixSet(false), analyticCenterComputed(false), + analyticCenterFailed(false), analyticCenterStatus(HighsModelStatus::kNotset), detectSymmetries(false), numRestarts(0), From 2eb63d3df1ec6b26adac5f75729dcd0550459da3 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 16:26:30 +0000 Subject: [PATCH 02/17] Now not attempring AC calculation in sub-MIP if it failed (or was not attempted) in parent MIP --- src/lp_data/HighsOptions.h | 8 ++++++++ src/mip/HighsMipSolver.cpp | 11 +++++++---- src/mip/HighsMipSolverData.cpp | 13 +++++++------ src/mip/HighsMipSolverData.h | 2 ++ src/mip/HighsPrimalHeuristics.cpp | 6 +++++- 5 files changed, 29 insertions(+), 11 deletions(-) diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 82424a0adc..3d6cc8d4f3 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -421,6 +421,7 @@ struct HighsOptionsStruct { HighsInt mip_pool_soft_limit; HighsInt mip_pscost_minreliable; HighsInt mip_min_cliquetable_entries_for_parallelism; + HighsInt mip_compute_analytic_centre; HighsInt mip_report_level; double mip_feasibility_tolerance; double mip_rel_gap; @@ -555,6 +556,7 @@ struct HighsOptionsStruct { mip_pool_soft_limit(0), mip_pscost_minreliable(0), mip_min_cliquetable_entries_for_parallelism(0), + mip_compute_analytic_centre(0), mip_report_level(0), mip_feasibility_tolerance(0.0), mip_rel_gap(0.0), @@ -1032,6 +1034,12 @@ class HighsOptions : public HighsOptionsStruct { kHighsIInf); records.push_back(record_int); + record_int = new OptionRecordInt( + "mip_compute_analytic_centre", + "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 => true", + advanced, &mip_compute_analytic_centre, 0, 1, 2); + records.push_back(record_int); + record_int = new OptionRecordInt("mip_report_level", "MIP solver reporting level", now_advanced, &mip_report_level, 0, 1, 2); diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index 2b5b93e6e6..ead83e320a 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -800,25 +800,28 @@ void HighsMipSolver::cleanupSolve() { " %.2f (presolve)\n" " %.2f (solve)\n" " %.2f (postsolve)\n" - " Max sub-MIP depth %d\n" " Nodes %llu\n" " Repair LPs %llu (%llu feasible; %llu iterations)\n" " LP iterations %llu (total)\n" " %llu (strong br.)\n" " %llu (separation)\n" - " %llu (heuristics)\n", + " %llu (heuristics)\n" + " AC failures %d\n" + " Max sub-MIP depth %d\n", timer_.read(timer_.total_clock), analysis_.mipTimerRead(kMipClockPresolve), analysis_.mipTimerRead(kMipClockSolve), analysis_.mipTimerRead(kMipClockPostsolve), - int(max_submip_level), (long long unsigned)mipdata_->num_nodes, + (long long unsigned)mipdata_->num_nodes, (long long unsigned)mipdata_->total_repair_lp, (long long unsigned)mipdata_->total_repair_lp_feasible, (long long unsigned)mipdata_->total_repair_lp_iterations, (long long unsigned)mipdata_->total_lp_iterations, (long long unsigned)mipdata_->sb_lp_iterations, (long long unsigned)mipdata_->sepa_lp_iterations, - (long long unsigned)mipdata_->heuristic_lp_iterations); + (long long unsigned)mipdata_->heuristic_lp_iterations, + int(mipdata_->num_analytic_centre_failures), + int(max_submip_level > 0 ? max_submip_level : 0)); analysis_.reportMipTimer(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 3dd88092c9..3679b72c89 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -347,6 +347,7 @@ void HighsMipSolverData::finishAnalyticCenterComputation( } analyticCenterComputed = true; analyticCenterFailed = analyticCenterStatus == HighsModelStatus::kSolveError; + if (analyticCenterFailed) num_analytic_centre_failures++; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -636,6 +637,7 @@ void HighsMipSolverData::init() { heuristic_lp_iterations_before_run = 0; sepa_lp_iterations_before_run = 0; sb_lp_iterations_before_run = 0; + num_analytic_centre_failures = 0; num_disp_lines = 0; numCliqueEntriesAfterPresolve = 0; numCliqueEntriesAfterFirstPresolve = 0; @@ -1857,8 +1859,6 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() { } void HighsMipSolverData::evaluateRootNode() { - const bool compute_analytic_centre = true; - if (!compute_analytic_centre) printf("NOT COMPUTING ANALYTIC CENTRE!\n"); HighsInt maxSepaRounds = mipsolver.submip ? 5 : kHighsIInf; if (numRestarts == 0) maxSepaRounds = @@ -1871,7 +1871,8 @@ void HighsMipSolverData::evaluateRootNode() { startSymmetryDetection(tg, symData); mipsolver.analysis_.mipTimerStop(kMipClockStartSymmetryDetection); } - if (compute_analytic_centre && !analyticCenterComputed && !analyticCenterFailed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed && !analyticCenterFailed) { if (mipsolver.analysis_.analyse_mip_time) { highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, "MIP-Timing: %11.2g - starting analytic centre calculation\n", @@ -2055,7 +2056,7 @@ void HighsMipSolverData::evaluateRootNode() { return; } if (nseparounds >= 5 && !mipsolver.submip && - compute_analytic_centre && !analyticCenterComputed) { + mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { if (checkLimits()) { mipsolver.analysis_.mipTimerStop(kMipClockSeparation); return; @@ -2150,7 +2151,7 @@ void HighsMipSolverData::evaluateRootNode() { rootlpsolobj = lp.getObjective(); lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters))); - if (compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { if (checkLimits()) return; mipsolver.analysis_.mipTimerStart(kMipClockFinishAnalyticCentreComputation); @@ -2272,7 +2273,7 @@ void HighsMipSolverData::evaluateRootNode() { if (lower_bound <= upper_limit) { if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart && mipsolver.options_mip_->presolve != kHighsOffString) { - if (compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { mipsolver.analysis_.mipTimerStart( kMipClockFinishAnalyticCentreComputation); finishAnalyticCenterComputation(tg); diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 8917b8dee2..6e101af638 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -120,6 +120,7 @@ struct HighsMipSolverData { HighsInt numintegercols; HighsInt maxTreeSizeLog2; + HighsInt num_analytic_centre_failures; HighsCDouble pruned_treeweight; double avgrootlpiters; double last_disptime; @@ -186,6 +187,7 @@ struct HighsMipSolverData { rootlpsolobj(-kHighsInf), numintegercols(0), maxTreeSizeLog2(0), + num_analytic_centre_failures(0), pruned_treeweight(0), avgrootlpiters(0.0), last_disptime(0.0), diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 5908f4373c..6a9a2c1680 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -132,8 +132,12 @@ bool HighsPrimalHeuristics::solveSubMip( submipoptions.presolve = "on"; submipoptions.mip_detect_symmetry = false; submipoptions.mip_heuristic_effort = 0.8; - // setup solver and run it + // If the analytic centre calculation has failed for the parent MIP, + // then don't perform it for the sub-MIP + if (mipsolver.mipdata_->analyticCenterFailed) + submipoptions.mip_compute_analytic_centre = 0; + // setup solver and run it HighsSolution solution; solution.value_valid = false; solution.dual_valid = false; From b62b112d7e078ff57132ed621d0e241d83181673 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 17:55:08 +0000 Subject: [PATCH 03/17] Introduce kkt_iter_limit --- src/ipm/ipx/conjugate_residuals.cc | 12 ++++++------ src/lp_data/HConst.h | 8 ++++++++ src/lp_data/HighsOptions.h | 5 ++++- src/mip/HighsMipSolverData.cpp | 23 +++++++++++++---------- 4 files changed, 31 insertions(+), 17 deletions(-) diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index eb25be7118..87b6c1e6ca 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -13,6 +13,7 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs) { const Int m = rhs.size(); + printf("ConjugateResiduals::Solve On entry, maxiter = %d\n", int(maxiter)); Vector residual(m); // rhs - C*lhs Vector step(m); // update to lhs Vector Cresidual(m); // C * residual @@ -26,10 +27,9 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (maxiter < 0) maxiter = m+100; - Int putative_large_iter = 100;//std::max(1000, maxiter/100); - maxiter = putative_large_iter; - Int large_iter = 0;//std::max(1000, maxiter/100); - Int report_frequency = 1;//100; + Int large_iter = std::max(1000, maxiter/100); + maxiter = std::min(10*large_iter, m+100); + Int report_frequency = 100; Int report_iter = 0; // Initialize residual, step and Cstep. @@ -65,8 +65,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = IPX_ERROR_cr_matrix_not_posdef; break; } - if (iter_ > large_iter && iter_ > report_iter + report_frequency) { - printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d\n", int(putative_large_iter), int(maxiter), int(iter_)); + if (iter_ >= std::max(large_iter, report_iter + report_frequency)) { + printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d\n", int(large_iter), int(maxiter), int(iter_)); report_iter = iter_; report_frequency *= 2; } diff --git a/src/lp_data/HConst.h b/src/lp_data/HConst.h index 5e6a9d67fc..d830f9a4fc 100644 --- a/src/lp_data/HConst.h +++ b/src/lp_data/HConst.h @@ -45,6 +45,14 @@ const double kExcessivelySmallCostValue = 1e-4; const bool kAllowDeveloperAssert = false; const bool kExtendInvertWhenAddingRows = false; +enum MipAnalyticCentreCalulation { + kMipAnalyticCentreCalulationMin = 0, + kMipAnalyticCentreCalulationNo = kMipAnalyticCentreCalulationMin, + kMipAnalyticCentreCalulationOriginal, + kMipAnalyticCentreCalulationTrue, + kMipAnalyticCentreCalulationMax = kMipAnalyticCentreCalulationTrue +}; + enum class HighsLogType { kInfo = 1, kDetailed, kVerbose, kWarning, kError }; enum SimplexScaleStrategy { diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 3d6cc8d4f3..44886db306 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -1037,7 +1037,10 @@ class HighsOptions : public HighsOptionsStruct { record_int = new OptionRecordInt( "mip_compute_analytic_centre", "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 => true", - advanced, &mip_compute_analytic_centre, 0, 1, 2); + advanced, &mip_compute_analytic_centre, + kMipAnalyticCentreCalulationMin, + kMipAnalyticCentreCalulationOriginal, + kMipAnalyticCentreCalulationMax); records.push_back(record_int); record_int = diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 3679b72c89..f7438d7358 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -302,24 +302,27 @@ void HighsMipSolverData::startAnalyticCenterComputation( taskGroup.spawn([&]() { // first check if the analytic centre computation should be cancelled, e.g. // due to early return in the root node evaluation + assert(mipsolver.options_mip_->mip_compute_analytic_centre >= kMipAnalyticCentreCalulationOriginal); Highs ipm; - ipm.setOptionValue("solver", "ipm"); - ipm.setOptionValue("run_crossover", kHighsOffString); + if (mipsolver.options_mip_->mip_compute_analytic_centre == kMipAnalyticCentreCalulationOriginal) { + // Original calculation is just IPM with crossover off + ipm.setOptionValue("solver", "ipm"); + ipm.setOptionValue("run_crossover", kHighsOffString); + } else { + // True calculation performs centring properly + ipm.setOptionValue("run_centring", true); + ipm.setOptionValue("ipm_optimality_tolerance", 1e-2); + } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - // ipm.setOptionValue("output_flag", !mipsolver.submip); + ipm.setOptionValue("output_flag", !mipsolver.submip); ipm.setOptionValue("ipm_iteration_limit", 200); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); - // if (!mipsolver.submip) { - // const std::string file_name = mipsolver.model_->model_name_ + - // "_ipm.mps"; printf("Calling ipm.writeModel(%s)\n", - // file_name.c_str()); fflush(stdout); ipm.writeModel(file_name); - // } - printf("HighsMipSolverData::startAnalyticCenterComputation submip = %d analyticCenterFailed = %d; num_Col = %d; num_row = %d\n", - mipsolver.submip, analyticCenterFailed, int(ipm.getNumCol()), int(ipm.getNumRow())); + // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); + mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp); ipm.run(); mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp); From 0bb8cf08d2ebef2cdf908f9b83195c9474f11192 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 21:14:47 +0000 Subject: [PATCH 04/17] Now to get data on max CG iterations --- src/ipm/IpxWrapper.cpp | 1 + src/ipm/ipx/conjugate_residuals.cc | 13 ++++++++----- src/ipm/ipx/control.cc | 4 +++- src/ipm/ipx/control.h | 1 + src/ipm/ipx/ipx_internal.h | 1 + src/ipm/ipx/ipx_parameters.h | 1 + src/ipm/ipx/lp_solver.cc | 1 + src/lp_data/HighsOptions.h | 7 +++++++ src/mip/HighsPrimalHeuristics.cpp | 3 +-- 9 files changed, 24 insertions(+), 8 deletions(-) diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index eae2aff189..ed50953852 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -121,6 +121,7 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, parameters.time_limit = options.time_limit - timer.readRunHighsClock(); parameters.ipm_maxiter = options.ipm_iteration_limit - highs_info.ipm_iteration_count; + parameters.kkt_maxiter = options.kkt_iteration_limit; // Determine if crossover is to be run or not // // When doing analytic centring calculations, crossover must not be diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index 87b6c1e6ca..907f1a156b 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -13,7 +13,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs) { const Int m = rhs.size(); - printf("ConjugateResiduals::Solve On entry, maxiter = %d\n", int(maxiter)); Vector residual(m); // rhs - C*lhs Vector step(m); // update to lhs Vector Cresidual(m); // C * residual @@ -27,10 +26,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (maxiter < 0) maxiter = m+100; + Int default_maxiter = maxiter; Int large_iter = std::max(1000, maxiter/100); maxiter = std::min(10*large_iter, m+100); Int report_frequency = 100; Int report_iter = 0; + printf("ConjugateResiduals::Solve Default maxiter = %d; using %d\n", int(default_maxiter), int(maxiter)); // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { @@ -65,8 +66,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = IPX_ERROR_cr_matrix_not_posdef; break; } - if (iter_ >= std::max(large_iter, report_iter + report_frequency)) { - printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d\n", int(large_iter), int(maxiter), int(iter_)); + if (iter_ >= report_iter + report_frequency) { + printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d: time = %g\n", int(large_iter), int(maxiter), int(iter_), control_.Elapsed()); report_iter = iter_; report_frequency *= 2; } @@ -90,8 +91,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, cdot = cdotnew; iter_++; - if ((errflag_ = control_.InterruptCheck()) != 0) - break; + if ((errflag_ = control_.InterruptCheck()) != 0) { + printf("ConjugateResiduals::Solve control_.InterruptCheck()) != 0: time = %g\n", control_.Elapsed()); + break; + } } time_ = timer.Elapsed(); } diff --git a/src/ipm/ipx/control.cc b/src/ipm/ipx/control.cc index 9a149fab15..4b6a3c58e5 100644 --- a/src/ipm/ipx/control.cc +++ b/src/ipm/ipx/control.cc @@ -13,8 +13,10 @@ Control::Control() { Int Control::InterruptCheck(const Int ipm_iteration_count) const { HighsTaskExecutor::getThisWorkerDeque()->checkInterrupt(); if (parameters_.time_limit >= 0.0 && - parameters_.time_limit < timer_.Elapsed()) + parameters_.time_limit < timer_.Elapsed()) { + printf("Control::InterruptCheck Reached time limit of %g\n", parameters_.time_limit); return IPX_ERROR_time_interrupt; + } // The pointer callback_ should not be null, since that indicates // that it's not been set assert(callback_); diff --git a/src/ipm/ipx/control.h b/src/ipm/ipx/control.h index d322b0d383..99e96307f0 100644 --- a/src/ipm/ipx/control.h +++ b/src/ipm/ipx/control.h @@ -69,6 +69,7 @@ class Control { double ipm_drop_primal() const { return parameters_.ipm_drop_primal; } double ipm_drop_dual() const { return parameters_.ipm_drop_dual; } double kkt_tol() const { return parameters_.kkt_tol; } + ipxint kkt_maxiter() const { return parameters_.kkt_maxiter; } ipxint crash_basis() const { return parameters_.crash_basis; } double dependency_tol() const { return parameters_.dependency_tol; } double volume_tol() const { return parameters_.volume_tol; } diff --git a/src/ipm/ipx/ipx_internal.h b/src/ipm/ipx/ipx_internal.h index f51654b2ac..d9a2133099 100644 --- a/src/ipm/ipx/ipx_internal.h +++ b/src/ipm/ipx/ipx_internal.h @@ -32,6 +32,7 @@ struct Parameters : public ipx_parameters { ipm_drop_primal = 1e-9; ipm_drop_dual = 1e-9; kkt_tol = 0.3; + kkt_maxiter = -1; crash_basis = 1; dependency_tol = 1e-6; volume_tol = 2.0; diff --git a/src/ipm/ipx/ipx_parameters.h b/src/ipm/ipx/ipx_parameters.h index 2f9db89186..590f0b90f0 100644 --- a/src/ipm/ipx/ipx_parameters.h +++ b/src/ipm/ipx/ipx_parameters.h @@ -28,6 +28,7 @@ struct ipx_parameters { double ipm_drop_dual; /* Linear solver */ + ipxint kkt_maxiter; double kkt_tol; /* Basis construction in IPM */ diff --git a/src/ipm/ipx/lp_solver.cc b/src/ipm/ipx/lp_solver.cc index 6fb771587d..4ae972c92d 100644 --- a/src/ipm/ipx/lp_solver.cc +++ b/src/ipm/ipx/lp_solver.cc @@ -558,6 +558,7 @@ void LpSolver::RunMainIPM(IPM& ipm) { KKTSolverBasis kkt(control_, *basis_); Timer timer; ipm.maxiter(control_.ipm_maxiter()); + kkt.maxiter(control_.kkt_maxiter()); ipm.Driver(&kkt, iterate_.get(), &info_); info_.time_ipm2 = timer.Elapsed(); } diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 44886db306..7976c0bd8c 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -333,6 +333,7 @@ struct HighsOptionsStruct { // Options for IPM solver HighsInt ipm_iteration_limit; + HighsInt kkt_iteration_limit; // Options for PDLP solver bool pdlp_native_termination; @@ -482,6 +483,7 @@ struct HighsOptionsStruct { output_flag(false), log_to_console(false), ipm_iteration_limit(0), + kkt_iteration_limit(0), pdlp_native_termination(false), pdlp_scaling(false), pdlp_iteration_limit(0), @@ -1082,6 +1084,11 @@ class HighsOptions : public HighsOptionsStruct { &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf); records.push_back(record_int); + record_int = new OptionRecordInt( + "kkt_iteration_limit", "Iteration limit for PCG in IPX IPM solver: default is -1, so set internally", advanced, + &kkt_iteration_limit, -1, -1, kHighsIInf); + records.push_back(record_int); + record_bool = new OptionRecordBool( "pdlp_native_termination", "Use native termination for PDLP solver: Default = false", advanced, diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 6a9a2c1680..3729858ec1 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -135,8 +135,7 @@ bool HighsPrimalHeuristics::solveSubMip( // If the analytic centre calculation has failed for the parent MIP, // then don't perform it for the sub-MIP - if (mipsolver.mipdata_->analyticCenterFailed) - submipoptions.mip_compute_analytic_centre = 0; + // if (mipsolver.mipdata_->analyticCenterFailed) submipoptions.mip_compute_analytic_centre = 0; // setup solver and run it HighsSolution solution; solution.value_valid = false; From 2e180e742e56b43b9a9ad1b09c2dc650a3f5b224 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 22:08:04 +0000 Subject: [PATCH 05/17] changed iter to iterSum in IPX --- src/ipm/ipx/ipm.cc | 2 +- src/ipm/ipx/kkt_solver.cc | 3 ++- src/ipm/ipx/kkt_solver.h | 9 +++++++-- src/ipm/ipx/kkt_solver_basis.cc | 6 ++++-- src/ipm/ipx/kkt_solver_basis.h | 6 ++++-- src/ipm/ipx/kkt_solver_diag.cc | 6 ++++-- src/ipm/ipx/kkt_solver_diag.h | 6 ++++-- src/lp_data/HighsInfo.cpp | 1 + src/lp_data/HighsInfo.h | 6 ++++++ 9 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/ipm/ipx/ipm.cc b/src/ipm/ipx/ipm.cc index 4a21364951..d8860f19fc 100644 --- a/src/ipm/ipx/ipm.cc +++ b/src/ipm/ipx/ipm.cc @@ -856,7 +856,7 @@ void IPM::PrintOutput() { control_.Debug() << " " << Fixed(step_primal_, 4, 2) << " " << Fixed(step_dual_, 4, 2) << " " << Format(kkt_->basis_changes(), 7) - << " " << Format(kkt_->iter(), 7); + << " " << Format(kkt_->iterSum(), 7); control_.Debug() << " " << Format(info_->dual_dropped, 7) << " " << Format(info_->primal_dropped, 7); diff --git a/src/ipm/ipx/kkt_solver.cc b/src/ipm/ipx/kkt_solver.cc index dde0967e52..1aaee27133 100644 --- a/src/ipm/ipx/kkt_solver.cc +++ b/src/ipm/ipx/kkt_solver.cc @@ -16,7 +16,8 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol, info->time_kkt_solve += timer.Elapsed(); } -Int KKTSolver::iter() const { return _iter(); } +Int KKTSolver::iterSum() const { return _iterSum(); } + //Int KKTSolver::iterMax() const { return _iterMax(); } Int KKTSolver::basis_changes() const { return _basis_changes(); } const Basis* KKTSolver::basis() const { return _basis(); } diff --git a/src/ipm/ipx/kkt_solver.h b/src/ipm/ipx/kkt_solver.h index 3889f1b97c..5de82440db 100644 --- a/src/ipm/ipx/kkt_solver.h +++ b/src/ipm/ipx/kkt_solver.h @@ -46,7 +46,12 @@ class KKTSolver { // If an iterative method is used, returns the # iterations in all Solve() // calls since the last call to Factorize(). A direct solver returns the # // iterative refinement steps. - Int iter() const; + Int iterSum() const; + + // If an iterative method is used, returns the max # iterations in + // all Solve() calls since the last call to Factorize(). A direct + // solver returns the max # iterative refinement steps. + // Int iterMax() const; // If a basis matrix is maintained, returns the # basis changes in the last // call to Factorize(). Otherwise returns 0. @@ -60,7 +65,7 @@ class KKTSolver { virtual void _Factorize(Iterate* iterate, Info* info) = 0; virtual void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) = 0; - virtual Int _iter() const = 0; + virtual Int _iterSum() const = 0; virtual Int _basis_changes() const { return 0; } virtual const Basis* _basis() const { return nullptr; } }; diff --git a/src/ipm/ipx/kkt_solver_basis.cc b/src/ipm/ipx/kkt_solver_basis.cc index dddb35ec8b..e335ec9dba 100644 --- a/src/ipm/ipx/kkt_solver_basis.cc +++ b/src/ipm/ipx/kkt_solver_basis.cc @@ -20,7 +20,8 @@ void KKTSolverBasis::_Factorize(Iterate* iterate, Info* info) { const Int n = model_.cols(); info->errflag = 0; factorized_ = false; - iter_ = 0; + iter_sum_ = 0; + iter_max_ = 0; basis_changes_ = 0; for (Int j = 0; j < n+m; j++) @@ -152,7 +153,8 @@ void KKTSolverBasis::_Solve(const Vector& a, const Vector& b, double tol, info->time_cr2_NNt += splitted_normal_matrix_.time_NNt(); info->time_cr2_B += splitted_normal_matrix_.time_B(); info->time_cr2_Bt += splitted_normal_matrix_.time_Bt(); - iter_ += cr.iter(); + iter_sum_ += cr.iter(); + iter_max_ = std::max(cr.iter(), iter_max_); // Permute back solution to normal equations. for (Int k = 0; k < m; k++) diff --git a/src/ipm/ipx/kkt_solver_basis.h b/src/ipm/ipx/kkt_solver_basis.h index db7e1af131..da637324bd 100644 --- a/src/ipm/ipx/kkt_solver_basis.h +++ b/src/ipm/ipx/kkt_solver_basis.h @@ -33,7 +33,8 @@ class KKTSolverBasis : public KKTSolver { void _Factorize(Iterate* iterate, Info* info) override; void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; - Int _iter() const override { return iter_; } + Int _iterSum() const override { return iter_sum_; } + // Int _iterMax() const override { return iter_max_; } Int _basis_changes() const override { return basis_changes_; } const Basis* _basis() const override { return &basis_; } @@ -57,7 +58,8 @@ class KKTSolverBasis : public KKTSolver { Vector colscale_; // interior point column scaling factors bool factorized_{false}; // preconditioner prepared? Int maxiter_{-1}; - Int iter_{0}; + Int iter_sum_{0}; + Int iter_max_{0}; Int basis_changes_{0}; }; diff --git a/src/ipm/ipx/kkt_solver_diag.cc b/src/ipm/ipx/kkt_solver_diag.cc index 382c9d06cc..4ff530cb37 100644 --- a/src/ipm/ipx/kkt_solver_diag.cc +++ b/src/ipm/ipx/kkt_solver_diag.cc @@ -18,7 +18,8 @@ KKTSolverDiag::KKTSolverDiag(const Control& control, const Model& model) : void KKTSolverDiag::_Factorize(Iterate* pt, Info* info) { const Int m = model_.rows(); const Int n = model_.cols(); - iter_ = 0; + iter_sum_ = 0; + iter_max_ = 0; factorized_ = false; if (pt) { @@ -102,7 +103,8 @@ void KKTSolverDiag::_Solve(const Vector& a, const Vector& b, double tol, info->time_cr1 += cr.time(); info->time_cr1_AAt += normal_matrix_.time(); info->time_cr1_pre += precond_.time(); - iter_ += cr.iter(); + iter_sum_ += cr.iter(); + iter_max_ = std::max(cr.iter(), iter_max_); // Recover solution to KKT system. for (Int i = 0; i < m; i++) diff --git a/src/ipm/ipx/kkt_solver_diag.h b/src/ipm/ipx/kkt_solver_diag.h index 31a2b708ba..66029e0a6a 100644 --- a/src/ipm/ipx/kkt_solver_diag.h +++ b/src/ipm/ipx/kkt_solver_diag.h @@ -29,7 +29,8 @@ class KKTSolverDiag : public KKTSolver { void _Factorize(Iterate* iterate, Info* info) override; void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; - Int _iter() const override { return iter_; }; + Int _iterSum() const override { return iter_sum_; }; + // Int _iterMax() const override { return iter_max_; }; const Control& control_; const Model& model_; @@ -40,7 +41,8 @@ class KKTSolverDiag : public KKTSolver { Vector resscale_; // residual scaling factors for CR termination test bool factorized_{false}; // KKT matrix factorized? Int maxiter_{-1}; - Int iter_{0}; // # CR iterations since last Factorize() + Int iter_sum_{0}; + Int iter_max_{0}; }; } // namespace ipx diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp index 2792f64550..987c18cc70 100644 --- a/src/lp_data/HighsInfo.cpp +++ b/src/lp_data/HighsInfo.cpp @@ -41,6 +41,7 @@ void HighsInfo::invalidate() { max_complementarity_violation = kHighsIllegalComplementarityViolation; sum_complementarity_violations = kHighsIllegalComplementarityViolation; primal_dual_integral = -kHighsInf; + max_cr_iteration_count = -1; } static std::string infoEntryTypeToString(const HighsInfoType type) { diff --git a/src/lp_data/HighsInfo.h b/src/lp_data/HighsInfo.h index 3292b8281e..294b4350ed 100644 --- a/src/lp_data/HighsInfo.h +++ b/src/lp_data/HighsInfo.h @@ -116,6 +116,7 @@ struct HighsInfoStruct { double max_complementarity_violation; double sum_complementarity_violations; double primal_dual_integral; + HighsInt max_cr_iteration_count; }; class HighsInfo : public HighsInfoStruct { @@ -279,6 +280,11 @@ class HighsInfo : public HighsInfoStruct { new InfoRecordDouble("primal_dual_integral", "Primal-dual integral", advanced, &primal_dual_integral, 0); records.push_back(record_double); + + record_int = new InfoRecordInt("max_cr_iteration_count", + "Maximum number of CG iterations in IPX", advanced, + &max_cr_iteration_count, -1); + records.push_back(record_int); } public: From ddd52b30beb06f3db57f4ccd6a6f6dfa1484d258 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Mon, 25 Nov 2024 22:53:03 +0000 Subject: [PATCH 06/17] Ready to do some testing to identify limit on KKT iterations --- src/ipm/IpxWrapper.cpp | 6 ++++++ src/ipm/ipx/conjugate_residuals.cc | 13 ------------- src/ipm/ipx/info.cc | 2 ++ src/ipm/ipx/ipx_info.h | 2 ++ src/ipm/ipx/kkt_solver.cc | 2 +- src/ipm/ipx/kkt_solver.h | 6 +++--- src/ipm/ipx/kkt_solver_basis.cc | 2 +- src/ipm/ipx/kkt_solver_basis.h | 2 +- src/ipm/ipx/kkt_solver_diag.cc | 2 +- src/ipm/ipx/kkt_solver_diag.h | 2 +- src/lp_data/HighsInfo.cpp | 3 ++- src/lp_data/HighsInfo.h | 16 +++++++++++---- src/lp_data/HighsOptions.h | 19 +++++++++--------- src/mip/HighsMipSolver.cpp | 6 +++--- src/mip/HighsMipSolverData.cpp | 31 +++++++++++++++++++++--------- src/mip/HighsMipSolverData.h | 4 ++-- src/mip/HighsPrimalHeuristics.cpp | 3 ++- 17 files changed, 71 insertions(+), 50 deletions(-) diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index ed50953852..db48e9c7b7 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -184,6 +184,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, if (report_solve_data) reportSolveData(options.log_options, ipx_info); highs_info.ipm_iteration_count += (HighsInt)ipx_info.iter; highs_info.crossover_iteration_count += (HighsInt)ipx_info.updates_crossover; + highs_info.max_cr_iteration_count1 = ipx_info.kkt_iter_max1; + highs_info.max_cr_iteration_count2 = ipx_info.kkt_iter_max2; // If not solved... if (solve_status != IPX_STATUS_solved) { @@ -961,6 +963,10 @@ void reportSolveData(const HighsLogOptions& log_options, (int)ipx_info.kktiter1); highsLogDev(log_options, HighsLogType::kInfo, " KKT iter 2 = %d\n", (int)ipx_info.kktiter2); + highsLogDev(log_options, HighsLogType::kInfo, " KKT iter max 1 = %d\n", + (int)ipx_info.kkt_iter_max1); + highsLogDev(log_options, HighsLogType::kInfo, " KKT iter max 2 = %d\n", + (int)ipx_info.kkt_iter_max2); highsLogDev(log_options, HighsLogType::kInfo, " Basis repairs = %d\n", (int)ipx_info.basis_repairs); highsLogDev(log_options, HighsLogType::kInfo, " Updates start = %d\n", diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index 907f1a156b..cb16878743 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -26,13 +26,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (maxiter < 0) maxiter = m+100; - Int default_maxiter = maxiter; - Int large_iter = std::max(1000, maxiter/100); - maxiter = std::min(10*large_iter, m+100); - Int report_frequency = 100; - Int report_iter = 0; - printf("ConjugateResiduals::Solve Default maxiter = %d; using %d\n", int(default_maxiter), int(maxiter)); - // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -66,12 +59,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = IPX_ERROR_cr_matrix_not_posdef; break; } - if (iter_ >= report_iter + report_frequency) { - printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d: time = %g\n", int(large_iter), int(maxiter), int(iter_), control_.Elapsed()); - report_iter = iter_; - report_frequency *= 2; - } - // Update lhs, residual and Cresidual. const double denom = Dot(Cstep,Cstep); const double alpha = cdot/denom; diff --git a/src/ipm/ipx/info.cc b/src/ipm/ipx/info.cc index fa09ab6cff..a64dea6a58 100644 --- a/src/ipm/ipx/info.cc +++ b/src/ipm/ipx/info.cc @@ -53,6 +53,8 @@ std::ostream& operator<<(std::ostream& os, const Info& info) { dump(os, "iter", info.iter); dump(os, "kktiter1", info.kktiter1); dump(os, "kktiter2", info.kktiter2); + dump(os, "kkt_iter_max1", info.kkt_iter_max1); + dump(os, "kkt_iter_max2", info.kkt_iter_max2); dump(os, "basis_repairs", info.basis_repairs); dump(os, "updates_start", info.updates_start); dump(os, "updates_ipm", info.updates_ipm); diff --git a/src/ipm/ipx/ipx_info.h b/src/ipm/ipx/ipx_info.h index ed43d6718e..03cda96090 100644 --- a/src/ipm/ipx/ipx_info.h +++ b/src/ipm/ipx/ipx_info.h @@ -58,6 +58,8 @@ struct ipx_info { ipxint iter; /* # interior point iterations */ ipxint kktiter1; /* # linear solver iterations before switch */ ipxint kktiter2; /* # linear solver iterations after switch */ + ipxint kkt_iter_max1; /* # max linear solver iterations before switch */ + ipxint kkt_iter_max2; /* # max linear solver iterations after switch */ ipxint basis_repairs; /* # basis repairs after crash, < 0 discarded */ ipxint updates_start; /* # basis updates for starting basis */ ipxint updates_ipm; /* # basis updates in IPM */ diff --git a/src/ipm/ipx/kkt_solver.cc b/src/ipm/ipx/kkt_solver.cc index 1aaee27133..a828fcb0a3 100644 --- a/src/ipm/ipx/kkt_solver.cc +++ b/src/ipm/ipx/kkt_solver.cc @@ -17,7 +17,7 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol, } Int KKTSolver::iterSum() const { return _iterSum(); } - //Int KKTSolver::iterMax() const { return _iterMax(); } +Int KKTSolver::iterMax() const { return _iterMax(); } Int KKTSolver::basis_changes() const { return _basis_changes(); } const Basis* KKTSolver::basis() const { return _basis(); } diff --git a/src/ipm/ipx/kkt_solver.h b/src/ipm/ipx/kkt_solver.h index 5de82440db..e5ab15bc22 100644 --- a/src/ipm/ipx/kkt_solver.h +++ b/src/ipm/ipx/kkt_solver.h @@ -49,9 +49,8 @@ class KKTSolver { Int iterSum() const; // If an iterative method is used, returns the max # iterations in - // all Solve() calls since the last call to Factorize(). A direct - // solver returns the max # iterative refinement steps. - // Int iterMax() const; + // _all_ Solve() calls. + Int iterMax() const; // If a basis matrix is maintained, returns the # basis changes in the last // call to Factorize(). Otherwise returns 0. @@ -66,6 +65,7 @@ class KKTSolver { virtual void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) = 0; virtual Int _iterSum() const = 0; + virtual Int _iterMax() const = 0; virtual Int _basis_changes() const { return 0; } virtual const Basis* _basis() const { return nullptr; } }; diff --git a/src/ipm/ipx/kkt_solver_basis.cc b/src/ipm/ipx/kkt_solver_basis.cc index e335ec9dba..46d28e2766 100644 --- a/src/ipm/ipx/kkt_solver_basis.cc +++ b/src/ipm/ipx/kkt_solver_basis.cc @@ -21,7 +21,6 @@ void KKTSolverBasis::_Factorize(Iterate* iterate, Info* info) { info->errflag = 0; factorized_ = false; iter_sum_ = 0; - iter_max_ = 0; basis_changes_ = 0; for (Int j = 0; j < n+m; j++) @@ -149,6 +148,7 @@ void KKTSolverBasis::_Solve(const Vector& a, const Vector& b, double tol, cr.Solve(splitted_normal_matrix_, work, tol, nullptr, maxiter_, lhs); info->errflag = cr.errflag(); info->kktiter2 += cr.iter(); + info->kkt_iter_max2 = std::max(cr.iter(), info->kkt_iter_max2); info->time_cr2 += cr.time(); info->time_cr2_NNt += splitted_normal_matrix_.time_NNt(); info->time_cr2_B += splitted_normal_matrix_.time_B(); diff --git a/src/ipm/ipx/kkt_solver_basis.h b/src/ipm/ipx/kkt_solver_basis.h index da637324bd..b46a89e5c2 100644 --- a/src/ipm/ipx/kkt_solver_basis.h +++ b/src/ipm/ipx/kkt_solver_basis.h @@ -34,7 +34,7 @@ class KKTSolverBasis : public KKTSolver { void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; Int _iterSum() const override { return iter_sum_; } - // Int _iterMax() const override { return iter_max_; } + Int _iterMax() const override { return iter_max_; } Int _basis_changes() const override { return basis_changes_; } const Basis* _basis() const override { return &basis_; } diff --git a/src/ipm/ipx/kkt_solver_diag.cc b/src/ipm/ipx/kkt_solver_diag.cc index 4ff530cb37..9e66ad6848 100644 --- a/src/ipm/ipx/kkt_solver_diag.cc +++ b/src/ipm/ipx/kkt_solver_diag.cc @@ -19,7 +19,6 @@ void KKTSolverDiag::_Factorize(Iterate* pt, Info* info) { const Int m = model_.rows(); const Int n = model_.cols(); iter_sum_ = 0; - iter_max_ = 0; factorized_ = false; if (pt) { @@ -100,6 +99,7 @@ void KKTSolverDiag::_Solve(const Vector& a, const Vector& b, double tol, cr.Solve(normal_matrix_, precond_, rhs, tol, &resscale_[0], maxiter_, y); info->errflag = cr.errflag(); info->kktiter1 += cr.iter(); + info->kkt_iter_max1 = std::max(cr.iter(), info->kkt_iter_max1); info->time_cr1 += cr.time(); info->time_cr1_AAt += normal_matrix_.time(); info->time_cr1_pre += precond_.time(); diff --git a/src/ipm/ipx/kkt_solver_diag.h b/src/ipm/ipx/kkt_solver_diag.h index 66029e0a6a..aaa18cfb85 100644 --- a/src/ipm/ipx/kkt_solver_diag.h +++ b/src/ipm/ipx/kkt_solver_diag.h @@ -30,7 +30,7 @@ class KKTSolverDiag : public KKTSolver { void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; Int _iterSum() const override { return iter_sum_; }; - // Int _iterMax() const override { return iter_max_; }; + Int _iterMax() const override { return iter_max_; }; const Control& control_; const Model& model_; diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp index 987c18cc70..bbb07eb170 100644 --- a/src/lp_data/HighsInfo.cpp +++ b/src/lp_data/HighsInfo.cpp @@ -41,7 +41,8 @@ void HighsInfo::invalidate() { max_complementarity_violation = kHighsIllegalComplementarityViolation; sum_complementarity_violations = kHighsIllegalComplementarityViolation; primal_dual_integral = -kHighsInf; - max_cr_iteration_count = -1; + max_cr_iteration_count1 = -1; + max_cr_iteration_count2 = -1; } static std::string infoEntryTypeToString(const HighsInfoType type) { diff --git a/src/lp_data/HighsInfo.h b/src/lp_data/HighsInfo.h index 294b4350ed..70b56c20ef 100644 --- a/src/lp_data/HighsInfo.h +++ b/src/lp_data/HighsInfo.h @@ -116,7 +116,8 @@ struct HighsInfoStruct { double max_complementarity_violation; double sum_complementarity_violations; double primal_dual_integral; - HighsInt max_cr_iteration_count; + HighsInt max_cr_iteration_count1; + HighsInt max_cr_iteration_count2; }; class HighsInfo : public HighsInfoStruct { @@ -281,9 +282,16 @@ class HighsInfo : public HighsInfoStruct { advanced, &primal_dual_integral, 0); records.push_back(record_double); - record_int = new InfoRecordInt("max_cr_iteration_count", - "Maximum number of CG iterations in IPX", advanced, - &max_cr_iteration_count, -1); + record_int = + new InfoRecordInt("max_cr_iteration_count1", + "Maximum number of diag CR iterations in IPX", + advanced, &max_cr_iteration_count1, -1); + records.push_back(record_int); + + record_int = + new InfoRecordInt("max_cr_iteration_count2", + "Maximum number of basic CR iterations in IPX", + advanced, &max_cr_iteration_count2, -1); records.push_back(record_int); } diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 7976c0bd8c..73c261a889 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -558,7 +558,7 @@ struct HighsOptionsStruct { mip_pool_soft_limit(0), mip_pscost_minreliable(0), mip_min_cliquetable_entries_for_parallelism(0), - mip_compute_analytic_centre(0), + mip_compute_analytic_centre(0), mip_report_level(0), mip_feasibility_tolerance(0.0), mip_rel_gap(0.0), @@ -1038,11 +1038,10 @@ class HighsOptions : public HighsOptionsStruct { record_int = new OptionRecordInt( "mip_compute_analytic_centre", - "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 => true", - advanced, &mip_compute_analytic_centre, - kMipAnalyticCentreCalulationMin, - kMipAnalyticCentreCalulationOriginal, - kMipAnalyticCentreCalulationMax); + "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 " + "=> true", + advanced, &mip_compute_analytic_centre, kMipAnalyticCentreCalulationMin, + kMipAnalyticCentreCalulationOriginal, kMipAnalyticCentreCalulationMax); records.push_back(record_int); record_int = @@ -1084,9 +1083,11 @@ class HighsOptions : public HighsOptionsStruct { &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf); records.push_back(record_int); - record_int = new OptionRecordInt( - "kkt_iteration_limit", "Iteration limit for PCG in IPX IPM solver: default is -1, so set internally", advanced, - &kkt_iteration_limit, -1, -1, kHighsIInf); + record_int = + new OptionRecordInt("kkt_iteration_limit", + "Iteration limit for PCG in IPX IPM solver: " + "default is -1, so set internally", + advanced, &kkt_iteration_limit, -1, -1, kHighsIInf); records.push_back(record_int); record_bool = new OptionRecordBool( diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index ead83e320a..66520a4470 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -806,13 +806,13 @@ void HighsMipSolver::cleanupSolve() { " %llu (strong br.)\n" " %llu (separation)\n" " %llu (heuristics)\n" - " AC failures %d\n" + " AC failures %d\n" " Max sub-MIP depth %d\n", timer_.read(timer_.total_clock), analysis_.mipTimerRead(kMipClockPresolve), analysis_.mipTimerRead(kMipClockSolve), analysis_.mipTimerRead(kMipClockPostsolve), - (long long unsigned)mipdata_->num_nodes, + (long long unsigned)mipdata_->num_nodes, (long long unsigned)mipdata_->total_repair_lp, (long long unsigned)mipdata_->total_repair_lp_feasible, (long long unsigned)mipdata_->total_repair_lp_iterations, @@ -820,7 +820,7 @@ void HighsMipSolver::cleanupSolve() { (long long unsigned)mipdata_->sb_lp_iterations, (long long unsigned)mipdata_->sepa_lp_iterations, (long long unsigned)mipdata_->heuristic_lp_iterations, - int(mipdata_->num_analytic_centre_failures), + int(mipdata_->num_analytic_centre_failures), int(max_submip_level > 0 ? max_submip_level : 0)); analysis_.reportMipTimer(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index f7438d7358..5fe9bff4b9 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -302,9 +302,11 @@ void HighsMipSolverData::startAnalyticCenterComputation( taskGroup.spawn([&]() { // first check if the analytic centre computation should be cancelled, e.g. // due to early return in the root node evaluation - assert(mipsolver.options_mip_->mip_compute_analytic_centre >= kMipAnalyticCentreCalulationOriginal); + assert(mipsolver.options_mip_->mip_compute_analytic_centre >= + kMipAnalyticCentreCalulationOriginal); Highs ipm; - if (mipsolver.options_mip_->mip_compute_analytic_centre == kMipAnalyticCentreCalulationOriginal) { + if (mipsolver.options_mip_->mip_compute_analytic_centre == + kMipAnalyticCentreCalulationOriginal) { // Original calculation is just IPM with crossover off ipm.setOptionValue("solver", "ipm"); ipm.setOptionValue("run_crossover", kHighsOffString); @@ -321,11 +323,18 @@ void HighsMipSolverData::startAnalyticCenterComputation( lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); - // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); + // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + + // "_ipm.mps"); mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp); ipm.run(); mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp); + if (!mipsolver.submip) { + printf("grepAC: model; num_row; max CR counts are, %s, %d, %d, %d\n", + mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), + int(ipm.getInfo().max_cr_iteration_count1), + int(ipm.getInfo().max_cr_iteration_count2)); + } const std::vector& sol = ipm.getSolution().col_value; if (HighsInt(sol.size()) != mipsolver.numCol()) return; analyticCenterStatus = ipm.getModelStatus(); @@ -1877,9 +1886,10 @@ void HighsMipSolverData::evaluateRootNode() { if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed && !analyticCenterFailed) { if (mipsolver.analysis_.analyse_mip_time) { - highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, - "MIP-Timing: %11.2g - starting analytic centre calculation\n", - mipsolver.timer_.read()); + highsLogUser( + mipsolver.options_mip_->log_options, HighsLogType::kInfo, + "MIP-Timing: %11.2g - starting analytic centre calculation\n", + mipsolver.timer_.read()); fflush(stdout); } mipsolver.analysis_.mipTimerStart(kMipClockStartAnalyticCentreComputation); @@ -2059,7 +2069,8 @@ void HighsMipSolverData::evaluateRootNode() { return; } if (nseparounds >= 5 && !mipsolver.submip && - mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { if (checkLimits()) { mipsolver.analysis_.mipTimerStop(kMipClockSeparation); return; @@ -2154,7 +2165,8 @@ void HighsMipSolverData::evaluateRootNode() { rootlpsolobj = lp.getObjective(); lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters))); - if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { if (checkLimits()) return; mipsolver.analysis_.mipTimerStart(kMipClockFinishAnalyticCentreComputation); @@ -2276,7 +2288,8 @@ void HighsMipSolverData::evaluateRootNode() { if (lower_bound <= upper_limit) { if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart && mipsolver.options_mip_->presolve != kHighsOffString) { - if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { mipsolver.analysis_.mipTimerStart( kMipClockFinishAnalyticCentreComputation); finishAnalyticCenterComputation(tg); diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 6e101af638..915d7c9761 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -172,7 +172,7 @@ struct HighsMipSolverData { cliquesExtracted(false), rowMatrixSet(false), analyticCenterComputed(false), - analyticCenterFailed(false), + analyticCenterFailed(false), analyticCenterStatus(HighsModelStatus::kNotset), detectSymmetries(false), numRestarts(0), @@ -187,7 +187,7 @@ struct HighsMipSolverData { rootlpsolobj(-kHighsInf), numintegercols(0), maxTreeSizeLog2(0), - num_analytic_centre_failures(0), + num_analytic_centre_failures(0), pruned_treeweight(0), avgrootlpiters(0.0), last_disptime(0.0), diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 3729858ec1..d3a90ef8b7 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -135,7 +135,8 @@ bool HighsPrimalHeuristics::solveSubMip( // If the analytic centre calculation has failed for the parent MIP, // then don't perform it for the sub-MIP - // if (mipsolver.mipdata_->analyticCenterFailed) submipoptions.mip_compute_analytic_centre = 0; + // if (mipsolver.mipdata_->analyticCenterFailed) + // submipoptions.mip_compute_analytic_centre = 0; // setup solver and run it HighsSolution solution; solution.value_valid = false; From f3ff0b4da0aceeae74dfb606d0a5b3aa83a3f8b8 Mon Sep 17 00:00:00 2001 From: jajhall Date: Tue, 26 Nov 2024 15:53:35 +0000 Subject: [PATCH 07/17] Now to test AC, gathering more data --- src/ipm/ipx/conjugate_residuals.cc | 31 +++++++++++++++++-------- src/mip/HighsMipSolver.cpp | 18 +++++++++++---- src/mip/HighsMipSolverData.cpp | 36 +++++++++++++++++++++++------- src/mip/HighsMipSolverData.h | 10 +++++++-- src/mip/HighsPrimalHeuristics.cpp | 1 + 5 files changed, 73 insertions(+), 23 deletions(-) diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index cb16878743..f7345def7c 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -9,6 +9,8 @@ namespace ipx { ConjugateResiduals::ConjugateResiduals(const Control& control) : control_(control) {} +const bool cr_logging = false; + void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs) { @@ -23,9 +25,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = 0; iter_ = 0; time_ = 0.0; + Int use_maxiter = maxiter; if (maxiter < 0) - maxiter = m+100; + use_maxiter = m+100; + if (cr_logging) printf("\nCR(C) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout); + maxiter = use_maxiter; // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -37,9 +42,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, step = residual; Cstep = Cresidual; + double resnorm; while (true) { // Termination check. - double resnorm = 0.0; + resnorm = 0.0; if (resscale) for (Int i = 0; i < m; i++) resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i])); @@ -48,6 +54,7 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (resnorm <= tol) break; if (iter_ == maxiter) { + if (cr_logging) printf("CR(C) reached maxiter!\n\n"); fflush(stdout); control_.Debug(3) << " CR method not converged in " << maxiter << " iterations." << " residual = " << sci2(resnorm) << ',' @@ -78,11 +85,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, cdot = cdotnew; iter_++; - if ((errflag_ = control_.InterruptCheck()) != 0) { - printf("ConjugateResiduals::Solve control_.InterruptCheck()) != 0: time = %g\n", control_.Elapsed()); - break; - } + if ((errflag_ = control_.InterruptCheck()) != 0) break; } + if (errflag_ != IPX_ERROR_cr_iter_limit) + if (cr_logging) printf("CR(C) iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout); time_ = timer.Elapsed(); } @@ -110,9 +116,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, errflag_ = 0; iter_ = 0; time_ = 0.0; + Int use_maxiter = maxiter; if (maxiter < 0) - maxiter = m+100; + use_maxiter = m+100; + if (cr_logging) printf("\nCR(C,P) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout); + maxiter = use_maxiter; // Initialize residual, sresidual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -125,9 +134,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, step = sresidual; Cstep = Csresidual; + double resnorm; while (true) { // Termination check. - double resnorm = 0.0; + resnorm = 0.0; if (resscale) for (Int i = 0; i < m; i++) resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i])); @@ -136,7 +146,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, if (resnorm <= tol) break; if (iter_ == maxiter) { - control_.Debug(3) + if (cr_logging) printf("CR(C,P) reached maxiter!\n\n"); fflush(stdout); + control_.Debug(3) << " PCR method not converged in " << maxiter << " iterations." << " residual = " << sci2(resnorm) << ',' << " tolerance = " << sci2(tol) << '\n'; @@ -208,6 +219,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, if ((errflag_ = control_.InterruptCheck()) != 0) break; } + if (errflag_ != IPX_ERROR_cr_iter_limit) + if (cr_logging) printf("CR(C,P) iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout); time_ = timer.Elapsed(); } diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index 66520a4470..6585c0ddb5 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -795,6 +795,11 @@ void HighsMipSolver::cleanupSolve() { " %.12g (row viol.)\n", solution_objective_, bound_violation_, integrality_violation_, row_violation_); + int ac_start = mipdata_->num_analytic_centre_start; + int ac_opt = mipdata_->num_analytic_centre_opt; + int ac_other = mipdata_->num_analytic_centre_other; + int ac_fail = mipdata_->num_analytic_centre_fail; + int ac_abort = ac_start - ac_opt - ac_other - ac_fail; highsLogUser(options_mip_->log_options, HighsLogType::kInfo, " Timing %.2f (total)\n" " %.2f (presolve)\n" @@ -806,8 +811,9 @@ void HighsMipSolver::cleanupSolve() { " %llu (strong br.)\n" " %llu (separation)\n" " %llu (heuristics)\n" - " AC failures %d\n" - " Max sub-MIP depth %d\n", + " Max sub-MIP depth %d\n" + " Analytic centre calculations (Optimal / other / abort / " + "fail) = (%d / %d / %d / %d) of %d\n", timer_.read(timer_.total_clock), analysis_.mipTimerRead(kMipClockPresolve), analysis_.mipTimerRead(kMipClockSolve), @@ -820,8 +826,12 @@ void HighsMipSolver::cleanupSolve() { (long long unsigned)mipdata_->sb_lp_iterations, (long long unsigned)mipdata_->sepa_lp_iterations, (long long unsigned)mipdata_->heuristic_lp_iterations, - int(mipdata_->num_analytic_centre_failures), - int(max_submip_level > 0 ? max_submip_level : 0)); + int(max_submip_level > 0 ? max_submip_level : 0), ac_opt, + ac_other, ac_abort, ac_fail, ac_start); + + highsLogUser(options_mip_->log_options, HighsLogType::kInfo, + "grepAcStats,%s,%d,%d,%d,%d,%d\n", model_->model_name_.c_str(), + ac_opt, ac_other, ac_abort, ac_fail, ac_start); analysis_.reportMipTimer(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 5fe9bff4b9..f4af08b640 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -304,6 +304,7 @@ void HighsMipSolverData::startAnalyticCenterComputation( // due to early return in the root node evaluation assert(mipsolver.options_mip_->mip_compute_analytic_centre >= kMipAnalyticCentreCalulationOriginal); + const bool ac_logging = !mipsolver.submip; Highs ipm; if (mipsolver.options_mip_->mip_compute_analytic_centre == kMipAnalyticCentreCalulationOriginal) { @@ -317,8 +318,11 @@ void HighsMipSolverData::startAnalyticCenterComputation( } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - ipm.setOptionValue("output_flag", !mipsolver.submip); + // ipm.setOptionValue("output_flag", !mipsolver.submip); ipm.setOptionValue("ipm_iteration_limit", 200); + HighsInt kkt_iteration_limit = + std::max(100, mipsolver.model_->num_row_ / 1000); + ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); @@ -326,14 +330,28 @@ void HighsMipSolverData::startAnalyticCenterComputation( // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + // "_ipm.mps"); + if (ac_logging) num_analytic_centre_start++; + double tt0 = mipsolver.timer_.read(mipsolver.timer_.total_clock); mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp); ipm.run(); mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp); - if (!mipsolver.submip) { - printf("grepAC: model; num_row; max CR counts are, %s, %d, %d, %d\n", - mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), - int(ipm.getInfo().max_cr_iteration_count1), - int(ipm.getInfo().max_cr_iteration_count2)); + double tt1 = mipsolver.timer_.read(mipsolver.timer_.total_clock); + if (ac_logging) { + HighsModelStatus ac_status = ipm.getModelStatus(); + printf( + "grepAcLocal: model; num_row; CR limit; max CR counts; time and " + "status, %s, %d, %d, %d, %d, %g, %s\n", + mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), + int(kkt_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1), + int(ipm.getInfo().max_cr_iteration_count2), tt1 - tt0, + ipm.modelStatusToString(ac_status).c_str()); + if (ac_status == HighsModelStatus::kSolveError) { + num_analytic_centre_fail++; + } else if (ac_status == HighsModelStatus::kOptimal) { + num_analytic_centre_opt++; + } else { + num_analytic_centre_other++; + } } const std::vector& sol = ipm.getSolution().col_value; if (HighsInt(sol.size()) != mipsolver.numCol()) return; @@ -359,7 +377,6 @@ void HighsMipSolverData::finishAnalyticCenterComputation( } analyticCenterComputed = true; analyticCenterFailed = analyticCenterStatus == HighsModelStatus::kSolveError; - if (analyticCenterFailed) num_analytic_centre_failures++; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -649,7 +666,10 @@ void HighsMipSolverData::init() { heuristic_lp_iterations_before_run = 0; sepa_lp_iterations_before_run = 0; sb_lp_iterations_before_run = 0; - num_analytic_centre_failures = 0; + num_analytic_centre_start = 0; + num_analytic_centre_opt = 0; + num_analytic_centre_other = 0; + num_analytic_centre_fail = 0; num_disp_lines = 0; numCliqueEntriesAfterPresolve = 0; numCliqueEntriesAfterFirstPresolve = 0; diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 915d7c9761..ba8b8c3323 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -120,7 +120,10 @@ struct HighsMipSolverData { HighsInt numintegercols; HighsInt maxTreeSizeLog2; - HighsInt num_analytic_centre_failures; + HighsInt num_analytic_centre_start; + HighsInt num_analytic_centre_opt; + HighsInt num_analytic_centre_other; + HighsInt num_analytic_centre_fail; HighsCDouble pruned_treeweight; double avgrootlpiters; double last_disptime; @@ -187,7 +190,10 @@ struct HighsMipSolverData { rootlpsolobj(-kHighsInf), numintegercols(0), maxTreeSizeLog2(0), - num_analytic_centre_failures(0), + num_analytic_centre_start(0), + num_analytic_centre_opt(0), + num_analytic_centre_other(0), + num_analytic_centre_fail(0), pruned_treeweight(0), avgrootlpiters(0.0), last_disptime(0.0), diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index d3a90ef8b7..91f53a0ea2 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -137,6 +137,7 @@ bool HighsPrimalHeuristics::solveSubMip( // if (mipsolver.mipdata_->analyticCenterFailed) // submipoptions.mip_compute_analytic_centre = 0; + // setup solver and run it HighsSolution solution; solution.value_valid = false; From 150f4f31ab1397197c486ed771a6e1537bffaec5 Mon Sep 17 00:00:00 2001 From: jajhall Date: Tue, 26 Nov 2024 16:04:36 +0000 Subject: [PATCH 08/17] Fixed std::max error in HighsMipSolverDat.cpp:324 --- src/mip/HighsMipSolverData.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index f4af08b640..c1d6c8a45f 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -320,8 +320,8 @@ void HighsMipSolverData::startAnalyticCenterComputation( ipm.setOptionValue("output_flag", false); // ipm.setOptionValue("output_flag", !mipsolver.submip); ipm.setOptionValue("ipm_iteration_limit", 200); - HighsInt kkt_iteration_limit = - std::max(100, mipsolver.model_->num_row_ / 1000); + HighsInt kkt_iteration_limit = mipsolver.model_->num_row_ / 1000; + kkt_iteration_limit = std::max(100, kkt_iteration_limit); ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); From 7feb35417e1afe164485838bb30b76ba2ef2049e Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 26 Nov 2024 19:56:58 +0000 Subject: [PATCH 09/17] Cast 100 to HighsInt in HighsMipSolverData.cpp --- src/mip/HighsMipSolverData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index c1d6c8a45f..b0c0c4bf1c 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -321,7 +321,7 @@ void HighsMipSolverData::startAnalyticCenterComputation( // ipm.setOptionValue("output_flag", !mipsolver.submip); ipm.setOptionValue("ipm_iteration_limit", 200); HighsInt kkt_iteration_limit = mipsolver.model_->num_row_ / 1000; - kkt_iteration_limit = std::max(100, kkt_iteration_limit); + kkt_iteration_limit = std::max(HighsInt(100), kkt_iteration_limit); ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); From 57c5f4cc1512535399a879f87f0e61bd20dae1cf Mon Sep 17 00:00:00 2001 From: JAJHall Date: Wed, 27 Nov 2024 11:24:27 +0000 Subject: [PATCH 10/17] Now to introduce control of kkt iteration limit for phase 1 --- src/ipm/IpxWrapper.cpp | 1 + src/ipm/ipx/conjugate_residuals.cc | 7 +++-- src/ipm/ipx/control.h | 1 + src/ipm/ipx/ipx_internal.h | 1 + src/ipm/ipx/ipx_parameters.h | 1 + src/lp_data/HighsOptions.h | 7 +++++ src/mip/HighsMipSolverData.cpp | 43 +++++++++++++++++++++++------- src/mip/HighsPrimalHeuristics.cpp | 2 ++ 8 files changed, 51 insertions(+), 12 deletions(-) diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index db48e9c7b7..4cd8f1e3ac 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -122,6 +122,7 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, parameters.ipm_maxiter = options.ipm_iteration_limit - highs_info.ipm_iteration_count; parameters.kkt_maxiter = options.kkt_iteration_limit; + parameters.kkt_logging = options.kkt_logging; // Determine if crossover is to be run or not // // When doing analytic centring calculations, crossover must not be diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index f7345def7c..26b89411b1 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -9,11 +9,12 @@ namespace ipx { ConjugateResiduals::ConjugateResiduals(const Control& control) : control_(control) {} -const bool cr_logging = false; - void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs) { + + const bool cr_logging = control_.kkt_logging(); + const Int m = rhs.size(); Vector residual(m); // rhs - C*lhs Vector step(m); // update to lhs @@ -95,6 +96,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs){ + const bool cr_logging = control_.kkt_logging(); + const Int m = rhs.size(); Vector residual(m); // rhs - C*lhs Vector sresidual(m); // preconditioned residual diff --git a/src/ipm/ipx/control.h b/src/ipm/ipx/control.h index 99e96307f0..0ea446b289 100644 --- a/src/ipm/ipx/control.h +++ b/src/ipm/ipx/control.h @@ -68,6 +68,7 @@ class Control { double ipm_optimality_tol() const { return parameters_.ipm_optimality_tol; } double ipm_drop_primal() const { return parameters_.ipm_drop_primal; } double ipm_drop_dual() const { return parameters_.ipm_drop_dual; } + bool kkt_logging() const { return parameters_.kkt_logging; } double kkt_tol() const { return parameters_.kkt_tol; } ipxint kkt_maxiter() const { return parameters_.kkt_maxiter; } ipxint crash_basis() const { return parameters_.crash_basis; } diff --git a/src/ipm/ipx/ipx_internal.h b/src/ipm/ipx/ipx_internal.h index d9a2133099..4b5a152016 100644 --- a/src/ipm/ipx/ipx_internal.h +++ b/src/ipm/ipx/ipx_internal.h @@ -31,6 +31,7 @@ struct Parameters : public ipx_parameters { ipm_optimality_tol = 1e-8; ipm_drop_primal = 1e-9; ipm_drop_dual = 1e-9; + kkt_logging = false; kkt_tol = 0.3; kkt_maxiter = -1; crash_basis = 1; diff --git a/src/ipm/ipx/ipx_parameters.h b/src/ipm/ipx/ipx_parameters.h index 590f0b90f0..388f35e7c8 100644 --- a/src/ipm/ipx/ipx_parameters.h +++ b/src/ipm/ipx/ipx_parameters.h @@ -28,6 +28,7 @@ struct ipx_parameters { double ipm_drop_dual; /* Linear solver */ + bool kkt_logging; ipxint kkt_maxiter; double kkt_tol; diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 73c261a889..c2ffd33e2e 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -334,6 +334,7 @@ struct HighsOptionsStruct { // Options for IPM solver HighsInt ipm_iteration_limit; HighsInt kkt_iteration_limit; + bool kkt_logging; // Options for PDLP solver bool pdlp_native_termination; @@ -484,6 +485,7 @@ struct HighsOptionsStruct { log_to_console(false), ipm_iteration_limit(0), kkt_iteration_limit(0), + kkt_logging(false), pdlp_native_termination(false), pdlp_scaling(false), pdlp_iteration_limit(0), @@ -1083,6 +1085,11 @@ class HighsOptions : public HighsOptionsStruct { &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf); records.push_back(record_int); + record_bool = new OptionRecordBool( + "kkt_logging", "Perform logging in CR solver for KKT in IPX: Default = false", + advanced, &kkt_logging, false); + records.push_back(record_bool); + record_int = new OptionRecordInt("kkt_iteration_limit", "Iteration limit for PCG in IPX IPM solver: " diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index b0c0c4bf1c..aa10e4a2d1 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -318,11 +318,22 @@ void HighsMipSolverData::startAnalyticCenterComputation( } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - // ipm.setOptionValue("output_flag", !mipsolver.submip); + // ipm.setOptionValue("output_flag", !mipsolver.submip); // 2049 unset this ultimately ipm.setOptionValue("ipm_iteration_limit", 200); - HighsInt kkt_iteration_limit = mipsolver.model_->num_row_ / 1000; - kkt_iteration_limit = std::max(HighsInt(100), kkt_iteration_limit); - ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit); + // ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately + // + // kkt_iteration_limit1 is what's set internal to IPX to limit the + // CR iterations before the initial basis is computed, and should + // not be changed + HighsInt kkt_iteration_limit1 = 10 + mipsolver.model_->num_row_ / 20; + kkt_iteration_limit1 = std::max(HighsInt(500), kkt_iteration_limit1); + // kkt_iteration_limit2 is not set internal to IPX, so the default + // value is m+100, which is OK if you're desperate to solve an LP, + // but can make the use of the AC in MIP prohibitively expensive + HighsInt kkt_iteration_limit2 = mipsolver.model_->num_row_ / 1000; + kkt_iteration_limit2 = std::max(HighsInt(100), kkt_iteration_limit2); + // 2049 Set this ultimately + // ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit2); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); @@ -339,13 +350,20 @@ void HighsMipSolverData::startAnalyticCenterComputation( if (ac_logging) { HighsModelStatus ac_status = ipm.getModelStatus(); printf( - "grepAcLocal: model; num_row; CR limit; max CR counts; time and " - "status, %s, %d, %d, %d, %d, %g, %s\n", + "grepAcLocal: model; num_row; CR limits+counts; time and status," + "%s,%d," + "%d,%d," + // "%d,%d," + "%g,%s\n", mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), - int(kkt_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1), - int(ipm.getInfo().max_cr_iteration_count2), tt1 - tt0, + // int(kkt_iteration_limit1), + int(ipm.getInfo().max_cr_iteration_count1), + // int(kkt_iteration_limit2), + int(ipm.getInfo().max_cr_iteration_count2), + tt1 - tt0, ipm.modelStatusToString(ac_status).c_str()); - if (ac_status == HighsModelStatus::kSolveError) { + if (ac_status == HighsModelStatus::kSolveError || + ac_status == HighsModelStatus::kUnknown) { num_analytic_centre_fail++; } else if (ac_status == HighsModelStatus::kOptimal) { num_analytic_centre_opt++; @@ -376,7 +394,9 @@ void HighsMipSolverData::finishAnalyticCenterComputation( fflush(stdout); } analyticCenterComputed = true; - analyticCenterFailed = analyticCenterStatus == HighsModelStatus::kSolveError; + analyticCenterFailed = + analyticCenterStatus == HighsModelStatus::kSolveError || + analyticCenterStatus == HighsModelStatus::kUnknown; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -413,6 +433,9 @@ void HighsMipSolverData::finishAnalyticCenterComputation( nfixed, nintfixed); mipsolver.mipdata_->domain.propagate(); if (mipsolver.mipdata_->domain.infeasible()) return; + } else if (!analyticCenterFailed) { + printf("HighsMipSolverData::finishAnalyticCenterComputation: analyticCenterStatus = %s\n", + lp.getLpSolver().modelStatusToString().c_str()); } } diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 91f53a0ea2..50d6e066f7 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -135,6 +135,8 @@ bool HighsPrimalHeuristics::solveSubMip( // If the analytic centre calculation has failed for the parent MIP, // then don't perform it for the sub-MIP + // 2049 Set this ultimately + // // if (mipsolver.mipdata_->analyticCenterFailed) // submipoptions.mip_compute_analytic_centre = 0; From 39f712d532079a9a11da33bcefcd3d11653d3c43 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Wed, 27 Nov 2024 13:34:07 +0000 Subject: [PATCH 11/17] Test 100 MIP with proposed cr1 and cr2 limits --- src/ipm/IpxWrapper.cpp | 3 ++- src/ipm/ipx/control.h | 3 ++- src/ipm/ipx/ipx_internal.h | 3 ++- src/ipm/ipx/ipx_parameters.h | 3 ++- src/ipm/ipx/lp_solver.cc | 18 +++++++++++++---- src/lp_data/HighsOptions.h | 19 +++++++++++++----- src/mip/HighsMipSolverData.cpp | 33 +++++++++++++++++-------------- src/mip/HighsPrimalHeuristics.cpp | 3 +-- 8 files changed, 55 insertions(+), 30 deletions(-) diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index 4cd8f1e3ac..4a4e3d4c76 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -121,7 +121,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, parameters.time_limit = options.time_limit - timer.readRunHighsClock(); parameters.ipm_maxiter = options.ipm_iteration_limit - highs_info.ipm_iteration_count; - parameters.kkt_maxiter = options.kkt_iteration_limit; + parameters.cr1_maxiter = options.cr1_iteration_limit; + parameters.cr2_maxiter = options.cr2_iteration_limit; parameters.kkt_logging = options.kkt_logging; // Determine if crossover is to be run or not // diff --git a/src/ipm/ipx/control.h b/src/ipm/ipx/control.h index 0ea446b289..c5b19969e0 100644 --- a/src/ipm/ipx/control.h +++ b/src/ipm/ipx/control.h @@ -70,7 +70,8 @@ class Control { double ipm_drop_dual() const { return parameters_.ipm_drop_dual; } bool kkt_logging() const { return parameters_.kkt_logging; } double kkt_tol() const { return parameters_.kkt_tol; } - ipxint kkt_maxiter() const { return parameters_.kkt_maxiter; } + ipxint cr1_maxiter() const { return parameters_.cr1_maxiter; } + ipxint cr2_maxiter() const { return parameters_.cr2_maxiter; } ipxint crash_basis() const { return parameters_.crash_basis; } double dependency_tol() const { return parameters_.dependency_tol; } double volume_tol() const { return parameters_.volume_tol; } diff --git a/src/ipm/ipx/ipx_internal.h b/src/ipm/ipx/ipx_internal.h index 4b5a152016..4fe05ea86f 100644 --- a/src/ipm/ipx/ipx_internal.h +++ b/src/ipm/ipx/ipx_internal.h @@ -33,7 +33,8 @@ struct Parameters : public ipx_parameters { ipm_drop_dual = 1e-9; kkt_logging = false; kkt_tol = 0.3; - kkt_maxiter = -1; + cr1_maxiter = -1; + cr2_maxiter = -1; crash_basis = 1; dependency_tol = 1e-6; volume_tol = 2.0; diff --git a/src/ipm/ipx/ipx_parameters.h b/src/ipm/ipx/ipx_parameters.h index 388f35e7c8..029b89f65c 100644 --- a/src/ipm/ipx/ipx_parameters.h +++ b/src/ipm/ipx/ipx_parameters.h @@ -29,7 +29,8 @@ struct ipx_parameters { /* Linear solver */ bool kkt_logging; - ipxint kkt_maxiter; + ipxint cr1_maxiter; + ipxint cr2_maxiter; double kkt_tol; /* Basis construction in IPM */ diff --git a/src/ipm/ipx/lp_solver.cc b/src/ipm/ipx/lp_solver.cc index 4ae972c92d..0d422e8c65 100644 --- a/src/ipm/ipx/lp_solver.cc +++ b/src/ipm/ipx/lp_solver.cc @@ -484,10 +484,20 @@ void LpSolver::RunInitialIPM(IPM& ipm) { Int switchiter = control_.switchiter(); if (switchiter < 0) { - // Switch iteration not specified by user. Run as long as KKT solver - // converges within min(500,10+m/20) iterations. + // Switch iteration not specified by user. Run as long as KKT + // solver converges within + // + // min(control_.cr1_maxiter(), 500,10+m/20) + // + // iterations, ignoring control_.cr1_maxiter() if negative Int m = model_.rows(); - kkt.maxiter(std::min(500l, (long) (10+m/20) )); + ipxint df_cr1_maxiter = std::min(500l, (long) (10+m/20)); + ipxint cr1_maxiter = df_cr1_maxiter; + ipxint control_cr1_maxiter = control_.cr1_maxiter(); + if (control_cr1_maxiter > 0) cr1_maxiter = std::min(control_cr1_maxiter, cr1_maxiter); + printf("LpSolver::RunInitialIPM Using kkt.maxiter = %d, from df_cr1_maxiter = %d and control_cr1_maxiter = %d\n", + int(cr1_maxiter), int(df_cr1_maxiter), int(control_cr1_maxiter)); + kkt.maxiter(cr1_maxiter); ipm.maxiter(control_.ipm_maxiter()); } else { ipm.maxiter(std::min(switchiter, control_.ipm_maxiter())); @@ -558,7 +568,7 @@ void LpSolver::RunMainIPM(IPM& ipm) { KKTSolverBasis kkt(control_, *basis_); Timer timer; ipm.maxiter(control_.ipm_maxiter()); - kkt.maxiter(control_.kkt_maxiter()); + kkt.maxiter(control_.cr2_maxiter()); ipm.Driver(&kkt, iterate_.get(), &info_); info_.time_ipm2 = timer.Elapsed(); } diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index c2ffd33e2e..359f179da1 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -333,7 +333,8 @@ struct HighsOptionsStruct { // Options for IPM solver HighsInt ipm_iteration_limit; - HighsInt kkt_iteration_limit; + HighsInt cr1_iteration_limit; + HighsInt cr2_iteration_limit; bool kkt_logging; // Options for PDLP solver @@ -484,7 +485,8 @@ struct HighsOptionsStruct { output_flag(false), log_to_console(false), ipm_iteration_limit(0), - kkt_iteration_limit(0), + cr1_iteration_limit(0), + cr2_iteration_limit(0), kkt_logging(false), pdlp_native_termination(false), pdlp_scaling(false), @@ -1091,10 +1093,17 @@ class HighsOptions : public HighsOptionsStruct { records.push_back(record_bool); record_int = - new OptionRecordInt("kkt_iteration_limit", - "Iteration limit for PCG in IPX IPM solver: " + new OptionRecordInt("cr1_iteration_limit", + "Iteration limit for initial CR in IPX IPM solver: " "default is -1, so set internally", - advanced, &kkt_iteration_limit, -1, -1, kHighsIInf); + advanced, &cr1_iteration_limit, -1, -1, kHighsIInf); + records.push_back(record_int); + + record_int = + new OptionRecordInt("cr2_iteration_limit", + "Iteration limit for main CR in IPX IPM solver: " + "default is -1, so set internally", + advanced, &cr2_iteration_limit, -1, -1, kHighsIInf); records.push_back(record_int); record_bool = new OptionRecordBool( diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index aa10e4a2d1..05f05cb46d 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -318,22 +318,25 @@ void HighsMipSolverData::startAnalyticCenterComputation( } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - // ipm.setOptionValue("output_flag", !mipsolver.submip); // 2049 unset this ultimately + ipm.setOptionValue("output_flag", !mipsolver.submip); // 2049 unset this ultimately ipm.setOptionValue("ipm_iteration_limit", 200); - // ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately + ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately // - // kkt_iteration_limit1 is what's set internal to IPX to limit the - // CR iterations before the initial basis is computed, and should - // not be changed - HighsInt kkt_iteration_limit1 = 10 + mipsolver.model_->num_row_ / 20; - kkt_iteration_limit1 = std::max(HighsInt(500), kkt_iteration_limit1); - // kkt_iteration_limit2 is not set internal to IPX, so the default + // cr1_iteration_limit is what's set internal to IPX to limit the + // CR iterations before the initial basis is computed, and perhaps + // should not be changed, but maybe 500 is too high + HighsInt cr1_iteration_limit = 10 + mipsolver.model_->num_row_ / 20; + cr1_iteration_limit = std::min(HighsInt(500), cr1_iteration_limit); // IPX internal default + cr1_iteration_limit = std::min(HighsInt(100), cr1_iteration_limit); //2049 set this ultimately + // 2049 set this ultimately + ipm.setOptionValue("cr1_iteration_limit", cr1_iteration_limit); + // cr2_iteration_limit is not set internal to IPX, so the default // value is m+100, which is OK if you're desperate to solve an LP, // but can make the use of the AC in MIP prohibitively expensive - HighsInt kkt_iteration_limit2 = mipsolver.model_->num_row_ / 1000; - kkt_iteration_limit2 = std::max(HighsInt(100), kkt_iteration_limit2); + HighsInt cr2_iteration_limit = 50 + mipsolver.model_->num_row_ / 1000; + cr2_iteration_limit = std::min(HighsInt(100), cr2_iteration_limit); // 2049 Set this ultimately - // ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit2); + ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); @@ -353,12 +356,12 @@ void HighsMipSolverData::startAnalyticCenterComputation( "grepAcLocal: model; num_row; CR limits+counts; time and status," "%s,%d," "%d,%d," - // "%d,%d," + "%d,%d," "%g,%s\n", mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), - // int(kkt_iteration_limit1), + int(cr1_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1), - // int(kkt_iteration_limit2), + int(cr2_iteration_limit), int(ipm.getInfo().max_cr_iteration_count2), tt1 - tt0, ipm.modelStatusToString(ac_status).c_str()); @@ -435,7 +438,7 @@ void HighsMipSolverData::finishAnalyticCenterComputation( if (mipsolver.mipdata_->domain.infeasible()) return; } else if (!analyticCenterFailed) { printf("HighsMipSolverData::finishAnalyticCenterComputation: analyticCenterStatus = %s\n", - lp.getLpSolver().modelStatusToString().c_str()); + lp.getLpSolver().modelStatusToString(analyticCenterStatus).c_str()); } } diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 50d6e066f7..d05ebefe0b 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -137,8 +137,7 @@ bool HighsPrimalHeuristics::solveSubMip( // 2049 Set this ultimately // - // if (mipsolver.mipdata_->analyticCenterFailed) - // submipoptions.mip_compute_analytic_centre = 0; + if (mipsolver.mipdata_->analyticCenterFailed) submipoptions.mip_compute_analytic_centre = 0; // setup solver and run it HighsSolution solution; From b981b614dcf056810845e25b4af1620df742ab4e Mon Sep 17 00:00:00 2001 From: JAJHall Date: Wed, 27 Nov 2024 23:10:09 +0000 Subject: [PATCH 12/17] Extract current fill from IPX basis --- src/mip/HighsMipSolverData.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 05f05cb46d..0e857d4315 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -333,8 +333,8 @@ void HighsMipSolverData::startAnalyticCenterComputation( // cr2_iteration_limit is not set internal to IPX, so the default // value is m+100, which is OK if you're desperate to solve an LP, // but can make the use of the AC in MIP prohibitively expensive - HighsInt cr2_iteration_limit = 50 + mipsolver.model_->num_row_ / 1000; - cr2_iteration_limit = std::min(HighsInt(100), cr2_iteration_limit); + HighsInt cr2_iteration_limit = 50 + mipsolver.model_->num_row_ / 50; + cr2_iteration_limit = std::min(HighsInt(500), cr2_iteration_limit); // 2049 Set this ultimately ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit); HighsLp lpmodel(*mipsolver.model_); From 45545b755d7ecc5d62fb5bdd13a5b0c672450c1c Mon Sep 17 00:00:00 2001 From: JAJHall Date: Wed, 27 Nov 2024 23:43:27 +0000 Subject: [PATCH 13/17] Now extracting fill factors for IPX basis factorization --- src/ipm/ipx/basis.cc | 7 +++++++ src/ipm/ipx/basis.h | 1 + src/ipm/ipx/ipm.cc | 1 + src/ipm/ipx/kkt_solver.cc | 1 + src/ipm/ipx/kkt_solver.h | 4 ++++ src/ipm/ipx/kkt_solver_basis.h | 2 ++ src/ipm/ipx/kkt_solver_diag.h | 3 ++- 7 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/ipm/ipx/basis.cc b/src/ipm/ipx/basis.cc index 607de407d1..994335163a 100644 --- a/src/ipm/ipx/basis.cc +++ b/src/ipm/ipx/basis.cc @@ -449,6 +449,13 @@ double Basis::time_update() const { return time_update_; } +double Basis::current_fill() const { + if (fill_factors_.empty()) + return 0.0; + Int num_factors = fill_factors_.size(); + return fill_factors_[num_factors-1]; +} + double Basis::mean_fill() const { if (fill_factors_.empty()) return 0.0; diff --git a/src/ipm/ipx/basis.h b/src/ipm/ipx/basis.h index 4ae358c171..19c8a6adca 100644 --- a/src/ipm/ipx/basis.h +++ b/src/ipm/ipx/basis.h @@ -219,6 +219,7 @@ class Basis { double time_ftran() const; // time FTRAN, including partial double time_btran() const; // time BTRAN, including partial double time_update() const; // time LU update + double current_fill() const; // Current LU fill factors double mean_fill() const; // geom. mean of LU fill factors double max_fill() const; // max LU fill factor diff --git a/src/ipm/ipx/ipm.cc b/src/ipm/ipx/ipm.cc index d8860f19fc..d417534e15 100644 --- a/src/ipm/ipx/ipm.cc +++ b/src/ipm/ipx/ipm.cc @@ -89,6 +89,7 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) { kkt->Factorize(iterate, info); if (info->errflag) break; + printf("IPM::Driver fill factor = %g\n", kkt_->current_fill()); Predictor(step); if (info->errflag) break; diff --git a/src/ipm/ipx/kkt_solver.cc b/src/ipm/ipx/kkt_solver.cc index a828fcb0a3..c0d4da2397 100644 --- a/src/ipm/ipx/kkt_solver.cc +++ b/src/ipm/ipx/kkt_solver.cc @@ -19,6 +19,7 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol, Int KKTSolver::iterSum() const { return _iterSum(); } Int KKTSolver::iterMax() const { return _iterMax(); } Int KKTSolver::basis_changes() const { return _basis_changes(); } +double KKTSolver::current_fill() const {return _current_fill(); } const Basis* KKTSolver::basis() const { return _basis(); } } // namespace ipx diff --git a/src/ipm/ipx/kkt_solver.h b/src/ipm/ipx/kkt_solver.h index e5ab15bc22..7d3fe8f395 100644 --- a/src/ipm/ipx/kkt_solver.h +++ b/src/ipm/ipx/kkt_solver.h @@ -56,6 +56,9 @@ class KKTSolver { // call to Factorize(). Otherwise returns 0. Int basis_changes() const; + // + double current_fill() const; + // If a basis matrix is maintained, returns a pointer to it. // Otherwise returns NULL. const Basis* basis() const; @@ -67,6 +70,7 @@ class KKTSolver { virtual Int _iterSum() const = 0; virtual Int _iterMax() const = 0; virtual Int _basis_changes() const { return 0; } + virtual double _current_fill() const { return 0.0; } virtual const Basis* _basis() const { return nullptr; } }; diff --git a/src/ipm/ipx/kkt_solver_basis.h b/src/ipm/ipx/kkt_solver_basis.h index b46a89e5c2..b155211926 100644 --- a/src/ipm/ipx/kkt_solver_basis.h +++ b/src/ipm/ipx/kkt_solver_basis.h @@ -35,7 +35,9 @@ class KKTSolverBasis : public KKTSolver { Vector& x, Vector& y, Info* info) override; Int _iterSum() const override { return iter_sum_; } Int _iterMax() const override { return iter_max_; } + Int _basis_changes() const override { return basis_changes_; } + double _current_fill() const override { return _basis()->current_fill(); } const Basis* _basis() const override { return &basis_; } // Processes basic variables that are close to a bound by either pivoting diff --git a/src/ipm/ipx/kkt_solver_diag.h b/src/ipm/ipx/kkt_solver_diag.h index aaa18cfb85..6de0e94523 100644 --- a/src/ipm/ipx/kkt_solver_diag.h +++ b/src/ipm/ipx/kkt_solver_diag.h @@ -24,7 +24,8 @@ class KKTSolverDiag : public KKTSolver { Int maxiter() const { return maxiter_; } void maxiter(Int new_maxiter) { maxiter_ = new_maxiter; } - + // double current_fill() const { return 0; } + private: void _Factorize(Iterate* iterate, Info* info) override; void _Solve(const Vector& a, const Vector& b, double tol, From b6f29b84ff06049768f0c2c44c6640963e6d6278 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Thu, 28 Nov 2024 11:32:18 +0000 Subject: [PATCH 14/17] Reporting IPX fill factor, and added timeout to AC solve --- src/ipm/ipx/ipm.cc | 9 ++++++++- src/mip/HighsMipSolverData.cpp | 20 ++++++++++++-------- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/ipm/ipx/ipm.cc b/src/ipm/ipx/ipm.cc index d417534e15..b67aa27193 100644 --- a/src/ipm/ipx/ipm.cc +++ b/src/ipm/ipx/ipm.cc @@ -89,7 +89,6 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) { kkt->Factorize(iterate, info); if (info->errflag) break; - printf("IPM::Driver fill factor = %g\n", kkt_->current_fill()); Predictor(step); if (info->errflag) break; @@ -97,6 +96,14 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) { if (info->errflag) break; MakeStep(step); + // + std::stringstream h_logging_stream; + h_logging_stream << "IPM::Driver fill factor = " << + kkt_->current_fill() << "; iter (sum = " << + kkt_->iterSum() << ", max = " << + kkt_->iterMax() << ")\n"; + control_.hLog(h_logging_stream); + // info->iter++; PrintOutput(); } diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 0e857d4315..539495e8be 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -320,7 +320,12 @@ void HighsMipSolverData::startAnalyticCenterComputation( ipm.setOptionValue("output_flag", false); ipm.setOptionValue("output_flag", !mipsolver.submip); // 2049 unset this ultimately ipm.setOptionValue("ipm_iteration_limit", 200); - ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately + double time_available = + std::max(mipsolver.options_mip_->time_limit - + mipsolver.timer_.read(mipsolver.timer_.total_clock), + 0.1); + ipm.setOptionValue("time_limit", time_available); + // ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately // // cr1_iteration_limit is what's set internal to IPX to limit the // CR iterations before the initial basis is computed, and perhaps @@ -329,20 +334,19 @@ void HighsMipSolverData::startAnalyticCenterComputation( cr1_iteration_limit = std::min(HighsInt(500), cr1_iteration_limit); // IPX internal default cr1_iteration_limit = std::min(HighsInt(100), cr1_iteration_limit); //2049 set this ultimately // 2049 set this ultimately - ipm.setOptionValue("cr1_iteration_limit", cr1_iteration_limit); + // ipm.setOptionValue("cr1_iteration_limit", cr1_iteration_limit); // cr2_iteration_limit is not set internal to IPX, so the default // value is m+100, which is OK if you're desperate to solve an LP, // but can make the use of the AC in MIP prohibitively expensive HighsInt cr2_iteration_limit = 50 + mipsolver.model_->num_row_ / 50; cr2_iteration_limit = std::min(HighsInt(500), cr2_iteration_limit); // 2049 Set this ultimately - ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit); + // ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); - // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + - // "_ipm.mps"); + if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); if (ac_logging) num_analytic_centre_start++; double tt0 = mipsolver.timer_.read(mipsolver.timer_.total_clock); @@ -1161,9 +1165,9 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution( } this->total_repair_lp++; double time_available = - std::max(mipsolver.options_mip_->time_limit - - mipsolver.timer_.read(mipsolver.timer_.total_clock), - 0.1); + std::max(mipsolver.options_mip_->time_limit - + mipsolver.timer_.read(mipsolver.timer_.total_clock), + 0.1); Highs tmpSolver; const bool debug_report = false; if (debug_report) { From d05a2d64f711b964634cb252eff7572335ab1964 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Thu, 28 Nov 2024 20:29:39 +0000 Subject: [PATCH 15/17] Added lp_data/HighsSolutionStats.h --- cmake/sources-python.cmake | 1 + cmake/sources.cmake | 1 + src/lp_data/HStruct.h | 14 -------------- src/simplex/HEkk.h | 1 + 4 files changed, 3 insertions(+), 14 deletions(-) diff --git a/cmake/sources-python.cmake b/cmake/sources-python.cmake index df2ff5a223..b6d3349a3e 100644 --- a/cmake/sources-python.cmake +++ b/cmake/sources-python.cmake @@ -311,6 +311,7 @@ set(highs_headers_python src/lp_data/HighsRanging.h src/lp_data/HighsSolution.h src/lp_data/HighsSolutionDebug.h + src/lp_data/HighsSolutionStats.h src/lp_data/HighsSolve.h src/lp_data/HighsStatus.h src/lp_data/HStruct.h diff --git a/cmake/sources.cmake b/cmake/sources.cmake index 96e84b7f5b..7335a2e83d 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -315,6 +315,7 @@ set(highs_headers lp_data/HighsRanging.h lp_data/HighsSolution.h lp_data/HighsSolutionDebug.h + lp_data/HighsSolutionStats.h lp_data/HighsSolve.h lp_data/HighsStatus.h lp_data/HStruct.h diff --git a/src/lp_data/HStruct.h b/src/lp_data/HStruct.h index 35a0b0ee42..3963d9df1c 100644 --- a/src/lp_data/HStruct.h +++ b/src/lp_data/HStruct.h @@ -154,18 +154,4 @@ struct HighsLinearObjective { void clear(); }; -struct HighsSimplexStats { - bool valid; - HighsInt iteration_count; - HighsInt num_invert; - HighsInt last_invert_num_el; - HighsInt last_factored_basis_num_el; - double col_aq_density; - double row_ep_density; - double row_ap_density; - double row_DSE_density; - void report(FILE* file, const std::string message = "") const; - void initialise(const HighsInt iteration_count_ = 0); -}; - #endif /* LP_DATA_HSTRUCT_H_ */ diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h index 0b253bebc0..b3b163a769 100644 --- a/src/simplex/HEkk.h +++ b/src/simplex/HEkk.h @@ -17,6 +17,7 @@ #include "lp_data/HighsCallback.h" #include "simplex/HSimplexNla.h" #include "simplex/HighsSimplexAnalysis.h" +#include "lp_data/HighsSolutionStats.h" #include "util/HSet.h" #include "util/HighsHash.h" #include "util/HighsRandom.h" From 4f393aaa8369a8e525156a8fcd7ad12c287d8c1f Mon Sep 17 00:00:00 2001 From: JAJHall Date: Fri, 29 Nov 2024 09:40:49 +0000 Subject: [PATCH 16/17] Now including simplex stats in IPX --- src/Highs.h | 3 +++ src/ipm/IpxWrapper.cpp | 16 ++++++++++----- src/ipm/IpxWrapper.h | 6 ++++-- src/ipm/ipx/ipx_parameters.h | 5 ++++- src/lp_data/HighsSolutionStats.h | 35 ++++++++++++++++++++++++++++++++ src/mip/HighsLpRelaxation.cpp | 8 ++++++-- src/mip/HighsMipSolverData.cpp | 4 ++++ src/simplex/HEkk.cpp | 4 ++++ src/simplex/HEkk.h | 1 + 9 files changed, 72 insertions(+), 10 deletions(-) create mode 100644 src/lp_data/HighsSolutionStats.h diff --git a/src/Highs.h b/src/Highs.h index 70d800a5de..286eb002ed 100644 --- a/src/Highs.h +++ b/src/Highs.h @@ -1219,6 +1219,9 @@ class Highs { const HighsSimplexStats& getSimplexStats() const { return ekk_instance_.getSimplexStats(); } + void passSimplexStats(const HighsSimplexStats simplex_stats) { + return ekk_instance_.passSimplexStats(simplex_stats); + } void reportSimplexStats(FILE* file) const { ekk_instance_.reportSimplexStats(file); } diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index 4a4e3d4c76..0ab4e58ebc 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -23,15 +23,19 @@ using std::min; HighsStatus solveLpIpx(HighsLpSolverObject& solver_object) { return solveLpIpx(solver_object.options_, solver_object.timer_, - solver_object.lp_, solver_object.basis_, - solver_object.solution_, solver_object.model_status_, - solver_object.highs_info_, solver_object.callback_); + solver_object.lp_, solver_object.ekk_instance_, + solver_object.basis_, solver_object.solution_, + solver_object.model_status_, + solver_object.highs_info_, + solver_object.callback_); } HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, - const HighsLp& lp, HighsBasis& highs_basis, + const HighsLp& lp, const HEkk& ekk_instance, + HighsBasis& highs_basis, HighsSolution& highs_solution, - HighsModelStatus& model_status, HighsInfo& highs_info, + HighsModelStatus& model_status, + HighsInfo& highs_info, HighsCallback& callback) { // Use IPX to try to solve the LP // @@ -149,6 +153,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, parameters.max_centring_steps = options.max_centring_steps; parameters.centring_ratio_tolerance = options.centring_ratio_tolerance; + parameters.simplex_stats = ekk_instance.getSimplexStats(); + // Set the internal IPX parameters lps.SetParameters(parameters); diff --git a/src/ipm/IpxWrapper.h b/src/ipm/IpxWrapper.h index 660165d82d..6f8a75d8b3 100644 --- a/src/ipm/IpxWrapper.h +++ b/src/ipm/IpxWrapper.h @@ -25,9 +25,11 @@ HighsStatus solveLpIpx(HighsLpSolverObject& solver_object); HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, - const HighsLp& lp, HighsBasis& highs_basis, + const HighsLp& lp, const HEkk& ekk_instance, + HighsBasis& highs_basis, HighsSolution& highs_solution, - HighsModelStatus& model_status, HighsInfo& highs_info, + HighsModelStatus& model_status, + HighsInfo& highs_info, HighsCallback& callback); void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row, diff --git a/src/ipm/ipx/ipx_parameters.h b/src/ipm/ipx/ipx_parameters.h index 029b89f65c..1fc4af238a 100644 --- a/src/ipm/ipx/ipx_parameters.h +++ b/src/ipm/ipx/ipx_parameters.h @@ -2,6 +2,7 @@ #define IPX_PARAMETERS_H_ #include "io/HighsIO.h" +#include "lp_data/HighsSolutionStats.h" #include "ipm/ipx/ipx_config.h" #include @@ -68,7 +69,9 @@ struct ipx_parameters { /* HiGHS logging parameters */ bool highs_logging; const HighsLogOptions* log_options; - + + /* Simplex solution stats */ + HighsSimplexStats simplex_stats; }; #ifdef __cplusplus diff --git a/src/lp_data/HighsSolutionStats.h b/src/lp_data/HighsSolutionStats.h new file mode 100644 index 0000000000..ea43a9522a --- /dev/null +++ b/src/lp_data/HighsSolutionStats.h @@ -0,0 +1,35 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* */ +/* This file is part of the HiGHS linear optimization suite */ +/* */ +/* Written and engineered 2008-2024 by Julian Hall, Ivet Galabova, */ +/* Leona Gottwald and Michael Feldmeier */ +/* */ +/* Available as open-source under the MIT License */ +/* */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/**@file lp_data/HighsSolutionStats.h + * @brief Structs for HiGHS + */ +#ifndef LP_DATA_HIGHSSOLUTIONSTATS_H_ +#define LP_DATA_HIGHSSOLUTIONSTATS_H_ + +#include + +#include "lp_data/HConst.h" + +struct HighsSimplexStats { + bool valid; + HighsInt iteration_count; + HighsInt num_invert; + HighsInt last_invert_num_el; + HighsInt last_factored_basis_num_el; + double col_aq_density; + double row_ep_density; + double row_ap_density; + double row_DSE_density; + void report(FILE* file, const std::string message = "") const; + void initialise(const HighsInt iteration_count_ = 0); +}; + +#endif /* LP_DATA_HIGHSSOLUTIONSTATS_H_ */ diff --git a/src/mip/HighsLpRelaxation.cpp b/src/mip/HighsLpRelaxation.cpp index 4288f99437..ff2fc7d3b1 100644 --- a/src/mip/HighsLpRelaxation.cpp +++ b/src/mip/HighsLpRelaxation.cpp @@ -1064,9 +1064,11 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) { int(lpsolver.getNumCol()), int(lpsolver.getNumRow())); } } - const bool solver_logging = false; + const bool solver_logging = false;//!mipsolver.submip && !valid_basis; const bool detailed_simplex_logging = false; - if (solver_logging) lpsolver.setOptionValue("output_flag", true); + if (solver_logging) { + lpsolver.setOptionValue("output_flag", true); + } if (detailed_simplex_logging) { lpsolver.setOptionValue("output_flag", true); lpsolver.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose); @@ -1076,6 +1078,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) { mipsolver.analysis_.mipTimerStart(simplex_solve_clock); HighsStatus callstatus = lpsolver.run(); + lpsolver.setOptionValue("output_flag", false); // !fix-2049 + // if (solver_logging) lpsolver.reportSimplexStats(stdout); // !fix-2049 mipsolver.analysis_.mipTimerStop(simplex_solve_clock); const HighsInfo& info = lpsolver.getInfo(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 1f1b2830cf..8a07cfb5c4 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -342,6 +342,10 @@ void HighsMipSolverData::startAnalyticCenterComputation( cr2_iteration_limit = std::min(HighsInt(500), cr2_iteration_limit); // 2049 Set this ultimately // ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit); + + HighsSimplexStats simplex_stats = lp.getLpSolver().getSimplexStats(); + if (!mipsolver.submip) simplex_stats.report(stdout); + ipm.passSimplexStats(simplex_stats); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp index caa4e641e0..ca743c1f98 100644 --- a/src/simplex/HEkk.cpp +++ b/src/simplex/HEkk.cpp @@ -4421,6 +4421,10 @@ void HEkk::unitBtranResidual(const HighsInt row_out, const HVector& row_ep, } } +void HEkk::passSimplexStats(const HighsSimplexStats simplex_stats) { + this->simplex_stats_ = simplex_stats; +} + void HighsSimplexStats::report(FILE* file, std::string message) const { fprintf(file, "\nSimplex stats: %s\n", message.c_str()); fprintf(file, " valid = %d\n", this->valid); diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h index b3b163a769..d62a4e9fe9 100644 --- a/src/simplex/HEkk.h +++ b/src/simplex/HEkk.h @@ -158,6 +158,7 @@ class HEkk { const vector& rowUpper); const HighsSimplexStats& getSimplexStats() const { return simplex_stats_; } + void passSimplexStats(const HighsSimplexStats simplex_stats); void initialiseSimplexStats() { simplex_stats_.initialise(iteration_count_); } void reportSimplexStats(FILE* file, const std::string message = "") const { simplex_stats_.report(file, message); From 54ce3ae06ca8e0f7cc036eaa92a20154e4c4c3be Mon Sep 17 00:00:00 2001 From: JAJHall Date: Fri, 29 Nov 2024 09:41:11 +0000 Subject: [PATCH 17/17] Formatted --- src/ipm/IpxWrapper.cpp | 13 +++----- src/ipm/IpxWrapper.h | 6 ++-- src/lp_data/HighsOptions.h | 3 +- src/mip/HighsLpRelaxation.cpp | 4 +-- src/mip/HighsMipSolverData.cpp | 52 +++++++++++++++++-------------- src/mip/HighsPrimalHeuristics.cpp | 3 +- src/simplex/HEkk.h | 2 +- 7 files changed, 42 insertions(+), 41 deletions(-) diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index 0ab4e58ebc..7a073eeba5 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -24,18 +24,15 @@ using std::min; HighsStatus solveLpIpx(HighsLpSolverObject& solver_object) { return solveLpIpx(solver_object.options_, solver_object.timer_, solver_object.lp_, solver_object.ekk_instance_, - solver_object.basis_, solver_object.solution_, - solver_object.model_status_, - solver_object.highs_info_, - solver_object.callback_); + solver_object.basis_, solver_object.solution_, + solver_object.model_status_, solver_object.highs_info_, + solver_object.callback_); } HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, const HighsLp& lp, const HEkk& ekk_instance, - HighsBasis& highs_basis, - HighsSolution& highs_solution, - HighsModelStatus& model_status, - HighsInfo& highs_info, + HighsBasis& highs_basis, HighsSolution& highs_solution, + HighsModelStatus& model_status, HighsInfo& highs_info, HighsCallback& callback) { // Use IPX to try to solve the LP // diff --git a/src/ipm/IpxWrapper.h b/src/ipm/IpxWrapper.h index 6f8a75d8b3..792cf05346 100644 --- a/src/ipm/IpxWrapper.h +++ b/src/ipm/IpxWrapper.h @@ -26,10 +26,8 @@ HighsStatus solveLpIpx(HighsLpSolverObject& solver_object); HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, const HighsLp& lp, const HEkk& ekk_instance, - HighsBasis& highs_basis, - HighsSolution& highs_solution, - HighsModelStatus& model_status, - HighsInfo& highs_info, + HighsBasis& highs_basis, HighsSolution& highs_solution, + HighsModelStatus& model_status, HighsInfo& highs_info, HighsCallback& callback); void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row, diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 359f179da1..6ae321368d 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -1088,7 +1088,8 @@ class HighsOptions : public HighsOptionsStruct { records.push_back(record_int); record_bool = new OptionRecordBool( - "kkt_logging", "Perform logging in CR solver for KKT in IPX: Default = false", + "kkt_logging", + "Perform logging in CR solver for KKT in IPX: Default = false", advanced, &kkt_logging, false); records.push_back(record_bool); diff --git a/src/mip/HighsLpRelaxation.cpp b/src/mip/HighsLpRelaxation.cpp index ff2fc7d3b1..73b675917c 100644 --- a/src/mip/HighsLpRelaxation.cpp +++ b/src/mip/HighsLpRelaxation.cpp @@ -1064,7 +1064,7 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) { int(lpsolver.getNumCol()), int(lpsolver.getNumRow())); } } - const bool solver_logging = false;//!mipsolver.submip && !valid_basis; + const bool solver_logging = false; //! mipsolver.submip && !valid_basis; const bool detailed_simplex_logging = false; if (solver_logging) { lpsolver.setOptionValue("output_flag", true); @@ -1078,7 +1078,7 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) { mipsolver.analysis_.mipTimerStart(simplex_solve_clock); HighsStatus callstatus = lpsolver.run(); - lpsolver.setOptionValue("output_flag", false); // !fix-2049 + lpsolver.setOptionValue("output_flag", false); // !fix-2049 // if (solver_logging) lpsolver.reportSimplexStats(stdout); // !fix-2049 mipsolver.analysis_.mipTimerStop(simplex_solve_clock); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 8a07cfb5c4..f03114f2ab 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -318,22 +318,26 @@ void HighsMipSolverData::startAnalyticCenterComputation( } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - ipm.setOptionValue("output_flag", !mipsolver.submip); // 2049 unset this ultimately + ipm.setOptionValue("output_flag", + !mipsolver.submip); // 2049 unset this ultimately ipm.setOptionValue("ipm_iteration_limit", 200); double time_available = - std::max(mipsolver.options_mip_->time_limit - - mipsolver.timer_.read(mipsolver.timer_.total_clock), - 0.1); + std::max(mipsolver.options_mip_->time_limit - + mipsolver.timer_.read(mipsolver.timer_.total_clock), + 0.1); ipm.setOptionValue("time_limit", time_available); - // ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this ultimately + // ipm.setOptionValue("kkt_logging", !mipsolver.submip); // 2049 unset this + // ultimately // // cr1_iteration_limit is what's set internal to IPX to limit the // CR iterations before the initial basis is computed, and perhaps // should not be changed, but maybe 500 is too high HighsInt cr1_iteration_limit = 10 + mipsolver.model_->num_row_ / 20; - cr1_iteration_limit = std::min(HighsInt(500), cr1_iteration_limit); // IPX internal default - cr1_iteration_limit = std::min(HighsInt(100), cr1_iteration_limit); //2049 set this ultimately - // 2049 set this ultimately + cr1_iteration_limit = + std::min(HighsInt(500), cr1_iteration_limit); // IPX internal default + cr1_iteration_limit = std::min( + HighsInt(100), cr1_iteration_limit); // 2049 set this ultimately + // 2049 set this ultimately // ipm.setOptionValue("cr1_iteration_limit", cr1_iteration_limit); // cr2_iteration_limit is not set internal to IPX, so the default // value is m+100, which is OK if you're desperate to solve an LP, @@ -350,7 +354,8 @@ void HighsMipSolverData::startAnalyticCenterComputation( lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); - if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); + if (!mipsolver.submip) + ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); if (ac_logging) num_analytic_centre_start++; double tt0 = mipsolver.timer_.read(mipsolver.timer_.total_clock); @@ -364,17 +369,14 @@ void HighsMipSolverData::startAnalyticCenterComputation( "grepAcLocal: model; num_row; CR limits+counts; time and status," "%s,%d," "%d,%d," - "%d,%d," + "%d,%d," "%g,%s\n", mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), - int(cr1_iteration_limit), - int(ipm.getInfo().max_cr_iteration_count1), - int(cr2_iteration_limit), - int(ipm.getInfo().max_cr_iteration_count2), - tt1 - tt0, - ipm.modelStatusToString(ac_status).c_str()); + int(cr1_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1), + int(cr2_iteration_limit), int(ipm.getInfo().max_cr_iteration_count2), + tt1 - tt0, ipm.modelStatusToString(ac_status).c_str()); if (ac_status == HighsModelStatus::kSolveError || - ac_status == HighsModelStatus::kUnknown) { + ac_status == HighsModelStatus::kUnknown) { num_analytic_centre_fail++; } else if (ac_status == HighsModelStatus::kOptimal) { num_analytic_centre_opt++; @@ -406,8 +408,8 @@ void HighsMipSolverData::finishAnalyticCenterComputation( } analyticCenterComputed = true; analyticCenterFailed = - analyticCenterStatus == HighsModelStatus::kSolveError || - analyticCenterStatus == HighsModelStatus::kUnknown; + analyticCenterStatus == HighsModelStatus::kSolveError || + analyticCenterStatus == HighsModelStatus::kUnknown; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -445,8 +447,10 @@ void HighsMipSolverData::finishAnalyticCenterComputation( mipsolver.mipdata_->domain.propagate(); if (mipsolver.mipdata_->domain.infeasible()) return; } else if (!analyticCenterFailed) { - printf("HighsMipSolverData::finishAnalyticCenterComputation: analyticCenterStatus = %s\n", - lp.getLpSolver().modelStatusToString(analyticCenterStatus).c_str()); + printf( + "HighsMipSolverData::finishAnalyticCenterComputation: " + "analyticCenterStatus = %s\n", + lp.getLpSolver().modelStatusToString(analyticCenterStatus).c_str()); } } @@ -1169,9 +1173,9 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution( } this->total_repair_lp++; double time_available = - std::max(mipsolver.options_mip_->time_limit - - mipsolver.timer_.read(mipsolver.timer_.total_clock), - 0.1); + std::max(mipsolver.options_mip_->time_limit - + mipsolver.timer_.read(mipsolver.timer_.total_clock), + 0.1); Highs tmpSolver; const bool debug_report = false; if (debug_report) { diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 1812803af0..8b147fd49d 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -137,7 +137,8 @@ bool HighsPrimalHeuristics::solveSubMip( // 2049 Set this ultimately // - if (mipsolver.mipdata_->analyticCenterFailed) submipoptions.mip_compute_analytic_centre = 0; + if (mipsolver.mipdata_->analyticCenterFailed) + submipoptions.mip_compute_analytic_centre = 0; // setup solver and run it HighsSolution solution; diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h index d62a4e9fe9..62594dd822 100644 --- a/src/simplex/HEkk.h +++ b/src/simplex/HEkk.h @@ -15,9 +15,9 @@ #define SIMPLEX_HEKK_H_ #include "lp_data/HighsCallback.h" +#include "lp_data/HighsSolutionStats.h" #include "simplex/HSimplexNla.h" #include "simplex/HighsSimplexAnalysis.h" -#include "lp_data/HighsSolutionStats.h" #include "util/HSet.h" #include "util/HighsHash.h" #include "util/HighsRandom.h"