diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index c6d0fa05e6..bb15ba6a75 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_loader", + ":cp_model_mapping", ":cp_model_utils", ":integer", ":linear_programming_constraint", @@ -129,6 +129,7 @@ cc_library( hdrs = ["cp_model_search.h"], deps = [ ":cp_model_cc_proto", + ":cp_model_mapping", ":cp_model_utils", ":integer", ":integer_search", @@ -214,6 +215,29 @@ cc_library( ], ) +cc_library( + name = "cp_model_mapping", + hdrs = ["cp_model_mapping.h"], + visibility = ["//visibility:public"], + deps = [ + ":cp_model_cc_proto", + ":cp_model_utils", + ":integer", + ":intervals", + ":linear_constraint", + ":model", + ":sat_base", + "@com_google_absl//absl/base:core_headers", + "@com_google_absl//absl/container:flat_hash_map", + "@com_google_absl//absl/container:flat_hash_set", + "//ortools/base", + "//ortools/base:int_type", + "//ortools/base:strong_vector", + "//ortools/util:logging", + "//ortools/util:sorted_interval_list", + ], +) + cc_library( name = "cp_model_loader", srcs = ["cp_model_loader.cc"], @@ -224,6 +248,7 @@ cc_library( ":circuit", ":cp_constraints", ":cp_model_cc_proto", + ":cp_model_mapping", ":cp_model_symmetries", ":cp_model_utils", ":cumulative", @@ -233,6 +258,7 @@ cc_library( ":integer", ":integer_expr", ":intervals", + ":linear_relaxation", ":model", ":pb_constraint", ":precedences", @@ -304,9 +330,11 @@ cc_library( srcs = ["cp_model_presolve.cc"], hdrs = ["cp_model_presolve.h"], deps = [ + ":circuit", ":cp_model_cc_proto", ":cp_model_checker", ":cp_model_expand", + ":cp_model_mapping", ":cp_model_loader", ":cp_model_objective", ":cp_model_symmetries", @@ -635,6 +663,7 @@ cc_library( srcs = ["integer_search.cc"], hdrs = ["integer_search.h"], deps = [ + ":cp_model_mapping", ":integer", ":linear_programming_constraint", ":probing", @@ -653,7 +682,7 @@ cc_library( srcs = ["lb_tree_search.cc"], hdrs = ["lb_tree_search.h"], deps = [ - ":cp_model_loader", + ":cp_model_mapping", ":integer", ":integer_search", ":sat_base", @@ -921,7 +950,8 @@ cc_library( srcs = ["linear_relaxation.cc"], hdrs = ["linear_relaxation.h"], deps = [ - ":cp_model_loader", + ":circuit", + ":cp_model_mapping", ":integer", ":linear_programming_constraint", ":model", @@ -1056,6 +1086,7 @@ cc_library( deps = [ ":boolean_problem", ":boolean_problem_cc_proto", + ":cp_model_mapping", ":encoding", ":integer", ":integer_expr", @@ -1225,7 +1256,7 @@ cc_library( visibility = ["//visibility:public"], deps = [ ":cp_model_cc_proto", - ":cp_model_loader", + ":cp_model_mapping", ":cp_model_utils", ":linear_programming_constraint", ":model", @@ -1249,7 +1280,7 @@ cc_library( visibility = ["//visibility:public"], deps = [ ":cp_model_cc_proto", - ":cp_model_loader", + ":cp_model_mapping", ":integer", ":linear_constraint", ":synchronization", @@ -1269,7 +1300,7 @@ cc_library( hdrs = ["rins.h"], visibility = ["//visibility:public"], deps = [ - ":cp_model_loader", + ":cp_model_mapping", ":integer", ":linear_programming_constraint", ":model", diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 605fca5419..5c53649506 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -21,7 +21,7 @@ #include "absl/container/flat_hash_set.h" #include "absl/synchronization/mutex.h" #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_programming_constraint.h" @@ -256,11 +256,35 @@ std::vector NeighborhoodGeneratorHelper::GetActiveIntervals( // We filter out fixed intervals. Because of presolve, if there is an // enforcement literal, it cannot be fixed. - if (interval_ct.enforcement_literal().empty() && - IsConstant(PositiveRef(interval_ct.interval().start())) && - IsConstant(PositiveRef(interval_ct.interval().size())) && - IsConstant(PositiveRef(interval_ct.interval().end()))) { - continue; + if (interval_ct.enforcement_literal().empty()) { + if (interval_ct.interval().has_start_view()) { + bool is_constant = true; + for (const int v : interval_ct.interval().start_view().vars()) { + if (!IsConstant(v)) { + is_constant = false; + break; + } + } + for (const int v : interval_ct.interval().size_view().vars()) { + if (!IsConstant(v)) { + is_constant = false; + break; + } + } + for (const int v : interval_ct.interval().end_view().vars()) { + if (!IsConstant(v)) { + is_constant = false; + break; + } + } + if (is_constant) continue; + } else { + if (IsConstant(PositiveRef(interval_ct.interval().start())) && + IsConstant(PositiveRef(interval_ct.interval().size())) && + IsConstant(PositiveRef(interval_ct.interval().end()))) { + continue; + } + } } active_intervals.push_back(i); @@ -645,6 +669,51 @@ Neighborhood ConstraintGraphNeighborhoodGenerator::Generate( } namespace { + +LinearExpressionProto GetStart(const IntervalConstraintProto& interval) { + if (interval.has_start_view()) return interval.start_view(); + LinearExpressionProto result; + result.add_vars(interval.start()); + result.add_coeffs(1); + return result; +} + +LinearExpressionProto GetSize(const IntervalConstraintProto& interval) { + if (interval.has_size_view()) return interval.size_view(); + LinearExpressionProto result; + result.add_vars(interval.size()); + result.add_coeffs(1); + return result; +} + +LinearExpressionProto GetEnd(const IntervalConstraintProto& interval) { + if (interval.has_end_view()) return interval.end_view(); + LinearExpressionProto result; + result.add_vars(interval.end()); + result.add_coeffs(1); + return result; +} + +int64_t GetLinearExpressionValue(const LinearExpressionProto& expr, + const CpSolverResponse& initial_solution) { + int64_t result = expr.offset(); + for (int i = 0; i < expr.vars_size(); ++i) { + result += expr.coeffs(i) * initial_solution.solution(expr.vars(i)); + } + return result; +} + +void AddLinearExpressionToConstraint(const int64_t coeff, + const LinearExpressionProto& expr, + LinearConstraintProto* constraint, + int64_t* rhs_offset) { + *rhs_offset -= coeff * expr.offset(); + for (int i = 0; i < expr.vars_size(); ++i) { + constraint->add_vars(expr.vars(i)); + constraint->add_coeffs(expr.coeffs(i) * coeff); + } +} + void AddPrecedenceConstraints(const absl::Span intervals, const absl::flat_hash_set& ignored_intervals, const CpSolverResponse& initial_solution, @@ -658,50 +727,63 @@ void AddPrecedenceConstraints(const absl::Span intervals, const ConstraintProto& interval_ct = helper.ModelProto().constraints(i); // TODO(user): we ignore size zero for now. - const int size_var = interval_ct.interval().size(); - if (initial_solution.solution(size_var) == 0) continue; + const LinearExpressionProto size_var = GetSize(interval_ct.interval()); + if (GetLinearExpressionValue(size_var, initial_solution) == 0) continue; + + const LinearExpressionProto start_var = GetStart(interval_ct.interval()); + const int64_t start_value = + GetLinearExpressionValue(start_var, initial_solution); - const int start_var = interval_ct.interval().start(); - const int64_t start_value = initial_solution.solution(start_var); start_interval_pairs.push_back({start_value, i}); } std::sort(start_interval_pairs.begin(), start_interval_pairs.end()); // Add precedence between the remaining intervals, forcing their order. for (int i = 0; i + 1 < start_interval_pairs.size(); ++i) { - const int before_start = helper.ModelProto() - .constraints(start_interval_pairs[i].second) - .interval() - .start(); - const int before_end = helper.ModelProto() - .constraints(start_interval_pairs[i].second) - .interval() - .end(); - const int after_start = helper.ModelProto() - .constraints(start_interval_pairs[i + 1].second) - .interval() - .start(); + const LinearExpressionProto before_start = + GetStart(helper.ModelProto() + .constraints(start_interval_pairs[i].second) + .interval()); + const LinearExpressionProto before_end = + GetEnd(helper.ModelProto() + .constraints(start_interval_pairs[i].second) + .interval()); + const LinearExpressionProto after_start = + GetStart(helper.ModelProto() + .constraints(start_interval_pairs[i + 1].second) + .interval()); // If the end was smaller we keep it that way, otherwise we just order the // start variables. LinearConstraintProto* linear = neighborhood->delta.add_constraints()->mutable_linear(); linear->add_domain(std::numeric_limits::min()); - if (initial_solution.solution(before_end) <= - initial_solution.solution(after_start)) { + int64_t rhs_offset = 0; + if (GetLinearExpressionValue(before_end, initial_solution) <= + GetLinearExpressionValue(after_start, initial_solution)) { // If the end was smaller than the next start, keep it that way. - linear->add_domain(0); - linear->add_vars(before_end); + AddLinearExpressionToConstraint(1, before_end, linear, &rhs_offset); } else { // Otherwise, keep the same minimum separation. This is done in order // to "simplify" the neighborhood. - linear->add_domain(initial_solution.solution(before_start) - - initial_solution.solution(after_start)); - linear->add_vars(before_start); + rhs_offset = GetLinearExpressionValue(before_start, initial_solution) - + GetLinearExpressionValue(after_start, initial_solution); + AddLinearExpressionToConstraint(1, before_start, linear, &rhs_offset); + } + + AddLinearExpressionToConstraint(-1, after_start, linear, &rhs_offset); + linear->add_domain(rhs_offset); + + // The linear constraint should be satisfied by the current solution. + if (DEBUG_MODE) { + int64_t activity = 0; + for (int i = 0; i < linear->vars().size(); ++i) { + activity += + linear->coeffs(i) * initial_solution.solution(linear->vars(i)); + } + CHECK_GE(activity, linear->domain(0)); + CHECK_LE(activity, linear->domain(1)); } - linear->add_coeffs(1); - linear->add_vars(after_start); - linear->add_coeffs(-1); } } } // namespace @@ -797,8 +879,9 @@ Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate( for (const int i : active_intervals) { const ConstraintProto& interval_ct = helper_.ModelProto().constraints(i); - const int start_var = interval_ct.interval().start(); - const int64_t start_value = initial_solution.solution(start_var); + const LinearExpressionProto start_var = GetStart(interval_ct.interval()); + const int64_t start_value = + GetLinearExpressionValue(start_var, initial_solution); start_interval_pairs.push_back({start_value, i}); } std::sort(start_interval_pairs.begin(), start_interval_pairs.end()); diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 7256aa2459..5d984866f7 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -115,9 +115,9 @@ bool ConstraintIsNEq(const LinearConstraintProto& proto, } // namespace -void CpModelMapping::CreateVariables(const CpModelProto& model_proto, - bool view_all_booleans_as_integers, - Model* m) { +void LoadVariables(const CpModelProto& model_proto, + bool view_all_booleans_as_integers, Model* m) { + auto* mapping = m->GetOrCreate(); const int num_proto_variables = model_proto.variables_size(); // All [0, 1] variables always have a corresponding Boolean, even if it is @@ -130,14 +130,14 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto, std::vector false_variables; std::vector true_variables; - booleans_.resize(num_proto_variables, kNoBooleanVariable); - reverse_boolean_map_.resize(num_proto_variables, -1); + mapping->booleans_.resize(num_proto_variables, kNoBooleanVariable); + mapping->reverse_boolean_map_.resize(num_proto_variables, -1); for (int i = 0; i < num_proto_variables; ++i) { const auto& domain = model_proto.variables(i).domain(); if (domain.size() != 2) continue; if (domain[0] >= 0 && domain[1] <= 1) { - booleans_[i] = new_var; - reverse_boolean_map_[new_var] = i; + mapping->booleans_[i] = new_var; + mapping->reverse_boolean_map_[new_var] = i; if (domain[1] == 0) { false_variables.push_back(new_var); } else if (domain[0] == 1) { @@ -188,7 +188,7 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto, // Make sure any unused variable, that is not already a Boolean is // considered "used". for (int i = 0; i < num_proto_variables; ++i) { - if (booleans_[i] == kNoBooleanVariable) { + if (mapping->booleans_[i] == kNoBooleanVariable) { used_variables.insert(i); } } @@ -198,18 +198,19 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto, used_variables.end()); gtl::STLSortAndRemoveDuplicates(&var_to_instantiate_as_integer); } - integers_.resize(num_proto_variables, kNoIntegerVariable); + mapping->integers_.resize(num_proto_variables, kNoIntegerVariable); auto* integer_trail = m->GetOrCreate(); integer_trail->ReserveSpaceForNumVariables( var_to_instantiate_as_integer.size()); - reverse_integer_map_.resize(2 * var_to_instantiate_as_integer.size(), -1); + mapping->reverse_integer_map_.resize(2 * var_to_instantiate_as_integer.size(), + -1); for (const int i : var_to_instantiate_as_integer) { const auto& var_proto = model_proto.variables(i); - integers_[i] = + mapping->integers_[i] = integer_trail->AddIntegerVariable(ReadDomainFromProto(var_proto)); - DCHECK_LT(integers_[i], reverse_integer_map_.size()); - reverse_integer_map_[integers_[i]] = i; + DCHECK_LT(mapping->integers_[i], mapping->reverse_integer_map_.size()); + mapping->reverse_integer_map_[mapping->integers_[i]] = i; } auto* encoder = m->GetOrCreate(); @@ -217,16 +218,18 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto, // Link any variable that has both views. for (int i = 0; i < num_proto_variables; ++i) { - if (integers_[i] == kNoIntegerVariable) continue; - if (booleans_[i] == kNoBooleanVariable) continue; + if (mapping->integers_[i] == kNoIntegerVariable) continue; + if (mapping->booleans_[i] == kNoBooleanVariable) continue; // Associate with corresponding integer variable. - encoder->AssociateToIntegerEqualValue(sat::Literal(booleans_[i], true), - integers_[i], IntegerValue(1)); + encoder->AssociateToIntegerEqualValue( + sat::Literal(mapping->booleans_[i], true), mapping->integers_[i], + IntegerValue(1)); } // Create the interval variables. - intervals_.resize(model_proto.constraints_size(), kNoIntervalVariable); + mapping->intervals_.resize(model_proto.constraints_size(), + kNoIntervalVariable); for (int c = 0; c < model_proto.constraints_size(); ++c) { const ConstraintProto& ct = model_proto.constraints(c); if (ct.constraint_case() != ConstraintProto::ConstraintCase::kInterval) { @@ -234,41 +237,43 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto, } if (HasEnforcementLiteral(ct)) { const sat::Literal enforcement_literal = - Literal(ct.enforcement_literal(0)); + mapping->Literal(ct.enforcement_literal(0)); // TODO(user): Fix the constant variable situation. An optional interval // with constant start/end or size cannot share the same constant // variable if it is used in non-optional situation. if (ct.interval().has_start_view()) { - intervals_[c] = intervals_repository->CreateInterval( - LoadAffineView(ct.interval().start_view()), - LoadAffineView(ct.interval().end_view()), - LoadAffineView(ct.interval().size_view()), + mapping->intervals_[c] = intervals_repository->CreateInterval( + mapping->LoadAffineView(ct.interval().start_view()), + mapping->LoadAffineView(ct.interval().end_view()), + mapping->LoadAffineView(ct.interval().size_view()), enforcement_literal.Index(), /*add_linear_relation=*/false); } else { - intervals_[c] = m->Add(NewOptionalInterval( - Integer(ct.interval().start()), Integer(ct.interval().end()), - Integer(ct.interval().size()), enforcement_literal)); + mapping->intervals_[c] = m->Add(NewOptionalInterval( + mapping->Integer(ct.interval().start()), + mapping->Integer(ct.interval().end()), + mapping->Integer(ct.interval().size()), enforcement_literal)); } } else { if (ct.interval().has_start_view()) { - intervals_[c] = intervals_repository->CreateInterval( - LoadAffineView(ct.interval().start_view()), - LoadAffineView(ct.interval().end_view()), - LoadAffineView(ct.interval().size_view()), kNoLiteralIndex, + mapping->intervals_[c] = intervals_repository->CreateInterval( + mapping->LoadAffineView(ct.interval().start_view()), + mapping->LoadAffineView(ct.interval().end_view()), + mapping->LoadAffineView(ct.interval().size_view()), kNoLiteralIndex, /*add_linear_relation=*/false); } else { - intervals_[c] = m->Add(NewInterval(Integer(ct.interval().start()), - Integer(ct.interval().end()), - Integer(ct.interval().size()))); + mapping->intervals_[c] = + m->Add(NewInterval(mapping->Integer(ct.interval().start()), + mapping->Integer(ct.interval().end()), + mapping->Integer(ct.interval().size()))); } } - already_loaded_ct_.insert(&ct); + mapping->already_loaded_ct_.insert(&ct); } } -void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto, - Model* m) { +void LoadBooleanSymmetries(const CpModelProto& model_proto, Model* m) { + auto* mapping = m->GetOrCreate(); const SymmetryProto symmetry = model_proto.symmetry(); if (symmetry.permutations().empty()) return; @@ -280,7 +285,7 @@ void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto, for (const SparsePermutationProto& perm : symmetry.permutations()) { bool all_bool = true; for (const int var : perm.support()) { - if (!IsBoolean(var)) { + if (!mapping->IsBoolean(var)) { all_bool = false; break; } @@ -297,7 +302,8 @@ void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto, const int saved_support_index = support_index; for (int j = 0; j < size; ++j) { const int var = perm.support(support_index++); - literal_permutation->AddToCurrentCycle(Literal(var).Index().value()); + literal_permutation->AddToCurrentCycle( + mapping->Literal(var).Index().value()); } literal_permutation->CloseCurrentCycle(); @@ -307,7 +313,7 @@ void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto, for (int j = 0; j < size; ++j) { const int var = perm.support(support_index++); literal_permutation->AddToCurrentCycle( - Literal(var).NegatedIndex().value()); + mapping->Literal(var).NegatedIndex().value()); } literal_permutation->CloseCurrentCycle(); } @@ -325,8 +331,8 @@ void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto, // // TODO(user): Regroup/presolve two encoding like b => x > 2 and the same // Boolean b => x > 5. These shouldn't happen if we merge linear constraints. -void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, - Model* m) { +void ExtractEncoding(const CpModelProto& model_proto, Model* m) { + auto* mapping = m->GetOrCreate(); auto* encoder = m->GetOrCreate(); auto* integer_trail = m->GetOrCreate(); auto* sat_solver = m->GetOrCreate(); @@ -386,7 +392,8 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, if (ct.linear().vars_size() != 1) continue; // ct is a linear constraint with one term and one enforcement literal. - const sat::Literal enforcement_literal = Literal(ct.enforcement_literal(0)); + const sat::Literal enforcement_literal = + mapping->Literal(ct.enforcement_literal(0)); const int ref = ct.linear().vars(0); const int var = PositiveRef(ref); @@ -400,17 +407,17 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, if (domain_if_enforced.NumIntervals() == 1) { if (domain_if_enforced.Max() >= domain.Max() && domain_if_enforced.Min() > domain.Min()) { - inequalities.push_back( - {&ct, enforcement_literal, - IntegerLiteral::GreaterOrEqual( - Integer(var), IntegerValue(domain_if_enforced.Min()))}); + inequalities.push_back({&ct, enforcement_literal, + IntegerLiteral::GreaterOrEqual( + mapping->Integer(var), + IntegerValue(domain_if_enforced.Min()))}); implied_bounds->Add(enforcement_literal, inequalities.back().i_lit); } else if (domain_if_enforced.Min() <= domain.Min() && domain_if_enforced.Max() < domain.Max()) { - inequalities.push_back( - {&ct, enforcement_literal, - IntegerLiteral::LowerOrEqual( - Integer(var), IntegerValue(domain_if_enforced.Max()))}); + inequalities.push_back({&ct, enforcement_literal, + IntegerLiteral::LowerOrEqual( + mapping->Integer(var), + IntegerValue(domain_if_enforced.Max()))}); implied_bounds->Add(enforcement_literal, inequalities.back().i_lit); } } @@ -426,7 +433,7 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, var_to_equalities[var].push_back( {&ct, enforcement_literal, inter.Min(), true}); if (domain.Contains(inter.Min())) { - variables_to_encoded_values_[var].insert(inter.Min()); + mapping->variables_to_encoded_values_[var].insert(inter.Min()); } } } @@ -437,7 +444,7 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, var_to_equalities[var].push_back( {&ct, enforcement_literal, inter.Min(), false}); if (domain.Contains(inter.Min())) { - variables_to_encoded_values_[var].insert(inter.Min()); + mapping->variables_to_encoded_values_[var].insert(inter.Min()); } } } @@ -469,23 +476,23 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, ++num_inequalities; encoder->AssociateToIntegerLiteral(inequalities[i].literal, inequalities[i].i_lit); - already_loaded_ct_.insert(inequalities[i].ct); - already_loaded_ct_.insert(inequalities[i + 1].ct); + mapping->already_loaded_ct_.insert(inequalities[i].ct); + mapping->already_loaded_ct_.insert(inequalities[i + 1].ct); } } // Encode the half-inequalities. int num_half_inequalities = 0; for (const auto inequality : inequalities) { - if (ConstraintIsAlreadyLoaded(inequality.ct)) continue; + if (mapping->ConstraintIsAlreadyLoaded(inequality.ct)) continue; m->Add( Implication(inequality.literal, encoder->GetOrCreateAssociatedLiteral(inequality.i_lit))); if (sat_solver->IsModelUnsat()) return; ++num_half_inequalities; - already_loaded_ct_.insert(inequality.ct); - is_half_encoding_ct_.insert(inequality.ct); + mapping->already_loaded_ct_.insert(inequality.ct); + mapping->is_half_encoding_ct_.insert(inequality.ct); } if (!inequalities.empty()) { @@ -517,10 +524,11 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, } ++num_equalities; - encoder->AssociateToIntegerEqualValue(encoding[j].literal, integers_[i], + encoder->AssociateToIntegerEqualValue(encoding[j].literal, + mapping->integers_[i], IntegerValue(encoding[j].value)); - already_loaded_ct_.insert(encoding[j].ct); - already_loaded_ct_.insert(encoding[j + 1].ct); + mapping->already_loaded_ct_.insert(encoding[j].ct); + mapping->already_loaded_ct_.insert(encoding[j + 1].ct); values.insert(encoding[j].value); } @@ -536,9 +544,9 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, // however, that in the presolve, we should only use the "representative" in // linear constraints, so we should be fine. for (const auto equality : encoding) { - if (ConstraintIsAlreadyLoaded(equality.ct)) continue; + if (mapping->ConstraintIsAlreadyLoaded(equality.ct)) continue; const class Literal eq = encoder->GetOrCreateLiteralAssociatedToEquality( - integers_[i], IntegerValue(equality.value)); + mapping->integers_[i], IntegerValue(equality.value)); if (equality.is_equality) { m->Add(Implication(equality.literal, eq)); } else { @@ -546,13 +554,13 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, } ++num_half_equalities; - already_loaded_ct_.insert(equality.ct); - is_half_encoding_ct_.insert(equality.ct); + mapping->already_loaded_ct_.insert(equality.ct); + mapping->is_half_encoding_ct_.insert(equality.ct); } // Update stats. if (VLOG_IS_ON(1)) { - if (encoder->VariableIsFullyEncoded(integers_[i])) { + if (encoder->VariableIsFullyEncoded(mapping->integers_[i])) { ++num_fully_encoded; } else { ++num_partially_encoded; @@ -572,8 +580,9 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto, } } -void CpModelMapping::PropagateEncodingFromEquivalenceRelations( - const CpModelProto& model_proto, Model* m) { +void PropagateEncodingFromEquivalenceRelations(const CpModelProto& model_proto, + Model* m) { + auto* mapping = m->GetOrCreate(); auto* encoder = m->GetOrCreate(); auto* sat_solver = m->GetOrCreate(); @@ -589,8 +598,8 @@ void CpModelMapping::PropagateEncodingFromEquivalenceRelations( const IntegerValue rhs(ct.linear().domain(0)); // Make sure the coefficient are positive. - IntegerVariable var1 = Integer(ct.linear().vars(0)); - IntegerVariable var2 = Integer(ct.linear().vars(1)); + IntegerVariable var1 = mapping->Integer(ct.linear().vars(0)); + IntegerVariable var2 = mapping->Integer(ct.linear().vars(1)); IntegerValue coeff1(ct.linear().coeffs(0)); IntegerValue coeff2(ct.linear().coeffs(1)); if (coeff1 < 0) { @@ -655,8 +664,8 @@ void CpModelMapping::PropagateEncodingFromEquivalenceRelations( } } -void CpModelMapping::DetectOptionalVariables(const CpModelProto& model_proto, - Model* m) { +void DetectOptionalVariables(const CpModelProto& model_proto, Model* m) { + auto* mapping = m->GetOrCreate(); const SatParameters& parameters = *(m->GetOrCreate()); if (!parameters.use_optional_variables()) return; if (parameters.enumerate_all_solutions()) return; @@ -724,7 +733,8 @@ void CpModelMapping::DetectOptionalVariables(const CpModelProto& model_proto, ++num_optionals; integer_trail->MarkIntegerVariableAsOptional( - Integer(var), Literal(enforcement_intersection[var].front())); + mapping->Integer(var), + mapping->Literal(enforcement_intersection[var].front())); } VLOG(2) << "Auto-detected " << num_optionals << " optional variables."; } @@ -1359,26 +1369,14 @@ void LoadIntMinConstraint(const ConstraintProto& ct, Model* m) { m->Add(IsEqualToMinOf(min, vars)); } -LinearExpression GetExprFromProto(const LinearExpressionProto& expr_proto, - const CpModelMapping& mapping) { - LinearExpression expr; - expr.vars = mapping.Integers(expr_proto.vars()); - for (int j = 0; j < expr_proto.coeffs_size(); ++j) { - expr.coeffs.push_back(IntegerValue(expr_proto.coeffs(j))); - } - expr.offset = IntegerValue(expr_proto.offset()); - return CanonicalizeExpr(expr); -} - void LoadLinMaxConstraint(const ConstraintProto& ct, Model* m) { auto* mapping = m->GetOrCreate(); - const LinearExpression max = - GetExprFromProto(ct.lin_max().target(), *mapping); + const LinearExpression max = mapping->GetExprFromProto(ct.lin_max().target()); std::vector negated_exprs; negated_exprs.reserve(ct.lin_max().exprs_size()); for (int i = 0; i < ct.lin_max().exprs_size(); ++i) { negated_exprs.push_back( - NegationOf(GetExprFromProto(ct.lin_max().exprs(i), *mapping))); + NegationOf(mapping->GetExprFromProto(ct.lin_max().exprs(i)))); } // TODO(user): Consider replacing the min propagator by max. m->Add(IsEqualToMinOf(NegationOf(max), negated_exprs)); diff --git a/ortools/sat/cp_model_loader.h b/ortools/sat/cp_model_loader.h index 73efefbd0a..d6815e8254 100644 --- a/ortools/sat/cp_model_loader.h +++ b/ortools/sat/cp_model_loader.h @@ -25,231 +25,55 @@ #include "ortools/base/map_util.h" #include "ortools/base/strong_vector.h" #include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/integer.h" #include "ortools/sat/intervals.h" +#include "ortools/sat/linear_relaxation.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" 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. +// Extracts all the used variables in the CpModelProto and creates a +// sat::Model representation for them. More precisely +// - All Boolean variables will be mapped. +// - All Interval variables will be mapped. +// - All non-Boolean variable will have a corresponding IntegerVariable, and +// depending on the view_all_booleans_as_integers, some or all of the +// BooleanVariable will also have an IntegerVariable corresponding to its +// "integer view". // -// This also holds some information used when loading a CpModel proto. -class CpModelMapping { - public: - // Extracts all the used variables in the CpModelProto and creates a - // sat::Model representation for them. More precisely - // - All Boolean variables will be mapped. - // - All Interval variables will be mapped. - // - All non-Boolean variable will have a corresponding IntegerVariable, and - // depending on the view_all_booleans_as_integers, some or all of the - // BooleanVariable will also have an IntegerVariable corresponding to its - // "integer view". - // - // Note(user): We could create IntegerVariable on the fly as they are needed, - // but that loose the original variable order which might be useful in - // heuristics later. - void CreateVariables(const CpModelProto& model_proto, - bool view_all_booleans_as_integers, Model* m); +// Note(user): We could create IntegerVariable on the fly as they are needed, +// but that loose the original variable order which might be useful in +// heuristics later. +void LoadVariables(const CpModelProto& model_proto, + bool view_all_booleans_as_integers, Model* m); - // Automatically detect optional variables. - void DetectOptionalVariables(const CpModelProto& model_proto, Model* m); +// Automatically detect optional variables. +void DetectOptionalVariables(const CpModelProto& model_proto, Model* m); - // Experimental. Loads the symmetry form the proto symmetry field, as long as - // they only involve Booleans. - // - // TODO(user): We currently only have the code for Booleans, it is why we - // currently ignore symmetries involving integer variables. - void LoadBooleanSymmetries(const CpModelProto& model_proto, Model* m); +// Experimental. Loads the symmetry form the proto symmetry field, as long as +// they only involve Booleans. +// +// TODO(user): We currently only have the code for Booleans, it is why we +// currently ignore symmetries involving integer variables. +void LoadBooleanSymmetries(const CpModelProto& model_proto, Model* m); - // Extract the encodings (IntegerVariable <-> Booleans) present in the model. - // This effectively load some linear constraints of size 1 that will be marked - // as already loaded. - void ExtractEncoding(const CpModelProto& model_proto, Model* m); +// Extract the encodings (IntegerVariable <-> Booleans) present in the model. +// This effectively load some linear constraints of size 1 that will be marked +// as already loaded. +void ExtractEncoding(const CpModelProto& model_proto, Model* m); - // Process all affine relations of the form a*X + b*Y == cte. For each - // literals associated to (X >= bound) or (X == value) associate it to its - // corresponding relation on Y. Also do the other side. - // - // TODO(user): In an ideal world, all affine relations like this should be - // removed in the presolve. - void PropagateEncodingFromEquivalenceRelations( - const CpModelProto& model_proto, Model* m); - - // Returns true if the given CpModelProto variable reference refers to a - // Boolean variable. Such variable will always have an associated Literal(), - // but not always an associated Integer(). - bool IsBoolean(int ref) const { - DCHECK_LT(PositiveRef(ref), booleans_.size()); - return booleans_[PositiveRef(ref)] != kNoBooleanVariable; - } - - bool IsInteger(int ref) const { - DCHECK_LT(PositiveRef(ref), integers_.size()); - return integers_[PositiveRef(ref)] != kNoIntegerVariable; - } - - sat::Literal Literal(int ref) const { - DCHECK(IsBoolean(ref)); - return sat::Literal(booleans_[PositiveRef(ref)], RefIsPositive(ref)); - } - - IntegerVariable Integer(int ref) const { - DCHECK(IsInteger(ref)); - const IntegerVariable var = integers_[PositiveRef(ref)]; - return RefIsPositive(ref) ? var : NegationOf(var); - } - - // TODO(user): We could "easily" create an intermediate variable for more - // complex linear expression. We could also identify duplicate expressions to - // not create two identical integer variable. - AffineExpression LoadAffineView(const LinearExpressionProto& exp) const { - CHECK_LE(exp.vars().size(), 1); - if (exp.vars().empty()) { - return AffineExpression(IntegerValue(exp.offset())); - } - return AffineExpression(Integer(exp.vars(0)), IntegerValue(exp.coeffs(0)), - IntegerValue(exp.offset())); - } - - IntervalVariable Interval(int i) const { - CHECK_GE(i, 0); - CHECK_LT(i, intervals_.size()); - CHECK_NE(intervals_[i], kNoIntervalVariable); - return intervals_[i]; - } - - template - std::vector Integers(const List& list) const { - std::vector result; - for (const auto i : list) result.push_back(Integer(i)); - return result; - } - - template - std::vector Literals(const ProtoIndices& indices) const { - std::vector result; - for (const int i : indices) result.push_back(CpModelMapping::Literal(i)); - return result; - } - - template - std::vector Intervals(const ProtoIndices& indices) const { - std::vector result; - for (const int i : indices) result.push_back(Interval(i)); - return result; - } - - // Depending on the option, we will load constraints in stages. This is used - // to detect constraints that are already loaded. For instance the interval - // constraints and the linear constraint of size 1 (encodings) are usually - // loaded first. - bool ConstraintIsAlreadyLoaded(const ConstraintProto* ct) const { - return already_loaded_ct_.contains(ct); - } - - // Returns true if the given constraint is a "half-encoding" constraint. That - // is, if it is of the form (b => size 1 linear) but there is no (<=) side in - // the model. Such constraint are detected while we extract integer encoding - // and are cached here so that we can deal properly with them during the - // linear relaxation. - bool IsHalfEncodingConstraint(const ConstraintProto* ct) const { - return is_half_encoding_ct_.contains(ct); - } - - // Note that both these functions returns positive reference or -1. - int GetProtoVariableFromBooleanVariable(BooleanVariable var) const { - if (var.value() >= reverse_boolean_map_.size()) return -1; - return reverse_boolean_map_[var]; - } - int GetProtoVariableFromIntegerVariable(IntegerVariable var) const { - if (var.value() >= reverse_integer_map_.size()) return -1; - return reverse_integer_map_[var]; - } - - const std::vector& GetVariableMapping() const { - return integers_; - } - - // For logging only, these are not super efficient. - int NumIntegerVariables() const { - int result = 0; - for (const IntegerVariable var : integers_) { - if (var != kNoIntegerVariable) result++; - } - return result; - } - int NumBooleanVariables() const { - int result = 0; - for (const BooleanVariable var : booleans_) { - if (var != kNoBooleanVariable) result++; - } - return result; - } - - // Returns a heuristic set of values that could be created for the given - // variable when the constraints will be loaded. - // Note that the pointer is not stable across calls. - // It returns nullptr if the set is empty. - const absl::flat_hash_set& PotentialEncodedValues(int var) { - const auto& it = variables_to_encoded_values_.find(var); - if (it != variables_to_encoded_values_.end()) { - return it->second; - } - return empty_set_; - } - - // Returns the number of variables in the loaded proto. - int NumProtoVariables() const { return integers_.size(); } - - private: - // Note that only the variables used by at least one constraint will be - // created, the other will have a kNo[Integer,Interval,Boolean]VariableValue. - std::vector integers_; - std::vector intervals_; - std::vector booleans_; - - // Recover from a IntervalVariable/BooleanVariable its associated CpModelProto - // index. The value of -1 is used to indicate that there is no correspondence - // (i.e. this variable is only used internally). - absl::StrongVector reverse_boolean_map_; - absl::StrongVector reverse_integer_map_; - - // Set of constraints to ignore because they were already dealt with by - // ExtractEncoding(). - absl::flat_hash_set already_loaded_ct_; - absl::flat_hash_set is_half_encoding_ct_; - - absl::flat_hash_map> - variables_to_encoded_values_; - const absl::flat_hash_set empty_set_; -}; +// Process all affine relations of the form a*X + b*Y == cte. For each +// literals associated to (X >= bound) or (X == value) associate it to its +// corresponding relation on Y. Also do the other side. +// +// TODO(user): In an ideal world, all affine relations like this should be +// removed in the presolve. +void PropagateEncodingFromEquivalenceRelations(const CpModelProto& model_proto, + Model* m); // Inspects the model and use some heuristic to decide which variable, if any, // should be fully encoded. Note that some constraints like the element or table @@ -297,9 +121,6 @@ void LoadRoutesConstraint(const ConstraintProto& ct, Model* m); void LoadCircuitCoveringConstraint(const ConstraintProto& ct, Model* m); void LoadInverseConstraint(const ConstraintProto& ct, Model* m); -LinearExpression GetExprFromProto(const LinearExpressionProto& expr_proto, - const CpModelMapping& mapping); - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 8c2384bab5..dcc1d82374 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -44,6 +44,7 @@ #include "ortools/sat/cp_model_checker.h" #include "ortools/sat/cp_model_expand.h" #include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/cp_model_objective.h" #include "ortools/sat/cp_model_symmetries.h" #include "ortools/sat/cp_model_utils.h" @@ -2341,18 +2342,22 @@ bool CpModelPresolver::PresolveLinearOnBooleans(ConstraintProto* ct) { namespace { -void AddLinearExpression(int64_t coeff, const LinearExpressionProto& exp, - LinearConstraintProto* linear_ct) { - const int size = exp.vars().size(); - for (int i = 0; i < size; ++i) { - linear_ct->add_vars(exp.vars(i)); - linear_ct->add_coeffs(coeff * exp.coeffs(i)); - } - if (exp.offset() != 0) { - FillDomainInProto(ReadDomainFromProto(*linear_ct) - .AdditionWith(Domain(-coeff * exp.offset())), - linear_ct); - } +void AddLinearConstraintFromInterval(const ConstraintProto& ct, + PresolveContext* context) { + const int start = ct.interval().start(); + const int end = ct.interval().end(); + const int size = ct.interval().size(); + ConstraintProto* new_ct = context->working_model->add_constraints(); + *(new_ct->mutable_enforcement_literal()) = ct.enforcement_literal(); + new_ct->mutable_linear()->add_domain(0); + new_ct->mutable_linear()->add_domain(0); + new_ct->mutable_linear()->add_vars(start); + new_ct->mutable_linear()->add_coeffs(1); + new_ct->mutable_linear()->add_vars(size); + new_ct->mutable_linear()->add_coeffs(1); + new_ct->mutable_linear()->add_vars(end); + new_ct->mutable_linear()->add_coeffs(-1); + context->UpdateNewConstraintsVariableUsage(); } } // namespace @@ -2392,21 +2397,7 @@ bool CpModelPresolver::PresolveInterval(int c, ConstraintProto* ct) { if (context_->IntervalUsage(c) == 0) { if (!ct->interval().has_start_view()) { - // Convert to linear. - const int start = ct->interval().start(); - const int end = ct->interval().end(); - const int size = ct->interval().size(); - ConstraintProto* new_ct = context_->working_model->add_constraints(); - *(new_ct->mutable_enforcement_literal()) = ct->enforcement_literal(); - new_ct->mutable_linear()->add_domain(0); - new_ct->mutable_linear()->add_domain(0); - new_ct->mutable_linear()->add_vars(start); - new_ct->mutable_linear()->add_coeffs(1); - new_ct->mutable_linear()->add_vars(size); - new_ct->mutable_linear()->add_coeffs(1); - new_ct->mutable_linear()->add_vars(end); - new_ct->mutable_linear()->add_coeffs(-1); - context_->UpdateNewConstraintsVariableUsage(); + AddLinearConstraintFromInterval(*ct, context_); } context_->UpdateRuleStats("interval: unused, converted to linear"); return RemoveConstraint(ct); @@ -2423,6 +2414,10 @@ bool CpModelPresolver::PresolveInterval(int c, ConstraintProto* ct) { if (!ct->interval().has_start_view()) { changed = true; + // Add a linear constraint. Our new format require a separate linear + // constraint which allow us to reuse all the propagation code. + AddLinearConstraintFromInterval(*ct, context_); + // Fill the view fields. interval->mutable_start_view()->add_vars(interval->start()); interval->mutable_start_view()->add_coeffs(1); @@ -2434,20 +2429,6 @@ bool CpModelPresolver::PresolveInterval(int c, ConstraintProto* ct) { interval->mutable_end_view()->add_coeffs(1); interval->mutable_end_view()->set_offset(0); - // Create a new linear constraint to detect affine relation and propagate - // the domain properly. - // - // Note(user): Aving an extra constraint is not super clean, but reusing - // the code is a must. - ConstraintProto* new_ct = context_->working_model->add_constraints(); - *(new_ct->mutable_enforcement_literal()) = ct->enforcement_literal(); - new_ct->mutable_linear()->add_domain(0); - new_ct->mutable_linear()->add_domain(0); - AddLinearExpression(1, interval->start_view(), new_ct->mutable_linear()); - AddLinearExpression(1, interval->size_view(), new_ct->mutable_linear()); - AddLinearExpression(-1, interval->end_view(), new_ct->mutable_linear()); - context_->UpdateNewConstraintsVariableUsage(); - // Set the old fields to their default. Not really needed. interval->set_start(0); interval->set_size(0); @@ -2461,17 +2442,6 @@ bool CpModelPresolver::PresolveInterval(int c, ConstraintProto* ct) { return changed; } - // If the interval is of fixed size, we can add the corresponsing affine - // relation to our pool. - // - // TODO(user): This will currently add another linear relation to the proto - // in addition to the interval at the end of the presolve though. - if (/* DISABLES CODE */ (false) && ct->enforcement_literal().empty() && - context_->IsFixed(ct->interval().size())) { - context_->StoreAffineRelation(ct->interval().end(), ct->interval().start(), - 1, context_->MinOf(ct->interval().size())); - } - // This never change the constraint-variable graph. return false; } @@ -3973,9 +3943,8 @@ void CpModelPresolver::Probe() { // Important: Because the model_proto do not contains affine relation or the // objective, we cannot call DetectOptionalVariables() ! This might wrongly // detect optionality and derive bad conclusion. - mapping->CreateVariables(model_proto, /*view_all_booleans_as_integers=*/false, - &model); - mapping->ExtractEncoding(model_proto, &model); + LoadVariables(model_proto, /*view_all_booleans_as_integers=*/false, &model); + ExtractEncoding(model_proto, &model); auto* sat_solver = model.GetOrCreate(); for (const ConstraintProto& ct : model_proto.constraints()) { if (mapping->ConstraintIsAlreadyLoaded(&ct)) continue; diff --git a/ortools/sat/cp_model_search.h b/ortools/sat/cp_model_search.h index 29851fc5df..314c1f7b5d 100644 --- a/ortools/sat/cp_model_search.h +++ b/ortools/sat/cp_model_search.h @@ -20,7 +20,7 @@ #include "ortools/base/integral_types.h" #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_search.h" #include "ortools/sat/model.h" diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 9ec80c3ea6..1c3e7efb8b 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -453,316 +453,13 @@ IntegerVariable GetOrCreateVariableGreaterOrEqualToSumOf( return new_var; } -void TryToAddCutGenerators(const CpModelProto& model_proto, - const ConstraintProto& ct, Model* m, - LinearRelaxation* relaxation) { - const int linearization_level = - m->GetOrCreate()->linearization_level(); - auto* mapping = m->GetOrCreate(); - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kCircuit && - linearization_level > 1) { - std::vector tails(ct.circuit().tails().begin(), - ct.circuit().tails().end()); - std::vector heads(ct.circuit().heads().begin(), - ct.circuit().heads().end()); - std::vector literals = mapping->Literals(ct.circuit().literals()); - const int num_nodes = ReindexArcs(&tails, &heads); - - relaxation->cut_generators.push_back( - CreateStronglyConnectedGraphCutGenerator(num_nodes, tails, heads, - literals, m)); - } - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kRoutes && - linearization_level > 1) { - std::vector tails(ct.routes().tails().begin(), - ct.routes().tails().end()); - std::vector heads(ct.routes().heads().begin(), - ct.routes().heads().end()); - std::vector literals = mapping->Literals(ct.routes().literals()); - - int num_nodes = 0; - for (int i = 0; i < ct.routes().tails_size(); ++i) { - num_nodes = std::max(num_nodes, 1 + ct.routes().tails(i)); - num_nodes = std::max(num_nodes, 1 + ct.routes().heads(i)); - } - if (ct.routes().demands().empty() || ct.routes().capacity() == 0) { - relaxation->cut_generators.push_back( - CreateStronglyConnectedGraphCutGenerator(num_nodes, tails, heads, - literals, m)); - } else { - const std::vector demands(ct.routes().demands().begin(), - ct.routes().demands().end()); - relaxation->cut_generators.push_back( - CreateCVRPCutGenerator(num_nodes, tails, heads, literals, demands, - ct.routes().capacity(), m)); - } - } - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntProd) { - if (HasEnforcementLiteral(ct)) return; - if (ct.int_prod().vars_size() != 2) return; - - // Constraint is z == x * y. - - IntegerVariable z = mapping->Integer(ct.int_prod().target()); - IntegerVariable x = mapping->Integer(ct.int_prod().vars(0)); - IntegerVariable y = mapping->Integer(ct.int_prod().vars(1)); - - IntegerTrail* const integer_trail = m->GetOrCreate(); - IntegerValue x_lb = integer_trail->LowerBound(x); - IntegerValue x_ub = integer_trail->UpperBound(x); - IntegerValue y_lb = integer_trail->LowerBound(y); - IntegerValue y_ub = integer_trail->UpperBound(y); - - if (x == y) { - // We currently only support variables with non-negative domains. - if (x_lb < 0 && x_ub > 0) return; - - // Change the sigh of x if its domain is non-positive. - if (x_ub <= 0) { - x = NegationOf(x); - } - - relaxation->cut_generators.push_back(CreateSquareCutGenerator(z, x, m)); - } else { - // We currently only support variables with non-negative domains. - if (x_lb < 0 && x_ub > 0) return; - if (y_lb < 0 && y_ub > 0) return; - - // Change signs to return to the case where all variables are a domain - // with non negative values only. - if (x_ub <= 0) { - x = NegationOf(x); - z = NegationOf(z); - } - if (y_ub <= 0) { - y = NegationOf(y); - z = NegationOf(z); - } - - relaxation->cut_generators.push_back( - CreatePositiveMultiplicationCutGenerator(z, x, y, m)); - } - } - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kAllDiff) { - if (linearization_level < 2) return; - if (HasEnforcementLiteral(ct)) return; - const int num_vars = ct.all_diff().vars_size(); - if (num_vars <= m->GetOrCreate()->max_all_diff_cut_size()) { - std::vector vars = - mapping->Integers(ct.all_diff().vars()); - relaxation->cut_generators.push_back( - CreateAllDifferentCutGenerator(vars, m)); - } - } - - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kCumulative) { - if (linearization_level < 2) return; - if (HasEnforcementLiteral(ct)) return; - - std::vector demands = - mapping->Integers(ct.cumulative().demands()); - std::vector intervals = - mapping->Intervals(ct.cumulative().intervals()); - const IntegerVariable capacity = - mapping->Integer(ct.cumulative().capacity()); - relaxation->cut_generators.push_back( - CreateCumulativeOverlappingCutGenerator(intervals, capacity, demands, - m)); - relaxation->cut_generators.push_back( - CreateCumulativeEnergyCutGenerator(intervals, capacity, demands, m)); - relaxation->cut_generators.push_back( - CreateCumulativeCompletionTimeCutGenerator(intervals, capacity, demands, - m)); - } - - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kNoOverlap) { - if (linearization_level < 2) return; - if (HasEnforcementLiteral(ct)) return; - std::vector intervals = - mapping->Intervals(ct.no_overlap().intervals()); - relaxation->cut_generators.push_back( - CreateNoOverlapEnergyCutGenerator(intervals, m)); - relaxation->cut_generators.push_back( - CreateNoOverlapPrecedenceCutGenerator(intervals, m)); - relaxation->cut_generators.push_back( - CreateNoOverlapCompletionTimeCutGenerator(intervals, m)); - } - - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kNoOverlap2D) { - if (linearization_level < 2) return; - if (HasEnforcementLiteral(ct)) return; - std::vector x_intervals = - mapping->Intervals(ct.no_overlap_2d().x_intervals()); - std::vector y_intervals = - mapping->Intervals(ct.no_overlap_2d().y_intervals()); - relaxation->cut_generators.push_back( - CreateNoOverlap2dCompletionTimeCutGenerator(x_intervals, y_intervals, - m)); - } - - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinMax) { - if (!m->GetOrCreate()->add_lin_max_cuts()) return; - if (HasEnforcementLiteral(ct)) return; - - // TODO(user): Support linearization of general target expression. - if (ct.lin_max().target().vars_size() != 1) return; - if (ct.lin_max().target().coeffs(0) != 1) return; - - const IntegerVariable target = - mapping->Integer(ct.lin_max().target().vars(0)); - std::vector exprs; - exprs.reserve(ct.lin_max().exprs_size()); - for (int i = 0; i < ct.lin_max().exprs_size(); ++i) { - // Note: Cut generator requires all expressions to contain only positive - // vars. - exprs.push_back( - PositiveVarExpr(GetExprFromProto(ct.lin_max().exprs(i), *mapping))); - } - - // Add initial big-M linear relaxation. - // z_vars[i] == 1 <=> target = exprs[i]. - const std::vector z_vars = - AppendLinMaxRelaxation(target, exprs, m, relaxation); - - if (linearization_level >= 2) { - relaxation->cut_generators.push_back( - CreateLinMaxCutGenerator(target, exprs, z_vars, m)); - } - } -} - } // namespace -LinearRelaxation ComputeLinearRelaxation(const CpModelProto& model_proto, - int linearization_level, Model* m) { - LinearRelaxation relaxation; - - // Linearize the constraints. - absl::flat_hash_set used_integer_variable; - - auto* mapping = m->GetOrCreate(); - auto* encoder = m->GetOrCreate(); - auto* trail = m->GetOrCreate(); - for (const auto& ct : model_proto.constraints()) { - // Make sure the literals from a circuit constraint always have a view. - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kCircuit) { - for (const int ref : ct.circuit().literals()) { - const Literal l = mapping->Literal(ref); - if (!encoder->LiteralOrNegationHasView(l)) { - m->Add(NewIntegerVariableFromLiteral(l)); - } - } - } - - // For now, we skip any constraint with literals that do not have an integer - // view. Ideally it should be up to the constraint to decide if creating a - // view is worth it. - // - // TODO(user): It should be possible to speed this up if needed. - const IndexReferences refs = GetReferencesUsedByConstraint(ct); - bool ok = true; - for (const int literal_ref : refs.literals) { - const Literal literal = mapping->Literal(literal_ref); - if (trail->Assignment().LiteralIsAssigned(literal)) { - // Create a view to the constant 0 or 1. - m->Add(NewIntegerVariableFromLiteral(literal)); - } else if (!encoder->LiteralOrNegationHasView(literal)) { - ok = false; - break; - } - } - if (!ok) continue; - - TryToLinearizeConstraint(model_proto, ct, m, linearization_level, - &relaxation); - TryToAddCutGenerators(model_proto, ct, m, &relaxation); - } - - // Linearize the encoding of variable that are fully encoded in the proto. - int num_full_encoding_relaxations = 0; - int num_partial_encoding_relaxations = 0; - for (int i = 0; i < model_proto.variables_size(); ++i) { - if (mapping->IsBoolean(i)) continue; - - const IntegerVariable var = mapping->Integer(i); - if (m->Get(IsFixed(var))) continue; - - // TODO(user): This different encoding for the partial variable might be - // better (less LP constraints), but we do need more investigation to - // decide. - if (/* DISABLES CODE */ (false)) { - AppendPartialEncodingRelaxation(var, *m, &relaxation); - continue; - } - - if (encoder->VariableIsFullyEncoded(var)) { - if (AppendFullEncodingRelaxation(var, *m, &relaxation)) { - ++num_full_encoding_relaxations; - continue; - } - } - - // Even if the variable is fully encoded, sometimes not all its associated - // literal have a view (if they are not part of the original model for - // instance). - // - // TODO(user): Should we add them to the LP anyway? this isn't clear as - // we can sometimes create a lot of Booleans like this. - const int old = relaxation.linear_constraints.size(); - AppendPartialGreaterThanEncodingRelaxation(var, *m, &relaxation); - if (relaxation.linear_constraints.size() > old) { - ++num_partial_encoding_relaxations; - } - } - - if (!m->GetOrCreate()->FinishPropagation()) return relaxation; - - // Linearize the at most one constraints. Note that we transform them - // into maximum "at most one" first and we removes redundant ones. - m->GetOrCreate()->TransformIntoMaxCliques( - &relaxation.at_most_ones); - for (const std::vector& at_most_one : relaxation.at_most_ones) { - if (at_most_one.empty()) continue; - - LinearConstraintBuilder lc(m, kMinIntegerValue, IntegerValue(1)); - for (const Literal literal : at_most_one) { - // Note that it is okay to simply ignore the literal if it has no - // integer view. - const bool unused ABSL_ATTRIBUTE_UNUSED = - lc.AddLiteralTerm(literal, IntegerValue(1)); - } - relaxation.linear_constraints.push_back(lc.Build()); - } - - // We converted all at_most_one to LP constraints, so we need to clear them - // so that we don't do extra work in the connected component computation. - relaxation.at_most_ones.clear(); - - // Remove size one LP constraints, they are not useful. - { - int new_size = 0; - for (int i = 0; i < relaxation.linear_constraints.size(); ++i) { - if (relaxation.linear_constraints[i].vars.size() <= 1) continue; - std::swap(relaxation.linear_constraints[new_size++], - relaxation.linear_constraints[i]); - } - relaxation.linear_constraints.resize(new_size); - } - - VLOG(3) << "num_full_encoding_relaxations: " << num_full_encoding_relaxations; - VLOG(3) << "num_partial_encoding_relaxations: " - << num_partial_encoding_relaxations; - VLOG(3) << relaxation.linear_constraints.size() - << " constraints in the LP relaxation."; - VLOG(3) << relaxation.cut_generators.size() << " cuts generators."; - return relaxation; -} - // Adds one LinearProgrammingConstraint per connected component of the model. IntegerVariable AddLPConstraints(const CpModelProto& model_proto, int linearization_level, Model* m) { - const LinearRelaxation relaxation = - ComputeLinearRelaxation(model_proto, linearization_level, m); + LinearRelaxation relaxation; + ComputeLinearRelaxation(model_proto, linearization_level, m, &relaxation); // The bipartite graph of LP constraints might be disconnected: // make a partition of the variables into connected components. @@ -1239,19 +936,19 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { (parameters.linearization_level() >= 2) || (parameters.search_branching() == SatParameters::FIXED_SEARCH && model_proto.search_strategy().empty()); - mapping->CreateVariables(model_proto, view_all_booleans_as_integers, model); - mapping->DetectOptionalVariables(model_proto, model); + LoadVariables(model_proto, view_all_booleans_as_integers, model); + DetectOptionalVariables(model_proto, model); // TODO(user): The core algo and symmetries seems to be problematic in some // cases. See for instance: neos-691058.mps.gz. This is probably because as // we modify the model, our symmetry might be wrong? investigate. if (!parameters.optimize_with_core() && parameters.symmetry_level() > 1 && !parameters.enumerate_all_solutions()) { - mapping->LoadBooleanSymmetries(model_proto, model); + LoadBooleanSymmetries(model_proto, model); } - mapping->ExtractEncoding(model_proto, model); - mapping->PropagateEncodingFromEquivalenceRelations(model_proto, model); + ExtractEncoding(model_proto, model); + PropagateEncodingFromEquivalenceRelations(model_proto, model); // Check the model is still feasible before continuing. if (sat_solver->IsModelUnsat()) return unsat(); @@ -1322,8 +1019,9 @@ void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) { if (parameters.linearization_level() == 0) return; // Add linear constraints to Feasibility Pump. - const LinearRelaxation relaxation = ComputeLinearRelaxation( - model_proto, parameters.linearization_level(), model); + LinearRelaxation relaxation; + ComputeLinearRelaxation(model_proto, parameters.linearization_level(), model, + &relaxation); const int num_lp_constraints = relaxation.linear_constraints.size(); if (num_lp_constraints == 0) return; auto* feasibility_pump = model->GetOrCreate(); diff --git a/ortools/sat/cumulative_energy.cc b/ortools/sat/cumulative_energy.cc index e17ed76dde..378a002ca8 100644 --- a/ortools/sat/cumulative_energy.cc +++ b/ortools/sat/cumulative_energy.cc @@ -105,7 +105,7 @@ void CumulativeEnergyConstraint::RegisterWith(GenericLiteralWatcher* watcher) { bool CumulativeEnergyConstraint::Propagate() { // This only uses one time direction, but the helper might be used elsewhere. // TODO(user): just keep the current direction? - helper_->SynchronizeAndSetTimeDirection(true); + if (!helper_->SynchronizeAndSetTimeDirection(true)) return false; const IntegerValue capacity_max = integer_trail_->UpperBound(capacity_); // TODO(user): force capacity_max >= 0, fail/remove optionals when 0. diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 385473f3a5..791f74c2f5 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -477,7 +477,7 @@ CutGenerator CreateKnapsackCoverCutGenerator( // TODO(user): When we use implied-bound substitution, we might still infer // an interesting cut even if all variables are integer. See if we still // want to skip all such constraints. - if (AllVarsTakeIntegerValue(vars, lp_values)) return; + if (AllVarsTakeIntegerValue(vars, lp_values)) return true; KnapsackSolverForCuts knapsack_solver( "Knapsack on demand cover cut generator"); @@ -611,6 +611,7 @@ CutGenerator CreateKnapsackCoverCutGenerator( if (skipped_constraints > 0) { VLOG(2) << "Skipped constraints: " << skipped_constraints; } + return true; }; return result; @@ -683,7 +684,7 @@ std::function GetSuperAdditiveRoundingFunction( // as low as 2 could lead to the better cut (this is exactly the Letchford & // Lodi function). // - // Another intersting fact, is that if we want to compute the maximum alpha + // Another interesting fact, is that if we want to compute the maximum alpha // for a constraint with 2 terms like: // divisor * Y + (ratio * divisor + remainder) * X // <= rhs_ratio * divisor + rhs_remainder @@ -1229,7 +1230,7 @@ bool CoverCutHelper::TrySimpleKnapsack( // Note(user): past this point, now that a given "base" cover has been chosen, // we basically compute the cut (of the form sum X <= bound) with the maximum // possible violation. Note also that we lift as much as possible, so we don't - // necessarilly optimize for the cut efficacity though. But we do get a + // necessarily optimize for the cut efficacity though. But we do get a // stronger cut. if (rhs >= 0) return false; if (new_size == 0) return false; @@ -1356,7 +1357,7 @@ CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z, if (CapProd(x_ub, y_ub) >= kMaxSafeInteger) { VLOG(3) << "Potential overflow in PositiveMultiplicationCutGenerator"; - return; + return true; } const double x_lp_value = lp_values[x]; @@ -1425,6 +1426,7 @@ CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z, try_add_above_cut(y_ub, x_ub, x_ub * y_ub); try_add_below_cut(y_ub, x_lb, x_lb * y_ub); try_add_below_cut(y_lb, x_ub, x_ub * y_lb); + return true; }; return result; @@ -1443,10 +1445,10 @@ CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x, const int64_t x_ub = integer_trail->LevelZeroUpperBound(x).value(); const int64_t x_lb = integer_trail->LevelZeroLowerBound(x).value(); - if (x_lb == x_ub) return; + if (x_lb == x_ub) return true; // Check for potential overflows. - if (x_ub > (int64_t{1} << 31)) return; + if (x_ub > (int64_t{1} << 31)) return true; DCHECK_GE(x_lb, 0); const double y_lp_value = lp_values[y]; @@ -1492,6 +1494,7 @@ CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x, below_cut.ub = kMaxIntegerValue; manager->AddCut(below_cut, "SquareLower", lp_values); } + return true; }; return result; @@ -1829,7 +1832,7 @@ CutGenerator CreateAllDifferentCutGenerator( // These cuts work at all levels but the generator adds too many cuts on // some instances and degrade the performance so we only use it at level // 0. - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; std::vector> sorted_vars; for (const IntegerVariable var : vars) { if (integer_trail->LevelZeroLowerBound(var) == @@ -1845,6 +1848,7 @@ CutGenerator CreateAllDifferentCutGenerator( std::reverse(sorted_vars.begin(), sorted_vars.end()); TryToGenerateAllDiffCut(sorted_vars, *integer_trail, lp_values, manager); + return true; }; VLOG(1) << "Created all_diff cut generator of size: " << vars.size(); return result; @@ -1978,6 +1982,7 @@ CutGenerator CreateLinMaxCutGenerator( if (violation > 1e-2) { manager->AddCut(cut.Build(), "LinMax", lp_values); } + return true; }; return result; } @@ -2015,7 +2020,7 @@ void AddIntegerVariableFromIntervals(SchedulingConstraintHelper* helper, gtl::STLSortAndRemoveDuplicates(vars); } -std::function&, +std::function&, LinearConstraintManager*)> GenerateCumulativeCut(const std::string& cut_name, SchedulingConstraintHelper* helper, @@ -2028,7 +2033,7 @@ GenerateCumulativeCut(const std::string& cut_name, return [capacity, demands, trail, integer_trail, helper, model, cut_name, encoder](const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; const auto demand_is_fixed = [integer_trail, &demands](int i) { return demands.empty() || integer_trail->IsFixed(demands[i]); @@ -2049,7 +2054,7 @@ GenerateCumulativeCut(const std::string& cut_name, } } - if (active_intervals.size() < 2) return; + if (active_intervals.size() < 2) return true; std::sort(active_intervals.begin(), active_intervals.end(), [helper](int a, int b) { @@ -2193,6 +2198,7 @@ GenerateCumulativeCut(const std::string& cut_name, manager->AddCut(cut.Build(), cut_name, lp_values); } } + return true; }; } @@ -2243,7 +2249,7 @@ CutGenerator CreateCumulativeOverlappingCutGenerator( [helper, capacity, demands, trail, integer_trail, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; std::vector events; // Iterate through the intervals. If start_max < end_min, the demand @@ -2325,6 +2331,7 @@ CutGenerator CreateCumulativeOverlappingCutGenerator( cut_events.resize(new_size); added_positive_event = false; } + return true; }; return result; } @@ -2362,7 +2369,7 @@ CutGenerator CreateNoOverlapPrecedenceCutGenerator( [trail, helper, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; // TODO(user): We can do much better in term of complexity: // Sort all tasks by min start time, loop other them 1 by 1, @@ -2405,6 +2412,7 @@ CutGenerator CreateNoOverlapPrecedenceCutGenerator( } } } + return true; }; return result; @@ -2525,7 +2533,7 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( [trail, helper, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; auto generate_cuts = [&lp_values, model, manager, helper](const std::string& cut_name) { @@ -2544,10 +2552,11 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( model, manager); }; - helper->SynchronizeAndSetTimeDirection(true); + if (!helper->SynchronizeAndSetTimeDirection(true)) return false; generate_cuts("NoOverlapCompletionTime"); - helper->SynchronizeAndSetTimeDirection(false); + if (!helper->SynchronizeAndSetTimeDirection(false)) return false; generate_cuts("NoOverlapCompletionTimeMirror"); + return true; }; return result; } @@ -2573,7 +2582,7 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( [trail, integer_trail, helper, demands, capacity, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; const IntegerValue capacity_max = integer_trail->UpperBound(capacity); auto generate_cuts = [&lp_values, model, manager, helper, capacity_max, @@ -2596,10 +2605,11 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( model, manager); }; - helper->SynchronizeAndSetTimeDirection(true); + if (!helper->SynchronizeAndSetTimeDirection(true)) return false; generate_cuts("CumulativeCompletionTime"); - helper->SynchronizeAndSetTimeDirection(false); + if (!helper->SynchronizeAndSetTimeDirection(false)) return false; generate_cuts("CumulativeCompletionTimeMirror"); + return true; }; return result; } @@ -2625,10 +2635,10 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( [trail, x_helper, y_helper, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { - if (trail->CurrentDecisionLevel() > 0) return; + if (trail->CurrentDecisionLevel() > 0) return true; - x_helper->SynchronizeAndSetTimeDirection(true); - y_helper->SynchronizeAndSetTimeDirection(true); + if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false; + if (!y_helper->SynchronizeAndSetTimeDirection(true)) return false; const int num_boxes = x_helper->NumTasks(); std::vector active_boxes; @@ -2639,17 +2649,20 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( if (cached_areas[box] == 0) continue; if (!y_helper->IsPresent(box) || !y_helper->IsPresent(box)) continue; - // TODO(user): Also consider shifted end max. + // TODO(user): It might be possible/better to use some shifted value + // here, but for now this code is not in the hot spot, so better be + // defensive and only do connected components on really disjoint + // boxes. Rectangle& rectangle = cached_rectangles[box]; - rectangle.x_min = x_helper->ShiftedStartMin(box); + rectangle.x_min = x_helper->StartMin(box); rectangle.x_max = x_helper->EndMax(box); - rectangle.y_min = y_helper->ShiftedStartMin(box); + rectangle.y_min = y_helper->StartMin(box); rectangle.y_max = y_helper->EndMax(box); active_boxes.push_back(box); } - if (active_boxes.size() <= 1) return; + if (active_boxes.size() <= 1) return true; std::vector> components = GetOverlappingRectangleComponents(cached_rectangles, @@ -2677,15 +2690,16 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( model, manager); }; - x_helper->SynchronizeAndSetTimeDirection(true); - y_helper->SynchronizeAndSetTimeDirection(true); + if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false; + if (!y_helper->SynchronizeAndSetTimeDirection(true)) return false; generate_cuts("NoOverlap2dXCompletionTime", x_helper, y_helper); generate_cuts("NoOverlap2dYCompletionTime", y_helper, x_helper); - x_helper->SynchronizeAndSetTimeDirection(false); - y_helper->SynchronizeAndSetTimeDirection(false); + if (!x_helper->SynchronizeAndSetTimeDirection(false)) return false; + if (!y_helper->SynchronizeAndSetTimeDirection(false)) return false; generate_cuts("NoOverlap2dXCompletionTimeMirror", x_helper, y_helper); generate_cuts("NoOverlap2dYCompletionTimeMirror", y_helper, x_helper); } + return true; }; return result; } @@ -2746,6 +2760,7 @@ CutGenerator CreateCliqueCutGenerator( manager->AddCut(builder.Build(), "clique", lp_values); } + return true; }; return result; } diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 7310298c4b..1661425e21 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -40,7 +40,7 @@ namespace sat { // - Only add cuts in term of the same variables or their negation. struct CutGenerator { std::vector vars; - std::function& lp_values, LinearConstraintManager* manager)> generate_cuts; diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index a4aaced1de..50dd7a6a6b 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -116,8 +116,8 @@ NonOverlappingRectanglesEnergyPropagator:: bool NonOverlappingRectanglesEnergyPropagator::Propagate() { const int num_boxes = x_.NumTasks(); - x_.SynchronizeAndSetTimeDirection(true); - y_.SynchronizeAndSetTimeDirection(true); + if (!x_.SynchronizeAndSetTimeDirection(true)) return false; + if (!y_.SynchronizeAndSetTimeDirection(true)) return false; active_boxes_.clear(); cached_areas_.resize(num_boxes); @@ -127,12 +127,16 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { if (cached_areas_[box] == 0) continue; if (!x_.IsPresent(box) || !y_.IsPresent(box)) continue; - // TODO(user): Also consider shifted end max. + // The code needs the size min to be larger or equal to the mandatory part + // for it to works correctly. This is always enforced by the helper. + DCHECK_GE(x_.SizeMin(box), x_.EndMin(box) - x_.StartMax(box)); + DCHECK_GE(y_.SizeMin(box), y_.EndMin(box) - y_.StartMax(box)); + Rectangle& rectangle = cached_rectangles_[box]; rectangle.x_min = x_.ShiftedStartMin(box); - rectangle.x_max = x_.EndMax(box); + rectangle.x_max = x_.ShiftedEndMax(box); rectangle.y_min = y_.ShiftedStartMin(box); - rectangle.y_max = y_.EndMax(box); + rectangle.y_max = y_.ShiftedEndMax(box); active_boxes_.push_back(box); } @@ -181,9 +185,9 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { int NonOverlappingRectanglesEnergyPropagator::RegisterWith( GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - x_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, + x_.WatchAllTasks(id, watcher, /*watch_start_max=*/true, /*watch_end_max=*/true); - y_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, + y_.WatchAllTasks(id, watcher, /*watch_start_max=*/true, /*watch_end_max=*/true); return id; } @@ -220,10 +224,8 @@ bool NonOverlappingRectanglesEnergyPropagator::FailWhenEnergyIsTooLarge( IntegerValue sum_of_areas = cached_areas_[box]; const auto add_box_energy_in_rectangle_reason = [&](int b) { - x_.AddEnergyAfterReason(b, x_.SizeMin(b), area.x_min); - x_.AddEndMaxReason(b, area.x_max); - y_.AddEnergyAfterReason(b, y_.SizeMin(b), area.y_min); - y_.AddEndMaxReason(b, area.y_max); + x_.AddEnergyMinInIntervalReason(b, area.x_min, area.x_max); + y_.AddEnergyMinInIntervalReason(b, area.y_min, area.y_max); }; for (int i = 0; i < neighbors_.size(); ++i) { @@ -457,8 +459,10 @@ bool NonOverlappingRectanglesDisjunctivePropagator:: // And finally propagate. // TODO(user): Sorting of boxes seems influential on the performance. Test. for (const absl::Span boxes : boxes_to_propagate_) { - x_.ResetFromSubset(x, boxes); - y_.ResetFromSubset(y, boxes); + x_.ClearOtherHelper(); + y_.ClearOtherHelper(); + if (!x_.ResetFromSubset(x, boxes)) return false; + if (!y_.ResetFromSubset(y, boxes)) return false; // Collect the common overlapping coordinates of all boxes. IntegerValue lb(std::numeric_limits::min()); @@ -489,8 +493,8 @@ bool NonOverlappingRectanglesDisjunctivePropagator:: } bool NonOverlappingRectanglesDisjunctivePropagator::Propagate() { - global_x_.SynchronizeAndSetTimeDirection(true); - global_y_.SynchronizeAndSetTimeDirection(true); + if (!global_x_.SynchronizeAndSetTimeDirection(true)) return false; + if (!global_y_.SynchronizeAndSetTimeDirection(true)) return false; std::function inner_propagate; if (watcher_->GetCurrentId() == fast_id_) { diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc index 19e0e373b5..83929c2bd3 100644 --- a/ortools/sat/diffn_util.cc +++ b/ortools/sat/diffn_util.cc @@ -13,7 +13,6 @@ #include "ortools/sat/diffn_util.h" -#include "absl/algorithm/container.h" #include "ortools/base/stl_util.h" namespace operations_research { @@ -58,16 +57,14 @@ bool ReportEnergyConflict(Rectangle bounding_box, absl::Span boxes, IntegerValue total_energy(0); for (const int b : boxes) { const IntegerValue x_min = x->ShiftedStartMin(b); - const IntegerValue x_max = x->EndMax(b); + const IntegerValue x_max = x->ShiftedEndMax(b); if (x_min < bounding_box.x_min || x_max > bounding_box.x_max) continue; const IntegerValue y_min = y->ShiftedStartMin(b); - const IntegerValue y_max = y->EndMax(b); + const IntegerValue y_max = y->ShiftedEndMax(b); if (y_min < bounding_box.y_min || y_max > bounding_box.y_max) continue; - x->AddEnergyAfterReason(b, x->SizeMin(b), bounding_box.x_min); - x->AddEndMaxReason(b, bounding_box.x_max); - y->AddEnergyAfterReason(b, y->SizeMin(b), bounding_box.y_min); - y->AddEndMaxReason(b, bounding_box.y_max); + x->AddEnergyMinInIntervalReason(b, bounding_box.x_min, bounding_box.x_max); + y->AddEnergyMinInIntervalReason(b, bounding_box.y_min, bounding_box.y_max); total_energy += x->SizeMin(b) * y->SizeMin(b); diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 8223711f51..84e68ecf73 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -251,7 +251,7 @@ IntegerValue TaskSet::ComputeEndMin(int task_to_ignore, bool DisjunctiveWithTwoItems::Propagate() { DCHECK_EQ(helper_->NumTasks(), 2); - helper_->SynchronizeAndSetTimeDirection(true); + if (!helper_->SynchronizeAndSetTimeDirection(true)) return false; // We can't propagate anything if one of the interval is absent for sure. if (helper_->IsAbsent(0) || helper_->IsAbsent(1)) return true; @@ -342,7 +342,7 @@ void CombinedDisjunctive::AddNoOverlap( template bool CombinedDisjunctive::Propagate() { - helper_->SynchronizeAndSetTimeDirection(time_direction); + if (!helper_->SynchronizeAndSetTimeDirection(time_direction)) return false; const auto& task_by_increasing_end_min = helper_->TaskByIncreasingEndMin(); const auto& task_by_decreasing_start_max = helper_->TaskByDecreasingStartMax(); @@ -458,7 +458,8 @@ bool CombinedDisjunctive::Propagate() { } bool DisjunctiveOverloadChecker::Propagate() { - helper_->SynchronizeAndSetTimeDirection(/*is_forward=*/true); + if (!helper_->SynchronizeAndSetTimeDirection(/*is_forward=*/true)) + return false; // Split problem into independent part. // @@ -631,14 +632,14 @@ bool DisjunctiveOverloadChecker::PropagateSubwindow( int DisjunctiveOverloadChecker::RegisterWith(GenericLiteralWatcher* watcher) { // This propagator reach the fix point in one pass. const int id = watcher->Register(this); - helper_->SynchronizeAndSetTimeDirection(/*is_forward=*/true); + helper_->SetTimeDirection(/*is_forward=*/true); helper_->WatchAllTasks(id, watcher, /*watch_start_max=*/false, /*watch_end_max=*/true); return id; } bool DisjunctiveDetectablePrecedences::Propagate() { - helper_->SynchronizeAndSetTimeDirection(time_direction_); + if (!helper_->SynchronizeAndSetTimeDirection(time_direction_)) return false; to_propagate_.clear(); processed_.assign(helper_->NumTasks(), false); @@ -727,6 +728,7 @@ bool DisjunctiveDetectablePrecedences::PropagateSubwindow() { // TODO(user): Maybe it is just faster to merge ComputeEndMin() with // AddEntry(). task_set_.Clear(); + to_propagate_.clear(); bool need_update = false; IntegerValue task_set_end_min = kMinIntegerValue; @@ -862,7 +864,7 @@ bool DisjunctiveDetectablePrecedences::PropagateSubwindow() { int DisjunctiveDetectablePrecedences::RegisterWith( GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - helper_->SynchronizeAndSetTimeDirection(time_direction_); + helper_->SetTimeDirection(time_direction_); helper_->WatchAllTasks(id, watcher, /*watch_start_max=*/true, /*watch_end_max=*/false); watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); @@ -870,7 +872,7 @@ int DisjunctiveDetectablePrecedences::RegisterWith( } bool DisjunctivePrecedences::Propagate() { - helper_->SynchronizeAndSetTimeDirection(time_direction_); + if (!helper_->SynchronizeAndSetTimeDirection(time_direction_)) return false; window_.clear(); IntegerValue window_end = kMinIntegerValue; for (const TaskTime task_time : helper_->TaskByIncreasingShiftedStartMin()) { @@ -983,14 +985,14 @@ bool DisjunctivePrecedences::PropagateSubwindow() { int DisjunctivePrecedences::RegisterWith(GenericLiteralWatcher* watcher) { // This propagator reach the fixed point in one go. const int id = watcher->Register(this); - helper_->SynchronizeAndSetTimeDirection(time_direction_); + helper_->SetTimeDirection(time_direction_); helper_->WatchAllTasks(id, watcher, /*watch_start_max=*/false, /*watch_end_max=*/false); return id; } bool DisjunctiveNotLast::Propagate() { - helper_->SynchronizeAndSetTimeDirection(time_direction_); + if (!helper_->SynchronizeAndSetTimeDirection(time_direction_)) return false; const auto& task_by_decreasing_start_max = helper_->TaskByDecreasingStartMax(); @@ -1171,7 +1173,7 @@ int DisjunctiveNotLast::RegisterWith(GenericLiteralWatcher* watcher) { bool DisjunctiveEdgeFinding::Propagate() { const int num_tasks = helper_->NumTasks(); - helper_->SynchronizeAndSetTimeDirection(time_direction_); + if (!helper_->SynchronizeAndSetTimeDirection(time_direction_)) return false; is_gray_.resize(num_tasks, false); non_gray_task_to_event_.resize(num_tasks); @@ -1385,7 +1387,7 @@ bool DisjunctiveEdgeFinding::PropagateSubwindow(IntegerValue window_end_min) { int DisjunctiveEdgeFinding::RegisterWith(GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - helper_->SynchronizeAndSetTimeDirection(time_direction_); + helper_->SetTimeDirection(time_direction_); helper_->WatchAllTasks(id, watcher, /*watch_start_max=*/false, /*watch_end_max=*/true); watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); diff --git a/ortools/sat/feasibility_pump.h b/ortools/sat/feasibility_pump.h index ec7a4adc8b..90f91d1999 100644 --- a/ortools/sat/feasibility_pump.h +++ b/ortools/sat/feasibility_pump.h @@ -21,7 +21,7 @@ #include "ortools/lp_data/lp_data.h" #include "ortools/lp_data/lp_data_utils.h" #include "ortools/lp_data/lp_types.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_constraint.h" #include "ortools/sat/sat_solver.h" diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index e629f757cc..8f72a355ba 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -484,7 +484,7 @@ class IntegerEncoder { // Arguments: // - map is just encoding_by_var_[associated_lit.var] and is passed as a // slight optimization. - // - 'it' is the current position of associated_lit in map, i.e we must have + // - 'it' is the current position of associated_lit in map, i.e. we must have // it->second == associated_lit. void AddImplications(const std::map& map, std::map::const_iterator it, @@ -941,7 +941,7 @@ class IntegerTrail : public SatPropagator { // The current bound on this variable. IntegerValue current_bound; - // Trail index of the last TrailEntry in the trail refering to this var. + // Trail index of the last TrailEntry in the trail referring to this var. int current_trail_index; }; absl::StrongVector vars_; diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index b6614de0dc..8d98a03099 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -20,7 +20,7 @@ #include "absl/container/flat_hash_set.h" #include "absl/types/span.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/probing.h" diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc index 51ad0ec5f3..6413b1f4b8 100644 --- a/ortools/sat/intervals.cc +++ b/ortools/sat/intervals.cc @@ -90,7 +90,9 @@ SchedulingConstraintHelper::SchedulingConstraintHelper( RegisterWith(model->GetOrCreate()); InitSortedVectors(); - SynchronizeAndSetTimeDirection(true); + if (!SynchronizeAndSetTimeDirection(true)) { + model->GetOrCreate()->NotifyThatModelIsUnsat(); + } } SchedulingConstraintHelper::SchedulingConstraintHelper(int num_tasks, @@ -141,24 +143,50 @@ void SchedulingConstraintHelper::RegisterWith(GenericLiteralWatcher* watcher) { integer_trail_->RegisterReversibleClass(this); } -void SchedulingConstraintHelper::UpdateCachedValues(int t) { +bool SchedulingConstraintHelper::UpdateCachedValues(int t) { recompute_cache_[t] = false; + if (IsAbsent(t)) return true; - const IntegerValue dmin = integer_trail_->LowerBound(sizes_[t]); - const IntegerValue smin = integer_trail_->LowerBound(starts_[t]); - const IntegerValue smax = integer_trail_->UpperBound(starts_[t]); - const IntegerValue emin = integer_trail_->LowerBound(ends_[t]); - const IntegerValue emax = integer_trail_->UpperBound(ends_[t]); + IntegerValue smin = integer_trail_->LowerBound(starts_[t]); + IntegerValue smax = integer_trail_->UpperBound(starts_[t]); + IntegerValue dmin = integer_trail_->LowerBound(sizes_[t]); + IntegerValue dmax = integer_trail_->UpperBound(sizes_[t]); + IntegerValue emin = integer_trail_->LowerBound(ends_[t]); + IntegerValue emax = integer_trail_->UpperBound(ends_[t]); - cached_duration_min_[t] = dmin; - cached_start_min_[t] = smin; - cached_negated_end_max_[t] = -emax; + // Detect first if we have a conflict using the relation start + size = end. + if (smin + dmin - emax > 0) { + ClearReason(); + AddStartMinReason(t, smin); + AddSizeMinReason(t, dmin); + AddEndMaxReason(t, emax); + return PushTaskAbsence(t); + } + if (smax + dmax - emin < 0) { + ClearReason(); + AddStartMaxReason(t, smax); + AddSizeMaxReason(t, dmax); + AddEndMinReason(t, emin); + return PushTaskAbsence(t); + } - // Sometimes, for optional interval with non-optional bounds, the two - // part of the max here is not the same. We always consider the value assuming + // Sometimes, for optional interval with non-optional bounds, this propagation + // give tighter bounds. We always consider the value assuming // the interval is present. - cached_end_min_[t] = std::max(emin, smin + dmin); - cached_negated_start_max_[t] = -std::min(smax, emax - dmin); + // + // Note that this is also useful in case not everything was propagated. Note + // also that since there is no conflict, we reach the fix point in one pass. + smin = std::max(smin, emin - dmax); + smax = std::min(smax, emax - dmin); + dmin = std::max(dmin, emin - smax); + emin = std::max(emin, smin + dmin); + emax = std::min(emax, smax + dmax); + + cached_start_min_[t] = smin; + cached_end_min_[t] = emin; + cached_negated_start_max_[t] = -smax; + cached_negated_end_max_[t] = -emax; + cached_size_min_[t] = dmin; // Note that we use the cached value here for EndMin()/StartMax(). const IntegerValue new_shifted_start_min = EndMin(t) - dmin; @@ -171,9 +199,10 @@ void SchedulingConstraintHelper::UpdateCachedValues(int t) { recompute_negated_shifted_end_max_ = true; cached_negated_shifted_end_max_[t] = new_negated_shifted_end_max; } + return true; } -void SchedulingConstraintHelper::ResetFromSubset( +bool SchedulingConstraintHelper::ResetFromSubset( const SchedulingConstraintHelper& other, absl::Span tasks) { current_time_direction_ = other.current_time_direction_; @@ -195,7 +224,7 @@ void SchedulingConstraintHelper::ResetFromSubset( } InitSortedVectors(); - SynchronizeAndSetTimeDirection(true); + return SynchronizeAndSetTimeDirection(true); } void SchedulingConstraintHelper::InitSortedVectors() { @@ -206,7 +235,7 @@ void SchedulingConstraintHelper::InitSortedVectors() { cached_shifted_start_min_.resize(num_tasks); cached_negated_shifted_end_max_.resize(num_tasks); - cached_duration_min_.resize(num_tasks); + cached_size_min_.resize(num_tasks); cached_start_min_.resize(num_tasks); cached_end_min_.resize(num_tasks); cached_negated_start_max_.resize(num_tasks); @@ -231,8 +260,7 @@ void SchedulingConstraintHelper::InitSortedVectors() { recompute_negated_shifted_end_max_ = true; } -void SchedulingConstraintHelper::SynchronizeAndSetTimeDirection( - bool is_forward) { +void SchedulingConstraintHelper::SetTimeDirection(bool is_forward) { if (current_time_direction_ != is_forward) { current_time_direction_ = is_forward; @@ -249,16 +277,23 @@ void SchedulingConstraintHelper::SynchronizeAndSetTimeDirection( std::swap(cached_shifted_start_min_, cached_negated_shifted_end_max_); std::swap(recompute_shifted_start_min_, recompute_negated_shifted_end_max_); } +} + +bool SchedulingConstraintHelper::SynchronizeAndSetTimeDirection( + bool is_forward) { + SetTimeDirection(is_forward); if (recompute_all_cache_) { for (int t = 0; t < recompute_cache_.size(); ++t) { - UpdateCachedValues(t); + if (!UpdateCachedValues(t)) return false; } } else { for (int t = 0; t < recompute_cache_.size(); ++t) { - if (recompute_cache_[t]) UpdateCachedValues(t); + if (recompute_cache_[t]) + if (!UpdateCachedValues(t)) return false; } } recompute_all_cache_ = false; + return true; } const std::vector& @@ -398,14 +433,14 @@ bool SchedulingConstraintHelper::PushIntervalBound(int t, IntegerLiteral lit) { if (!PushIntegerLiteralIfTaskPresent(t, lit)) return false; if (IsAbsent(t)) return true; if (!precedences_->PropagateOutgoingArcs(lit.var)) return false; - UpdateCachedValues(t); + if (!UpdateCachedValues(t)) return false; return true; } bool SchedulingConstraintHelper::IncreaseStartMin(int t, IntegerValue new_start_min) { if (starts_[t].var == kNoIntegerVariable) { - if (new_start_min > starts_[t].constant) return ReportConflict(); + if (new_start_min > starts_[t].constant) return PushTaskAbsence(t); return true; } return PushIntervalBound(t, starts_[t].GreaterOrEqual(new_start_min)); @@ -414,15 +449,15 @@ bool SchedulingConstraintHelper::IncreaseStartMin(int t, bool SchedulingConstraintHelper::DecreaseEndMax(int t, IntegerValue new_end_max) { if (ends_[t].var == kNoIntegerVariable) { - if (new_end_max < ends_[t].constant) return ReportConflict(); + if (new_end_max < ends_[t].constant) return PushTaskAbsence(t); return true; } return PushIntervalBound(t, ends_[t].LowerOrEqual(new_end_max)); } bool SchedulingConstraintHelper::PushTaskAbsence(int t) { - DCHECK_NE(reason_for_presence_[t], kNoLiteralIndex); - DCHECK(!IsAbsent(t)); + if (IsAbsent(t)) return true; + if (!IsOptional(t)) return ReportConflict(); AddOtherReason(t); diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index 0456557369..d454c5c802 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -196,8 +196,8 @@ class SchedulingConstraintHelper : public PropagatorInterface, // Resets the class to the same state as if it was constructed with // the given subset of tasks from other. - void ResetFromSubset(const SchedulingConstraintHelper& other, - absl::Span tasks); + ABSL_MUST_USE_RESULT bool ResetFromSubset( + const SchedulingConstraintHelper& other, absl::Span tasks); // Returns the number of task. int NumTasks() const { return starts_.size(); } @@ -206,7 +206,8 @@ class SchedulingConstraintHelper : public PropagatorInterface, // either forward/backward. This will impact all the functions below. This // MUST be called at the beginning of all Propagate() call that uses this // helper. - void SynchronizeAndSetTimeDirection(bool is_forward); + void SetTimeDirection(bool is_forward); + ABSL_MUST_USE_RESULT bool SynchronizeAndSetTimeDirection(bool is_forward); // Helpers for the current bounds on the current task time window. // [ (size-min) ... (size-min) ] @@ -221,7 +222,7 @@ class SchedulingConstraintHelper : public PropagatorInterface, // cases where pushing the start of one task will change values for many // others. This is fine as the new values will be picked up as we reach the // propagation fixed point. - IntegerValue SizeMin(int t) const { return cached_duration_min_[t]; } + IntegerValue SizeMin(int t) const { return cached_size_min_[t]; } IntegerValue SizeMax(int t) const { // This one is "rare" so we don't cache it. return integer_trail_->UpperBound(sizes_[t]); @@ -287,11 +288,14 @@ class SchedulingConstraintHelper : public PropagatorInterface, void AddAbsenceReason(int t); void AddSizeMinReason(int t); void AddSizeMinReason(int t, IntegerValue lower_bound); + void AddSizeMaxReason(int t, IntegerValue upper_bound); void AddStartMinReason(int t, IntegerValue lower_bound); void AddStartMaxReason(int t, IntegerValue upper_bound); void AddEndMinReason(int t, IntegerValue lower_bound); void AddEndMaxReason(int t, IntegerValue upper_bound); + void AddEnergyAfterReason(int t, IntegerValue energy_min, IntegerValue time); + void AddEnergyMinInIntervalReason(int t, IntegerValue min, IntegerValue max); // Adds the reason why task "before" must be before task "after". // That is StartMax(before) < EndMin(after). @@ -363,8 +367,13 @@ class SchedulingConstraintHelper : public PropagatorInterface, bool InPropagationLoop() const { return integer_trail_->InPropagationLoop(); } private: + // Generic reason for a <= upper_bound, given that a = b + c in case the + // current upper bound of a is not good enough. + void AddGenericReason(const AffineExpression& a, IntegerValue upper_bound, + const AffineExpression& b, const AffineExpression& c); + void InitSortedVectors(); - void UpdateCachedValues(int t); + ABSL_MUST_USE_RESULT bool UpdateCachedValues(int t); // Internal function for IncreaseStartMin()/DecreaseEndMax(). bool PushIntervalBound(int t, IntegerLiteral lit); @@ -401,7 +410,7 @@ class SchedulingConstraintHelper : public PropagatorInterface, int previous_level_ = 0; // The caches of all relevant interval values. - std::vector cached_duration_min_; + std::vector cached_size_min_; std::vector cached_start_min_; std::vector cached_end_min_; std::vector cached_negated_start_max_; @@ -493,98 +502,106 @@ inline void SchedulingConstraintHelper::AddAbsenceReason(int t) { } inline void SchedulingConstraintHelper::AddSizeMinReason(int t) { - AddOtherReason(t); - if (sizes_[t].var != kNoIntegerVariable) { - integer_reason_.push_back( - integer_trail_->LowerBoundAsLiteral(sizes_[t].var)); + AddSizeMinReason(t, SizeMin(t)); +} + +inline void SchedulingConstraintHelper::AddGenericReason( + const AffineExpression& a, IntegerValue upper_bound, + const AffineExpression& b, const AffineExpression& c) { + if (integer_trail_->UpperBound(a) <= upper_bound) { + if (a.var != kNoIntegerVariable) { + integer_reason_.push_back(a.LowerOrEqual(upper_bound)); + } + return; + } + CHECK_NE(a.var, kNoIntegerVariable); + + // Here we assume that the upper_bound on a comes from the lower bound of b + + // c. + const IntegerValue slack = upper_bound - integer_trail_->UpperBound(b) - + integer_trail_->UpperBound(c); + CHECK_GE(slack, 0); + if (b.var == kNoIntegerVariable && c.var == kNoIntegerVariable) return; + if (b.var == kNoIntegerVariable) { + integer_reason_.push_back(c.LowerOrEqual(upper_bound - b.constant)); + } else if (c.var == kNoIntegerVariable) { + integer_reason_.push_back(b.LowerOrEqual(upper_bound - c.constant)); + } else { + integer_trail_->AppendRelaxedLinearReason( + slack, {IntegerValue(1), IntegerValue(1)}, + {NegationOf(b.var), NegationOf(c.var)}, &integer_reason_); } } inline void SchedulingConstraintHelper::AddSizeMinReason( int t, IntegerValue lower_bound) { AddOtherReason(t); - if (sizes_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(sizes_[t].GreaterOrEqual(lower_bound)); - } + DCHECK(!IsAbsent(t)); + AddGenericReason(sizes_[t].Negated(), -lower_bound, minus_ends_[t], + starts_[t]); +} + +inline void SchedulingConstraintHelper::AddSizeMaxReason( + int t, IntegerValue upper_bound) { + AddOtherReason(t); + CHECK(!IsAbsent(t)); + AddGenericReason(sizes_[t], upper_bound, ends_[t], minus_starts_[t]); } inline void SchedulingConstraintHelper::AddStartMinReason( int t, IntegerValue lower_bound) { - DCHECK_GE(StartMin(t), lower_bound); AddOtherReason(t); - if (starts_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(starts_[t].GreaterOrEqual(lower_bound)); - } + DCHECK(!IsAbsent(t)); + AddGenericReason(minus_starts_[t], -lower_bound, minus_ends_[t], sizes_[t]); } inline void SchedulingConstraintHelper::AddStartMaxReason( int t, IntegerValue upper_bound) { AddOtherReason(t); - - // Note that we cannot use the cache here! - if (integer_trail_->UpperBound(starts_[t]) <= upper_bound) { - if (starts_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(starts_[t].LowerOrEqual(upper_bound)); - } - } else { - // This might happen if we used StartMax() <= EndMax() - SizeMin(). - if (sizes_[t].var != kNoIntegerVariable) { - integer_reason_.push_back( - integer_trail_->LowerBoundAsLiteral(sizes_[t].var)); - } - if (ends_[t].var != kNoIntegerVariable) { - integer_reason_.push_back( - ends_[t].LowerOrEqual(upper_bound + SizeMin(t))); - } - } + DCHECK(!IsAbsent(t)); + AddGenericReason(starts_[t], upper_bound, ends_[t], sizes_[t].Negated()); } inline void SchedulingConstraintHelper::AddEndMinReason( int t, IntegerValue lower_bound) { AddOtherReason(t); - - // Note that we cannot use the cache here! - if (integer_trail_->LowerBound(ends_[t]) >= lower_bound) { - if (ends_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(ends_[t].GreaterOrEqual(lower_bound)); - } - } else { - // This might happen if we used EndMin() >= StartMin() + SizeMin(). - if (sizes_[t].var != kNoIntegerVariable) { - integer_reason_.push_back( - integer_trail_->LowerBoundAsLiteral(sizes_[t].var)); - } - if (starts_[t].var != kNoIntegerVariable) { - integer_reason_.push_back( - starts_[t].GreaterOrEqual(lower_bound - SizeMin(t))); - } - } + DCHECK(!IsAbsent(t)); + AddGenericReason(minus_ends_[t], -lower_bound, minus_starts_[t], + sizes_[t].Negated()); } inline void SchedulingConstraintHelper::AddEndMaxReason( int t, IntegerValue upper_bound) { - DCHECK_LE(EndMax(t), upper_bound); AddOtherReason(t); - if (ends_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(ends_[t].LowerOrEqual(upper_bound)); - } + DCHECK(!IsAbsent(t)); + AddGenericReason(ends_[t], upper_bound, starts_[t], sizes_[t]); } inline void SchedulingConstraintHelper::AddEnergyAfterReason( int t, IntegerValue energy_min, IntegerValue time) { - AddOtherReason(t); if (StartMin(t) >= time) { - if (starts_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(starts_[t].GreaterOrEqual(time)); - } + AddStartMinReason(t, time); } else { - if (ends_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(ends_[t].GreaterOrEqual(time + energy_min)); - } + AddEndMinReason(t, time + energy_min); } - if (sizes_[t].var != kNoIntegerVariable) { - integer_reason_.push_back(sizes_[t].GreaterOrEqual(energy_min)); + AddSizeMinReason(t, energy_min); +} + +inline void SchedulingConstraintHelper::AddEnergyMinInIntervalReason( + int t, IntegerValue time_min, IntegerValue time_max) { + const IntegerValue energy_min = SizeMin(t); + CHECK_LE(time_min + energy_min, time_max); + if (StartMin(t) >= time_min) { + AddStartMinReason(t, time_min); + } else { + AddEndMinReason(t, time_min + energy_min); } + if (EndMax(t) <= time_max) { + AddEndMaxReason(t, time_max); + } else { + AddStartMaxReason(t, time_max - energy_min); + } + AddSizeMinReason(t, energy_min); } // ============================================================================= diff --git a/ortools/sat/lb_tree_search.cc b/ortools/sat/lb_tree_search.cc index 56c8222aa1..a8dae02ca8 100644 --- a/ortools/sat/lb_tree_search.cc +++ b/ortools/sat/lb_tree_search.cc @@ -13,7 +13,9 @@ #include "ortools/sat/lb_tree_search.h" -#include "ortools/sat/cp_model_loader.h" +#include + +#include "ortools/sat/cp_model_mapping.h" namespace operations_research { namespace sat { diff --git a/ortools/sat/lb_tree_search.h b/ortools/sat/lb_tree_search.h index b3e98b2106..7081bb7101 100644 --- a/ortools/sat/lb_tree_search.h +++ b/ortools/sat/lb_tree_search.h @@ -14,6 +14,7 @@ #ifndef OR_TOOLS_SAT_LB_TREE_SEARCH_H_ #define OR_TOOLS_SAT_LB_TREE_SEARCH_H_ +#include #include #include "ortools/sat/integer.h" @@ -81,8 +82,8 @@ class LbTreeSearch { IntegerValue false_objective; // Points to adjacent nodes in the tree. Large if no connection. - int true_child = kint32max; - int false_child = kint32max; + int true_child = std::numeric_limits::max(); + int false_child = std::numeric_limits::max(); // Instead of storing the full reason for an objective LB increase in one // the branches (which can lead to a quadratic memory usage), we stores the @@ -95,8 +96,8 @@ class LbTreeSearch { // and skip all the nodes in-between by connecting directly the correct // ancestor to this node. Note that when we do that, the level of the nodes // in the sub-branch change, but this still work. - int true_level = kint32max; - int false_level = kint32max; + int true_level = std::numeric_limits::max(); + int false_level = std::numeric_limits::max(); }; // Returns true if this node objective lb is greater than the root level diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 6122e23273..280a54061b 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -1436,7 +1436,10 @@ bool LinearProgrammingConstraint::Propagate() { (trail_->CurrentDecisionLevel() == 0 || !sat_parameters_.only_add_cuts_at_level_zero())) { for (const CutGenerator& generator : cut_generators_) { - generator.generate_cuts(expanded_lp_solution_, &constraint_manager_); + if (!generator.generate_cuts(expanded_lp_solution_, + &constraint_manager_)) { + return false; + } } } @@ -1484,32 +1487,34 @@ bool LinearProgrammingConstraint::Propagate() { if (objective_is_defined_ && (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL || simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_FEASIBLE)) { - // Try to filter optimal objective value. Note that GetObjectiveValue() - // already take care of the scaling so that it returns an objective in the - // CP world. - const double relaxed_optimal_objective = simplex_.GetObjectiveValue(); - const IntegerValue approximate_new_lb(static_cast( - std::ceil(relaxed_optimal_objective - kCpEpsilon))); - // TODO(user): Maybe do a bit less computation when we cannot propagate // anything. if (sat_parameters_.use_exact_lp_reason()) { if (!ExactLpReasonning()) return false; // Display when the inexact bound would have propagated more. - const IntegerValue propagated_lb = - integer_trail_->LowerBound(objective_cp_); - if (approximate_new_lb > propagated_lb) { - VLOG(2) << "LP objective [ " << ToDouble(propagated_lb) << ", " - << ToDouble(integer_trail_->UpperBound(objective_cp_)) - << " ] approx_lb += " - << ToDouble(approximate_new_lb - propagated_lb) << " gap: " - << integer_trail_->UpperBound(objective_cp_) - propagated_lb; + if (VLOG_IS_ON(2)) { + const double relaxed_optimal_objective = simplex_.GetObjectiveValue(); + const IntegerValue approximate_new_lb(static_cast( + std::ceil(relaxed_optimal_objective - kCpEpsilon))); + const IntegerValue propagated_lb = + integer_trail_->LowerBound(objective_cp_); + if (approximate_new_lb > propagated_lb) { + VLOG(2) << "LP objective [ " << ToDouble(propagated_lb) << ", " + << ToDouble(integer_trail_->UpperBound(objective_cp_)) + << " ] approx_lb += " + << ToDouble(approximate_new_lb - propagated_lb) << " gap: " + << integer_trail_->UpperBound(objective_cp_) - propagated_lb; + } } } else { + // Try to filter optimal objective value. Note that GetObjectiveValue() + // already take care of the scaling so that it returns an objective in the + // CP world. FillReducedCostReasonIn(simplex_.GetReducedCosts(), &integer_reason_); const double objective_cp_ub = ToDouble(integer_trail_->UpperBound(objective_cp_)); + const double relaxed_optimal_objective = simplex_.GetObjectiveValue(); ReducedCostStrengtheningDeductions(objective_cp_ub - relaxed_optimal_objective); if (!deductions_.empty()) { @@ -1519,6 +1524,8 @@ bool LinearProgrammingConstraint::Propagate() { } // Push new objective lb. + const IntegerValue approximate_new_lb(static_cast( + std::ceil(relaxed_optimal_objective - kCpEpsilon))); if (approximate_new_lb > integer_trail_->LowerBound(objective_cp_)) { const IntegerLiteral deduction = IntegerLiteral::GreaterOrEqual(objective_cp_, approximate_new_lb); @@ -2547,6 +2554,7 @@ CutGenerator CreateStronglyConnectedGraphCutGenerator( SeparateSubtourInequalities( num_nodes, tails, heads, literals, lp_values, /*demands=*/{}, /*capacity=*/0, manager, model); + return true; }; return result; } @@ -2566,6 +2574,7 @@ CutGenerator CreateCVRPCutGenerator(int num_nodes, SeparateSubtourInequalities(num_nodes, tails, heads, literals, lp_values, demands, capacity, manager, model); + return true; }; return result; } diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index d8159b3169..e5c62dd438 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -20,8 +20,9 @@ #include "absl/container/flat_hash_set.h" #include "ortools/base/iterator_adaptors.h" #include "ortools/base/stl_util.h" +#include "ortools/sat/circuit.h" // for ReindexArcs. #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_expr.h" #include "ortools/sat/linear_constraint.h" @@ -314,208 +315,265 @@ bool AllLiteralsHaveViews(const IntegerEncoder& encoder, } // namespace -// Add a linear relaxation of the CP constraint to the set of linear -// constraints. The highest linearization_level is, the more types of constraint -// we encode. This method should be called only for linearization_level > 0. -// -// Note: IntProd is linearized dynamically using the cut generators. -// -// TODO(user): In full generality, we could encode all the constraint as an LP. -// TODO(user,user): Add unit tests for this method. -void TryToLinearizeConstraint(const CpModelProto& model_proto, - const ConstraintProto& ct, Model* model, - int linearization_level, - LinearRelaxation* relaxation) { - CHECK_EQ(model->GetOrCreate()->CurrentDecisionLevel(), 0); - DCHECK_GT(linearization_level, 0); +void AppendBoolOrRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { auto* mapping = model->GetOrCreate(); - const auto& encoder = *model->GetOrCreate(); - if (ct.constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) { - if (linearization_level < 2) return; - LinearConstraintBuilder lc(model, IntegerValue(1), kMaxIntegerValue); - for (const int enforcement_ref : ct.enforcement_literal()) { - CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)), - IntegerValue(1))); - } - for (const int ref : ct.bool_or().literals()) { - CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1))); - } - relaxation->linear_constraints.push_back(lc.Build()); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kBoolAnd) { - // TODO(user): These constraints can be many, and if they are not regrouped - // in big at most ones, then they should probably only added lazily as cuts. - // Regroup this with future clique-cut separation logic. - if (linearization_level < 2) return; - if (!HasEnforcementLiteral(ct)) return; - if (ct.enforcement_literal().size() == 1) { - const Literal enforcement = mapping->Literal(ct.enforcement_literal(0)); - for (const int ref : ct.bool_and().literals()) { - relaxation->at_most_ones.push_back( - {enforcement, mapping->Literal(ref).Negated()}); - } - return; - } + LinearConstraintBuilder lc(model, IntegerValue(1), kMaxIntegerValue); + for (const int enforcement_ref : ct.enforcement_literal()) { + CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)), + IntegerValue(1))); + } + for (const int ref : ct.bool_or().literals()) { + CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1))); + } + relaxation->linear_constraints.push_back(lc.Build()); +} - // Andi(e_i) => Andj(x_j) - // <=> num_rhs_terms <= Sum_j(x_j) + num_rhs_terms * Sum_i(~e_i) - int num_literals = ct.bool_and().literals_size(); - LinearConstraintBuilder lc(model, IntegerValue(num_literals), - kMaxIntegerValue); +void AppendBoolAndRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { + // TODO(user): These constraints can be many, and if they are not regrouped + // in big at most ones, then they should probably only added lazily as cuts. + // Regroup this with future clique-cut separation logic. + if (!HasEnforcementLiteral(ct)) return; + + auto* mapping = model->GetOrCreate(); + if (ct.enforcement_literal().size() == 1) { + const Literal enforcement = mapping->Literal(ct.enforcement_literal(0)); for (const int ref : ct.bool_and().literals()) { - CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1))); + relaxation->at_most_ones.push_back( + {enforcement, mapping->Literal(ref).Negated()}); } - for (const int enforcement_ref : ct.enforcement_literal()) { - CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)), - IntegerValue(num_literals))); + return; + } + + // Andi(e_i) => Andj(x_j) + // <=> num_rhs_terms <= Sum_j(x_j) + num_rhs_terms * Sum_i(~e_i) + int num_literals = ct.bool_and().literals_size(); + LinearConstraintBuilder lc(model, IntegerValue(num_literals), + kMaxIntegerValue); + for (const int ref : ct.bool_and().literals()) { + CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1))); + } + for (const int enforcement_ref : ct.enforcement_literal()) { + CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)), + IntegerValue(num_literals))); + } + relaxation->linear_constraints.push_back(lc.Build()); +} + +void AppendAtMostOneRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + + auto* mapping = model->GetOrCreate(); + relaxation->at_most_ones.push_back( + mapping->Literals(ct.at_most_one().literals())); +} + +void AppendExactlyOneRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + auto* mapping = model->GetOrCreate(); + auto* encoder = model->GetOrCreate(); + + const std::vector literals = + mapping->Literals(ct.exactly_one().literals()); + if (AllLiteralsHaveViews(*encoder, literals)) { + LinearConstraintBuilder lc(model, IntegerValue(1), IntegerValue(1)); + for (const Literal lit : literals) { + CHECK(lc.AddLiteralTerm(lit, IntegerValue(1))); } relaxation->linear_constraints.push_back(lc.Build()); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kAtMostOne) { - if (HasEnforcementLiteral(ct)) return; - relaxation->at_most_ones.push_back( - mapping->Literals(ct.at_most_one().literals())); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kExactlyOne) { - if (HasEnforcementLiteral(ct)) return; - const std::vector literals = - mapping->Literals(ct.exactly_one().literals()); - if (linearization_level >= 2 || AllLiteralsHaveViews(encoder, literals)) { - LinearConstraintBuilder lc(model, IntegerValue(1), IntegerValue(1)); - for (const Literal lit : literals) { - CHECK(lc.AddLiteralTerm(lit, IntegerValue(1))); - } - relaxation->linear_constraints.push_back(lc.Build()); - } else if (linearization_level == 1) { - // We just encode the at most one part that might be partially linearized - // later. - relaxation->at_most_ones.push_back(literals); - } - } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMax) { - if (HasEnforcementLiteral(ct)) return; - const IntegerVariable target = mapping->Integer(ct.int_max().target()); - const std::vector vars = - mapping->Integers(ct.int_max().vars()); - AppendMaxRelaxation(target, vars, linearization_level, model, relaxation); + } else { + // We just encode the at most one part that might be partially linearized + // later. + relaxation->at_most_ones.push_back(literals); + } +} - } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMin) { - if (HasEnforcementLiteral(ct)) return; - const IntegerVariable negative_target = - NegationOf(mapping->Integer(ct.int_min().target())); - const std::vector negative_vars = - NegationOf(mapping->Integers(ct.int_min().vars())); - AppendMaxRelaxation(negative_target, negative_vars, linearization_level, - model, relaxation); - } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear) { - AppendLinearConstraintRelaxation(ct, linearization_level, *model, - relaxation); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kCircuit) { - if (HasEnforcementLiteral(ct)) return; - const int num_arcs = ct.circuit().literals_size(); - CHECK_EQ(num_arcs, ct.circuit().tails_size()); - CHECK_EQ(num_arcs, ct.circuit().heads_size()); +namespace { +// Adds linearization of int max constraints. This can also be used to linearize +// int min with negated variables. +void AppendMaxRelaxationHelper(IntegerVariable target, + const std::vector& vars, + bool encode_other_direction, Model* model, + LinearRelaxation* relaxation) { + // Case X = max(X_1, X_2, ..., X_N) + // Part 1: Encode X >= max(X_1, X_2, ..., X_N) + for (const IntegerVariable var : vars) { + // This deal with the corner case X = max(X, Y, Z, ..) ! + // Note that this can be presolved into X >= Y, X >= Z, ... + if (target == var) continue; + LinearConstraintBuilder lc(model, kMinIntegerValue, IntegerValue(0)); + lc.AddTerm(var, IntegerValue(1)); + lc.AddTerm(target, IntegerValue(-1)); + relaxation->linear_constraints.push_back(lc.Build()); + } - // Each node must have exactly one incoming and one outgoing arc (note that - // it can be the unique self-arc of this node too). - std::map> incoming_arc_constraints; - std::map> outgoing_arc_constraints; - for (int i = 0; i < num_arcs; i++) { - const Literal arc = mapping->Literal(ct.circuit().literals(i)); - const int tail = ct.circuit().tails(i); - const int head = ct.circuit().heads(i); + // Part 2: Encode upper bound on X. + if (!encode_other_direction) return; + GenericLiteralWatcher* watcher = model->GetOrCreate(); + // For size = 2, we do this with 1 less variable. + IntegerEncoder* encoder = model->GetOrCreate(); + if (vars.size() == 2) { + IntegerVariable y = model->Add(NewIntegerVariable(0, 1)); + const Literal y_lit = + encoder->GetOrCreateLiteralAssociatedToEquality(y, IntegerValue(1)); + AppendEnforcedUpperBound(y_lit, target, vars[0], model, relaxation); - // Make sure this literal has a view. - model->Add(NewIntegerVariableFromLiteral(arc)); - outgoing_arc_constraints[tail].push_back(arc); - incoming_arc_constraints[head].push_back(arc); - } - for (const auto* node_map : - {&outgoing_arc_constraints, &incoming_arc_constraints}) { - for (const auto& entry : *node_map) { - const std::vector& exactly_one = entry.second; - if (exactly_one.size() > 1) { - LinearConstraintBuilder at_least_one_lc(model, IntegerValue(1), - kMaxIntegerValue); - for (const Literal l : exactly_one) { - CHECK(at_least_one_lc.AddLiteralTerm(l, IntegerValue(1))); - } + // TODO(user,user): It makes more sense to use ConditionalLowerOrEqual() + // here, but that degrades perf on the road*.fzn problem. Understand why. + IntegerSumLE* upper_bound1 = new IntegerSumLE( + {y_lit}, {target, vars[0]}, {IntegerValue(1), IntegerValue(-1)}, + IntegerValue(0), model); + upper_bound1->RegisterWith(watcher); + model->TakeOwnership(upper_bound1); + AppendEnforcedUpperBound(y_lit.Negated(), target, vars[1], model, + relaxation); + IntegerSumLE* upper_bound2 = new IntegerSumLE( + {y_lit.Negated()}, {target, vars[1]}, + {IntegerValue(1), IntegerValue(-1)}, IntegerValue(0), model); + upper_bound2->RegisterWith(watcher); + model->TakeOwnership(upper_bound2); + return; + } + // For each X_i, we encode y_i => X <= X_i. And at least one of the y_i is + // true. Note that the correct y_i will be chosen because of the first part in + // linearlization (X >= X_i). + // TODO(user,user): Only lower bound is needed, experiment. + LinearConstraintBuilder lc_exactly_one(model, IntegerValue(1), + IntegerValue(1)); + std::vector exactly_one_literals; + exactly_one_literals.reserve(vars.size()); + for (const IntegerVariable var : vars) { + if (target == var) continue; + // y => X <= X_i. + // <=> max_term_value * y + X - X_i <= max_term_value. + // where max_tern_value is X_ub - X_i_lb. + IntegerVariable y = model->Add(NewIntegerVariable(0, 1)); + const Literal y_lit = + encoder->GetOrCreateLiteralAssociatedToEquality(y, IntegerValue(1)); - // We separate the two constraints. - relaxation->at_most_ones.push_back(exactly_one); - relaxation->linear_constraints.push_back(at_least_one_lc.Build()); + AppendEnforcedUpperBound(y_lit, target, var, model, relaxation); + IntegerSumLE* upper_bound_constraint = new IntegerSumLE( + {y_lit}, {target, var}, {IntegerValue(1), IntegerValue(-1)}, + IntegerValue(0), model); + upper_bound_constraint->RegisterWith(watcher); + model->TakeOwnership(upper_bound_constraint); + exactly_one_literals.push_back(y_lit); + + CHECK(lc_exactly_one.AddLiteralTerm(y_lit, IntegerValue(1))); + } + model->Add(ExactlyOneConstraint(exactly_one_literals)); + relaxation->linear_constraints.push_back(lc_exactly_one.Build()); +} +} // namespace + +void AppendIntMaxRelaxation(const ConstraintProto& ct, + bool encode_other_direction, Model* model, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + + auto* mapping = model->GetOrCreate(); + const IntegerVariable target = mapping->Integer(ct.int_max().target()); + const std::vector vars = + mapping->Integers(ct.int_max().vars()); + AppendMaxRelaxationHelper(target, vars, encode_other_direction, model, + relaxation); +} + +void AppendIntMinRelaxation(const ConstraintProto& ct, + bool encode_other_direction, Model* model, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + auto* mapping = model->GetOrCreate(); + const IntegerVariable negative_target = + NegationOf(mapping->Integer(ct.int_min().target())); + const std::vector negative_vars = + NegationOf(mapping->Integers(ct.int_min().vars())); + AppendMaxRelaxationHelper(negative_target, negative_vars, + encode_other_direction, model, relaxation); +} + +void AppendCircuitRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + auto* mapping = model->GetOrCreate(); + const int num_arcs = ct.circuit().literals_size(); + CHECK_EQ(num_arcs, ct.circuit().tails_size()); + CHECK_EQ(num_arcs, ct.circuit().heads_size()); + + // Each node must have exactly one incoming and one outgoing arc (note + // that it can be the unique self-arc of this node too). + std::map> incoming_arc_constraints; + std::map> outgoing_arc_constraints; + for (int i = 0; i < num_arcs; i++) { + const Literal arc = mapping->Literal(ct.circuit().literals(i)); + const int tail = ct.circuit().tails(i); + const int head = ct.circuit().heads(i); + + // Make sure this literal has a view. + model->Add(NewIntegerVariableFromLiteral(arc)); + outgoing_arc_constraints[tail].push_back(arc); + incoming_arc_constraints[head].push_back(arc); + } + for (const auto* node_map : + {&outgoing_arc_constraints, &incoming_arc_constraints}) { + for (const auto& entry : *node_map) { + const std::vector& exactly_one = entry.second; + if (exactly_one.size() > 1) { + LinearConstraintBuilder at_least_one_lc(model, IntegerValue(1), + kMaxIntegerValue); + for (const Literal l : exactly_one) { + CHECK(at_least_one_lc.AddLiteralTerm(l, IntegerValue(1))); } + + // We separate the two constraints. + relaxation->at_most_ones.push_back(exactly_one); + relaxation->linear_constraints.push_back(at_least_one_lc.Build()); } } - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kElement) { - const IntegerVariable index = mapping->Integer(ct.element().index()); - const IntegerVariable target = mapping->Integer(ct.element().target()); - const std::vector vars = - mapping->Integers(ct.element().vars()); + } +} - // We only relax the case where all the vars are constant. - // target = sum (index == i) * fixed_vars[i]. - LinearConstraintBuilder constraint(model, IntegerValue(0), IntegerValue(0)); - constraint.AddTerm(target, IntegerValue(-1)); - IntegerTrail* integer_trail = model->GetOrCreate(); - for (const auto literal_value : model->Add(FullyEncodeVariable((index)))) { - const IntegerVariable var = vars[literal_value.value.value()]; - if (!model->Get(IsFixed(var))) return; +void AppendIntervalRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation) { + // If the interval is using views, then the linear equation is already + // present in the model. + if (ct.interval().has_start_view()) return; - // Make sure this literal has a view. - model->Add(NewIntegerVariableFromLiteral(literal_value.literal)); - CHECK(constraint.AddLiteralTerm(literal_value.literal, - integer_trail->LowerBound(var))); - } - - relaxation->linear_constraints.push_back(constraint.Build()); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kInterval) { - // If the interval is using views, then the linear equation is already - // present in the model. - if (linearization_level < 2) return; - if (ct.interval().has_start_view()) return; - const IntegerVariable start = mapping->Integer(ct.interval().start()); - const IntegerVariable size = mapping->Integer(ct.interval().size()); - const IntegerVariable end = mapping->Integer(ct.interval().end()); - IntegerTrail* integer_trail = model->GetOrCreate(); - const bool size_is_fixed = integer_trail->IsFixed(size); - const IntegerValue rhs = - size_is_fixed ? -integer_trail->LowerBound(size) : IntegerValue(0); - LinearConstraintBuilder lc(model, rhs, rhs); - lc.AddTerm(start, IntegerValue(1)); - if (!size_is_fixed) { - lc.AddTerm(size, IntegerValue(1)); - } - lc.AddTerm(end, IntegerValue(-1)); - if (HasEnforcementLiteral(ct)) { - LinearConstraint tmp_lc = lc.Build(); - LinearExpression expr; - expr.coeffs = tmp_lc.coeffs; - expr.vars = tmp_lc.vars; - AppendEnforcedLinearExpression( - mapping->Literals(ct.enforcement_literal()), expr, tmp_lc.ub, - tmp_lc.ub, *model, relaxation); - } else { - relaxation->linear_constraints.push_back(lc.Build()); - } - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kNoOverlap) { - AppendNoOverlapRelaxation(model_proto, ct, linearization_level, model, - relaxation); - } else if (ct.constraint_case() == - ConstraintProto::ConstraintCase::kCumulative) { - AppendCumulativeRelaxation(model_proto, ct, linearization_level, model, - relaxation); + auto* mapping = model->GetOrCreate(); + const IntegerVariable start = mapping->Integer(ct.interval().start()); + const IntegerVariable size = mapping->Integer(ct.interval().size()); + const IntegerVariable end = mapping->Integer(ct.interval().end()); + IntegerTrail* integer_trail = model->GetOrCreate(); + const bool size_is_fixed = integer_trail->IsFixed(size); + const IntegerValue rhs = + size_is_fixed ? -integer_trail->LowerBound(size) : IntegerValue(0); + LinearConstraintBuilder lc(model, rhs, rhs); + lc.AddTerm(start, IntegerValue(1)); + if (!size_is_fixed) { + lc.AddTerm(size, IntegerValue(1)); + } + lc.AddTerm(end, IntegerValue(-1)); + if (HasEnforcementLiteral(ct)) { + LinearConstraint tmp_lc = lc.Build(); + LinearExpression expr; + expr.coeffs = tmp_lc.coeffs; + expr.vars = tmp_lc.vars; + AppendEnforcedLinearExpression(mapping->Literals(ct.enforcement_literal()), + expr, tmp_lc.ub, tmp_lc.ub, *model, + relaxation); + } else { + relaxation->linear_constraints.push_back(lc.Build()); } } // TODO(user): Use affine demand. -void AddCumulativeCut(const std::vector& intervals, - const std::vector& demands, - IntegerValue capacity_upper_bound, Model* model, - LinearRelaxation* relaxation) { +void AddCumulativeRelaxation(const std::vector& intervals, + const std::vector& demands, + IntegerValue capacity_upper_bound, Model* model, + LinearRelaxation* relaxation) { // TODO(user): Keep a map intervals -> helper, or ct_index->helper to avoid // creating many helpers for the same constraint. auto* helper = new SchedulingConstraintHelper(intervals, model); @@ -605,11 +663,9 @@ void AddCumulativeCut(const std::vector& intervals, } void AppendCumulativeRelaxation(const CpModelProto& model_proto, - const ConstraintProto& ct, - int linearization_level, Model* model, + const ConstraintProto& ct, Model* model, LinearRelaxation* relaxation) { CHECK(ct.has_cumulative()); - if (linearization_level < 2) return; if (HasEnforcementLiteral(ct)) return; auto* mapping = model->GetOrCreate(); @@ -620,96 +676,22 @@ void AppendCumulativeRelaxation(const CpModelProto& model_proto, const IntegerValue capacity_upper_bound = model->GetOrCreate()->UpperBound( mapping->Integer(ct.cumulative().capacity())); - AddCumulativeCut(intervals, demands, capacity_upper_bound, model, relaxation); + AddCumulativeRelaxation(intervals, demands, capacity_upper_bound, model, + relaxation); } void AppendNoOverlapRelaxation(const CpModelProto& model_proto, - const ConstraintProto& ct, - int linearization_level, Model* model, + const ConstraintProto& ct, Model* model, LinearRelaxation* relaxation) { CHECK(ct.has_no_overlap()); - if (linearization_level < 2) return; if (HasEnforcementLiteral(ct)) return; auto* mapping = model->GetOrCreate(); std::vector intervals = mapping->Intervals(ct.no_overlap().intervals()); - AddCumulativeCut(intervals, /*demands=*/{}, - /*capacity_upper_bound=*/IntegerValue(1), model, relaxation); -} - -void AppendMaxRelaxation(IntegerVariable target, - const std::vector& vars, - int linearization_level, Model* model, - LinearRelaxation* relaxation) { - // Case X = max(X_1, X_2, ..., X_N) - // Part 1: Encode X >= max(X_1, X_2, ..., X_N) - for (const IntegerVariable var : vars) { - // This deal with the corner case X = max(X, Y, Z, ..) ! - // Note that this can be presolved into X >= Y, X >= Z, ... - if (target == var) continue; - LinearConstraintBuilder lc(model, kMinIntegerValue, IntegerValue(0)); - lc.AddTerm(var, IntegerValue(1)); - lc.AddTerm(target, IntegerValue(-1)); - relaxation->linear_constraints.push_back(lc.Build()); - } - - // Part 2: Encode upper bound on X. - if (linearization_level < 2) return; - GenericLiteralWatcher* watcher = model->GetOrCreate(); - // For size = 2, we do this with 1 less variable. - IntegerEncoder* encoder = model->GetOrCreate(); - if (vars.size() == 2) { - IntegerVariable y = model->Add(NewIntegerVariable(0, 1)); - const Literal y_lit = - encoder->GetOrCreateLiteralAssociatedToEquality(y, IntegerValue(1)); - AppendEnforcedUpperBound(y_lit, target, vars[0], model, relaxation); - - // TODO(user,user): It makes more sense to use ConditionalLowerOrEqual() - // here, but that degrades perf on the road*.fzn problem. Understand why. - IntegerSumLE* upper_bound1 = new IntegerSumLE( - {y_lit}, {target, vars[0]}, {IntegerValue(1), IntegerValue(-1)}, - IntegerValue(0), model); - upper_bound1->RegisterWith(watcher); - model->TakeOwnership(upper_bound1); - AppendEnforcedUpperBound(y_lit.Negated(), target, vars[1], model, - relaxation); - IntegerSumLE* upper_bound2 = new IntegerSumLE( - {y_lit.Negated()}, {target, vars[1]}, - {IntegerValue(1), IntegerValue(-1)}, IntegerValue(0), model); - upper_bound2->RegisterWith(watcher); - model->TakeOwnership(upper_bound2); - return; - } - // For each X_i, we encode y_i => X <= X_i. And at least one of the y_i is - // true. Note that the correct y_i will be chosen because of the first part in - // linearlization (X >= X_i). - // TODO(user): Only lower bound is needed, experiment. - LinearConstraintBuilder lc_exactly_one(model, IntegerValue(1), - IntegerValue(1)); - std::vector exactly_one_literals; - exactly_one_literals.reserve(vars.size()); - for (const IntegerVariable var : vars) { - if (target == var) continue; - // y => X <= X_i. - // <=> max_term_value * y + X - X_i <= max_term_value. - // where max_tern_value is X_ub - X_i_lb. - IntegerVariable y = model->Add(NewIntegerVariable(0, 1)); - const Literal y_lit = - encoder->GetOrCreateLiteralAssociatedToEquality(y, IntegerValue(1)); - - AppendEnforcedUpperBound(y_lit, target, var, model, relaxation); - IntegerSumLE* upper_bound_constraint = new IntegerSumLE( - {y_lit}, {target, var}, {IntegerValue(1), IntegerValue(-1)}, - IntegerValue(0), model); - upper_bound_constraint->RegisterWith(watcher); - model->TakeOwnership(upper_bound_constraint); - exactly_one_literals.push_back(y_lit); - - CHECK(lc_exactly_one.AddLiteralTerm(y_lit, IntegerValue(1))); - } - model->Add(ExactlyOneConstraint(exactly_one_literals)); - relaxation->linear_constraints.push_back(lc_exactly_one.Build()); + AddCumulativeRelaxation(intervals, /*demands=*/{}, + /*capacity_upper_bound=*/IntegerValue(1), model, + relaxation); } std::vector AppendLinMaxRelaxation( @@ -734,7 +716,8 @@ std::vector AppendLinMaxRelaxation( IntegerEncoder* encoder = model->GetOrCreate(); GenericLiteralWatcher* watcher = model->GetOrCreate(); - // TODO(user): For the case where num_exprs = 2, Create only one z var. + // TODO(user,user): For the case where num_exprs = 2, Create only one z + // var. std::vector z_vars; std::vector z_lits; z_vars.reserve(num_exprs); @@ -764,7 +747,7 @@ std::vector AppendLinMaxRelaxation( // For the relaxation, we use different constraints with a stronger linear // relaxation as explained in the .h - // TODO(user): Consider passing the x_vars to this method instead of + // TODO(user,user): Consider passing the x_vars to this method instead of // computing it here. std::vector x_vars; for (int i = 0; i < num_exprs; ++i) { @@ -809,11 +792,11 @@ std::vector AppendLinMaxRelaxation( return z_vars; } -void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto, - const int linearization_level, - const Model& model, +void AppendLinearConstraintRelaxation(const ConstraintProto& ct, + bool linearize_enforced_constraints, + Model* model, LinearRelaxation* relaxation) { - auto* mapping = model.Get(); + auto* mapping = model->Get(); // Note that we ignore the holes in the domain. // @@ -822,20 +805,18 @@ void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto, // possible. // // TODO(user): process the "at most one" part of a == 1 separately? - const IntegerValue rhs_domain_min = - IntegerValue(constraint_proto.linear().domain(0)); + const IntegerValue rhs_domain_min = IntegerValue(ct.linear().domain(0)); const IntegerValue rhs_domain_max = - IntegerValue(constraint_proto.linear().domain( - constraint_proto.linear().domain_size() - 1)); + IntegerValue(ct.linear().domain(ct.linear().domain_size() - 1)); if (rhs_domain_min == std::numeric_limits::min() && rhs_domain_max == std::numeric_limits::max()) return; - if (!HasEnforcementLiteral(constraint_proto)) { - LinearConstraintBuilder lc(&model, rhs_domain_min, rhs_domain_max); - for (int i = 0; i < constraint_proto.linear().vars_size(); i++) { - const int ref = constraint_proto.linear().vars(i); - const int64_t coeff = constraint_proto.linear().coeffs(i); + if (!HasEnforcementLiteral(ct)) { + LinearConstraintBuilder lc(model, rhs_domain_min, rhs_domain_max); + for (int i = 0; i < ct.linear().vars_size(); i++) { + const int ref = ct.linear().vars(i); + const int64_t coeff = ct.linear().coeffs(i); lc.AddTerm(mapping->Integer(ref), IntegerValue(coeff)); } relaxation->linear_constraints.push_back(lc.Build()); @@ -843,26 +824,25 @@ void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto, } // Reified version. - if (linearization_level < 2) return; + if (!linearize_enforced_constraints) return; // We linearize fully reified constraints of size 1 all together for a given // variable. But we need to process half-reified ones. - if (!mapping->IsHalfEncodingConstraint(&constraint_proto) && - constraint_proto.linear().vars_size() <= 1) { + if (!mapping->IsHalfEncodingConstraint(&ct) && ct.linear().vars_size() <= 1) { return; } std::vector enforcing_literals; - enforcing_literals.reserve(constraint_proto.enforcement_literal_size()); - for (const int enforcement_ref : constraint_proto.enforcement_literal()) { + enforcing_literals.reserve(ct.enforcement_literal_size()); + for (const int enforcement_ref : ct.enforcement_literal()) { enforcing_literals.push_back(mapping->Literal(enforcement_ref)); } LinearExpression expr; - expr.vars.reserve(constraint_proto.linear().vars_size()); - expr.coeffs.reserve(constraint_proto.linear().vars_size()); - for (int i = 0; i < constraint_proto.linear().vars_size(); i++) { - int ref = constraint_proto.linear().vars(i); - IntegerValue coeff(constraint_proto.linear().coeffs(i)); + expr.vars.reserve(ct.linear().vars_size()); + expr.coeffs.reserve(ct.linear().vars_size()); + for (int i = 0; i < ct.linear().vars_size(); i++) { + int ref = ct.linear().vars(i); + IntegerValue coeff(ct.linear().coeffs(i)); if (!RefIsPositive(ref)) { ref = PositiveRef(ref); coeff = -coeff; @@ -872,7 +852,416 @@ void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto, expr.coeffs.push_back(coeff); } AppendEnforcedLinearExpression(enforcing_literals, expr, rhs_domain_min, - rhs_domain_max, model, relaxation); + rhs_domain_max, *model, relaxation); +} + +// Add a linear relaxation of the CP constraint to the set of linear +// constraints. The highest linearization_level is, the more types of constraint +// we encode. This method should be called only for linearization_level > 0. +// +// Note: IntProd is linearized dynamically using the cut generators. +// +// TODO(user): In full generality, we could encode all the constraint as an LP. +// TODO(user,user): Add unit tests for this method. +void TryToLinearizeConstraint(const CpModelProto& model_proto, + const ConstraintProto& ct, + int linearization_level, Model* model, + LinearRelaxation* relaxation) { + CHECK_EQ(model->GetOrCreate()->CurrentDecisionLevel(), 0); + DCHECK_GT(linearization_level, 0); + + switch (ct.constraint_case()) { + case ConstraintProto::ConstraintCase::kBoolOr: { + if (linearization_level > 1) { + AppendBoolOrRelaxation(ct, model, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kBoolAnd: { + if (linearization_level > 1) { + AppendBoolAndRelaxation(ct, model, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kAtMostOne: { + AppendAtMostOneRelaxation(ct, model, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kExactlyOne: { + AppendExactlyOneRelaxation(ct, model, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kIntMax: { + AppendIntMaxRelaxation(ct, + /*encode_other_direction=*/linearization_level > 1, + model, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kIntMin: { + AppendIntMinRelaxation(ct, + /*encode_other_direction=*/linearization_level > 1, + model, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kLinear: { + AppendLinearConstraintRelaxation( + ct, /*linearize_enforced_constraints=*/linearization_level > 1, model, + relaxation); + break; + } + case ConstraintProto::ConstraintCase::kCircuit: { + AppendCircuitRelaxation(ct, model, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kInterval: { + if (linearization_level > 1) { + AppendIntervalRelaxation(ct, model, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kNoOverlap: { + if (linearization_level > 1) { + AppendNoOverlapRelaxation(model_proto, ct, model, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kCumulative: { + if (linearization_level > 1) { + AppendCumulativeRelaxation(model_proto, ct, model, relaxation); + } + break; + } + default: { + } + } +} + +// Cut generators. + +void AddCircuitCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + std::vector tails(ct.circuit().tails().begin(), + ct.circuit().tails().end()); + std::vector heads(ct.circuit().heads().begin(), + ct.circuit().heads().end()); + auto* mapping = m->GetOrCreate(); + std::vector literals = mapping->Literals(ct.circuit().literals()); + const int num_nodes = ReindexArcs(&tails, &heads); + + relaxation->cut_generators.push_back(CreateStronglyConnectedGraphCutGenerator( + num_nodes, tails, heads, literals, m)); +} + +void AddRoutesCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + std::vector tails(ct.routes().tails().begin(), + ct.routes().tails().end()); + std::vector heads(ct.routes().heads().begin(), + ct.routes().heads().end()); + auto* mapping = m->GetOrCreate(); + std::vector literals = mapping->Literals(ct.routes().literals()); + + int num_nodes = 0; + for (int i = 0; i < ct.routes().tails_size(); ++i) { + num_nodes = std::max(num_nodes, 1 + ct.routes().tails(i)); + num_nodes = std::max(num_nodes, 1 + ct.routes().heads(i)); + } + if (ct.routes().demands().empty() || ct.routes().capacity() == 0) { + relaxation->cut_generators.push_back( + CreateStronglyConnectedGraphCutGenerator(num_nodes, tails, heads, + literals, m)); + } else { + const std::vector demands(ct.routes().demands().begin(), + ct.routes().demands().end()); + relaxation->cut_generators.push_back(CreateCVRPCutGenerator( + num_nodes, tails, heads, literals, demands, ct.routes().capacity(), m)); + } +} + +void AddIntProdCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + if (ct.int_prod().vars_size() != 2) return; + auto* mapping = m->GetOrCreate(); + + // Constraint is z == x * y. + + IntegerVariable z = mapping->Integer(ct.int_prod().target()); + IntegerVariable x = mapping->Integer(ct.int_prod().vars(0)); + IntegerVariable y = mapping->Integer(ct.int_prod().vars(1)); + + IntegerTrail* const integer_trail = m->GetOrCreate(); + IntegerValue x_lb = integer_trail->LowerBound(x); + IntegerValue x_ub = integer_trail->UpperBound(x); + IntegerValue y_lb = integer_trail->LowerBound(y); + IntegerValue y_ub = integer_trail->UpperBound(y); + + if (x == y) { + // We currently only support variables with non-negative domains. + if (x_lb < 0 && x_ub > 0) return; + + // Change the sigh of x if its domain is non-positive. + if (x_ub <= 0) { + x = NegationOf(x); + } + + relaxation->cut_generators.push_back(CreateSquareCutGenerator(z, x, m)); + } else { + // We currently only support variables with non-negative domains. + if (x_lb < 0 && x_ub > 0) return; + if (y_lb < 0 && y_ub > 0) return; + + // Change signs to return to the case where all variables are a domain + // with non negative values only. + if (x_ub <= 0) { + x = NegationOf(x); + z = NegationOf(z); + } + if (y_ub <= 0) { + y = NegationOf(y); + z = NegationOf(z); + } + + relaxation->cut_generators.push_back( + CreatePositiveMultiplicationCutGenerator(z, x, y, m)); + } +} + +void AddAllDiffCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + auto* mapping = m->GetOrCreate(); + const int num_vars = ct.all_diff().vars_size(); + if (num_vars <= m->GetOrCreate()->max_all_diff_cut_size()) { + std::vector vars = mapping->Integers(ct.all_diff().vars()); + relaxation->cut_generators.push_back( + CreateAllDifferentCutGenerator(vars, m)); + } +} + +void AddCumulativeCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + auto* mapping = m->GetOrCreate(); + + std::vector demands = + mapping->Integers(ct.cumulative().demands()); + std::vector intervals = + mapping->Intervals(ct.cumulative().intervals()); + const IntegerVariable capacity = mapping->Integer(ct.cumulative().capacity()); + relaxation->cut_generators.push_back( + CreateCumulativeOverlappingCutGenerator(intervals, capacity, demands, m)); + relaxation->cut_generators.push_back( + CreateCumulativeEnergyCutGenerator(intervals, capacity, demands, m)); + relaxation->cut_generators.push_back( + CreateCumulativeCompletionTimeCutGenerator(intervals, capacity, demands, + m)); +} + +void AddNoOverlapCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + + auto* mapping = m->GetOrCreate(); + std::vector intervals = + mapping->Intervals(ct.no_overlap().intervals()); + relaxation->cut_generators.push_back( + CreateNoOverlapEnergyCutGenerator(intervals, m)); + relaxation->cut_generators.push_back( + CreateNoOverlapPrecedenceCutGenerator(intervals, m)); + relaxation->cut_generators.push_back( + CreateNoOverlapCompletionTimeCutGenerator(intervals, m)); +} + +void AddNoOverlap2dCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (HasEnforcementLiteral(ct)) return; + + auto* mapping = m->GetOrCreate(); + std::vector x_intervals = + mapping->Intervals(ct.no_overlap_2d().x_intervals()); + std::vector y_intervals = + mapping->Intervals(ct.no_overlap_2d().y_intervals()); + relaxation->cut_generators.push_back( + CreateNoOverlap2dCompletionTimeCutGenerator(x_intervals, y_intervals, m)); +} + +void AddLinMaxCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation) { + if (!m->GetOrCreate()->add_lin_max_cuts()) return; + if (HasEnforcementLiteral(ct)) return; + + // TODO(user,user): Support linearization of general target expression. + auto* mapping = m->GetOrCreate(); + if (ct.lin_max().target().vars_size() != 1) return; + if (ct.lin_max().target().coeffs(0) != 1) return; + + const IntegerVariable target = + mapping->Integer(ct.lin_max().target().vars(0)); + std::vector exprs; + exprs.reserve(ct.lin_max().exprs_size()); + for (int i = 0; i < ct.lin_max().exprs_size(); ++i) { + // Note: Cut generator requires all expressions to contain only positive + // vars. + exprs.push_back( + PositiveVarExpr(mapping->GetExprFromProto(ct.lin_max().exprs(i)))); + } + + // FIXME: Add lin max relaxation at level 1. + // Add initial big-M linear relaxation. + // z_vars[i] == 1 <=> target = exprs[i]. + const std::vector z_vars = + AppendLinMaxRelaxation(target, exprs, m, relaxation); + + relaxation->cut_generators.push_back( + CreateLinMaxCutGenerator(target, exprs, z_vars, m)); +} + +// TODO(user): Remove and merge with model loading. +void TryToAddCutGenerators(const ConstraintProto& ct, int linearization_level, + Model* m, LinearRelaxation* relaxation) { + switch (ct.constraint_case()) { + case ConstraintProto::ConstraintCase::kCircuit: { + if (linearization_level > 1) { + AddCircuitCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kRoutes: { + if (linearization_level > 1) { + AddRoutesCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kIntProd: { + AddIntProdCutGenerator(ct, m, relaxation); + break; + } + case ConstraintProto::ConstraintCase::kAllDiff: { + if (linearization_level > 1) { + AddAllDiffCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kCumulative: { + if (linearization_level > 1) { + AddCumulativeCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kNoOverlap: { + if (linearization_level > 1) { + AddNoOverlapCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kNoOverlap2D: { + if (linearization_level > 1) { + AddNoOverlap2dCutGenerator(ct, m, relaxation); + } + break; + } + case ConstraintProto::ConstraintCase::kLinMax: { + if (linearization_level > 1) { + AddLinMaxCutGenerator(ct, m, relaxation); + } + break; + } + default: { + } + } +} + +void ComputeLinearRelaxation(const CpModelProto& model_proto, + int linearization_level, Model* m, + LinearRelaxation* relaxation) { + CHECK(relaxation != nullptr); + + // Linearize the constraints. + absl::flat_hash_set used_integer_variable; + + auto* mapping = m->GetOrCreate(); + auto* encoder = m->GetOrCreate(); + for (const auto& ct : model_proto.constraints()) { + TryToLinearizeConstraint(model_proto, ct, linearization_level, m, + relaxation); + TryToAddCutGenerators(ct, linearization_level, m, relaxation); + } + + // Linearize the encoding of variable that are fully encoded in the proto. + int num_full_encoding_relaxations = 0; + int num_partial_encoding_relaxations = 0; + for (int i = 0; i < model_proto.variables_size(); ++i) { + if (mapping->IsBoolean(i)) continue; + + const IntegerVariable var = mapping->Integer(i); + if (m->Get(IsFixed(var))) continue; + + // TODO(user): This different encoding for the partial variable might be + // better (less LP constraints), but we do need more investigation to + // decide. + if (/* DISABLES CODE */ (false)) { + AppendPartialEncodingRelaxation(var, *m, relaxation); + continue; + } + + if (encoder->VariableIsFullyEncoded(var)) { + if (AppendFullEncodingRelaxation(var, *m, relaxation)) { + ++num_full_encoding_relaxations; + continue; + } + } + + // Even if the variable is fully encoded, sometimes not all its associated + // literal have a view (if they are not part of the original model for + // instance). + // + // TODO(user): Should we add them to the LP anyway? this isn't clear as + // we can sometimes create a lot of Booleans like this. + const int old = relaxation->linear_constraints.size(); + AppendPartialGreaterThanEncodingRelaxation(var, *m, relaxation); + if (relaxation->linear_constraints.size() > old) { + ++num_partial_encoding_relaxations; + } + } + + if (!m->GetOrCreate()->FinishPropagation()) return; + + // Linearize the at most one constraints. Note that we transform them + // into maximum "at most one" first and we removes redundant ones. + m->GetOrCreate()->TransformIntoMaxCliques( + &relaxation->at_most_ones); + for (const std::vector& at_most_one : relaxation->at_most_ones) { + if (at_most_one.empty()) continue; + + LinearConstraintBuilder lc(m, kMinIntegerValue, IntegerValue(1)); + for (const Literal literal : at_most_one) { + // Note that it is okay to simply ignore the literal if it has no + // integer view. + const bool unused ABSL_ATTRIBUTE_UNUSED = + lc.AddLiteralTerm(literal, IntegerValue(1)); + } + relaxation->linear_constraints.push_back(lc.Build()); + } + + // We converted all at_most_one to LP constraints, so we need to clear them + // so that we don't do extra work in the connected component computation. + relaxation->at_most_ones.clear(); + + // Remove size one LP constraints, they are not useful. + relaxation->linear_constraints.erase( + std::remove_if( + relaxation->linear_constraints.begin(), + relaxation->linear_constraints.end(), + [](const LinearConstraint& lc) { return lc.vars.size() <= 1; }), + relaxation->linear_constraints.end()); + + VLOG(3) << "num_full_encoding_relaxations: " << num_full_encoding_relaxations; + VLOG(3) << "num_partial_encoding_relaxations: " + << num_partial_encoding_relaxations; + VLOG(3) << relaxation->linear_constraints.size() + << " constraints in the LP relaxation."; + VLOG(3) << relaxation->cut_generators.size() << " cuts generators."; } } // namespace sat diff --git a/ortools/sat/linear_relaxation.h b/ortools/sat/linear_relaxation.h index bb5a7d76ce..12acdab187 100644 --- a/ortools/sat/linear_relaxation.h +++ b/ortools/sat/linear_relaxation.h @@ -16,7 +16,7 @@ #include -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_constraint.h" #include "ortools/sat/linear_programming_constraint.h" @@ -70,35 +70,6 @@ void AppendPartialGreaterThanEncodingRelaxation(IntegerVariable var, const Model& model, LinearRelaxation* relaxation); -// Adds linearization of different types of constraints. -void TryToLinearizeConstraint(const CpModelProto& model_proto, - const ConstraintProto& ct, Model* model, - int linearization_level, - LinearRelaxation* relaxation); - -// Adds linearization of no overlap constraints. -// It adds an energetic equation linking the duration of all potential tasks to -// the actual span of the no overlap constraint. -void AppendNoOverlapRelaxation(const CpModelProto& model_proto, - const ConstraintProto& ct, - int linearization_level, Model* model, - LinearRelaxation* relaxation); - -// Adds linearization of cumulative constraints.The second part adds an -// energetic equation linking the duration of all potential tasks to the actual -// max span * capacity of the cumulative constraint. -void AppendCumulativeRelaxation(const CpModelProto& model_proto, - const ConstraintProto& ct, - int linearization_level, Model* model, - LinearRelaxation* relaxation); - -// Adds linearization of int max constraints. This can also be used to linearize -// int min with negated variables. -void AppendMaxRelaxation(IntegerVariable target, - const std::vector& vars, - int linearization_level, Model* model, - LinearRelaxation* relaxation); - // Adds linearization of int max constraints. Returns a vector of z vars such // that: z_vars[l] == 1 <=> target = exprs[l]. // @@ -122,6 +93,26 @@ std::vector AppendLinMaxRelaxation( IntegerVariable target, const std::vector& exprs, Model* model, LinearRelaxation* relaxation); +void AppendBoolOrRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +void AppendBoolAndRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +void AppendAtMostOneRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +void AppendExactlyOneRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +void AppendIntMaxRelaxation(const ConstraintProto& ct, + bool encode_other_direction, Model* model, + LinearRelaxation* relaxation); + +void AppendIntMinRelaxation(const ConstraintProto& ct, + bool encode_other_direction, Model* model, + LinearRelaxation* relaxation); + // Appends linear constraints to the relaxation. This also handles the // relaxation of linear constraints with enforcement literals. // A linear constraint lb <= ax <= ub with enforcement literals {ei} is relaxed @@ -130,11 +121,72 @@ std::vector AppendLinMaxRelaxation( // -inf <= (Sum Negated(ei) * (ub - implied_ub)) + ax <= ub // Where implied_lb and implied_ub are trivial lower and upper bounds of the // constraint. -void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto, - const int linearization_level, - const Model& model, +void AppendLinearConstraintRelaxation(const ConstraintProto& ct, + bool linearize_enforced_constraints, + Model* model, LinearRelaxation* relaxation); +void AppendCircuitRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +void AppendIntervalRelaxation(const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +// Adds linearization of no overlap constraints. +// It adds an energetic equation linking the duration of all potential tasks to +// the actual span of the no overlap constraint. +void AppendNoOverlapRelaxation(const CpModelProto& model_proto, + const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +// Adds linearization of cumulative constraints.The second part adds an +// energetic equation linking the duration of all potential tasks to the actual +// max span * capacity of the cumulative constraint. +void AppendCumulativeRelaxation(const CpModelProto& model_proto, + const ConstraintProto& ct, Model* model, + LinearRelaxation* relaxation); + +// Adds linearization of different types of constraints. +void TryToLinearizeConstraint(const CpModelProto& model_proto, + const ConstraintProto& ct, + int linearization_level, Model* model, + LinearRelaxation* relaxation); + +// Cut generators. +void AddCircuitCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddRoutesCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddIntProdCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddAllDiffCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddCumulativeCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddNoOverlapCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddNoOverlap2dCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +void AddLinMaxCutGenerator(const ConstraintProto& ct, Model* m, + LinearRelaxation* relaxation); + +// Scan the model and add cut generators. +void TryToAddCutGenerators(const ConstraintProto& ct, int linearization_level, + Model* m, LinearRelaxation* relaxation); + +// Builds the linear relaxaton of a CpModelProto and stores it in the +// LinearRelaxation container. +void ComputeLinearRelaxation(const CpModelProto& model_proto, + int linearization_level, Model* m, + LinearRelaxation* relaxation); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index d976ce2505..01c2d4a45b 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -105,7 +105,7 @@ void DeleteVectorIndices(const std::vector& indices, Vector* v) { } // In the Fu & Malik algorithm (or in WPM1), when two cores overlap, we -// artifically introduce symmetries. More precisely: +// artificially introduce symmetries. More precisely: // // The picture below shows two cores with index 0 and 1, with one blocking // variable per '-' and with the variables ordered from left to right (by their @@ -701,7 +701,7 @@ SatSolver::Status SolveWithWPM1(LogBehavior log, // If the soft clause protected by old_a has a cost greater than // min_cost then: // - its cost is disminished by min_cost. - // - an identical clause with cost min_cost is artifically added to + // - an identical clause with cost min_cost is artificially added to // the problem. CHECK_GE(costs[index], min_cost); if (costs[index] == min_cost) { @@ -1662,7 +1662,7 @@ SatSolver::Status CoreBasedOptimizer::Optimize() { // Solve under the assumptions. // - // TODO(user): If the "search" is interupted while computing cores, we + // TODO(user): If the "search" is interrupted while computing cores, we // currently do not resume it flawlessly. We however add any cores we found // before aborting. std::vector> cores; diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index b9c187e82a..fe110f8594 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -18,7 +18,7 @@ #include #include "ortools/sat/boolean_problem.pb.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_search.h" #include "ortools/sat/model.h" diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 1c7577008b..42ff91893d 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -1498,45 +1498,5 @@ void PresolveContext::ClearPrecedenceCache() { reified_precedences_cache_.clear(); } -void PresolveContext::BuildIntervalDictionary() { - interval_dictionary_.clear(); - optional_interval_dictionary_.clear(); - start_end_of_intervals_.clear(); - - for (int i = 0; i < working_model->constraints_size(); ++i) { - const ConstraintProto& ct = working_model->constraints(i); - if (ct.constraint_case() != ConstraintProto::ConstraintCase::kInterval) { - continue; - } - start_end_of_intervals_.insert( - std::make_pair(ct.interval().start(), ct.interval().end())); - if (ct.enforcement_literal().empty()) { - interval_dictionary_[std::make_tuple(ct.interval().start(), - ct.interval().end())] = i; - } else { - CHECK_EQ(ct.enforcement_literal_size(), 1); - optional_interval_dictionary_[std::make_tuple( - ct.interval().start(), ct.interval().end(), - ct.enforcement_literal(0))] = i; - } - } -} - -int PresolveContext::GetIntervalFromStartAndEnd(int start, int end) { - const auto& it = interval_dictionary_.find(std::make_tuple(start, end)); - return it == interval_dictionary_.end() ? -1 : it->second; -} - -int PresolveContext::GetIntervalFromStartEndAndPresence(int start, int end, - int presence) { - const auto& it = - optional_interval_dictionary_.find(std::make_tuple(start, end, presence)); - return it == optional_interval_dictionary_.end() ? -1 : it->second; -} - -bool PresolveContext::IsStartAndEndOfOneInterval(int start, int end) { - return start_end_of_intervals_.contains(std::make_pair(start, end)); -} - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 6a71c7f8e4..673a38a7ec 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -366,17 +366,6 @@ class PresolveContext { TimeLimit* time_limit() { return time_limit_; } ModelRandomGenerator* random() { return random_; } - // Scan all interval constraints in the model for later query. - // The index is temporary as it will not be updated after variable - // substitution. - void BuildIntervalDictionary(); - - // These getters perform a lookup in the interval dictionary. - // They return the index of the interval constraint if present, or -1 if not. - int GetIntervalFromStartAndEnd(int start, int end); - int GetIntervalFromStartEndAndPresence(int start, int end, int presence); - bool IsStartAndEndOfOneInterval(int start, int end); - // For each variables, list the constraints that just enforce a lower bound // (resp. upper bound) on that variable. If all the constraints in which a // variable appear are in the same direction, then we can usually fix a @@ -538,12 +527,6 @@ class PresolveContext { // phase, and is cleared afterwards. absl::flat_hash_map, int> reified_precedences_cache_; - - // Cache for intervals. - absl::flat_hash_map, int> interval_dictionary_; - absl::flat_hash_map, int> - optional_interval_dictionary_; - absl::flat_hash_set> start_end_of_intervals_; }; } // namespace sat diff --git a/ortools/sat/rins.cc b/ortools/sat/rins.cc index 4ee247c391..bb5f516a00 100644 --- a/ortools/sat/rins.cc +++ b/ortools/sat/rins.cc @@ -16,7 +16,7 @@ #include #include -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_programming_constraint.h" diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index fc630b0c2f..7a3e698e15 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -984,8 +984,8 @@ message SatParameters { // Temporary flag util the feature is more mature. This convert intervals to // the newer proto format that support affine start/var/end instead of just - // variables. It changes a bit the search and is not always better currently. - optional bool convert_intervals = 177 [default = false]; + // variables. + optional bool convert_intervals = 177 [default = true]; // Whether we try to automatically detect the symmetries in a model and // exploit them. Currently, at level 1 we detect them in presolve and try diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 409fb635ca..beeaea707e 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -18,7 +18,7 @@ #if !defined(__PORTABLE_PLATFORM__) #include "ortools/base/file.h" -#include "ortools/sat/cp_model_loader.h" +#include "ortools/sat/cp_model_mapping.h" #endif // __PORTABLE_PLATFORM__ #include "absl/container/flat_hash_set.h" diff --git a/ortools/sat/timetable.cc b/ortools/sat/timetable.cc index cd1369ca70..fc647a5572 100644 --- a/ortools/sat/timetable.cc +++ b/ortools/sat/timetable.cc @@ -336,6 +336,7 @@ bool TimeTablingPerTask::Propagate() { if (!SweepAllTasks(/*is_forward=*/true)) return false; // We reuse the same profile, but reversed, to update the maximum end times. + if (!helper_->SynchronizeAndSetTimeDirection(false)) return false; ReverseProfile(); // Update the maximum end times (reversed problem). @@ -345,7 +346,7 @@ bool TimeTablingPerTask::Propagate() { } bool TimeTablingPerTask::BuildProfile() { - helper_->SynchronizeAndSetTimeDirection(true); // forward + if (!helper_->SynchronizeAndSetTimeDirection(true)) return false; // Update the set of tasks that contribute to the profile. Tasks that were // contributing are still part of the profile so we only need to check the @@ -430,8 +431,6 @@ bool TimeTablingPerTask::BuildProfile() { } void TimeTablingPerTask::ReverseProfile() { - helper_->SynchronizeAndSetTimeDirection(false); // backward - // We keep the sentinels inchanged. for (int i = 1; i + 1 < profile_.size(); ++i) { profile_[i].start = -profile_[i + 1].start; diff --git a/ortools/sat/timetable_edgefinding.cc b/ortools/sat/timetable_edgefinding.cc index 177160cda4..6da8844746 100644 --- a/ortools/sat/timetable_edgefinding.cc +++ b/ortools/sat/timetable_edgefinding.cc @@ -58,10 +58,10 @@ bool TimeTableEdgeFinding::Propagate() { while (true) { const int64_t old_timestamp = integer_trail_->num_enqueues(); - helper_->SynchronizeAndSetTimeDirection(true); + if (!helper_->SynchronizeAndSetTimeDirection(true)) return false; if (!TimeTableEdgeFindingPass()) return false; - helper_->SynchronizeAndSetTimeDirection(false); + if (!helper_->SynchronizeAndSetTimeDirection(false)) return false; if (!TimeTableEdgeFindingPass()) return false; // Stop if no propagation.