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

[CP-SAT] improvements to feasibility_jump, better presolve on lin_max; tuning

Browse Source
This commit is contained in:
Laurent Perron
2023-04-21 12:48:03 +02:00
parent f4f9b622a5
commit dbaef62eb5
19 changed files with 1675 additions and 497 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
ortools/sat/BUILD.bazel
View File

@@ -176,6 +176,7 @@ cc_library(
":util",
"//ortools/base",
"//ortools/base:hash",
"//ortools/graph:strongly_connected_components",
"//ortools/port:proto_utils",
"//ortools/util:saturated_arithmetic",
"//ortools/util:sorted_interval_list",

974
ortools/sat/constraint_violation.cc
View File

File diff suppressed because it is too large Load Diff

271
ortools/sat/constraint_violation.h
View File

@@ -17,6 +17,7 @@
#include <cstddef>
#include <cstdint>
#include <memory>
#include <utility>
#include <vector>
#include "absl/types/span.h"
@@ -51,17 +52,39 @@ class LinearIncrementalEvaluator {
// and before the class starts to be used. This is DCHECKed.
void PrecomputeCompactView();
// Compute activities and query violations.
// Compute activities and update them.
void ComputeInitialActivities(absl::Span<const int64_t> solution);
void Update(int var, int64_t delta);
// Function specific to the feasibility jump heuristic.
// Note that the score of the changed variable will not be updated correctly!
void UpdateVariableAndScores(int var, int64_t delta,
absl::Span<const int64_t> solution,
absl::Span<const double> weights,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
void UpdateScoreOnWeightUpdate(int c, double weight_delta,
absl::Span<const int64_t> solution,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
// Variables whose score was updated since the last clear.
// Note that we also include variable whose score might be the same but for
// which different jump values might have a changed score.
void ClearAffectedVariables();
const std::vector<int>& VariablesAffectedByLastUpdate() const {
return last_affected_variables_;
}
// Query violation.
int64_t Activity(int c) const;
int64_t Violation(int c) const;
bool IsViolated(int c) const;
bool AppearsInViolatedConstraints(int var) const;
// Used to DCHECK the state of the evaluator.
bool VarIsConsistent(int var) const;
bool AppearsInViolatedConstraints(int var) const;
// Intersect constraint bounds with [lb..ub].
void ReduceBounds(int c, int64_t lb, int64_t ub);
@@ -77,8 +100,8 @@ class LinearIncrementalEvaluator {
double WeightedViolationDelta(absl::Span<const double> weights, int var,
int64_t delta) const;
// The violation for each constraints is a piecewise linear function. This
// computes and aggregate all the breakpoints for the given variable and its
// The violation for each constraint is a piecewise linear function. This
// computes and aggregates all the breakpoints for the given variable and its
// domain.
//
// Note that if the domain contains less than two values, we return an empty
@@ -86,13 +109,8 @@ class LinearIncrementalEvaluator {
std::vector<int64_t> SlopeBreakpoints(int var, int64_t current_value,
const Domain& var_domain) const;
// AffectedVariableOnUpdate() return the set of variables that have at least
// one constraint in common with the given variable. One need to call
// PrecomputeCompactView() after all the constraint have been added for
// this to work.
std::vector<int> AffectedVariableOnUpdate(int var);
std::vector<int> ConstraintsToAffectedVariables(
absl::Span<const int> ct_indices);
// Checks if the jump value of a variable is always optimal.
bool ViolationChangeIsConvex(int var) const;
double DeterministicTime() const {
return 5e-9 * static_cast<double>(dtime_);
@@ -111,6 +129,50 @@ class LinearIncrementalEvaluator {
bool positive; // bool_var or its negation.
};
struct SpanData {
int start = 0;
int num_pos_literal = 0;
int num_neg_literal = 0;
int linear_start = 0;
int num_linear_entries = 0;
};
absl::Span<const int> VarToConstraints(int var) const {
if (var >= columns_.size()) return {};
const SpanData& data = columns_[var];
const int size =
data.num_pos_literal + data.num_neg_literal + data.num_linear_entries;
if (size == 0) return {};
return absl::MakeSpan(&ct_buffer_[data.start], size);
}
absl::Span<const int> ConstraintToVars(int c) const {
const SpanData& data = rows_[c];
const int size =
data.num_pos_literal + data.num_neg_literal + data.num_linear_entries;
if (size == 0) return {};
return absl::MakeSpan(&row_var_buffer_[data.start], size);
}
void ComputeAndCacheDistance(int ct_index);
// Incremental row-based update.
void UpdateScoreOnNewlyEnforced(int c, double weight,
absl::Span<const int64_t> solution,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
void UpdateScoreOnNewlyUnenforced(int c, double weight,
absl::Span<const int64_t> solution,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
void UpdateScoreOfEnforcementIncrease(int c, double score_change,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
void UpdateScoreOnActivityChange(int c, double weight, int64_t activity_delta,
absl::Span<const int64_t> solution,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
// Constraint indexed data (static).
int num_constraints_ = 0;
std::vector<Domain> domains_;
@@ -119,41 +181,40 @@ class LinearIncrementalEvaluator {
// Variable indexed data.
// Note that this is just used at construction and is replaced by a compact
// view when PrecomputeCompactView() is called.
bool creation_phase_ = true;
std::vector<std::vector<Entry>> var_entries_;
std::vector<std::vector<LiteralEntry>> literal_entries_;
// Memory efficient column based data.
struct VarData {
int start = 0;
int num_pos_literal = 0;
int num_neg_literal = 0;
int linear_start = 0;
int num_linear_entries = 0;
};
absl::Span<const int> VarToConstraints(int var) const {
const VarData& data = columns_[var];
const int size =
data.num_pos_literal + data.num_neg_literal + data.num_linear_entries;
return absl::MakeSpan(&ct_buffer_[data.start], size);
}
std::vector<VarData> columns_;
// Memory efficient column based data (static).
std::vector<SpanData> columns_;
std::vector<int> ct_buffer_;
std::vector<int64_t> coeff_buffer_;
// Transposed view data (only variables).
bool creation_phase_ = true;
std::vector<int> row_buffer_;
std::vector<absl::Span<const int>> rows_;
// Memory efficient row based data (static).
std::vector<SpanData> rows_;
std::vector<int> row_var_buffer_;
std::vector<int64_t> row_coeff_buffer_;
// Temporary data.
std::vector<int> tmp_row_sizes_;
std::vector<bool> tmp_in_result_;
std::vector<int> tmp_row_num_positive_literals_;
std::vector<int> tmp_row_num_negative_literals_;
std::vector<int> tmp_row_num_linear_entries_;
// Constraint indexed data (dynamic).
std::vector<bool> is_violated_;
std::vector<int64_t> activities_;
std::vector<int64_t> distances_;
std::vector<int> num_false_enforcement_;
std::vector<bool> row_update_will_not_impact_deltas_;
// Code to update the scores on a variable change.
std::vector<int64_t> old_distances_;
std::vector<int> old_num_false_enforcement_;
std::vector<int64_t> cached_deltas_;
std::vector<double> cached_scores_;
std::vector<bool> in_last_affected_variables_;
std::vector<int> last_affected_variables_;
mutable size_t dtime_ = 0;
};
@@ -167,22 +228,31 @@ class CompiledConstraint {
explicit CompiledConstraint(const ConstraintProto& ct_proto);
virtual ~CompiledConstraint() = default;
// Recomputes the violation of the constraint from scratch.
void InitializeViolation(absl::Span<const int64_t> solution);
// Update the violation with the new value.
void PerformMove(int var, int64_t old_value,
absl::Span<const int64_t> solution_with_new_value);
// Computes the violation of a constraint.
//
// A violation is a positive integer value. A zero value means the constraint
// is not violated..
virtual int64_t Violation(absl::Span<const int64_t> solution) = 0;
// is not violated.
virtual int64_t ComputeViolation(absl::Span<const int64_t> solution) = 0;
void CacheViolation(absl::Span<const int64_t> solution);
// Returns the delta if var changes from old_value to solution[var].
virtual int64_t ViolationDelta(
int var, int64_t old_value,
absl::Span<const int64_t> solution_with_new_value);
// Cache the violation of a constraint.
int64_t violation() const;
const ConstraintProto& ct_proto() const;
// Getters.
const ConstraintProto& ct_proto() const { return ct_proto_; }
int64_t violation() const { return violation_; }
private:
const ConstraintProto& ct_proto_;
int64_t violation_ = -1;
int64_t violation_;
};
// Evaluation container for the local search.
@@ -205,15 +275,24 @@ class LsEvaluator {
// Intersects the domain of the objective with [lb..ub].
void ReduceObjectiveBounds(int64_t lb, int64_t ub);
// Sets the current solution, and compute violations for each constraints.
// Sets the current solution, and computes violations for each constraint.
void ComputeInitialViolations(absl::Span<const int64_t> solution);
// Recompute the violations of non linear constraints.
void UpdateNonLinearViolations();
void UpdateAllNonLinearViolations();
// Update the value of one variable and recompute violations.
void UpdateVariableValueAndRecomputeViolations(int var, int64_t value,
bool focus_on_linear = false);
// Sets the value of the variable in the current solution.
// It must be called after UpdateLinearScores().
void UpdateVariableValue(int var, int64_t new_value);
// Recomputes the violations of all impacted non linear constraints.
void UpdateNonLinearViolations(int var, int64_t new_value);
// Function specific to the linear only feasibility jump.
void UpdateLinearScores(int var, int64_t value,
absl::Span<const double> weights,
absl::Span<const int64_t> jump_values,
absl::Span<double> jump_scores);
// Simple summation metric for the constraint and objective violations.
int64_t SumOfViolations();
@@ -231,9 +310,13 @@ class LsEvaluator {
// Returns the weighted sum of violation. The weights must have the same
// size as NumEvaluatorConstraints().
int64_t Violation(int c) const;
bool IsViolated(int c) const;
double WeightedViolation(absl::Span<const double> weights) const;
double WeightedViolationDelta(absl::Span<const double> weights, int var,
int64_t delta) const;
// Ignores the violations of the linear constraints.
double WeightedNonLinearViolationDelta(absl::Span<const double> weights,
int var, int64_t delta) const;
LinearIncrementalEvaluator* MutableLinearEvaluator() {
return &linear_evaluator_;
@@ -242,11 +325,10 @@ class LsEvaluator {
// Returns the set of variables appearing in a violated constraint.
std::vector<int> VariablesInViolatedConstraints() const;
// Counters: number of times UpdateVariableValueAndRecomputeViolations() has
// been called.
int64_t num_deltas_computed() const { return num_deltas_computed_; }
// Indicates if the computed jump value is always the best choice.
bool VariableOnlyInLinearConstraintWithConvexViolationChange(int var) const;
// Access the solution stored. Useful for debugging.
// Access the solution stored.
const std::vector<int64_t>& current_solution() const {
return current_solution_;
}
@@ -259,10 +341,97 @@ class LsEvaluator {
LinearIncrementalEvaluator linear_evaluator_;
std::vector<std::unique_ptr<CompiledConstraint>> constraints_;
std::vector<std::vector<int>> var_to_constraint_graph_;
std::vector<bool> jump_value_optimal_;
// We need the mutable to evaluate a move.
mutable std::vector<int64_t> current_solution_;
bool model_is_supported_ = true;
mutable int64_t num_deltas_computed_ = 0;
};
// Individual compiled constraints.
// The violation of a bool_xor constraint is 0 or 1.
class CompiledBoolXorConstraint : public CompiledConstraint {
public:
explicit CompiledBoolXorConstraint(const ConstraintProto& ct_proto);
~CompiledBoolXorConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
int64_t ViolationDelta(
int /*var*/, int64_t /*old_value*/,
absl::Span<const int64_t> solution_with_new_value) override;
};
// The violation of a lin_max constraint is:
// - the sum(max(0, expr_value - target_value) forall expr). This part will be
// maintained by the linear part.
// - target_value - max(expressions) if positive.
class CompiledLinMaxConstraint : public CompiledConstraint {
public:
explicit CompiledLinMaxConstraint(const ConstraintProto& ct_proto);
~CompiledLinMaxConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
};
// The violation of an int_prod constraint is
// abs(value(target) - prod(value(expr)).
class CompiledIntProdConstraint : public CompiledConstraint {
public:
explicit CompiledIntProdConstraint(const ConstraintProto& ct_proto);
~CompiledIntProdConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
};
// The violation of an int_div constraint is
// abs(value(target) - value(expr0) / value(expr1)).
class CompiledIntDivConstraint : public CompiledConstraint {
public:
explicit CompiledIntDivConstraint(const ConstraintProto& ct_proto);
~CompiledIntDivConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
};
// The violation of a all_diff is the number of unordered pairs of expressions
// with the same value.
class CompiledAllDiffConstraint : public CompiledConstraint {
public:
explicit CompiledAllDiffConstraint(const ConstraintProto& ct_proto);
~CompiledAllDiffConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
private:
std::vector<int64_t> values_;
};
// The violation of a no_overlap is the sum of overloads over time.
class CompiledNoOverlapConstraint : public CompiledConstraint {
public:
explicit CompiledNoOverlapConstraint(const ConstraintProto& ct_proto,
const CpModelProto& cp_model);
~CompiledNoOverlapConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
private:
const CpModelProto& cp_model_;
std::vector<std::pair<int64_t, int64_t>> events_;
};
// The violation of a cumulative is the sum of overloads over time.
class CompiledCumulativeConstraint : public CompiledConstraint {
public:
explicit CompiledCumulativeConstraint(const ConstraintProto& ct_proto,
const CpModelProto& cp_model);
~CompiledCumulativeConstraint() override = default;
int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
private:
const CpModelProto& cp_model_;
std::vector<std::pair<int64_t, int64_t>> events_;
};
} // namespace sat

10
ortools/sat/cp_model_postsolve.cc
View File

@@ -207,9 +207,8 @@ int64_t EvaluateLinearExpression(const LinearExpressionProto& expr,
} // namespace
// Compute the max of each expression, and assign it to the target expr (which
// must be of the form +ref or -ref);
// We only support post-solving the case were the target is unassigned,
// but everything else is fixed.
// must be of the form +ref or -ref); We only support post-solving the case
// where all expression are fixed and correct.
void PostsolveLinMax(const ConstraintProto& ct, std::vector<Domain>* domains) {
int64_t max_value = std::numeric_limits<int64_t>::min();
for (const LinearExpressionProto& expr : ct.lin_max().exprs()) {
@@ -217,9 +216,8 @@ void PostsolveLinMax(const ConstraintProto& ct, std::vector<Domain>* domains) {
}
const int target_ref = GetSingleRefFromExpression(ct.lin_max().target());
const int target_var = PositiveRef(target_ref);
(*domains)[target_var] = (*domains)[target_var].IntersectionWith(
Domain(RefIsPositive(target_ref) ? max_value : -max_value));
CHECK(!(*domains)[target_var].IsEmpty());
(*domains)[target_var] =
Domain(RefIsPositive(target_ref) ? max_value : -max_value);
}
// We only support 3 cases in the presolve currently.

62
ortools/sat/cp_model_presolve.cc
View File

@@ -14,11 +14,9 @@
#include "ortools/sat/cp_model_presolve.h"
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <cstdlib>
#include <deque>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
@@ -32,14 +30,15 @@
#include "absl/container/flat_hash_map.h"
#include "absl/container/flat_hash_set.h"
#include "absl/hash/hash.h"
#include "absl/log/check.h"
#include "absl/numeric/int128.h"
#include "absl/strings/str_cat.h"
#include "absl/types/span.h"
#include "ortools/base/integral_types.h"
#include "ortools/base/logging.h"
#include "ortools/base/mathutil.h"
#include "ortools/base/stl_util.h"
#include "ortools/base/timer.h"
#include "ortools/base/vlog_is_on.h"
#include "ortools/graph/topologicalsorter.h"
#include "ortools/sat/circuit.h"
#include "ortools/sat/clause.h"
@@ -869,6 +868,61 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) {
if (abort) return changed;
}
// Checks if the affine target domain is constraining.
bool affine_target_domain_contains_max_domain = false;
if (ExpressionContainsSingleRef(target)) { // target = +/- var.
int64_t infered_min = std::numeric_limits<int64_t>::min();
int64_t infered_max = std::numeric_limits<int64_t>::min();
for (const LinearExpressionProto& expr : ct->lin_max().exprs()) {
infered_min = std::max(infered_min, context_->MinOf(expr));
infered_max = std::max(infered_max, context_->MaxOf(expr));
}
Domain rhs_domain;
for (const LinearExpressionProto& expr : ct->lin_max().exprs()) {
rhs_domain = rhs_domain.UnionWith(
context_->DomainSuperSetOf(expr).IntersectionWith(
{infered_min, infered_max}));
}
// Checks if all values from the max(exprs) belong in the domain of the
// target.
// Note that the target is +/-var.
DCHECK_EQ(std::abs(target.coeffs(0)), 1);
const Domain target_domain =
target.coeffs(0) == 1 ? context_->DomainOf(target.vars(0))
: context_->DomainOf(target.vars(0)).Negation();
affine_target_domain_contains_max_domain =
rhs_domain.IsIncludedIn(target_domain);
}
// If the target is not used, and safe, we can remove the constraint.
if (affine_target_domain_contains_max_domain &&
context_->VariableIsUniqueAndRemovable(target.vars(0))) {
context_->UpdateRuleStats("lin_max: unused affine target");
context_->MarkVariableAsRemoved(target.vars(0));
*context_->mapping_model->add_constraints() = *ct;
return RemoveConstraint(ct);
}
// If the target is only used in the objective, and safe, we can simplify the
// constraint.
if (affine_target_domain_contains_max_domain &&
context_->VariableWithCostIsUniqueAndRemovable(target.vars(0)) &&
(target.coeffs(0) > 0) ==
(context_->ObjectiveCoeff(target.vars(0)) > 0)) {
context_->UpdateRuleStats("lin_max: rewrite with precedences");
for (const LinearExpressionProto& expr : ct->lin_max().exprs()) {
LinearConstraintProto* prec =
context_->working_model->add_constraints()->mutable_linear();
prec->add_domain(0);
prec->add_domain(std::numeric_limits<int64_t>::max());
AddLinearExpressionToLinearConstraint(target, 1, prec);
AddLinearExpressionToLinearConstraint(expr, -1, prec);
}
*context_->mapping_model->add_constraints() = *ct;
return RemoveConstraint(ct);
}
// Deal with fixed target case.
if (target_min == target_max) {
bool all_booleans = true;
@@ -5111,7 +5165,7 @@ bool CpModelPresolver::PresolveNoOverlap(ConstraintProto* ct) {
return changed;
}
bool CpModelPresolver::PresolveNoOverlap2D(int c, ConstraintProto* ct) {
bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
if (context_->ModelIsUnsat()) {
return false;
}

15
ortools/sat/cp_model_solver.cc
View File

@@ -40,7 +40,6 @@
#include "absl/container/btree_set.h"
#include "absl/container/flat_hash_set.h"
#include "absl/flags/flag.h"
#include "absl/status/status.h"
#include "absl/strings/str_cat.h"
#include "absl/strings/str_format.h"
#include "absl/strings/str_join.h"
@@ -3202,7 +3201,11 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
absl::StrCat("rens_lns_", local_params.name())),
local_params, helper, &shared));
}
if (params.use_violation_ls() && !params.interleave_search()) {
const bool feasibility_jump_possible =
!params.interleave_search() &&
helper->TypeToConstraints(ConstraintProto::kNoOverlap2D).empty();
if (params.use_violation_ls() && feasibility_jump_possible) {
SatParameters local_params = params;
local_params.set_random_seed(params.random_seed());
local_params.set_feasibility_jump_decay(0.95);
@@ -3241,10 +3244,10 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
// schedule more than the available number of threads. They will just be
// interleaved. We will get an higher diversity, but use more memory.
const int num_feasibility_jump =
params.interleave_search()
? 0
: (params.test_feasibility_jump() ? num_available
: num_available / 2);
feasibility_jump_possible
? (params.test_feasibility_jump() ? num_available
: num_available / 2)
: 0;
const int num_first_solution_subsolvers =
num_available - num_feasibility_jump;
for (int i = 0; i < num_feasibility_jump; ++i) {

592
ortools/sat/feasibility_jump.cc
View File

@@ -16,6 +16,7 @@
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <limits>
#include <memory>
@@ -23,6 +24,8 @@
#include <utility>
#include <vector>
#include "absl/log/check.h"
#include "absl/random/bit_gen_ref.h"
#include "absl/random/random.h"
#include "absl/strings/str_cat.h"
#include "absl/types/span.h"
@@ -32,10 +35,31 @@
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_checker.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/restart.h"
#include "ortools/sat/subsolver.h"
#include "ortools/sat/synchronization.h"
#include "ortools/sat/util.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research::sat {
FeasibilityJumpSolver::~FeasibilityJumpSolver() {
if (!VLOG_IS_ON(1)) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back({"fs_jump/num_general_moves_computed", num_general_evals_});
stats.push_back({"fs_jump/num_general_moves_done", num_general_moves_});
stats.push_back({"fs_jump/num_linear_moves_computed", num_linear_evals_});
stats.push_back({"fs_jump/num_linear_moves_repaired_with_full_scan",
num_repairs_with_full_scan_});
stats.push_back({"fs_jump/num_linear_moves_done", num_linear_moves_});
stats.push_back({"fs_jump/num_solutions_imported", num_solutions_imported_});
stats.push_back(
{"fs_jump/num_variables_partially_scanned", num_partial_scans_});
stats.push_back({"fs_jump/num_weight_updates", num_weight_updates_});
shared_stats_->AddStats(stats);
}
void FeasibilityJumpSolver::Initialize() {
is_initialized_ = true;
evaluator_ = std::make_unique<LsEvaluator>(cp_model_);
@@ -46,18 +70,85 @@ void FeasibilityJumpSolver::Initialize() {
// validation.
const int num_variables = cp_model_.variables().size();
var_domains_.resize(num_variables);
var_has_two_values_.resize(num_variables);
for (int var = 0; var < num_variables; ++var) {
var_domains_[var] = ReadDomainFromProto(cp_model_.variables(var));
var_has_two_values_[var] = var_domains_[var].HasTwoValues();
}
}
namespace {
int64_t ComputeRange(int64_t range, double range_ratio) {
return static_cast<int64_t>(
std::ceil(static_cast<double>(range) * range_ratio));
}
// TODO(user): Optimize and move to the Domain class.
// TODO(user): Improve entropy on non continuous domains.
int64_t RandomValueNearMin(const Domain& domain, double range_ratio,
absl::BitGenRef random) {
if (domain.Size() == 1) return domain.FixedValue();
if (domain.Size() == 2) {
return absl::Bernoulli(random, 1 - range_ratio) ? domain.Min()
: domain.Max();
}
const int64_t range = ComputeRange(domain.Max() - domain.Min(), range_ratio);
return domain.ValueAtOrBefore(domain.Min() +
absl::LogUniform<int64_t>(random, 0, range));
}
int64_t RandomValueNearMax(const Domain& domain, double range_ratio,
absl::BitGenRef random) {
if (domain.Size() == 1) return domain.FixedValue();
if (domain.Size() == 2) {
return absl::Bernoulli(random, 1 - range_ratio) ? domain.Max()
: domain.Min();
}
const int64_t range = ComputeRange(domain.Max() - domain.Min(), range_ratio);
return domain.ValueAtOrAfter(domain.Max() -
absl::LogUniform<int64_t>(random, 0, range));
}
int64_t RandomValueNearZero(const Domain& domain, double range_ratio,
absl::BitGenRef random) {
if (domain.Min() >= 0) {
return RandomValueNearMin(domain, range_ratio, random);
}
if (domain.Max() <= 0) {
return RandomValueNearMax(domain, range_ratio, random);
}
const Domain positive_domain = domain.IntersectionWith({0, domain.Max()});
const double choose_positive_probability =
static_cast<double>(positive_domain.Size()) /
static_cast<double>(domain.Size());
if (absl::Bernoulli(random, choose_positive_probability)) {
return RandomValueNearMin(positive_domain, range_ratio, random);
} else {
return RandomValueNearMax(domain.IntersectionWith({domain.Min(), 0}),
range_ratio, random);
}
}
} // namespace
void FeasibilityJumpSolver::RestartFromDefaultSolution() {
// Starts with solution closest to zero.
// TODO(user): randomize the start instead?
const int num_variables = cp_model_.variables().size();
const double default_value_probability =
1.0 - params_.feasibility_jump_var_randomization_ratio();
const double pertubation_ratio =
params_.feasibility_jump_var_perburbation_range_ratio();
// Starts with values closest to zero.
std::vector<int64_t> solution(num_variables);
for (int var = 0; var < num_variables; ++var) {
solution[var] = var_domains_[var].SmallestValue();
if (num_batches_ == 0 ||
absl::Bernoulli(random_, default_value_probability)) {
solution[var] = var_domains_[var].SmallestValue();
} else {
solution[var] =
RandomValueNearZero(var_domains_[var], pertubation_ratio, random_);
}
}
if (cp_model_.has_objective() &&
@@ -66,9 +157,21 @@ void FeasibilityJumpSolver::RestartFromDefaultSolution() {
for (int i = 0; i < num_terms; ++i) {
const int var = cp_model_.objective().vars(i);
if (cp_model_.objective().coeffs(i) > 0) {
solution[var] = var_domains_[var].Min();
if (num_batches_ == 0 ||
absl::Bernoulli(random_, default_value_probability)) {
solution[var] = var_domains_[var].Min();
} else {
solution[var] =
RandomValueNearMin(var_domains_[var], pertubation_ratio, random_);
}
} else {
solution[var] = var_domains_[var].Max();
if (num_batches_ == 0 ||
absl::Bernoulli(random_, default_value_probability)) {
solution[var] = var_domains_[var].Max();
} else {
solution[var] =
RandomValueNearMax(var_domains_[var], pertubation_ratio, random_);
}
}
}
}
@@ -79,13 +182,33 @@ void FeasibilityJumpSolver::RestartFromDefaultSolution() {
}
std::string FeasibilityJumpSolver::OneLineStats() const {
return absl::StrCat(
"batch:", num_batches_, " restart:", FormatCounter(num_restarts_),
" inf:", FormatCounter(evaluator_->NumInfeasibleConstraints()),
" improving:", FormatCounter(good_jumps_.size()),
" jumps:", FormatCounter(num_jumps_),
" updates:", FormatCounter(num_weight_updates_),
" non_linear_evals:", FormatCounter(evaluator_->num_deltas_computed()));
// Restarts or solutions imported.
const std::string restart_str =
num_solutions_imported_ == 0
? absl::StrCat(" #restarts:", num_restarts_)
: absl::StrCat(" #solutions_imported:", num_solutions_imported_);
// Moves and evaluations in the general iterations.
const std::string general_str =
num_general_evals_ == 0 && num_general_moves_ == 0
? ""
: absl::StrCat(" #gen_moves:", FormatCounter(num_general_moves_), "/",
FormatCounter(num_general_evals_));
// Improving jumps and infeasible constraints.
const int num_infeasible_cts = evaluator_->NumInfeasibleConstraints();
const std::string non_solution_str =
good_jumps_.empty() && num_infeasible_cts == 0
? ""
: absl::StrCat(" #good_lin_moves:", FormatCounter(good_jumps_.size()),
" #inf_cts:",
FormatCounter(evaluator_->NumInfeasibleConstraints()));
return absl::StrCat("batch:", num_batches_, restart_str,
" #lin_moves:", FormatCounter(num_linear_moves_), "/",
FormatCounter(num_linear_evals_), general_str,
non_solution_str,
" #weight_updates:", FormatCounter(num_weight_updates_));
}
std::function<void()> FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) {
@@ -105,21 +228,36 @@ std::function<void()> FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) {
CHECK_GT(repo.NumSolutions(), 0);
const SharedSolutionRepository<int64_t>::Solution solution =
repo.GetRandomBiasedSolution(random_);
evaluator_->ComputeInitialViolations(solution.variable_values);
if (solution.rank < last_solution_rank_) {
evaluator_->ComputeInitialViolations(solution.variable_values);
// Reset weights for each new solution.
weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0);
// Reset weights for each new solution.
weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0);
// Update last solution rank.
last_solution_rank_ = solution.rank;
VLOG(2) << name() << " import a solution with value " << solution.rank;
++num_solutions_imported_;
}
} else {
// Restart? Note that we always "restart" the first time.
const double dtime =
evaluator_->MutableLinearEvaluator()->DeterministicTime();
if (dtime >= dtime_restart_threshold_) {
if (dtime >= dtime_restart_threshold_ &&
num_weight_updates_ >= update_restart_threshold_) {
++num_restarts_;
RestartFromDefaultSolution();
// We use luby restart with a base of 0.5s.
// We use luby restart with a base of 1 deterministic unit.
// We also block the restart if there was not enough weight update.
// Note that we only restart between batches too.
dtime_restart_threshold_ = dtime + 0.5 * SUniv(num_restarts_);
//
// TODO(user): Ideally batch should use deterministic time too so we
// can just use number of batch for the luby restart.
// TODO(user): Maybe have one worker with very low restart
// rate.
dtime_restart_threshold_ = dtime + 1.0 * SUniv(num_restarts_);
update_restart_threshold_ = num_weight_updates_ + 10;
}
}
@@ -133,6 +271,9 @@ std::function<void()> FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) {
// chunks.
if (shared_bounds_ != nullptr) {
shared_bounds_->UpdateDomains(&var_domains_);
for (int var = 0; var < var_domains_.size(); ++var) {
var_has_two_values_[var] = var_domains_[var].HasTwoValues();
}
}
// Checks the current solution is compatible with updated domains.
@@ -191,21 +332,21 @@ bool IsGood(double delta) { return delta < 0.0; }
void FeasibilityJumpSolver::RecomputeJump(int var) {
const std::vector<int64_t>& solution = evaluator_->current_solution();
++num_computed_jumps_;
++num_linear_evals_;
jump_need_recomputation_[var] = false;
if (var_domains_[var].IsFixed()) {
jump_values_[var] = solution[var];
jump_deltas_[var] = 0.0;
jump_scores_[var] = 0.0;
return;
}
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
if (var_domains_[var].Size() == 2) {
if (var_has_two_values_[var]) {
jump_values_[var] = solution[var] == var_domains_[var].Min()
? var_domains_[var].Max()
: var_domains_[var].Min();
jump_deltas_[var] = linear_evaluator->WeightedViolationDelta(
jump_scores_[var] = linear_evaluator->WeightedViolationDelta(
weights_, var, jump_values_[var] - solution[var]);
} else {
// In practice, after a few iterations, the chance of finding an improving
@@ -274,21 +415,10 @@ void FeasibilityJumpSolver::RecomputeJump(int var) {
DCHECK_NE(best_jump.first, solution[var]);
jump_values_[var] = best_jump.first;
jump_deltas_[var] = best_jump.second;
jump_scores_[var] = best_jump.second;
}
if (IsGood(jump_deltas_[var]) && !in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
}
void FeasibilityJumpSolver::MarkJumpForRecomputation(int var) {
jump_need_recomputation_[var] = true;
// This jump might be good, so we need to add it to the queue so it can be
// evaluated when choosing the next jump.
if (!in_good_jumps_[var]) {
if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
@@ -297,7 +427,7 @@ void FeasibilityJumpSolver::MarkJumpForRecomputation(int var) {
void FeasibilityJumpSolver::RecomputeAllJumps() {
const int num_variables = static_cast<int>(var_domains_.size());
jump_values_.resize(num_variables);
jump_deltas_.resize(num_variables);
jump_scores_.resize(num_variables);
jump_need_recomputation_.assign(num_variables, true);
in_good_jumps_.assign(num_variables, false);
@@ -313,69 +443,91 @@ int FeasibilityJumpSolver::UpdateConstraintWeights(bool linear_mode) {
const double kMaxWeight = 1e50;
const double kBumpFactor = 1.0 / params_.feasibility_jump_decay();
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
bump_value_ *= kBumpFactor;
bool rescale = false;
int num_bumps = 0;
tmp_constraints_.clear();
const int num_constraints = linear_mode
? evaluator_->NumLinearConstraints()
: evaluator_->NumEvaluatorConstraints();
linear_evaluator->ClearAffectedVariables();
for (int c = 0; c < num_constraints; ++c) {
if (evaluator_->Violation(c) > 0) {
++num_bumps;
weights_[c] += bump_value_;
if (weights_[c] > kMaxWeight) rescale = true;
tmp_constraints_.push_back(c);
if (!evaluator_->IsViolated(c)) continue;
++num_bumps;
weights_[c] += bump_value_;
if (weights_[c] > kMaxWeight) rescale = true;
if (!rescale) {
linear_evaluator->UpdateScoreOnWeightUpdate(
c, bump_value_, evaluator_->current_solution(), jump_values_,
absl::MakeSpan(jump_scores_));
}
}
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
if (rescale) {
const double factor = 1.0 / kMaxWeight;
bump_value_ *= factor;
for (int c = 0; c < num_constraints; ++c) {
weights_[c] *= factor;
}
if (linear_mode) {
solution_score_ = linear_evaluator->WeightedViolation(weights_);
RecomputeAllJumps();
return num_bumps;
}
RecomputeAllJumps();
return num_bumps;
}
// Recompute the affected jumps and overall score.
// Recompute the affected jumps.
// Note that the constraint violations are unaffected.
if (linear_mode) {
solution_score_ = linear_evaluator->WeightedViolation(weights_);
for (const int var :
linear_evaluator->ConstraintsToAffectedVariables(tmp_constraints_)) {
MarkJumpForRecomputation(var);
for (const int var : linear_evaluator->VariablesAffectedByLastUpdate()) {
// We don't need to recompute score of binary variable, it should
// already be correct.
if (!jump_need_recomputation_[var] && var_has_two_values_[var]) {
DCHECK(JumpIsUpToDate(var));
if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
continue;
}
// This jump might be good, so we need to add it to the queue so it can be
// evaluated when choosing the next jump.
jump_need_recomputation_[var] = true;
if (!in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
} else {
solution_score_ = evaluator_->WeightedViolation(weights_);
}
return num_bumps;
}
// Important: This is for debugging, but unfortunately it currently change the
// deterministic time and hence the overall algo behavior.
bool FeasibilityJumpSolver::JumpIsUpToDate(int var) {
const int64_t old_jump = jump_values_[var];
const double old_score = jump_scores_[var];
RecomputeJump(var);
CHECK_EQ(jump_values_[var], old_jump); // No change
const double relative =
std::max({std::abs(jump_scores_[var]), std::abs(old_score), 1.0});
return std::abs(jump_scores_[var] - old_score) / relative < 1e-6;
}
bool FeasibilityJumpSolver::DoSomeLinearIterations() {
++num_batches_;
// Recompute at the beginning of each batch.
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
solution_score_ = linear_evaluator->WeightedViolation(weights_);
const std::vector<int64_t>& solution = evaluator_->current_solution();
RecomputeAllJumps();
if (VLOG_IS_ON(1)) {
if (num_batches_ < 10) {
// The first few batches are more informative to understand how the
// heuristic behave.
// heuristic behaves.
shared_response_->LogMessage(name(), OneLineStats());
} else {
shared_response_->LogPeriodicMessage(name(), OneLineStats(),
@@ -403,7 +555,9 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() {
// Take the best jump score amongst some random candidates.
// It is okay if we pick twice the same, we don't really care.
int best_var = -1;
double best_delta = 0.0;
int best_index = -1;
int64_t best_value = 0;
double best_score = 0.0;
int num_improving_jump_tested = 0;
while (!good_jumps_.empty() && num_improving_jump_tested < 5) {
const int index = absl::Uniform<int>(
@@ -411,52 +565,78 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() {
const int var = good_jumps_[index];
// We lazily update the jump value.
if (jump_need_recomputation_[var]) RecomputeJump(var);
if (jump_need_recomputation_[var]) {
RecomputeJump(var);
} else {
DCHECK(JumpIsUpToDate(var));
}
if (!IsGood(jump_deltas_[var])) {
if (!IsGood(jump_scores_[var])) {
// Lazily remove.
in_good_jumps_[var] = false;
good_jumps_[index] = good_jumps_.back();
good_jumps_.pop_back();
if (best_index == good_jumps_.size()) best_index = index;
continue;
}
++num_improving_jump_tested;
if (jump_deltas_[var] <= best_delta) {
if (jump_scores_[var] <= best_score) {
best_var = var;
best_delta = jump_deltas_[var];
best_index = index;
best_value = jump_values_[var];
best_score = jump_scores_[var];
}
}
if (good_jumps_.empty()) {
break;
// TODO(user): Re-enable the code. It can be a bit slow currently
// especially since in many situation it doesn't achieve anything.
if (true) break;
bool time_limit_crossed = false;
if (ScanAllVariables(solution, /*linear_mode=*/true, &best_var,
&best_value, &best_score, &time_limit_crossed)) {
VLOG(2) << name()
<< " scans and finds improving solution (var:" << best_var
<< " value:" << best_value << ", delta:" << best_score << ")";
++num_repairs_with_full_scan_;
} else if (time_limit_crossed) {
return false;
} else {
break;
}
} else {
const int64_t var_change = best_value - solution[best_var];
DCHECK_EQ(best_score, linear_evaluator->WeightedViolationDelta(
weights_, best_var, var_change));
}
const int64_t var_change = jump_values_[best_var] - solution[best_var];
DCHECK_EQ(best_delta, linear_evaluator->WeightedViolationDelta(
weights_, best_var, var_change));
CHECK_NE(best_var, -1);
CHECK_NE(best_index, -1);
// Perform the move.
++num_jumps_;
CHECK_NE(best_var, -1);
evaluator_->UpdateVariableValueAndRecomputeViolations(
best_var, jump_values_[best_var], /*focus_on_linear=*/true);
solution_score_ += best_delta;
DCHECK_EQ(solution_score_, linear_evaluator->WeightedViolation(weights_));
++num_linear_moves_;
const int64_t old_value = evaluator_->current_solution()[best_var];
linear_evaluator->ClearAffectedVariables();
evaluator_->UpdateLinearScores(best_var, best_value, weights_,
jump_values_,
absl::MakeSpan(jump_scores_));
evaluator_->UpdateVariableValue(best_var, best_value);
// For now we recompute all the affected best jump value.
//
// TODO(user): In the paper, they just recompute the scores and only
// change the jump values when the constraint weight change. This should
// be faster.
//
// TODO(user): We can reduce the complexity if we just recompute the
// score. we don't need to rescan each column of the affected variables
// again, but could update their scores directly.
for (const int var :
linear_evaluator->AffectedVariableOnUpdate(best_var)) {
MarkJumpForRecomputation(var);
// We already know the score of undoing the move we just did, and we know
// this move is bad, so we can remove it from good_jumps_ right away.
jump_values_[best_var] = old_value;
jump_scores_[best_var] = -best_score;
if (var_has_two_values_[best_var]) {
CHECK_EQ(good_jumps_[best_index], best_var);
in_good_jumps_[best_var] = false;
good_jumps_[best_index] = good_jumps_.back();
good_jumps_.pop_back();
} else {
jump_need_recomputation_[best_var] = true;
}
MarkJumpsThatNeedsToBeRecomputed(best_var);
}
// We will update the weight unless the queue is non-empty.
@@ -468,52 +648,90 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() {
return false;
}
// Update the jump scores.
//
// We incrementally maintain the score (except for best_var).
// However for non-Boolean, we still need to recompute the jump value.
// We will do that in a lazy fashion.
//
// TODO(user): In the paper, they just recompute the scores and only
// change the jump values when the constraint weight changes. Experiment?
// Note however that the current code is quite fast.
//
// TODO(user): For non-Boolean, we could easily detect if a non-improving
// score cannot become improving. We don't need to add such variable to
// the queue.
void FeasibilityJumpSolver::MarkJumpsThatNeedsToBeRecomputed(int changed_var) {
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
for (const int var : linear_evaluator->VariablesAffectedByLastUpdate()) {
if (var == changed_var) continue;
if (jump_need_recomputation_[var]) {
DCHECK(in_good_jumps_[var]);
continue;
}
// We don't need to recompute score of binary variable, it should
// already be correct.
if (var_has_two_values_[var]) {
DCHECK(JumpIsUpToDate(var));
if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
continue;
}
jump_need_recomputation_[var] = true;
if (!in_good_jumps_[var]) {
in_good_jumps_[var] = true;
good_jumps_.push_back(var);
}
}
}
bool FeasibilityJumpSolver::DoSomeGeneralIterations() {
if (evaluator_->NumNonLinearConstraints() == 0) {
return true;
}
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
const std::vector<int64_t>& solution = evaluator_->current_solution();
// Non linear constraints are not evaluated in the linear phase.
evaluator_->UpdateNonLinearViolations();
solution_score_ = evaluator_->WeightedViolation(weights_);
// Non-linear constraints are not evaluated in the linear phase.
evaluator_->UpdateAllNonLinearViolations();
const int kLoops = 10000;
std::vector<int> to_scan;
for (int loop = 0; loop < kLoops; ++loop) {
to_scan = evaluator_->VariablesInViolatedConstraints();
std::shuffle(to_scan.begin(), to_scan.end(), random_);
bool improving_move_found = false;
int var;
int64_t value;
double score;
bool time_limit_crossed = false;
if (ScanAllVariables(solution, /*linear_mode=*/false, &var, &value, &score,
&time_limit_crossed)) {
// Perform the move.
num_general_moves_++;
for (const int current_var : to_scan) {
const int64_t initial_value = solution[current_var];
int64_t best_value = initial_value;
double best_delta = 0.0;
for (const int64_t new_value : var_domains_[current_var].Values()) {
if (new_value == initial_value) continue;
const double delta = evaluator_->WeightedViolationDelta(
weights_, current_var, new_value - initial_value);
if (evaluator_->num_deltas_computed() % 1000 == 0 &&
shared_time_limit_ != nullptr &&
shared_time_limit_->LimitReached()) {
return false;
}
if (delta < best_delta) {
improving_move_found = true;
best_value = new_value;
best_delta = delta;
// Update the linear part.
linear_evaluator->ClearAffectedVariables();
evaluator_->UpdateLinearScores(var, value, weights_, jump_values_,
absl::MakeSpan(jump_scores_));
jump_values_[var] = solution[var];
jump_scores_[var] = -score;
for (const int v : linear_evaluator->VariablesAffectedByLastUpdate()) {
if (!var_has_two_values_[var]) {
jump_need_recomputation_[v] = true;
}
}
if (improving_move_found) { // Accept the move with the best delta.
solution_score_ += best_delta;
evaluator_->UpdateVariableValueAndRecomputeViolations(current_var,
best_value);
break;
}
// Update the non-linear part. Note it also commits the move.
evaluator_->UpdateNonLinearViolations(var, value);
evaluator_->UpdateVariableValue(var, value);
continue;
} else if (time_limit_crossed) {
return false;
}
if (UpdateConstraintWeights(/*linear_mode=*/false) == 0) {
return true;
}
@@ -521,4 +739,142 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() {
return false;
}
bool FeasibilityJumpSolver::ScanAllVariables(absl::Span<const int64_t> solution,
bool linear_mode,
int* improving_var,
int64_t* improving_value,
double* improving_delta,
bool* time_limit_crossed) {
LinearIncrementalEvaluator* linear_evaluator =
evaluator_->MutableLinearEvaluator();
// We stop at the first improving move that does not degrade the linear score.
// Otherwise, we return the best improving move.
tmp_to_scan_ = evaluator_->VariablesInViolatedConstraints();
std::shuffle(tmp_to_scan_.begin(), tmp_to_scan_.end(), random_);
// We maintain the best solution found, regardless of the linear violations.
int best_var = -1;
int64_t best_value = 0;
double best_delta = -1e-4;
// We will favor moves that keep the linear assignment feasible for as long as
// possible (if feasible to starts with).
//
// TODO(user): We just need to check if the assignment is linear
// feasible, not compute its score.
const bool accept_moves_with_linear_violation =
!params_.feasibility_jump_protect_linear_feasibility() ||
linear_evaluator->WeightedViolation(weights_) > 0.0;
for (const int var : tmp_to_scan_) {
if (evaluator_->VariableOnlyInLinearConstraintWithConvexViolationChange(
var)) {
// We lazily update the jump value.
if (jump_need_recomputation_[var]) {
RecomputeJump(var);
} else {
DCHECK(JumpIsUpToDate(var));
}
const double delta = jump_scores_[var];
if (IsGood(delta)) {
*improving_var = var;
*improving_value = jump_values_[var];
*improving_delta = delta;
return true;
}
continue;
}
// Compute the values to scan.
// If the domain is large we do not scan the full variable domain.
Domain scan_range;
const int64_t initial_value = solution[var];
const Domain& var_domain = var_domains_[var];
if (var_domains_[var].Size() <=
params_.feasibility_jump_max_num_values_scanned()) {
scan_range = var_domain;
} else {
// TODO(user): Partial scan of the rest of the domain.
++num_partial_scans_;
const int64_t half_span =
params_.feasibility_jump_max_num_values_scanned() / 2;
Domain scan_range =
Domain(initial_value).AdditionWith({-half_span, half_span});
if (scan_range.Min() < var_domain.Min()) {
scan_range = scan_range.AdditionWith(
Domain(var_domain.Min() - scan_range.Min()));
} else if (scan_range.Max() > var_domain.Max()) {
scan_range = scan_range.AdditionWith(
Domain(scan_range.Max() - var_domain.Max()));
}
scan_range = scan_range.IntersectionWith(var_domain);
}
// Take the best amongst all the values to scan for the current variable if
// it does not degrade the linear score (when
// accept_moves_with_linear_violation is false).
int best_var_non_increasing_linear = -1;
int64_t best_value_non_increasing_linear = 0;
double best_delta_non_increasing_linear = -1e-4;
for (const int64_t new_value : scan_range.Values()) {
if (new_value == initial_value) continue;
// Checks the time limit periodically.
if (++num_general_evals_ % 1000 == 0 && shared_time_limit_ != nullptr &&
shared_time_limit_->LimitReached()) {
*time_limit_crossed = true;
return false;
}
// Computes both linear and general delta.
const double linear_delta = linear_evaluator->WeightedViolationDelta(
weights_, var, new_value - initial_value);
const double non_linear_delta =
linear_mode ? 0.0
: evaluator_->WeightedNonLinearViolationDelta(
weights_, var, new_value - initial_value);
const double delta = linear_delta + non_linear_delta;
// Do we have a valid improving move ? We select the best value for the
// first variable found with such a move so we delay the exit until after
// the loop on values is finished.
if ((accept_moves_with_linear_violation || linear_delta <= 0.0) &&
delta < best_delta_non_increasing_linear) {
best_var_non_increasing_linear = var;
best_value_non_increasing_linear = new_value;
best_delta_non_increasing_linear = delta;
}
// We keep the best move so far.
if (delta < best_delta) {
best_var = var;
best_value = new_value;
best_delta = delta;
}
} // End scan values.
// Early accept the best improving move that does not degrade the linear
// score, or the best move if accept_moves_with_linear_violation is true.
if (best_var_non_increasing_linear != -1) {
*improving_var = best_var_non_increasing_linear;
*improving_value = best_value_non_increasing_linear;
*improving_delta = best_delta_non_increasing_linear;
return true;
}
}
// Accept the best improving move.
if (best_var != -1) {
*improving_var = best_var;
*improving_value = best_value;
*improving_delta = best_delta;
return true;
}
// No solution found.
return false;
}
} // namespace operations_research::sat

53
ortools/sat/feasibility_jump.h
View File

@@ -17,6 +17,7 @@
#include <atomic>
#include <cstdint>
#include <functional>
#include <limits>
#include <memory>
#include <string>
#include <utility>
@@ -24,9 +25,9 @@
#include "absl/time/time.h"
#include "absl/types/span.h"
#include "ortools/base/vlog_is_on.h"
#include "ortools/sat/constraint_violation.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/sat/subsolver.h"
#include "ortools/sat/synchronization.h"
@@ -64,14 +65,8 @@ class FeasibilityJumpSolver : public SubSolver {
shared_stats_(shared_stats),
random_(params_) {}
~FeasibilityJumpSolver() override {
if (!VLOG_IS_ON(1)) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back({"fs_jump/num_jumps", num_jumps_});
stats.push_back({"fs_jump/num_computed_jumps", num_computed_jumps_});
stats.push_back({"fs_jump/num_updates", num_weight_updates_});
shared_stats_->AddStats(stats);
}
// If VLOG_IS_ON(1), it will export a bunch of statistics.
~FeasibilityJumpSolver() override;
// No synchronization needed for TaskIsAvailable().
void Synchronize() final {}
@@ -105,14 +100,20 @@ class FeasibilityJumpSolver : public SubSolver {
std::string OneLineStats() const;
// Linear only.
bool JumpIsUpToDate(int var); // Debug.
void RecomputeJump(int var);
void MarkJumpForRecomputation(int var);
void RecomputeAllJumps();
void MarkJumpsThatNeedsToBeRecomputed(int changed_var);
// Moves.
bool DoSomeLinearIterations();
bool DoSomeGeneralIterations();
// Returns true if an improving move was found.
bool ScanAllVariables(absl::Span<const int64_t> solution, bool linear_mode,
int* improving_var, int64_t* improving_value,
double* improving_delta, bool* time_limit_crossed);
// Return the number of infeasible constraints.
int UpdateConstraintWeights(bool linear_mode);
@@ -134,6 +135,7 @@ class FeasibilityJumpSolver : public SubSolver {
std::unique_ptr<LsEvaluator> evaluator_;
std::vector<Domain> var_domains_;
std::vector<bool> var_has_two_values_;
// For each variable, we store:
// - A jump value, which is different from the current solution, except for
@@ -142,37 +144,46 @@ class FeasibilityJumpSolver : public SubSolver {
// jump.
std::vector<bool> jump_need_recomputation_;
std::vector<int64_t> jump_values_;
std::vector<double> jump_deltas_;
std::vector<double> jump_scores_;
// The score of a solution is just the sum of infeasibility of each
// constraints weighted by these scores.
// constraint weighted by these scores.
std::vector<double> weights_;
// The current weighted violation (lower is better).
double solution_score_;
// Depending on the options, we use an exponentially decaying constraint
// weight like for SAT activities.
double bump_value_ = 1.0;
// Sparse list of jump with a potential delta <= 0.0.
// We lazily recompute the exact delta as we randomly pick variable from here.
// Sparse list of jumps with a potential delta < 0.0.
// If jump_need_recomputation_[var] is true, we lazily recompute the exact
// delta as we randomly pick variables from here.
std::vector<bool> in_good_jumps_;
std::vector<int> good_jumps_;
// We restart each time our local deterministic time crosses this.
double dtime_restart_threshold_ = 0.0;
int64_t update_restart_threshold_ = 0;
std::vector<int> tmp_constraints_;
std::vector<int64_t> tmp_breakpoints_;
// Statistics
absl::Time last_logging_time_;
int64_t num_restarts_ = 0;
int64_t num_batches_ = 0;
int64_t num_jumps_ = 0;
int64_t num_computed_jumps_ = 0;
int64_t num_linear_evals_ = 0;
int64_t num_general_evals_ = 0;
int64_t num_general_moves_ = 0;
int64_t num_linear_moves_ = 0;
int64_t num_partial_scans_ = 0;
int64_t num_repairs_with_full_scan_ = 0;
int64_t num_restarts_ = 0;
int64_t num_solutions_imported_ = 0;
int64_t num_weight_updates_ = 0;
// Temporary storage.
std::vector<int> tmp_to_scan_;
// Info on the last solution loaded.
int64_t last_solution_rank_ = std::numeric_limits<int64_t>::max();
};
} // namespace operations_research::sat

43
ortools/sat/linear_constraint_manager.cc
View File

@@ -107,7 +107,7 @@ void LinearConstraintManager::RescaleActiveCounts(const double scaling_factor) {
constraint_infos_[i].active_count *= scaling_factor;
}
constraint_active_count_increase_ *= scaling_factor;
VLOG(2) << "Rescaled active counts by " << scaling_factor;
VLOG(3) << "Rescaled active counts by " << scaling_factor;
}
bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
@@ -148,7 +148,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
lp_constraints_.resize(new_size);
solution_state->statuses.resize(num_cols + glop::ColIndex(new_size));
if (num_removed_constraints > 0) {
VLOG(2) << "Removed " << num_removed_constraints << " constraints";
VLOG(3) << "Removed " << num_removed_constraints << " constraints";
}
return num_removed_constraints > 0;
}
@@ -171,18 +171,22 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
const ConstraintIndex ct_index = equiv_constraints_[key];
if (constraint_infos_[ct_index].constraint.vars == ct.vars &&
constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) {
if (added != nullptr) *added = false;
bool tightened = false;
if (ct.lb > constraint_infos_[ct_index].constraint.lb) {
tightened = true;
if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
constraint_infos_[ct_index].constraint.lb = ct.lb;
if (added != nullptr) *added = true;
}
if (ct.ub < constraint_infos_[ct_index].constraint.ub) {
tightened = true;
if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
constraint_infos_[ct_index].constraint.ub = ct.ub;
if (added != nullptr) *added = true;
}
++num_merged_constraints_;
if (added != nullptr) *added = tightened;
if (tightened) {
++num_merged_constraints_;
FillActivityRange(&constraint_infos_[ct_index]);
}
return ct_index;
}
}
@@ -192,6 +196,7 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
ConstraintInfo ct_info;
ct_info.constraint = std::move(ct);
ct_info.l2_norm = ComputeL2Norm(ct_info.constraint);
FillActivityRange(&ct_info);
ct_info.hash = key;
equiv_constraints_[key] = ct_index;
ct_info.active_count = constraint_active_count_increase_;
@@ -246,7 +251,7 @@ bool LinearConstraintManager::AddCut(
// Only add cut with sufficient efficacy.
if (violation / l2_norm < 1e-5) {
VLOG(2) << "BAD Cut '" << type_name << "'"
VLOG(3) << "BAD Cut '" << type_name << "'"
<< " size=" << ct.vars.size()
<< " max_magnitude=" << ComputeInfinityNorm(ct)
<< " norm=" << l2_norm << " violation=" << violation
@@ -328,7 +333,7 @@ void LinearConstraintManager::PermanentlyRemoveSomeConstraints() {
}
if (num_deleted_constraints > 0) {
VLOG(2) << "Constraint manager cleanup: #deleted:"
VLOG(3) << "Constraint manager cleanup: #deleted:"
<< num_deleted_constraints;
}
num_deletable_constraints_ -= num_deleted_constraints;
@@ -509,6 +514,28 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
return term_changed;
}
void LinearConstraintManager::FillActivityRange(ConstraintInfo* info) {
IntegerValue min_sum(0);
IntegerValue max_sum(0);
const int num_terms = info->constraint.vars.size();
for (int i = 0; i < num_terms; ++i) {
const IntegerVariable var = info->constraint.vars[i];
const IntegerValue coeff = info->constraint.coeffs[i];
const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
if (coeff > 0.0) {
min_sum += coeff * lb;
max_sum += coeff * ub;
} else {
min_sum += coeff * ub;
max_sum += coeff * lb;
}
}
const IntegerValue tight_lb = std::max(min_sum, info->constraint.lb);
const IntegerValue tight_ub = std::min(max_sum, info->constraint.ub);
info->activity_range = CapSub(tight_ub.value(), tight_lb.value());
}
bool LinearConstraintManager::ChangeLp(
const absl::StrongVector<IntegerVariable, double>& lp_solution,
glop::BasisState* solution_state, int* num_new_constraints) {

10
ortools/sat/linear_constraint_manager.h
View File

@@ -52,6 +52,12 @@ class LinearConstraintManager {
public:
struct ConstraintInfo {
LinearConstraint constraint;
// For a boxed constraint, this is the constraint rhs - lhs. For a one sided
// constraint, we first use trivial activity bound to compute tight rhs and
// lhs and then use the same formula.
IntegerValue activity_range;
double l2_norm = 0.0;
int64_t inactive_count = 0;
double objective_parallelism = 0.0;
@@ -169,7 +175,7 @@ class LinearConstraintManager {
// Helper method to compute objective parallelism for a given constraint. This
// also lazily computes objective norm.
void ComputeObjectiveParallelism(const ConstraintIndex ct_index);
void ComputeObjectiveParallelism(ConstraintIndex ct_index);
// Multiplies all active counts and the increment counter by the given
// 'scaling_factor'. This should be called when at least one of the active
@@ -180,6 +186,8 @@ class LinearConstraintManager {
// don't remove any constraints which are already in LP.
void PermanentlyRemoveSomeConstraints();
void FillActivityRange(ConstraintInfo* info);
const SatParameters& sat_parameters_;
const IntegerTrail& integer_trail_;

29
ortools/sat/linear_programming_constraint.cc
View File

@@ -292,6 +292,7 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
// Fill integer_lp_.
integer_lp_.clear();
infinity_norms_.clear();
ct_bound_diff_.clear();
const auto& all_constraints = constraint_manager_.AllConstraints();
for (const auto index : constraint_manager_.LpConstraints()) {
const LinearConstraint& ct = all_constraints[index].constraint;
@@ -321,8 +322,9 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
new_ct.terms.push_back({GetMirrorVariable(var), coeff});
}
infinity_norms_.push_back(infinity_norm);
ct_bound_diff_.push_back(all_constraints[index].activity_range);
// Important to keep lp_data_ "clean".
// It is important to keep lp_data_ "clean".
DCHECK(std::is_sorted(new_ct.terms.begin(), new_ct.terms.end()));
}
@@ -889,8 +891,10 @@ bool LinearProgrammingConstraint::PreprocessCut(LinearConstraint* cut) {
const IntegerValue slack{CapSub(cut->ub.value(), min_sum.value())};
if (!min_sum_overflow && !AtMinOrMaxInt64(slack.value())) {
// TODO(user): raise conflict or report UNSAT.
if (slack < 0) return false; // Always false.
if (slack < 0) {
problem_proven_infeasible_by_cuts_ = true;
return false;
}
if (trail_->CurrentDecisionLevel() == 0 && max_range > slack) {
bool newly_fixed = false;
@@ -936,6 +940,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
IntegerValue cut_ub;
if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
&cut_ub)) {
++num_cut_overflows_;
VLOG(1) << "Issue, overflow!";
return false;
}
@@ -988,6 +993,8 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
// TODO(user): Keep track of the potential overflow here.
if (!base_ct_.FillFromLinearConstraint(cut_, expanded_lp_solution_,
integer_trail_)) {
++num_cut_overflows_;
VLOG(1) << "Issue, overflow!";
return false;
}
@@ -1019,8 +1026,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
CutTerm entry;
entry.coeff = IntTypeAbs(coeff);
entry.lp_value = 0.0;
entry.bound_diff =
CapSub(integer_lp_[row].ub.value(), integer_lp_[row].lb.value());
entry.bound_diff = ct_bound_diff_[row];
entry.expr_vars[0] =
first_slack + 2 * IntegerVariable(tmp_slack_rows_.size());
entry.expr_coeffs[1] = 0;
@@ -1650,10 +1656,14 @@ bool LinearProgrammingConstraint::Propagate() {
// TODO(user): Refactor so that they are just normal cut generators?
const int level = trail_->CurrentDecisionLevel();
if (trail_->CurrentDecisionLevel() == 0) {
problem_proven_infeasible_by_cuts_ = false;
if (parameters_.add_objective_cut()) AddObjectiveCut();
if (parameters_.add_mir_cuts()) AddMirCuts();
if (parameters_.add_cg_cuts()) AddCGCuts();
if (parameters_.add_zero_half_cuts()) AddZeroHalfCuts();
if (problem_proven_infeasible_by_cuts_) {
return integer_trail_->ReportConflict({});
}
}
// Try to add cuts.
@@ -2316,6 +2326,11 @@ bool LinearProgrammingConstraint::ExactLpReasonning() {
tmp_lp_multipliers_.push_back({row, value});
}
// In this case, the LP lower bound match the basic objective "constraint"
// propagation. That is there is an LP solution with all objective variable at
// their current best bound. There is no need to do more work here.
if (tmp_lp_multipliers_.empty()) return true;
IntegerValue scaling = 0;
tmp_integer_multipliers_ = ScaleLpMultiplier(
/*take_objective_into_account=*/true, tmp_lp_multipliers_, &scaling);
@@ -2613,8 +2628,8 @@ std::string LinearProgrammingConstraint::Statistics() const {
FormatCounter(total_num_simplex_iterations_), "\n");
absl::StrAppend(&result, " total num cut propagation: ",
FormatCounter(total_num_cut_propagations_), "\n");
absl::StrAppend(&result, " total num prevent overflow: ",
FormatCounter(num_prevent_overflows_), "\n");
absl::StrAppend(&result, " total num cut overflow: ",
FormatCounter(num_cut_overflows_), "\n");
absl::StrAppend(&result, " total num adjust: ", FormatCounter(num_adjusts_),
"\n");
absl::StrAppend(&result, " total num scaling issues: ",

4
ortools/sat/linear_programming_constraint.h
View File

@@ -438,6 +438,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
IntegerValue objective_infinity_norm_ = IntegerValue(0);
absl::StrongVector<glop::RowIndex, LinearConstraintInternal> integer_lp_;
absl::StrongVector<glop::RowIndex, IntegerValue> infinity_norms_;
absl::StrongVector<glop::RowIndex, IntegerValue> ct_bound_diff_;
// Underlying LP solver API.
glop::GlopParameters simplex_params_;
@@ -455,6 +456,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
FlowCoverCutHelper flow_cover_cut_helper_;
IntegerRoundingCutHelper integer_rounding_cut_helper_;
bool problem_proven_infeasible_by_cuts_ = false;
CutData base_ct_;
LinearConstraint cut_;
LinearConstraint saved_cut_;
@@ -580,7 +582,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
// Some stats on the LP statuses encountered.
int64_t num_solves_ = 0;
mutable int64_t num_adjusts_ = 0;
mutable int64_t num_prevent_overflows_ = 0;
mutable int64_t num_cut_overflows_ = 0;
mutable int64_t num_scaling_issues_ = 0;
std::vector<int64_t> num_solves_by_status_;
};

11
ortools/sat/parameters_validation.cc
View File

@@ -13,6 +13,8 @@
#include "ortools/sat/parameters_validation.h"
#include <cmath>
#include <limits>
#include <string>
#include <vector>
@@ -89,7 +91,6 @@ std::string ValidateParameters(const SatParameters& params) {
TEST_NOT_NAN(max_time_in_seconds);
TEST_NOT_NAN(max_deterministic_time);
TEST_NOT_NAN(feasibility_jump_decay);
// TODO(user): Consider using annotations directly in the proto for these
// validation. It is however not open sourced.
@@ -97,7 +98,15 @@ std::string ValidateParameters(const SatParameters& params) {
TEST_IN_RANGE(mip_max_bound, 0, 1e17);
TEST_IN_RANGE(solution_pool_size, 1, std::numeric_limits<int32_t>::max());
TEST_IN_RANGE(shared_tree_worker_objective_split_probability, 0.0, 1.0);
// Feasibility jump.
TEST_NOT_NAN(feasibility_jump_decay);
TEST_NOT_NAN(feasibility_jump_var_randomization_ratio);
TEST_NOT_NAN(feasibility_jump_var_perburbation_range_ratio);
TEST_IN_RANGE(feasibility_jump_max_num_values_scanned, 2, 1'000'000'000);
TEST_IN_RANGE(feasibility_jump_decay, 0.0, 1.0);
TEST_IN_RANGE(feasibility_jump_var_randomization_ratio, 0.0, 1.0);
TEST_IN_RANGE(feasibility_jump_var_perburbation_range_ratio, 0.0, 1.0);
TEST_POSITIVE(glucose_decay_increment_period);
TEST_POSITIVE(shared_tree_max_nodes_per_worker);

12
ortools/sat/presolve_context.cc
View File

@@ -1916,6 +1916,18 @@ bool PresolveContext::ShiftCostInExactlyOne(absl::Span<const int> exactly_one,
// Note that the domain never include the offset, so we need to update it.
if (offset != 0) AddToObjectiveOffset(offset);
// When we shift the cost using an exactly one, our objective implied bounds
// might be more or less precise. If the objective domain is not constraining
// (and thus just constraining the upper bound), we relax it to make sure its
// stay "non constraining".
//
// TODO(user): This is a bit hacky, find a nicer way.
if (!objective_domain_is_constraining_) {
objective_domain_ =
Domain(std::numeric_limits<int64_t>::min(), objective_domain_.Max());
}
return true;
}

5
ortools/sat/presolve_context.h
View File

@@ -24,6 +24,7 @@
#include "absl/base/attributes.h"
#include "absl/container/flat_hash_map.h"
#include "absl/container/flat_hash_set.h"
#include "absl/log/check.h"
#include "absl/strings/str_cat.h"
#include "absl/types/span.h"
#include "ortools/base/logging.h"
@@ -55,7 +56,7 @@ class PresolveContext;
// bunch of places.
class SavedLiteral {
public:
SavedLiteral() {}
SavedLiteral() = default;
explicit SavedLiteral(int ref) : ref_(ref) {}
int Get(PresolveContext* context) const;
@@ -70,7 +71,7 @@ class SavedLiteral {
// general affine for the linear1 involving an absolute value.
class SavedVariable {
public:
SavedVariable() {}
SavedVariable() = default;
explicit SavedVariable(int ref) : ref_(ref) {}
int Get() const;

30
ortools/sat/sat_parameters.proto
View File

@@ -23,7 +23,7 @@ option csharp_namespace = "Google.OrTools.Sat";
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 245
// NEXT TAG: 249
message SatParameters {
// In some context, like in a portfolio of search, it makes sense to name a
// given parameters set for logging purpose.
@@ -1069,12 +1069,32 @@ message SatParameters {
//
// The test_feasibility_jump is used to only enable this for benchmarking.
optional bool test_feasibility_jump = 240 [default = false];
// Max number of values scanned when scanning the domain of a variable.
// This number is checked to be >= 2 and <= 1e9.
optional int64 feasibility_jump_max_num_values_scanned = 245 [default = 4096];
// Feasibility jump focuses on the linear model, then switches to the full
// model when a linear feasible assignment is found. In the second phase, this
// parameters favors moves that keep the linear assignment feasible over move
// that will break it.
optional bool feasibility_jump_protect_linear_feasibility = 246
[default = true];
// Use feasibility_jump to find improving solutions.
optional bool use_violation_ls = 244 [default = false];
// The following feasibility_jump parameters meant to be used by the solver on
// the different workers that use feasibility_jump.
optional double feasibility_jump_decay = 242 [default = 1.0];
optional bool feasibility_jump_start_with_objective = 243 [default = false];
// Use a local search worker based on constraint violations to find improving
// solutions.
optional bool use_violation_ls = 244 [default = false];
// Ratio of variables that will not have their default value upon restart.
optional double feasibility_jump_var_randomization_ratio = 247
[default = 0.0];
// Max distance between the default value and the pertubated value relative to
// the range of the domain of the variable.
optional double feasibility_jump_var_perburbation_range_ratio = 248
[default = 0.2];
// Enables experimental workstealing-like shared tree search.
// If non-zero, start this many complete worker threads to explore a shared

21
ortools/sat/synchronization.cc
View File

@@ -35,13 +35,14 @@
#include "absl/container/flat_hash_map.h"
#include "absl/container/flat_hash_set.h"
#include "absl/flags/flag.h"
#include "absl/status/status.h"
#include "absl/log/check.h"
#include "absl/strings/str_cat.h"
#include "absl/strings/str_format.h"
#include "absl/strings/string_view.h"
#include "absl/synchronization/mutex.h"
#include "absl/time/clock.h"
#include "absl/time/time.h"
#include "absl/types/span.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/integer.h"
@@ -54,7 +55,6 @@
#include "ortools/util/logging.h"
#include "ortools/util/sorted_interval_list.h"
#include "ortools/util/strong_integers.h"
#include "ortools/util/time_limit.h"
ABSL_FLAG(bool, cp_model_dump_solutions, false,
"DEBUG ONLY. If true, all the intermediate solution will be dumped "
@@ -665,7 +665,7 @@ void SharedResponseManager::NewSolution(
if (obj.scaling_factor() < 0) {
std::swap(lb, ub);
}
RegisterSolutionFound(solution_message);
RegisterSolutionFound(solution_message, num_solutions_);
SOLVER_LOG(logger_, ProgressMessage(absl::StrCat(num_solutions_),
wall_timer_.Get(), best, lb, ub,
solution_message));
@@ -733,9 +733,12 @@ std::string ExtractSubSolverName(const std::string& improvement_info) {
}
void SharedResponseManager::RegisterSolutionFound(
const std::string& improvement_info) {
const std::string& improvement_info, int solution_rank) {
if (improvement_info.empty()) return;
primal_improvements_count_[ExtractSubSolverName(improvement_info)]++;
const std::string subsolver_name = ExtractSubSolverName(improvement_info);
primal_improvements_count_[subsolver_name]++;
primal_improvements_min_rank_.insert({subsolver_name, solution_rank});
primal_improvements_max_rank_[subsolver_name] = solution_rank;
}
void SharedResponseManager::RegisterObjectiveBoundImprovement(
@@ -748,9 +751,13 @@ void SharedResponseManager::DisplayImprovementStatistics() {
absl::MutexLock mutex_lock(&mutex_);
if (!primal_improvements_count_.empty()) {
SOLVER_LOG(logger_, "");
SOLVER_LOG(logger_, "Solutions found per subsolver:");
SOLVER_LOG(logger_, "Solutions found per subsolver (", num_solutions_,
"):");
for (const auto& entry : primal_improvements_count_) {
SOLVER_LOG(logger_, " '", entry.first, "': ", entry.second);
const int min_rank = primal_improvements_min_rank_[entry.first];
const int max_rank = primal_improvements_max_rank_[entry.first];
SOLVER_LOG(logger_, " '", entry.first, "':", entry.second, " rank:[",
min_rank, ",", max_rank, "]");
}
}
if (!dual_improvements_count_.empty()) {

12
ortools/sat/synchronization.h
View File

@@ -31,18 +31,18 @@
#include "absl/random/random.h"
#include "absl/synchronization/mutex.h"
#include "absl/time/time.h"
#include "ortools/base/integral_types.h"
#include "absl/types/span.h"
#include "ortools/base/logging.h"
#include "ortools/base/stl_util.h"
#include "ortools/base/timer.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/model.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/sat/util.h"
#include "ortools/util/bitset.h"
#include "ortools/util/logging.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research {
namespace sat {
@@ -358,7 +358,8 @@ class SharedResponseManager {
ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_);
void UpdateGapIntegralInternal() ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_);
void RegisterSolutionFound(const std::string& improvement_info)
void RegisterSolutionFound(const std::string& improvement_info,
int solution_rank)
ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_);
void RegisterObjectiveBoundImprovement(const std::string& improvement_info)
ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_);
@@ -423,6 +424,11 @@ class SharedResponseManager {
// Used for statistics of the improvements found by workers.
absl::btree_map<std::string, int> primal_improvements_count_
ABSL_GUARDED_BY(mutex_);
absl::btree_map<std::string, int> primal_improvements_min_rank_
ABSL_GUARDED_BY(mutex_);
absl::btree_map<std::string, int> primal_improvements_max_rank_
ABSL_GUARDED_BY(mutex_);
absl::btree_map<std::string, int> dual_improvements_count_
ABSL_GUARDED_BY(mutex_);

17
ortools/sat/work_assignment.cc
View File

@@ -187,10 +187,14 @@ absl::Span<const ProtoLiteral> ProtoTrail::Implications(int level) const {
SharedTreeManager::SharedTreeManager(Model* model)
: params_(*model->GetOrCreate<SatParameters>()),
num_workers_(params_.shared_tree_num_workers()),
num_workers_(std::max(1, params_.shared_tree_num_workers())),
shared_response_manager_(model->GetOrCreate<SharedResponseManager>()),
num_splits_wanted_(num_workers_ - 1),
max_nodes_(num_workers_ * params_.shared_tree_max_nodes_per_worker()) {
max_nodes_(params_.shared_tree_max_nodes_per_worker() >=
std::numeric_limits<int>::max() / num_workers_
? std::numeric_limits<int>::max()
: num_workers_ *
params_.shared_tree_max_nodes_per_worker()) {
// Create the root node with a fake literal.
nodes_.push_back(
{.literal = ProtoLiteral(),
@@ -597,7 +601,10 @@ LiteralIndex SharedTreeWorker::NextDecision() {
CHECK(!sat_solver_->Assignment().LiteralIsTrue(decision));
return decision.Index();
}
if (objective_ == nullptr) return helper_->GetDecision(decision_policy);
if (objective_ == nullptr ||
objective_->objective_var == kNoIntegerVariable) {
return helper_->GetDecision(decision_policy);
}
// If the current node is close to the global lower bound, maybe try to
// improve it.
const IntegerValue root_obj_lb =
@@ -631,6 +638,8 @@ SatSolver::Status SharedTreeWorker::Search(
sat_solver_->Backtrack(0);
encoder_->GetTrueLiteral();
encoder_->GetFalseLiteral();
const bool has_objective =
objective_ != nullptr && objective_->objective_var != kNoIntegerVariable;
std::vector<Literal> clause;
while (!time_limit_->LimitReached()) {
if (!sat_solver_->FinishPropagation()) {
@@ -667,7 +676,7 @@ SatSolver::Status SharedTreeWorker::Search(
if (time_limit_->LimitReached()) return SatSolver::LIMIT_REACHED;
if (decision.Index() == kNoLiteralIndex) {
feasible_solution_observer();
if (objective_ == nullptr) return SatSolver::FEASIBLE;
if (!has_objective) return SatSolver::FEASIBLE;
const IntegerValue objective =
integer_trail_->LowerBound(objective_->objective_var);
sat_solver_->Backtrack(0);