diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index 65a491b4af..e05abdcc80 100644 --- a/ortools/sat/BUILD +++ b/ortools/sat/BUILD @@ -89,7 +89,7 @@ cc_library( visibility = ["//visibility:public"], deps = [ ":cp_model_cc_proto", - ":cp_model_search", + ":cp_model_loader", ":cp_model_utils", ":integer", ":model", @@ -1108,6 +1108,7 @@ cc_library( ":integer", ":linear_programming_constraint", ":model", + ":synchronization", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/synchronization", "@com_google_absl//absl/types:optional", diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc index f931bd56e8..b033038067 100644 --- a/ortools/sat/cp_model_checker.cc +++ b/ortools/sat/cp_model_checker.cc @@ -49,6 +49,7 @@ bool DomainInProtoIsValid(const ProtoWithDomain& proto) { if (proto.domain().size() % 2) return false; std::vector domain; for (int i = 0; i < proto.domain_size(); i += 2) { + if (proto.domain(i) > proto.domain(i + 1)) return false; domain.push_back({proto.domain(i), proto.domain(i + 1)}); } return IntervalsAreSortedAndNonAdjacent(domain); diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index a322e8b31a..bcd518e968 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -922,6 +922,333 @@ void ExpandLinMin(ConstraintProto* ct, PresolveContext* context) { ct->Clear(); } +// Add the implications and clauses to link one variable of a table to the +// literals controling if the tuples are possible or not. The parallel vectors +// (tuple_literals, values) contains all valid projected tuples. The +// tuples_with_any vector provides a list of tuple_literals that will support +// any value. +void ProcessOneVariable(const std::vector& tuple_literals, + const std::vector& values, int variable, + const std::vector& tuples_with_any, + PresolveContext* context) { + VLOG(2) << "Process var(" << variable << ") with domain " + << context->DomainOf(variable) << " and " << values.size() + << " active tuples, and " << tuples_with_any.size() << " any tuples"; + CHECK_EQ(tuple_literals.size(), values.size()); + std::vector> pairs; + + // Collect pairs of value-literal. + for (int i = 0; i < values.size(); ++i) { + const int64 value = values[i]; + CHECK(context->DomainContains(variable, value)); + pairs.emplace_back(value, tuple_literals[i]); + } + + // Regroup literal with the same value and add for each the clause: If all the + // tuples containing a value are false, then this value must be false too. + std::vector selected; + std::sort(pairs.begin(), pairs.end()); + for (int i = 0; i < pairs.size();) { + selected.clear(); + const int64 value = pairs[i].first; + for (; i < pairs.size() && pairs[i].first == value; ++i) { + selected.push_back(pairs[i].second); + } + + CHECK(!selected.empty() || !tuples_with_any.empty()); + if (selected.size() == 1 && tuples_with_any.empty()) { + context->InsertVarValueEncoding(selected.front(), variable, value); + } else { + const int value_literal = + context->GetOrCreateVarValueEncoding(variable, value); + BoolArgumentProto* no_support = + context->working_model->add_constraints()->mutable_bool_or(); + for (const int lit : selected) { + no_support->add_literals(lit); + context->AddImplication(lit, value_literal); + } + for (const int lit : tuples_with_any) { + no_support->add_literals(lit); + } + + // And the "value" literal. + no_support->add_literals(NegatedRef(value_literal)); + } + } +} + +// Simpler encoding for table constraints with 2 variables. +void AddSizeTwoTable( + const std::vector& vars, const std::vector>& tuples, + const std::vector>& values_per_var, + PresolveContext* context) { + CHECK_EQ(vars.size(), 2); + const int left_var = vars[0]; + const int right_var = vars[1]; + if (context->DomainOf(left_var).Size() == 1 || + context->DomainOf(right_var).Size() == 1) { + // A table constraint with at most one variable not fixed is trivially + // enforced after domain reduction. + return; + } + + std::map> left_to_right; + std::map> right_to_left; + + for (const auto& tuple : tuples) { + const int64 left_value(tuple[0]); + const int64 right_value(tuple[1]); + CHECK(context->DomainContains(left_var, left_value)); + CHECK(context->DomainContains(right_var, right_value)); + + const int left_literal = + context->GetOrCreateVarValueEncoding(left_var, left_value); + const int right_literal = + context->GetOrCreateVarValueEncoding(right_var, right_value); + left_to_right[left_literal].push_back(right_literal); + right_to_left[right_literal].push_back(left_literal); + } + + int num_implications = 0; + int num_clause_added = 0; + int num_large_clause_added = 0; + auto add_support_constraint = + [context, &num_clause_added, &num_large_clause_added, &num_implications]( + int lit, const std::vector& support_literals, + int max_support_size) { + if (support_literals.size() == max_support_size) return; + if (support_literals.size() == 1) { + context->AddImplication(lit, support_literals.front()); + num_implications++; + } else { + BoolArgumentProto* bool_or = + context->working_model->add_constraints()->mutable_bool_or(); + for (const int support_literal : support_literals) { + bool_or->add_literals(support_literal); + } + bool_or->add_literals(NegatedRef(lit)); + num_clause_added++; + if (support_literals.size() > max_support_size / 2) { + num_large_clause_added++; + } + } + }; + + for (const auto& it : left_to_right) { + add_support_constraint(it.first, it.second, values_per_var[1].size()); + } + for (const auto& it : right_to_left) { + add_support_constraint(it.first, it.second, values_per_var[0].size()); + } + VLOG(2) << "Table: 2 variables, " << tuples.size() << " tuples encoded using " + << num_clause_added << " clauses, including " + << num_large_clause_added << " large clauses, " << num_implications + << " implications"; +} + +void ExpandPositiveTable(ConstraintProto* ct, PresolveContext* context) { + const TableConstraintProto& table = ct->table(); + const std::vector vars(table.vars().begin(), table.vars().end()); + const int num_vars = table.vars_size(); + const int num_original_tuples = table.values_size() / num_vars; + + // Read tuples flat array and recreate the vector of tuples. + std::vector> tuples(num_original_tuples); + int count = 0; + for (int tuple_index = 0; tuple_index < num_original_tuples; ++tuple_index) { + for (int var_index = 0; var_index < num_vars; ++var_index) { + tuples[tuple_index].push_back(table.values(count++)); + } + } + + // Compute the set of possible values for each variable (from the table). + // Remove invalid tuples along the way. + std::vector> values_per_var(num_vars); + int new_size = 0; + for (int tuple_index = 0; tuple_index < num_original_tuples; ++tuple_index) { + bool keep = true; + for (int var_index = 0; var_index < num_vars; ++var_index) { + const int64 value = tuples[tuple_index][var_index]; + if (!context->DomainContains(vars[var_index], value)) { + keep = false; + break; + } + } + if (keep) { + for (int var_index = 0; var_index < num_vars; ++var_index) { + values_per_var[var_index].insert(tuples[tuple_index][var_index]); + } + std::swap(tuples[tuple_index], tuples[new_size]); + new_size++; + } + } + tuples.resize(new_size); + const int num_valid_tuples = tuples.size(); + + if (tuples.empty()) { + context->UpdateRuleStats("table: empty"); + return (void)context->NotifyThatModelIsUnsat(); + } + + // Update variable domains. It is redundant with presolve, but we could be + // here with presolve = false. + // Also counts the number of fixed variables. + int num_fixed_variables = 0; + for (int var_index = 0; var_index < num_vars; ++var_index) { + CHECK(context->IntersectDomainWith( + vars[var_index], + Domain::FromValues({values_per_var[var_index].begin(), + values_per_var[var_index].end()}))); + if (context->DomainOf(vars[var_index]).Size() == 1) { + num_fixed_variables++; + } + } + + if (num_fixed_variables == num_vars - 1) { + context->UpdateRuleStats("table: one variable not fixed"); + ct->Clear(); + return; + } else if (num_fixed_variables == num_vars) { + context->UpdateRuleStats("table: all variables fixed"); + ct->Clear(); + return; + } + + // Tables with two variables do not need tuple literals. + if (num_vars == 2) { + AddSizeTwoTable(vars, tuples, values_per_var, context); + context->UpdateRuleStats( + "table: expanded positive constraint with two variables"); + ct->Clear(); + return; + } + + // It is easier to compute this before compression, as compression will merge + // tuples. + int num_prefix_tuples = 0; + { + absl::flat_hash_set> prefixes; + for (const std::vector& tuple : tuples) { + prefixes.insert(absl::MakeSpan(tuple.data(), num_vars - 1)); + } + num_prefix_tuples = prefixes.size(); + } + + // TODO(user): reinvestigate ExploreSubsetOfVariablesAndAddNegatedTables. + + // Compress tuples. + const int64 any_value = kint64min; + std::vector domain_sizes; + for (int i = 0; i < num_vars; ++i) { + domain_sizes.push_back(values_per_var[i].size()); + } + CompressTuples(domain_sizes, any_value, &tuples); + const int num_compressed_tuples = tuples.size(); + + if (num_compressed_tuples == 1) { + // Domains are propagated. We can remove the constraint. + context->UpdateRuleStats("table: one tuple"); + ct->Clear(); + return; + } + + // Detect if prefix tuples are all different. + const bool prefixes_are_all_different = num_prefix_tuples == num_valid_tuples; + if (prefixes_are_all_different) { + context->UpdateRuleStats( + "TODO table: last value implied by previous values"); + } + // TODO(user): if 2 table constraints share the same valid prefix, the + // tuple literals can be reused. + // TODO(user): investigate different encoding for prefix tables. Maybe + // we can remove the need to create tuple literals. + + // Debug message to log the status of the expansion. + if (VLOG_IS_ON(2)) { + // Compute the maximum number of prefix tuples. + int64 max_num_prefix_tuples = 1; + for (int var_index = 0; var_index + 1 < num_vars; ++var_index) { + max_num_prefix_tuples = + CapProd(max_num_prefix_tuples, values_per_var[var_index].size()); + } + + std::string message = + absl::StrCat("Table: ", num_vars, + " variables, original tuples = ", num_original_tuples); + if (num_valid_tuples != num_original_tuples) { + absl::StrAppend(&message, ", valid tuples = ", num_valid_tuples); + } + if (prefixes_are_all_different) { + if (num_prefix_tuples < max_num_prefix_tuples) { + absl::StrAppend(&message, ", partial prefix = ", num_prefix_tuples, "/", + max_num_prefix_tuples); + } else { + absl::StrAppend(&message, ", full prefix = true"); + } + } else { + absl::StrAppend(&message, ", num prefix tuples = ", num_prefix_tuples); + } + if (num_compressed_tuples != num_valid_tuples) { + absl::StrAppend(&message, + ", compressed tuples = ", num_compressed_tuples); + } + VLOG(2) << message; + } + + // Log if we have only two tuples. + if (num_compressed_tuples == 2) { + context->UpdateRuleStats("TODO table: two tuples"); + } + + // Create one Boolean variable per tuple to indicate if it can still be + // selected or not. Note that we don't enforce exactly one tuple to be + // selected as this is costly. + // + // TODO(user): Investigate adding the at_most_one if the number of tuples + // is small. + // TODO(user): Investigate it we could recover a linear constraint: + // var = sum(tuple_literals[i] * values[i]) + // It could be done here or along the deductions grouping. + std::vector tuple_literals(num_compressed_tuples); + BoolArgumentProto* at_least_one_tuple = + context->working_model->add_constraints()->mutable_bool_or(); + + for (int var_index = 0; var_index < num_compressed_tuples; ++var_index) { + tuple_literals[var_index] = context->NewBoolVar(); + at_least_one_tuple->add_literals(tuple_literals[var_index]); + } + + std::vector active_tuple_literals; + std::vector active_values; + std::vector any_tuple_literals; + for (int var_index = 0; var_index < num_vars; ++var_index) { + if (values_per_var[var_index].size() == 1) continue; + + active_tuple_literals.clear(); + active_values.clear(); + any_tuple_literals.clear(); + for (int tuple_index = 0; tuple_index < tuple_literals.size(); + ++tuple_index) { + const int64 value = tuples[tuple_index][var_index]; + const int tuple_literal = tuple_literals[tuple_index]; + + if (value == any_value) { + any_tuple_literals.push_back(tuple_literal); + } else { + active_tuple_literals.push_back(tuple_literal); + active_values.push_back(value); + } + } + + if (!active_tuple_literals.empty()) { + ProcessOneVariable(active_tuple_literals, active_values, vars[var_index], + any_tuple_literals, context); + } + } + + context->UpdateRuleStats("table: expanded positive constraint"); + ct->Clear(); +} } // namespace void ExpandCpModel(PresolveOptions options, PresolveContext* context) { @@ -963,6 +1290,8 @@ void ExpandCpModel(PresolveOptions options, PresolveContext* context) { case ConstraintProto::ConstraintCase::kTable: if (ct->table().negated()) { ExpandNegativeTable(ct, context); + } else if (options.parameters.expand_table_constraints()) { + ExpandPositiveTable(ct, context); } break; default: diff --git a/ortools/sat/cp_model_loader.h b/ortools/sat/cp_model_loader.h index 9d05848d38..90bc26c7af 100644 --- a/ortools/sat/cp_model_loader.h +++ b/ortools/sat/cp_model_loader.h @@ -33,6 +33,30 @@ namespace operations_research { namespace sat { +// For an optimization problem, this contains the internal integer objective +// to minimize and information on how to display it correctly in the logs. +struct ObjectiveDefinition { + double scaling_factor = 1.0; + double offset = 0.0; + IntegerVariable objective_var = kNoIntegerVariable; + + // The objective linear expression that should be equal to objective_var. + // If not all proto variable have an IntegerVariable view, then some vars + // will be set to kNoIntegerVariable. In practice, when this is used, we make + // sure there is a view though. + std::vector vars; + std::vector coeffs; + + // List of variable that when set to their lower bound should help getting a + // better objective. This is used by some search heuristic to preferably + // assign any of the variable here to their lower bound first. + absl::flat_hash_set objective_impacting_variables; + + double ScaleIntegerObjective(IntegerValue value) const { + return (ToDouble(value) + offset) * scaling_factor; + } +}; + // Holds the mapping between CpModel proto indices and the sat::model ones. // // This also holds some information used when loading a CpModel proto. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index ded4fabddd..7cf2bed043 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -1808,6 +1808,10 @@ class FullProblemSolver : public SubSolver { shared_->rins_manager); } + if (shared->response != nullptr) { + local_model_->Register(shared->response); + } + // Level zero variable bounds sharing. if (shared_->bounds != nullptr) { RegisterVariableBoundsLevelZeroExport( diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 02076b1366..86c536f82e 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -462,6 +462,7 @@ CutGenerator CreateKnapsackCoverCutGenerator( LinearConstraint mutable_constraint; // Iterate through all knapsack constraints. + implied_bounds_processor.ClearCache(); for (const LinearConstraint& constraint : knapsack_constraints) { if (model->GetOrCreate()->LimitReached()) break; VLOG(2) << "Processing constraint: " << constraint.DebugString(); @@ -1255,6 +1256,79 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( cut, nullptr, nullptr); } +ImpliedBoundsProcessor::BestImpliedBoundInfo +ImpliedBoundsProcessor::ComputeBestImpliedBound( + IntegerVariable var, + const gtl::ITIVector& lp_values, + std::vector* implied_bound_cuts) const { + auto it = cache_.find(var); + if (it != cache_.end()) return it->second; + BestImpliedBoundInfo result; + + const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var); + for (const ImpliedBoundEntry& entry : + implied_bounds_->GetImpliedBounds(var)) { + // Only process entries with a Boolean variable currently part of the LP + // we are considering for this cut. + // + // TODO(user): the more we use cuts, the less it make sense to have a + // lot of small independent LPs. + if (!lp_vars_.contains(PositiveVariable(entry.literal_view))) { + continue; + } + + // The equation is X = lb + diff * Bool + Slack where Bool is in [0, 1] + // and slack in [0, ub - lb]. + const IntegerValue diff = entry.lower_bound - lb; + CHECK_GE(diff, 0); + const double bool_lp_value = entry.is_positive + ? lp_values[entry.literal_view] + : 1.0 - lp_values[entry.literal_view]; + const double slack_lp_value = + lp_values[var] - ToDouble(lb) - bool_lp_value * ToDouble(diff); + + // If the implied bound equation is not respected, we just add it + // to implied_bound_cuts, and skip the entry for now. + if (slack_lp_value < -1e-6) { + if (implied_bound_cuts != nullptr) { + LinearConstraint ib_cut; + std::vector> terms; + ib_cut.lb = kMinIntegerValue; // Not relevant. + ib_cut.ub = IntegerValue(0); + if (entry.is_positive) { + // X >= Indicator * (bound - lb) + lb + terms.push_back({entry.literal_view, diff}); + terms.push_back({var, IntegerValue(-1)}); + ib_cut.ub = -lb; + } else { + // X >= -Indicator * (bound - lb) + bound + terms.push_back({entry.literal_view, -diff}); + terms.push_back({var, IntegerValue(-1)}); + ib_cut.ub = -entry.lower_bound; + } + CleanTermsAndFillConstraint(&terms, &ib_cut); + implied_bound_cuts->push_back(std::move(ib_cut)); + } + continue; + } + + // We look for tight implied bounds, and amongst the tightest one, we + // prefer larger coefficient in front of the Boolean. + if (slack_lp_value + 1e-4 < result.slack_lp_value || + (slack_lp_value < result.slack_lp_value + 1e-4 && + diff > result.bound_diff)) { + result.bool_lp_value = bool_lp_value; + result.slack_lp_value = slack_lp_value; + + result.bound_diff = diff; + result.is_positive = entry.is_positive; + result.bool_var = entry.literal_view; + } + } + cache_[var] = result; + return result; +} + void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( IntegerVariable first_slack, const gtl::ITIVector& lp_values, @@ -1279,89 +1353,19 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( var = NegationOf(var); } - double best_bool_lp_value = 0.0; - double best_slack_lp_value = std::numeric_limits::infinity(); - bool best_is_positive; - IntegerValue best_diff; - IntegerVariable best_bool_var = kNoIntegerVariable; - IntegerVariable best_x = kNoIntegerVariable; - // Find the best implied bound to use. - // - // TODO(user): We could also use implied upper bound, it is why the code is - // this way so it is easy to experiment with - // for (const IntegerVariable x : {var, NegationOf(var)}) { - // instead of just trying var. - for (const IntegerVariable x : {var}) { - const IntegerValue lb = integer_trail_->LevelZeroLowerBound(x); - for (const ImpliedBoundEntry& entry : - implied_bounds_->GetImpliedBounds(x)) { - // Only process entries with a Boolean variable currently part of the LP - // we are considering for this cut. - // - // TODO(user): the more we use cuts, the less it make sense to have a - // lot of small independent LPs. - if (!lp_vars_.contains(PositiveVariable(entry.literal_view))) { - continue; - } - - // The equation is X = lb + diff * Bool + Slack where Bool is in [0, 1] - // and slack in [0, ub - lb]. - const IntegerValue diff = entry.lower_bound - lb; - CHECK_GE(diff, 0); - const double bool_lp_value = entry.is_positive - ? lp_values[entry.literal_view] - : 1.0 - lp_values[entry.literal_view]; - const double slack_lp_value = - lp_values[x] - ToDouble(lb) - bool_lp_value * ToDouble(diff); - - // If the implied bound equation is not respected, we just add it - // to implied_bound_cuts, and skip the entry for now. - if (slack_lp_value < -1e-6) { - if (implied_bound_cuts != nullptr) { - LinearConstraint ib_cut; - std::vector> terms; - ib_cut.lb = kMinIntegerValue; // Not relevant. - ib_cut.ub = cut->ub; - if (entry.is_positive) { - // X >= Indicator * (bound - lb) + lb - terms.push_back({entry.literal_view, diff}); - terms.push_back({x, IntegerValue(-1)}); - ib_cut.ub = -lb; - } else { - // X >= -Indicator * (bound - lb) + bound - terms.push_back({entry.literal_view, -diff}); - terms.push_back({x, IntegerValue(-1)}); - ib_cut.ub = -entry.lower_bound; - } - CleanTermsAndFillConstraint(&terms, &ib_cut); - implied_bound_cuts->push_back(std::move(ib_cut)); - } - continue; - } - - // We look for tight implied bounds, and amongst the tightest one, we - // prefer larger coefficient in front of the Boolean. - if (slack_lp_value + 1e-4 < best_slack_lp_value || - (slack_lp_value < best_slack_lp_value + 1e-4 && diff > best_diff)) { - best_bool_lp_value = bool_lp_value; - best_slack_lp_value = slack_lp_value; - - best_diff = diff; - best_is_positive = entry.is_positive; - best_bool_var = entry.literal_view; - - best_x = x; - } - } - } + // TODO(user): We could also use implied upper bound, that is try with + // NegationOf(var). + const BestImpliedBoundInfo info = + ComputeBestImpliedBound(var, lp_values, implied_bound_cuts); const int old_size = tmp_terms_.size(); // Shall we keep the original term ? bool keep_term = false; - if (best_x == kNoIntegerVariable) keep_term = true; - if (CapProd(std::abs(coeff.value()), best_diff.value()) == kint64max) { + if (info.bool_var == kNoIntegerVariable) keep_term = true; + if (CapProd(std::abs(coeff.value()), info.bound_diff.value()) == + kint64max) { keep_term = true; } @@ -1369,34 +1373,24 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( if (slack_infos == nullptr) { // We do not want to loose anything, so we only replace if the slack lp is // zero. - if (best_slack_lp_value > 1e-6) keep_term = true; - - // Without slack, we do need the replacement to be a lower bound, and - // this is only the case if the coefficient is positive. - if ((best_x == var) != (coeff > 0)) keep_term = true; + if (info.slack_lp_value > 1e-6) keep_term = true; } if (keep_term) { tmp_terms_.push_back({var, coeff}); } else { // Substitute. - if (best_x != var) { - coeff = -coeff; - var = NegationOf(var); - } - CHECK_EQ(best_x, var); - const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var); const IntegerValue ub = integer_trail_->LevelZeroUpperBound(var); - SlackInfo info; - info.lp_value = best_slack_lp_value; - info.lb = 0; - info.ub = ub - lb; + SlackInfo slack_info; + slack_info.lp_value = info.slack_lp_value; + slack_info.lb = 0; + slack_info.ub = ub - lb; - if (best_is_positive) { + if (info.is_positive) { // X = Indicator * diff + lb + Slack - tmp_terms_.push_back({best_bool_var, coeff * best_diff}); + tmp_terms_.push_back({info.bool_var, coeff * info.bound_diff}); if (!AddProductTo(-coeff, lb, &new_ub)) { VLOG(2) << "Overflow"; return; @@ -1405,17 +1399,17 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( tmp_terms_.push_back({first_slack, coeff}); first_slack += 2; - // slack = X - Indicator * best_diff - lb; - info.terms.push_back({var, IntegerValue(1)}); - info.terms.push_back({best_bool_var, -best_diff}); - info.offset = -lb; - slack_infos->push_back(info); + // slack = X - Indicator * info.bound_diff - lb; + slack_info.terms.push_back({var, IntegerValue(1)}); + slack_info.terms.push_back({info.bool_var, -info.bound_diff}); + slack_info.offset = -lb; + slack_infos->push_back(slack_info); } } else { // X = (1 - Indicator) * (diff) + lb + Slack // X = -Indicator * (diff) + lb + diff + Slack - tmp_terms_.push_back({best_bool_var, -coeff * best_diff}); - if (!AddProductTo(-coeff, lb + best_diff, &new_ub)) { + tmp_terms_.push_back({info.bool_var, -coeff * info.bound_diff}); + if (!AddProductTo(-coeff, lb + info.bound_diff, &new_ub)) { VLOG(2) << "Overflow"; return; } @@ -1423,11 +1417,11 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( tmp_terms_.push_back({first_slack, coeff}); first_slack += 2; - // slack = X + Indicator * best_diff - lb - diff; - info.terms.push_back({var, IntegerValue(1)}); - info.terms.push_back({best_bool_var, +best_diff}); - info.offset = -lb - best_diff; - slack_infos->push_back(info); + // slack = X + Indicator * info.bound_diff - lb - diff; + slack_info.terms.push_back({var, IntegerValue(1)}); + slack_info.terms.push_back({info.bool_var, +info.bound_diff}); + slack_info.offset = -lb - info.bound_diff; + slack_infos->push_back(slack_info); } } changed = true; diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 652e7706f1..921dbb90fa 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -100,8 +100,25 @@ class ImpliedBoundsProcessor { // Add a new variable that could be used in the new cuts. void AddLpVariable(IntegerVariable var) { lp_vars_.insert(var); } + // Must be called before we process any constraints with a different + // lp_values or level zero bounds. + void ClearCache() const { cache_.clear(); } + private: + struct BestImpliedBoundInfo { + double bool_lp_value = 0.0; + double slack_lp_value = std::numeric_limits::infinity(); + bool is_positive; + IntegerValue bound_diff; + IntegerVariable bool_var = kNoIntegerVariable; + }; + BestImpliedBoundInfo ComputeBestImpliedBound( + IntegerVariable var, + const gtl::ITIVector& lp_values, + std::vector* implied_bound_cuts) const; + absl::flat_hash_set lp_vars_; + mutable absl::flat_hash_map cache_; // Data from the constructor. IntegerTrail* integer_trail_; diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index c39bc8815f..c64d10e118 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -17,6 +17,7 @@ #include #include +#include "ortools/sat/cp_model_loader.h" #include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_programming_constraint.h" @@ -778,9 +779,6 @@ SatSolver::Status SolveIntegerProblem(Model* model) { return sat_solver->UnsatStatus(); } if (model->Get() != nullptr) { - // If RINS is activated, we need to make sure the SolutionDetails is - // created. - model->GetOrCreate(); num_decisions_without_rins++; // TODO(user): Experiment more around dynamically changing the // threshold for trigerring RINS. Alternatively expose this as parameter diff --git a/ortools/sat/integer_search.h b/ortools/sat/integer_search.h index 657b51f91a..7bf3165b89 100644 --- a/ortools/sat/integer_search.h +++ b/ortools/sat/integer_search.h @@ -59,30 +59,6 @@ struct SearchHeuristics { // given model. void ConfigureSearchHeuristics(Model* model); -// For an optimization problem, this contains the internal integer objective -// to minimize and information on how to display it correctly in the logs. -struct ObjectiveDefinition { - double scaling_factor = 1.0; - double offset = 0.0; - IntegerVariable objective_var = kNoIntegerVariable; - - // The objective linear expression that should be equal to objective_var. - // If not all proto variable have an IntegerVariable view, then some vars - // will be set to kNoIntegerVariable. In practice, when this is used, we make - // sure there is a view though. - std::vector vars; - std::vector coeffs; - - // List of variable that when set to their lower bound should help getting a - // better objective. This is used by some search heuristic to preferably - // assign any of the variable here to their lower bound first. - absl::flat_hash_set objective_impacting_variables; - - double ScaleIntegerObjective(IntegerValue value) const { - return (ToDouble(value) + offset) * scaling_factor; - } -}; - // Callbacks that will be called when the search goes back to level 0. // Callbacks should return false if the propagation fails. struct LevelZeroCallbackHelper { diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 921eb793fd..9bfc116907 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -1211,6 +1211,7 @@ bool LinearProgrammingConstraint::Propagate() { // // TODO(user): Refactor so that they are just normal cut generators? if (trail_->CurrentDecisionLevel() == 0) { + implied_bounds_processor_.ClearCache(); if (sat_parameters_.add_mir_cuts()) AddMirCuts(); if (sat_parameters_.add_cg_cuts()) AddCGCuts(); } diff --git a/ortools/sat/rins.cc b/ortools/sat/rins.cc index 63d972a723..4409a3d6ad 100644 --- a/ortools/sat/rins.cc +++ b/ortools/sat/rins.cc @@ -16,6 +16,7 @@ #include "ortools/sat/cp_model_loader.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_programming_constraint.h" +#include "ortools/sat/synchronization.h" namespace operations_research { namespace sat { @@ -67,7 +68,7 @@ SharedRINSNeighborhoodManager::GetUnexploredNeighborhood() { void AddRINSNeighborhood(Model* model) { IntegerTrail* const integer_trail = model->GetOrCreate(); - auto* solution_details = model->Get(); + auto* response_manager = model->Get(); const RINSVariables& rins_vars = *model->GetOrCreate(); RINSNeighborhood rins_neighborhood; @@ -75,6 +76,7 @@ void AddRINSNeighborhood(Model* model) { const double tolerance = 1e-6; for (const RINSVariable& rins_var : rins_vars.vars) { const IntegerVariable positive_var = rins_var.positive_var; + const int model_var = rins_var.model_var; if (integer_trail->IsCurrentlyIgnored(positive_var)) continue; @@ -82,7 +84,8 @@ void AddRINSNeighborhood(Model* model) { if (lp == nullptr || !lp->HasSolution()) continue; const double lp_value = lp->GetSolutionValue(positive_var); - if (solution_details == nullptr || solution_details->solution_count == 0) { + if (response_manager == nullptr || + response_manager->SolutionsRepository().NumSolutions() == 0) { // The tolerance make sure that if the lp_values is close to an integer, // then we fix the variable to this integer value. // @@ -109,10 +112,11 @@ void AddRINSNeighborhood(Model* model) { << " UB: " << integer_trail->UpperBound(positive_var); } } else { - if (positive_var >= solution_details->best_solution.size()) continue; + const SharedSolutionRepository::Solution& solution = + response_manager->SolutionsRepository().GetSolution(0); const IntegerValue best_solution_value = - solution_details->best_solution[positive_var]; + IntegerValue(solution.variable_values[model_var]); if (std::abs(best_solution_value.value() - lp_value) < 1e-4) { if (best_solution_value >= integer_trail->LowerBound(positive_var) && best_solution_value <= integer_trail->UpperBound(positive_var)) { diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index b9f554cc5b..3860e8409a 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -21,7 +21,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 158 +// NEXT TAG: 159 message SatParameters { // ========================================================================== // Branching and polarity @@ -365,8 +365,8 @@ message SatParameters { // ========================================================================== // During presolve, only try to perform the bounded variable elimination (BVE) - // of a variable x if the number of occurences of x times the number of - // occurences of not(x) is not greater than this parameter. + // of a variable x if the number of occurrences of x times the number of + // occurrences of not(x) is not greater than this parameter. optional int32 presolve_bve_threshold = 54 [default = 500]; // During presolve, we apply BVE only if this weight times the number of @@ -416,6 +416,10 @@ message SatParameters { // If true, the automaton constraints are expanded. optional bool expand_automaton_constraints = 143 [default = true]; + // If true, the positive table constraints are expanded. + // Note that currently, negative table constraints are always expanded. + optional bool expand_table_constraints = 158 [default = true]; + // During presolve, we use a maximum clique heuristic to merge together // no-overlap constraints or at most one constraints. This code can be slow, // so we have a limit in place on the number of explored nodes in the diff --git a/ortools/sat/simplification.cc b/ortools/sat/simplification.cc index f18535b1f3..171834a909 100644 --- a/ortools/sat/simplification.cc +++ b/ortools/sat/simplification.cc @@ -301,7 +301,7 @@ bool SatPresolver::ProcessAllClauses() { clause_to_process_.pop_front(); if (!ProcessClauseToSimplifyOthers(ci)) return false; if (++num_skipped_checks >= kCheckFrequency) { - if (num_inspected_signatures_ > 1e9) { + if (num_inspected_signatures_ + num_inspected_literals_ > 1e9) { VLOG(1) << "Aborting ProcessAllClauses() because work limit has been " "reached"; return true; @@ -332,6 +332,7 @@ bool SatPresolver::Presolve(const std::vector& can_be_removed) { DisplayStats(timer.Get()); if (time_limit_ != nullptr && time_limit_->LimitReached()) return true; + if (num_inspected_signatures_ + num_inspected_literals_ > 1e9) return true; InitializePriorityQueue(); while (var_pq_.Size() > 0) { @@ -341,6 +342,8 @@ bool SatPresolver::Presolve(const std::vector& can_be_removed) { if (CrossProduct(Literal(var, true))) { if (!ProcessAllClauses()) return false; } + if (time_limit_ != nullptr && time_limit_->LimitReached()) return true; + if (num_inspected_signatures_ + num_inspected_literals_ > 1e9) return true; } DisplayStats(timer.Get()); @@ -557,7 +560,8 @@ bool SatPresolver::ProcessClauseToSimplifyOthersUsingLiteral( // 'a' are a subset of the one in 'b'. We use the signatures to abort // early as a speed optimization. if (ci != clause_index && (clause_signature & ~ci_signature) == 0 && - SimplifyClause(clause, &clauses_[ci], &opposite_literal)) { + SimplifyClause(clause, &clauses_[ci], &opposite_literal, + &num_inspected_literals_)) { if (opposite_literal == kNoLiteralIndex) { need_cleaning = true; Remove(ci); @@ -641,7 +645,8 @@ bool SatPresolver::ProcessClauseToSimplifyOthers(ClauseIndex clause_index) { // applies to the second call to // ProcessClauseToSimplifyOthersUsingLiteral() above too. if ((clause_signature & ~ci_signature) == 0 && - SimplifyClause(clause, &clauses_[ci], &opposite_literal)) { + SimplifyClause(clause, &clauses_[ci], &opposite_literal, + &num_inspected_literals_)) { DCHECK_EQ(opposite_literal, lit.NegatedIndex()); if (clauses_[ci].empty()) return false; // UNSAT. if (drat_proof_handler_ != nullptr) { @@ -910,10 +915,15 @@ void SatPresolver::DisplayStats(double elapsed_seconds) { } bool SimplifyClause(const std::vector& a, std::vector* b, - LiteralIndex* opposite_literal) { + LiteralIndex* opposite_literal, + int64* num_inspected_literals) { if (b->size() < a.size()) return false; DCHECK(std::is_sorted(a.begin(), a.end())); DCHECK(std::is_sorted(b->begin(), b->end())); + if (num_inspected_literals != nullptr) { + *num_inspected_literals += a.size(); + *num_inspected_literals += b->size(); + } *opposite_literal = LiteralIndex(-1); diff --git a/ortools/sat/simplification.h b/ortools/sat/simplification.h index 1cad2279bd..b874c06f39 100644 --- a/ortools/sat/simplification.h +++ b/ortools/sat/simplification.h @@ -331,6 +331,7 @@ class SatPresolver { // The cached value of ComputeSignatureOfClauseVariables() for each clause. std::vector signatures_; // Indexed by ClauseIndex int64 num_inspected_signatures_ = 0; + int64 num_inspected_literals_ = 0; // Occurrence list. For each literal, contains the ClauseIndex of the clause // that contains it (ordered by clause index). @@ -361,8 +362,15 @@ class SatPresolver { // the clause a with one of its literal negated is a subset of b, in which // case opposite_literal is set to this negated literal index. Moreover, this // opposite_literal is then removed from b. +// +// If num_inspected_literals_ is not nullptr, the "complexity" of this function +// will be added to it in order to track the amount of work done. +// +// TODO(user): when a.size() << b.size(), we should use binary search instead +// of scanning b linearly. bool SimplifyClause(const std::vector& a, std::vector* b, - LiteralIndex* opposite_literal); + LiteralIndex* opposite_literal, + int64* num_inspected_literals = nullptr); // Visible for testing. Returns kNoLiteralIndex except if: // - a and b differ in only one literal. diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index d6588624fd..b2eb6163dc 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -23,7 +23,6 @@ #include "ortools/base/integral_types.h" #include "ortools/base/stl_util.h" #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_search.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/integer.h" #include "ortools/sat/model.h"