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

add deterministic time limit to the graph symmetry finder; improve support for scheduling constraints in regular LNS for CP-SAT; improve symmetry detection for CP-SAT

Browse Source
This commit is contained in:
Laurent Perron
2021-02-03 15:39:48 +01:00
parent 780fc717bf
commit 111f43b139
9 changed files with 448 additions and 276 deletions
Download Patch File Download Diff File Expand all files Collapse all files

53
ortools/algorithms/find_graph_symmetries.cc
View File

@@ -105,7 +105,7 @@ GraphSymmetryFinder::GraphSymmetryFinder(const Graph& graph, bool is_undirected)
tmp_degree_(NumNodes(), 0),
tmp_nodes_with_degree_(NumNodes() + 1) {
// Set up an "unlimited" time limit by default.
time_limit_ = TimeLimit::Infinite();
time_limit_ = &dummy_time_limit_;
tmp_partition_.Reset(NumNodes());
if (is_undirected) {
DCHECK(GraphIsSymmetric(graph));
@@ -178,7 +178,9 @@ template <class T>
inline void IncrementCounterForNonSingletons(const T& nodes,
const DynamicPartition& partition,
std::vector<int>* node_count,
std::vector<int>* nodes_seen) {
std::vector<int>* nodes_seen,
int64* num_operations) {
*num_operations += nodes.end() - nodes.begin();
for (const int node : nodes) {
if (partition.ElementsInSamePartAs(node).size() == 1) continue;
const int count = ++(*node_count)[node];
@@ -192,6 +194,15 @@ void GraphSymmetryFinder::RecursivelyRefinePartitionByAdjacency(
// Rename, for readability of the code below.
std::vector<int>& tmp_nodes_with_nonzero_degree = tmp_stack_;
// This function is the main bottleneck of the whole algorithm. We count the
// number of blocks in the inner-most loops in num_operations. At the end we
// will multiply it by a factor to have some deterministic time that we will
// append to the deterministic time counter.
//
// TODO(user): We are really imprecise in our counting, but it is fine. We
// just need a way to enforce a deterministic limit on the computation effort.
int64 num_operations = 0;
// Assuming that the partition was refined based on the adjacency on
// parts [0 .. first_unrefined_part_index) already, we simply need to
// refine parts first_unrefined_part_index ... NumParts()-1, the latter bound
@@ -215,19 +226,20 @@ void GraphSymmetryFinder::RecursivelyRefinePartitionByAdjacency(
// come from/to the current part.
if (outgoing_adjacency) {
for (const int node : partition->ElementsInPart(part_index)) {
IncrementCounterForNonSingletons(graph_[node], *partition,
&tmp_degree_,
&tmp_nodes_with_nonzero_degree);
IncrementCounterForNonSingletons(
graph_[node], *partition, &tmp_degree_,
&tmp_nodes_with_nonzero_degree, &num_operations);
}
} else {
for (const int node : partition->ElementsInPart(part_index)) {
IncrementCounterForNonSingletons(TailsOfIncomingArcsTo(node),
*partition, &tmp_degree_,
&tmp_nodes_with_nonzero_degree);
IncrementCounterForNonSingletons(
TailsOfIncomingArcsTo(node), *partition, &tmp_degree_,
&tmp_nodes_with_nonzero_degree, &num_operations);
}
}
// Group the nodes by (nonzero) degree. Remember the maximum degree.
int max_degree = 0;
num_operations += 3 + tmp_nodes_with_nonzero_degree.size();
for (const int node : tmp_nodes_with_nonzero_degree) {
const int degree = tmp_degree_[node];
tmp_degree_[node] = 0; // To clean up after us.
@@ -238,11 +250,20 @@ void GraphSymmetryFinder::RecursivelyRefinePartitionByAdjacency(
// For each degree, refine the partition by the set of nodes with that
// degree.
for (int degree = 1; degree <= max_degree; ++degree) {
// We use a manually tuned factor 3 because Refine() does quite a bit of
// operations for each node in its argument.
num_operations += 1 + 3 * tmp_nodes_with_degree_[degree].size();
partition->Refine(tmp_nodes_with_degree_[degree]);
tmp_nodes_with_degree_[degree].clear(); // To clean up after us.
}
}
}
// The coefficient was manually tuned (only on a few instances) so that the
// time is roughly correlated with seconds on a fast desktop computer from
// 2020.
time_limit_->AdvanceDeterministicTime(1e-8 *
static_cast<double>(num_operations));
}
void GraphSymmetryFinder::DistinguishNodeInPartition(
@@ -353,11 +374,12 @@ void GetAllOtherRepresentativesInSamePartAs(
} // namespace
absl::Status GraphSymmetryFinder::FindSymmetries(
double time_limit_seconds, std::vector<int>* node_equivalence_classes_io,
std::vector<int>* node_equivalence_classes_io,
std::vector<std::unique_ptr<SparsePermutation>>* generators,
std::vector<int>* factorized_automorphism_group_size) {
std::vector<int>* factorized_automorphism_group_size,
TimeLimit* time_limit) {
// Initialization.
time_limit_ = absl::make_unique<TimeLimit>(time_limit_seconds);
time_limit_ = time_limit == nullptr ? &dummy_time_limit_ : time_limit;
IF_STATS_ENABLED(stats_.initialization_time.StartTimer());
generators->clear();
factorized_automorphism_group_size->clear();
@@ -744,8 +766,10 @@ GraphSymmetryFinder::FindOneSuitablePermutation(
tmp_dynamic_permutation_.UndoLastMappings(&base_singletons);
} else {
ScopedTimeDistributionUpdater u(&stats_.map_reelection_time);
// TODO(user): try to get the non-singleton part from
// DynamicPermutation in O(1), if it matters to the overall speed.
// TODO(user, viger): try to get the non-singleton part from
// DynamicPermutation in O(1). On some graphs like the symmetry of the
// mip problem lectsched-4-obj.mps.gz, this take the majority of the
// time!
int non_singleton_part = 0;
{
ScopedTimeDistributionUpdater u(&stats_.non_singleton_search_time);
@@ -754,6 +778,9 @@ GraphSymmetryFinder::FindOneSuitablePermutation(
DCHECK_LT(non_singleton_part, base_partition->NumParts());
}
}
time_limit_->AdvanceDeterministicTime(
1e-9 * static_cast<double>(non_singleton_part));
// The partitions are compatible, but we'll deepen the search on some
// non-singleton part. We can pick any base and image node in this case.
GetBestMapping(*base_partition, *image_partition, non_singleton_part,

17
ortools/algorithms/find_graph_symmetries.h
View File

@@ -92,8 +92,9 @@ class GraphSymmetryFinder {
// large as N!).
//
// DEADLINE AND PARTIAL COMPLETION:
// If the deadline passed as argument is reached, this method will return
// quickly (within a few milliseconds). The outputs may be partially filled:
// If the deadline passed as argument (via TimeLimit) is reached, this method
// will return quickly (within a few milliseconds of the limit). The outputs
// may be partially filled:
// - Each element of "generators", if non-empty, will be a valid permutation.
// - "node_equivalence_classes_io" will contain the equivalence classes
// corresponding to the orbits under all the generators in "generators".
@@ -101,9 +102,10 @@ class GraphSymmetryFinder {
// partially valid: its last element may be undervalued. But all prior
// elements are valid factors of the automorphism group size.
absl::Status FindSymmetries(
double time_limit_seconds, std::vector<int>* node_equivalence_classes_io,
std::vector<int>* node_equivalence_classes_io,
std::vector<std::unique_ptr<SparsePermutation> >* generators,
std::vector<int>* factorized_automorphism_group_size);
std::vector<int>* factorized_automorphism_group_size,
TimeLimit* time_limit = nullptr);
// Fully refine the partition of nodes, using the graph as symmetry breaker.
// This means applying the following steps on each part P of the partition:
@@ -148,8 +150,11 @@ class GraphSymmetryFinder {
util::BeginEndWrapper<std::vector<int>::const_iterator> TailsOfIncomingArcsTo(
int node) const;
// Deadline management. Populated upon FindSymmetries().
mutable std::unique_ptr<TimeLimit> time_limit_;
// Deadline management. Populated upon FindSymmetries(). If the passed
// time limit is nullptr, time_limit_ will point to dummy_time_limit_ which
// is an object with infinite limits by default.
TimeLimit dummy_time_limit_;
TimeLimit* time_limit_;
// Internal search code used in FindSymmetries(), split out for readability:
// find one permutation (if it exists) that maps root_node to root_image_node

5
ortools/sat/boolean_problem.cc
View File

@@ -705,9 +705,8 @@ void FindLinearBooleanProblemSymmetries(
/*is_undirected=*/true);
std::vector<int> factorized_automorphism_group_size;
// TODO(user): inject the appropriate time limit here.
CHECK_OK(symmetry_finder.FindSymmetries(
/*time_limit_seconds=*/std::numeric_limits<double>::infinity(),
&equivalence_classes, generators, &factorized_automorphism_group_size));
CHECK_OK(symmetry_finder.FindSymmetries(&equivalence_classes, generators,
&factorized_automorphism_group_size));
// Remove from the permutations the part not concerning the literals.
// Note that some permutation may becomes empty, which means that we had

130
ortools/sat/cp_model_lns.cc
View File

@@ -111,14 +111,30 @@ void NeighborhoodGeneratorHelper::RecomputeHelperData() {
// presolve?
var_to_constraint_.assign(model_proto_.variables_size(), {});
constraint_to_var_.assign(model_proto_.constraints_size(), {});
const auto register_var = [&](int var, int ct_index) {
DCHECK(RefIsPositive(var));
if (IsConstant(var)) return;
var_to_constraint_[var].push_back(ct_index);
constraint_to_var_[ct_index].push_back(var);
CHECK_GE(var, 0);
CHECK_LT(var, model_proto_.variables_size());
};
for (int ct_index = 0; ct_index < model_proto_.constraints_size();
++ct_index) {
for (const int var : UsedVariables(model_proto_.constraints(ct_index))) {
if (IsConstant(var)) continue;
var_to_constraint_[var].push_back(ct_index);
constraint_to_var_[ct_index].push_back(var);
CHECK_GE(var, 0);
CHECK_LT(var, model_proto_.variables_size());
register_var(var, ct_index);
}
// We replace intervals by their underlying integer variables.
if (parameters_.lns_expand_intervals_in_constraint_graph()) {
for (const int interval :
UsedIntervals(model_proto_.constraints(ct_index))) {
for (const int var :
UsedVariables(model_proto_.constraints(interval))) {
register_var(var, ct_index);
}
}
}
}
@@ -179,6 +195,37 @@ Neighborhood NeighborhoodGeneratorHelper::FullNeighborhood() const {
return neighborhood;
}
std::vector<int> NeighborhoodGeneratorHelper::GetActiveIntervals(
const CpSolverResponse& initial_solution) const {
std::vector<int> active_intervals;
for (const int i : TypeToConstraints(ConstraintProto::kInterval)) {
const ConstraintProto& interval_ct = ModelProto().constraints(i);
// We only look at intervals that are performed in the solution. The
// unperformed intervals should be automatically freed during the generation
// phase.
if (interval_ct.enforcement_literal().size() == 1) {
const int enforcement_ref = interval_ct.enforcement_literal(0);
const int enforcement_var = PositiveRef(enforcement_ref);
const int value = initial_solution.solution(enforcement_var);
if (RefIsPositive(enforcement_ref) == (value == 0)) {
continue;
}
}
// We filter out fixed intervals. Because of presolve, if there is an
// enforcement literal, it cannot be fixed.
if (interval_ct.enforcement_literal().empty() &&
IsConstant(PositiveRef(interval_ct.interval().start())) &&
IsConstant(PositiveRef(interval_ct.interval().size())) &&
IsConstant(PositiveRef(interval_ct.interval().end()))) {
continue;
}
active_intervals.push_back(i);
}
return active_intervals;
}
Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables(
const CpSolverResponse& initial_solution,
const std::vector<int>& variables_to_fix) const {
@@ -363,6 +410,47 @@ Neighborhood SimpleNeighborhoodGenerator::Generate(
return helper_.FixGivenVariables(initial_solution, fixed_variables);
}
Neighborhood SimpleConstraintNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, double difficulty,
absl::BitGenRef random) {
std::vector<int> active_constraints;
for (int ct = 0; ct < helper_.ModelProto().constraints_size(); ++ct) {
if (helper_.ModelProto().constraints(ct).constraint_case() ==
ConstraintProto::CONSTRAINT_NOT_SET) {
continue;
}
active_constraints.push_back(ct);
}
const int num_active_vars = helper_.ActiveVariables().size();
const int num_model_vars = helper_.ModelProto().variables_size();
const int target_size = std::ceil(difficulty * num_active_vars);
const int num_constraints = helper_.ConstraintToVar().size();
if (num_constraints == 0 || target_size == num_active_vars) {
return helper_.FullNeighborhood();
}
CHECK_GT(target_size, 0);
std::shuffle(active_constraints.begin(), active_constraints.end(), random);
std::vector<bool> visited_variables_set(num_model_vars, false);
std::vector<int> relaxed_variables;
for (const int constraint_index : active_constraints) {
CHECK_LT(constraint_index, num_constraints);
for (const int var : helper_.ConstraintToVar()[constraint_index]) {
if (visited_variables_set[var]) continue;
visited_variables_set[var] = true;
if (helper_.IsActive(var)) {
relaxed_variables.push_back(var);
if (relaxed_variables.size() == target_size) break;
}
}
if (relaxed_variables.size() == target_size) break;
}
return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
}
Neighborhood VariableGraphNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, double difficulty,
absl::BitGenRef random) {
@@ -442,14 +530,14 @@ Neighborhood ConstraintGraphNeighborhoodGenerator::Generate(
// Pick a random unprocessed constraint.
const int i = absl::Uniform<int>(random, 0, next_constraints.size());
const int contraint_index = next_constraints[i];
const int constraint_index = next_constraints[i];
std::swap(next_constraints[i], next_constraints.back());
next_constraints.pop_back();
// Add all the variable of this constraint and increase the set of next
// possible constraints.
CHECK_LT(contraint_index, num_constraints);
random_variables = helper_.ConstraintToVar()[contraint_index];
CHECK_LT(constraint_index, num_constraints);
random_variables = helper_.ConstraintToVar()[constraint_index];
std::shuffle(random_variables.begin(), random_variables.end(), random);
for (const int var : random_variables) {
if (visited_variables_set[var]) continue;
@@ -494,20 +582,25 @@ Neighborhood GenerateSchedulingNeighborhoodForRelaxation(
const int enforcement_var = PositiveRef(enforcement_ref);
const int value = initial_solution.solution(enforcement_var);
// If the interval is not enforced, we just relax it. If it belongs to an
// exactly one constraint, and the enforced interval is not relaxed, then
// propagation will force this interval to stay not enforced. Otherwise,
// LNS will be able to change which interval will be enforced among all
// alternatives.
if (RefIsPositive(enforcement_ref) == (value == 0)) {
ignored_intervals.insert(i);
continue;
}
// Fix the value.
neighborhood.cp_model.mutable_variables(enforcement_var)->clear_domain();
neighborhood.cp_model.mutable_variables(enforcement_var)->add_domain(value);
neighborhood.cp_model.mutable_variables(enforcement_var)->add_domain(value);
// If the interval is ignored, skip for the loop below as there is no
// point adding precedence on it.
if (RefIsPositive(enforcement_ref) == (value == 0)) {
ignored_intervals.insert(i);
}
}
for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap)) {
// Sort all non-relaxed intervals of this constraint by current start time.
// Sort all non-relaxed intervals of this constraint by current start
// time.
std::vector<std::pair<int64, int>> start_interval_pairs;
for (const int i :
neighborhood.cp_model.constraints(c).no_overlap().intervals()) {
@@ -565,8 +658,8 @@ Neighborhood GenerateSchedulingNeighborhoodForRelaxation(
Neighborhood SchedulingNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, double difficulty,
absl::BitGenRef random) {
const auto span = helper_.TypeToConstraints(ConstraintProto::kInterval);
std::vector<int> intervals_to_relax(span.begin(), span.end());
std::vector<int> intervals_to_relax =
helper_.GetActiveIntervals(initial_solution);
GetRandomSubset(difficulty, &intervals_to_relax, random);
return GenerateSchedulingNeighborhoodForRelaxation(intervals_to_relax,
@@ -577,9 +670,8 @@ Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, double difficulty,
absl::BitGenRef random) {
std::vector<std::pair<int64, int>> start_interval_pairs;
for (const int i : helper_.TypeToConstraints(ConstraintProto::kInterval)) {
for (const int i : helper_.GetActiveIntervals(initial_solution)) {
const ConstraintProto& interval_ct = helper_.ModelProto().constraints(i);
const int start_var = interval_ct.interval().start();
const int64 start_value = initial_solution.solution(start_var);
start_interval_pairs.push_back({start_value, i});

21
ortools/sat/cp_model_lns.h
View File

@@ -121,6 +121,13 @@ class NeighborhoodGeneratorHelper : public SubSolver {
return absl::MakeSpan(type_to_constraints_[type]);
}
// Returns the list of indices of active interval constraints according
// to the initial_solution and the parameter lns_focus_on_performed_intervals.
// If true, this method returns the list of performed intervals in the
// solution. If false, it returns all intervals of the model.
std::vector<int> GetActiveIntervals(
const CpSolverResponse& initial_solution) const;
// The initial problem.
// Note that the domain of the variables are not updated here.
const CpModelProto& ModelProto() const { return model_proto_; }
@@ -348,6 +355,18 @@ class SimpleNeighborhoodGenerator : public NeighborhoodGenerator {
double difficulty, absl::BitGenRef random) final;
};
// Pick a random subset of constraints and relax all the variables of these
// constraints. Note that to satisfy the difficulty, we might not relax all the
// variable of the "last" constraint.
class SimpleConstraintNeighborhoodGenerator : public NeighborhoodGenerator {
public:
explicit SimpleConstraintNeighborhoodGenerator(
NeighborhoodGeneratorHelper const* helper, const std::string& name)
: NeighborhoodGenerator(name, helper) {}
Neighborhood Generate(const CpSolverResponse& initial_solution,
double difficulty, absl::BitGenRef random) final;
};
// Pick a random subset of variables that are constructed by a BFS in the
// variable <-> constraint graph. That is, pick a random variable, then all the
// variable connected by some constraint to the first one, and so on. The
@@ -362,7 +381,7 @@ class VariableGraphNeighborhoodGenerator : public NeighborhoodGenerator {
};
// Pick a random subset of constraint and relax all of their variables. We are a
// bit smarter than this because after the first contraint is selected, we only
// bit smarter than this because after the first constraint is selected, we only
// select constraints that share at least one variable with the already selected
// constraints. The variable from the "last" constraint are selected randomly.
class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator {

19
ortools/sat/cp_model_solver.cc
View File

@@ -1197,7 +1197,14 @@ void LoadBaseModel(const CpModelProto& model_proto,
model_proto.search_strategy().empty());
mapping->CreateVariables(model_proto, view_all_booleans_as_integers, model);
mapping->DetectOptionalVariables(model_proto, model);
mapping->DetectAndLoadBooleanSymmetries(model_proto, model);
// TODO(user): The core algo and symmetries seems to be problematic in some
// cases. See for instance: neos-691058.mps.gz. This is probably because as
// we modify the model, our symmetry might be wrong? investigate.
if (!parameters.optimize_with_core()) {
mapping->DetectAndLoadBooleanSymmetries(model_proto, model);
}
mapping->ExtractEncoding(model_proto, model);
mapping->PropagateEncodingFromEquivalenceRelations(model_proto, model);
@@ -2760,15 +2767,19 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
// Each will have their own metrics.
subsolvers.push_back(absl::make_unique<LnsSolver>(
absl::make_unique<SimpleNeighborhoodGenerator>(
helper, absl::StrCat("rnd_lns_", local_params.name())),
helper, absl::StrCat("rnd_var_lns_", local_params.name())),
local_params, helper, &shared));
subsolvers.push_back(absl::make_unique<LnsSolver>(
absl::make_unique<SimpleConstraintNeighborhoodGenerator>(
helper, absl::StrCat("rnd_cst_lns_", local_params.name())),
local_params, helper, &shared));
subsolvers.push_back(absl::make_unique<LnsSolver>(
absl::make_unique<VariableGraphNeighborhoodGenerator>(
helper, absl::StrCat("var_lns_", local_params.name())),
helper, absl::StrCat("graph_var_lns_", local_params.name())),
local_params, helper, &shared));
subsolvers.push_back(absl::make_unique<LnsSolver>(
absl::make_unique<ConstraintGraphNeighborhoodGenerator>(
helper, absl::StrCat("cst_lns_", local_params.name())),
helper, absl::StrCat("graph_cst_lns_", local_params.name())),
local_params, helper, &shared));
if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap).empty()) {

468
ortools/sat/cp_model_symmetries.cc
View File

@@ -44,10 +44,10 @@ class IdGenerator {
public:
IdGenerator() {}
// If the key was never seen before, then generate a new id, otherwise return
// the previously generated id.
int GetId(const std::vector<int64>& key) {
return gtl::LookupOrInsert(&id_map_, key, id_map_.size());
// If the color was never seen before, then generate a new id, otherwise
// return the previously generated id.
int GetId(const std::vector<int64>& color) {
return gtl::LookupOrInsert(&id_map_, color, id_map_.size());
}
int NextFreeId() const { return id_map_.size(); }
@@ -91,173 +91,177 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
const int num_variables = problem.variables_size();
auto graph = absl::make_unique<Graph>();
initial_equivalence_classes->clear();
// Each node will be created with a given color. Two nodes of different color
// can never be send one into another by a symmetry. The first element of
// the color vector will always be the NodeType.
//
// TODO(user): Using a full int64 for storing 3 values is not great. We
// can optimize this at the price of a bit more code.
enum NodeType {
VARIABLE_NODE,
VAR_COEFFICIENT_NODE,
CONSTRAINT_NODE,
CONSTRAINT_COEFFICIENT_NODE,
ENFORCEMENT_LITERAL
};
IdGenerator id_generator;
IdGenerator color_id_generator;
initial_equivalence_classes->clear();
auto new_node = [&initial_equivalence_classes,
&color_id_generator](const std::vector<int64>& color) {
// Since we add nodes one by one, initial_equivalence_classes->size() gives
// the number of nodes at any point, which we use as the next node index.
const int node = initial_equivalence_classes->size();
initial_equivalence_classes->push_back(color_id_generator.GetId(color));
return node;
};
// For two variables to be in the same equivalence class, they need to have
// the same objective coefficient, and the same possible bounds.
//
// TODO(user): We could ignore the objective coefficients, and just make sure
// that when we break symmetry amongst variables, we choose the possibility
// with the smallest cost?
std::vector<int64> objective_by_var(num_variables, 0);
for (int i = 0; i < problem.objective().vars_size(); ++i) {
objective_by_var[problem.objective().vars(i)] =
problem.objective().coeffs(i);
}
auto new_node = [&initial_equivalence_classes,
&id_generator](const std::vector<int64>& key) {
// Since we add nodes one by one, initial_equivalence_classes->size() gives
// the number of nodes at any point, which we use as the next node index.
const int node = initial_equivalence_classes->size();
initial_equivalence_classes->push_back(id_generator.GetId(key));
return node;
};
std::vector<int64> tmp_key;
// Create one node for each variable. Note that the code rely on the fact that
// the index of a VARIABLE_NODE type is the same as the variable index.
std::vector<int64> tmp_color;
for (int v = 0; v < num_variables; ++v) {
tmp_key = {VARIABLE_NODE, objective_by_var[v]};
Append(problem.variables(v).domain(), &tmp_key);
// Note that the code rely on the fact that the index of a VARIABLE_NODE
// type is the same as the variable index.
CHECK_EQ(v, new_node(tmp_key));
tmp_color = {VARIABLE_NODE, objective_by_var[v]};
Append(problem.variables(v).domain(), &tmp_color);
CHECK_EQ(v, new_node(tmp_color));
// Make sure the graph contains all the variable nodes, even if no edges are
// attached to them through constraints.
graph->AddNode(v);
}
auto add_edge = [&graph](int node_1, int node_2) {
graph->AddArc(node_1, node_2);
graph->AddArc(node_2, node_1);
};
// We will lazily create "coefficient nodes" that correspond to a variable
// with a given coefficient.
absl::flat_hash_map<std::pair<int64, int64>, int> coefficient_nodes;
auto get_coefficient_node = [&new_node, &graph, &coefficient_nodes,
&tmp_color](int var, int64 coeff) {
const int var_node = var;
DCHECK(RefIsPositive(var));
// We will create a bunch of nodes linked to a variable node. Only one node
// per (var, type) will be required so we cache them to not create more node
// than necessary.
absl::flat_hash_map<std::vector<int64>, int, VectorHash> secondary_var_nodes;
auto get_secondary_var_node = [&new_node, &add_edge, &secondary_var_nodes,
&tmp_key](int var_node,
absl::Span<const int64> type) {
tmp_key.assign(type.begin(), type.end());
tmp_key.push_back(var_node);
const auto insert = secondary_var_nodes.insert({tmp_key, 0});
// For a coefficient of one, which are the most common, we can optimize the
// size of the graph by omitting the coefficient node altogether and using
// directly the var_node in this case.
if (coeff == 1) return var_node;
const auto insert =
coefficient_nodes.insert({std::make_pair(var, coeff), 0});
if (!insert.second) return insert.first->second;
tmp_key.pop_back();
const int secondary_node = new_node(tmp_key);
add_edge(var_node, secondary_node);
tmp_color = {VAR_COEFFICIENT_NODE, coeff};
const int secondary_node = new_node(tmp_color);
graph->AddArc(var_node, secondary_node);
insert.first->second = secondary_node;
return secondary_node;
};
auto add_literal_edge = [&add_edge, &get_secondary_var_node](
int ref, int constraint_node) {
const int variable_node = PositiveRef(ref);
if (RefIsPositive(ref)) {
// For all coefficients equal to one, which are the most common, we
// can optimize the size of the graph by omitting the coefficient
// node altogether.
add_edge(variable_node, constraint_node);
} else {
const int coefficient_node = get_secondary_var_node(
variable_node, {CONSTRAINT_COEFFICIENT_NODE, -1});
add_edge(coefficient_node, constraint_node);
}
// For a literal we use the same as a coefficient 1 or -1. We can do that
// because literal and (var, coefficient) never appear together in the same
// constraint.
auto get_literal_node = [&get_coefficient_node](int ref) {
return get_coefficient_node(PositiveRef(ref), RefIsPositive(ref) ? 1 : -1);
};
// Because the implication can be numerous, we encode them without constraints
// node by using an arc from the lhs to the rhs. Note that we also always add
// the other direction. We use a set to remove duplicates both for efficiency
// and to not artificially break symmetries by using multi-arcs.
absl::flat_hash_set<std::pair<int, int>> implications;
auto add_implication = [&get_literal_node, &graph, &implications](int ref_a,
int ref_b) {
const auto insert = implications.insert({ref_a, ref_b});
if (!insert.second) return;
graph->AddArc(get_literal_node(ref_a), get_literal_node(ref_b));
// Always add the other side.
implications.insert({NegatedRef(ref_b), NegatedRef(ref_a)});
graph->AddArc(get_literal_node(NegatedRef(ref_b)),
get_literal_node(NegatedRef(ref_a)));
};
// Add constraints to the graph.
for (const ConstraintProto& constraint : problem.constraints()) {
const int constraint_node = initial_equivalence_classes->size();
std::vector<int64> key = {CONSTRAINT_NODE, constraint.constraint_case()};
std::vector<int64> color = {CONSTRAINT_NODE, constraint.constraint_case()};
switch (constraint.constraint_case()) {
case ConstraintProto::CONSTRAINT_NOT_SET:
break;
case ConstraintProto::kLinear: {
Append(constraint.linear().domain(), &key);
CHECK_EQ(constraint_node, new_node(key));
// TODO(user): We can use the same trick as for the implications to
// encode relations of the form coeff * var_a <= coeff * var_b without
// creating a constraint node by directly adding an arc between the two
// var coefficient nodes.
Append(constraint.linear().domain(), &color);
CHECK_EQ(constraint_node, new_node(color));
for (int i = 0; i < constraint.linear().vars_size(); ++i) {
if (constraint.linear().coeffs(i) == 0) continue;
const int ref = constraint.linear().vars(i);
const int variable_node = PositiveRef(ref);
const int64 coeff = RefIsPositive(ref)
? constraint.linear().coeffs(i)
: -constraint.linear().coeffs(i);
if (coeff == 1) {
// For all coefficients equal to one, which are the most common, we
// can optimize the size of the graph by omitting the coefficient
// node altogether.
add_edge(variable_node, constraint_node);
} else {
const int coefficient_node = get_secondary_var_node(
variable_node, {CONSTRAINT_COEFFICIENT_NODE, coeff});
add_edge(coefficient_node, constraint_node);
}
graph->AddArc(get_coefficient_node(variable_node, coeff),
constraint_node);
}
break;
}
case ConstraintProto::kBoolOr: {
CHECK_EQ(constraint_node, new_node(key));
CHECK_EQ(constraint_node, new_node(color));
for (const int ref : constraint.bool_or().literals()) {
add_literal_edge(ref, constraint_node);
graph->AddArc(get_literal_node(ref), constraint_node);
}
break;
}
case ConstraintProto::kAtMostOne: {
CHECK_EQ(constraint_node, new_node(key));
if (constraint.at_most_one().literals().size() == 2) {
// Treat it as an implication to avoid creating a node.
add_implication(constraint.at_most_one().literals(0),
NegatedRef(constraint.at_most_one().literals(1)));
break;
}
CHECK_EQ(constraint_node, new_node(color));
for (const int ref : constraint.at_most_one().literals()) {
add_literal_edge(ref, constraint_node);
graph->AddArc(get_literal_node(ref), constraint_node);
}
break;
}
case ConstraintProto::kExactlyOne: {
CHECK_EQ(constraint_node, new_node(key));
CHECK_EQ(constraint_node, new_node(color));
for (const int ref : constraint.exactly_one().literals()) {
add_literal_edge(ref, constraint_node);
graph->AddArc(get_literal_node(ref), constraint_node);
}
break;
}
case ConstraintProto::kBoolXor: {
CHECK_EQ(constraint_node, new_node(key));
CHECK_EQ(constraint_node, new_node(color));
for (const int ref : constraint.bool_xor().literals()) {
add_literal_edge(ref, constraint_node);
graph->AddArc(get_literal_node(ref), constraint_node);
}
break;
}
// TODO(user): We could directly connect variable node together to
// deal more efficiently with this constraint. Make sure not to create
// multi-arc since I am not sure the symmetry code works with these
// though.
case ConstraintProto::kBoolAnd: {
if (constraint.enforcement_literal_size() == 0) {
// All literals are true in this case.
CHECK_EQ(constraint_node, new_node(key));
for (const int ref : constraint.bool_and().literals()) {
add_literal_edge(ref, constraint_node);
}
break;
// The other cases should be presolved before this is called.
// TODO(user): not 100% true, this happen on rmatr200-p5, Fix.
if (constraint.enforcement_literal_size() != 1) {
VLOG(1) << "BoolAnd with multiple enforcement literal are not "
"supported in symmetry code:"
<< constraint.ShortDebugString();
return nullptr;
}
// To make the BoolAnd constraint more generic in the graph, we expand
// it into a set of BoolOr constraints where
// not(enforcement_literal) OR literal = true
// for all constraint's literals. This is equivalent to
// enforcement_literal => literal
// for all literals.
std::vector<int64> key = {CONSTRAINT_NODE, ConstraintProto::kBoolOr};
const int non_enforcement_literal =
NegatedRef(constraint.enforcement_literal(0));
for (const int literal : constraint.bool_and().literals()) {
const int constraint_node = new_node(key);
add_literal_edge(non_enforcement_literal, constraint_node);
add_literal_edge(literal, constraint_node);
CHECK_EQ(constraint.enforcement_literal_size(), 1);
const int ref_a = constraint.enforcement_literal(0);
for (const int ref_b : constraint.bool_and().literals()) {
add_implication(ref_a, ref_b);
}
break;
}
@@ -274,11 +278,13 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
}
}
// For enforcement, we use a similar trick than for the implications.
// Because all our constraint arcs are in the direction var_node to
// constraint_node, we just use the reverse direction for the enforcement
// part. This way we can reuse the same get_literal_node() function.
if (constraint.constraint_case() != ConstraintProto::kBoolAnd) {
for (const int ref : constraint.enforcement_literal()) {
const int enforcement_node = get_secondary_var_node(
PositiveRef(ref), {ENFORCEMENT_LITERAL, RefIsPositive(ref)});
add_edge(constraint_node, enforcement_node);
graph->AddArc(constraint_node, get_literal_node(ref));
}
}
}
@@ -294,7 +300,7 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
//
// TODO(user): It will probably be more efficient to not even create these
// nodes, but we will need a mapping to know the variable <-> node index.
int next_id = id_generator.NextFreeId();
int next_id = color_id_generator.NextFreeId();
for (int i = 0; i < num_variables; ++i) {
if ((*graph)[i].empty()) {
(*initial_equivalence_classes)[i] = next_id++;
@@ -319,7 +325,7 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
void FindCpModelSymmetries(
const SatParameters& params, const CpModelProto& problem,
std::vector<std::unique_ptr<SparsePermutation>>* generators,
double time_limit_seconds) {
double deterministic_limit) {
const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
CHECK(generators != nullptr);
generators->clear();
@@ -333,15 +339,17 @@ void FindCpModelSymmetries(
if (log_info) {
LOG(INFO) << "Graph for symmetry has " << graph->num_nodes()
<< " nodes and " << graph->num_arcs() / 2 << " edges.";
<< " nodes and " << graph->num_arcs() << " arcs.";
}
if (graph->num_nodes() == 0) return;
GraphSymmetryFinder symmetry_finder(*graph, /*is_undirected=*/true);
GraphSymmetryFinder symmetry_finder(*graph, /*is_undirected=*/false);
std::vector<int> factorized_automorphism_group_size;
std::unique_ptr<TimeLimit> time_limit =
TimeLimit::FromDeterministicTime(deterministic_limit);
const absl::Status status = symmetry_finder.FindSymmetries(
time_limit_seconds, &equivalence_classes, generators,
&factorized_automorphism_group_size);
&equivalence_classes, generators, &factorized_automorphism_group_size,
time_limit.get());
// TODO(user): Change the API to not return an error when the time limit is
// reached.
@@ -384,12 +392,17 @@ void FindCpModelSymmetries(
}
generators->resize(num_generators);
average_support_size /= num_generators;
if (log_info && num_generators > 0) {
LOG(INFO) << "# of generators: " << num_generators;
LOG(INFO) << "Average support size: " << average_support_size;
if (num_duplicate_constraints > 0) {
LOG(INFO) << "The model contains " << num_duplicate_constraints
<< " duplicate constraints !";
if (log_info) {
LOG(INFO) << "Symmetry computation done. time: "
<< time_limit->GetElapsedTime()
<< " dtime: " << time_limit->GetElapsedDeterministicTime();
if (num_generators > 0) {
LOG(INFO) << "# of generators: " << num_generators;
LOG(INFO) << "Average support size: " << average_support_size;
if (num_duplicate_constraints > 0) {
LOG(INFO) << "The model contains " << num_duplicate_constraints
<< " duplicate constraints !";
}
}
}
}
@@ -463,7 +476,7 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
// Collect the at most ones.
//
// Note(user): This relies on the fact that the pointer remain stable when
// Note(user): This relies on the fact that the pointers remain stable when
// we adds new constraints. It should be the case, but it is a bit unsafe.
// On the other hand it is annoying to deal with both cases below.
std::vector<const google::protobuf::RepeatedField<int32>*> at_most_ones;
@@ -496,7 +509,8 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
}
}
while (!orbitope.empty() && !orbitope[0].empty()) {
while (!orbitope.empty() && orbitope[0].size() > 1) {
const int num_cols = orbitope[0].size();
const std::vector<int> orbits = GetOrbitopeOrbits(num_vars, orbitope);
// Because in the orbitope case, we have a full symmetry group of the
@@ -528,25 +542,39 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
// Note(user): On the miplib, only 1/ happens currently. Not sure with LNS
// though.
std::vector<bool> all_equivalent_rows(orbitope.size(), false);
std::vector<bool> at_most_one_one(orbitope.size(), false);
std::vector<bool> at_most_one_zero(orbitope.size(), false);
// The result described above can be generalized if an at most one intersect
// many of the orbitope rows, each in at leat two positions. We will track
// the set of best rows on which we have an at most one (or at most one
// zero) on all their entries.
bool at_most_one_in_best_rows; // The alternative is at most one zero.
int64 best_score = 0;
std::vector<int> best_rows;
std::vector<int> rows_in_at_most_one;
for (const google::protobuf::RepeatedField<int32>* literals :
at_most_ones) {
tmp_to_clear.clear();
int num_in_intersections = 0;
for (const int literal : *literals) {
if (context->IsFixed(literal)) continue;
const int var = PositiveRef(literal);
const int rep = orbits[var];
if (rep == -1) continue;
++num_in_intersections;
if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
tmp_sizes[rep]++;
if (RefIsPositive(literal)) tmp_num_positive[rep]++;
}
int num_positive_direction = 0;
int num_negative_direction = 0;
// An at most one touching two positions in an orbitope row can possibly
// be extended, depending if it has singleton intersection swith other
// rows and where.
bool possible_extension = false;
rows_in_at_most_one.clear();
for (const int row : tmp_to_clear) {
const int size = tmp_sizes[row];
const int num_positive = tmp_num_positive[row];
@@ -555,37 +583,55 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
tmp_num_positive[row] = 0;
if (num_positive > 1 && num_negative == 0) {
at_most_one_one[row] = true;
if (size < num_cols) possible_extension = true;
rows_in_at_most_one.push_back(row);
++num_positive_direction;
} else if (num_positive == 0 && num_negative > 1) {
at_most_one_zero[row] = true;
if (size < num_cols) possible_extension = true;
rows_in_at_most_one.push_back(row);
++num_negative_direction;
} else if (num_positive > 0 && num_negative > 0) {
all_equivalent_rows[row] = true;
}
}
// We might be able to presolve more in these cases.
if (at_most_one_zero[row] || at_most_one_one[row] ||
all_equivalent_rows[row]) {
if (tmp_to_clear.size() > 1) {
context->UpdateRuleStats(
"TODO symmetry: at most one across orbits.");
} else if (size < orbitope[0].size()) {
context->UpdateRuleStats(
"TODO symmetry: at most one can be extended");
}
}
if (possible_extension) {
context->UpdateRuleStats(
"TODO symmetry: possible at most one extension.");
}
if (num_positive_direction > 0 && num_negative_direction > 0) {
return context->NotifyThatModelIsUnsat("Symmetry and at most ones");
}
const bool direction = num_positive_direction > 0;
// Because of symmetry, the choice of the column shouldn't matter (they
// will all appear in the same number of constraints of the same types),
// however we prefer to fix the variables that seems to touch more
// constraints.
//
// TODO(user): maybe we should simplify the constraint using the variable
// we fix before choosing the next row to break symmetry on. If there are
// multiple row involved, we could also take the intersection instead of
// probably counting the same constraints more than once.
int64 score = 0;
for (const int row : rows_in_at_most_one) {
score +=
context->VarToConstraints(PositiveRef(orbitope[row][0])).size();
}
if (score > best_score) {
at_most_one_in_best_rows = direction;
best_score = score;
best_rows = rows_in_at_most_one;
}
}
// Heuristically choose a "best" row/col to "distinguish" and break the
// symmetry on.
int best_row = 0;
int best_col = 0;
int best_score = 0;
bool fix_others_to_zero = true;
// Mark all the equivalence.
// Note that this operation do not change the symmetry group.
//
// TODO(user): We could remove these rows from the orbitope. Note that
// currently this never happen on the miplib (maybe in LNS though).
for (int i = 0; i < all_equivalent_rows.size(); ++i) {
const int num_cols = orbitope[i].size();
// Note that this operation do not change the symmetry group.
if (all_equivalent_rows[i]) {
for (int j = 1; j < num_cols; ++j) {
context->StoreBooleanEqualityRelation(orbitope[i][0], orbitope[i][j]);
@@ -593,94 +639,68 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
if (context->ModelIsUnsat()) return false;
}
}
// Because of symmetry, the choice of the column shouldn't matter (they
// will all appear in the same number of constraints of the same types),
// however we prefer to fix a variable that seems to touch more
// constraints.
//
// TODO(user): maybe we should simplify the constraint using the variable
// we fix before choosing the next row to break symmetry on.
const int row_score =
context->VarToConstraints(PositiveRef(orbitope[i][0])).size();
// TODO(user): If one variable make the line already fixed, we should just
// ignore this row. Not too important as actually this shouldn't happen
// because we never compute symmetries involving fixed variables. But in
// the future, fixing some literal might have some side effects and fix
// others.
if (at_most_one_one[i] && row_score > best_score) {
best_col = 0;
for (int j = 0; j < num_cols; ++j) {
if (context->LiteralIsTrue(orbitope[i][j])) {
best_col = j;
break;
}
}
best_row = i;
best_score = row_score;
fix_others_to_zero = true;
}
if (at_most_one_zero[i] && row_score > best_score) {
best_col = 0;
for (int j = 0; j < num_cols; ++j) {
if (context->LiteralIsFalse(orbitope[i][j])) {
best_col = j;
break;
}
}
best_row = i;
best_score = row_score;
fix_others_to_zero = false;
}
}
if (best_score == 0) break;
for (int j = 1; j < orbitope[best_row].size(); ++j) {
if (fix_others_to_zero) {
context->UpdateRuleStats("symmetry: fixed to false");
if (!context->SetLiteralToFalse(orbitope[best_row][j])) return false;
} else {
context->UpdateRuleStats("symmetry: fixed to true");
if (!context->SetLiteralToTrue(orbitope[best_row][j])) return false;
}
}
// We add the symmetry breaking inequalities: best_var >= all other var
// in orbit. That is not(best_var) => not(other) for Booleans. We only add
// them if we didn't fix any variable just above.
// Break the symmetry on our set of best rows by picking one columns
// and setting all the other entries to zero or one. Note that the at most
// one applies to all entries in all rows.
//
// TODO(user): Add the inequality for non-Boolean too? Also note that this
// code only run if the code above is disabled. It is here for testing
// alternatives. In particular, if there is no at most one, we cannot fix
// n-1 variables, but we can still add inequalities.
const int best_var = orbitope[best_row][best_col];
const bool maximize_best_var = fix_others_to_zero;
if (context->CanBeUsedAsLiteral(best_var) && !context->IsFixed(best_var)) {
ConstraintProto* ct = context->working_model->add_constraints();
ct->add_enforcement_literal(maximize_best_var ? NegatedRef(best_var)
: best_var);
for (const int other : orbitope[best_row]) {
if (other == best_var) continue;
if (context->IsFixed(other)) continue;
ct->mutable_bool_and()->add_literals(
maximize_best_var ? NegatedRef(other) : other);
context->UpdateRuleStats("symmetry: added implication");
}
context->UpdateNewConstraintsVariableUsage();
// TODO(user): We don't have any at most one relation on this orbitope,
// but we could still add symmetry breaking inequality by picking any matrix
// entry and making it the largest/lowest value on its row. This also work
// for non-Booleans.
if (best_score == 0) {
context->UpdateRuleStats(
"TODO symmetry: add symmetry breaking inequalities?");
break;
}
// Remove the column of best_var.
CHECK_NE(best_col, -1);
// If our symmetry group is valid, they cannot be any variable already
// fixed to one (or zero if !at_most_one_in_best_rows). Otherwise all would
// be fixed to one and the problem would be unsat.
for (const int i : best_rows) {
for (int j = 0; j < num_cols; ++j) {
const int var = orbitope[i][j];
if ((at_most_one_in_best_rows && context->LiteralIsTrue(var)) ||
(!at_most_one_in_best_rows && context->LiteralIsFalse(var))) {
return context->NotifyThatModelIsUnsat("Symmetry and at most one");
}
}
}
// We have an at most one on a set of rows, we will pick a column, and set
// all other entries on these rows to zero.
//
// TODO(user): All choices should be equivalent, but double check?
const int best_col = 0;
for (const int i : best_rows) {
for (int j = 0; j < num_cols; ++j) {
if (j == best_col) continue;
const int var = orbitope[i][j];
if (at_most_one_in_best_rows) {
context->UpdateRuleStats("symmetry: fixed to false");
if (!context->SetLiteralToFalse(var)) return false;
} else {
context->UpdateRuleStats("symmetry: fixed to true");
if (!context->SetLiteralToTrue(var)) return false;
}
}
}
// Remove all best rows.
for (const int i : best_rows) orbitope[i].clear();
int new_size = 0;
for (int i = 0; i < orbitope.size(); ++i) {
if (!orbitope[i].empty()) orbitope[new_size++] = orbitope[i];
}
CHECK_LT(new_size, orbitope.size());
orbitope.resize(new_size);
// Remove best_col.
for (int i = 0; i < orbitope.size(); ++i) {
std::swap(orbitope[i][best_col], orbitope[i].back());
orbitope[i].pop_back();
}
// We also remove the line of best_var since heuristicially, it is better to
// not add symmetries involving any of the variable on this line.
std::swap(orbitope[best_row], orbitope.back());
orbitope.pop_back();
}
return true;

4
ortools/sat/cp_model_symmetries.h
View File

@@ -41,12 +41,10 @@ namespace sat {
// TODO(user): As long as we only exploit symmetry involving only Boolean
// variables we can make this code more efficient by not detecting symmetries
// involing integer variable.
//
// TODO(user): Make the limit deterministic.
void FindCpModelSymmetries(
const SatParameters& params, const CpModelProto& problem,
std::vector<std::unique_ptr<SparsePermutation>>* generators,
double time_limit_seconds = std::numeric_limits<double>::infinity());
double deterministic_limit = std::numeric_limits<double>::infinity());
// Basic implementation of some symmetry breaking during presolve.
//

7
ortools/sat/sat_parameters.proto
View File

@@ -21,7 +21,7 @@ option java_multiple_files = true;
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 184
// NEXT TAG: 185
message SatParameters {
// In some context, like in a portfolio of search, it makes sense to name a
// given parameters set for logging purpose.
@@ -686,8 +686,8 @@ message SatParameters {
// Target number of constraints to remove during cleanup.
optional int32 cut_cleanup_target = 157 [default = 1000];
// Add that many lazy contraints (or cuts) at once in the LP. Note that at the
// beginning of the solve, we do add more than this.
// Add that many lazy constraints (or cuts) at once in the LP. Note that at
// the beginning of the solve, we do add more than this.
optional int32 new_constraints_batch_size = 122 [default = 50];
// The search branching will be used to decide how to branch on unfixed nodes.
@@ -872,6 +872,7 @@ message SatParameters {
// LNS parameters.
optional bool use_lns_only = 101 [default = false];
optional bool lns_focus_on_decision_variables = 105 [default = false];
optional bool lns_expand_intervals_in_constraint_graph = 184 [default = true];
// Turns on relaxation induced neighborhood generator.
optional bool use_rins_lns = 129 [default = true];