Skip to content

Commit

Permalink
Now to test AC, gathering more data
Browse files Browse the repository at this point in the history
  • Loading branch information
@jajhall
jajhall committed Nov 26, 2024
1 parent ddd52b3 commit f3ff0b4
Show file tree
Hide file tree
Showing 5 changed files with 73 additions and 23 deletions.
31 changes: 22 additions & 9 deletions src/ipm/ipx/conjugate_residuals.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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
Expand All @@ -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]));
Expand All @@ -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) << ','
Expand Down Expand Up @@ -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();
}

Expand Down Expand Up @@ -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
Expand All @@ -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]));
Expand All @@ -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';
Expand Down Expand Up @@ -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();
}

Expand Down
18 changes: 14 additions & 4 deletions src/mip/HighsMipSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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),
Expand All @@ -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();

Expand Down
36 changes: 28 additions & 8 deletions src/mip/HighsMipSolverData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -317,23 +318,40 @@ 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));

// 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<double>& sol = ipm.getSolution().col_value;
if (HighsInt(sol.size()) != mipsolver.numCol()) return;
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
10 changes: 8 additions & 2 deletions src/mip/HighsMipSolverData.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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),
Expand Down
1 change: 1 addition & 0 deletions src/mip/HighsPrimalHeuristics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit f3ff0b4

Please sign in to comment.