From 9d916b7f554c505fdbe48928acad1ebfd8eef939 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Sat, 29 May 2021 17:56:46 +0200 Subject: [PATCH] [CP-SAT] rewrite diffn propagation; refactor completion time cut generation; add diffn cuts; fix rare bug with cores --- makefiles/Makefile.gen.mk | 27 +++ ortools/sat/BUILD | 15 +- ortools/sat/cp_model_lns.cc | 118 +++++---- ortools/sat/cp_model_solver.cc | 16 +- ortools/sat/cuts.cc | 208 +++++++++++----- ortools/sat/cuts.h | 6 + ortools/sat/diffn.cc | 425 ++++++++++++--------------------- ortools/sat/diffn.h | 37 +-- ortools/sat/diffn_util.cc | 320 +++++++++++++++++++++++++ ortools/sat/diffn_util.h | 102 ++++++++ ortools/sat/intervals.h | 6 + ortools/sat/optimization.cc | 3 +- 12 files changed, 884 insertions(+), 399 deletions(-) create mode 100644 ortools/sat/diffn_util.cc create mode 100644 ortools/sat/diffn_util.h diff --git a/makefiles/Makefile.gen.mk b/makefiles/Makefile.gen.mk index bdfa8d04b2..a4be2d176a 100644 --- a/makefiles/Makefile.gen.mk +++ b/makefiles/Makefile.gen.mk @@ -1335,6 +1335,7 @@ SAT_DEPS = \ $(SRC_DIR)/ortools/sat/cumulative.h \ $(SRC_DIR)/ortools/sat/cuts.h \ $(SRC_DIR)/ortools/sat/diffn.h \ + $(SRC_DIR)/ortools/sat/diffn_util.h \ $(SRC_DIR)/ortools/sat/disjunctive.h \ $(SRC_DIR)/ortools/sat/drat_checker.h \ $(SRC_DIR)/ortools/sat/drat_proof_handler.h \ @@ -1407,6 +1408,7 @@ SAT_LIB_OBJS = \ $(OBJ_DIR)/sat/cumulative_energy.$O \ $(OBJ_DIR)/sat/cuts.$O \ $(OBJ_DIR)/sat/diffn.$O \ + $(OBJ_DIR)/sat/diffn_util.$O \ $(OBJ_DIR)/sat/disjunctive.$O \ $(OBJ_DIR)/sat/drat_checker.$O \ $(OBJ_DIR)/sat/drat_proof_handler.$O \ @@ -1978,6 +1980,31 @@ objs/sat/diffn.$O: ortools/sat/diffn.cc ortools/sat/diffn.h \ ortools/base/stl_util.h ortools/sat/cumulative.h ortools/util/sort.h | $(OBJ_DIR)/sat $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sdiffn.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sdiffn.$O +objs/sat/diffn_util.$O: ortools/sat/diffn_util.cc ortools/sat/diffn_util.h \ + ortools/base/int_type.h ortools/base/macros.h \ + ortools/base/integral_types.h ortools/base/logging.h \ + ortools/base/commandlineflags.h ortools/base/log_severity.h \ + ortools/base/logging_export.h ortools/base/vlog_is_on.h \ + ortools/sat/disjunctive.h ortools/sat/integer.h ortools/base/hash.h \ + ortools/base/basictypes.h ortools/base/map_util.h \ + ortools/base/strong_vector.h ortools/graph/iterators.h \ + ortools/sat/model.h ortools/base/typeid.h ortools/sat/sat_base.h \ + ortools/util/bitset.h ortools/sat/sat_solver.h ortools/base/timer.h \ + ortools/sat/clause.h ortools/sat/drat_proof_handler.h \ + ortools/sat/drat_checker.h ortools/sat/drat_writer.h ortools/base/file.h \ + ortools/gen/ortools/sat/sat_parameters.pb.h ortools/sat/util.h \ + ortools/util/random_engine.h ortools/util/time_limit.h \ + ortools/util/running_stat.h ortools/util/stats.h \ + ortools/sat/pb_constraint.h ortools/sat/restart.h \ + ortools/sat/sat_decision.h ortools/util/integer_pq.h ortools/util/rev.h \ + ortools/util/saturated_arithmetic.h ortools/util/sorted_interval_list.h \ + ortools/sat/intervals.h ortools/sat/cp_constraints.h \ + ortools/sat/integer_expr.h ortools/base/mathutil.h \ + ortools/sat/linear_constraint.h ortools/sat/precedences.h \ + ortools/sat/theta_tree.h ortools/base/iterator_adaptors.h \ + ortools/base/stl_util.h ortools/sat/cumulative.h ortools/util/sort.h | $(OBJ_DIR)/sat + $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sdiffn_util.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sdiffn_util.$O + objs/sat/disjunctive.$O: ortools/sat/disjunctive.cc \ ortools/sat/disjunctive.h ortools/base/int_type.h ortools/base/macros.h \ ortools/sat/integer.h ortools/base/hash.h ortools/base/basictypes.h \ diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index eec5314d50..c6d0fa05e6 100644 --- a/ortools/sat/BUILD +++ b/ortools/sat/BUILD @@ -664,7 +664,6 @@ cc_library( ], ) - cc_library( name = "pseudo_costs", srcs = ["pseudo_costs.cc"], @@ -999,6 +998,7 @@ cc_library( srcs = ["cuts.cc"], hdrs = ["cuts.h"], deps = [ + ":diffn_util", ":implied_bounds", ":integer", ":intervals", @@ -1152,12 +1152,25 @@ cc_library( ], ) +cc_library( + name = "diffn_util", + srcs = ["diffn_util.cc"], + hdrs = ["diffn_util.h"], + deps = [ + ":integer", + ":intervals", + "@com_google_absl//absl/strings:str_format", + "//ortools/graph:connected_components", + ], +) + cc_library( name = "diffn", srcs = ["diffn.cc"], hdrs = ["diffn.h"], deps = [ ":cumulative", + ":diffn_util", ":disjunctive", ":integer", ":intervals", diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index fe00827354..605fca5419 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -644,6 +644,68 @@ Neighborhood ConstraintGraphNeighborhoodGenerator::Generate( return helper_.RelaxGivenVariables(initial_solution, relaxed_variables); } +namespace { +void AddPrecedenceConstraints(const absl::Span intervals, + const absl::flat_hash_set& ignored_intervals, + const CpSolverResponse& initial_solution, + const NeighborhoodGeneratorHelper& helper, + Neighborhood* neighborhood) { + // Sort all non-relaxed intervals of this constraint by current start + // time. + std::vector> start_interval_pairs; + for (const int i : intervals) { + if (ignored_intervals.contains(i)) continue; + 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 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(); + + // 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)) { + // If the end was smaller than the next start, keep it that way. + linear->add_domain(0); + linear->add_vars(before_end); + } 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); + } + linear->add_coeffs(1); + linear->add_vars(after_start); + linear->add_coeffs(-1); + } +} +} // namespace + Neighborhood GenerateSchedulingNeighborhoodForRelaxation( const absl::Span intervals_to_relax, const CpSolverResponse& initial_solution, @@ -687,46 +749,22 @@ Neighborhood GenerateSchedulingNeighborhoodForRelaxation( } for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap)) { - // Sort all non-relaxed intervals of this constraint by current start - // time. - std::vector> start_interval_pairs; - for (const int i : - helper.ModelProto().constraints(c).no_overlap().intervals()) { - if (ignored_intervals.contains(i)) continue; - 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 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_var = helper.ModelProto() - .constraints(start_interval_pairs[i].second) - .interval() - .end(); - const int after_var = helper.ModelProto() - .constraints(start_interval_pairs[i + 1].second) - .interval() - .start(); - CHECK_LE(initial_solution.solution(before_var), - initial_solution.solution(after_var)); - - LinearConstraintProto* linear = - neighborhood.delta.add_constraints()->mutable_linear(); - linear->add_domain(std::numeric_limits::min()); - linear->add_domain(0); - linear->add_vars(before_var); - linear->add_coeffs(1); - linear->add_vars(after_var); - linear->add_coeffs(-1); - } + AddPrecedenceConstraints( + helper.ModelProto().constraints(c).no_overlap().intervals(), + ignored_intervals, initial_solution, helper, &neighborhood); + } + for (const int c : helper.TypeToConstraints(ConstraintProto::kCumulative)) { + AddPrecedenceConstraints( + helper.ModelProto().constraints(c).cumulative().intervals(), + ignored_intervals, initial_solution, helper, &neighborhood); + } + for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap2D)) { + AddPrecedenceConstraints( + helper.ModelProto().constraints(c).no_overlap_2d().x_intervals(), + ignored_intervals, initial_solution, helper, &neighborhood); + AddPrecedenceConstraints( + helper.ModelProto().constraints(c).no_overlap_2d().y_intervals(), + ignored_intervals, initial_solution, helper, &neighborhood); } // Set the current solution as a hint. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 7bfae3e98f..9ec80c3ea6 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -588,6 +588,18 @@ void TryToAddCutGenerators(const CpModelProto& model_proto, 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; @@ -2870,7 +2882,9 @@ void SolveCpModelParallel(const CpModelProto& model_proto, helper, absl::StrCat("graph_cst_lns_", local_params.name())), local_params, helper, &shared)); - if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap).empty()) { + if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap).empty() || + !helper->TypeToConstraints(ConstraintProto::kNoOverlap2D).empty() || + !helper->TypeToConstraints(ConstraintProto::kCumulative).empty()) { subsolvers.push_back(absl::make_unique( absl::make_unique( helper, absl::StrCat("scheduling_time_window_lns_", diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 5dd0366db4..385473f3a5 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -26,11 +26,13 @@ #include "ortools/base/integral_types.h" #include "ortools/base/stl_util.h" #include "ortools/base/strong_vector.h" +#include "ortools/sat/diffn_util.h" #include "ortools/sat/integer.h" #include "ortools/sat/intervals.h" #include "ortools/sat/linear_constraint.h" #include "ortools/sat/linear_constraint_manager.h" #include "ortools/sat/sat_base.h" +#include "ortools/sat/util.h" #include "ortools/util/time_limit.h" namespace operations_research { @@ -2413,6 +2415,8 @@ struct CtEvent { IntegerValue start_min; IntegerValue size_min; double lp_end; + IntegerValue y_min; + IntegerValue y_max; }; // We generate the cut from the Smith's rule from: @@ -2427,14 +2431,10 @@ struct CtEvent { // // A second difference is that we look at a set of intervals starting // after a given start_min, sorted by relative (end_lp - start_min). -// -// We actually compute the cuts with all the sizes actually equal to (size / -// size_divisor), but to keep the computation in the integer domain, we multiply -// by size_divisor where needed instead. void GenerateCompletionTimeCut( const std::string& cut_name, const absl::StrongVector& lp_values, - IntegerValue size_divisor, std::vector events, Model* model, + std::vector events, Model* model, LinearConstraintManager* manager) { TopNCuts top_n_cuts(15); @@ -2461,36 +2461,47 @@ void GenerateCompletionTimeCut( IntegerValue best_min_contrib(0); IntegerValue sum_duration(0); IntegerValue sum_square_duration(0); - double lp_contrib = 0; + IntegerValue best_size_divisor(0); + double unscaled_lp_contrib = 0; IntegerValue current_start_min(kMaxIntegerValue); + IntegerValue y_min = kMaxIntegerValue; + IntegerValue y_max = kMinIntegerValue; for (int i = 0; i < residual_tasks.size(); ++i) { DCHECK_GE(residual_tasks[i].start_min, sequence_start_min); const IntegerValue duration = residual_tasks[i].size_min; sum_duration += duration; sum_square_duration += duration * duration; - lp_contrib += - residual_tasks[i].lp_end * ToDouble(duration * size_divisor); + unscaled_lp_contrib += residual_tasks[i].lp_end * ToDouble(duration); current_start_min = std::min(current_start_min, residual_tasks[i].start_min); + y_min = std::min(y_min, residual_tasks[i].y_min); + y_max = std::max(y_max, residual_tasks[i].y_max); + const IntegerValue size_divisor = y_max - y_min; + + // We compute the cuts with all the sizes actually equal to + // (size / size_divisor), but to keep the computation in the integer + // domain, we multiply by size_divisor where needed instead. const IntegerValue min_contrib = (sum_duration * sum_duration + sum_square_duration) / 2 + current_start_min * sum_duration * size_divisor; - const double efficacy = (ToDouble(min_contrib) - lp_contrib) / + const double efficacy = (ToDouble(min_contrib) - + unscaled_lp_contrib * ToDouble(size_divisor)) / std::sqrt(ToDouble(sum_square_duration)); // TODO(user): Check overflow and ignore if too big. if (efficacy > best_efficacy) { best_efficacy = efficacy; best_end = i; best_min_contrib = min_contrib; + best_size_divisor = size_divisor; } } if (best_end != -1) { LinearConstraintBuilder cut(model, best_min_contrib, kMaxIntegerValue); for (int i = 0; i <= best_end; ++i) { cut.AddTerm(residual_tasks[i].end, - residual_tasks[i].size_min * size_divisor); + residual_tasks[i].size_min * best_size_divisor); } top_n_cuts.AddCut(cut.Build(), cut_name, lp_values); } @@ -2516,29 +2527,27 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return; - std::vector events; - std::vector reverse_events; - - for (int index = 0; index < helper->NumTasks(); ++index) { - if (!helper->IsPresent(index)) continue; - const IntegerValue size_min = helper->SizeMin(index); - if (size_min > 0) { - const AffineExpression end_expr = helper->Ends()[index]; - events.push_back({end_expr, helper->StartMin(index), size_min, - end_expr.LpValue(lp_values)}); - const AffineExpression r_start_expr = - helper->Starts()[index].Negated(); - reverse_events.push_back({r_start_expr, -helper->EndMax(index), - size_min, - r_start_expr.LpValue(lp_values)}); + auto generate_cuts = [&lp_values, model, manager, + helper](const std::string& cut_name) { + std::vector events; + for (int index = 0; index < helper->NumTasks(); ++index) { + if (!helper->IsPresent(index)) continue; + const IntegerValue size_min = helper->SizeMin(index); + if (size_min > 0) { + const AffineExpression end_expr = helper->Ends()[index]; + events.push_back({end_expr, helper->StartMin(index), size_min, + end_expr.LpValue(lp_values), IntegerValue(0), + IntegerValue(1)}); + } } - } - GenerateCompletionTimeCut("NoOverlapCompletionTime", lp_values, - IntegerValue(1), std::move(events), model, - manager); - GenerateCompletionTimeCut("NoOverlapCompletionTimeMirror", lp_values, - IntegerValue(1), std::move(reverse_events), - model, manager); + GenerateCompletionTimeCut(cut_name, lp_values, std::move(events), + model, manager); + }; + + helper->SynchronizeAndSetTimeDirection(true); + generate_cuts("NoOverlapCompletionTime"); + helper->SynchronizeAndSetTimeDirection(false); + generate_cuts("NoOverlapCompletionTimeMirror"); }; return result; } @@ -2566,32 +2575,117 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return; - std::vector events; - std::vector reverse_events; - - for (int index = 0; index < helper->NumTasks(); ++index) { - if (!helper->IsPresent(index)) continue; - const IntegerValue area_min = - helper->SizeMin(index) * - integer_trail->LowerBound(demands[index]); - if (area_min > 0) { - const AffineExpression end_expr = helper->Ends()[index]; - events.push_back({end_expr, helper->StartMin(index), area_min, - end_expr.LpValue(lp_values)}); - const AffineExpression r_start_expr = - helper->Starts()[index].Negated(); - reverse_events.push_back({r_start_expr, -helper->EndMax(index), - area_min, - r_start_expr.LpValue(lp_values)}); - } - } const IntegerValue capacity_max = integer_trail->UpperBound(capacity); - GenerateCompletionTimeCut("CumulativeCompletionTime", lp_values, - capacity_max, std::move(events), model, - manager); - GenerateCompletionTimeCut("CumulativeCompletionTimeMirror", lp_values, - capacity_max, std::move(reverse_events), - model, manager); + auto generate_cuts = [&lp_values, model, manager, helper, capacity_max, + integer_trail, + &demands](const std::string& cut_name) { + std::vector events; + for (int index = 0; index < helper->NumTasks(); ++index) { + if (!helper->IsPresent(index)) continue; + const IntegerValue area_min = + helper->SizeMin(index) * + integer_trail->LowerBound(demands[index]); + if (area_min > 0) { + const AffineExpression end_expr = helper->Ends()[index]; + events.push_back({end_expr, helper->StartMin(index), area_min, + end_expr.LpValue(lp_values), IntegerValue(0), + capacity_max}); + } + } + GenerateCompletionTimeCut(cut_name, lp_values, std::move(events), + model, manager); + }; + + helper->SynchronizeAndSetTimeDirection(true); + generate_cuts("CumulativeCompletionTime"); + helper->SynchronizeAndSetTimeDirection(false); + generate_cuts("CumulativeCompletionTimeMirror"); + }; + return result; +} + +CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( + const std::vector& x_intervals, + const std::vector& y_intervals, Model* model) { + CutGenerator result; + + SchedulingConstraintHelper* x_helper = + new SchedulingConstraintHelper(x_intervals, model); + model->TakeOwnership(x_helper); + + SchedulingConstraintHelper* y_helper = + new SchedulingConstraintHelper(y_intervals, model); + model->TakeOwnership(y_helper); + AddIntegerVariableFromIntervals(x_helper, model, &result.vars); + AddIntegerVariableFromIntervals(y_helper, model, &result.vars); + + Trail* trail = model->GetOrCreate(); + + result.generate_cuts = + [trail, x_helper, y_helper, model]( + const absl::StrongVector& lp_values, + LinearConstraintManager* manager) { + if (trail->CurrentDecisionLevel() > 0) return; + + x_helper->SynchronizeAndSetTimeDirection(true); + y_helper->SynchronizeAndSetTimeDirection(true); + + const int num_boxes = x_helper->NumTasks(); + std::vector active_boxes; + std::vector cached_areas(num_boxes); + std::vector cached_rectangles(num_boxes); + for (int box = 0; box < num_boxes; ++box) { + cached_areas[box] = x_helper->SizeMin(box) * y_helper->SizeMin(box); + if (cached_areas[box] == 0) continue; + if (!y_helper->IsPresent(box) || !y_helper->IsPresent(box)) continue; + + // TODO(user): Also consider shifted end max. + Rectangle& rectangle = cached_rectangles[box]; + rectangle.x_min = x_helper->ShiftedStartMin(box); + rectangle.x_max = x_helper->EndMax(box); + rectangle.y_min = y_helper->ShiftedStartMin(box); + rectangle.y_max = y_helper->EndMax(box); + + active_boxes.push_back(box); + } + + if (active_boxes.size() <= 1) return; + + std::vector> components = + GetOverlappingRectangleComponents(cached_rectangles, + absl::MakeSpan(active_boxes)); + for (absl::Span boxes : components) { + if (boxes.size() <= 1) continue; + + auto generate_cuts = [&lp_values, model, manager, &boxes, + &cached_areas]( + const std::string& cut_name, + SchedulingConstraintHelper* x_helper, + SchedulingConstraintHelper* y_helper) { + std::vector events; + + for (const int box : boxes) { + const AffineExpression x_end_expr = x_helper->Ends()[box]; + events.push_back({x_end_expr, x_helper->ShiftedStartMin(box), + cached_areas[box], + x_end_expr.LpValue(lp_values), + y_helper->ShiftedStartMin(box), + y_helper->ShiftedEndMax(box)}); + } + + GenerateCompletionTimeCut(cut_name, lp_values, std::move(events), + model, manager); + }; + + x_helper->SynchronizeAndSetTimeDirection(true); + y_helper->SynchronizeAndSetTimeDirection(true); + generate_cuts("NoOverlap2dXCompletionTime", x_helper, y_helper); + generate_cuts("NoOverlap2dYCompletionTime", y_helper, x_helper); + x_helper->SynchronizeAndSetTimeDirection(false); + y_helper->SynchronizeAndSetTimeDirection(false); + generate_cuts("NoOverlap2dXCompletionTimeMirror", x_helper, y_helper); + generate_cuts("NoOverlap2dYCompletionTimeMirror", y_helper, x_helper); + } }; return result; } diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 7fdd00ae4e..7310298c4b 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -515,6 +515,12 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( const IntegerVariable capacity, const std::vector& demands, Model* model); +// Completion time cuts for the no_overlap_2d constraint. It actually generates +// the completion time cumulative cuts in both axis. +CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( + const std::vector& x_intervals, + const std::vector& y_intervals, Model* model); + // For a given set of intervals, we first compute the min and max of all // intervals. Then we create a cut that indicates that all intervals must fit // in that span. diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 16e47c898c..a4aaced1de 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -23,6 +23,7 @@ #include "ortools/base/map_util.h" #include "ortools/base/stl_util.h" #include "ortools/sat/cumulative.h" +#include "ortools/sat/diffn_util.h" #include "ortools/sat/disjunctive.h" #include "ortools/sat/intervals.h" #include "ortools/sat/sat_solver.h" @@ -107,6 +108,156 @@ void AddCumulativeRelaxation(const std::vector& x_intervals, model->Add(Cumulative(x_intervals, sizes, capacity, x)); } +#define RETURN_IF_FALSE(f) \ + if (!(f)) return false; + +NonOverlappingRectanglesEnergyPropagator:: + ~NonOverlappingRectanglesEnergyPropagator() {} + +bool NonOverlappingRectanglesEnergyPropagator::Propagate() { + const int num_boxes = x_.NumTasks(); + x_.SynchronizeAndSetTimeDirection(true); + y_.SynchronizeAndSetTimeDirection(true); + + active_boxes_.clear(); + cached_areas_.resize(num_boxes); + cached_rectangles_.resize(num_boxes); + for (int box = 0; box < num_boxes; ++box) { + cached_areas_[box] = x_.SizeMin(box) * y_.SizeMin(box); + if (cached_areas_[box] == 0) continue; + if (!x_.IsPresent(box) || !y_.IsPresent(box)) continue; + + // TODO(user): Also consider shifted end max. + Rectangle& rectangle = cached_rectangles_[box]; + rectangle.x_min = x_.ShiftedStartMin(box); + rectangle.x_max = x_.EndMax(box); + rectangle.y_min = y_.ShiftedStartMin(box); + rectangle.y_max = y_.EndMax(box); + + active_boxes_.push_back(box); + } + + if (active_boxes_.size() <= 1) return true; + + Rectangle conflicting_rectangle; + std::vector> components = GetOverlappingRectangleComponents( + cached_rectangles_, absl::MakeSpan(active_boxes_)); + for (absl::Span boxes : components) { + // Computes the size on x and y past which there is no point doing any + // energetic reasonning. We do a few iterations since as we filter one size, + // we can potentially filter more on the transposed dimension. + threshold_x_ = kMaxIntegerValue; + threshold_y_ = kMaxIntegerValue; + for (int i = 0; i < 3; ++i) { + if (!AnalyzeIntervals(/*transpose=*/i == 1, boxes, cached_rectangles_, + cached_areas_, &threshold_x_, &threshold_y_, + &conflicting_rectangle)) { + return ReportEnergyConflict(conflicting_rectangle, boxes, &x_, &y_); + } + boxes = FilterBoxesAndRandomize(cached_rectangles_, boxes, threshold_x_, + threshold_y_, *random_); + if (boxes.size() <= 1) break; + } + if (boxes.size() <= 1) continue; + + // Now that we removed boxes with a large domain, resplit everything. + const std::vector> local_components = + GetOverlappingRectangleComponents(cached_rectangles_, boxes); + for (absl::Span local_boxes : local_components) { + IntegerValue total_sum_of_areas(0); + for (const int box : local_boxes) { + total_sum_of_areas += cached_areas_[box]; + } + for (const int box : local_boxes) { + RETURN_IF_FALSE( + FailWhenEnergyIsTooLarge(box, local_boxes, total_sum_of_areas)); + } + } + } + + return true; +} + +int NonOverlappingRectanglesEnergyPropagator::RegisterWith( + GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + x_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, + /*watch_end_max=*/true); + y_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, + /*watch_end_max=*/true); + return id; +} + +void NonOverlappingRectanglesEnergyPropagator::SortBoxesIntoNeighbors( + int box, absl::Span local_boxes, + IntegerValue total_sum_of_areas) { + const Rectangle& box_dim = cached_rectangles_[box]; + + neighbors_.clear(); + for (const int other_box : local_boxes) { + if (other_box == box) continue; + const Rectangle& other_dim = cached_rectangles_[other_box]; + const IntegerValue span_x = std::max(box_dim.x_max, other_dim.x_max) - + std::min(box_dim.x_min, other_dim.x_min); + if (span_x > threshold_x_) continue; + const IntegerValue span_y = std::max(box_dim.y_max, other_dim.y_max) - + std::min(box_dim.y_min, other_dim.y_min); + if (span_y > threshold_y_) continue; + const IntegerValue bounding_area = span_x * span_y; + if (bounding_area < total_sum_of_areas) { + neighbors_.push_back({other_box, bounding_area}); + } + } + std::sort(neighbors_.begin(), neighbors_.end()); +} + +bool NonOverlappingRectanglesEnergyPropagator::FailWhenEnergyIsTooLarge( + int box, absl::Span local_boxes, + IntegerValue total_sum_of_areas) { + SortBoxesIntoNeighbors(box, local_boxes, total_sum_of_areas); + + Rectangle area = cached_rectangles_[box]; + 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); + }; + + for (int i = 0; i < neighbors_.size(); ++i) { + const int other_box = neighbors_[i].box; + CHECK_GT(cached_areas_[other_box], 0); + + // Update Bounding box. + area.TakeUnionWith(cached_rectangles_[other_box]); + if (area.x_max - area.x_min > threshold_x_) break; + if (area.y_max - area.y_min > threshold_y_) break; + + // Update sum of areas. + sum_of_areas += cached_areas_[other_box]; + const IntegerValue bounding_area = + (area.x_max - area.x_min) * (area.y_max - area.y_min); + if (bounding_area >= total_sum_of_areas) { + // Nothing will be deduced. Exiting. + return true; + } + + if (sum_of_areas > bounding_area) { + x_.ClearReason(); + y_.ClearReason(); + add_box_energy_in_rectangle_reason(box); + for (int j = 0; j <= i; ++j) { + add_box_energy_in_rectangle_reason(neighbors_[j].box); + } + x_.ImportOtherReasons(y_); + return x_.ReportConflict(); + } + } + return true; +} + namespace { // We want for different propagation to reuse as much as possible the same @@ -138,15 +289,16 @@ void SplitDisjointBoxes(const SchedulingConstraintHelper& x, absl::Span boxes, std::vector>* result) { result->clear(); - std::sort(boxes.begin(), boxes.end(), - [&x](int a, int b) { return x.StartMin(a) < x.StartMin(b); }); + std::sort(boxes.begin(), boxes.end(), [&x](int a, int b) { + return x.ShiftedStartMin(a) < x.ShiftedStartMin(b); + }); int current_start = 0; std::size_t current_length = 1; IntegerValue current_max_end = x.EndMax(boxes[0]); for (int b = 1; b < boxes.size(); ++b) { const int box = boxes[b]; - if (x.StartMin(box) < current_max_end) { + if (x.ShiftedStartMin(box) < current_max_end) { // Merge. current_length++; current_max_end = std::max(current_max_end, x.EndMax(box)); @@ -168,273 +320,6 @@ void SplitDisjointBoxes(const SchedulingConstraintHelper& x, } // namespace -#define RETURN_IF_FALSE(f) \ - if (!(f)) return false; - -NonOverlappingRectanglesEnergyPropagator:: - ~NonOverlappingRectanglesEnergyPropagator() {} - -bool NonOverlappingRectanglesEnergyPropagator::AnalyzeIntervals( - SchedulingConstraintHelper* x, SchedulingConstraintHelper* y, - IntegerValue x_threshold, IntegerValue* y_threshold) { - // First, we compute the possible x_min values (removing duplicates). - std::vector starts; - const int num_tasks = x->NumTasks(); - std::vector task_is_active(num_tasks, false); - for (int t = 0; t < num_tasks; ++t) { - if (!x->IsPresent(t) || !y->IsPresent(t)) continue; - - // Ignore all the box with a domain that is too large. - const IntegerValue x_min = x->ShiftedStartMin(t); - const IntegerValue x_max = x->EndMax(t); - if (x_max - x_min > x_threshold) continue; - - const IntegerValue energy = x->SizeMin(t) * y->SizeMin(t); - if (energy == 0) continue; - - const IntegerValue y_min = y->ShiftedStartMin(t); - const IntegerValue y_max = y->EndMax(t); - if (y_max - y_min > *y_threshold) continue; - - task_is_active[t] = true; - starts.push_back(x_min); - } - gtl::STLSortAndRemoveDuplicates(&starts); - - // The maximum y dimension of a bounding area for which there is a potential - // conflict. - IntegerValue max_conflict_height(0); - - // This is currently only used for logging. - absl::flat_hash_set> stripes; - - // All quantities at index j correspond to the interval [starts[j], x_max]. - std::vector energies(starts.size(), IntegerValue(0)); - std::vector min_height(starts.size(), kMaxIntegerValue); - std::vector y_mins(starts.size(), kMaxIntegerValue); - std::vector y_maxs(starts.size(), -kMaxIntegerValue); - std::vector energy_at_max_y(starts.size(), IntegerValue(0)); - std::vector energy_at_min_y(starts.size(), IntegerValue(0)); - - // Iterate over all boxes by increasing x_max values. - const std::vector& by_end_max = x->TaskByDecreasingEndMax(); - for (int i = by_end_max.size(); --i >= 0;) { - const int t = by_end_max[i].task_index; - if (!task_is_active[t]) continue; - - const IntegerValue x_min = x->ShiftedStartMin(t); - const IntegerValue x_max = x->EndMax(t); - const IntegerValue energy = x->SizeMin(t) * y->SizeMin(t); - const IntegerValue y_min = y->ShiftedStartMin(t); - const IntegerValue y_max = y->EndMax(t); - - // Add this box contribution to all the [starts[j], x_max] intervals. - for (int j = 0; j < starts.size(); ++j) { - if (x_min < starts[j]) continue; - if (x_max - starts[j] > x_threshold) continue; - - energies[j] += energy; - min_height[j] = std::min(min_height[j], y_max - y_min); - - if (y_min < y_mins[j]) { - y_mins[j] = y_min; - energy_at_min_y[j] = energy; - } else if (y_min == y_mins[j]) { - energy_at_min_y[j] += energy; - } - - if (y_max > y_maxs[j]) { - y_maxs[j] = y_max; - energy_at_max_y[j] = energy; - } else if (y_max == y_maxs[j]) { - energy_at_max_y[j] += energy; - } - - const IntegerValue y_height = y_maxs[j] - y_mins[j]; - IntegerValue conflict_height = - CeilRatio(energies[j], x_max - starts[j]) - 1; - if (conflict_height >= y_height) { - // We have a conflict. - x->ClearReason(); - y->ClearReason(); - for (int k = i; k < by_end_max.size(); ++k) { - const int task_index = by_end_max[k].task_index; - if (!task_is_active[task_index]) continue; - if (x->ShiftedStartMin(task_index) < starts[j]) continue; - - x->AddEnergyAfterReason(task_index, x->SizeMin(task_index), - starts[j]); - x->AddEndMaxReason(task_index, x_max); - - y->AddEnergyAfterReason(task_index, y->SizeMin(task_index), - y_mins[j]); - y->AddEndMaxReason(task_index, y_maxs[j]); - } - x->ImportOtherReasons(*y); - return x->ReportConflict(); - } - - // The only way to have a conflict is to reduce the y_domain, thus at - // least removing either the boxes at the max or at the min. - conflict_height = CeilRatio(energies[j] - std::min(energy_at_max_y[j], - energy_at_min_y[j]), - x_max - starts[j]) - - 1; - - // If this is true, there can never be a conflict for any subset of the - // boxes in the current [starts[j], x_max] interval. - if (min_height[j] > conflict_height) continue; - - if (VLOG_IS_ON(2)) stripes.insert({starts[j], x_max}); - max_conflict_height = std::max(max_conflict_height, conflict_height); - } - } - - VLOG(2) << starts.size() << " " << x->NumTasks() - << " conflict_height: " << max_conflict_height - << " num_stripes:" << stripes.size() << " (<= " << x_threshold << ")"; - - *y_threshold = std::min(*y_threshold, max_conflict_height); - return true; -} - -bool NonOverlappingRectanglesEnergyPropagator::Propagate() { - const int num_boxes = x_.NumTasks(); - x_.SynchronizeAndSetTimeDirection(true); - y_.SynchronizeAndSetTimeDirection(true); - - // Computes the size on x and y past which there is no point doing any - // energetic reasonning. We do two iterations since as we filter one size, - // we can potentially filter more on the other. - // - // TODO(user): This is visible in the cpu profile, optimize more. - threshold_x_ = kMaxIntegerValue; - threshold_y_ = kMaxIntegerValue; - RETURN_IF_FALSE(AnalyzeIntervals(&x_, &y_, threshold_x_, &threshold_y_)); - RETURN_IF_FALSE(AnalyzeIntervals(&y_, &x_, threshold_y_, &threshold_x_)); - RETURN_IF_FALSE(AnalyzeIntervals(&x_, &y_, threshold_x_, &threshold_y_)); - - active_boxes_.clear(); - cached_areas_.resize(num_boxes); - cached_dimensions_.resize(num_boxes); - for (int box = 0; box < num_boxes; ++box) { - cached_areas_[box] = x_.SizeMin(box) * y_.SizeMin(box); - if (cached_areas_[box] == 0) continue; - if (!x_.IsPresent(box) || !y_.IsPresent(box)) continue; - - // TODO(user): Also consider shifted end max. - Dimension& dimension = cached_dimensions_[box]; - dimension.x_min = x_.ShiftedStartMin(box); - dimension.x_max = x_.EndMax(box); - dimension.y_min = y_.ShiftedStartMin(box); - dimension.y_max = y_.EndMax(box); - - // We can skip box with a domain that is too large. - if (dimension.x_max - dimension.x_min > threshold_x_) continue; - if (dimension.y_max - dimension.y_min > threshold_y_) continue; - - active_boxes_.push_back(box); - } - if (active_boxes_.size() <= 1) return true; - - SplitDisjointBoxes(x_, absl::MakeSpan(active_boxes_), &x_split_); - for (absl::Span x_boxes : x_split_) { - SplitDisjointBoxes(y_, x_boxes, &y_split_); - for (absl::Span y_boxes : y_split_) { - IntegerValue total_sum_of_areas(0); - for (const int box : y_boxes) { - total_sum_of_areas += cached_areas_[box]; - } - for (const int box : y_boxes) { - RETURN_IF_FALSE( - FailWhenEnergyIsTooLarge(box, y_boxes, total_sum_of_areas)); - } - } - } - - return true; -} - -int NonOverlappingRectanglesEnergyPropagator::RegisterWith( - GenericLiteralWatcher* watcher) { - const int id = watcher->Register(this); - x_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, - /*watch_end_max=*/true); - y_.WatchAllTasks(id, watcher, /*watch_start_max=*/false, - /*watch_end_max=*/true); - return id; -} - -void NonOverlappingRectanglesEnergyPropagator::SortBoxesIntoNeighbors( - int box, absl::Span local_boxes, - IntegerValue total_sum_of_areas) { - const Dimension& box_dim = cached_dimensions_[box]; - - neighbors_.clear(); - for (const int other_box : local_boxes) { - if (other_box == box) continue; - const Dimension& other_dim = cached_dimensions_[other_box]; - const IntegerValue span_x = std::max(box_dim.x_max, other_dim.x_max) - - std::min(box_dim.x_min, other_dim.x_min); - if (span_x > threshold_x_) continue; - const IntegerValue span_y = std::max(box_dim.y_max, other_dim.y_max) - - std::min(box_dim.y_min, other_dim.y_min); - if (span_y > threshold_y_) continue; - const IntegerValue bounding_area = span_x * span_y; - if (bounding_area < total_sum_of_areas) { - neighbors_.push_back({other_box, bounding_area}); - } - } - std::sort(neighbors_.begin(), neighbors_.end()); -} - -bool NonOverlappingRectanglesEnergyPropagator::FailWhenEnergyIsTooLarge( - int box, absl::Span local_boxes, - IntegerValue total_sum_of_areas) { - SortBoxesIntoNeighbors(box, local_boxes, total_sum_of_areas); - - Dimension area = cached_dimensions_[box]; - 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); - }; - - for (int i = 0; i < neighbors_.size(); ++i) { - const int other_box = neighbors_[i].box; - CHECK_GT(cached_areas_[other_box], 0); - - // Update Bounding box. - area.TakeUnionWith(cached_dimensions_[other_box]); - if (area.x_max - area.x_min > threshold_x_) break; - if (area.y_max - area.y_min > threshold_y_) break; - - // Update sum of areas. - sum_of_areas += cached_areas_[other_box]; - const IntegerValue bounding_area = - (area.x_max - area.x_min) * (area.y_max - area.y_min); - if (bounding_area >= total_sum_of_areas) { - // Nothing will be deduced. Exiting. - return true; - } - - if (sum_of_areas > bounding_area) { - x_.ClearReason(); - y_.ClearReason(); - add_box_energy_in_rectangle_reason(box); - for (int j = 0; j <= i; ++j) { - add_box_energy_in_rectangle_reason(neighbors_[j].box); - } - x_.ImportOtherReasons(y_); - return x_.ReportConflict(); - } - } - return true; -} - // Note that x_ and y_ must be initialized with enough intervals when passed // to the disjunctive propagators. NonOverlappingRectanglesDisjunctivePropagator:: diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h index b8475dda07..224684f8e1 100644 --- a/ortools/sat/diffn.h +++ b/ortools/sat/diffn.h @@ -20,6 +20,7 @@ #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/macros.h" +#include "ortools/sat/diffn_util.h" #include "ortools/sat/disjunctive.h" #include "ortools/sat/integer.h" #include "ortools/sat/intervals.h" @@ -36,21 +37,15 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { // boxes. If strict is true, these boxes must not 'cross' another box, and are // pushed by the other boxes. NonOverlappingRectanglesEnergyPropagator(SchedulingConstraintHelper* x, - SchedulingConstraintHelper* y) - : x_(*x), y_(*y) {} + SchedulingConstraintHelper* y, + Model* model) + : x_(*x), y_(*y), random_(model->GetOrCreate()) {} ~NonOverlappingRectanglesEnergyPropagator() override; bool Propagate() final; int RegisterWith(GenericLiteralWatcher* watcher); private: - // A O(n^2) algorithm to analyze all the relevant [x_min, x_max] intervals - // and infer a threshold of the y size of a bounding box after which there - // is no point checking for energy overload. - bool AnalyzeIntervals(SchedulingConstraintHelper* x, - SchedulingConstraintHelper* y, IntegerValue x_threshold, - IntegerValue* y_threshold); - void SortBoxesIntoNeighbors(int box, absl::Span local_boxes, IntegerValue total_sum_of_areas); bool FailWhenEnergyIsTooLarge(int box, absl::Span local_boxes, @@ -58,32 +53,16 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { SchedulingConstraintHelper& x_; SchedulingConstraintHelper& y_; + ModelRandomGenerator* random_; // When the size of the bounding box is greater than any of the corresponding - // dimension, then there is no point checking for overload. + // rectangles, then there is no point checking for overload. IntegerValue threshold_x_; IntegerValue threshold_y_; - std::vector> x_split_; - std::vector> y_split_; - std::vector active_boxes_; std::vector cached_areas_; - - struct Dimension { - IntegerValue x_min; - IntegerValue x_max; - IntegerValue y_min; - IntegerValue y_max; - - void TakeUnionWith(const Dimension& other) { - x_min = std::min(x_min, other.x_min); - y_min = std::min(y_min, other.y_min); - x_max = std::max(x_max, other.x_max); - y_max = std::max(y_max, other.y_max); - } - }; - std::vector cached_dimensions_; + std::vector cached_rectangles_; struct Neighbor { int box; @@ -177,7 +156,7 @@ inline std::function NonOverlappingRectangles( model->TakeOwnership(y_helper); NonOverlappingRectanglesEnergyPropagator* energy_constraint = - new NonOverlappingRectanglesEnergyPropagator(x_helper, y_helper); + new NonOverlappingRectanglesEnergyPropagator(x_helper, y_helper, model); GenericLiteralWatcher* const watcher = model->GetOrCreate(); watcher->SetPropagatorPriority(energy_constraint->RegisterWith(watcher), 3); diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc new file mode 100644 index 0000000000..19e0e373b5 --- /dev/null +++ b/ortools/sat/diffn_util.cc @@ -0,0 +1,320 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/sat/diffn_util.h" + +#include "absl/algorithm/container.h" +#include "ortools/base/stl_util.h" + +namespace operations_research { +namespace sat { + +bool Rectangle::IsDisjoint(const Rectangle& other) const { + return x_min >= other.x_max || other.x_min >= x_max || y_min >= other.y_max || + other.y_min >= y_max; +} + +std::vector> GetOverlappingRectangleComponents( + const std::vector& rectangles, + absl::Span active_rectangles) { + if (active_rectangles.empty()) return {}; + + std::vector> result; + const int size = active_rectangles.size(); + for (int start = 0; start < size;) { + // Find the component of active_rectangles[start]. + int end = start + 1; + for (int i = start; i < end; i++) { + for (int j = end; j < size; ++j) { + if (!rectangles[active_rectangles[i]].IsDisjoint( + rectangles[active_rectangles[j]])) { + std::swap(active_rectangles[end++], active_rectangles[j]); + } + } + } + if (end > start + 1) { + result.push_back(active_rectangles.subspan(start, end - start)); + } + start = end; + } + return result; +} + +bool ReportEnergyConflict(Rectangle bounding_box, absl::Span boxes, + SchedulingConstraintHelper* x, + SchedulingConstraintHelper* y) { + x->ClearReason(); + y->ClearReason(); + IntegerValue total_energy(0); + for (const int b : boxes) { + const IntegerValue x_min = x->ShiftedStartMin(b); + const IntegerValue x_max = x->EndMax(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); + 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); + + total_energy += x->SizeMin(b) * y->SizeMin(b); + + // We abort early if a subset of boxes is enough. + // TODO(user): Also relax the box if possible. + if (total_energy > bounding_box.Area()) break; + } + + CHECK_GT(total_energy, bounding_box.Area()); + x->ImportOtherReasons(*y); + return x->ReportConflict(); +} + +bool BoxesAreInEnergyConflict(const std::vector& rectangles, + const std::vector& energies, + absl::Span boxes, + Rectangle* conflict) { + // First consider all relevant intervals along the x axis. + std::vector x_starts; + std::vector boxes_by_increasing_x_max; + for (const int b : boxes) { + x_starts.push_back(rectangles[b].x_min); + boxes_by_increasing_x_max.push_back({b, rectangles[b].x_max}); + } + gtl::STLSortAndRemoveDuplicates(&x_starts); + std::sort(boxes_by_increasing_x_max.begin(), boxes_by_increasing_x_max.end()); + + std::vector y_starts; + std::vector energy_sum; + std::vector boxes_by_increasing_y_max; + + std::vector> stripes(x_starts.size()); + for (int i = 0; i < boxes_by_increasing_x_max.size(); ++i) { + const int b = boxes_by_increasing_x_max[i].task_index; + const IntegerValue x_min = rectangles[b].x_min; + const IntegerValue x_max = rectangles[b].x_max; + for (int j = 0; j < x_starts.size(); ++j) { + if (x_starts[j] > x_min) break; + stripes[j].push_back(b); + + // Redo the same on the y coordinate for the current x interval + // which is [starts[j], x_max]. + y_starts.clear(); + boxes_by_increasing_y_max.clear(); + for (const int b : stripes[j]) { + y_starts.push_back(rectangles[b].y_min); + boxes_by_increasing_y_max.push_back({b, rectangles[b].y_max}); + } + gtl::STLSortAndRemoveDuplicates(&y_starts); + std::sort(boxes_by_increasing_y_max.begin(), + boxes_by_increasing_y_max.end()); + + const IntegerValue x_size = x_max - x_starts[j]; + energy_sum.assign(y_starts.size(), IntegerValue(0)); + for (int i = 0; i < boxes_by_increasing_y_max.size(); ++i) { + const int b = boxes_by_increasing_y_max[i].task_index; + const IntegerValue y_min = rectangles[b].y_min; + const IntegerValue y_max = rectangles[b].y_max; + for (int j = 0; j < y_starts.size(); ++j) { + if (y_starts[j] > y_min) break; + energy_sum[j] += energies[b]; + if (energy_sum[j] > x_size * (y_max - y_starts[j])) { + if (conflict != nullptr) { + *conflict = rectangles[b]; + for (int k = 0; k < i; ++k) { + const int task_index = boxes_by_increasing_y_max[k].task_index; + if (rectangles[task_index].y_min >= y_starts[j]) { + conflict->TakeUnionWith(rectangles[task_index]); + } + } + } + return true; + } + } + } + } + } + return false; +} + +bool AnalyzeIntervals(bool transpose, absl::Span local_boxes, + const std::vector& rectangles, + const std::vector& rectangle_energies, + IntegerValue* x_threshold, IntegerValue* y_threshold, + Rectangle* conflict) { + // First, we compute the possible x_min values (removing duplicates). + // We also sort the relevant tasks by their x_max. + // + // TODO(user): If the number of unique x_max is smaller than the number of + // unique x_min, it is better to do it the other way around. + std::vector starts; + std::vector task_by_increasing_x_max; + for (const int t : local_boxes) { + const IntegerValue x_min = + transpose ? rectangles[t].y_min : rectangles[t].x_min; + const IntegerValue x_max = + transpose ? rectangles[t].y_max : rectangles[t].x_max; + starts.push_back(x_min); + task_by_increasing_x_max.push_back({t, x_max}); + } + gtl::STLSortAndRemoveDuplicates(&starts); + + // Note that for the same end_max, the order change our heuristic to + // evaluate the max_conflict_height. + std::sort(task_by_increasing_x_max.begin(), task_by_increasing_x_max.end()); + + // The maximum y dimension of a bounding area for which there is a potential + // conflict. + IntegerValue max_conflict_height(0); + + // This is currently only used for logging. + absl::flat_hash_set> stripes; + + // All quantities at index j correspond to the interval [starts[j], x_max]. + std::vector energies(starts.size(), IntegerValue(0)); + std::vector y_mins(starts.size(), kMaxIntegerValue); + std::vector y_maxs(starts.size(), -kMaxIntegerValue); + std::vector energy_at_max_y(starts.size(), IntegerValue(0)); + std::vector energy_at_min_y(starts.size(), IntegerValue(0)); + + // Sentinel. + starts.push_back(kMaxIntegerValue); + + // Iterate over all boxes by increasing x_max values. + int first_j = 0; + const IntegerValue threshold = transpose ? *y_threshold : *x_threshold; + for (int i = 0; i < task_by_increasing_x_max.size(); ++i) { + const int t = task_by_increasing_x_max[i].task_index; + + const IntegerValue energy = rectangle_energies[t]; + IntegerValue x_min = rectangles[t].x_min; + IntegerValue x_max = rectangles[t].x_max; + IntegerValue y_min = rectangles[t].y_min; + IntegerValue y_max = rectangles[t].y_max; + if (transpose) { + std::swap(x_min, y_min); + std::swap(x_max, y_max); + } + + // Add this box contribution to all the [starts[j], x_max] intervals. + while (first_j + 1 < starts.size() && x_max - starts[first_j] > threshold) { + ++first_j; + } + for (int j = first_j; starts[j] <= x_min; ++j) { + const IntegerValue old_energy_at_max = energy_at_max_y[j]; + const IntegerValue old_energy_at_min = energy_at_min_y[j]; + + energies[j] += energy; + + const bool is_disjoint = y_min >= y_maxs[j] || y_max <= y_mins[j]; + + if (y_min <= y_mins[j]) { + if (y_min < y_mins[j]) { + y_mins[j] = y_min; + energy_at_min_y[j] = energy; + } else { + energy_at_min_y[j] += energy; + } + } + + if (y_max >= y_maxs[j]) { + if (y_max > y_maxs[j]) { + y_maxs[j] = y_max; + energy_at_max_y[j] = energy; + } else { + energy_at_max_y[j] += energy; + } + } + + // If the new box is disjoint in y from the ones added so far, there + // cannot be a new conflict involving this box, so we skip until we add + // new boxes. + if (is_disjoint) continue; + + const IntegerValue width = x_max - starts[j]; + IntegerValue conflict_height = CeilRatio(energies[j], width) - 1; + if (y_max - y_min > conflict_height) continue; + if (conflict_height >= y_maxs[j] - y_mins[j]) { + // We have a conflict. + if (conflict != nullptr) { + *conflict = rectangles[t]; + for (int k = 0; k < i; ++k) { + const int task_index = task_by_increasing_x_max[k].task_index; + const IntegerValue task_x_min = transpose + ? rectangles[task_index].y_min + : rectangles[task_index].x_min; + if (task_x_min < starts[j]) continue; + conflict->TakeUnionWith(rectangles[task_index]); + } + } + return false; + } + + // Because we currently do not have a conflict involving the new box, the + // only way to have one is to remove enough energy to reduce the y domain. + IntegerValue can_remove = std::min(old_energy_at_min, old_energy_at_max); + if (old_energy_at_min < old_energy_at_max) { + if (y_maxs[j] - y_min >= + CeilRatio(energies[j] - old_energy_at_min, width)) { + // In this case, we need to remove at least old_energy_at_max to have + // a conflict. + can_remove = old_energy_at_max; + } + } else if (old_energy_at_max < old_energy_at_min) { + if (y_max - y_mins[j] >= + CeilRatio(energies[j] - old_energy_at_max, width)) { + can_remove = old_energy_at_min; + } + } + conflict_height = CeilRatio(energies[j] - can_remove, width) - 1; + + // If the new box height is above the conflict_height, do not count + // it now. We only need to consider conflict involving the new box. + if (y_max - y_min > conflict_height) continue; + + if (VLOG_IS_ON(2)) stripes.insert({starts[j], x_max}); + max_conflict_height = std::max(max_conflict_height, conflict_height); + } + } + + VLOG(2) << " num_starts: " << starts.size() - 1 << "/" << local_boxes.size() + << " conflict_height: " << max_conflict_height + << " num_stripes:" << stripes.size() << " (<= " << threshold << ")"; + + if (transpose) { + *x_threshold = std::min(*x_threshold, max_conflict_height); + } else { + *y_threshold = std::min(*y_threshold, max_conflict_height); + } + return true; +} + +absl::Span FilterBoxesAndRandomize( + const std::vector& cached_rectangles, absl::Span boxes, + IntegerValue threshold_x, IntegerValue threshold_y, + absl::BitGenRef random) { + size_t new_size = 0; + for (const int b : boxes) { + const Rectangle& dim = cached_rectangles[b]; + if (dim.x_max - dim.x_min > threshold_x) continue; + if (dim.y_max - dim.y_min > threshold_y) continue; + boxes[new_size++] = b; + } + if (new_size == 0) return {}; + std::shuffle(&boxes[0], &boxes[0] + new_size, random); + return {&boxes[0], new_size}; +} + +} // namespace sat +} // namespace operations_research diff --git a/ortools/sat/diffn_util.h b/ortools/sat/diffn_util.h new file mode 100644 index 0000000000..2dc6f25d0f --- /dev/null +++ b/ortools/sat/diffn_util.h @@ -0,0 +1,102 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_SAT_DIFFN_UTIL_H_ +#define OR_TOOLS_SAT_DIFFN_UTIL_H_ + +#include +#include + +#include "absl/container/flat_hash_set.h" +#include "absl/strings/str_format.h" +#include "ortools/graph/connected_components.h" +#include "ortools/sat/integer.h" +#include "ortools/sat/intervals.h" + +namespace operations_research { +namespace sat { + +struct Rectangle { + IntegerValue x_min; + IntegerValue x_max; + IntegerValue y_min; + IntegerValue y_max; + + void TakeUnionWith(const Rectangle& other) { + x_min = std::min(x_min, other.x_min); + y_min = std::min(y_min, other.y_min); + x_max = std::max(x_max, other.x_max); + y_max = std::max(y_max, other.y_max); + } + + IntegerValue Area() const { return (x_max - x_min) * (y_max - y_min); } + + bool IsDisjoint(const Rectangle& other) const; + + std::string DebugString() const { + return absl::StrFormat("rectangle(x(%i..%i), y(%i..%i))", x_min.value(), + x_max.value(), y_min.value(), y_max.value()); + } +}; + +// Creates a graph when two nodes are connected iif their rectangles overlap. +// Then partition into connected components. +// +// This method removes all singleton components. It will modify the +// active_rectangle span in place. +std::vector> GetOverlappingRectangleComponents( + const std::vector& rectangles, + absl::Span active_rectangles); + +// Visible for testing. The algo is in O(n^4) so shouldn't be used directly. +// Returns true if there exist a bounding box with too much energy. +bool BoxesAreInEnergyConflict(const std::vector& rectangles, + const std::vector& energies, + absl::Span boxes, + Rectangle* conflict = nullptr); + +// Checks that there is indeed a conflict for the given bounding_box and +// report it. This returns false for convenience as we usually want to return +// false on a conflict. +// +// TODO(user): relax the bounding box dimension to have a relaxed explanation. +// We can also minimize the number of required intervals. +bool ReportEnergyConflict(Rectangle bounding_box, absl::Span boxes, + SchedulingConstraintHelper* x, + SchedulingConstraintHelper* y); + +// A O(n^2) algorithm to analyze all the relevant X intervals and infer a +// threshold of the y size of a bounding box after which there is no point +// checking for energy overload. +// +// Returns false on conflict, and fill the bounding box that caused the +// conflict. +// +// If transpose is true, we analyze the relevant Y intervals instead. +bool AnalyzeIntervals(bool transpose, absl::Span boxes, + const std::vector& rectangles, + const std::vector& rectangle_energies, + IntegerValue* x_threshold, IntegerValue* y_threshold, + Rectangle* conflict = nullptr); + +// Removes boxes with a size above the thresholds. Also randomize the order. +// Because we rely on various heuristic, this allow to change the order from +// one call to the next. +absl::Span FilterBoxesAndRandomize( + const std::vector& cached_rectangles, absl::Span boxes, + IntegerValue threshold_x, IntegerValue threshold_y, absl::BitGenRef random); + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_DIFFN_UTIL_H_ diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index a0aae7dd03..0456557369 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -249,6 +249,12 @@ class SchedulingConstraintHelper : public PropagatorInterface, return cached_shifted_start_min_[t]; } + // As with ShiftedStartMin(), we can compute the shifted end max (that is + // start_max + size_min. + IntegerValue ShiftedEndMax(int t) const { + return -cached_negated_shifted_end_max_[t]; + } + bool StartIsFixed(int t) const; bool EndIsFixed(int t) const; bool SizeIsFixed(int t) const; diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index 0d8d952814..d976ce2505 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -224,6 +224,7 @@ void MinimizeCoreWithPropagation(TimeLimit* limit, SatSolver* solver, solver->Backtrack(0); solver->SetAssumptionLevel(0); + if (!solver->FinishPropagation()) return; while (!limit->LimitReached()) { // We want each literal in candidate to appear last once in our propagation // order. We want to do that while maximizing the reutilization of the @@ -1223,7 +1224,7 @@ SatSolver::Status FindCores(std::vector assumptions, if (sat_solver->parameters().minimize_core()) { MinimizeCoreWithPropagation(limit, sat_solver, &core); } - CHECK(!core.empty()); + if (core.empty()) return sat_solver->UnsatStatus(); cores->push_back(core); if (!sat_solver->parameters().find_multiple_cores()) break;