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

[CP-SAT] cleaning

Browse Source
This commit is contained in:
Laurent Perron
2022-10-17 17:11:59 +02:00
parent 1f164099e6
commit 8948fcc24b
21 changed files with 570 additions and 452 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
ortools/sat/BUILD.bazel
View File

@@ -814,6 +814,7 @@ cc_library(
":linear_constraint",
":model",
":sat_base",
":synchronization",
"//ortools/base",
"//ortools/base:strong_vector",
"//ortools/util:bitset",
@@ -1031,7 +1032,7 @@ cc_library(
":linear_constraint",
":linear_programming_constraint",
":model",
":presolve_util",
":presolve_util",
":routing_cuts",
":sat_base",
":sat_parameters_cc_proto",

8
ortools/sat/cp_model_checker.cc
View File

@@ -1075,8 +1075,8 @@ namespace {
class ConstraintChecker {
public:
explicit ConstraintChecker(const std::vector<int64_t>& variable_values)
: variable_values_(variable_values) {}
explicit ConstraintChecker(absl::Span<const int64_t> variable_values)
: variable_values_(variable_values.begin(), variable_values.end()) {}
bool LiteralIsTrue(int l) const {
if (l >= 0) return variable_values_[l] != 0;
@@ -1524,13 +1524,13 @@ class ConstraintChecker {
}
private:
std::vector<int64_t> variable_values_;
const std::vector<int64_t> variable_values_;
};
} // namespace
bool SolutionIsFeasible(const CpModelProto& model,
const std::vector<int64_t>& variable_values,
absl::Span<const int64_t> variable_values,
const CpModelProto* mapping_proto,
const std::vector<int>* postsolve_mapping) {
if (variable_values.size() != model.variables_size()) {

2
ortools/sat/cp_model_checker.h
View File

@@ -63,7 +63,7 @@ bool PossibleIntegerOverflow(const CpModelProto& model,
// The last two arguments are optional and help debugging a failing constraint
// due to presolve.
bool SolutionIsFeasible(const CpModelProto& model,
const std::vector<int64_t>& variable_values,
absl::Span<const int64_t> variable_values,
const CpModelProto* mapping_proto = nullptr,
const std::vector<int>* postsolve_mapping = nullptr);

169
ortools/sat/cp_model_lns.cc
View File

@@ -1161,12 +1161,8 @@ void NeighborhoodGenerator::Synchronize() {
// to a better "new objective" if the base solution wasn't the best one.
//
// This might not be a final solution, but it does work ok for now.
const IntegerValue best_objective_improvement =
IsRelaxationGenerator()
? IntegerValue(CapSub(data.new_objective_bound.value(),
data.initial_best_objective_bound.value()))
: IntegerValue(CapSub(data.initial_best_objective.value(),
data.new_objective.value()));
const IntegerValue best_objective_improvement = IntegerValue(CapSub(
data.initial_best_objective.value(), data.new_objective.value()));
if (best_objective_improvement > 0) {
num_consecutive_non_improving_calls_ = 0;
} else {
@@ -1867,166 +1863,5 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate(
return neighborhood;
}
Neighborhood ConsecutiveConstraintsRelaxationNeighborhoodGenerator::Generate(
const CpSolverResponse& /*initial_solution*/, double difficulty,
absl::BitGenRef random) {
std::vector<int> removable_constraints;
const int num_constraints = helper_.ModelProto().constraints_size();
removable_constraints.reserve(num_constraints);
for (int c = 0; c < num_constraints; ++c) {
// Removing intervals is not easy because other constraint might require
// them, so for now, we don't remove them.
if (helper_.ModelProto().constraints(c).constraint_case() ==
ConstraintProto::kInterval) {
continue;
}
removable_constraints.push_back(c);
}
const int target_size =
std::round((1.0 - difficulty) * removable_constraints.size());
const int random_start_index =
absl::Uniform<int>(random, 0, removable_constraints.size());
std::vector<int> removed_constraints;
removed_constraints.reserve(target_size);
int c = random_start_index;
while (removed_constraints.size() < target_size) {
removed_constraints.push_back(removable_constraints[c]);
++c;
if (c == removable_constraints.size()) {
c = 0;
}
}
return helper_.RemoveMarkedConstraints(removed_constraints);
}
WeightedRandomRelaxationNeighborhoodGenerator::
WeightedRandomRelaxationNeighborhoodGenerator(
NeighborhoodGeneratorHelper const* helper, const std::string& name)
: NeighborhoodGenerator(name, helper) {
std::vector<int> removable_constraints;
const int num_constraints = helper_.ModelProto().constraints_size();
constraint_weights_.reserve(num_constraints);
// TODO(user): Experiment with different starting weights.
for (int c = 0; c < num_constraints; ++c) {
switch (helper_.ModelProto().constraints(c).constraint_case()) {
case ConstraintProto::kCumulative:
case ConstraintProto::kAllDiff:
case ConstraintProto::kElement:
case ConstraintProto::kRoutes:
case ConstraintProto::kCircuit:
constraint_weights_.push_back(3.0);
num_removable_constraints_++;
break;
case ConstraintProto::kBoolOr:
case ConstraintProto::kBoolAnd:
case ConstraintProto::kBoolXor:
case ConstraintProto::kIntProd:
case ConstraintProto::kIntDiv:
case ConstraintProto::kIntMod:
case ConstraintProto::kLinMax:
case ConstraintProto::kNoOverlap:
case ConstraintProto::kNoOverlap2D:
constraint_weights_.push_back(2.0);
num_removable_constraints_++;
break;
case ConstraintProto::kLinear:
case ConstraintProto::kTable:
case ConstraintProto::kAutomaton:
case ConstraintProto::kInverse:
case ConstraintProto::kReservoir:
case ConstraintProto::kAtMostOne:
case ConstraintProto::kExactlyOne:
constraint_weights_.push_back(1.0);
num_removable_constraints_++;
break;
case ConstraintProto::CONSTRAINT_NOT_SET:
case ConstraintProto::kDummyConstraint:
case ConstraintProto::kInterval:
// Removing intervals is not easy because other constraint might require
// them, so for now, we don't remove them.
constraint_weights_.push_back(0.0);
break;
}
}
}
void WeightedRandomRelaxationNeighborhoodGenerator::
AdditionalProcessingOnSynchronize(const SolveData& solve_data) {
const IntegerValue best_objective_improvement =
solve_data.new_objective_bound - solve_data.initial_best_objective_bound;
const std::vector<int>& removed_constraints =
removed_constraints_[solve_data.neighborhood_id];
// Heuristic: We change the weights of the removed constraints if the
// neighborhood is solved (status is OPTIMAL or INFEASIBLE) or we observe an
// improvement in objective bounds. Otherwise we assume that the
// difficulty/time wasn't right for us to record feedbacks.
//
// If the objective bounds are improved, we bump up the weights. If the
// objective bounds are worse and the problem status is OPTIMAL, we bump down
// the weights. Otherwise if the new objective bounds are same as current
// bounds (which happens a lot on some instances), we do not update the
// weights as we do not have a clear signal whether the constraints removed
// were good choices or not.
// TODO(user): We can improve this heuristic with more experiments.
if (best_objective_improvement > 0) {
// Bump up the weights of all removed constraints.
for (int c : removed_constraints) {
if (constraint_weights_[c] <= 90.0) {
constraint_weights_[c] += 10.0;
} else {
constraint_weights_[c] = 100.0;
}
}
} else if (solve_data.status == CpSolverStatus::OPTIMAL &&
best_objective_improvement < 0) {
// Bump down the weights of all removed constraints.
for (int c : removed_constraints) {
if (constraint_weights_[c] > 0.5) {
constraint_weights_[c] -= 0.5;
}
}
}
removed_constraints_.erase(solve_data.neighborhood_id);
}
Neighborhood WeightedRandomRelaxationNeighborhoodGenerator::Generate(
const CpSolverResponse& /*initial_solution*/, double difficulty,
absl::BitGenRef random) {
const int target_size =
std::round((1.0 - difficulty) * num_removable_constraints_);
std::vector<int> removed_constraints;
// Generate a random number between (0,1) = u[i] and use score[i] =
// u[i]^(1/w[i]) and then select top k items with largest scores.
// Reference: https://utopia.duth.gr/~pefraimi/research/data/2007EncOfAlg.pdf
std::vector<std::pair<double, int>> constraint_removal_scores;
std::uniform_real_distribution<double> random_var(0.0, 1.0);
for (int c = 0; c < constraint_weights_.size(); ++c) {
if (constraint_weights_[c] <= 0) continue;
const double u = random_var(random);
const double score = std::pow(u, (1 / constraint_weights_[c]));
constraint_removal_scores.push_back({score, c});
}
std::sort(constraint_removal_scores.rbegin(),
constraint_removal_scores.rend());
for (int i = 0; i < target_size; ++i) {
removed_constraints.push_back(constraint_removal_scores[i].second);
}
Neighborhood result = helper_.RemoveMarkedConstraints(removed_constraints);
absl::MutexLock mutex_lock(&generator_mutex_);
result.id = next_available_id_;
next_available_id_++;
removed_constraints_.insert({result.id, removed_constraints});
return result;
}
} // namespace sat
} // namespace operations_research

74
ortools/sat/cp_model_lns.h
View File

@@ -334,10 +334,6 @@ class NeighborhoodGenerator {
// Returns true if the neighborhood generator can generate a neighborhood.
virtual bool ReadyToGenerate() const;
// Returns true if the neighborhood generator generates relaxation of the
// given problem.
virtual bool IsRelaxationGenerator() const { return false; }
// Uses UCB1 algorithm to compute the score (Multi armed bandit problem).
// Details are at
// https://lilianweng.github.io/lil-log/2018/01/23/the-multi-armed-bandit-problem-and-its-solutions.html.
@@ -350,10 +346,6 @@ class NeighborhoodGenerator {
// Adds solve data about one "solved" neighborhood.
struct SolveData {
// Neighborhood Id. Used to identify the neighborhood by a generator.
// Currently only used by WeightedRandomRelaxationNeighborhoodGenerator.
int64_t neighborhood_id = 0;
// The status of the sub-solve.
CpSolverStatus status = CpSolverStatus::UNKNOWN;
@@ -379,22 +371,14 @@ class NeighborhoodGenerator {
IntegerValue base_objective = IntegerValue(0);
IntegerValue new_objective = IntegerValue(0);
// Bounds data is only used by relaxation neighborhoods.
IntegerValue initial_best_objective_bound = IntegerValue(0);
IntegerValue new_objective_bound = IntegerValue(0);
// This is just used to construct a deterministic order for the updates.
bool operator<(const SolveData& o) const {
return std::tie(status, difficulty, deterministic_limit,
deterministic_time, initial_best_objective,
base_objective, new_objective,
initial_best_objective_bound, new_objective_bound,
neighborhood_id) <
base_objective, new_objective) <
std::tie(o.status, o.difficulty, o.deterministic_limit,
o.deterministic_time, o.initial_best_objective,
o.base_objective, o.new_objective,
o.initial_best_objective_bound, o.new_objective_bound,
o.neighborhood_id);
o.base_objective, o.new_objective);
}
};
void AddSolveData(SolveData data) {
@@ -690,60 +674,6 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator {
SharedIncompleteSolutionManager* incomplete_solutions_;
};
// Generates a relaxation of the original model by removing a consecutive span
// of constraints starting at a random index. The number of constraints removed
// is in sync with the difficulty passed to the generator.
class ConsecutiveConstraintsRelaxationNeighborhoodGenerator
: public NeighborhoodGenerator {
public:
explicit ConsecutiveConstraintsRelaxationNeighborhoodGenerator(
NeighborhoodGeneratorHelper const* helper, const std::string& name)
: NeighborhoodGenerator(name, helper) {}
Neighborhood Generate(const CpSolverResponse& initial_solution,
double difficulty, absl::BitGenRef random) final;
bool IsRelaxationGenerator() const override { return true; }
bool ReadyToGenerate() const override { return true; }
};
// Generates a relaxation of the original model by removing some constraints
// randomly with a given weight for each constraint that controls the
// probability of constraint getting removed. The number of constraints removed
// is in sync with the difficulty passed to the generator. Higher weighted
// constraints are more likely to get removed.
class WeightedRandomRelaxationNeighborhoodGenerator
: public NeighborhoodGenerator {
public:
WeightedRandomRelaxationNeighborhoodGenerator(
NeighborhoodGeneratorHelper const* helper, const std::string& name);
// Generates the neighborhood as described above. Also stores the removed
// constraints indices for adjusting the weights.
Neighborhood Generate(const CpSolverResponse& initial_solution,
double difficulty, absl::BitGenRef random) final;
bool IsRelaxationGenerator() const override { return true; }
bool ReadyToGenerate() const override { return true; }
private:
// Adjusts the weights of the constraints removed to get the neighborhood
// based on the solve_data.
void AdditionalProcessingOnSynchronize(const SolveData& solve_data) override
ABSL_EXCLUSIVE_LOCKS_REQUIRED(generator_mutex_);
// Higher weighted constraints are more likely to get removed.
std::vector<double> constraint_weights_;
int num_removable_constraints_ = 0;
// Indices of the removed constraints per generated neighborhood.
absl::flat_hash_map<int64_t, std::vector<int>> removed_constraints_
ABSL_GUARDED_BY(generator_mutex_);
// TODO(user): Move this to parent class if other generators start using
// feedbacks.
int64_t next_available_id_ ABSL_GUARDED_BY(generator_mutex_) = 0;
};
} // namespace sat
} // namespace operations_research

43
ortools/sat/cp_model_loader.cc
View File

@@ -383,6 +383,7 @@ void ExtractEncoding(const CpModelProto& model_proto, Model* m) {
// if some of the literal view used in the LP are created later, but that
// should be fixable via calls to implied_bounds->NotifyNewIntegerView().
auto* implied_bounds = m->GetOrCreate<ImpliedBounds>();
auto* detector = m->GetOrCreate<ProductDetector>();
// Detection of literal equivalent to (i_var >= bound). We also collect
// all the half-refied part and we will sort the vector for detection of the
@@ -468,6 +469,10 @@ void ExtractEncoding(const CpModelProto& model_proto, Model* m) {
{
const Domain inter = domain.IntersectionWith(domain_if_enforced);
if (!inter.IsEmpty() && inter.Min() == inter.Max()) {
if (inter.Min() == 0) {
detector->ProcessConditionalZero(enforcement_literal,
mapping->Integer(var));
}
var_to_equalities[var].push_back(
{&ct, enforcement_literal, inter.Min(), true});
if (domain.Contains(inter.Min())) {
@@ -869,6 +874,9 @@ void LoadBoolOrConstraint(const ConstraintProto& ct, Model* m) {
literals.push_back(mapping->Literal(ref).Negated());
}
m->Add(ClauseConstraint(literals));
if (literals.size() == 3) {
m->GetOrCreate<ProductDetector>()->ProcessTernaryClause(literals);
}
}
void LoadBoolAndConstraint(const ConstraintProto& ct, Model* m) {
@@ -894,7 +902,11 @@ void LoadAtMostOneConstraint(const ConstraintProto& ct, Model* m) {
void LoadExactlyOneConstraint(const ConstraintProto& ct, Model* m) {
auto* mapping = m->GetOrCreate<CpModelMapping>();
CHECK(!HasEnforcementLiteral(ct)) << "Not supported.";
m->Add(ExactlyOneConstraint(mapping->Literals(ct.exactly_one().literals())));
const auto& literals = mapping->Literals(ct.exactly_one().literals());
m->Add(ExactlyOneConstraint(literals));
if (literals.size() == 3) {
m->GetOrCreate<ProductDetector>()->ProcessTernaryExactlyOne(literals);
}
}
void LoadBoolXorConstraint(const ConstraintProto& ct, Model* m) {
@@ -975,6 +987,18 @@ void LoadEquivalenceNeqAC(const std::vector<Literal> enforcement_literal,
}
}
bool IsPartOfProductEncoding(const ConstraintProto& ct) {
if (ct.enforcement_literal().size() != 1) return false;
if (ct.linear().vars().size() > 2) return false;
if (ct.linear().domain().size() != 2) return false;
if (ct.linear().domain(0) != 0) return false;
if (ct.linear().domain(1) != 0) return false;
for (const int64_t coeff : ct.linear().coeffs()) {
if (std::abs(coeff) != 1) return false;
}
return true;
}
} // namespace
void LoadLinearConstraint(const ConstraintProto& ct, Model* m) {
@@ -995,6 +1019,23 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) {
return;
}
if (IsPartOfProductEncoding(ct)) {
const Literal l = mapping->Literal(ct.enforcement_literal(0));
auto* detector = m->GetOrCreate<ProductDetector>();
if (ct.linear().vars().size() == 1) {
// TODO(user): Actually this should never be called since we process
// linear1 in ExtractEncoding().
detector->ProcessConditionalZero(l,
mapping->Integer(ct.linear().vars(0)));
} else if (ct.linear().vars().size() == 2) {
const IntegerVariable x = mapping->Integer(ct.linear().vars(0));
const IntegerVariable y = mapping->Integer(ct.linear().vars(1));
detector->ProcessConditionalEquality(
l, x,
ct.linear().coeffs(0) == ct.linear().coeffs(1) ? NegationOf(y) : y);
}
}
auto* integer_trail = m->GetOrCreate<IntegerTrail>();
const std::vector<IntegerVariable> vars =
mapping->Integers(ct.linear().vars());

32
ortools/sat/cp_model_presolve.cc
View File

@@ -2483,6 +2483,15 @@ bool CpModelPresolver::PresolveDiophantine(ConstraintProto* ct) {
// This tries to decompose the constraint into coeff * part1 + part2 and show
// that the value that part2 take is not important, thus the constraint can
// only be transformed on a constraint on the first part.
//
// TODO(user): Improve !! we miss simple case like x + 47 y + 50 z >= 50
// for positive variables. We should remove x, and ideally we should rewrite
// this as y + 2z >= 2 if we can show that its relaxation is just better?
// We should at least see that it is the same as 47y + 50 z >= 48.
//
// TODO(user): One easy algo is to first remove all enforcement term (even
// non-Boolean one) before applying the algo here and then re-linearize the
// non-Boolean terms.
void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint(
int c, ConstraintProto* ct) {
if (ct->constraint_case() != ConstraintProto::kLinear) return;
@@ -3282,7 +3291,7 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint(
// Note that for now the default is false, and also there are problem calling
// GetOrCreateVarValueEncoding() after expansion because we might have removed
// the variable used in the encoding.
const bool only_booleans =
const bool only_extract_booleans =
!context_->params().presolve_extract_integer_enforcement() ||
context_->ModelIsExpanded();
@@ -3302,10 +3311,25 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint(
const bool is_boolean = context_->CanBeUsedAsLiteral(ref);
if (context_->IsFixed(ref) || coeff < threshold ||
(only_booleans && !is_boolean)) {
// We keep this term.
(only_extract_booleans && !is_boolean)) {
mutable_arg->set_vars(new_size, mutable_arg->vars(i));
mutable_arg->set_coeffs(new_size, mutable_arg->coeffs(i));
// We keep this term but reduces its coeff.
// This is only for the case where only_extract_booleans == true.
if (coeff > threshold) {
context_->UpdateRuleStats("linear: coefficient strenghtening.");
if (lower_bounded) {
// coeff * (X - LB + LB) -> threshold * (X - LB) + coeff * LB
rhs_offset -= (coeff - threshold) * context_->MinOf(ref);
} else {
// coeff * (X - UB + UB) -> threshold * (X - UB) + coeff * UB
rhs_offset -= (coeff - threshold) * context_->MaxOf(ref);
}
mutable_arg->set_coeffs(new_size,
arg.coeffs(i) > 0 ? threshold : -threshold);
} else {
mutable_arg->set_coeffs(new_size, arg.coeffs(i));
}
++new_size;
continue;
}

3
ortools/sat/cp_model_search.cc
View File

@@ -693,8 +693,7 @@ std::vector<SatParameters> GetDiverseSetOfParameters(
// worker for them.
int target = base_params.num_workers();
if (!base_params.interleave_search() &&
(base_params.use_rins_lns() || base_params.use_relaxation_lns() ||
base_params.use_feasibility_pump())) {
(base_params.use_rins_lns() || base_params.use_feasibility_pump())) {
target = std::max(1, base_params.num_workers() - 1);
}
if (!base_params.interleave_search() && result.size() > target) {

229
ortools/sat/cp_model_solver.cc
View File

@@ -67,6 +67,7 @@
#include "ortools/sat/drat_checker.h"
#include "ortools/sat/drat_proof_handler.h"
#include "ortools/sat/feasibility_pump.h"
#include "ortools/sat/implied_bounds.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/integer_expr.h"
#include "ortools/sat/integer_search.h"
@@ -598,16 +599,14 @@ void LoadDebugSolution(const CpModelProto& model_proto, Model* model) {
const IntegerVariable objective_var = objective_def->objective_var;
const int64_t objective_value =
ComputeInnerObjective(model_proto.objective(), response);
ComputeInnerObjective(model_proto.objective(), response.solution());
debug_solution[objective_var] = objective_value;
debug_solution[NegationOf(objective_var)] = -objective_value;
#endif // __PORTABLE_PLATFORM__
}
void FillSolutionInResponse(const CpModelProto& model_proto, const Model& model,
CpSolverResponse* response) {
response->clear_solution();
std::vector<int64_t> GetSolutionValues(const CpModelProto& model_proto,
const Model& model) {
auto* mapping = model.Get<CpModelMapping>();
auto* trail = model.Get<Trail>();
@@ -636,9 +635,7 @@ void FillSolutionInResponse(const CpModelProto& model_proto, const Model& model,
// TODO(user): Checks against initial model.
CHECK(SolutionIsFeasible(model_proto, solution));
}
for (const int64_t value : solution) {
response->add_solution(value);
}
return solution;
}
namespace {
@@ -1269,6 +1266,9 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) {
model->GetOrCreate<IntegerEncoder>()
->AddAllImplicationsBetweenAssociatedLiterals();
if (!sat_solver->FinishPropagation()) return unsat();
model->GetOrCreate<ProductDetector>()->ProcessImplicationGraph(
model->GetOrCreate<BinaryImplicationGraph>());
}
void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) {
@@ -1522,10 +1522,9 @@ void LoadCpModel(const CpModelProto& model_proto, Model* model) {
const std::string solution_info = model->Name();
const auto solution_observer = [&model_proto, model, solution_info,
shared_response_manager]() {
CpSolverResponse response;
FillSolutionInResponse(model_proto, *model, &response);
response.set_solution_info(solution_info);
shared_response_manager->NewSolution(response, model);
const std::vector<int64_t> solution =
GetSolutionValues(model_proto, *model);
shared_response_manager->NewSolution(solution, solution_info, model);
};
const auto& objective = *model->GetOrCreate<ObjectiveDefinition>();
@@ -1557,10 +1556,9 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, Model* model) {
const std::string& solution_info = model->Name();
const auto solution_observer = [&model_proto, &model, &solution_info,
&shared_response_manager]() {
CpSolverResponse response;
FillSolutionInResponse(model_proto, *model, &response);
response.set_solution_info(solution_info);
shared_response_manager->NewSolution(response, model);
const std::vector<int64_t> solution =
GetSolutionValues(model_proto, *model);
shared_response_manager->NewSolution(solution, solution_info, model);
};
// Reconfigure search heuristic if it was changed.
@@ -1694,10 +1692,10 @@ void QuickSolveWithHint(const CpModelProto& model_proto, Model* model) {
const std::string& solution_info = model->Name();
if (status == SatSolver::Status::FEASIBLE) {
CpSolverResponse response;
FillSolutionInResponse(model_proto, *model, &response);
response.set_solution_info(absl::StrCat(solution_info, " [hint]"));
shared_response_manager->NewSolution(response, model);
const std::vector<int64_t> solution =
GetSolutionValues(model_proto, *model);
shared_response_manager->NewSolution(
solution, absl::StrCat(solution_info, " [hint]"), model);
if (!model_proto.has_objective()) {
if (parameters->enumerate_all_solutions()) {
@@ -1828,18 +1826,17 @@ void MinimizeL1DistanceWithHint(const CpModelProto& model_proto, Model* model) {
const std::string& solution_info = model->Name();
if (status == SatSolver::Status::FEASIBLE) {
CpSolverResponse response;
FillSolutionInResponse(model_proto, local_model, &response);
const std::vector<int64_t> solution =
GetSolutionValues(model_proto, local_model);
if (DEBUG_MODE) {
CpSolverResponse updated_response;
FillSolutionInResponse(updated_model_proto, local_model,
&updated_response);
const std::vector<int64_t> updated_solution =
GetSolutionValues(updated_model_proto, local_model);
LOG(INFO) << "Found solution with repaired hint penalty = "
<< ComputeInnerObjective(updated_model_proto.objective(),
updated_response);
updated_solution);
}
response.set_solution_info(absl::StrCat(solution_info, " [repaired]"));
shared_response_manager->NewSolution(response, &local_model);
shared_response_manager->NewSolution(
solution, absl::StrCat(solution_info, " [repaired]"), &local_model);
}
}
@@ -2512,8 +2509,6 @@ class LnsSolver : public SubSolver {
if (!neighborhood.is_generated) return;
data.neighborhood_id = neighborhood.id;
const int64_t num_calls = std::max(int64_t{1}, generator_->num_calls());
const double fully_solved_proportion =
static_cast<double>(generator_->num_fully_solved_calls()) /
@@ -2522,7 +2517,7 @@ class LnsSolver : public SubSolver {
if (!neighborhood.source_info.empty()) {
absl::StrAppend(&source_info, "_", neighborhood.source_info);
}
const std::string solution_info = absl::StrFormat(
const std::string lns_info = absl::StrFormat(
"%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", source_info, data.difficulty,
task_id, data.deterministic_limit, fully_solved_proportion);
@@ -2534,7 +2529,7 @@ class LnsSolver : public SubSolver {
local_params.set_symmetry_level(0);
local_params.set_solution_pool_size(0);
Model local_model(solution_info);
Model local_model(lns_info);
*(local_model.GetOrCreate<SatParameters>()) = local_params;
TimeLimit* local_time_limit = local_model.GetOrCreate<TimeLimit>();
local_time_limit->ResetLimitFromParameters(local_params);
@@ -2616,106 +2611,49 @@ class LnsSolver : public SubSolver {
local_response_manager->SetResponseStatus(presolve_status);
}
CpSolverResponse local_response = local_response_manager->GetResponse();
const std::string solution_info = local_response.solution_info();
std::vector<int64_t> solution_values(local_response.solution().begin(),
local_response.solution().end());
data.status = local_response.status();
// TODO(user): we actually do not need to postsolve if the solution is
// not going to be used...
if (local_params.cp_model_presolve() &&
(local_response.status() == CpSolverStatus::OPTIMAL ||
local_response.status() == CpSolverStatus::FEASIBLE)) {
std::vector<int64_t> solution(local_response.solution().begin(),
local_response.solution().end());
PostsolveResponseWrapper(local_params,
helper_->ModelProto().variables_size(),
mapping_proto, postsolve_mapping, &solution);
local_response.mutable_solution()->Assign(solution.begin(),
solution.end());
(data.status == CpSolverStatus::OPTIMAL ||
data.status == CpSolverStatus::FEASIBLE)) {
PostsolveResponseWrapper(
local_params, helper_->ModelProto().variables_size(), mapping_proto,
postsolve_mapping, &solution_values);
local_response.mutable_solution()->Assign(solution_values.begin(),
solution_values.end());
}
data.status = local_response.status();
data.deterministic_time = local_time_limit->GetElapsedDeterministicTime();
bool new_solution = false;
bool display_lns_info = VLOG_IS_ON(2);
if (generator_->IsRelaxationGenerator()) {
bool has_feasible_solution = false;
if (local_response.status() == CpSolverStatus::OPTIMAL ||
local_response.status() == CpSolverStatus::FEASIBLE) {
has_feasible_solution = true;
}
const std::vector<int64_t> solution(local_response.solution().begin(),
local_response.solution().end());
if (local_response.status() == CpSolverStatus::INFEASIBLE) {
shared_->response->NotifyThatImprovingProblemIsInfeasible(
local_response.solution_info());
}
if (shared_->model_proto->has_objective()) {
// TODO(user): This is not deterministic since it is updated without
// synchronization! So we shouldn't base the LNS score out of that.
const IntegerValue current_obj_lb =
shared_->response->GetInnerObjectiveLowerBound();
const IntegerValue local_obj_lb =
local_response_manager->GetInnerObjectiveLowerBound();
const double scaled_local_obj_bound = ScaleObjectiveValue(
lns_fragment.objective(), local_obj_lb.value());
// Update the bound.
const IntegerValue new_inner_obj_lb = IntegerValue(
std::ceil(UnscaleObjectiveValue(shared_->model_proto->objective(),
scaled_local_obj_bound) -
1e-6));
data.new_objective_bound = new_inner_obj_lb;
data.initial_best_objective_bound = current_obj_lb;
if (new_inner_obj_lb > current_obj_lb) {
shared_->response->UpdateInnerObjectiveBounds(
solution_info, new_inner_obj_lb, kMaxIntegerValue);
}
}
// If we have a solution of the relaxed problem, we check if it is also
// a valid solution of the non-relaxed one.
if (has_feasible_solution) {
if (SolutionIsFeasible(
*shared_->model_proto,
std::vector<int64_t>(local_response.solution().begin(),
local_response.solution().end()))) {
shared_->response->NewSolution(local_response,
/*model=*/nullptr);
// Mark the solution optimal if the relaxation status is optimal.
if (local_response.status() == CpSolverStatus::OPTIMAL) {
shared_->response->NotifyThatImprovingProblemIsInfeasible(
local_response.solution_info());
shared_->time_limit->Stop();
}
}
shared_->relaxation_solutions->NewRelaxationSolution(local_response);
}
} else {
const std::vector<int64_t> solution(local_response.solution().begin(),
local_response.solution().end());
if (!solution.empty()) {
// A solution that does not pass our validator indicates a bug. We
// abort and dump the problematic model to facilitate debugging.
//
// TODO(user): In a production environment, we should probably just
// ignore this fragment and continue.
const bool feasible =
SolutionIsFeasible(*shared_->model_proto, solution);
if (!feasible) {
if (absl::GetFlag(FLAGS_cp_model_dump_problematic_lns)) {
const std::string name =
absl::StrCat(absl::GetFlag(FLAGS_cp_model_dump_prefix),
debug_copy.name(), ".pb.txt");
LOG(INFO) << "Dumping problematic LNS model to '" << name << "'.";
CHECK_OK(file::SetTextProto(name, debug_copy, file::Defaults()));
}
LOG(FATAL) << "Infeasible LNS solution! " << solution_info
<< " solved with params "
<< local_params.ShortDebugString();
if (!solution.empty()) {
// A solution that does not pass our validator indicates a bug. We
// abort and dump the problematic model to facilitate debugging.
//
// TODO(user): In a production environment, we should probably just
// ignore this fragment and continue.
const bool feasible =
SolutionIsFeasible(*shared_->model_proto, solution_values);
if (!feasible) {
if (absl::GetFlag(FLAGS_cp_model_dump_problematic_lns)) {
const std::string name =
absl::StrCat(absl::GetFlag(FLAGS_cp_model_dump_prefix),
debug_copy.name(), ".pb.txt");
LOG(INFO) << "Dumping problematic LNS model to '" << name << "'.";
CHECK_OK(file::SetTextProto(name, debug_copy, file::Defaults()));
}
LOG(FATAL) << "Infeasible LNS solution! " << solution_info
<< " solved with params "
<< local_params.ShortDebugString();
}
// Special case if we solved a part of the full problem!
@@ -2729,41 +2667,42 @@ class LnsSolver : public SubSolver {
// - It should be fine if all our generator are fully or not at
// all included in the variable we are fixing. So we can relax the
// test here. Try on z26.mps or spot5_1401.fzn.
if (local_response.status() == CpSolverStatus::OPTIMAL &&
!shared_->model_proto->has_symmetry() && !solution.empty() &&
if (data.status == CpSolverStatus::OPTIMAL &&
!shared_->model_proto->has_symmetry() && !solution_values.empty() &&
neighborhood.is_simple &&
!neighborhood.variables_that_can_be_fixed_to_local_optimum
.empty()) {
display_lns_info = true;
shared_->bounds->FixVariablesFromPartialSolution(
solution,
solution_values,
neighborhood.variables_that_can_be_fixed_to_local_optimum);
}
// Finish to fill the SolveData now that the local solve is done.
data.new_objective = data.base_objective;
if (local_response.status() == CpSolverStatus::OPTIMAL ||
local_response.status() == CpSolverStatus::FEASIBLE) {
if (data.status == CpSolverStatus::OPTIMAL ||
data.status == CpSolverStatus::FEASIBLE) {
data.new_objective = IntegerValue(ComputeInnerObjective(
shared_->model_proto->objective(), local_response));
shared_->model_proto->objective(), solution_values));
}
// Report any feasible solution we have. Optimization: We don't do that
// if we just recovered the base solution.
if (local_response.status() == CpSolverStatus::OPTIMAL ||
local_response.status() == CpSolverStatus::FEASIBLE) {
if (data.status == CpSolverStatus::OPTIMAL ||
data.status == CpSolverStatus::FEASIBLE) {
const std::vector<int64_t> base_solution(
base_response.solution().begin(), base_response.solution().end());
if (solution != base_solution) {
if (solution_values != base_solution) {
new_solution = true;
shared_->response->NewSolution(local_response, /*model=*/nullptr);
shared_->response->NewSolution(solution_values, solution_info,
/*model=*/nullptr);
}
}
if (!neighborhood.is_reduced &&
(local_response.status() == CpSolverStatus::OPTIMAL ||
local_response.status() == CpSolverStatus::INFEASIBLE)) {
(data.status == CpSolverStatus::OPTIMAL ||
data.status == CpSolverStatus::INFEASIBLE)) {
shared_->response->NotifyThatImprovingProblemIsInfeasible(
local_response.solution_info());
solution_info);
shared_->time_limit->Stop();
}
}
@@ -2777,11 +2716,11 @@ class LnsSolver : public SubSolver {
const double base_obj = ScaleObjectiveValue(
shared_->model_proto->objective(),
ComputeInnerObjective(shared_->model_proto->objective(),
base_response));
base_response.solution()));
const double new_obj = ScaleObjectiveValue(
shared_->model_proto->objective(),
ComputeInnerObjective(shared_->model_proto->objective(),
local_response));
solution_values));
absl::StrAppend(&s, " [new_sol:", base_obj, " -> ", new_obj, "]");
}
if (neighborhood.is_simple) {
@@ -2831,13 +2770,6 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
std::unique_ptr<SharedRelaxationSolutionRepository>
shared_relaxation_solutions;
if (params.use_relaxation_lns()) {
shared_relaxation_solutions =
std::make_unique<SharedRelaxationSolutionRepository>(
/*num_solutions_to_keep=*/10);
global_model->Register<SharedRelaxationSolutionRepository>(
shared_relaxation_solutions.get());
}
auto shared_lp_solutions = std::make_unique<SharedLPSolutionRepository>(
/*num_solutions_to_keep=*/10);
@@ -3041,19 +2973,6 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
absl::StrCat("rens_lns_", local_params.name())),
local_params, helper, &shared));
}
if (params.use_relaxation_lns()) {
subsolvers.push_back(std::make_unique<LnsSolver>(
std::make_unique<
ConsecutiveConstraintsRelaxationNeighborhoodGenerator>(
helper, absl::StrCat("rnd_rel_lns_", local_params.name())),
local_params, helper, &shared));
subsolvers.push_back(std::make_unique<LnsSolver>(
std::make_unique<WeightedRandomRelaxationNeighborhoodGenerator>(
helper, absl::StrCat("wgt_rel_lns_", local_params.name())),
local_params, helper, &shared));
}
}
// Add a synchronization point for the gap integral that is executed last.

4
ortools/sat/cp_model_utils.cc
View File

@@ -517,14 +517,14 @@ std::vector<int> UsedIntervals(const ConstraintProto& ct) {
}
int64_t ComputeInnerObjective(const CpObjectiveProto& objective,
const CpSolverResponse& response) {
absl::Span<const int64_t> solution) {
int64_t objective_value = 0;
for (int i = 0; i < objective.vars_size(); ++i) {
int64_t coeff = objective.coeffs(i);
const int ref = objective.vars(i);
const int var = PositiveRef(ref);
if (!RefIsPositive(ref)) coeff = -coeff;
objective_value += coeff * response.solution()[var];
objective_value += coeff * solution[var];
}
return objective_value;
}

2
ortools/sat/cp_model_utils.h
View File

@@ -164,7 +164,7 @@ inline double UnscaleObjectiveValue(const CpObjectiveProto& proto,
// This is the objective without offset and scaling. Call ScaleObjectiveValue()
// to get the user facing objective.
int64_t ComputeInnerObjective(const CpObjectiveProto& objective,
const CpSolverResponse& response);
absl::Span<const int64_t> solution);
// Returns true if a linear expression can be reduced to a single ref.
bool ExpressionContainsSingleRef(const LinearExpressionProto& expr);

263
ortools/sat/implied_bounds.cc
View File

@@ -16,6 +16,8 @@
#include <stdint.h>
#include <algorithm>
#include <array>
#include <bitset>
#include <limits>
#include <optional>
#include <string>
@@ -41,10 +43,14 @@ namespace sat {
// Just display some global statistics on destruction.
ImpliedBounds::~ImpliedBounds() {
VLOG(1) << num_deductions_ << " enqueued deductions.";
VLOG(1) << bounds_.size() << " implied bounds stored.";
VLOG(1) << num_enqueued_in_var_to_bounds_
<< " implied bounds with view stored.";
if (!VLOG_IS_ON(1)) return;
if (shared_stats_ == nullptr) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back({"implied_bound/num_deductions", num_deductions_});
stats.push_back({"implied_bound/num_stored", bounds_.size()});
stats.push_back(
{"implied_bound/num_stored_with_view", num_enqueued_in_var_to_bounds_});
shared_stats_->AddStats(stats);
}
void ImpliedBounds::Add(Literal literal, IntegerLiteral integer_literal) {
@@ -124,7 +130,7 @@ void ImpliedBounds::Add(Literal literal, IntegerLiteral integer_literal) {
level_zero_lower_bounds_[var] = deduction;
new_level_zero_bounds_.Set(var);
VLOG(1) << "Deduction old: "
VLOG(2) << "Deduction old: "
<< IntegerLiteral::GreaterOrEqual(
var, integer_trail_->LevelZeroLowerBound(var))
<< " new: " << IntegerLiteral::GreaterOrEqual(var, deduction);
@@ -339,7 +345,7 @@ std::vector<LiteralValueValue> TryToReconcileSize2Encodings(
const std::vector<ValueLiteralPair>& right_enc =
integer_encoder->FullDomainEncoding(right.var);
if (left_enc.size() != 2 || right_enc.size() != 2) {
VLOG(3) << "encodings are not fully propagated";
VLOG(2) << "encodings are not fully propagated";
return terms;
}
@@ -426,7 +432,7 @@ std::vector<LiteralValueValue> TryToDecomposeProduct(
}
if (compatible_keys.size() > 1) {
VLOG(1) << "More than one exactly_one involved in the encoding of the two "
VLOG(3) << "More than one exactly_one involved in the encoding of the two "
"variables";
}
@@ -515,5 +521,248 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left,
return true;
}
ProductDetector::ProductDetector(Model* model)
: enabled_(model->GetOrCreate<SatParameters>()->linearization_level() > 1),
sat_solver_(model->GetOrCreate<SatSolver>()),
trail_(model->GetOrCreate<Trail>()),
integer_trail_(model->GetOrCreate<IntegerTrail>()),
integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
shared_stats_(model->GetOrCreate<SharedStatistics>()) {}
ProductDetector::~ProductDetector() {
if (!VLOG_IS_ON(1)) return;
if (shared_stats_ == nullptr) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back(
{"product_detector/num_processed_binary", num_processed_binary_});
stats.push_back(
{"product_detector/num_processed_exactly_one", num_processed_exo_});
stats.push_back(
{"product_detector/num_processed_ternary", num_processed_ternary_});
stats.push_back({"product_detector/num_trail_updates", num_trail_updates_});
stats.push_back({"product_detector/num_products", num_products_});
stats.push_back({"product_detector/num_conditional_equalities",
num_conditional_equalities_});
stats.push_back(
{"product_detector/num_conditional_zeros", num_conditional_zeros_});
stats.push_back({"product_detector/num_int_products", num_int_products_});
shared_stats_->AddStats(stats);
}
void ProductDetector::ProcessTernaryClause(
absl::Span<const Literal> ternary_clause) {
if (!enabled_) return;
if (ternary_clause.size() != 3) return;
++num_processed_ternary_;
candidates_[GetKey(ternary_clause[0].Index(), ternary_clause[1].Index())]
.push_back(ternary_clause[2].Index());
candidates_[GetKey(ternary_clause[0].Index(), ternary_clause[2].Index())]
.push_back(ternary_clause[1].Index());
candidates_[GetKey(ternary_clause[1].Index(), ternary_clause[2].Index())]
.push_back(ternary_clause[0].Index());
// We mark the literal of the ternary clause as seen.
// Only a => b with a seen need to be looked at.
for (const Literal l : ternary_clause) {
if (l.Index() >= seen_.size()) seen_.resize(l.Index() + 1);
seen_[l.Index()] = true;
}
}
void ProductDetector::ProcessTernaryExactlyOne(
absl::Span<const Literal> ternary_exo) {
if (!enabled_) return;
if (ternary_exo.size() != 3) return;
++num_processed_exo_;
ProcessNewProduct(ternary_exo[0].Index(), ternary_exo[1].NegatedIndex(),
ternary_exo[2].NegatedIndex());
ProcessNewProduct(ternary_exo[1].Index(), ternary_exo[0].NegatedIndex(),
ternary_exo[2].NegatedIndex());
ProcessNewProduct(ternary_exo[2].Index(), ternary_exo[0].NegatedIndex(),
ternary_exo[1].NegatedIndex());
}
// TODO(user): As product are discovered, we could remove entries from our
// hash maps!
void ProductDetector::ProcessBinaryClause(
absl::Span<const Literal> binary_clause) {
if (!enabled_) return;
if (binary_clause.size() != 2) return;
++num_processed_binary_;
const std::array<LiteralIndex, 2> key =
GetKey(binary_clause[0].NegatedIndex(), binary_clause[1].NegatedIndex());
std::array<LiteralIndex, 3> ternary;
for (const LiteralIndex l : candidates_[key]) {
ternary[0] = key[0];
ternary[1] = key[1];
ternary[2] = l;
std::sort(ternary.begin(), ternary.end());
const int l_index = ternary[0] == l ? 0 : ternary[1] == l ? 1 : 2;
std::bitset<3>& bs = detector_[ternary];
if (bs[l_index]) continue;
bs[l_index] = true;
if (bs[0] && bs[1] && l_index != 2) {
ProcessNewProduct(ternary[2], Literal(ternary[0]).NegatedIndex(),
Literal(ternary[1]).NegatedIndex());
}
if (bs[0] && bs[2] && l_index != 1) {
ProcessNewProduct(ternary[1], Literal(ternary[0]).NegatedIndex(),
Literal(ternary[2]).NegatedIndex());
}
if (bs[1] && bs[2] && l_index != 0) {
ProcessNewProduct(ternary[0], Literal(ternary[1]).NegatedIndex(),
Literal(ternary[2]).NegatedIndex());
}
}
}
void ProductDetector::ProcessImplicationGraph(BinaryImplicationGraph* graph) {
if (!enabled_) return;
for (LiteralIndex a(0); a < seen_.size(); ++a) {
if (!seen_[a]) continue;
if (trail_->Assignment().LiteralIsAssigned(Literal(a))) continue;
const Literal not_a = Literal(a).Negated();
for (const Literal b : graph->DirectImplications(Literal(a))) {
ProcessBinaryClause({not_a, b}); // a => b;
}
}
}
void ProductDetector::ProcessTrailAtLevelOne() {
if (!enabled_) return;
if (trail_->CurrentDecisionLevel() != 1) return;
++num_trail_updates_;
const SatSolver::Decision decision = sat_solver_->Decisions()[0];
if (decision.literal.Index() >= seen_.size() ||
!seen_[decision.literal.Index()]) {
return;
}
const Literal not_a = decision.literal.Negated();
const int current_index = trail_->Index();
for (int i = decision.trail_index + 1; i < current_index; ++i) {
const Literal b = (*trail_)[i];
ProcessBinaryClause({not_a, b});
}
}
LiteralIndex ProductDetector::GetProduct(Literal a, Literal b) const {
const auto it = products_.find(GetKey(a.Index(), b.Index()));
if (it == products_.end()) return kNoLiteralIndex;
return it->second;
}
std::array<LiteralIndex, 2> ProductDetector::GetKey(LiteralIndex a,
LiteralIndex b) const {
std::array<LiteralIndex, 2> key{a, b};
if (key[0] > key[1]) std::swap(key[0], key[1]);
return key;
}
void ProductDetector::ProcessNewProduct(LiteralIndex p, LiteralIndex a,
LiteralIndex b) {
// If many literal correspond to the same product, we just keep one.
++num_products_;
products_[GetKey(a, b)] = p;
// This is used by ProductIsLinearizable().
has_product_.insert(
GetKey(Literal(a).IsPositive() ? a : Literal(a).NegatedIndex(),
Literal(b).IsPositive() ? b : Literal(b).NegatedIndex()));
}
bool ProductDetector::ProductIsLinearizable(IntegerVariable a,
IntegerVariable b) const {
if (a == b) return true;
if (a == NegationOf(b)) return true;
// Otherwise, we need both a and b to be expressible as linear expression
// involving Booleans whose product is also expressible.
if (integer_trail_->InitialVariableDomain(a).Size() != 2) return false;
if (integer_trail_->InitialVariableDomain(b).Size() != 2) return false;
const LiteralIndex la =
integer_encoder_->GetAssociatedLiteral(IntegerLiteral::GreaterOrEqual(
a, integer_trail_->LevelZeroUpperBound(a)));
if (la == kNoLiteralIndex) return false;
const LiteralIndex lb =
integer_encoder_->GetAssociatedLiteral(IntegerLiteral::GreaterOrEqual(
b, integer_trail_->LevelZeroUpperBound(b)));
if (lb == kNoLiteralIndex) return false;
// Any product involving la/not(la) * lb/not(lb) can be used.
return has_product_.contains(
GetKey(Literal(la).IsPositive() ? la : Literal(la).NegatedIndex(),
Literal(lb).IsPositive() ? lb : Literal(lb).NegatedIndex()));
}
IntegerVariable ProductDetector::GetProduct(Literal a,
IntegerVariable b) const {
const auto it = int_products_.find({a.Index(), PositiveVariable(b)});
if (it == int_products_.end()) return kNoIntegerVariable;
return VariableIsPositive(b) ? it->second : NegationOf(it->second);
}
void ProductDetector::ProcessNewProduct(IntegerVariable p, Literal l,
IntegerVariable x) {
if (!VariableIsPositive(x)) {
x = NegationOf(x);
p = NegationOf(p);
}
// We only store one product if there are many.
++num_int_products_;
int_products_[{l.Index(), x}] = p;
}
void ProductDetector::ProcessConditionalEquality(Literal l, IntegerVariable x,
IntegerVariable y) {
++num_conditional_equalities_;
if (x == y) return;
// We process both possibilities (product = x or product = y).
for (int i = 0; i < 2; ++i) {
if (!VariableIsPositive(x)) {
x = NegationOf(x);
y = NegationOf(y);
}
bool seen = false;
// TODO(user): Linear scan can be bad if b => X = many other variables.
// Hopefully this will not be common.
std::vector<IntegerVariable>& others =
conditional_equalities_[{l.Index(), x}];
for (const IntegerVariable o : others) {
if (o == y) {
seen = true;
break;
}
}
if (!seen) {
others.push_back(y);
if (conditional_zeros_.contains({l.NegatedIndex(), x})) {
ProcessNewProduct(/*p=*/x, l, y);
}
}
std::swap(x, y);
}
}
void ProductDetector::ProcessConditionalZero(Literal l, IntegerVariable p) {
++num_conditional_zeros_;
p = PositiveVariable(p);
auto [_, inserted] = conditional_zeros_.insert({l.Index(), p});
if (inserted) {
const auto it = conditional_equalities_.find({l.NegatedIndex(), p});
if (it != conditional_equalities_.end()) {
for (const IntegerVariable x : it->second) {
ProcessNewProduct(p, l.Negated(), x);
}
}
}
}
} // namespace sat
} // namespace operations_research

117
ortools/sat/implied_bounds.h
View File

@@ -15,6 +15,8 @@
#define OR_TOOLS_SAT_IMPLIED_BOUNDS_H_
#include <algorithm>
#include <array>
#include <bitset>
#include <cstdint>
#include <utility>
#include <vector>
@@ -29,6 +31,7 @@
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/sat/sat_solver.h"
#include "ortools/sat/synchronization.h"
#include "ortools/util/bitset.h"
#include "ortools/util/strong_integers.h"
@@ -85,7 +88,8 @@ class ImpliedBounds {
: parameters_(*model->GetOrCreate<SatParameters>()),
sat_solver_(model->GetOrCreate<SatSolver>()),
integer_trail_(model->GetOrCreate<IntegerTrail>()),
integer_encoder_(model->GetOrCreate<IntegerEncoder>()) {}
integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
shared_stats_(model->GetOrCreate<SharedStatistics>()) {}
~ImpliedBounds();
// Adds literal => integer_literal to the repository.
@@ -157,6 +161,7 @@ class ImpliedBounds {
SatSolver* sat_solver_;
IntegerTrail* integer_trail_;
IntegerEncoder* integer_encoder_;
SharedStatistics* shared_stats_;
// TODO(user): Remove the need for this.
std::vector<IntegerLiteral> tmp_integer_literals_;
@@ -225,6 +230,116 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left,
const AffineExpression& right, Model* model,
LinearConstraintBuilder* builder);
// Class used to detect and hold all the information about a variable beeing the
// product of two others. This class is meant to be used by LP relaxation and
// cuts.
class ProductDetector {
public:
explicit ProductDetector(Model* model);
~ProductDetector();
// Internally, a Boolean product is encoded in a linear fashion:
// p = a * b become:
// 1/ a and b => p, i.e. a clause (not(a), not(b), p).
// 2/ p => a and p => b, which is a clause (not(p), a) and (not(p), b).
//
// In particular if we have a+b+c==1 then we have a=b*c, b=a*c, and c=a*b !!
//
// For the detection to work, we must load all ternary clause first, then the
// implication.
void ProcessTernaryClause(absl::Span<const Literal> ternary_clause);
void ProcessTernaryExactlyOne(absl::Span<const Literal> ternary_exo);
void ProcessBinaryClause(absl::Span<const Literal> binary_clause);
// Utility function to process a bunch of implication all at once.
void ProcessImplicationGraph(BinaryImplicationGraph* graph);
void ProcessTrailAtLevelOne();
// We also detect product of Boolean with IntegerVariable.
// After presolve, a product P = l * X should be encoded with:
// l => P = X
// not(l) => P = 0
//
// TODO(user): Generalize to a * X + b = l * (Y + c) since these are also
// easy to linearize if we see l * Y.
void ProcessConditionalEquality(Literal l, IntegerVariable x,
IntegerVariable y);
void ProcessConditionalZero(Literal l, IntegerVariable p);
// Query function mainly used for testing.
LiteralIndex GetProduct(Literal a, Literal b) const;
IntegerVariable GetProduct(Literal a, IntegerVariable b) const;
// Query Functions. LinearizeProduct() should only be called if
// ProductIsLinearizable() is true.
bool ProductIsLinearizable(IntegerVariable a, IntegerVariable b) const;
// TODO(user): Implement!
LinearExpression LinearizeProduct(IntegerVariable a, IntegerVariable b);
// Returns an expression that is always lower or equal to the product a * b.
// This use the exact LinearizeProduct() if ProductIsLinearizable() otherwise
// it uses the simple McCormick lower bound.
//
// TODO(user): Implement!
LinearExpression ProductLowerBound(IntegerVariable a, IntegerVariable b);
private:
std::array<LiteralIndex, 2> GetKey(LiteralIndex a, LiteralIndex b) const;
void ProcessNewProduct(LiteralIndex p, LiteralIndex a, LiteralIndex b);
void ProcessNewProduct(IntegerVariable p, Literal l, IntegerVariable x);
// Fixed at creation time.
bool enabled_;
SatSolver* sat_solver_;
Trail* trail_;
IntegerTrail* integer_trail_;
IntegerEncoder* integer_encoder_;
SharedStatistics* shared_stats_;
// No need to process implication a => b if a was never seen.
absl::StrongVector<LiteralIndex, bool> seen_;
// For each clause of size 3 (l0, l1, l2) and a permutation of index (i, j, k)
// we bitset[i] to true if lj => not(lk) and lk => not(lj).
//
// The key is sorted.
absl::flat_hash_map<std::array<LiteralIndex, 3>, std::bitset<3>> detector_;
// For each (l0, l1) we list all the l2 such that (l0, l1, l2) is a 3 clause.
absl::flat_hash_map<std::array<LiteralIndex, 2>, std::vector<LiteralIndex>>
candidates_;
// Products (a, b) -> p such that p == a * b. They key is sorted.
absl::flat_hash_map<std::array<LiteralIndex, 2>, LiteralIndex> products_;
// Same keys has in products_ but canonicalized so we capture all 4 products
// a * b, (1 - a) * b, a * (1 - b) and (1 - a) * (1 - b) with one query.
absl::flat_hash_set<std::array<LiteralIndex, 2>> has_product_;
// For bool * int detection. Note that we only use positive IntegerVariable
// in the key part.
absl::flat_hash_set<std::pair<LiteralIndex, IntegerVariable>>
conditional_zeros_;
absl::flat_hash_map<std::pair<LiteralIndex, IntegerVariable>,
std::vector<IntegerVariable>>
conditional_equalities_;
// Stores l * X = P.
absl::flat_hash_map<std::pair<LiteralIndex, IntegerVariable>, IntegerVariable>
int_products_;
// Stats.
int64_t num_products_ = 0;
int64_t num_int_products_ = 0;
int64_t num_trail_updates_ = 0;
int64_t num_processed_binary_ = 0;
int64_t num_processed_ternary_ = 0;
int64_t num_processed_exo_ = 0;
int64_t num_conditional_zeros_ = 0;
int64_t num_conditional_equalities_ = 0;
};
} // namespace sat
} // namespace operations_research

2
ortools/sat/integer_search.cc
View File

@@ -791,6 +791,7 @@ IntegerSearchHelper::IntegerSearchHelper(Model* model)
integer_trail_(model->GetOrCreate<IntegerTrail>()),
encoder_(model->GetOrCreate<IntegerEncoder>()),
implied_bounds_(model->GetOrCreate<ImpliedBounds>()),
product_detector_(model->GetOrCreate<ProductDetector>()),
time_limit_(model->GetOrCreate<TimeLimit>()),
pseudo_costs_(model->GetOrCreate<PseudoCosts>()) {
// This is needed for recording the pseudo-costs.
@@ -907,6 +908,7 @@ bool IntegerSearchHelper::TakeDecision(Literal decision) {
// This is "almost free", so we might as well do it.
if (old_level == 0 && sat_solver_->CurrentDecisionLevel() == 1) {
implied_bounds_->ProcessIntegerTrail(decision);
product_detector_->ProcessTrailAtLevelOne();
}
// Update the pseudo costs.

1
ortools/sat/integer_search.h
View File

@@ -273,6 +273,7 @@ class IntegerSearchHelper {
IntegerTrail* integer_trail_;
IntegerEncoder* encoder_;
ImpliedBounds* implied_bounds_;
ProductDetector* product_detector_;
TimeLimit* time_limit_;
PseudoCosts* pseudo_costs_;
IntegerVariable objective_var_ = kNoIntegerVariable;

5
ortools/sat/linear_propagation.h
View File

@@ -195,10 +195,11 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface {
int rev_size; // The size of the non-fixed terms.
IntegerValue rev_rhs; // The current rhs, updated on fixed terms.
};
#if !defined(_MSC_VER)
#if !defined(_MSC_VER)
static_assert(sizeof(ConstraintInfo) == 24,
"ERROR_ConstraintInfo_is_not_well_compacted");
#endif
#endif // !defined(_MSC_VER)
absl::Span<IntegerValue> GetCoeffs(const ConstraintInfo& info);
absl::Span<IntegerVariable> GetVariables(const ConstraintInfo& info);

2
ortools/sat/probing.cc
View File

@@ -45,6 +45,7 @@ Prober::Prober(Model* model)
assignment_(model->GetOrCreate<SatSolver>()->Assignment()),
integer_trail_(model->GetOrCreate<IntegerTrail>()),
implied_bounds_(model->GetOrCreate<ImpliedBounds>()),
product_detector_(model->GetOrCreate<ProductDetector>()),
sat_solver_(model->GetOrCreate<SatSolver>()),
time_limit_(model->GetOrCreate<TimeLimit>()),
implication_graph_(model->GetOrCreate<BinaryImplicationGraph>()),
@@ -83,6 +84,7 @@ bool Prober::ProbeOneVariableInternal(BooleanVariable b) {
}
implied_bounds_->ProcessIntegerTrail(decision);
product_detector_->ProcessTrailAtLevelOne();
integer_trail_->AppendNewBounds(&new_integer_bounds_);
for (int i = saved_index + 1; i < trail_.Index(); ++i) {
const Literal l = trail_[i];

1
ortools/sat/probing.h
View File

@@ -102,6 +102,7 @@ class Prober {
const VariablesAssignment& assignment_;
IntegerTrail* integer_trail_;
ImpliedBounds* implied_bounds_;
ProductDetector* product_detector_;
SatSolver* sat_solver_;
TimeLimit* time_limit_;
BinaryImplicationGraph* implication_graph_;

4
ortools/sat/sat_parameters.proto
View File

@@ -1119,10 +1119,6 @@ message SatParameters {
}
optional FPRoundingMethod fp_rounding = 165 [default = PROPAGATION_ASSISTED];
// Turns on a lns worker which solves relaxed version of the original problem
// by removing constraints from the problem in order to get better bounds.
optional bool use_relaxation_lns = 150 [default = false];
// If true, registers more lns subsolvers with different parameters.
optional bool diversify_lns_params = 137 [default = false];

47
ortools/sat/synchronization.cc
View File

@@ -64,21 +64,22 @@ namespace operations_research {
namespace sat {
void SharedRelaxationSolutionRepository::NewRelaxationSolution(
const CpSolverResponse& response) {
absl::Span<const int64_t> solution_values,
IntegerValue inner_objective_value) {
// Note that the Add() method already applies mutex lock. So we don't need it
// here.
if (response.solution().empty()) return;
if (solution_values.empty()) return;
// Add this solution to the pool.
SharedSolutionRepository<int64_t>::Solution solution;
solution.variable_values.assign(response.solution().begin(),
response.solution().end());
solution.variable_values.assign(solution_values.begin(),
solution_values.end());
// For now we use the negated lower bound as the "internal objective" to
// prefer solution with an higher bound.
//
// Note: If the model doesn't have objective, the best_objective_bound is set
// to default value 0.
solution.rank = -response.best_objective_bound();
solution.rank = -inner_objective_value.value();
Add(solution);
}
@@ -507,28 +508,29 @@ void SharedResponseManager::FillObjectiveValuesInBestResponse() {
best_response_.set_gap_integral(gap_integral_);
}
void SharedResponseManager::NewSolution(const CpSolverResponse& response,
Model* model) {
void SharedResponseManager::NewSolution(
absl::Span<const int64_t> solution_values, const std::string& solution_info,
Model* model) {
absl::MutexLock mutex_lock(&mutex_);
// Special case if the user asked to keep solutions in the pool.
if (objective_or_null_ == nullptr && parameters_.enumerate_all_solutions() &&
parameters_.fill_additional_solutions_in_response()) {
SharedSolutionRepository<int64_t>::Solution solution;
solution.variable_values.assign(response.solution().begin(),
response.solution().end());
solution.variable_values.assign(solution_values.begin(),
solution_values.end());
solutions_.Add(solution);
}
if (objective_or_null_ != nullptr) {
const int64_t objective_value =
ComputeInnerObjective(*objective_or_null_, response);
ComputeInnerObjective(*objective_or_null_, solution_values);
// Add this solution to the pool, even if it is not improving.
if (!response.solution().empty()) {
if (!solution_values.empty()) {
SharedSolutionRepository<int64_t>::Solution solution;
solution.variable_values.assign(response.solution().begin(),
response.solution().end());
solution.variable_values.assign(solution_values.begin(),
solution_values.end());
solution.rank = objective_value;
solutions_.Add(solution);
}
@@ -563,8 +565,9 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response,
best_response_.set_status(CpSolverStatus::FEASIBLE);
}
best_response_.set_solution_info(response.solution_info());
*best_response_.mutable_solution() = response.solution();
best_response_.set_solution_info(solution_info);
best_response_.mutable_solution()->Assign(solution_values.begin(),
solution_values.end());
// Mark model as OPTIMAL if the inner bound crossed.
if (objective_or_null_ != nullptr &&
@@ -574,12 +577,13 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response,
// Logging.
++num_solutions_;
// TODO(user): Remove this code and the need for model in this function.
if (logger_->LoggingIsEnabled()) {
std::string solution_info = response.solution_info();
std::string solution_message = solution_info;
if (model != nullptr) {
const int64_t num_bool = model->Get<Trail>()->NumVariables();
const int64_t num_fixed = model->Get<SatSolver>()->NumFixedVariables();
absl::StrAppend(&solution_info, " fixed_bools:", num_fixed, "/",
absl::StrAppend(&solution_message, " fixed_bools:", num_fixed, "/",
num_bool);
}
@@ -592,13 +596,14 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response,
if (obj.scaling_factor() < 0) {
std::swap(lb, ub);
}
RegisterSolutionFound(solution_info);
RegisterSolutionFound(solution_message);
SOLVER_LOG(logger_, ProgressMessage(absl::StrCat(num_solutions_),
wall_timer_.Get(), best, lb, ub,
solution_info));
solution_message));
} else {
SOLVER_LOG(logger_, SatProgressMessage(absl::StrCat(num_solutions_),
wall_timer_.Get(), solution_info));
SOLVER_LOG(logger_,
SatProgressMessage(absl::StrCat(num_solutions_),
wall_timer_.Get(), solution_message));
}
}

11
ortools/sat/synchronization.h
View File

@@ -138,7 +138,8 @@ class SharedRelaxationSolutionRepository
explicit SharedRelaxationSolutionRepository(int num_solutions_to_keep)
: SharedSolutionRepository<int64_t>(num_solutions_to_keep) {}
void NewRelaxationSolution(const CpSolverResponse& response);
void NewRelaxationSolution(absl::Span<const int64_t> solution_values,
IntegerValue inner_objective_value);
};
class SharedLPSolutionRepository : public SharedSolutionRepository<double> {
@@ -291,12 +292,8 @@ class SharedResponseManager {
// Reads the new solution from the response and update our state. For an
// optimization problem, we only do something if the solution is strictly
// improving.
//
// TODO(user): Only the following fields from response are accessed here, we
// might want a tighter API:
// - solution_info
// - solution
void NewSolution(const CpSolverResponse& response, Model* model);
void NewSolution(absl::Span<const int64_t> solution_values,
const std::string& solution_info, Model* model = nullptr);
// Changes the solution to reflect the fact that the "improving" problem is
// infeasible. This means that if we have a solution, we have proven