From dbaef62eb5e2cbc75c5852c37a8a92d46d878c8b Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 21 Apr 2023 12:48:03 +0200 Subject: [PATCH] [CP-SAT] improvements to feasibility_jump, better presolve on lin_max; tuning --- ortools/sat/BUILD.bazel | 1 + ortools/sat/constraint_violation.cc | 974 ++++++++++++++----- ortools/sat/constraint_violation.h | 271 +++++- ortools/sat/cp_model_postsolve.cc | 10 +- ortools/sat/cp_model_presolve.cc | 62 +- ortools/sat/cp_model_solver.cc | 15 +- ortools/sat/feasibility_jump.cc | 592 ++++++++--- ortools/sat/feasibility_jump.h | 53 +- ortools/sat/linear_constraint_manager.cc | 43 +- ortools/sat/linear_constraint_manager.h | 10 +- ortools/sat/linear_programming_constraint.cc | 29 +- ortools/sat/linear_programming_constraint.h | 4 +- ortools/sat/parameters_validation.cc | 11 +- ortools/sat/presolve_context.cc | 12 + ortools/sat/presolve_context.h | 5 +- ortools/sat/sat_parameters.proto | 30 +- ortools/sat/synchronization.cc | 21 +- ortools/sat/synchronization.h | 12 +- ortools/sat/work_assignment.cc | 17 +- 19 files changed, 1675 insertions(+), 497 deletions(-) diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index a7adce7dc9..67c78b885a 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -176,6 +176,7 @@ cc_library( ":util", "//ortools/base", "//ortools/base:hash", + "//ortools/graph:strongly_connected_components", "//ortools/port:proto_utils", "//ortools/util:saturated_arithmetic", "//ortools/util:sorted_interval_list", diff --git a/ortools/sat/constraint_violation.cc b/ortools/sat/constraint_violation.cc index e8c861089d..68642c6b6d 100644 --- a/ortools/sat/constraint_violation.cc +++ b/ortools/sat/constraint_violation.cc @@ -17,13 +17,17 @@ #include #include #include +#include #include #include +#include "absl/algorithm/container.h" +#include "absl/container/flat_hash_set.h" #include "absl/log/check.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" +#include "ortools/graph/strongly_connected_components.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/util.h" @@ -85,6 +89,7 @@ int LinearIncrementalEvaluator::NewConstraint(Domain domain) { activities_.push_back(0); num_false_enforcement_.push_back(0); distances_.push_back(0); + is_violated_.push_back(false); return num_constraints_++; } @@ -156,12 +161,12 @@ void LinearIncrementalEvaluator::ComputeInitialActivities( // Resets the activity as the offset and the number of false enforcement to 0. activities_ = offsets_; + in_last_affected_variables_.resize(columns_.size(), false); num_false_enforcement_.assign(num_constraints_, 0); - row_update_will_not_impact_deltas_.assign(num_constraints_, false); // Update these numbers for all columns. for (int var = 0; var < columns_.size(); ++var) { - const VarData& data = columns_[var]; + const SpanData& data = columns_[var]; const int64_t value = solution[var]; int i = data.start; @@ -186,6 +191,7 @@ void LinearIncrementalEvaluator::ComputeInitialActivities( // Cache violations (not counting enforcement). for (int c = 0; c < num_constraints_; ++c) { distances_[c] = domains_[c].Distance(activities_[c]); + is_violated_[c] = Violation(c) > 0; } } @@ -195,61 +201,351 @@ void LinearIncrementalEvaluator::Update(int var, int64_t delta) { DCHECK_NE(delta, 0); if (var >= columns_.size()) return; - const VarData& data = columns_[var]; + const SpanData& data = columns_[var]; int i = data.start; for (int k = 0; k < data.num_pos_literal; ++k, ++i) { const int c = ct_buffer_[i]; if (delta == 1) { num_false_enforcement_[c]--; DCHECK_GE(num_false_enforcement_[c], 0); - - // A transition 3 -> 2 or 2 -> 1 with a distance of zero can be ignored. - row_update_will_not_impact_deltas_[c] = - num_false_enforcement_[c] >= 2 || - (num_false_enforcement_[c] == 1 && distances_[c] == 0); } else { num_false_enforcement_[c]++; - - // A transition 2 -> 3 or 1 -> 2 with a distance of zero can be ignored. - row_update_will_not_impact_deltas_[c] = - num_false_enforcement_[c] >= 3 || - (num_false_enforcement_[c] == 2 && distances_[c] == 0); } + is_violated_[c] = Violation(c) > 0; } for (int k = 0; k < data.num_neg_literal; ++k, ++i) { const int c = ct_buffer_[i]; if (delta == -1) { num_false_enforcement_[c]--; DCHECK_GE(num_false_enforcement_[c], 0); - - // A transition 3 -> 2 or 2 -> 1 with a distance of zero can be ignored. - row_update_will_not_impact_deltas_[c] = - num_false_enforcement_[c] >= 2 || - (num_false_enforcement_[c] == 1 && distances_[c] == 0); } else { num_false_enforcement_[c]++; - - // A transition 2 -> 3 or 1 -> 2 with a distance of zero can be ignored. - row_update_will_not_impact_deltas_[c] = - num_false_enforcement_[c] >= 3 || - (num_false_enforcement_[c] == 2 && distances_[c] == 0); } + is_violated_[c] = Violation(c) > 0; } int j = data.linear_start; for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { const int c = ct_buffer_[i]; const int64_t coeff = coeff_buffer_[j]; activities_[c] += coeff * delta; - - const int64_t old_distance = distances_[c]; distances_[c] = domains_[c].Distance(activities_[c]); + is_violated_[c] = Violation(c) > 0; + } +} - // If this constraint enforcement was >= 2, then the change - // in activity cannot impact any deltas. Or if enforcement - // is >= 1 and the distance before and after is the same. - row_update_will_not_impact_deltas_[c] = - num_false_enforcement_[c] >= 2 || - (num_false_enforcement_[c] >= 1 && distances_[c] == old_distance); +void LinearIncrementalEvaluator::ClearAffectedVariables() { + in_last_affected_variables_.resize(columns_.size(), false); + for (const int var : last_affected_variables_) { + in_last_affected_variables_[var] = false; + } + last_affected_variables_.clear(); +} + +void LinearIncrementalEvaluator::UpdateScoreOnWeightUpdate( + int c, double weight_delta, absl::Span solution, + absl::Span jump_values, absl::Span jump_scores) { + if (c >= rows_.size()) return; + + DCHECK_EQ(num_false_enforcement_[c], 0); + const SpanData& data = rows_[c]; + + // Update enforcement part, all change are 0 -> 1 transition and decrease + // by this. + const double enforcement_change = + weight_delta * static_cast(distances_[c]); + if (enforcement_change != 0.0) { + int i = data.start; + const int end = data.num_pos_literal + data.num_neg_literal; + dtime_ += end; + for (int k = 0; k < end; ++k, ++i) { + const int var = row_var_buffer_[i]; + jump_scores[var] -= enforcement_change; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } + + // Update linear part. + int i = data.start + data.num_pos_literal + data.num_neg_literal; + int j = data.linear_start; + dtime_ += 2 * data.num_linear_entries; + const int64_t old_distance = distances_[c]; + for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { + const int var = row_var_buffer_[i]; + const int64_t coeff = row_coeff_buffer_[j]; + const int64_t delta = jump_values[var] - solution[var]; + const int64_t new_distance = + domains_[c].Distance(activities_[c] + coeff * delta); + jump_scores[var] += + weight_delta * static_cast(new_distance - old_distance); + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } +} + +void LinearIncrementalEvaluator::UpdateScoreOnNewlyEnforced( + int c, double weight, absl::Span solution, + absl::Span jump_values, absl::Span jump_scores) { + const SpanData& data = rows_[c]; + + // Everyone else had a zero cost transition that now become enforced -> + // unenforced. + const double weight_time_violation = + weight * static_cast(distances_[c]); + if (weight_time_violation > 0.0) { + int i = data.start; + const int end = data.num_pos_literal + data.num_neg_literal; + dtime_ += end; + for (int k = 0; k < end; ++k, ++i) { + const int var = row_var_buffer_[i]; + jump_scores[var] -= weight_time_violation; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } + + // Update linear part! It was zero and is now a diff. + { + int i = data.start + data.num_pos_literal + data.num_neg_literal; + int j = data.linear_start; + dtime_ += 2 * data.num_linear_entries; + const int64_t old_distance = distances_[c]; + for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { + const int var = row_var_buffer_[i]; + const int64_t coeff = row_coeff_buffer_[j]; + const int64_t delta = jump_values[var] - solution[var]; + const int64_t new_distance = + domains_[c].Distance(activities_[c] + coeff * delta); + jump_scores[var] += + weight * static_cast(new_distance - old_distance); + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } +} + +void LinearIncrementalEvaluator::UpdateScoreOnNewlyUnenforced( + int c, double weight, absl::Span solution, + absl::Span jump_values, absl::Span jump_scores) { + const SpanData& data = rows_[c]; + + // Everyone else had a enforced -> unenforced transition that now become zero. + const double weight_time_violation = + weight * static_cast(distances_[c]); + if (weight_time_violation > 0.0) { + int i = data.start; + const int end = data.num_pos_literal + data.num_neg_literal; + dtime_ += end; + for (int k = 0; k < end; ++k, ++i) { + const int var = row_var_buffer_[i]; + jump_scores[var] += weight_time_violation; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } + + // Update linear part! It had a diff and is now zero. + { + int i = data.start + data.num_pos_literal + data.num_neg_literal; + int j = data.linear_start; + dtime_ += 2 * data.num_linear_entries; + const int64_t old_distance = distances_[c]; + for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { + const int var = row_var_buffer_[i]; + const int64_t coeff = row_coeff_buffer_[j]; + const int64_t delta = jump_values[var] - solution[var]; + const int64_t new_distance = + domains_[c].Distance(activities_[c] + coeff * delta); + jump_scores[var] -= + weight * static_cast(new_distance - old_distance); + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } +} + +// We just need to modifie the old/new transition that decrease the number of +// enforcement literal at false. +void LinearIncrementalEvaluator::UpdateScoreOfEnforcementIncrease( + int c, double score_change, absl::Span jump_values, + absl::Span jump_scores) { + if (score_change == 0.0) return; + + const SpanData& data = rows_[c]; + int i = data.start; + dtime_ += data.num_pos_literal; + for (int k = 0; k < data.num_pos_literal; ++k, ++i) { + const int var = row_var_buffer_[i]; + if (jump_values[var] == 1) { + jump_scores[var] += score_change; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } + dtime_ += data.num_neg_literal; + for (int k = 0; k < data.num_neg_literal; ++k, ++i) { + const int var = row_var_buffer_[i]; + if (jump_values[var] == 0) { + jump_scores[var] += score_change; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } +} + +void LinearIncrementalEvaluator::UpdateScoreOnActivityChange( + int c, double weight, int64_t activity_delta, + absl::Span solution, absl::Span jump_values, + absl::Span jump_scores) { + if (activity_delta == 0) return; + const SpanData& data = rows_[c]; + + // Enforcement is always enforced -> unforced. + // So it was -weight_time_distance and is now -weight_time_new_distance. + const double delta = + -weight * static_cast( + domains_[c].Distance(activities_[c] + activity_delta) - + distances_[c]); + if (delta != 0.0) { + int i = data.start; + const int end = data.num_pos_literal + data.num_neg_literal; + dtime_ += end; + for (int k = 0; k < end; ++k, ++i) { + const int var = row_var_buffer_[i]; + jump_scores[var] += delta; + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } + + // Update linear part. + { + int i = data.start + data.num_pos_literal + data.num_neg_literal; + int j = data.linear_start; + dtime_ += 2 * data.num_linear_entries; + const int64_t old_a_minus_new_a = + distances_[c] - domains_[c].Distance(activities_[c] + activity_delta); + for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { + const int var = row_var_buffer_[i]; + const int64_t coeff = row_coeff_buffer_[j]; + const int64_t delta = jump_values[var] - solution[var]; + const int64_t old_b = + domains_[c].Distance(activities_[c] + coeff * delta); + const int64_t new_b = + domains_[c].Distance(activities_[c] + activity_delta + coeff * delta); + + // The old score was: + // weight * static_cast(old_b - old_a); + // the new score is + // weight * static_cast(new_b - new_a); so the diff is: + jump_scores[var] += + weight * static_cast(old_a_minus_new_a + new_b - old_b); + if (!in_last_affected_variables_[var]) { + last_affected_variables_.push_back(var); + } + } + } +} + +void LinearIncrementalEvaluator::UpdateVariableAndScores( + int var, int64_t delta, absl::Span solution, + absl::Span weights, absl::Span jump_values, + absl::Span jump_scores) { + DCHECK(!creation_phase_); + DCHECK_NE(delta, 0); + if (var >= columns_.size()) return; + + const SpanData& data = columns_[var]; + int i = data.start; + for (int k = 0; k < data.num_pos_literal; ++k, ++i) { + const int c = ct_buffer_[i]; + if (delta == 1) { + num_false_enforcement_[c]--; + DCHECK_GE(num_false_enforcement_[c], 0); + if (num_false_enforcement_[c] == 0) { + UpdateScoreOnNewlyEnforced(c, weights[c], solution, jump_values, + jump_scores); + } else if (num_false_enforcement_[c] == 1) { + const double enforcement_change = + weights[c] * static_cast(distances_[c]); + UpdateScoreOfEnforcementIncrease(c, enforcement_change, jump_values, + jump_scores); + } + } else { + num_false_enforcement_[c]++; + if (num_false_enforcement_[c] == 1) { + UpdateScoreOnNewlyUnenforced(c, weights[c], solution, jump_values, + jump_scores); + } else if (num_false_enforcement_[c] == 2) { + const double enforcement_change = + weights[c] * static_cast(distances_[c]); + UpdateScoreOfEnforcementIncrease(c, -enforcement_change, jump_values, + jump_scores); + } + } + is_violated_[c] = Violation(c) > 0; + } + for (int k = 0; k < data.num_neg_literal; ++k, ++i) { + const int c = ct_buffer_[i]; + if (delta == -1) { + num_false_enforcement_[c]--; + DCHECK_GE(num_false_enforcement_[c], 0); + if (num_false_enforcement_[c] == 0) { + UpdateScoreOnNewlyEnforced(c, weights[c], solution, jump_values, + jump_scores); + } else if (num_false_enforcement_[c] == 1) { + const double enforcement_change = + weights[c] * static_cast(distances_[c]); + UpdateScoreOfEnforcementIncrease(c, enforcement_change, jump_values, + jump_scores); + } + } else { + num_false_enforcement_[c]++; + if (num_false_enforcement_[c] == 1) { + UpdateScoreOnNewlyUnenforced(c, weights[c], solution, jump_values, + jump_scores); + } else if (num_false_enforcement_[c] == 2) { + const double enforcement_change = + weights[c] * static_cast(distances_[c]); + UpdateScoreOfEnforcementIncrease(c, -enforcement_change, jump_values, + jump_scores); + } + } + is_violated_[c] = Violation(c) > 0; + } + int j = data.linear_start; + for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { + const int c = ct_buffer_[i]; + const int64_t coeff = coeff_buffer_[j]; + + if (num_false_enforcement_[c] == 1) { + // Only the 1 -> 0 are impacted. + // This is the same as the 1->2 transition, but the old 1->0 needs to + // be changed from - weight * distance to - weight * new_distance. + const int64_t new_distance = + domains_[c].Distance(activities_[c] + coeff * delta); + if (new_distance != distances_[c]) { + UpdateScoreOfEnforcementIncrease( + c, -weights[c] * static_cast(distances_[c] - new_distance), + jump_values, jump_scores); + } + } else if (num_false_enforcement_[c] == 0) { + UpdateScoreOnActivityChange(c, weights[c], coeff * delta, solution, + jump_values, jump_scores); + } + + activities_[c] += coeff * delta; + distances_[c] = domains_[c].Distance(activities_[c]); + is_violated_[c] = Violation(c) > 0; } } @@ -261,6 +557,10 @@ int64_t LinearIncrementalEvaluator::Violation(int c) const { return num_false_enforcement_[c] > 0 ? 0 : distances_[c]; } +bool LinearIncrementalEvaluator::IsViolated(int c) const { + return is_violated_[c]; +} + void LinearIncrementalEvaluator::ReduceBounds(int c, int64_t lb, int64_t ub) { domains_[c] = domains_[c].IntersectionWith(Domain(lb, ub)); distances_[c] = domains_[c].Distance(activities_[c]); @@ -286,7 +586,7 @@ double LinearIncrementalEvaluator::WeightedViolationDelta( absl::Span weights, int var, int64_t delta) const { DCHECK_NE(delta, 0); if (var >= columns_.size()) return 0.0; - const VarData& data = columns_[var]; + const SpanData& data = columns_[var]; int i = data.start; double result = 0.0; @@ -319,7 +619,7 @@ double LinearIncrementalEvaluator::WeightedViolationDelta( } int j = data.linear_start; - dtime_ += data.num_linear_entries; + dtime_ += 2 * data.num_linear_entries; for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { const int c = ct_buffer_[i]; if (num_false_enforcement_[c] > 0) continue; @@ -346,7 +646,7 @@ std::vector LinearIncrementalEvaluator::SlopeBreakpoints( std::vector result = var_domain.FlattenedIntervals(); if (var_domain.Size() <= 2 || var >= columns_.size()) return result; - const VarData& data = columns_[var]; + const SpanData& data = columns_[var]; int i = data.start + data.num_pos_literal + data.num_neg_literal; int j = data.linear_start; for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) { @@ -390,20 +690,50 @@ void LinearIncrementalEvaluator::PrecomputeCompactView() { creation_phase_ = false; if (num_constraints_ == 0) return; + // Compute the total size. + // Note that at this point the constraint indices are not "encoded" yet. + int total_size = 0; + int total_linear_size = 0; + tmp_row_sizes_.assign(num_constraints_, 0); + tmp_row_num_positive_literals_.assign(num_constraints_, 0); + tmp_row_num_negative_literals_.assign(num_constraints_, 0); + tmp_row_num_linear_entries_.assign(num_constraints_, 0); + for (const auto& column : literal_entries_) { + total_size += column.size(); + for (const auto [c, is_positive] : column) { + tmp_row_sizes_[c]++; + if (is_positive) { + tmp_row_num_positive_literals_[c]++; + } else { + tmp_row_num_negative_literals_[c]++; + } + } + } + for (const auto& column : var_entries_) { + total_size += column.size(); + total_linear_size += column.size(); + for (const auto& entry : column) { + tmp_row_sizes_[entry.ct_index]++; + tmp_row_num_linear_entries_[entry.ct_index]++; + } + } + // Compactify for faster WeightedViolationDelta(). + ct_buffer_.reserve(total_size); + coeff_buffer_.reserve(total_linear_size); columns_.resize(std::max(literal_entries_.size(), var_entries_.size())); for (int var = 0; var < columns_.size(); ++var) { columns_[var].start = static_cast(ct_buffer_.size()); columns_[var].linear_start = static_cast(coeff_buffer_.size()); if (var < literal_entries_.size()) { - for (const auto [c, positive] : literal_entries_[var]) { - if (positive) { + for (const auto [c, is_positive] : literal_entries_[var]) { + if (is_positive) { columns_[var].num_pos_literal++; ct_buffer_.push_back(c); } } - for (const auto [c, positive] : literal_entries_[var]) { - if (!positive) { + for (const auto [c, is_positive] : literal_entries_[var]) { + if (!is_positive) { columns_[var].num_neg_literal++; ct_buffer_.push_back(c); } @@ -418,101 +748,74 @@ void LinearIncrementalEvaluator::PrecomputeCompactView() { } } - // Compute the total size. - // Note that at this point the constraint indices are not "encoded" yet. - int total_size = 0; - tmp_row_sizes_.assign(num_constraints_, 0); - for (auto& column : literal_entries_) { - for (auto& entry : column) { - total_size++; - tmp_row_sizes_[entry.ct_index]++; - } - } - for (auto& column : var_entries_) { - for (auto& entry : column) { - total_size++; - tmp_row_sizes_[entry.ct_index]++; - } - } - row_buffer_.resize(total_size); - - // Corner case, some constraints but no entry !? we need size 1 for MakeSpan() - // below not to crash. - // - // TODO(user): we should probably filter empty constraint before here (this - // only happens in test). - if (total_size == 0) row_buffer_.resize(1); - - // transform tmp_row_sizes_ to starts in the buffer. - // Initialize the spans. - int offset = 0; - rows_.resize(num_constraints_); - for (int c = 0; c < num_constraints_; ++c) { - const int start = offset; - offset += tmp_row_sizes_[c]; - rows_[c] = absl::MakeSpan(&row_buffer_[start], tmp_row_sizes_[c]); - tmp_row_sizes_[c] = start; - } - DCHECK_EQ(offset, total_size); - - // Copy data. - for (int var = 0; var < literal_entries_.size(); ++var) { - for (const auto [c, _] : literal_entries_[var]) { - row_buffer_[tmp_row_sizes_[c]++] = var; - } - } - for (int var = 0; var < var_entries_.size(); ++var) { - for (const auto [c, _] : var_entries_[var]) { - row_buffer_[tmp_row_sizes_[c]++] = var; - } - } - // We do not need var_entries_ or literal_entries_ anymore. // // TODO(user): We could delete them before. But at the time of this // optimization, I didn't want to change the behavior of the algorithm at all. gtl::STLClearObject(&var_entries_); gtl::STLClearObject(&literal_entries_); + + // Initialize the SpanData. + // Transform tmp_row_sizes_ to starts in the row_var_buffer_. + // Transform tmp_row_num_linear_entries_ to starts in the row_coeff_buffer_. + int offset = 0; + int linear_offset = 0; + rows_.resize(num_constraints_); + for (int c = 0; c < num_constraints_; ++c) { + rows_[c].num_pos_literal = tmp_row_num_positive_literals_[c]; + rows_[c].num_neg_literal = tmp_row_num_negative_literals_[c]; + rows_[c].num_linear_entries = tmp_row_num_linear_entries_[c]; + + rows_[c].start = offset; + offset += tmp_row_sizes_[c]; + tmp_row_sizes_[c] = rows_[c].start; + + rows_[c].linear_start = linear_offset; + linear_offset += tmp_row_num_linear_entries_[c]; + tmp_row_num_linear_entries_[c] = rows_[c].linear_start; + } + DCHECK_EQ(offset, total_size); + DCHECK_EQ(linear_offset, total_linear_size); + + // Copy data. + row_var_buffer_.resize(total_size); + row_coeff_buffer_.resize(total_linear_size); + for (int var = 0; var < columns_.size(); ++var) { + const SpanData& data = columns_[var]; + int i = data.start; + for (int k = 0; k < data.num_pos_literal; ++i, ++k) { + const int c = ct_buffer_[i]; + row_var_buffer_[tmp_row_sizes_[c]++] = var; + } + } + for (int var = 0; var < columns_.size(); ++var) { + const SpanData& data = columns_[var]; + int i = data.start + data.num_pos_literal; + for (int k = 0; k < data.num_neg_literal; ++i, ++k) { + const int c = ct_buffer_[i]; + row_var_buffer_[tmp_row_sizes_[c]++] = var; + } + } + for (int var = 0; var < columns_.size(); ++var) { + const SpanData& data = columns_[var]; + int i = data.start + data.num_pos_literal + data.num_neg_literal; + int j = data.linear_start; + for (int k = 0; k < data.num_linear_entries; ++i, ++j, ++k) { + const int c = ct_buffer_[i]; + row_var_buffer_[tmp_row_sizes_[c]++] = var; + row_coeff_buffer_[tmp_row_num_linear_entries_[c]++] = coeff_buffer_[j]; + } + } + + cached_deltas_.assign(columns_.size(), 0); + cached_scores_.assign(columns_.size(), 0); } -// TODO(user): Depending on the status of an enforced constraint, we might not -// need to recompute data. We do a bit of that but we could do more. -std::vector LinearIncrementalEvaluator::AffectedVariableOnUpdate(int var) { - std::vector result; - if (var >= columns_.size()) return result; - tmp_in_result_.resize(columns_.size(), false); +bool LinearIncrementalEvaluator::ViolationChangeIsConvex(int var) const { for (const int c : VarToConstraints(var)) { - if (row_update_will_not_impact_deltas_[c]) continue; - for (const int affected : rows_[c]) { - if (!tmp_in_result_[affected]) { - tmp_in_result_[affected] = true; - result.push_back(affected); - } - } + if (domains_[c].intervals().size() > 2) return false; } - for (const int var : result) tmp_in_result_[var] = false; - return result; -} - -// Note that this is meant to be called on weight update and as such only on -// enforced constraint. No need to check row_update_will_not_impact_deltas_. -std::vector LinearIncrementalEvaluator::ConstraintsToAffectedVariables( - absl::Span ct_indices) { - std::vector result; - tmp_in_result_.resize(columns_.size(), false); - for (const int c : ct_indices) { - // This is needed since constraint indices can refer to more general - // constraint than linear ones. - if (c >= rows_.size()) continue; - for (const int affected : rows_[c]) { - if (!tmp_in_result_[affected]) { - tmp_in_result_[affected] = true; - result.push_back(affected); - } - } - } - for (const int var : result) tmp_in_result_[var] = false; - return result; + return true; } // ----- CompiledConstraint ----- @@ -520,38 +823,50 @@ std::vector LinearIncrementalEvaluator::ConstraintsToAffectedVariables( CompiledConstraint::CompiledConstraint(const ConstraintProto& ct_proto) : ct_proto_(ct_proto) {} -void CompiledConstraint::CacheViolation(absl::Span solution) { - violation_ = Violation(solution); +void CompiledConstraint::InitializeViolation( + absl::Span solution) { + violation_ = ComputeViolation(solution); } -int64_t CompiledConstraint::violation() const { - DCHECK_GE(violation_, 0); - return violation_; +void CompiledConstraint::PerformMove( + int var, int64_t old_value, + absl::Span solution_with_new_value) { + violation_ += ViolationDelta(var, old_value, solution_with_new_value); } -const ConstraintProto& CompiledConstraint::ct_proto() const { - return ct_proto_; +int64_t CompiledConstraint::ViolationDelta(int /*var*/, int64_t /*old_value*/, + absl::Span solution) { + return ComputeViolation(solution) - violation_; +} + +// ----- CompiledBoolXorConstraint ----- + +CompiledBoolXorConstraint::CompiledBoolXorConstraint( + const ConstraintProto& ct_proto) + : CompiledConstraint(ct_proto) {} + +int64_t CompiledBoolXorConstraint::ComputeViolation( + absl::Span solution) { + int64_t sum_of_literals = 0; + for (const int lit : ct_proto().bool_xor().literals()) { + sum_of_literals += LiteralValue(lit, solution); + } + return 1 - (sum_of_literals % 2); +} + +int64_t CompiledBoolXorConstraint::ViolationDelta( + int /*var*/, int64_t /*old_value*/, + absl::Span /*solution_with_new_value*/) { + return violation() == 0 ? 1 : -1; } // ----- CompiledLinMaxConstraint ----- -// The violation of a lin_max constraint is: -// - the sum(max(0, expr_value - target_value) forall expr). This part will be -// maintained by the linear part. -// - target_value - max(expressions) if positive. -class CompiledLinMaxConstraint : public CompiledConstraint { - public: - explicit CompiledLinMaxConstraint(const ConstraintProto& ct_proto); - ~CompiledLinMaxConstraint() override = default; - - int64_t Violation(absl::Span solution) override; -}; - CompiledLinMaxConstraint::CompiledLinMaxConstraint( const ConstraintProto& ct_proto) : CompiledConstraint(ct_proto) {} -int64_t CompiledLinMaxConstraint::Violation( +int64_t CompiledLinMaxConstraint::ComputeViolation( absl::Span solution) { const int64_t target_value = ExprValue(ct_proto().lin_max().target(), solution); @@ -565,21 +880,11 @@ int64_t CompiledLinMaxConstraint::Violation( // ----- CompiledIntProdConstraint ----- -// The violation of an int_prod constraint is -// abs(value(target) - prod(value(expr)). -class CompiledIntProdConstraint : public CompiledConstraint { - public: - explicit CompiledIntProdConstraint(const ConstraintProto& ct_proto); - ~CompiledIntProdConstraint() override = default; - - int64_t Violation(absl::Span solution) override; -}; - CompiledIntProdConstraint::CompiledIntProdConstraint( const ConstraintProto& ct_proto) : CompiledConstraint(ct_proto) {} -int64_t CompiledIntProdConstraint::Violation( +int64_t CompiledIntProdConstraint::ComputeViolation( absl::Span solution) { const int64_t target_value = ExprValue(ct_proto().int_prod().target(), solution); @@ -592,21 +897,11 @@ int64_t CompiledIntProdConstraint::Violation( // ----- CompiledIntDivConstraint ----- -// The violation of an int_div constraint is -// abs(value(target) - value(expr0) / value(expr1)). -class CompiledIntDivConstraint : public CompiledConstraint { - public: - explicit CompiledIntDivConstraint(const ConstraintProto& ct_proto); - ~CompiledIntDivConstraint() override = default; - - int64_t Violation(absl::Span solution) override; -}; - CompiledIntDivConstraint::CompiledIntDivConstraint( const ConstraintProto& ct_proto) : CompiledConstraint(ct_proto) {} -int64_t CompiledIntDivConstraint::Violation( +int64_t CompiledIntDivConstraint::ComputeViolation( absl::Span solution) { const int64_t target_value = ExprValue(ct_proto().int_div().target(), solution); @@ -618,24 +913,11 @@ int64_t CompiledIntDivConstraint::Violation( // ----- CompiledAllDiffConstraint ----- -// The violation of a all_diff is the number of unordered pairs of expressions -// with the same value. -class CompiledAllDiffConstraint : public CompiledConstraint { - public: - explicit CompiledAllDiffConstraint(const ConstraintProto& ct_proto); - ~CompiledAllDiffConstraint() override = default; - - int64_t Violation(absl::Span solution) override; - - private: - std::vector values_; -}; - CompiledAllDiffConstraint::CompiledAllDiffConstraint( const ConstraintProto& ct_proto) : CompiledConstraint(ct_proto) {} -int64_t CompiledAllDiffConstraint::Violation( +int64_t CompiledAllDiffConstraint::ComputeViolation( absl::Span solution) { values_.clear(); for (const LinearExpressionProto& expr : ct_proto().all_diff().exprs()) { @@ -688,62 +970,45 @@ int64_t ComputeOverloadArea( const int64_t max_end = std::max(start + size, end); if (start >= max_end) continue; - events.push_back({start, demand}); - events.push_back({max_end, -demand}); + events.emplace_back(start, demand); + events.emplace_back(max_end, -demand); } if (events.empty()) return 0; std::sort(events.begin(), events.end(), [](const std::pair& e1, const std::pair& e2) { - return e1.first < e2.first || - (e1.first == e2.first && e1.second < e2.second); + return e1.first < e2.first; }); + int64_t overload = 0; int64_t current_load = 0; - int64_t previous_load = 0; int64_t previous_time = events.front().first; - for (int i = 0; i < events.size();) { - if (previous_load > max_capacity) { - overload = - CapAdd(overload, CapProd(CapSub(previous_load, max_capacity), - CapSub(events[i].first, previous_time))); + // At this point, current_load is the load at previous_time. + const int64_t time = events[i].first; + if (current_load > max_capacity) { + overload = CapAdd( + overload, CapProd(current_load - max_capacity, time - previous_time)); } - const int start_index = i; - while (i < events.size() && events[i].first == events[start_index].first) { + while (i < events.size() && events[i].first == time) { current_load += events[i].second; i++; } - CHECK_GE(current_load, 0); - previous_time = events[start_index].first; - previous_load = current_load; + DCHECK_GE(current_load, 0); + previous_time = time; } - CHECK_EQ(current_load, 0); + DCHECK_EQ(current_load, 0); return overload; } } // namespace -// The violation of a no_overlap is the sum of overloads over time. -class CompiledNoOverlapConstraint : public CompiledConstraint { - public: - explicit CompiledNoOverlapConstraint(const ConstraintProto& ct_proto, - const CpModelProto& cp_model); - ~CompiledNoOverlapConstraint() override = default; - - int64_t Violation(absl::Span solution) override; - - private: - const CpModelProto& cp_model_; - std::vector> events_; -}; - CompiledNoOverlapConstraint::CompiledNoOverlapConstraint( const ConstraintProto& ct_proto, const CpModelProto& cp_model) : CompiledConstraint(ct_proto), cp_model_(cp_model) {} -int64_t CompiledNoOverlapConstraint::Violation( +int64_t CompiledNoOverlapConstraint::ComputeViolation( absl::Span solution) { return ComputeOverloadArea(ct_proto().no_overlap().intervals(), {}, cp_model_, solution, 1, events_); @@ -751,25 +1016,11 @@ int64_t CompiledNoOverlapConstraint::Violation( // ----- CompiledCumulativeConstraint ----- -// The violation of a cumulative is the sum of overloads over time. -class CompiledCumulativeConstraint : public CompiledConstraint { - public: - explicit CompiledCumulativeConstraint(const ConstraintProto& ct_proto, - const CpModelProto& cp_model); - ~CompiledCumulativeConstraint() override = default; - - int64_t Violation(absl::Span solution) override; - - private: - const CpModelProto& cp_model_; - std::vector> events_; -}; - CompiledCumulativeConstraint::CompiledCumulativeConstraint( const ConstraintProto& ct_proto, const CpModelProto& cp_model) : CompiledConstraint(ct_proto), cp_model_(cp_model) {} -int64_t CompiledCumulativeConstraint::Violation( +int64_t CompiledCumulativeConstraint::ComputeViolation( absl::Span solution) { return ComputeOverloadArea( ct_proto().cumulative().intervals(), ct_proto().cumulative().demands(), @@ -777,10 +1028,172 @@ int64_t CompiledCumulativeConstraint::Violation( ExprValue(ct_proto().cumulative().capacity(), solution), events_); } +// ----- CompiledCircuitConstraint ----- + +// The violation of a circuit has three parts: +// 1. Flow imbalance, maintained by the linear part. +// 2. The number of non-skipped SCCs in the graph minus 1. +// 3. The number of non-skipped SCCs that cannot be reached from any other +// component minus 1. +// +// #3 is not necessary for correctness, but makes the function much smoother. +// +// The only difference between single and multi circuit is flow balance at the +// depot, so we use the same compiled constraint for both. +class CompiledCircuitConstraint : public CompiledConstraint { + public: + explicit CompiledCircuitConstraint(const ConstraintProto& ct_proto); + ~CompiledCircuitConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; + + private: + struct SccOutput { + void emplace_back(const int* start, const int* end); + void reset(int num_nodes); + + int num_components = 0; + std::vector skipped; + std::vector root; + }; + void UpdateGraph(absl::Span solution); + absl::Span literals_; + absl::Span tails_; + absl::Span heads_; + // Stores the currently active arcs per tail node. + std::vector> graph_; + SccOutput sccs_; + std::vector has_in_arc_; + StronglyConnectedComponentsFinder>, + SccOutput> + scc_finder_; +}; + +void CompiledCircuitConstraint::SccOutput::emplace_back(int const* start, + int const* end) { + const int root_node = *start; + const int size = end - start; + if (size == 1) { + skipped[root_node] = true; + } else { + ++num_components; + } + for (; start != end; ++start) { + root[*start] = root_node; + } +} +void CompiledCircuitConstraint::SccOutput::reset(int num_nodes) { + num_components = 0; + root.clear(); + root.resize(num_nodes); + skipped.clear(); + skipped.resize(num_nodes); +} + +CompiledCircuitConstraint::CompiledCircuitConstraint( + const ConstraintProto& ct_proto) + : CompiledConstraint(ct_proto) { + const bool routes = ct_proto.has_routes(); + tails_ = routes ? ct_proto.routes().tails() : ct_proto.circuit().tails(); + heads_ = absl::MakeConstSpan(routes ? ct_proto.routes().heads() + : ct_proto.circuit().heads()); + literals_ = absl::MakeConstSpan(routes ? ct_proto.routes().literals() + : ct_proto.circuit().literals()); + graph_.resize(*absl::c_max_element(tails_) + 1); +} + +void CompiledCircuitConstraint::UpdateGraph( + absl::Span solution) { + for (std::vector& edges : graph_) { + edges.clear(); + } + for (int i = 0; i < tails_.size(); ++i) { + if (!LiteralValue(literals_[i], solution)) continue; + graph_[tails_[i]].push_back(heads_[i]); + } +} +int64_t CompiledCircuitConstraint::ComputeViolation( + absl::Span solution) { + const int num_nodes = graph_.size(); + sccs_.reset(num_nodes); + UpdateGraph(solution); + scc_finder_.FindStronglyConnectedComponents(num_nodes, graph_, &sccs_); + // Skipping all nodes causes off-by-one errors below, so it's simpler to + // handle explicitly. + if (sccs_.num_components == 0) return 0; + // Count the number of SCCs that have inbound cross-component arcs + // as a smoother measure of progress towards strong connectivity. + int num_half_connected_components = 0; + has_in_arc_.clear(); + has_in_arc_.resize(num_nodes, false); + for (int tail = 0; tail < graph_.size(); ++tail) { + if (sccs_.skipped[tail]) continue; + for (const int head : graph_[tail]) { + const int head_root = sccs_.root[head]; + if (sccs_.root[tail] == head_root) continue; + if (has_in_arc_[head_root]) continue; + if (sccs_.skipped[head_root]) continue; + has_in_arc_[head_root] = true; + ++num_half_connected_components; + } + } + const int64_t violation = sccs_.num_components - 1 + sccs_.num_components - + num_half_connected_components - 1 + + (ct_proto().has_routes() ? sccs_.skipped[0] : 0); + VLOG(2) << "#SCCs=" << sccs_.num_components << " #nodes=" << num_nodes + << " #half_connected_components=" << num_half_connected_components + << " violation=" << violation; + return violation; +} + +void AddCircuitFlowConstraints(LinearIncrementalEvaluator& linear_evaluator, + const ConstraintProto& ct_proto) { + const bool routes = ct_proto.has_routes(); + auto heads = routes ? ct_proto.routes().heads() : ct_proto.circuit().heads(); + auto tails = routes ? ct_proto.routes().tails() : ct_proto.circuit().tails(); + auto literals = + routes ? ct_proto.routes().literals() : ct_proto.circuit().literals(); + + std::vector> inflow_lits; + std::vector> outflow_lits; + for (int i = 0; i < heads.size(); ++i) { + if (heads[i] >= inflow_lits.size()) { + inflow_lits.resize(heads[i] + 1); + } + inflow_lits[heads[i]].push_back(literals[i]); + if (tails[i] >= outflow_lits.size()) { + outflow_lits.resize(tails[i] + 1); + } + outflow_lits[tails[i]].push_back(literals[i]); + } + if (routes) { + const int depot_net_flow = linear_evaluator.NewConstraint({0, 0}); + for (const int lit : inflow_lits[0]) { + linear_evaluator.AddTerm(depot_net_flow, lit, 1); + } + for (const int lit : outflow_lits[0]) { + linear_evaluator.AddTerm(depot_net_flow, lit, -1); + } + } + for (int i = routes ? 1 : 0; i < inflow_lits.size(); ++i) { + const int inflow_ct = linear_evaluator.NewConstraint({1, 1}); + for (const int lit : inflow_lits[i]) { + linear_evaluator.AddLiteral(inflow_ct, lit); + } + } + for (int i = routes ? 1 : 0; i < outflow_lits.size(); ++i) { + const int outflow_ct = linear_evaluator.NewConstraint({1, 1}); + for (const int lit : outflow_lits[i]) { + linear_evaluator.AddLiteral(outflow_ct, lit); + } + } +} + // ----- LsEvaluator ----- LsEvaluator::LsEvaluator(const CpModelProto& model) : model_(model) { var_to_constraint_graph_.resize(model_.variables_size()); + jump_value_optimal_.resize(model_.variables_size(), true); CompileConstraintsAndObjective(); BuildVarConstraintGraph(); } @@ -808,6 +1221,25 @@ void LsEvaluator::BuildVarConstraintGraph() { for (std::vector& deps : var_to_constraint_graph_) { gtl::STLSortAndRemoveDuplicates(&deps); } + + // Scan the model to decide if a variable is linked to a convex evaluation. + jump_value_optimal_.resize(model_.variables_size()); + for (int i = 0; i < model_.variables_size(); ++i) { + if (!var_to_constraint_graph_[i].empty()) { + jump_value_optimal_[i] = false; + continue; + } + + const IntegerVariableProto& var_proto = model_.variables(i); + if (var_proto.domain_size() == 2 && var_proto.domain(0) == 0 && + var_proto.domain(1) == 1) { + // Boolean variables violation change is always convex. + jump_value_optimal_[i] = true; + continue; + } + + jump_value_optimal_[i] = linear_evaluator_.ViolationChangeIsConvex(i); + } } void LsEvaluator::CompileConstraintsAndObjective() { @@ -868,6 +1300,10 @@ void LsEvaluator::CompileConstraintsAndObjective() { } break; } + case ConstraintProto::ConstraintCase::kBoolXor: { + constraints_.emplace_back(new CompiledBoolXorConstraint(ct)); + break; + } case ConstraintProto::ConstraintCase::kAllDiff: { constraints_.emplace_back(new CompiledAllDiffConstraint(ct)); break; @@ -920,6 +1356,11 @@ void LsEvaluator::CompileConstraintsAndObjective() { constraints_.emplace_back(new CompiledCumulativeConstraint(ct, model_)); break; } + case ConstraintProto::ConstraintCase::kCircuit: + case ConstraintProto::ConstraintCase::kRoutes: + constraints_.emplace_back(new CompiledCircuitConstraint(ct)); + AddCircuitFlowConstraints(linear_evaluator_, ct); + break; default: VLOG(1) << "Not implemented: " << ct.constraint_case(); model_is_supported_ = false; @@ -933,7 +1374,7 @@ void LsEvaluator::CompileConstraintsAndObjective() { void LsEvaluator::ReduceObjectiveBounds(int64_t lb, int64_t ub) { if (!model_.has_objective()) return; - linear_evaluator_.ReduceBounds(/*ct_index=*/0, lb, ub); + linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub); } void LsEvaluator::ComputeInitialViolations(absl::Span solution) { @@ -944,32 +1385,43 @@ void LsEvaluator::ComputeInitialViolations(absl::Span solution) { // Generic constraints. for (const auto& ct : constraints_) { - ct->CacheViolation(current_solution_); + ct->InitializeViolation(current_solution_); } } -void LsEvaluator::UpdateNonLinearViolations() { +void LsEvaluator::UpdateAllNonLinearViolations() { // Generic constraints. for (const auto& ct : constraints_) { - ct->CacheViolation(current_solution_); + ct->InitializeViolation(current_solution_); } } -void LsEvaluator::UpdateVariableValueAndRecomputeViolations( - int var, int64_t new_value, bool focus_on_linear) { - DCHECK(RefIsPositive(var)); +void LsEvaluator::UpdateNonLinearViolations(int var, int64_t new_value) { const int64_t old_value = current_solution_[var]; if (old_value == new_value) return; - // Linear part. - linear_evaluator_.Update(var, new_value - old_value); current_solution_[var] = new_value; - if (focus_on_linear) return; - - // Non linear part. for (const int ct_index : var_to_constraint_graph_[var]) { - constraints_[ct_index]->CacheViolation(current_solution_); + constraints_[ct_index]->PerformMove(var, old_value, current_solution_); } + current_solution_[var] = old_value; +} + +void LsEvaluator::UpdateLinearScores(int var, int64_t value, + absl::Span weights, + absl::Span jump_values, + absl::Span jump_scores) { + DCHECK(RefIsPositive(var)); + const int64_t old_value = current_solution_[var]; + if (old_value == value) return; + + linear_evaluator_.UpdateVariableAndScores(var, value - old_value, + current_solution_, weights, + jump_values, jump_scores); +} + +void LsEvaluator::UpdateVariableValue(int var, int64_t new_value) { + current_solution_[var] = new_value; } int64_t LsEvaluator::SumOfViolations() { @@ -978,18 +1430,20 @@ int64_t LsEvaluator::SumOfViolations() { // Process the linear part. for (int i = 0; i < linear_evaluator_.num_constraints(); ++i) { evaluation += linear_evaluator_.Violation(i); + DCHECK_GE(linear_evaluator_.Violation(i), 0); } // Process the generic constraint part. for (const auto& ct : constraints_) { evaluation += ct->violation(); + DCHECK_GE(ct->violation(), 0); } return evaluation; } int64_t LsEvaluator::ObjectiveActivity() const { DCHECK(model_.has_objective()); - return linear_evaluator_.Activity(/*ct_index=*/0); + return linear_evaluator_.Activity(/*c=*/0); } int LsEvaluator::NumLinearConstraints() const { @@ -1028,6 +1482,15 @@ int64_t LsEvaluator::Violation(int c) const { } } +bool LsEvaluator::IsViolated(int c) const { + if (c < linear_evaluator_.num_constraints()) { + return linear_evaluator_.IsViolated(c); + } else { + return constraints_[c - linear_evaluator_.num_constraints()]->violation() > + 0; + } +} + double LsEvaluator::WeightedViolation(absl::Span weights) const { DCHECK_EQ(weights.size(), NumEvaluatorConstraints()); double violations = linear_evaluator_.WeightedViolation(weights); @@ -1040,30 +1503,32 @@ double LsEvaluator::WeightedViolation(absl::Span weights) const { return violations; } -double LsEvaluator::WeightedViolationDelta(absl::Span weights, - int var, int64_t delta) const { - DCHECK_LT(var, current_solution_.size()); - ++num_deltas_computed_; - double violation_delta = - linear_evaluator_.WeightedViolationDelta(weights, var, delta); - +double LsEvaluator::WeightedNonLinearViolationDelta( + absl::Span weights, int var, int64_t delta) const { + const int64_t old_value = current_solution_[var]; + double violation_delta = 0; // We change the mutable solution here, are restore it after the evaluation. current_solution_[var] += delta; const int num_linear_constraints = linear_evaluator_.num_constraints(); for (const int ct_index : var_to_constraint_graph_[var]) { DCHECK_LT(ct_index, constraints_.size()); - const int64_t ct_eval_before = constraints_[ct_index]->violation(); - const int64_t ct_eval_after = - constraints_[ct_index]->Violation(current_solution_); - violation_delta += static_cast(ct_eval_after - ct_eval_before) * - weights[ct_index + num_linear_constraints]; + const int64_t delta = constraints_[ct_index]->ViolationDelta( + var, old_value, current_solution_); + violation_delta += + static_cast(delta) * weights[ct_index + num_linear_constraints]; } // Restore. current_solution_[var] -= delta; - return violation_delta; } +double LsEvaluator::WeightedViolationDelta(absl::Span weights, + int var, int64_t delta) const { + DCHECK_LT(var, current_solution_.size()); + return linear_evaluator_.WeightedViolationDelta(weights, var, delta) + + WeightedNonLinearViolationDelta(weights, var, delta); +} + // TODO(user): Speed-up: // - Use a row oriented representation of the model (could reuse the Apply // methods on proto). @@ -1087,5 +1552,10 @@ std::vector LsEvaluator::VariablesInViolatedConstraints() const { return variables; } +bool LsEvaluator::VariableOnlyInLinearConstraintWithConvexViolationChange( + int var) const { + return jump_value_optimal_[var]; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/constraint_violation.h b/ortools/sat/constraint_violation.h index 99dc7f24c0..b4a43d0503 100644 --- a/ortools/sat/constraint_violation.h +++ b/ortools/sat/constraint_violation.h @@ -17,6 +17,7 @@ #include #include #include +#include #include #include "absl/types/span.h" @@ -51,17 +52,39 @@ class LinearIncrementalEvaluator { // and before the class starts to be used. This is DCHECKed. void PrecomputeCompactView(); - // Compute activities and query violations. + // Compute activities and update them. void ComputeInitialActivities(absl::Span solution); void Update(int var, int64_t delta); + + // Function specific to the feasibility jump heuristic. + // Note that the score of the changed variable will not be updated correctly! + void UpdateVariableAndScores(int var, int64_t delta, + absl::Span solution, + absl::Span weights, + absl::Span jump_values, + absl::Span jump_scores); + void UpdateScoreOnWeightUpdate(int c, double weight_delta, + absl::Span solution, + absl::Span jump_values, + absl::Span jump_scores); + + // Variables whose score was updated since the last clear. + // Note that we also include variable whose score might be the same but for + // which different jump values might have a changed score. + void ClearAffectedVariables(); + const std::vector& VariablesAffectedByLastUpdate() const { + return last_affected_variables_; + } + + // Query violation. int64_t Activity(int c) const; int64_t Violation(int c) const; + bool IsViolated(int c) const; + bool AppearsInViolatedConstraints(int var) const; // Used to DCHECK the state of the evaluator. bool VarIsConsistent(int var) const; - bool AppearsInViolatedConstraints(int var) const; - // Intersect constraint bounds with [lb..ub]. void ReduceBounds(int c, int64_t lb, int64_t ub); @@ -77,8 +100,8 @@ class LinearIncrementalEvaluator { double WeightedViolationDelta(absl::Span weights, int var, int64_t delta) const; - // The violation for each constraints is a piecewise linear function. This - // computes and aggregate all the breakpoints for the given variable and its + // The violation for each constraint is a piecewise linear function. This + // computes and aggregates all the breakpoints for the given variable and its // domain. // // Note that if the domain contains less than two values, we return an empty @@ -86,13 +109,8 @@ class LinearIncrementalEvaluator { std::vector SlopeBreakpoints(int var, int64_t current_value, const Domain& var_domain) const; - // AffectedVariableOnUpdate() return the set of variables that have at least - // one constraint in common with the given variable. One need to call - // PrecomputeCompactView() after all the constraint have been added for - // this to work. - std::vector AffectedVariableOnUpdate(int var); - std::vector ConstraintsToAffectedVariables( - absl::Span ct_indices); + // Checks if the jump value of a variable is always optimal. + bool ViolationChangeIsConvex(int var) const; double DeterministicTime() const { return 5e-9 * static_cast(dtime_); @@ -111,6 +129,50 @@ class LinearIncrementalEvaluator { bool positive; // bool_var or its negation. }; + struct SpanData { + int start = 0; + int num_pos_literal = 0; + int num_neg_literal = 0; + int linear_start = 0; + int num_linear_entries = 0; + }; + + absl::Span VarToConstraints(int var) const { + if (var >= columns_.size()) return {}; + const SpanData& data = columns_[var]; + const int size = + data.num_pos_literal + data.num_neg_literal + data.num_linear_entries; + if (size == 0) return {}; + return absl::MakeSpan(&ct_buffer_[data.start], size); + } + + absl::Span ConstraintToVars(int c) const { + const SpanData& data = rows_[c]; + const int size = + data.num_pos_literal + data.num_neg_literal + data.num_linear_entries; + if (size == 0) return {}; + return absl::MakeSpan(&row_var_buffer_[data.start], size); + } + + void ComputeAndCacheDistance(int ct_index); + + // Incremental row-based update. + void UpdateScoreOnNewlyEnforced(int c, double weight, + absl::Span solution, + absl::Span jump_values, + absl::Span jump_scores); + void UpdateScoreOnNewlyUnenforced(int c, double weight, + absl::Span solution, + absl::Span jump_values, + absl::Span jump_scores); + void UpdateScoreOfEnforcementIncrease(int c, double score_change, + absl::Span jump_values, + absl::Span jump_scores); + void UpdateScoreOnActivityChange(int c, double weight, int64_t activity_delta, + absl::Span solution, + absl::Span jump_values, + absl::Span jump_scores); + // Constraint indexed data (static). int num_constraints_ = 0; std::vector domains_; @@ -119,41 +181,40 @@ class LinearIncrementalEvaluator { // Variable indexed data. // Note that this is just used at construction and is replaced by a compact // view when PrecomputeCompactView() is called. + bool creation_phase_ = true; std::vector> var_entries_; std::vector> literal_entries_; - // Memory efficient column based data. - struct VarData { - int start = 0; - int num_pos_literal = 0; - int num_neg_literal = 0; - int linear_start = 0; - int num_linear_entries = 0; - }; - absl::Span VarToConstraints(int var) const { - const VarData& data = columns_[var]; - const int size = - data.num_pos_literal + data.num_neg_literal + data.num_linear_entries; - return absl::MakeSpan(&ct_buffer_[data.start], size); - } - std::vector columns_; + // Memory efficient column based data (static). + std::vector columns_; std::vector ct_buffer_; std::vector coeff_buffer_; - // Transposed view data (only variables). - bool creation_phase_ = true; - std::vector row_buffer_; - std::vector> rows_; + // Memory efficient row based data (static). + std::vector rows_; + std::vector row_var_buffer_; + std::vector row_coeff_buffer_; // Temporary data. std::vector tmp_row_sizes_; - std::vector tmp_in_result_; + std::vector tmp_row_num_positive_literals_; + std::vector tmp_row_num_negative_literals_; + std::vector tmp_row_num_linear_entries_; // Constraint indexed data (dynamic). + std::vector is_violated_; std::vector activities_; std::vector distances_; std::vector num_false_enforcement_; - std::vector row_update_will_not_impact_deltas_; + + // Code to update the scores on a variable change. + std::vector old_distances_; + std::vector old_num_false_enforcement_; + std::vector cached_deltas_; + std::vector cached_scores_; + + std::vector in_last_affected_variables_; + std::vector last_affected_variables_; mutable size_t dtime_ = 0; }; @@ -167,22 +228,31 @@ class CompiledConstraint { explicit CompiledConstraint(const ConstraintProto& ct_proto); virtual ~CompiledConstraint() = default; + // Recomputes the violation of the constraint from scratch. + void InitializeViolation(absl::Span solution); + + // Update the violation with the new value. + void PerformMove(int var, int64_t old_value, + absl::Span solution_with_new_value); + // Computes the violation of a constraint. // // A violation is a positive integer value. A zero value means the constraint - // is not violated.. - virtual int64_t Violation(absl::Span solution) = 0; + // is not violated. + virtual int64_t ComputeViolation(absl::Span solution) = 0; - void CacheViolation(absl::Span solution); + // Returns the delta if var changes from old_value to solution[var]. + virtual int64_t ViolationDelta( + int var, int64_t old_value, + absl::Span solution_with_new_value); - // Cache the violation of a constraint. - int64_t violation() const; - - const ConstraintProto& ct_proto() const; + // Getters. + const ConstraintProto& ct_proto() const { return ct_proto_; } + int64_t violation() const { return violation_; } private: const ConstraintProto& ct_proto_; - int64_t violation_ = -1; + int64_t violation_; }; // Evaluation container for the local search. @@ -205,15 +275,24 @@ class LsEvaluator { // Intersects the domain of the objective with [lb..ub]. void ReduceObjectiveBounds(int64_t lb, int64_t ub); - // Sets the current solution, and compute violations for each constraints. + // Sets the current solution, and computes violations for each constraint. void ComputeInitialViolations(absl::Span solution); // Recompute the violations of non linear constraints. - void UpdateNonLinearViolations(); + void UpdateAllNonLinearViolations(); - // Update the value of one variable and recompute violations. - void UpdateVariableValueAndRecomputeViolations(int var, int64_t value, - bool focus_on_linear = false); + // Sets the value of the variable in the current solution. + // It must be called after UpdateLinearScores(). + void UpdateVariableValue(int var, int64_t new_value); + + // Recomputes the violations of all impacted non linear constraints. + void UpdateNonLinearViolations(int var, int64_t new_value); + + // Function specific to the linear only feasibility jump. + void UpdateLinearScores(int var, int64_t value, + absl::Span weights, + absl::Span jump_values, + absl::Span jump_scores); // Simple summation metric for the constraint and objective violations. int64_t SumOfViolations(); @@ -231,9 +310,13 @@ class LsEvaluator { // Returns the weighted sum of violation. The weights must have the same // size as NumEvaluatorConstraints(). int64_t Violation(int c) const; + bool IsViolated(int c) const; double WeightedViolation(absl::Span weights) const; double WeightedViolationDelta(absl::Span weights, int var, int64_t delta) const; + // Ignores the violations of the linear constraints. + double WeightedNonLinearViolationDelta(absl::Span weights, + int var, int64_t delta) const; LinearIncrementalEvaluator* MutableLinearEvaluator() { return &linear_evaluator_; @@ -242,11 +325,10 @@ class LsEvaluator { // Returns the set of variables appearing in a violated constraint. std::vector VariablesInViolatedConstraints() const; - // Counters: number of times UpdateVariableValueAndRecomputeViolations() has - // been called. - int64_t num_deltas_computed() const { return num_deltas_computed_; } + // Indicates if the computed jump value is always the best choice. + bool VariableOnlyInLinearConstraintWithConvexViolationChange(int var) const; - // Access the solution stored. Useful for debugging. + // Access the solution stored. const std::vector& current_solution() const { return current_solution_; } @@ -259,10 +341,97 @@ class LsEvaluator { LinearIncrementalEvaluator linear_evaluator_; std::vector> constraints_; std::vector> var_to_constraint_graph_; + std::vector jump_value_optimal_; // We need the mutable to evaluate a move. mutable std::vector current_solution_; bool model_is_supported_ = true; - mutable int64_t num_deltas_computed_ = 0; +}; + +// Individual compiled constraints. + +// The violation of a bool_xor constraint is 0 or 1. +class CompiledBoolXorConstraint : public CompiledConstraint { + public: + explicit CompiledBoolXorConstraint(const ConstraintProto& ct_proto); + ~CompiledBoolXorConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; + int64_t ViolationDelta( + int /*var*/, int64_t /*old_value*/, + absl::Span solution_with_new_value) override; +}; + +// The violation of a lin_max constraint is: +// - the sum(max(0, expr_value - target_value) forall expr). This part will be +// maintained by the linear part. +// - target_value - max(expressions) if positive. +class CompiledLinMaxConstraint : public CompiledConstraint { + public: + explicit CompiledLinMaxConstraint(const ConstraintProto& ct_proto); + ~CompiledLinMaxConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; +}; + +// The violation of an int_prod constraint is +// abs(value(target) - prod(value(expr)). +class CompiledIntProdConstraint : public CompiledConstraint { + public: + explicit CompiledIntProdConstraint(const ConstraintProto& ct_proto); + ~CompiledIntProdConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; +}; + +// The violation of an int_div constraint is +// abs(value(target) - value(expr0) / value(expr1)). +class CompiledIntDivConstraint : public CompiledConstraint { + public: + explicit CompiledIntDivConstraint(const ConstraintProto& ct_proto); + ~CompiledIntDivConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; +}; + +// The violation of a all_diff is the number of unordered pairs of expressions +// with the same value. +class CompiledAllDiffConstraint : public CompiledConstraint { + public: + explicit CompiledAllDiffConstraint(const ConstraintProto& ct_proto); + ~CompiledAllDiffConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; + + private: + std::vector values_; +}; + +// The violation of a no_overlap is the sum of overloads over time. +class CompiledNoOverlapConstraint : public CompiledConstraint { + public: + explicit CompiledNoOverlapConstraint(const ConstraintProto& ct_proto, + const CpModelProto& cp_model); + ~CompiledNoOverlapConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; + + private: + const CpModelProto& cp_model_; + std::vector> events_; +}; + +// The violation of a cumulative is the sum of overloads over time. +class CompiledCumulativeConstraint : public CompiledConstraint { + public: + explicit CompiledCumulativeConstraint(const ConstraintProto& ct_proto, + const CpModelProto& cp_model); + ~CompiledCumulativeConstraint() override = default; + + int64_t ComputeViolation(absl::Span solution) override; + + private: + const CpModelProto& cp_model_; + std::vector> events_; }; } // namespace sat diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc index 61673125a0..c8a962d45f 100644 --- a/ortools/sat/cp_model_postsolve.cc +++ b/ortools/sat/cp_model_postsolve.cc @@ -207,9 +207,8 @@ int64_t EvaluateLinearExpression(const LinearExpressionProto& expr, } // namespace // Compute the max of each expression, and assign it to the target expr (which -// must be of the form +ref or -ref); -// We only support post-solving the case were the target is unassigned, -// but everything else is fixed. +// must be of the form +ref or -ref); We only support post-solving the case +// where all expression are fixed and correct. void PostsolveLinMax(const ConstraintProto& ct, std::vector* domains) { int64_t max_value = std::numeric_limits::min(); for (const LinearExpressionProto& expr : ct.lin_max().exprs()) { @@ -217,9 +216,8 @@ void PostsolveLinMax(const ConstraintProto& ct, std::vector* domains) { } const int target_ref = GetSingleRefFromExpression(ct.lin_max().target()); const int target_var = PositiveRef(target_ref); - (*domains)[target_var] = (*domains)[target_var].IntersectionWith( - Domain(RefIsPositive(target_ref) ? max_value : -max_value)); - CHECK(!(*domains)[target_var].IsEmpty()); + (*domains)[target_var] = + Domain(RefIsPositive(target_ref) ? max_value : -max_value); } // We only support 3 cases in the presolve currently. diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 071ff84260..b440680d17 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -14,11 +14,9 @@ #include "ortools/sat/cp_model_presolve.h" #include -#include #include #include #include -#include #include #include #include @@ -32,14 +30,15 @@ #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" #include "absl/hash/hash.h" +#include "absl/log/check.h" #include "absl/numeric/int128.h" #include "absl/strings/str_cat.h" #include "absl/types/span.h" -#include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/base/stl_util.h" #include "ortools/base/timer.h" +#include "ortools/base/vlog_is_on.h" #include "ortools/graph/topologicalsorter.h" #include "ortools/sat/circuit.h" #include "ortools/sat/clause.h" @@ -869,6 +868,61 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { if (abort) return changed; } + // Checks if the affine target domain is constraining. + bool affine_target_domain_contains_max_domain = false; + if (ExpressionContainsSingleRef(target)) { // target = +/- var. + int64_t infered_min = std::numeric_limits::min(); + int64_t infered_max = std::numeric_limits::min(); + for (const LinearExpressionProto& expr : ct->lin_max().exprs()) { + infered_min = std::max(infered_min, context_->MinOf(expr)); + infered_max = std::max(infered_max, context_->MaxOf(expr)); + } + Domain rhs_domain; + for (const LinearExpressionProto& expr : ct->lin_max().exprs()) { + rhs_domain = rhs_domain.UnionWith( + context_->DomainSuperSetOf(expr).IntersectionWith( + {infered_min, infered_max})); + } + + // Checks if all values from the max(exprs) belong in the domain of the + // target. + // Note that the target is +/-var. + DCHECK_EQ(std::abs(target.coeffs(0)), 1); + const Domain target_domain = + target.coeffs(0) == 1 ? context_->DomainOf(target.vars(0)) + : context_->DomainOf(target.vars(0)).Negation(); + affine_target_domain_contains_max_domain = + rhs_domain.IsIncludedIn(target_domain); + } + + // If the target is not used, and safe, we can remove the constraint. + if (affine_target_domain_contains_max_domain && + context_->VariableIsUniqueAndRemovable(target.vars(0))) { + context_->UpdateRuleStats("lin_max: unused affine target"); + context_->MarkVariableAsRemoved(target.vars(0)); + *context_->mapping_model->add_constraints() = *ct; + return RemoveConstraint(ct); + } + + // If the target is only used in the objective, and safe, we can simplify the + // constraint. + if (affine_target_domain_contains_max_domain && + context_->VariableWithCostIsUniqueAndRemovable(target.vars(0)) && + (target.coeffs(0) > 0) == + (context_->ObjectiveCoeff(target.vars(0)) > 0)) { + context_->UpdateRuleStats("lin_max: rewrite with precedences"); + for (const LinearExpressionProto& expr : ct->lin_max().exprs()) { + LinearConstraintProto* prec = + context_->working_model->add_constraints()->mutable_linear(); + prec->add_domain(0); + prec->add_domain(std::numeric_limits::max()); + AddLinearExpressionToLinearConstraint(target, 1, prec); + AddLinearExpressionToLinearConstraint(expr, -1, prec); + } + *context_->mapping_model->add_constraints() = *ct; + return RemoveConstraint(ct); + } + // Deal with fixed target case. if (target_min == target_max) { bool all_booleans = true; @@ -5111,7 +5165,7 @@ bool CpModelPresolver::PresolveNoOverlap(ConstraintProto* ct) { return changed; } -bool CpModelPresolver::PresolveNoOverlap2D(int c, ConstraintProto* ct) { +bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) { if (context_->ModelIsUnsat()) { return false; } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 41efb156ea..37b0188537 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -40,7 +40,6 @@ #include "absl/container/btree_set.h" #include "absl/container/flat_hash_set.h" #include "absl/flags/flag.h" -#include "absl/status/status.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_format.h" #include "absl/strings/str_join.h" @@ -3202,7 +3201,11 @@ void SolveCpModelParallel(const CpModelProto& model_proto, absl::StrCat("rens_lns_", local_params.name())), local_params, helper, &shared)); } - if (params.use_violation_ls() && !params.interleave_search()) { + + const bool feasibility_jump_possible = + !params.interleave_search() && + helper->TypeToConstraints(ConstraintProto::kNoOverlap2D).empty(); + if (params.use_violation_ls() && feasibility_jump_possible) { SatParameters local_params = params; local_params.set_random_seed(params.random_seed()); local_params.set_feasibility_jump_decay(0.95); @@ -3241,10 +3244,10 @@ void SolveCpModelParallel(const CpModelProto& model_proto, // schedule more than the available number of threads. They will just be // interleaved. We will get an higher diversity, but use more memory. const int num_feasibility_jump = - params.interleave_search() - ? 0 - : (params.test_feasibility_jump() ? num_available - : num_available / 2); + feasibility_jump_possible + ? (params.test_feasibility_jump() ? num_available + : num_available / 2) + : 0; const int num_first_solution_subsolvers = num_available - num_feasibility_jump; for (int i = 0; i < num_feasibility_jump; ++i) { diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc index c738dc186f..221f8acc1f 100644 --- a/ortools/sat/feasibility_jump.cc +++ b/ortools/sat/feasibility_jump.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -23,6 +24,8 @@ #include #include +#include "absl/log/check.h" +#include "absl/random/bit_gen_ref.h" #include "absl/random/random.h" #include "absl/strings/str_cat.h" #include "absl/types/span.h" @@ -32,10 +35,31 @@ #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_checker.h" #include "ortools/sat/cp_model_utils.h" +#include "ortools/sat/integer.h" #include "ortools/sat/restart.h" +#include "ortools/sat/subsolver.h" +#include "ortools/sat/synchronization.h" +#include "ortools/sat/util.h" +#include "ortools/util/sorted_interval_list.h" namespace operations_research::sat { +FeasibilityJumpSolver::~FeasibilityJumpSolver() { + if (!VLOG_IS_ON(1)) return; + std::vector> stats; + stats.push_back({"fs_jump/num_general_moves_computed", num_general_evals_}); + stats.push_back({"fs_jump/num_general_moves_done", num_general_moves_}); + stats.push_back({"fs_jump/num_linear_moves_computed", num_linear_evals_}); + stats.push_back({"fs_jump/num_linear_moves_repaired_with_full_scan", + num_repairs_with_full_scan_}); + stats.push_back({"fs_jump/num_linear_moves_done", num_linear_moves_}); + stats.push_back({"fs_jump/num_solutions_imported", num_solutions_imported_}); + stats.push_back( + {"fs_jump/num_variables_partially_scanned", num_partial_scans_}); + stats.push_back({"fs_jump/num_weight_updates", num_weight_updates_}); + shared_stats_->AddStats(stats); +} + void FeasibilityJumpSolver::Initialize() { is_initialized_ = true; evaluator_ = std::make_unique(cp_model_); @@ -46,18 +70,85 @@ void FeasibilityJumpSolver::Initialize() { // validation. const int num_variables = cp_model_.variables().size(); var_domains_.resize(num_variables); + var_has_two_values_.resize(num_variables); for (int var = 0; var < num_variables; ++var) { var_domains_[var] = ReadDomainFromProto(cp_model_.variables(var)); + var_has_two_values_[var] = var_domains_[var].HasTwoValues(); } } +namespace { + +int64_t ComputeRange(int64_t range, double range_ratio) { + return static_cast( + std::ceil(static_cast(range) * range_ratio)); +} + +// TODO(user): Optimize and move to the Domain class. +// TODO(user): Improve entropy on non continuous domains. +int64_t RandomValueNearMin(const Domain& domain, double range_ratio, + absl::BitGenRef random) { + if (domain.Size() == 1) return domain.FixedValue(); + if (domain.Size() == 2) { + return absl::Bernoulli(random, 1 - range_ratio) ? domain.Min() + : domain.Max(); + } + const int64_t range = ComputeRange(domain.Max() - domain.Min(), range_ratio); + return domain.ValueAtOrBefore(domain.Min() + + absl::LogUniform(random, 0, range)); +} + +int64_t RandomValueNearMax(const Domain& domain, double range_ratio, + absl::BitGenRef random) { + if (domain.Size() == 1) return domain.FixedValue(); + if (domain.Size() == 2) { + return absl::Bernoulli(random, 1 - range_ratio) ? domain.Max() + : domain.Min(); + } + const int64_t range = ComputeRange(domain.Max() - domain.Min(), range_ratio); + return domain.ValueAtOrAfter(domain.Max() - + absl::LogUniform(random, 0, range)); +} + +int64_t RandomValueNearZero(const Domain& domain, double range_ratio, + absl::BitGenRef random) { + if (domain.Min() >= 0) { + return RandomValueNearMin(domain, range_ratio, random); + } + if (domain.Max() <= 0) { + return RandomValueNearMax(domain, range_ratio, random); + } + const Domain positive_domain = domain.IntersectionWith({0, domain.Max()}); + const double choose_positive_probability = + static_cast(positive_domain.Size()) / + static_cast(domain.Size()); + if (absl::Bernoulli(random, choose_positive_probability)) { + return RandomValueNearMin(positive_domain, range_ratio, random); + } else { + return RandomValueNearMax(domain.IntersectionWith({domain.Min(), 0}), + range_ratio, random); + } +} + +} // namespace + void FeasibilityJumpSolver::RestartFromDefaultSolution() { - // Starts with solution closest to zero. - // TODO(user): randomize the start instead? const int num_variables = cp_model_.variables().size(); + const double default_value_probability = + 1.0 - params_.feasibility_jump_var_randomization_ratio(); + const double pertubation_ratio = + params_.feasibility_jump_var_perburbation_range_ratio(); + + // Starts with values closest to zero. std::vector solution(num_variables); for (int var = 0; var < num_variables; ++var) { - solution[var] = var_domains_[var].SmallestValue(); + if (num_batches_ == 0 || + absl::Bernoulli(random_, default_value_probability)) { + solution[var] = var_domains_[var].SmallestValue(); + } else { + solution[var] = + RandomValueNearZero(var_domains_[var], pertubation_ratio, random_); + } } if (cp_model_.has_objective() && @@ -66,9 +157,21 @@ void FeasibilityJumpSolver::RestartFromDefaultSolution() { for (int i = 0; i < num_terms; ++i) { const int var = cp_model_.objective().vars(i); if (cp_model_.objective().coeffs(i) > 0) { - solution[var] = var_domains_[var].Min(); + if (num_batches_ == 0 || + absl::Bernoulli(random_, default_value_probability)) { + solution[var] = var_domains_[var].Min(); + } else { + solution[var] = + RandomValueNearMin(var_domains_[var], pertubation_ratio, random_); + } } else { - solution[var] = var_domains_[var].Max(); + if (num_batches_ == 0 || + absl::Bernoulli(random_, default_value_probability)) { + solution[var] = var_domains_[var].Max(); + } else { + solution[var] = + RandomValueNearMax(var_domains_[var], pertubation_ratio, random_); + } } } } @@ -79,13 +182,33 @@ void FeasibilityJumpSolver::RestartFromDefaultSolution() { } std::string FeasibilityJumpSolver::OneLineStats() const { - return absl::StrCat( - "batch:", num_batches_, " restart:", FormatCounter(num_restarts_), - " inf:", FormatCounter(evaluator_->NumInfeasibleConstraints()), - " improving:", FormatCounter(good_jumps_.size()), - " jumps:", FormatCounter(num_jumps_), - " updates:", FormatCounter(num_weight_updates_), - " non_linear_evals:", FormatCounter(evaluator_->num_deltas_computed())); + // Restarts or solutions imported. + const std::string restart_str = + num_solutions_imported_ == 0 + ? absl::StrCat(" #restarts:", num_restarts_) + : absl::StrCat(" #solutions_imported:", num_solutions_imported_); + + // Moves and evaluations in the general iterations. + const std::string general_str = + num_general_evals_ == 0 && num_general_moves_ == 0 + ? "" + : absl::StrCat(" #gen_moves:", FormatCounter(num_general_moves_), "/", + FormatCounter(num_general_evals_)); + + // Improving jumps and infeasible constraints. + const int num_infeasible_cts = evaluator_->NumInfeasibleConstraints(); + const std::string non_solution_str = + good_jumps_.empty() && num_infeasible_cts == 0 + ? "" + : absl::StrCat(" #good_lin_moves:", FormatCounter(good_jumps_.size()), + " #inf_cts:", + FormatCounter(evaluator_->NumInfeasibleConstraints())); + + return absl::StrCat("batch:", num_batches_, restart_str, + " #lin_moves:", FormatCounter(num_linear_moves_), "/", + FormatCounter(num_linear_evals_), general_str, + non_solution_str, + " #weight_updates:", FormatCounter(num_weight_updates_)); } std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { @@ -105,21 +228,36 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { CHECK_GT(repo.NumSolutions(), 0); const SharedSolutionRepository::Solution solution = repo.GetRandomBiasedSolution(random_); - evaluator_->ComputeInitialViolations(solution.variable_values); + if (solution.rank < last_solution_rank_) { + evaluator_->ComputeInitialViolations(solution.variable_values); - // Reset weights for each new solution. - weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0); + // Reset weights for each new solution. + weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0); + + // Update last solution rank. + last_solution_rank_ = solution.rank; + VLOG(2) << name() << " import a solution with value " << solution.rank; + ++num_solutions_imported_; + } } else { // Restart? Note that we always "restart" the first time. const double dtime = evaluator_->MutableLinearEvaluator()->DeterministicTime(); - if (dtime >= dtime_restart_threshold_) { + if (dtime >= dtime_restart_threshold_ && + num_weight_updates_ >= update_restart_threshold_) { ++num_restarts_; RestartFromDefaultSolution(); - // We use luby restart with a base of 0.5s. + // We use luby restart with a base of 1 deterministic unit. + // We also block the restart if there was not enough weight update. // Note that we only restart between batches too. - dtime_restart_threshold_ = dtime + 0.5 * SUniv(num_restarts_); + // + // TODO(user): Ideally batch should use deterministic time too so we + // can just use number of batch for the luby restart. + // TODO(user): Maybe have one worker with very low restart + // rate. + dtime_restart_threshold_ = dtime + 1.0 * SUniv(num_restarts_); + update_restart_threshold_ = num_weight_updates_ + 10; } } @@ -133,6 +271,9 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { // chunks. if (shared_bounds_ != nullptr) { shared_bounds_->UpdateDomains(&var_domains_); + for (int var = 0; var < var_domains_.size(); ++var) { + var_has_two_values_[var] = var_domains_[var].HasTwoValues(); + } } // Checks the current solution is compatible with updated domains. @@ -191,21 +332,21 @@ bool IsGood(double delta) { return delta < 0.0; } void FeasibilityJumpSolver::RecomputeJump(int var) { const std::vector& solution = evaluator_->current_solution(); - ++num_computed_jumps_; + ++num_linear_evals_; jump_need_recomputation_[var] = false; if (var_domains_[var].IsFixed()) { jump_values_[var] = solution[var]; - jump_deltas_[var] = 0.0; + jump_scores_[var] = 0.0; return; } LinearIncrementalEvaluator* linear_evaluator = evaluator_->MutableLinearEvaluator(); - if (var_domains_[var].Size() == 2) { + if (var_has_two_values_[var]) { jump_values_[var] = solution[var] == var_domains_[var].Min() ? var_domains_[var].Max() : var_domains_[var].Min(); - jump_deltas_[var] = linear_evaluator->WeightedViolationDelta( + jump_scores_[var] = linear_evaluator->WeightedViolationDelta( weights_, var, jump_values_[var] - solution[var]); } else { // In practice, after a few iterations, the chance of finding an improving @@ -274,21 +415,10 @@ void FeasibilityJumpSolver::RecomputeJump(int var) { DCHECK_NE(best_jump.first, solution[var]); jump_values_[var] = best_jump.first; - jump_deltas_[var] = best_jump.second; + jump_scores_[var] = best_jump.second; } - if (IsGood(jump_deltas_[var]) && !in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); - } -} - -void FeasibilityJumpSolver::MarkJumpForRecomputation(int var) { - jump_need_recomputation_[var] = true; - - // This jump might be good, so we need to add it to the queue so it can be - // evaluated when choosing the next jump. - if (!in_good_jumps_[var]) { + if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) { in_good_jumps_[var] = true; good_jumps_.push_back(var); } @@ -297,7 +427,7 @@ void FeasibilityJumpSolver::MarkJumpForRecomputation(int var) { void FeasibilityJumpSolver::RecomputeAllJumps() { const int num_variables = static_cast(var_domains_.size()); jump_values_.resize(num_variables); - jump_deltas_.resize(num_variables); + jump_scores_.resize(num_variables); jump_need_recomputation_.assign(num_variables, true); in_good_jumps_.assign(num_variables, false); @@ -313,69 +443,91 @@ int FeasibilityJumpSolver::UpdateConstraintWeights(bool linear_mode) { const double kMaxWeight = 1e50; const double kBumpFactor = 1.0 / params_.feasibility_jump_decay(); + LinearIncrementalEvaluator* linear_evaluator = + evaluator_->MutableLinearEvaluator(); bump_value_ *= kBumpFactor; bool rescale = false; int num_bumps = 0; - tmp_constraints_.clear(); const int num_constraints = linear_mode ? evaluator_->NumLinearConstraints() : evaluator_->NumEvaluatorConstraints(); + linear_evaluator->ClearAffectedVariables(); for (int c = 0; c < num_constraints; ++c) { - if (evaluator_->Violation(c) > 0) { - ++num_bumps; - weights_[c] += bump_value_; - if (weights_[c] > kMaxWeight) rescale = true; - tmp_constraints_.push_back(c); + if (!evaluator_->IsViolated(c)) continue; + + ++num_bumps; + weights_[c] += bump_value_; + if (weights_[c] > kMaxWeight) rescale = true; + if (!rescale) { + linear_evaluator->UpdateScoreOnWeightUpdate( + c, bump_value_, evaluator_->current_solution(), jump_values_, + absl::MakeSpan(jump_scores_)); } } - LinearIncrementalEvaluator* linear_evaluator = - evaluator_->MutableLinearEvaluator(); - if (rescale) { const double factor = 1.0 / kMaxWeight; bump_value_ *= factor; for (int c = 0; c < num_constraints; ++c) { weights_[c] *= factor; } - if (linear_mode) { - solution_score_ = linear_evaluator->WeightedViolation(weights_); - RecomputeAllJumps(); - return num_bumps; - } + RecomputeAllJumps(); + return num_bumps; } - // Recompute the affected jumps and overall score. + // Recompute the affected jumps. // Note that the constraint violations are unaffected. - if (linear_mode) { - solution_score_ = linear_evaluator->WeightedViolation(weights_); - for (const int var : - linear_evaluator->ConstraintsToAffectedVariables(tmp_constraints_)) { - MarkJumpForRecomputation(var); + for (const int var : linear_evaluator->VariablesAffectedByLastUpdate()) { + // We don't need to recompute score of binary variable, it should + // already be correct. + if (!jump_need_recomputation_[var] && var_has_two_values_[var]) { + DCHECK(JumpIsUpToDate(var)); + if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) { + in_good_jumps_[var] = true; + good_jumps_.push_back(var); + } + continue; + } + + // This jump might be good, so we need to add it to the queue so it can be + // evaluated when choosing the next jump. + jump_need_recomputation_[var] = true; + if (!in_good_jumps_[var]) { + in_good_jumps_[var] = true; + good_jumps_.push_back(var); } - } else { - solution_score_ = evaluator_->WeightedViolation(weights_); } return num_bumps; } +// Important: This is for debugging, but unfortunately it currently change the +// deterministic time and hence the overall algo behavior. +bool FeasibilityJumpSolver::JumpIsUpToDate(int var) { + const int64_t old_jump = jump_values_[var]; + const double old_score = jump_scores_[var]; + RecomputeJump(var); + CHECK_EQ(jump_values_[var], old_jump); // No change + const double relative = + std::max({std::abs(jump_scores_[var]), std::abs(old_score), 1.0}); + return std::abs(jump_scores_[var] - old_score) / relative < 1e-6; +} + bool FeasibilityJumpSolver::DoSomeLinearIterations() { ++num_batches_; // Recompute at the beginning of each batch. LinearIncrementalEvaluator* linear_evaluator = evaluator_->MutableLinearEvaluator(); - solution_score_ = linear_evaluator->WeightedViolation(weights_); const std::vector& solution = evaluator_->current_solution(); RecomputeAllJumps(); if (VLOG_IS_ON(1)) { if (num_batches_ < 10) { // The first few batches are more informative to understand how the - // heuristic behave. + // heuristic behaves. shared_response_->LogMessage(name(), OneLineStats()); } else { shared_response_->LogPeriodicMessage(name(), OneLineStats(), @@ -403,7 +555,9 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { // Take the best jump score amongst some random candidates. // It is okay if we pick twice the same, we don't really care. int best_var = -1; - double best_delta = 0.0; + int best_index = -1; + int64_t best_value = 0; + double best_score = 0.0; int num_improving_jump_tested = 0; while (!good_jumps_.empty() && num_improving_jump_tested < 5) { const int index = absl::Uniform( @@ -411,52 +565,78 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { const int var = good_jumps_[index]; // We lazily update the jump value. - if (jump_need_recomputation_[var]) RecomputeJump(var); + if (jump_need_recomputation_[var]) { + RecomputeJump(var); + } else { + DCHECK(JumpIsUpToDate(var)); + } - if (!IsGood(jump_deltas_[var])) { + if (!IsGood(jump_scores_[var])) { // Lazily remove. in_good_jumps_[var] = false; good_jumps_[index] = good_jumps_.back(); good_jumps_.pop_back(); + if (best_index == good_jumps_.size()) best_index = index; continue; } ++num_improving_jump_tested; - if (jump_deltas_[var] <= best_delta) { + if (jump_scores_[var] <= best_score) { best_var = var; - best_delta = jump_deltas_[var]; + best_index = index; + best_value = jump_values_[var]; + best_score = jump_scores_[var]; } } if (good_jumps_.empty()) { - break; + // TODO(user): Re-enable the code. It can be a bit slow currently + // especially since in many situation it doesn't achieve anything. + if (true) break; + + bool time_limit_crossed = false; + if (ScanAllVariables(solution, /*linear_mode=*/true, &best_var, + &best_value, &best_score, &time_limit_crossed)) { + VLOG(2) << name() + << " scans and finds improving solution (var:" << best_var + << " value:" << best_value << ", delta:" << best_score << ")"; + ++num_repairs_with_full_scan_; + } else if (time_limit_crossed) { + return false; + } else { + break; + } + } else { + const int64_t var_change = best_value - solution[best_var]; + DCHECK_EQ(best_score, linear_evaluator->WeightedViolationDelta( + weights_, best_var, var_change)); } - const int64_t var_change = jump_values_[best_var] - solution[best_var]; - DCHECK_EQ(best_delta, linear_evaluator->WeightedViolationDelta( - weights_, best_var, var_change)); + CHECK_NE(best_var, -1); + CHECK_NE(best_index, -1); // Perform the move. - ++num_jumps_; - CHECK_NE(best_var, -1); - evaluator_->UpdateVariableValueAndRecomputeViolations( - best_var, jump_values_[best_var], /*focus_on_linear=*/true); - solution_score_ += best_delta; - DCHECK_EQ(solution_score_, linear_evaluator->WeightedViolation(weights_)); + ++num_linear_moves_; + const int64_t old_value = evaluator_->current_solution()[best_var]; + linear_evaluator->ClearAffectedVariables(); + evaluator_->UpdateLinearScores(best_var, best_value, weights_, + jump_values_, + absl::MakeSpan(jump_scores_)); + evaluator_->UpdateVariableValue(best_var, best_value); - // For now we recompute all the affected best jump value. - // - // TODO(user): In the paper, they just recompute the scores and only - // change the jump values when the constraint weight change. This should - // be faster. - // - // TODO(user): We can reduce the complexity if we just recompute the - // score. we don't need to rescan each column of the affected variables - // again, but could update their scores directly. - for (const int var : - linear_evaluator->AffectedVariableOnUpdate(best_var)) { - MarkJumpForRecomputation(var); + // We already know the score of undoing the move we just did, and we know + // this move is bad, so we can remove it from good_jumps_ right away. + jump_values_[best_var] = old_value; + jump_scores_[best_var] = -best_score; + if (var_has_two_values_[best_var]) { + CHECK_EQ(good_jumps_[best_index], best_var); + in_good_jumps_[best_var] = false; + good_jumps_[best_index] = good_jumps_.back(); + good_jumps_.pop_back(); + } else { + jump_need_recomputation_[best_var] = true; } + MarkJumpsThatNeedsToBeRecomputed(best_var); } // We will update the weight unless the queue is non-empty. @@ -468,52 +648,90 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { return false; } +// Update the jump scores. +// +// We incrementally maintain the score (except for best_var). +// However for non-Boolean, we still need to recompute the jump value. +// We will do that in a lazy fashion. +// +// TODO(user): In the paper, they just recompute the scores and only +// change the jump values when the constraint weight changes. Experiment? +// Note however that the current code is quite fast. +// +// TODO(user): For non-Boolean, we could easily detect if a non-improving +// score cannot become improving. We don't need to add such variable to +// the queue. +void FeasibilityJumpSolver::MarkJumpsThatNeedsToBeRecomputed(int changed_var) { + LinearIncrementalEvaluator* linear_evaluator = + evaluator_->MutableLinearEvaluator(); + for (const int var : linear_evaluator->VariablesAffectedByLastUpdate()) { + if (var == changed_var) continue; + if (jump_need_recomputation_[var]) { + DCHECK(in_good_jumps_[var]); + continue; + } + + // We don't need to recompute score of binary variable, it should + // already be correct. + if (var_has_two_values_[var]) { + DCHECK(JumpIsUpToDate(var)); + if (IsGood(jump_scores_[var]) && !in_good_jumps_[var]) { + in_good_jumps_[var] = true; + good_jumps_.push_back(var); + } + continue; + } + + jump_need_recomputation_[var] = true; + if (!in_good_jumps_[var]) { + in_good_jumps_[var] = true; + good_jumps_.push_back(var); + } + } +} + bool FeasibilityJumpSolver::DoSomeGeneralIterations() { if (evaluator_->NumNonLinearConstraints() == 0) { return true; } + LinearIncrementalEvaluator* linear_evaluator = + evaluator_->MutableLinearEvaluator(); const std::vector& solution = evaluator_->current_solution(); - // Non linear constraints are not evaluated in the linear phase. - evaluator_->UpdateNonLinearViolations(); - solution_score_ = evaluator_->WeightedViolation(weights_); + // Non-linear constraints are not evaluated in the linear phase. + evaluator_->UpdateAllNonLinearViolations(); const int kLoops = 10000; - std::vector to_scan; for (int loop = 0; loop < kLoops; ++loop) { - to_scan = evaluator_->VariablesInViolatedConstraints(); - std::shuffle(to_scan.begin(), to_scan.end(), random_); - bool improving_move_found = false; + int var; + int64_t value; + double score; + bool time_limit_crossed = false; + if (ScanAllVariables(solution, /*linear_mode=*/false, &var, &value, &score, + &time_limit_crossed)) { + // Perform the move. + num_general_moves_++; - for (const int current_var : to_scan) { - const int64_t initial_value = solution[current_var]; - int64_t best_value = initial_value; - double best_delta = 0.0; - for (const int64_t new_value : var_domains_[current_var].Values()) { - if (new_value == initial_value) continue; - const double delta = evaluator_->WeightedViolationDelta( - weights_, current_var, new_value - initial_value); - if (evaluator_->num_deltas_computed() % 1000 == 0 && - shared_time_limit_ != nullptr && - shared_time_limit_->LimitReached()) { - return false; - } - if (delta < best_delta) { - improving_move_found = true; - best_value = new_value; - best_delta = delta; + // Update the linear part. + linear_evaluator->ClearAffectedVariables(); + evaluator_->UpdateLinearScores(var, value, weights_, jump_values_, + absl::MakeSpan(jump_scores_)); + jump_values_[var] = solution[var]; + jump_scores_[var] = -score; + for (const int v : linear_evaluator->VariablesAffectedByLastUpdate()) { + if (!var_has_two_values_[var]) { + jump_need_recomputation_[v] = true; } } - if (improving_move_found) { // Accept the move with the best delta. - solution_score_ += best_delta; - evaluator_->UpdateVariableValueAndRecomputeViolations(current_var, - best_value); - break; - } + // Update the non-linear part. Note it also commits the move. + evaluator_->UpdateNonLinearViolations(var, value); + evaluator_->UpdateVariableValue(var, value); + continue; + } else if (time_limit_crossed) { + return false; } - if (UpdateConstraintWeights(/*linear_mode=*/false) == 0) { return true; } @@ -521,4 +739,142 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { return false; } +bool FeasibilityJumpSolver::ScanAllVariables(absl::Span solution, + bool linear_mode, + int* improving_var, + int64_t* improving_value, + double* improving_delta, + bool* time_limit_crossed) { + LinearIncrementalEvaluator* linear_evaluator = + evaluator_->MutableLinearEvaluator(); + + // We stop at the first improving move that does not degrade the linear score. + // Otherwise, we return the best improving move. + tmp_to_scan_ = evaluator_->VariablesInViolatedConstraints(); + std::shuffle(tmp_to_scan_.begin(), tmp_to_scan_.end(), random_); + + // We maintain the best solution found, regardless of the linear violations. + int best_var = -1; + int64_t best_value = 0; + double best_delta = -1e-4; + // We will favor moves that keep the linear assignment feasible for as long as + // possible (if feasible to starts with). + // + // TODO(user): We just need to check if the assignment is linear + // feasible, not compute its score. + const bool accept_moves_with_linear_violation = + !params_.feasibility_jump_protect_linear_feasibility() || + linear_evaluator->WeightedViolation(weights_) > 0.0; + + for (const int var : tmp_to_scan_) { + if (evaluator_->VariableOnlyInLinearConstraintWithConvexViolationChange( + var)) { + // We lazily update the jump value. + if (jump_need_recomputation_[var]) { + RecomputeJump(var); + } else { + DCHECK(JumpIsUpToDate(var)); + } + + const double delta = jump_scores_[var]; + if (IsGood(delta)) { + *improving_var = var; + *improving_value = jump_values_[var]; + *improving_delta = delta; + return true; + } + continue; + } + + // Compute the values to scan. + // If the domain is large we do not scan the full variable domain. + Domain scan_range; + const int64_t initial_value = solution[var]; + const Domain& var_domain = var_domains_[var]; + if (var_domains_[var].Size() <= + params_.feasibility_jump_max_num_values_scanned()) { + scan_range = var_domain; + } else { + // TODO(user): Partial scan of the rest of the domain. + ++num_partial_scans_; + const int64_t half_span = + params_.feasibility_jump_max_num_values_scanned() / 2; + Domain scan_range = + Domain(initial_value).AdditionWith({-half_span, half_span}); + if (scan_range.Min() < var_domain.Min()) { + scan_range = scan_range.AdditionWith( + Domain(var_domain.Min() - scan_range.Min())); + } else if (scan_range.Max() > var_domain.Max()) { + scan_range = scan_range.AdditionWith( + Domain(scan_range.Max() - var_domain.Max())); + } + scan_range = scan_range.IntersectionWith(var_domain); + } + + // Take the best amongst all the values to scan for the current variable if + // it does not degrade the linear score (when + // accept_moves_with_linear_violation is false). + int best_var_non_increasing_linear = -1; + int64_t best_value_non_increasing_linear = 0; + double best_delta_non_increasing_linear = -1e-4; + + for (const int64_t new_value : scan_range.Values()) { + if (new_value == initial_value) continue; + + // Checks the time limit periodically. + if (++num_general_evals_ % 1000 == 0 && shared_time_limit_ != nullptr && + shared_time_limit_->LimitReached()) { + *time_limit_crossed = true; + return false; + } + + // Computes both linear and general delta. + const double linear_delta = linear_evaluator->WeightedViolationDelta( + weights_, var, new_value - initial_value); + const double non_linear_delta = + linear_mode ? 0.0 + : evaluator_->WeightedNonLinearViolationDelta( + weights_, var, new_value - initial_value); + const double delta = linear_delta + non_linear_delta; + + // Do we have a valid improving move ? We select the best value for the + // first variable found with such a move so we delay the exit until after + // the loop on values is finished. + if ((accept_moves_with_linear_violation || linear_delta <= 0.0) && + delta < best_delta_non_increasing_linear) { + best_var_non_increasing_linear = var; + best_value_non_increasing_linear = new_value; + best_delta_non_increasing_linear = delta; + } + + // We keep the best move so far. + if (delta < best_delta) { + best_var = var; + best_value = new_value; + best_delta = delta; + } + } // End scan values. + + // Early accept the best improving move that does not degrade the linear + // score, or the best move if accept_moves_with_linear_violation is true. + if (best_var_non_increasing_linear != -1) { + *improving_var = best_var_non_increasing_linear; + *improving_value = best_value_non_increasing_linear; + *improving_delta = best_delta_non_increasing_linear; + return true; + } + } + + // Accept the best improving move. + if (best_var != -1) { + *improving_var = best_var; + *improving_value = best_value; + *improving_delta = best_delta; + return true; + } + + // No solution found. + return false; +} + } // namespace operations_research::sat diff --git a/ortools/sat/feasibility_jump.h b/ortools/sat/feasibility_jump.h index 8e16b51226..b53b3a1751 100644 --- a/ortools/sat/feasibility_jump.h +++ b/ortools/sat/feasibility_jump.h @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -24,9 +25,9 @@ #include "absl/time/time.h" #include "absl/types/span.h" +#include "ortools/base/vlog_is_on.h" #include "ortools/sat/constraint_violation.h" #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/integer.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/subsolver.h" #include "ortools/sat/synchronization.h" @@ -64,14 +65,8 @@ class FeasibilityJumpSolver : public SubSolver { shared_stats_(shared_stats), random_(params_) {} - ~FeasibilityJumpSolver() override { - if (!VLOG_IS_ON(1)) return; - std::vector> stats; - stats.push_back({"fs_jump/num_jumps", num_jumps_}); - stats.push_back({"fs_jump/num_computed_jumps", num_computed_jumps_}); - stats.push_back({"fs_jump/num_updates", num_weight_updates_}); - shared_stats_->AddStats(stats); - } + // If VLOG_IS_ON(1), it will export a bunch of statistics. + ~FeasibilityJumpSolver() override; // No synchronization needed for TaskIsAvailable(). void Synchronize() final {} @@ -105,14 +100,20 @@ class FeasibilityJumpSolver : public SubSolver { std::string OneLineStats() const; // Linear only. + bool JumpIsUpToDate(int var); // Debug. void RecomputeJump(int var); - void MarkJumpForRecomputation(int var); void RecomputeAllJumps(); + void MarkJumpsThatNeedsToBeRecomputed(int changed_var); // Moves. bool DoSomeLinearIterations(); bool DoSomeGeneralIterations(); + // Returns true if an improving move was found. + bool ScanAllVariables(absl::Span solution, bool linear_mode, + int* improving_var, int64_t* improving_value, + double* improving_delta, bool* time_limit_crossed); + // Return the number of infeasible constraints. int UpdateConstraintWeights(bool linear_mode); @@ -134,6 +135,7 @@ class FeasibilityJumpSolver : public SubSolver { std::unique_ptr evaluator_; std::vector var_domains_; + std::vector var_has_two_values_; // For each variable, we store: // - A jump value, which is different from the current solution, except for @@ -142,37 +144,46 @@ class FeasibilityJumpSolver : public SubSolver { // jump. std::vector jump_need_recomputation_; std::vector jump_values_; - std::vector jump_deltas_; + std::vector jump_scores_; // The score of a solution is just the sum of infeasibility of each - // constraints weighted by these scores. + // constraint weighted by these scores. std::vector weights_; - // The current weighted violation (lower is better). - double solution_score_; - // Depending on the options, we use an exponentially decaying constraint // weight like for SAT activities. double bump_value_ = 1.0; - // Sparse list of jump with a potential delta <= 0.0. - // We lazily recompute the exact delta as we randomly pick variable from here. + // Sparse list of jumps with a potential delta < 0.0. + // If jump_need_recomputation_[var] is true, we lazily recompute the exact + // delta as we randomly pick variables from here. std::vector in_good_jumps_; std::vector good_jumps_; // We restart each time our local deterministic time crosses this. double dtime_restart_threshold_ = 0.0; + int64_t update_restart_threshold_ = 0; - std::vector tmp_constraints_; std::vector tmp_breakpoints_; // Statistics absl::Time last_logging_time_; - int64_t num_restarts_ = 0; int64_t num_batches_ = 0; - int64_t num_jumps_ = 0; - int64_t num_computed_jumps_ = 0; + int64_t num_linear_evals_ = 0; + int64_t num_general_evals_ = 0; + int64_t num_general_moves_ = 0; + int64_t num_linear_moves_ = 0; + int64_t num_partial_scans_ = 0; + int64_t num_repairs_with_full_scan_ = 0; + int64_t num_restarts_ = 0; + int64_t num_solutions_imported_ = 0; int64_t num_weight_updates_ = 0; + + // Temporary storage. + std::vector tmp_to_scan_; + + // Info on the last solution loaded. + int64_t last_solution_rank_ = std::numeric_limits::max(); }; } // namespace operations_research::sat diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 7a904e30be..8309bb8ac4 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -107,7 +107,7 @@ void LinearConstraintManager::RescaleActiveCounts(const double scaling_factor) { constraint_infos_[i].active_count *= scaling_factor; } constraint_active_count_increase_ *= scaling_factor; - VLOG(2) << "Rescaled active counts by " << scaling_factor; + VLOG(3) << "Rescaled active counts by " << scaling_factor; } bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints( @@ -148,7 +148,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints( lp_constraints_.resize(new_size); solution_state->statuses.resize(num_cols + glop::ColIndex(new_size)); if (num_removed_constraints > 0) { - VLOG(2) << "Removed " << num_removed_constraints << " constraints"; + VLOG(3) << "Removed " << num_removed_constraints << " constraints"; } return num_removed_constraints > 0; } @@ -171,18 +171,22 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add( const ConstraintIndex ct_index = equiv_constraints_[key]; if (constraint_infos_[ct_index].constraint.vars == ct.vars && constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) { - if (added != nullptr) *added = false; + bool tightened = false; if (ct.lb > constraint_infos_[ct_index].constraint.lb) { + tightened = true; if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true; constraint_infos_[ct_index].constraint.lb = ct.lb; - if (added != nullptr) *added = true; } if (ct.ub < constraint_infos_[ct_index].constraint.ub) { + tightened = true; if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true; constraint_infos_[ct_index].constraint.ub = ct.ub; - if (added != nullptr) *added = true; } - ++num_merged_constraints_; + if (added != nullptr) *added = tightened; + if (tightened) { + ++num_merged_constraints_; + FillActivityRange(&constraint_infos_[ct_index]); + } return ct_index; } } @@ -192,6 +196,7 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add( ConstraintInfo ct_info; ct_info.constraint = std::move(ct); ct_info.l2_norm = ComputeL2Norm(ct_info.constraint); + FillActivityRange(&ct_info); ct_info.hash = key; equiv_constraints_[key] = ct_index; ct_info.active_count = constraint_active_count_increase_; @@ -246,7 +251,7 @@ bool LinearConstraintManager::AddCut( // Only add cut with sufficient efficacy. if (violation / l2_norm < 1e-5) { - VLOG(2) << "BAD Cut '" << type_name << "'" + VLOG(3) << "BAD Cut '" << type_name << "'" << " size=" << ct.vars.size() << " max_magnitude=" << ComputeInfinityNorm(ct) << " norm=" << l2_norm << " violation=" << violation @@ -328,7 +333,7 @@ void LinearConstraintManager::PermanentlyRemoveSomeConstraints() { } if (num_deleted_constraints > 0) { - VLOG(2) << "Constraint manager cleanup: #deleted:" + VLOG(3) << "Constraint manager cleanup: #deleted:" << num_deleted_constraints; } num_deletable_constraints_ -= num_deleted_constraints; @@ -509,6 +514,28 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { return term_changed; } +void LinearConstraintManager::FillActivityRange(ConstraintInfo* info) { + IntegerValue min_sum(0); + IntegerValue max_sum(0); + const int num_terms = info->constraint.vars.size(); + for (int i = 0; i < num_terms; ++i) { + const IntegerVariable var = info->constraint.vars[i]; + const IntegerValue coeff = info->constraint.coeffs[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + if (coeff > 0.0) { + min_sum += coeff * lb; + max_sum += coeff * ub; + } else { + min_sum += coeff * ub; + max_sum += coeff * lb; + } + } + const IntegerValue tight_lb = std::max(min_sum, info->constraint.lb); + const IntegerValue tight_ub = std::min(max_sum, info->constraint.ub); + info->activity_range = CapSub(tight_ub.value(), tight_lb.value()); +} + bool LinearConstraintManager::ChangeLp( const absl::StrongVector& lp_solution, glop::BasisState* solution_state, int* num_new_constraints) { diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h index fc1d1c32ff..581b920596 100644 --- a/ortools/sat/linear_constraint_manager.h +++ b/ortools/sat/linear_constraint_manager.h @@ -52,6 +52,12 @@ class LinearConstraintManager { public: struct ConstraintInfo { LinearConstraint constraint; + + // For a boxed constraint, this is the constraint rhs - lhs. For a one sided + // constraint, we first use trivial activity bound to compute tight rhs and + // lhs and then use the same formula. + IntegerValue activity_range; + double l2_norm = 0.0; int64_t inactive_count = 0; double objective_parallelism = 0.0; @@ -169,7 +175,7 @@ class LinearConstraintManager { // Helper method to compute objective parallelism for a given constraint. This // also lazily computes objective norm. - void ComputeObjectiveParallelism(const ConstraintIndex ct_index); + void ComputeObjectiveParallelism(ConstraintIndex ct_index); // Multiplies all active counts and the increment counter by the given // 'scaling_factor'. This should be called when at least one of the active @@ -180,6 +186,8 @@ class LinearConstraintManager { // don't remove any constraints which are already in LP. void PermanentlyRemoveSomeConstraints(); + void FillActivityRange(ConstraintInfo* info); + const SatParameters& sat_parameters_; const IntegerTrail& integer_trail_; diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index b24ca3c945..9a52593719 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -292,6 +292,7 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() { // Fill integer_lp_. integer_lp_.clear(); infinity_norms_.clear(); + ct_bound_diff_.clear(); const auto& all_constraints = constraint_manager_.AllConstraints(); for (const auto index : constraint_manager_.LpConstraints()) { const LinearConstraint& ct = all_constraints[index].constraint; @@ -321,8 +322,9 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() { new_ct.terms.push_back({GetMirrorVariable(var), coeff}); } infinity_norms_.push_back(infinity_norm); + ct_bound_diff_.push_back(all_constraints[index].activity_range); - // Important to keep lp_data_ "clean". + // It is important to keep lp_data_ "clean". DCHECK(std::is_sorted(new_ct.terms.begin(), new_ct.terms.end())); } @@ -889,8 +891,10 @@ bool LinearProgrammingConstraint::PreprocessCut(LinearConstraint* cut) { const IntegerValue slack{CapSub(cut->ub.value(), min_sum.value())}; if (!min_sum_overflow && !AtMinOrMaxInt64(slack.value())) { - // TODO(user): raise conflict or report UNSAT. - if (slack < 0) return false; // Always false. + if (slack < 0) { + problem_proven_infeasible_by_cuts_ = true; + return false; + } if (trail_->CurrentDecisionLevel() == 0 && max_range > slack) { bool newly_fixed = false; @@ -936,6 +940,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( IntegerValue cut_ub; if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_, &cut_ub)) { + ++num_cut_overflows_; VLOG(1) << "Issue, overflow!"; return false; } @@ -988,6 +993,8 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( // TODO(user): Keep track of the potential overflow here. if (!base_ct_.FillFromLinearConstraint(cut_, expanded_lp_solution_, integer_trail_)) { + ++num_cut_overflows_; + VLOG(1) << "Issue, overflow!"; return false; } @@ -1019,8 +1026,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( CutTerm entry; entry.coeff = IntTypeAbs(coeff); entry.lp_value = 0.0; - entry.bound_diff = - CapSub(integer_lp_[row].ub.value(), integer_lp_[row].lb.value()); + entry.bound_diff = ct_bound_diff_[row]; entry.expr_vars[0] = first_slack + 2 * IntegerVariable(tmp_slack_rows_.size()); entry.expr_coeffs[1] = 0; @@ -1650,10 +1656,14 @@ bool LinearProgrammingConstraint::Propagate() { // TODO(user): Refactor so that they are just normal cut generators? const int level = trail_->CurrentDecisionLevel(); if (trail_->CurrentDecisionLevel() == 0) { + problem_proven_infeasible_by_cuts_ = false; if (parameters_.add_objective_cut()) AddObjectiveCut(); if (parameters_.add_mir_cuts()) AddMirCuts(); if (parameters_.add_cg_cuts()) AddCGCuts(); if (parameters_.add_zero_half_cuts()) AddZeroHalfCuts(); + if (problem_proven_infeasible_by_cuts_) { + return integer_trail_->ReportConflict({}); + } } // Try to add cuts. @@ -2316,6 +2326,11 @@ bool LinearProgrammingConstraint::ExactLpReasonning() { tmp_lp_multipliers_.push_back({row, value}); } + // In this case, the LP lower bound match the basic objective "constraint" + // propagation. That is there is an LP solution with all objective variable at + // their current best bound. There is no need to do more work here. + if (tmp_lp_multipliers_.empty()) return true; + IntegerValue scaling = 0; tmp_integer_multipliers_ = ScaleLpMultiplier( /*take_objective_into_account=*/true, tmp_lp_multipliers_, &scaling); @@ -2613,8 +2628,8 @@ std::string LinearProgrammingConstraint::Statistics() const { FormatCounter(total_num_simplex_iterations_), "\n"); absl::StrAppend(&result, " total num cut propagation: ", FormatCounter(total_num_cut_propagations_), "\n"); - absl::StrAppend(&result, " total num prevent overflow: ", - FormatCounter(num_prevent_overflows_), "\n"); + absl::StrAppend(&result, " total num cut overflow: ", + FormatCounter(num_cut_overflows_), "\n"); absl::StrAppend(&result, " total num adjust: ", FormatCounter(num_adjusts_), "\n"); absl::StrAppend(&result, " total num scaling issues: ", diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 3819698bd6..902861df5c 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -438,6 +438,7 @@ class LinearProgrammingConstraint : public PropagatorInterface, IntegerValue objective_infinity_norm_ = IntegerValue(0); absl::StrongVector integer_lp_; absl::StrongVector infinity_norms_; + absl::StrongVector ct_bound_diff_; // Underlying LP solver API. glop::GlopParameters simplex_params_; @@ -455,6 +456,7 @@ class LinearProgrammingConstraint : public PropagatorInterface, FlowCoverCutHelper flow_cover_cut_helper_; IntegerRoundingCutHelper integer_rounding_cut_helper_; + bool problem_proven_infeasible_by_cuts_ = false; CutData base_ct_; LinearConstraint cut_; LinearConstraint saved_cut_; @@ -580,7 +582,7 @@ class LinearProgrammingConstraint : public PropagatorInterface, // Some stats on the LP statuses encountered. int64_t num_solves_ = 0; mutable int64_t num_adjusts_ = 0; - mutable int64_t num_prevent_overflows_ = 0; + mutable int64_t num_cut_overflows_ = 0; mutable int64_t num_scaling_issues_ = 0; std::vector num_solves_by_status_; }; diff --git a/ortools/sat/parameters_validation.cc b/ortools/sat/parameters_validation.cc index 5c49179cbc..e2b087a975 100644 --- a/ortools/sat/parameters_validation.cc +++ b/ortools/sat/parameters_validation.cc @@ -13,6 +13,8 @@ #include "ortools/sat/parameters_validation.h" +#include +#include #include #include @@ -89,7 +91,6 @@ std::string ValidateParameters(const SatParameters& params) { TEST_NOT_NAN(max_time_in_seconds); TEST_NOT_NAN(max_deterministic_time); - TEST_NOT_NAN(feasibility_jump_decay); // TODO(user): Consider using annotations directly in the proto for these // validation. It is however not open sourced. @@ -97,7 +98,15 @@ std::string ValidateParameters(const SatParameters& params) { TEST_IN_RANGE(mip_max_bound, 0, 1e17); TEST_IN_RANGE(solution_pool_size, 1, std::numeric_limits::max()); TEST_IN_RANGE(shared_tree_worker_objective_split_probability, 0.0, 1.0); + + // Feasibility jump. + TEST_NOT_NAN(feasibility_jump_decay); + TEST_NOT_NAN(feasibility_jump_var_randomization_ratio); + TEST_NOT_NAN(feasibility_jump_var_perburbation_range_ratio); + TEST_IN_RANGE(feasibility_jump_max_num_values_scanned, 2, 1'000'000'000); TEST_IN_RANGE(feasibility_jump_decay, 0.0, 1.0); + TEST_IN_RANGE(feasibility_jump_var_randomization_ratio, 0.0, 1.0); + TEST_IN_RANGE(feasibility_jump_var_perburbation_range_ratio, 0.0, 1.0); TEST_POSITIVE(glucose_decay_increment_period); TEST_POSITIVE(shared_tree_max_nodes_per_worker); diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index e8a1a4fdbe..819b29ff7b 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -1916,6 +1916,18 @@ bool PresolveContext::ShiftCostInExactlyOne(absl::Span exactly_one, // Note that the domain never include the offset, so we need to update it. if (offset != 0) AddToObjectiveOffset(offset); + + // When we shift the cost using an exactly one, our objective implied bounds + // might be more or less precise. If the objective domain is not constraining + // (and thus just constraining the upper bound), we relax it to make sure its + // stay "non constraining". + // + // TODO(user): This is a bit hacky, find a nicer way. + if (!objective_domain_is_constraining_) { + objective_domain_ = + Domain(std::numeric_limits::min(), objective_domain_.Max()); + } + return true; } diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index cb45d50797..0a53e2316d 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -24,6 +24,7 @@ #include "absl/base/attributes.h" #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" +#include "absl/log/check.h" #include "absl/strings/str_cat.h" #include "absl/types/span.h" #include "ortools/base/logging.h" @@ -55,7 +56,7 @@ class PresolveContext; // bunch of places. class SavedLiteral { public: - SavedLiteral() {} + SavedLiteral() = default; explicit SavedLiteral(int ref) : ref_(ref) {} int Get(PresolveContext* context) const; @@ -70,7 +71,7 @@ class SavedLiteral { // general affine for the linear1 involving an absolute value. class SavedVariable { public: - SavedVariable() {} + SavedVariable() = default; explicit SavedVariable(int ref) : ref_(ref) {} int Get() const; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 9bb1206b37..7d48d86dd4 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -23,7 +23,7 @@ option csharp_namespace = "Google.OrTools.Sat"; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 245 +// NEXT TAG: 249 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -1069,12 +1069,32 @@ message SatParameters { // // The test_feasibility_jump is used to only enable this for benchmarking. optional bool test_feasibility_jump = 240 [default = false]; + + // Max number of values scanned when scanning the domain of a variable. + // This number is checked to be >= 2 and <= 1e9. + optional int64 feasibility_jump_max_num_values_scanned = 245 [default = 4096]; + + // Feasibility jump focuses on the linear model, then switches to the full + // model when a linear feasible assignment is found. In the second phase, this + // parameters favors moves that keep the linear assignment feasible over move + // that will break it. + optional bool feasibility_jump_protect_linear_feasibility = 246 + [default = true]; + + // Use feasibility_jump to find improving solutions. + optional bool use_violation_ls = 244 [default = false]; + + // The following feasibility_jump parameters meant to be used by the solver on + // the different workers that use feasibility_jump. optional double feasibility_jump_decay = 242 [default = 1.0]; optional bool feasibility_jump_start_with_objective = 243 [default = false]; - - // Use a local search worker based on constraint violations to find improving - // solutions. - optional bool use_violation_ls = 244 [default = false]; + // Ratio of variables that will not have their default value upon restart. + optional double feasibility_jump_var_randomization_ratio = 247 + [default = 0.0]; + // Max distance between the default value and the pertubated value relative to + // the range of the domain of the variable. + optional double feasibility_jump_var_perburbation_range_ratio = 248 + [default = 0.2]; // Enables experimental workstealing-like shared tree search. // If non-zero, start this many complete worker threads to explore a shared diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index a6a8141c4c..cceb2f0bd4 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -35,13 +35,14 @@ #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" #include "absl/flags/flag.h" -#include "absl/status/status.h" +#include "absl/log/check.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_format.h" #include "absl/strings/string_view.h" #include "absl/synchronization/mutex.h" #include "absl/time/clock.h" #include "absl/time/time.h" +#include "absl/types/span.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/integer.h" @@ -54,7 +55,6 @@ #include "ortools/util/logging.h" #include "ortools/util/sorted_interval_list.h" #include "ortools/util/strong_integers.h" -#include "ortools/util/time_limit.h" ABSL_FLAG(bool, cp_model_dump_solutions, false, "DEBUG ONLY. If true, all the intermediate solution will be dumped " @@ -665,7 +665,7 @@ void SharedResponseManager::NewSolution( if (obj.scaling_factor() < 0) { std::swap(lb, ub); } - RegisterSolutionFound(solution_message); + RegisterSolutionFound(solution_message, num_solutions_); SOLVER_LOG(logger_, ProgressMessage(absl::StrCat(num_solutions_), wall_timer_.Get(), best, lb, ub, solution_message)); @@ -733,9 +733,12 @@ std::string ExtractSubSolverName(const std::string& improvement_info) { } void SharedResponseManager::RegisterSolutionFound( - const std::string& improvement_info) { + const std::string& improvement_info, int solution_rank) { if (improvement_info.empty()) return; - primal_improvements_count_[ExtractSubSolverName(improvement_info)]++; + const std::string subsolver_name = ExtractSubSolverName(improvement_info); + primal_improvements_count_[subsolver_name]++; + primal_improvements_min_rank_.insert({subsolver_name, solution_rank}); + primal_improvements_max_rank_[subsolver_name] = solution_rank; } void SharedResponseManager::RegisterObjectiveBoundImprovement( @@ -748,9 +751,13 @@ void SharedResponseManager::DisplayImprovementStatistics() { absl::MutexLock mutex_lock(&mutex_); if (!primal_improvements_count_.empty()) { SOLVER_LOG(logger_, ""); - SOLVER_LOG(logger_, "Solutions found per subsolver:"); + SOLVER_LOG(logger_, "Solutions found per subsolver (", num_solutions_, + "):"); for (const auto& entry : primal_improvements_count_) { - SOLVER_LOG(logger_, " '", entry.first, "': ", entry.second); + const int min_rank = primal_improvements_min_rank_[entry.first]; + const int max_rank = primal_improvements_max_rank_[entry.first]; + SOLVER_LOG(logger_, " '", entry.first, "':", entry.second, " rank:[", + min_rank, ",", max_rank, "]"); } } if (!dual_improvements_count_.empty()) { diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 34fd8ce3df..f12e5ca923 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -31,18 +31,18 @@ #include "absl/random/random.h" #include "absl/synchronization/mutex.h" #include "absl/time/time.h" -#include "ortools/base/integral_types.h" +#include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" #include "ortools/base/timer.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/integer.h" #include "ortools/sat/model.h" -#include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/util.h" #include "ortools/util/bitset.h" #include "ortools/util/logging.h" +#include "ortools/util/sorted_interval_list.h" namespace operations_research { namespace sat { @@ -358,7 +358,8 @@ class SharedResponseManager { ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); void UpdateGapIntegralInternal() ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); - void RegisterSolutionFound(const std::string& improvement_info) + void RegisterSolutionFound(const std::string& improvement_info, + int solution_rank) ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); void RegisterObjectiveBoundImprovement(const std::string& improvement_info) ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); @@ -423,6 +424,11 @@ class SharedResponseManager { // Used for statistics of the improvements found by workers. absl::btree_map primal_improvements_count_ ABSL_GUARDED_BY(mutex_); + absl::btree_map primal_improvements_min_rank_ + ABSL_GUARDED_BY(mutex_); + absl::btree_map primal_improvements_max_rank_ + ABSL_GUARDED_BY(mutex_); + absl::btree_map dual_improvements_count_ ABSL_GUARDED_BY(mutex_); diff --git a/ortools/sat/work_assignment.cc b/ortools/sat/work_assignment.cc index 4c08656066..e6adda3a80 100644 --- a/ortools/sat/work_assignment.cc +++ b/ortools/sat/work_assignment.cc @@ -187,10 +187,14 @@ absl::Span ProtoTrail::Implications(int level) const { SharedTreeManager::SharedTreeManager(Model* model) : params_(*model->GetOrCreate()), - num_workers_(params_.shared_tree_num_workers()), + num_workers_(std::max(1, params_.shared_tree_num_workers())), shared_response_manager_(model->GetOrCreate()), num_splits_wanted_(num_workers_ - 1), - max_nodes_(num_workers_ * params_.shared_tree_max_nodes_per_worker()) { + max_nodes_(params_.shared_tree_max_nodes_per_worker() >= + std::numeric_limits::max() / num_workers_ + ? std::numeric_limits::max() + : num_workers_ * + params_.shared_tree_max_nodes_per_worker()) { // Create the root node with a fake literal. nodes_.push_back( {.literal = ProtoLiteral(), @@ -597,7 +601,10 @@ LiteralIndex SharedTreeWorker::NextDecision() { CHECK(!sat_solver_->Assignment().LiteralIsTrue(decision)); return decision.Index(); } - if (objective_ == nullptr) return helper_->GetDecision(decision_policy); + if (objective_ == nullptr || + objective_->objective_var == kNoIntegerVariable) { + return helper_->GetDecision(decision_policy); + } // If the current node is close to the global lower bound, maybe try to // improve it. const IntegerValue root_obj_lb = @@ -631,6 +638,8 @@ SatSolver::Status SharedTreeWorker::Search( sat_solver_->Backtrack(0); encoder_->GetTrueLiteral(); encoder_->GetFalseLiteral(); + const bool has_objective = + objective_ != nullptr && objective_->objective_var != kNoIntegerVariable; std::vector clause; while (!time_limit_->LimitReached()) { if (!sat_solver_->FinishPropagation()) { @@ -667,7 +676,7 @@ SatSolver::Status SharedTreeWorker::Search( if (time_limit_->LimitReached()) return SatSolver::LIMIT_REACHED; if (decision.Index() == kNoLiteralIndex) { feasible_solution_observer(); - if (objective_ == nullptr) return SatSolver::FEASIBLE; + if (!has_objective) return SatSolver::FEASIBLE; const IntegerValue objective = integer_trail_->LowerBound(objective_->objective_var); sat_solver_->Backtrack(0);