always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[CP-SAT] fix best bound callback non C++ code; more work on pseudo-costs

Browse Source
This commit is contained in:
Laurent Perron
2024-05-02 14:04:47 +02:00
committed by Corentin Le Molgat
parent b54f5d422a
commit ec453aafb6
11 changed files with 194 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
ortools/java/com/google/ortools/sat/CpSolver.java
View File

@@ -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<Double> 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);

3
ortools/sat/cp_model_solver.h
View File

@@ -106,8 +106,7 @@ std::function<void(Model*)> 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<void(Model*)> NewBestBoundCallback(
const std::function<void(double)>& callback);

10
ortools/sat/csharp/CpSolver.cs
View File

@@ -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 {

92
ortools/sat/integer_search.cc
View File

@@ -233,20 +233,14 @@ std::function<BooleanOrIntegerLiteral()> BoolPseudoCostHeuristic(Model* model) {
std::function<BooleanOrIntegerLiteral()> LpPseudoCostHeuristic(Model* model) {
auto* lp_values = model->GetOrCreate<ModelLpValues>();
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* pseudo_costs = model->GetOrCreate<PseudoCosts>();
auto* encoder = model->GetOrCreate<IntegerEncoder>();
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<BooleanOrIntegerLiteral()> 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<BooleanOrIntegerLiteral()>& 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;
}

4
ortools/sat/integer_search.h
View File

@@ -281,6 +281,10 @@ class IntegerSearchHelper {
bool GetDecision(const std::function<BooleanOrIntegerLiteral()>& 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() {

12
ortools/sat/linear_programming_constraint.cc
View File

@@ -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<double>::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++) {

15
ortools/sat/linear_programming_constraint.h
View File

@@ -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<int>& 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<double> lp_solution_;
std::vector<double> 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.

1
ortools/sat/linear_propagation.h
View File

@@ -18,6 +18,7 @@
#include <deque>
#include <functional>
#include <limits>
#include <ostream>
#include <string>
#include <utility>

111
ortools/sat/pseudo_costs.cc
View File

@@ -14,7 +14,9 @@
#include "ortools/sat/pseudo_costs.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <limits>
#include <string>
#include <tuple>
@@ -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<double, double> 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<const double> 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<const double> 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<const Literal> 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::VariableBoundChange> PseudoCosts::GetBoundChanges(
Literal decision) {
Literal decision, absl::Span<const double> lp_values) {
std::vector<PseudoCosts::VariableBoundChange> 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);
}

32
ortools/sat/pseudo_costs.h
View File

@@ -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<const double> 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<double, double> 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<const double> 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<VariableBoundChange> GetBoundChanges(Literal decision);
std::vector<VariableBoundChange> GetBoundChanges(
Literal decision, absl::Span<const double> lp_values);
private:
// Returns the current objective info.

3
ortools/sat/synchronization.h
View File

@@ -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<void(double)> callback);