From e5d09e1e12c5c89b907b9423b0238e37d2e1a81a Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 2 May 2024 14:04:47 +0200 Subject: [PATCH] [CP-SAT] fix best bound callback non C++ code; more work on pseudo-costs --- .../java/com/google/ortools/sat/CpSolver.java | 21 +++- ortools/sat/cp_model_solver.h | 3 +- ortools/sat/csharp/CpSolver.cs | 10 ++ ortools/sat/integer_search.cc | 92 ++++++--------- ortools/sat/integer_search.h | 4 + ortools/sat/linear_programming_constraint.cc | 12 +- ortools/sat/linear_programming_constraint.h | 15 ++- ortools/sat/linear_propagation.h | 1 + ortools/sat/pseudo_costs.cc | 111 ++++++++++++------ ortools/sat/pseudo_costs.h | 32 +++-- ortools/sat/synchronization.h | 3 - 11 files changed, 194 insertions(+), 110 deletions(-) diff --git a/ortools/java/com/google/ortools/sat/CpSolver.java b/ortools/java/com/google/ortools/sat/CpSolver.java index d5f67fbf7b..ae336d7195 100644 --- a/ortools/java/com/google/ortools/sat/CpSolver.java +++ b/ortools/java/com/google/ortools/sat/CpSolver.java @@ -36,7 +36,7 @@ public final class CpSolver { /** Solves the given model, and returns the solve status. */ public CpSolverStatus solve(CpModel model) { - return solveWithSolutionCallback(model, null); + return solve(model, null); } /** @@ -97,9 +97,9 @@ public final class CpSolver { public CpSolverStatus searchAllSolutions(CpModel model, CpSolverSolutionCallback cb) { boolean oldValue = solveParameters.getEnumerateAllSolutions(); solveParameters.setEnumerateAllSolutions(true); - solve(model, cb); + CpSolverStatus status = solve(model, cb); solveParameters.setEnumerateAllSolutions(oldValue); - return solveResponse.getStatus(); + return status; } private synchronized void createSolveWrapper() { @@ -189,6 +189,21 @@ public final class CpSolver { this.logCallback = cb; } + /** Clears the log callback. */ + public void clearLogCallback() { + this.logCallback = null; + } + + /** Sets the best bound callback for the solver. */ + public void setBestBoundCallback(Consumer cb) { + this.bestBoundCallback = cb; + } + + /** Clears the best bound callback. */ + public void clearBestBoundCallback() { + this.bestBoundCallback = null; + } + /** Returns some statistics on the solution found as a string. */ public String responseStats() { return CpSatHelper.solverResponseStats(solveResponse); diff --git a/ortools/sat/cp_model_solver.h b/ortools/sat/cp_model_solver.h index 24fd56cee2..c4362e2061 100644 --- a/ortools/sat/cp_model_solver.h +++ b/ortools/sat/cp_model_solver.h @@ -106,8 +106,7 @@ std::function NewFeasibleSolutionLogCallback( /** Creates a callbacks that will be called on each new best objective bound * found. * - * Note that adding a callback is not free since some computation will be done - * before this is called, and it will only be called on optimization models. + * Note that this function is called before the update takes place. */ std::function NewBestBoundCallback( const std::function& callback); diff --git a/ortools/sat/csharp/CpSolver.cs b/ortools/sat/csharp/CpSolver.cs index bdf763b4b2..25e9538778 100644 --- a/ortools/sat/csharp/CpSolver.cs +++ b/ortools/sat/csharp/CpSolver.cs @@ -145,11 +145,21 @@ public class CpSolver log_callback_ = new LogCallbackDelegate(del); } + public void ClearLogCallback() + { + log_callback_ = null; + } + public void SetBestBoundCallback(DoubleToVoidDelegate del) { best_bound_callback_ = new BestBoundCallbackDelegate(del); } + public void ClearBestBoundCallback() + { + best_bound_callback_ = null; + } + public CpSolverResponse Response { get { diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index f95de4c846..f7aa8dd8bc 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -233,20 +233,14 @@ std::function BoolPseudoCostHeuristic(Model* model) { std::function LpPseudoCostHeuristic(Model* model) { auto* lp_values = model->GetOrCreate(); - auto* integer_trail = model->GetOrCreate(); auto* pseudo_costs = model->GetOrCreate(); - auto* encoder = model->GetOrCreate(); - return [lp_values, pseudo_costs, integer_trail, encoder, model]() { + return [lp_values, pseudo_costs]() { double best_score = 0.0; BooleanOrIntegerLiteral decision; for (IntegerVariable var(0); var < lp_values->size(); var += 2) { - const IntegerValue lb = integer_trail->LowerBound(var); - const IntegerValue ub = integer_trail->UpperBound(var); - if (lb == ub) continue; - - const double lp_value = (*lp_values)[var]; - const bool is_reliable = pseudo_costs->LpReliability(var) >= 4; - const bool is_integer = std::abs(lp_value - std::round(lp_value)) < 1e-6; + const PseudoCosts::BranchingInfo info = + pseudo_costs->EvaluateVar(var, *lp_values); + if (info.is_fixed) continue; // When not reliable, we skip integer. // @@ -254,35 +248,18 @@ std::function LpPseudoCostHeuristic(Model* model) { // TODO(user): do not branch on integer lp? however it seems better to // do that !? Maybe this is because if it has a high pseudo cost // average, it is good anyway? - if (!is_reliable && is_integer) continue; - - // There are some corner cases if we are at the bound. Note that it is - // important to be in sync with the SplitAroundLpValue() below. - double down_fractionality = lp_value - std::floor(lp_value); - IntegerValue down_target = IntegerValue(std::floor(lp_value)); - if (lp_value >= ToDouble(ub)) { - down_fractionality = 1.0; - down_target = ub - 1; - } else if (lp_value <= ToDouble(lb)) { - down_fractionality = 0.0; - down_target = lb; - } - const auto [down_score, up_score] = - pseudo_costs->LpPseudoCost(var, down_fractionality); - const double score = pseudo_costs->CombineScores(down_score, up_score); + if (!info.is_reliable && info.is_integer) continue; // We delay to subsequent heuristic if the score is 0.0. - if (score > best_score) { - best_score = score; + if (info.score > best_score) { + best_score = info.score; // This direction works better than the inverse in the benchs. But // always branching up seems even better. TODO(user): investigate. - if (down_score > up_score) { - decision = BooleanOrIntegerLiteral( - IntegerLiteral::LowerOrEqual(var, down_target)); + if (info.down_score > info.up_score) { + decision = BooleanOrIntegerLiteral(info.down_branch); } else { - decision = BooleanOrIntegerLiteral( - IntegerLiteral::GreaterOrEqual(var, down_target + 1)); + decision = BooleanOrIntegerLiteral(info.down_branch.Negated()); } } } @@ -1380,6 +1357,31 @@ bool IntegerSearchHelper::BeforeTakingDecision() { return true; } +LiteralIndex IntegerSearchHelper::GetDecisionLiteral( + const BooleanOrIntegerLiteral& decision) { + DCHECK(decision.HasValue()); + + // Convert integer decision to literal one if needed. + // + // TODO(user): Ideally it would be cool to delay the creation even more + // until we have a conflict with these decisions, but it is currently + // hard to do so. + const Literal literal = + decision.boolean_literal_index != kNoLiteralIndex + ? Literal(decision.boolean_literal_index) + : encoder_->GetOrCreateAssociatedLiteral(decision.integer_literal); + if (sat_solver_->Assignment().LiteralIsAssigned(literal)) { + // TODO(user): It would be nicer if this can never happen. For now, it + // does because of the Propagate() not reaching the fixed point as + // mentioned in a TODO above. As a work-around, we display a message + // but do not crash and recall the decision heuristic. + VLOG(1) << "Trying to take a decision that is already assigned!" + << " Fix this. Continuing for now..."; + return kNoLiteralIndex; + } + return literal.Index(); +} + bool IntegerSearchHelper::GetDecision( const std::function& f, LiteralIndex* decision) { *decision = kNoLiteralIndex; @@ -1411,28 +1413,8 @@ bool IntegerSearchHelper::GetDecision( } if (!new_decision.HasValue()) break; - // Convert integer decision to literal one if needed. - // - // TODO(user): Ideally it would be cool to delay the creation even more - // until we have a conflict with these decisions, but it is currently - // hard to do so. - if (new_decision.boolean_literal_index != kNoLiteralIndex) { - *decision = new_decision.boolean_literal_index; - } else { - *decision = - encoder_->GetOrCreateAssociatedLiteral(new_decision.integer_literal) - .Index(); - } - if (sat_solver_->Assignment().LiteralIsAssigned(Literal(*decision))) { - // TODO(user): It would be nicer if this can never happen. For now, it - // does because of the Propagate() not reaching the fixed point as - // mentioned in a TODO above. As a work-around, we display a message - // but do not crash and recall the decision heuristic. - VLOG(1) << "Trying to take a decision that is already assigned!" - << " Fix this. Continuing for now..."; - continue; - } - break; + *decision = GetDecisionLiteral(new_decision); + if (*decision != kNoLiteralIndex) break; } return true; } diff --git a/ortools/sat/integer_search.h b/ortools/sat/integer_search.h index 511b3d554d..d4722d8e59 100644 --- a/ortools/sat/integer_search.h +++ b/ortools/sat/integer_search.h @@ -281,6 +281,10 @@ class IntegerSearchHelper { bool GetDecision(const std::function& f, LiteralIndex* decision); + // Inner function used by GetDecision(). + // It will create a new associated literal if needed. + LiteralIndex GetDecisionLiteral(const BooleanOrIntegerLiteral& decision); + // Functions passed to GetDecision() might call this to notify a conflict // was detected. void NotifyThatConflictWasFoundDuringGetDecision() { diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 18a810cdc9..da9812a817 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -557,6 +557,11 @@ void LinearProgrammingConstraint::SetLevel(int level) { if (lp_solution_is_set_ && level < lp_solution_level_) { lp_solution_is_set_ = false; } + if (level < previous_level_) { + lp_at_optimal_ = false; + lp_objective_lower_bound_ = -std::numeric_limits::infinity(); + } + previous_level_ = level; // Special case for level zero, we "reload" any previously known optimal // solution from that level. @@ -697,6 +702,12 @@ bool LinearProgrammingConstraint::SolveLp() { << objective_definition_->ScaleObjective( simplex_.GetObjectiveValue()); + if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL || + simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_FEASIBLE) { + lp_objective_lower_bound_ = simplex_.GetObjectiveValue(); + } + lp_at_optimal_ = simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL; + if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { lp_solution_is_set_ = true; lp_solution_level_ = trail_->CurrentDecisionLevel(); @@ -799,7 +810,6 @@ bool LinearProgrammingConstraint::AnalyzeLp() { if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { CHECK(lp_solution_is_set_); - lp_objective_ = simplex_.GetObjectiveValue(); lp_solution_is_integer_ = true; const int num_vars = integer_variables_.size(); for (int i = 0; i < num_vars; i++) { diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 125d3ca856..2a41be232e 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -164,11 +164,17 @@ class LinearProgrammingConstraint : public PropagatorInterface, // Note that this solution is always an OPTIMAL solution of an LP above or // at the current decision level. We "erase" it when we backtrack over it. bool HasSolution() const { return lp_solution_is_set_; } - double SolutionObjectiveValue() const { return lp_objective_; } double GetSolutionValue(IntegerVariable variable) const; double GetSolutionReducedCost(IntegerVariable variable) const; bool SolutionIsInteger() const { return lp_solution_is_integer_; } + // Returns a valid lp lower bound for the current branch, and indicates if + // the latest LP solve in that branch was solved to optimality or not. + // During normal operation, we will always propagate the LP, so this should + // not refer to an optimality status lower in that branch. + bool AtOptimal() const { return lp_at_optimal_; } + double ObjectiveLpLowerBound() const { return lp_objective_lower_bound_; } + // PropagatorInterface API. bool Propagate() override; bool IncrementalPropagate(const std::vector& watch_indices) override; @@ -526,10 +532,15 @@ class LinearProgrammingConstraint : public PropagatorInterface, int lp_solution_level_ = 0; bool lp_solution_is_set_ = false; bool lp_solution_is_integer_ = false; - double lp_objective_; std::vector lp_solution_; std::vector lp_reduced_cost_; + // Last objective lower bound found by the LP Solver. + // We erase this on backtrack. + int previous_level_ = 0; + bool lp_at_optimal_ = false; + double lp_objective_lower_bound_; + // If non-empty, this is the last known optimal lp solution at root-node. If // the variable bounds changed, or cuts where added, it is possible that this // solution is no longer optimal though. diff --git a/ortools/sat/linear_propagation.h b/ortools/sat/linear_propagation.h index 1932ab499c..ba255d4927 100644 --- a/ortools/sat/linear_propagation.h +++ b/ortools/sat/linear_propagation.h @@ -18,6 +18,7 @@ #include #include +#include #include #include #include diff --git a/ortools/sat/pseudo_costs.cc b/ortools/sat/pseudo_costs.cc index 79a17e58de..2afe38e790 100644 --- a/ortools/sat/pseudo_costs.cc +++ b/ortools/sat/pseudo_costs.cc @@ -14,7 +14,9 @@ #include "ortools/sat/pseudo_costs.h" #include +#include #include +#include #include #include #include @@ -85,30 +87,70 @@ PseudoCosts::ObjectiveInfo PseudoCosts::GetCurrentObjectiveInfo() { result.lp_bound = 0.0; result.lp_at_optimal = true; for (const auto* lp : *lps_) { - if (!lp->HasSolution()) result.lp_at_optimal = false; - result.lp_bound += lp->SolutionObjectiveValue(); + if (!lp->AtOptimal()) result.lp_at_optimal = false; + result.lp_bound += lp->ObjectiveLpLowerBound(); } return result; } -void PseudoCosts::BeforeTakingDecision(Literal decision) { - if (objective_var_ == kNoIntegerVariable) return; +bool PseudoCosts::SaveLpInfo() { saved_info_ = GetCurrentObjectiveInfo(); - bound_changes_ = GetBoundChanges(decision); + return saved_info_.lp_at_optimal; } -std::pair PseudoCosts::LpPseudoCost( - IntegerVariable var, double down_fractionality) const { - const int max_index = std::max(var.value(), NegationOf(var).value()); - if (max_index >= average_unit_objective_increase_.size()) return {0.0, 0.0}; +void PseudoCosts::SaveBoundChanges(Literal decision, + absl::Span lp_values) { + bound_changes_ = GetBoundChanges(decision, lp_values); +} - const double up_fractionality = 1.0 - down_fractionality; - const double up_branch = - up_fractionality * average_unit_objective_increase_[var].CurrentAverage(); - const double down_branch = - down_fractionality * - average_unit_objective_increase_[NegationOf(var)].CurrentAverage(); - return {down_branch, up_branch}; +void PseudoCosts::BeforeTakingDecision(Literal decision) { + if (objective_var_ == kNoIntegerVariable) return; + SaveLpInfo(); + SaveBoundChanges(decision, *lp_values_); +} + +PseudoCosts::BranchingInfo PseudoCosts::EvaluateVar( + IntegerVariable var, absl::Span lp_values) { + DCHECK_NE(var, kNoIntegerVariable); + BranchingInfo result; + const IntegerValue lb = integer_trail_->LowerBound(var); + const IntegerValue ub = integer_trail_->UpperBound(var); + if (lb == ub) { + result.is_fixed = true; + return result; + } + + const double lp_value = lp_values[var.value()]; + double down_fractionality = lp_value - std::floor(lp_value); + IntegerValue down_target = IntegerValue(std::floor(lp_value)); + if (lp_value >= ToDouble(ub)) { + down_fractionality = 1.0; + down_target = ub - 1; + } else if (lp_value <= ToDouble(lb)) { + down_fractionality = 0.0; + down_target = lb; + } + + result.is_integer = std::abs(lp_value - std::round(lp_value)) < 1e-6; + result.down_fractionality = down_fractionality; + result.down_branch = IntegerLiteral::LowerOrEqual(var, down_target); + + const int max_index = std::max(var.value(), NegationOf(var).value()); + if (max_index < average_unit_objective_increase_.size()) { + result.down_score = + down_fractionality * + average_unit_objective_increase_[NegationOf(var)].CurrentAverage(); + result.up_score = (1.0 - down_fractionality) * + average_unit_objective_increase_[var].CurrentAverage(); + result.score = CombineScores(result.down_score, result.up_score); + + const int reliablitity = std::min( + average_unit_objective_increase_[var].NumRecords(), + average_unit_objective_increase_[NegationOf(var)].NumRecords()); + result.is_reliable = reliablitity >= 4; + } + + return result; } void PseudoCosts::UpdateBoolPseudoCosts(absl::Span reason, @@ -136,13 +178,19 @@ double PseudoCosts::BoolPseudoCost(Literal lit, double lp_value) const { return CombineScores(down_branch, up_branch); } -int PseudoCosts::LpReliability(IntegerVariable var) const { - const int max_index = std::max(var.value(), NegationOf(var).value()); - if (max_index >= average_unit_objective_increase_.size()) return 0; +double PseudoCosts::ObjectiveIncrease(bool conflict) { + const ObjectiveInfo new_info = GetCurrentObjectiveInfo(); + const double obj_lp_diff = + std::max(0.0, new_info.lp_bound - saved_info_.lp_bound); + const IntegerValue obj_int_diff = new_info.lb - saved_info_.lb; - return std::min( - average_unit_objective_increase_[var].NumRecords(), - average_unit_objective_increase_[NegationOf(var)].NumRecords()); + double obj_increase = + obj_lp_diff > 0.0 ? obj_lp_diff : ToDouble(obj_int_diff); + if (conflict) { + // We count a conflict as a max increase + 1.0 + obj_increase = ToDouble(saved_info_.ub) - ToDouble(saved_info_.lb) + 1.0; + } + return obj_increase; } void PseudoCosts::AfterTakingDecision(bool conflict) { @@ -156,23 +204,14 @@ void PseudoCosts::AfterTakingDecision(bool conflict) { // just be the "artificial" continuing of the current lp solve that create the // increase. if (saved_info_.lp_at_optimal) { - // Compute the increase in objective. - const double obj_lp_diff = - std::max(0.0, new_info.lp_bound - saved_info_.lp_bound); - const IntegerValue obj_int_diff = new_info.lb - saved_info_.lb; - double obj_diff = obj_lp_diff > 0.0 ? obj_lp_diff : ToDouble(obj_int_diff); - if (conflict) { - // We count a conflict as a max increase + 1.0 - obj_diff = ToDouble(saved_info_.ub) - ToDouble(saved_info_.lb) + 1.0; - } - // Update the average unit increases. + const double obj_increase = ObjectiveIncrease(conflict); for (const auto [var, lb_change, lp_increase] : bound_changes_) { if (lp_increase < 1e-6) continue; if (var >= average_unit_objective_increase_.size()) { average_unit_objective_increase_.resize(var + 1); } - average_unit_objective_increase_[var].AddData(obj_diff / lp_increase); + average_unit_objective_increase_[var].AddData(obj_increase / lp_increase); } } @@ -243,16 +282,16 @@ IntegerVariable PseudoCosts::GetBestDecisionVar() { } std::vector PseudoCosts::GetBoundChanges( - Literal decision) { + Literal decision, absl::Span lp_values) { std::vector bound_changes; for (const IntegerLiteral l : encoder_->GetIntegerLiterals(decision)) { PseudoCosts::VariableBoundChange entry; entry.var = l.var; entry.lower_bound_change = l.bound - integer_trail_->LowerBound(l.var); - if (l.var < lp_values_->size()) { + if (l.var < lp_values.size()) { entry.lp_increase = - std::max(0.0, ToDouble(l.bound) - (*lp_values_)[l.var]); + std::max(0.0, ToDouble(l.bound) - lp_values[l.var.value()]); } bound_changes.push_back(entry); } diff --git a/ortools/sat/pseudo_costs.h b/ortools/sat/pseudo_costs.h index a7b8c0447c..0f73294921 100644 --- a/ortools/sat/pseudo_costs.h +++ b/ortools/sat/pseudo_costs.h @@ -44,6 +44,12 @@ class PseudoCosts { // BeforeTakingDecision(). void AfterTakingDecision(bool conflict = false); + // Advanced usage. Internal functions used by Before/AfterTakingDecision(), + // that are exposed for strong branching. + double ObjectiveIncrease(bool conflict); + bool SaveLpInfo(); + void SaveBoundChanges(Literal decision, absl::Span lp_values); + // Returns the variable with best reliable pseudo cost that is not fixed. IntegerVariable GetBestDecisionVar(); @@ -60,15 +66,24 @@ class PseudoCosts { return pseudo_costs_[var].NumRecords(); } - // Alternative pseudo-costs. This relies on the LP more heavily and is more - // in line with what a MIP solver would do. Return the (down, up) costs which - // can be combined with CombineScores(); + // Combines the score of the two branch into one score. double CombineScores(double down_branch, double up_branch) const; - std::pair LpPseudoCost(IntegerVariable var, - double down_fractionality) const; - // Returns the pseudo cost "reliability". - int LpReliability(IntegerVariable var) const; + // Alternative pseudo-cost. This relies on the LP more heavily and is more in + // line with what a MIP solver would do. Returns all the info about taking a + // branch around the current lp_value of var. + struct BranchingInfo { + bool is_fixed = false; + bool is_reliable = false; + bool is_integer = false; + double down_fractionality = 0.0; + double score = 0.0; + double down_score = 0.0; + double up_score = 0.0; + IntegerLiteral down_branch; + }; + BranchingInfo EvaluateVar(IntegerVariable var, + absl::Span lp_values); // Experimental alternative pseudo cost based on the explanation for bound // increases. @@ -83,7 +98,8 @@ class PseudoCosts { IntegerValue lower_bound_change = IntegerValue(0); double lp_increase = 0.0; }; - std::vector GetBoundChanges(Literal decision); + std::vector GetBoundChanges( + Literal decision, absl::Span lp_values); private: // Returns the current objective info. diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index b24a651eb4..2ba4999651 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -272,9 +272,6 @@ class SharedResponseManager { // Adds a callback that will be called on each new best objective bound // found. Returns its id so it can be unregistered if needed. // - // Note that adding a callback is not free since some computation will be done - // before this is called. - // // Note that currently the class is waiting for the callback to finish before // accepting any new updates. That could be changed if needed. int AddBestBoundCallback(std::function callback);