diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp index f2543ba1f5..b8f7e249cf 100644 --- a/src/highs_bindings.cpp +++ b/src/highs_bindings.cpp @@ -832,7 +832,9 @@ PYBIND11_MODULE(_core, m) { .def_readwrite("max_complementarity_violation", &HighsInfo::max_complementarity_violation) .def_readwrite("sum_complementarity_violations", - &HighsInfo::sum_complementarity_violations); + &HighsInfo::sum_complementarity_violations) + .def_readwrite("primal_dual_integral", + &HighsInfo::primal_dual_integral); py::class_(m, "HighsOptions") .def(py::init<>()) .def_readwrite("presolve", &HighsOptions::presolve) diff --git a/src/interfaces/highs_c_api.cpp b/src/interfaces/highs_c_api.cpp index 5dd76f62b4..8dad0e5216 100644 --- a/src/interfaces/highs_c_api.cpp +++ b/src/interfaces/highs_c_api.cpp @@ -529,6 +529,17 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray, return retcode; } +HighsInt Highs_getDualUnboundednessDirection( + const void* highs, HighsInt* has_dual_unboundedness_direction, + double* dual_unboundedness_direction_value) { + bool v; + HighsInt retcode = (HighsInt)((Highs*)highs) + ->getDualUnboundednessDirection( + v, dual_unboundedness_direction_value); + *has_dual_unboundedness_direction = (HighsInt)v; + return retcode; +} + HighsInt Highs_getPrimalRay(const void* highs, HighsInt* has_primal_ray, double* primal_ray_value) { bool v; diff --git a/src/interfaces/highs_c_api.h b/src/interfaces/highs_c_api.h index 6237dbac49..a59ec64a17 100644 --- a/src/interfaces/highs_c_api.h +++ b/src/interfaces/highs_c_api.h @@ -915,11 +915,13 @@ HighsInt Highs_getBasis(const void* highs, HighsInt* col_status, HighsInt Highs_getModelStatus(const void* highs); /** - * Get an unbounded dual ray that is a certificate of primal infeasibility. + * Indicates whether a dual ray that is a certificate of primal + * infeasibility currently exists, and (at the expense of solving an + * LP) gets it if it does not and dual_ray_value is not nullptr. * * @param highs A pointer to the Highs instance. - * @param has_dual_ray A pointer to an int to store 1 if the dual ray - * exists. + * @param has_dual_ray A pointer to an int to store 1 if a dual ray + * currently exists. * @param dual_ray_value An array of length [num_row] filled with the * unbounded ray. * @@ -929,7 +931,28 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray, double* dual_ray_value); /** - * Get an unbounded primal ray that is a certificate of dual infeasibility. + * Indicates whether a dual unboundedness direction (corresponding to a + * certificate of primal infeasibility) exists, and (at the expense of + * solving an LP) gets it if it does not and + * dual_unboundedness_direction is not nullptr + * + * @param highs A pointer to the Highs + * instance. + * @param has_dual_unboundedness_direction A pointer to an int to store 1 + * if the dual unboundedness + * direction exists. + * @param dual_unboundedness_direction_value An array of length [num_col] + * filled with the unboundedness + * direction. + */ +HighsInt getDualUnboundednessDirection( + const void* highs, HighsInt* has_dual_unboundedness_direction, + double* dual_unboundedness_direction_value); + +/** + * Indicates whether a primal ray that is a certificate of primal + * unboundedness currently exists, and (at the expense of solving an + * LP) gets it if it does not and primal_ray_value is not nullptr. * * @param highs A pointer to the Highs instance. * @param has_primal_ray A pointer to an int to store 1 if the primal ray diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index 5d75b493fe..823bff0047 100644 --- a/src/lp_data/Highs.cpp +++ b/src/lp_data/Highs.cpp @@ -3938,6 +3938,7 @@ HighsStatus Highs::callSolveMip() { info_.mip_node_count = solver.node_count_; info_.mip_dual_bound = solver.dual_bound_; info_.mip_gap = solver.gap_; + info_.primal_dual_integral = solver.primal_dual_integral_; // Get the number of LP iterations, avoiding overflow if the int64_t // value is too large int64_t mip_total_lp_iterations = solver.total_lp_iterations_; @@ -4502,6 +4503,7 @@ HighsStatus Highs::returnFromHighs(HighsStatus highs_return_status) { } void Highs::reportSolvedLpQpStats() { + if (!options_.output_flag) return; HighsLogOptions& log_options = options_.log_options; if (this->model_.lp_.model_name_.length()) highsLogUser(log_options, HighsLogType::kInfo, "Model name : %s\n", @@ -4541,11 +4543,7 @@ void Highs::reportSolvedLpQpStats() { std::fabs(info_.objective_function_value - dual_objective_value) / std::max(1.0, std::fabs(info_.objective_function_value)); highsLogUser(log_options, HighsLogType::kInfo, - "Highs::reportSolvedLpQpStats Objective for %s: primal = " - "%17.10e; dual = %17.10e; rel gap = %17.10e\n", - this->model_.lp_.model_name_.c_str(), - info_.objective_function_value, dual_objective_value, - relative_primal_dual_gap); + "Relative P-D gap : %17.10e\n", relative_primal_dual_gap); } double run_time = timer_.readRunHighsClock(); highsLogUser(log_options, HighsLogType::kInfo, diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp index 99fb4ad16c..ce01a585b8 100644 --- a/src/lp_data/HighsInfo.cpp +++ b/src/lp_data/HighsInfo.cpp @@ -40,6 +40,7 @@ void HighsInfo::invalidate() { sum_dual_infeasibilities = kHighsIllegalInfeasibilityMeasure; max_complementarity_violation = kHighsIllegalComplementarityViolation; sum_complementarity_violations = kHighsIllegalComplementarityViolation; + primal_dual_integral = -kHighsInf; } static std::string infoEntryTypeToString(const HighsInfoType type) { diff --git a/src/lp_data/HighsInfo.h b/src/lp_data/HighsInfo.h index 94766e7a5f..3292b8281e 100644 --- a/src/lp_data/HighsInfo.h +++ b/src/lp_data/HighsInfo.h @@ -115,6 +115,7 @@ struct HighsInfoStruct { double sum_dual_infeasibilities; double max_complementarity_violation; double sum_complementarity_violations; + double primal_dual_integral; }; class HighsInfo : public HighsInfoStruct { @@ -273,6 +274,11 @@ class HighsInfo : public HighsInfoStruct { "sum_complementarity_violations", "Sum of complementarity violations", advanced, &sum_complementarity_violations, 0); records.push_back(record_double); + + record_double = + new InfoRecordDouble("primal_dual_integral", "Primal-dual integral", + advanced, &primal_dual_integral, 0); + records.push_back(record_double); } public: diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index 8482442e35..2b5b93e6e6 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -225,8 +225,16 @@ void HighsMipSolver::run() { search.setLpRelaxation(&mipdata_->lp); sepa.setLpRelaxation(&mipdata_->lp); + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = mipdata_->nodequeue.getBestLowerBound(); + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral(prev_lower_bound, mipdata_->lower_bound, + mipdata_->upper_bound, + mipdata_->upper_bound); + mipdata_->printDisplayLine(); search.installNode(mipdata_->nodequeue.popBestBoundNode()); int64_t numStallNodes = 0; @@ -348,8 +356,16 @@ void HighsMipSolver::run() { search.flushStatistics(); if (limit_reached) { + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min(mipdata_->upper_bound, mipdata_->nodequeue.getBestLowerBound()); + + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); mipdata_->printDisplayLine(); break; } @@ -371,13 +387,29 @@ void HighsMipSolver::run() { if (mipdata_->domain.infeasible()) { mipdata_->nodequeue.clear(); mipdata_->pruned_treeweight = 1.0; + + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound); + + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); mipdata_->printDisplayLine(); break; } + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min(mipdata_->upper_bound, mipdata_->nodequeue.getBestLowerBound()); + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); mipdata_->printDisplayLine(); if (mipdata_->nodequeue.empty()) break; @@ -538,7 +570,16 @@ void HighsMipSolver::run() { if (mipdata_->domain.infeasible()) { mipdata_->nodequeue.clear(); mipdata_->pruned_treeweight = 1.0; + + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound); + + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); break; } @@ -547,9 +588,16 @@ void HighsMipSolver::run() { break; } + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min( mipdata_->upper_bound, mipdata_->nodequeue.getBestLowerBound()); + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); mipdata_->printDisplayLine(); if (!mipdata_->domain.getChangedCols().empty()) { @@ -579,7 +627,16 @@ void HighsMipSolver::run() { search.openNodesToQueue(mipdata_->nodequeue); mipdata_->nodequeue.clear(); mipdata_->pruned_treeweight = 1.0; + + double prev_lower_bound = mipdata_->lower_bound; + mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound); + + bool bound_change = mipdata_->lower_bound != prev_lower_bound; + if (!submip && bound_change) + mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound); break; } @@ -609,12 +666,20 @@ void HighsMipSolver::run() { } void HighsMipSolver::cleanupSolve() { + // Force a final logging line mipdata_->printDisplayLine(kSolutionSourceCleanup); // Stop the solve clock - which won't be running if presolve // determines the model status if (analysis_.mipTimerRunning(kMipClockSolve)) analysis_.mipTimerStop(kMipClockSolve); + + // Need to complete the calculation of P-D integral, checking for NO + // gap change + mipdata_->updatePrimalDualIntegral( + mipdata_->lower_bound, mipdata_->lower_bound, mipdata_->upper_bound, + mipdata_->upper_bound, false); analysis_.mipTimerStart(kMipClockPostsolve); + bool havesolution = solution_objective_ != kHighsInf; bool feasible; if (havesolution) @@ -639,6 +704,7 @@ void HighsMipSolver::cleanupSolve() { node_count_ = mipdata_->num_nodes; total_lp_iterations_ = mipdata_->total_lp_iterations; dual_bound_ = std::min(dual_bound_, primal_bound_); + primal_dual_integral_ = mipdata_->primal_dual_integral.value; // adjust objective sense in case of maximization problem if (orig_model_->sense_ == ObjSense::kMaximize) { @@ -716,9 +782,11 @@ void HighsMipSolver::cleanupSolve() { " Primal bound %.12g\n" " Dual bound %.12g\n" " Gap %s\n" + " P-D integral %.12g\n" " Solution status %s\n", utilModelStatusToString(modelstatus_).c_str(), primal_bound_, - dual_bound_, gapString.data(), solutionstatus.c_str()); + dual_bound_, gapString.data(), + mipdata_->primal_dual_integral.value, solutionstatus.c_str()); if (solutionstatus != "-") highsLogUser(options_mip_->log_options, HighsLogType::kInfo, " %.12g (objective)\n" diff --git a/src/mip/HighsMipSolver.h b/src/mip/HighsMipSolver.h index d50d7c68a8..0a1eaae0ca 100644 --- a/src/mip/HighsMipSolver.h +++ b/src/mip/HighsMipSolver.h @@ -41,6 +41,7 @@ class HighsMipSolver { double gap_; int64_t node_count_; int64_t total_lp_iterations_; + double primal_dual_integral_; FILE* improving_solution_file_; std::vector saved_objective_and_solution_; diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 7762d7495c..6516dc3e89 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -448,6 +448,26 @@ void HighsMipSolverData::finishSymmetryDetection( globalOrbits = symmetries.computeStabilizerOrbits(domain); } +double HighsMipSolverData::limitsToGap(const double use_lower_bound, + const double use_upper_bound, double& lb, + double& ub) const { + double offset = mipsolver.model_->offset_; + lb = use_lower_bound + offset; + if (std::abs(lb) <= epsilon) lb = 0; + ub = kHighsInf; + double gap = kHighsInf; + if (use_upper_bound != kHighsInf) { + ub = use_upper_bound + offset; + if (std::fabs(ub) <= epsilon) ub = 0; + lb = std::min(ub, lb); + if (ub == 0.0) + gap = lb == 0.0 ? 0.0 : kHighsInf; + else + gap = (ub - lb) / fabs(ub); + } + return gap; +} + double HighsMipSolverData::computeNewUpperLimit(double ub, double mip_abs_gap, double mip_rel_gap) const { double new_upper_limit; @@ -622,6 +642,7 @@ void HighsMipSolverData::init() { upper_bound = kHighsInf; upper_limit = mipsolver.options_mip_->objective_bound; optimality_limit = mipsolver.options_mip_->objective_bound; + primal_dual_integral.initialise(); if (mipsolver.options_mip_->mip_report_level == 0) dispfreq = 0; @@ -661,9 +682,14 @@ void HighsMipSolverData::runSetup() { last_disptime = -kHighsInf; - // transform the objective limit to the current model + // Transform the reference of the objective limit and lower/upper + // bounds from the original model to the current model, undoing the + // transformation done before restart so that the offset change due + // to presolve is incorporated. Bound changes are transitory, so no + // real gap change, and no update to P-D integral is necessary upper_limit -= mipsolver.model_->offset_; optimality_limit -= mipsolver.model_->offset_; + lower_bound -= mipsolver.model_->offset_; upper_bound -= mipsolver.model_->offset_; @@ -686,7 +712,15 @@ void HighsMipSolverData::runSetup() { mipsolver.solution_objective_); } if (feasible && solobj < upper_bound) { + double prev_upper_bound = upper_bound; + upper_bound = solobj; + + bool bound_change = upper_bound != prev_upper_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(lower_bound, lower_bound, prev_upper_bound, + upper_bound); + double new_upper_limit = computeNewUpperLimit(solobj, 0.0, 0.0); saveReportMipSolution(new_upper_limit); if (new_upper_limit < upper_limit) { @@ -792,7 +826,16 @@ void HighsMipSolverData::runSetup() { domain.propagate(); if (domain.infeasible()) { mipsolver.modelstatus_ = HighsModelStatus::kInfeasible; + + double prev_lower_bound = lower_bound; + lower_bound = kHighsInf; + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); + pruned_treeweight = 1.0; return; } @@ -836,7 +879,16 @@ void HighsMipSolverData::runSetup() { if (fractionality(domain.col_lower_[i]) > feastol) { // integer variable is fixed to a fractional value -> infeasible mipsolver.modelstatus_ = HighsModelStatus::kInfeasible; + + double prev_lower_bound = lower_bound; + lower_bound = kHighsInf; + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, + upper_bound, upper_bound); + pruned_treeweight = 1.0; return; } @@ -1137,7 +1189,6 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution( check_row_data += ")"; } highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kWarning, - // printf( "Solution with objective %g has untransformed violations: " "bound = %.4g%s; integrality = %.4g%s; row = %.4g%s\n", mipsolver_objective_value, bound_violation_, @@ -1232,10 +1283,13 @@ void HighsMipSolverData::performRestart() { mipsolver.rootbasis = &root_basis; } - // transform the objective upper bound into the original space, as it is - // expected during presolve + // Transform the reference of the objective limit and lower/upper + // bounds to the original model, since offset will generally changte + // in presolve. Bound changes are transitory, so no real gap change, + // and no update to P-D integral is necessary upper_limit += mipsolver.model_->offset_; optimality_limit += mipsolver.model_->offset_; + upper_bound += mipsolver.model_->offset_; lower_bound += mipsolver.model_->offset_; @@ -1280,15 +1334,34 @@ void HighsMipSolverData::performRestart() { mipsolver.mipdata_->upper_bound = 0; mipsolver.mipdata_->transformNewIntegerFeasibleSolution( std::vector()); - } else + } else { upper_bound -= mipsolver.model_->offset_; + } + + // lower_bound still relates to the original model, and the offset + // is never applied, since MIP solving is complete, and + // lower_bound is set to upper_bound, so apply the offset now, so + // that housekeeping in updatePrimalDualIntegral is correct + double prev_lower_bound = lower_bound - mipsolver.model_->offset_; lower_bound = upper_bound; + + // There must be a gap change, since it's now zero, so always call + // updatePrimalDualIntegral (unless solving a sub-MIP) + // + // Surely there must be a lower bound change + bool bound_change = lower_bound != prev_lower_bound; + assert(bound_change); + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); if (mipsolver.solution_objective_ != kHighsInf && mipsolver.modelstatus_ == HighsModelStatus::kInfeasible) mipsolver.modelstatus_ = HighsModelStatus::kOptimal; return; } + // Bounds are currently in the original space since presolve will have + // changed offset_ runSetup(); postSolveStack.removeCutsFromModel(numCuts); @@ -1358,7 +1431,16 @@ bool HighsMipSolverData::addIncumbent(const std::vector& sol, // solobj = transformNewIntegerFeasibleSolution(sol); if (solobj >= upper_bound) return false; + + double prev_upper_bound = upper_bound; + upper_bound = solobj; + + bool bound_change = upper_bound != prev_upper_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(lower_bound, lower_bound, prev_upper_bound, + upper_bound); + incumbent = sol; double new_upper_limit = computeNewUpperLimit(solobj, 0.0, 0.0); @@ -1549,23 +1631,15 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) { double explored = 100 * double(pruned_treeweight); - double offset = mipsolver.model_->offset_; - double lb = lower_bound + offset; - if (std::abs(lb) <= epsilon) lb = 0; - double ub = kHighsInf; - double gap = kHighsInf; + double lb; + double ub; + double gap = limitsToGap(lower_bound, upper_bound, lb, ub); + gap *= 1e2; + if (mipsolver.options_mip_->objective_bound < ub) + ub = mipsolver.options_mip_->objective_bound; auto print_lp_iters = convertToPrintString(total_lp_iterations); if (upper_bound != kHighsInf) { - ub = upper_bound + offset; - - if (std::fabs(ub) <= epsilon) ub = 0; - lb = std::min(ub, lb); - if (ub == 0.0) - gap = lb == 0.0 ? 0.0 : kHighsInf; - else - gap = 100. * (ub - lb) / fabs(ub); - std::array gap_string = {}; if (gap >= 9999.) std::strcpy(gap_string.data(), "Large"); @@ -1574,7 +1648,6 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) { std::array ub_string; if (mipsolver.options_mip_->objective_bound < ub) { - ub = mipsolver.options_mip_->objective_bound; ub_string = convertToPrintString((int)mipsolver.orig_model_->sense_ * ub, "*"); } else @@ -1596,7 +1669,6 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) { } else { std::array ub_string; if (mipsolver.options_mip_->objective_bound < ub) { - ub = mipsolver.options_mip_->objective_bound; ub_string = convertToPrintString((int)mipsolver.orig_model_->sense_ * ub, "*"); } else @@ -1623,9 +1695,11 @@ void HighsMipSolverData::printDisplayLine(const int solution_source) { double primal_bound; double mip_rel_gap; limitsToBounds(dual_bound, primal_bound, mip_rel_gap); + mip_rel_gap *= 1e2; assert(dual_bound == (int)mipsolver.orig_model_->sense_ * lb); assert(primal_bound == (int)mipsolver.orig_model_->sense_ * ub); - assert(mip_rel_gap == gap); + assert(gap == mip_rel_gap); + // Possibly interrupt from MIP logging callback mipsolver.callback_->clearHighsCallbackDataOut(); const bool interrupt = interruptFromCallbackWithData( @@ -1665,7 +1739,14 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() { globalOrbits->orbitalFixing(domain); if (domain.infeasible()) { + double prev_lower_bound = lower_bound; + lower_bound = std::min(kHighsInf, upper_bound); + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); pruned_treeweight = 1.0; num_nodes += 1; num_leaves += 1; @@ -1707,7 +1788,15 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() { addIncumbent(lp.getLpSolver().getSolution().col_value, lp.getObjective(), kSolutionSourceEvaluateNode)) { mipsolver.modelstatus_ = HighsModelStatus::kOptimal; + + double prev_lower_bound = lower_bound; + lower_bound = upper_bound; + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); pruned_treeweight = 1.0; num_nodes += 1; num_leaves += 1; @@ -1717,7 +1806,14 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() { status = lp.getStatus(); if (status == HighsLpRelaxation::Status::kInfeasible) { + double prev_lower_bound = lower_bound; + lower_bound = std::min(kHighsInf, upper_bound); + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); pruned_treeweight = 1.0; num_nodes += 1; num_leaves += 1; @@ -1725,7 +1821,15 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() { } if (lp.unscaledDualFeasible(lp.getStatus())) { + double prev_lower_bound = lower_bound; + lower_bound = std::max(lp.getObjective(), lower_bound); + + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); + if (lpWasSolved) { redcostfixing.addRootRedcost(mipsolver, lp.getLpSolver().getSolution().col_dual, @@ -1776,8 +1880,16 @@ void HighsMipSolverData::evaluateRootNode() { lp.loadModel(); domain.clearChangedCols(); lp.setObjectiveLimit(upper_limit); + + double prev_lower_bound = lower_bound; + lower_bound = std::max(lower_bound, domain.getObjectiveLowerBound()); + bool bound_change = lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + updatePrimalDualIntegral(prev_lower_bound, lower_bound, upper_bound, + upper_bound); + printDisplayLine(); if (firstrootbasis.valid) @@ -2352,30 +2464,11 @@ void HighsMipSolverData::saveReportMipSolution(const double new_upper_limit) { void HighsMipSolverData::limitsToBounds(double& dual_bound, double& primal_bound, double& mip_rel_gap) const { - const HighsLp* model = this->mipsolver.model_; - const HighsLp* orig_model = this->mipsolver.orig_model_; - - const double offset = model->offset_; - dual_bound = lower_bound + offset; - if (std::abs(dual_bound) <= epsilon) dual_bound = 0; - primal_bound = kHighsInf; - mip_rel_gap = kHighsInf; - - if (upper_bound != kHighsInf) { - primal_bound = upper_bound + offset; - - if (std::fabs(primal_bound) <= epsilon) primal_bound = 0; - dual_bound = std::min(dual_bound, primal_bound); - if (primal_bound == 0.0) - mip_rel_gap = dual_bound == 0.0 ? 0.0 : kHighsInf; - else - mip_rel_gap = 100. * (primal_bound - dual_bound) / fabs(primal_bound); - } + mip_rel_gap = limitsToGap(lower_bound, upper_bound, dual_bound, primal_bound); primal_bound = std::min(mipsolver.options_mip_->objective_bound, primal_bound); - // Adjust objective sense in case of maximization problem - if (orig_model->sense_ == ObjSense::kMaximize) { + if (this->mipsolver.orig_model_->sense_ == ObjSense::kMaximize) { dual_bound = -dual_bound; primal_bound = -primal_bound; } @@ -2404,9 +2497,141 @@ bool HighsMipSolverData::interruptFromCallbackWithData( mipsolver.mipdata_->total_lp_iterations; mipsolver.callback_->data_out.mip_primal_bound = primal_bound; mipsolver.callback_->data_out.mip_dual_bound = dual_bound; - // Option mip_rel_gap, and mip_gap in HighsInfo, are both fractions, - // whereas mip_rel_gap in logging output (mimicked by - // limitsToBounds) gives a percentage, so convert it a fraction - mipsolver.callback_->data_out.mip_gap = 1e-2 * mip_rel_gap; + mipsolver.callback_->data_out.mip_gap = mip_rel_gap; return mipsolver.callback_->callbackAction(callback_type, message); } + +double possInfRelDiff(const double v0, const double v1, const double den) { + double rel_diff; + if (std::fabs(v0) == kHighsInf) { + if (std::fabs(v1) == kHighsInf) { + rel_diff = 0; + } else { + rel_diff = kHighsInf; + } + } else { + if (std::fabs(v1) == kHighsInf) { + rel_diff = kHighsInf; + } else { + rel_diff = std::fabs(v1 - v0) / std::max(1.0, std::fabs(den)); + } + } + return rel_diff; +} + +void HighsMipSolverData::updatePrimalDualIntegral(const double from_lower_bound, + const double to_lower_bound, + const double from_upper_bound, + const double to_upper_bound, + const bool check_bound_change, + const bool check_prev_data) { + // Parameters to updatePrimalDualIntegral are lower and upper bounds + // before/after a change + // + // updatePrimalDualIntegral should only be called when there is a + // change in one of the bounds, except when the final update is + // made, in which case the bounds must NOT have changed. By default, + // a check for some bound change is made, unless check_bound_change + // is false, in which case there is a check for unchanged bounds. + // + HighsPrimaDualIntegral& pdi = this->primal_dual_integral; + // HighsPrimaDualIntegral struct contains the following data + // + // * value: Current value of the P-D integral + // + // * prev_lb: Value of lb that was computed from to_lower_bound in + // the previous call. Used as a check that the value of lb + // computed from from_lower_bound in this call is equal - to + // within bound_change_tolerance. If not true, then a change in lb + // has been missed. Only for checking/debugging + // + // * prev_ub: Ditto for upper_bound. Only for checking/debugging + // + // * prev_gap: Ditto for gap. Only for checking/debugging + // + // * prev_time: Used to determine the time spent at the previous gap + + double from_lb; + double from_ub; + const double from_gap = + this->limitsToGap(from_lower_bound, from_upper_bound, from_lb, from_ub); + double to_lb; + double to_ub; + const double to_gap = + this->limitsToGap(to_lower_bound, to_upper_bound, to_lb, to_ub); + + const double lb_difference = possInfRelDiff(from_lb, to_lb, to_lb); + const double ub_difference = possInfRelDiff(from_ub, to_ub, to_ub); + const double bound_change_tolerance = 0; + const bool bound_change = lb_difference > bound_change_tolerance || + ub_difference > bound_change_tolerance; + + if (check_bound_change) { + if (!bound_change) { + if (from_lower_bound == to_lower_bound && + from_upper_bound == to_upper_bound) { + const double lower_bound_difference = + possInfRelDiff(from_lower_bound, to_lower_bound, to_lower_bound); + const double upper_bound_difference = + possInfRelDiff(from_upper_bound, to_upper_bound, to_upper_bound); + assert(bound_change); + } + } + } else { + if (bound_change) { + if (from_lower_bound != to_lower_bound || + from_upper_bound != to_upper_bound) { + const double lower_bound_difference = + possInfRelDiff(from_lower_bound, to_lower_bound, to_lower_bound); + const double upper_bound_difference = + possInfRelDiff(from_upper_bound, to_upper_bound, to_upper_bound); + assert(!bound_change); + } + } + } + if (pdi.value > -kHighsInf) { + // updatePrimalDualIntegral has been called previously, so can + // usually test housekeeping, even if gap is still inf + // + // The one case where the checking can't be done comes after restart, where + // the + // + if (check_prev_data) { + // These housekeeping tests check that the previous saved + // lower/upper bounds and gap are very close to the "from" + // lower/upper bounds and corresponding gap. They are usually + // identical, but rounding error can occur when passing through + // reset, when the old/new offsets are added/subtracted from the + // bounds due to changes in offset during presolve. + const double lb_inconsistency = + possInfRelDiff(from_lb, pdi.prev_lb, pdi.prev_lb); + const bool lb_consistent = lb_inconsistency < 1e-12; + const double ub_inconsistency = + possInfRelDiff(from_ub, pdi.prev_ub, pdi.prev_ub); + const bool ub_consistent = ub_inconsistency < 1e-12; + const double gap_inconsistency = + possInfRelDiff(from_gap, pdi.prev_gap, 1.0); + const bool gap_consistent = gap_inconsistency < 1e-12; + assert(lb_consistent); + assert(ub_consistent); + assert(gap_consistent); + } + if (to_gap < kHighsInf) { + double time = mipsolver.timer_.read(mipsolver.timer_.solve_clock); + if (from_gap < kHighsInf) { + // Need to update the P-D integral + double time_diff = time - pdi.prev_time; + assert(time_diff >= 0); + pdi.value += time_diff * pdi.prev_gap; + } + pdi.prev_time = time; + } + } else { + pdi.value = 0; + } + pdi.prev_lb = to_lb; + pdi.prev_ub = to_ub; + pdi.prev_gap = to_gap; +} + +void HighsPrimaDualIntegral::initialise() { this->value = -kHighsInf; } diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index bf7522b086..38693bf361 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -33,6 +33,15 @@ #include "presolve/HighsSymmetry.h" #include "util/HighsTimer.h" +struct HighsPrimaDualIntegral { + double value; + double prev_lb; + double prev_ub; + double prev_gap; + double prev_time; + void initialise(); +}; + enum MipSolutionSource : int { kSolutionSourceNone = -1, kSolutionSourceMin = kSolutionSourceNone, @@ -140,6 +149,8 @@ struct HighsMipSolverData { HighsNodeQueue nodequeue; + HighsPrimaDualIntegral primal_dual_integral; + HighsDebugSol debugSolution; HighsMipSolverData(HighsMipSolver& mipsolver) @@ -222,6 +233,15 @@ struct HighsMipSolverData { void finishSymmetryDetection(const highs::parallel::TaskGroup& taskGroup, std::unique_ptr& symData); + void updatePrimalDualIntegral(const double from_lower_bound, + const double to_lower_bound, + const double from_upper_bound, + const double to_upper_bound, + const bool check_bound_change = true, + const bool check_prev_data = true); + double limitsToGap(const double use_lower_bound, const double use_upper_bound, + double& lb, double& ub) const; + double computeNewUpperLimit(double upper_bound, double mip_abs_gap, double mip_rel_gap) const; bool moreHeuristicsAllowed() const; diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 4a306c09b8..5908f4373c 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -280,8 +280,18 @@ void HighsPrimalHeuristics::rootReducedCost() { localdom.propagate(); if (localdom.infeasible()) { localdom.conflictAnalysis(mipsolver.mipdata_->conflictPool); + + double prev_lower_bound = mipsolver.mipdata_->lower_bound; + mipsolver.mipdata_->lower_bound = std::max(mipsolver.mipdata_->lower_bound, currCutoff); + + const bool bound_change = + mipsolver.mipdata_->lower_bound != prev_lower_bound; + if (!mipsolver.submip && bound_change) + mipsolver.mipdata_->updatePrimalDualIntegral( + prev_lower_bound, mipsolver.mipdata_->lower_bound, + mipsolver.mipdata_->upper_bound, mipsolver.mipdata_->upper_bound); localdom.backtrack(); if (localdom.getBranchDepth() == 0) break; neighbourhood.backtracked();