From 7e98e7b4d4b4d47ea27797d5fc5f685c2a9033be Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Mon, 10 Feb 2025 12:45:46 +0100 Subject: [PATCH] [CP-SAT] improve no_overlap_2d cuts; more no_overlap_2d cuts; bugfixes; add sat/c_api subdirectory --- ortools/base/macros.h | 2 +- ortools/graph/shortest_paths.h | 68 +----- ortools/init/BUILD.bazel | 2 +- ortools/init/init.cc | 33 ++- ortools/linear_solver/model_exporter.cc | 22 +- ortools/pdlp/primal_dual_hybrid_gradient.cc | 2 +- ortools/sat/2d_rectangle_presolve.cc | 250 ++++++++++++++------ ortools/sat/2d_rectangle_presolve.h | 64 +++++ ortools/sat/2d_rectangle_presolve_test.cc | 87 +++++++ ortools/sat/BUILD.bazel | 8 +- ortools/sat/CMakeLists.txt | 2 +- ortools/sat/all_different_test.cc | 6 +- ortools/sat/c_api/BUILD.bazel | 31 +++ ortools/sat/c_api/cp_solver_c.cc | 83 +++++++ ortools/sat/c_api/cp_solver_c.h | 41 ++++ ortools/sat/cp_model_expand.cc | 4 +- ortools/sat/cp_model_presolve.cc | 141 +++++++++++ ortools/sat/cp_model_presolve.h | 13 + ortools/sat/cp_model_solver.cc | 8 +- ortools/sat/cp_model_solver_helpers.cc | 3 +- ortools/sat/diffn_cuts.cc | 152 ++++++------ ortools/sat/diffn_cuts.h | 1 - ortools/sat/go/cpmodel/BUILD.bazel | 17 +- ortools/sat/go/cpmodel/cp_solver.go | 2 +- ortools/sat/linear_relaxation.cc | 24 +- ortools/sat/optimization.cc | 10 +- ortools/sat/sat_parameters.proto | 16 +- ortools/sat/scheduling_cuts.cc | 34 +-- ortools/sat/solution_crush.cc | 83 +++++++ ortools/sat/solution_crush.h | 17 ++ ortools/util/fixed_shape_binary_tree.h | 1 + 31 files changed, 927 insertions(+), 300 deletions(-) create mode 100644 ortools/sat/c_api/BUILD.bazel create mode 100644 ortools/sat/c_api/cp_solver_c.cc create mode 100644 ortools/sat/c_api/cp_solver_c.h diff --git a/ortools/base/macros.h b/ortools/base/macros.h index 85dc470feb..0a67d92e6d 100644 --- a/ortools/base/macros.h +++ b/ortools/base/macros.h @@ -16,7 +16,7 @@ #include // for size_t. -#include "ortools/base/base_export.h" // for OR_DLL +#include "ortools/base/base_export.h" // for OR_DLL #define COMPILE_ASSERT(x, msg) diff --git a/ortools/graph/shortest_paths.h b/ortools/graph/shortest_paths.h index 6fcdcbb095..c117e99d16 100644 --- a/ortools/graph/shortest_paths.h +++ b/ortools/graph/shortest_paths.h @@ -116,10 +116,6 @@ class GenericPathContainer { using NodeIndex = typename GraphType::NodeIndex; using Impl = internal::PathContainerImpl; - // TODO(b/385094969): Remove this when all clients are migrated, and use - // factory functions instead. - GenericPathContainer(); - // This type is neither copyable nor movable. GenericPathContainer(const GenericPathContainer&) = delete; GenericPathContainer& operator=(const GenericPathContainer&) = delete; @@ -153,9 +149,6 @@ class GenericPathContainer { // Builds a path container which only stores distances between path nodes. static GenericPathContainer BuildPathDistanceContainer(); - ABSL_DEPRECATED("Use factory function BuildPathDistanceContainer instead.") - static void BuildPathDistanceContainer(GenericPathContainer* path_container); - // Builds a path container which stores explicit paths and distances between // path nodes in a memory-compact representation. // In this case `GetPenultimateNodeInPath()` is `O(log(path_tree_size))`, @@ -168,11 +161,6 @@ class GenericPathContainer { // `O(log(path_tree_size) * path_size)`. static GenericPathContainer BuildInMemoryCompactPathContainer(); - ABSL_DEPRECATED( - "Use factory function BuildInMemoryCompactPathContainer instead.") - static void BuildInMemoryCompactPathContainer( - GenericPathContainer* path_container); - // TODO(user): Add save-to-disk container. // TODO(user): Add `BuildInMemoryFastPathContainer()`, which does // `GetPenultimateNodeInPath()` in `O(1)`. @@ -230,12 +218,12 @@ void ComputeOneToAllShortestPaths( // Computes shortest paths from the node `source` to nodes in `destinations`. // TODO(b/385094969): Remove second template parameter when all clients are // migrated. -template +template void ComputeOneToManyShortestPaths( const GraphType& graph, const std::vector& arc_lengths, typename GraphType::NodeIndex source, const std::vector& destinations, - GenericPathContainer* const path_container) { + GenericPathContainer* const path_container) { std::vector sources(1, source); ComputeManyToManyShortestPathsWithMultipleThreads( graph, arc_lengths, sources, destinations, 1, path_container); @@ -282,13 +270,10 @@ void ComputeManyToAllShortestPathsWithMultipleThreads( } // Computes shortest paths between all nodes of the graph. -// TODO(b/385094969): Remove second template parameter when all clients are -// migrated. -template +template void ComputeAllToAllShortestPathsWithMultipleThreads( const GraphType& graph, const std::vector& arc_lengths, - int num_threads, - GenericPathContainer* const path_container) { + int num_threads, GenericPathContainer* const path_container) { std::vector all_nodes; GetGraphNodesFromGraph(graph, &all_nodes); ComputeManyToManyShortestPathsWithMultipleThreads( @@ -635,15 +620,13 @@ bool InsertOrUpdateEntry( // using a binary heap-based Dijkstra algorithm. // TODO(user): Investigate alternate implementation which wouldn't use // AdjustablePriorityQueue. -// TODO(b/385094969): Remove second template parameter when all clients are -// migrated. -template +template void ComputeOneToManyOnGraph( const GraphType* const graph, const std::vector* const arc_lengths, typename GraphType::NodeIndex source, const std::vector* const destinations, - typename GenericPathContainer::Impl* const paths) { + typename GenericPathContainer::Impl* const paths) { using NodeIndex = typename GraphType::NodeIndex; using ArcIndex = typename GraphType::ArcIndex; using NodeEntryT = NodeEntry; @@ -715,9 +698,6 @@ void ComputeOneToManyOnGraph( } // namespace internal -template -GenericPathContainer::GenericPathContainer() = default; - template GenericPathContainer::~GenericPathContainer() = default; @@ -744,22 +724,6 @@ void GenericPathContainer::GetPath( container_->GetPath(from, to, path); } -template -void GenericPathContainer::BuildPathDistanceContainer( - GenericPathContainer* const path_container) { - CHECK(path_container != nullptr); - path_container->container_ = std::make_unique< - internal::DistanceContainer>(); -} - -template -void GenericPathContainer::BuildInMemoryCompactPathContainer( - GenericPathContainer* const path_container) { - CHECK(path_container != nullptr); - path_container->container_ = std::make_unique< - internal::InMemoryCompactPathContainer>(); -} - template GenericPathContainer GenericPathContainer::BuildPathDistanceContainer() { @@ -776,22 +740,12 @@ GenericPathContainer::BuildInMemoryCompactPathContainer() { NodeIndex, GraphType::kNilNode>>()); } -// TODO(b/385094969): Remove second template parameter when all clients are -// migrated. -template +template void ComputeManyToManyShortestPathsWithMultipleThreads( const GraphType& graph, const std::vector& arc_lengths, const std::vector& sources, const std::vector& destinations, - int num_threads, - GenericPathContainer* const paths) { - static_assert(std::is_same_v, - "use an explicit `GenericPathContainer` instead of using " - "`PathContainer`"); - static_assert(GraphType::kNilNode == PathContainerGraphType::kNilNode, - "use an explicit `GenericPathContainer` instead of using " - "`PathContainer`"); + int num_threads, GenericPathContainer* const paths) { if (graph.num_nodes() > 0) { CHECK_EQ(graph.num_arcs(), arc_lengths.size()) << "Number of arcs in graph must match arc length vector size"; @@ -812,10 +766,8 @@ void ComputeManyToManyShortestPathsWithMultipleThreads( pool->StartWorkers(); for (int i = 0; i < unique_sources.size(); ++i) { pool->Schedule(absl::bind_front( - &internal::ComputeOneToManyOnGraph, - &graph, &arc_lengths, unique_sources[i], &unique_destinations, - container)); + &internal::ComputeOneToManyOnGraph, &graph, &arc_lengths, + unique_sources[i], &unique_destinations, container)); } } container->Finalize(); diff --git a/ortools/init/BUILD.bazel b/ortools/init/BUILD.bazel index 95b55c918c..47a185fa3d 100644 --- a/ortools/init/BUILD.bazel +++ b/ortools/init/BUILD.bazel @@ -15,8 +15,8 @@ package(default_visibility = ["//visibility:public"]) cc_library( name = "init", - hdrs = ["init.h"], srcs = ["init.cc"], + hdrs = ["init.h"], deps = [ "//ortools/base", "//ortools/gurobi:environment", diff --git a/ortools/init/init.cc b/ortools/init/init.cc index c02a33f709..e80ce87a4f 100644 --- a/ortools/init/init.cc +++ b/ortools/init/init.cc @@ -22,25 +22,24 @@ #include "ortools/sat/cp_model_solver_helpers.h" namespace operations_research { - void CppBridge::InitLogging(const std::string& usage) { - absl::SetProgramUsageMessage(usage); - absl::InitializeLog(); - } +void CppBridge::InitLogging(const std::string& usage) { + absl::SetProgramUsageMessage(usage); + absl::InitializeLog(); +} - void CppBridge::SetFlags(const CppFlags& flags) { - absl::SetFlag(&FLAGS_stderrthreshold, flags.stderrthreshold); - absl::EnableLogPrefix(flags.log_prefix); - if (!flags.cp_model_dump_prefix.empty()) { - absl::SetFlag(&FLAGS_cp_model_dump_prefix, flags.cp_model_dump_prefix); - } - absl::SetFlag(&FLAGS_cp_model_dump_models, flags.cp_model_dump_models); - absl::SetFlag(&FLAGS_cp_model_dump_submodels, - flags.cp_model_dump_submodels); - absl::SetFlag(&FLAGS_cp_model_dump_response, flags.cp_model_dump_response); +void CppBridge::SetFlags(const CppFlags& flags) { + absl::SetFlag(&FLAGS_stderrthreshold, flags.stderrthreshold); + absl::EnableLogPrefix(flags.log_prefix); + if (!flags.cp_model_dump_prefix.empty()) { + absl::SetFlag(&FLAGS_cp_model_dump_prefix, flags.cp_model_dump_prefix); } + absl::SetFlag(&FLAGS_cp_model_dump_models, flags.cp_model_dump_models); + absl::SetFlag(&FLAGS_cp_model_dump_submodels, flags.cp_model_dump_submodels); + absl::SetFlag(&FLAGS_cp_model_dump_response, flags.cp_model_dump_response); +} - bool CppBridge::LoadGurobiSharedLibrary(const std::string& full_library_path) { - return LoadGurobiDynamicLibrary({full_library_path}).ok(); - } +bool CppBridge::LoadGurobiSharedLibrary(const std::string& full_library_path) { + return LoadGurobiDynamicLibrary({full_library_path}).ok(); +} } // namespace operations_research diff --git a/ortools/linear_solver/model_exporter.cc b/ortools/linear_solver/model_exporter.cc index cd4ecd5983..23fd23ff5b 100644 --- a/ortools/linear_solver/model_exporter.cc +++ b/ortools/linear_solver/model_exporter.cc @@ -140,7 +140,7 @@ class MPModelProtoExporter { // into two constraints, one for the left hand side (_lhs) and one for right // hand side (_rhs). bool AppendConstraint(const MPConstraintProto& ct_proto, - const std::string& name, LineBreaker& line_breaker, + absl::string_view name, LineBreaker& line_breaker, std::vector& show_variable, std::string* output); // Clears "output" and writes a term to it, in "LP" format. Returns false on @@ -159,8 +159,8 @@ class MPModelProtoExporter { // Same as AppendMpsLineHeader. Appends an extra new-line at the end the // string pointed to by output. - void AppendMpsLineHeaderWithNewLine(const std::string& id, - const std::string& name, + void AppendMpsLineHeaderWithNewLine(absl::string_view id, + absl::string_view name, std::string* output) const; // Appends an MPS term in various contexts. The term consists of a head name, @@ -186,7 +186,7 @@ class MPModelProtoExporter { // Appends a line describing the bound of a variablenew-line if two columns // are already present on the MPS line. // Used by and in complement to AppendMpsTermWithContext. - void AppendMpsBound(const std::string& bound_type, const std::string& name, + void AppendMpsBound(absl::string_view bound_type, absl::string_view name, double value, std::string* output) const; const MPModelProto& proto_; @@ -392,7 +392,7 @@ std::string DoubleToString(double d) { return absl::StrCat((d)); } } // namespace bool MPModelProtoExporter::AppendConstraint(const MPConstraintProto& ct_proto, - const std::string& name, + absl::string_view name, LineBreaker& line_breaker, std::vector& show_variable, std::string* output) { @@ -414,7 +414,7 @@ bool MPModelProtoExporter::AppendConstraint(const MPConstraintProto& ct_proto, absl::StrAppend(output, " ", name, ": ", line_breaker.GetOutput()); } else { if (ub != +kInfinity) { - std::string rhs_name = name; + std::string rhs_name(name); if (lb != -kInfinity) { absl::StrAppend(&rhs_name, "_rhs"); } @@ -427,7 +427,7 @@ bool MPModelProtoExporter::AppendConstraint(const MPConstraintProto& ct_proto, absl::StrAppend(output, relation); } if (lb != -kInfinity) { - std::string lhs_name = name; + std::string lhs_name(name); if (ub != +kInfinity) { absl::StrAppend(&lhs_name, "_lhs"); } @@ -462,7 +462,7 @@ bool IsBoolean(const MPVariableProto& var) { floor(var.upper_bound()) == 1.0; } -void UpdateMaxSize(const std::string& new_string, int* size) { +void UpdateMaxSize(absl::string_view new_string, int* size) { const int new_size = new_string.size(); if (new_size > *size) *size = new_size; } @@ -687,7 +687,7 @@ void MPModelProtoExporter::AppendMpsLineHeader(absl::string_view id, } void MPModelProtoExporter::AppendMpsLineHeaderWithNewLine( - const std::string& id, const std::string& name, std::string* output) const { + absl::string_view id, absl::string_view name, std::string* output) const { AppendMpsLineHeader(id, name, output); absl::StripTrailingAsciiWhitespace(output); absl::StrAppend(output, "\n"); @@ -704,8 +704,8 @@ void MPModelProtoExporter::AppendMpsTermWithContext(absl::string_view head_name, AppendNewLineIfTwoColumns(output); } -void MPModelProtoExporter::AppendMpsBound(const std::string& bound_type, - const std::string& name, double value, +void MPModelProtoExporter::AppendMpsBound(absl::string_view bound_type, + absl::string_view name, double value, std::string* output) const { AppendMpsLineHeader(bound_type, "BOUND", output); AppendMpsPair(name, value, output); diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc index ef0c9a18e5..a027b46936 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc @@ -1262,7 +1262,7 @@ std::optional PreprocessSolver::ApplyPresolveIfEnabled( // set it for completeness. presolved_qp->objective_scaling_factor = glop_lp.objective_scaling_factor(); sharded_qp_ = ShardedQuadraticProgram(std::move(*presolved_qp), num_threads_, - num_shards_); + num_shards_, params.scheduler_type()); // A status of `INIT` means the preprocessor created a (usually) smaller // problem that needs solving. Other statuses mean the preprocessor solved // the problem completely. diff --git a/ortools/sat/2d_rectangle_presolve.cc b/ortools/sat/2d_rectangle_presolve.cc index ebc940b3c0..b755bcd58c 100644 --- a/ortools/sat/2d_rectangle_presolve.cc +++ b/ortools/sat/2d_rectangle_presolve.cc @@ -40,85 +40,14 @@ namespace operations_research { namespace sat { -bool PresolveFixed2dRectangles( +namespace { +std::vector FindSpacesThatCannotBeOccupied( + absl::Span fixed_boxes, absl::Span non_fixed_boxes, - std::vector* fixed_boxes) { - // This implementation compiles a set of areas that cannot be occupied by any - // item, then calls ReduceNumberofBoxes() to use these areas to minimize - // `fixed_boxes`. - bool changed = false; - - DCHECK(FindPartialRectangleIntersections(*fixed_boxes).empty()); - IntegerValue original_area = 0; - std::vector fixed_boxes_copy; - if (VLOG_IS_ON(1)) { - for (const Rectangle& r : *fixed_boxes) { - original_area += r.Area(); - } - } - if (VLOG_IS_ON(2)) { - fixed_boxes_copy = *fixed_boxes; - } - - const int original_num_boxes = fixed_boxes->size(); - - // The greedy algorithm is really fast. Run it first since it might greatly - // reduce the size of large trivial instances. - std::vector empty_vec; - if (ReduceNumberofBoxesGreedy(fixed_boxes, &empty_vec)) { - changed = true; - } - - IntegerValue min_x_size = std::numeric_limits::max(); - IntegerValue min_y_size = std::numeric_limits::max(); - - CHECK(!non_fixed_boxes.empty()); - Rectangle bounding_box = non_fixed_boxes[0].bounding_area; - - for (const RectangleInRange& box : non_fixed_boxes) { - bounding_box.GrowToInclude(box.bounding_area); - min_x_size = std::min(min_x_size, box.x_size); - min_y_size = std::min(min_y_size, box.y_size); - } - DCHECK_GT(min_x_size, 0); - DCHECK_GT(min_y_size, 0); - - // Fixed items are only useful to constraint where the non-fixed items can be - // placed. This means in particular that any part of a fixed item outside the - // bounding box of the non-fixed items is useless. Clip them. - int new_size = 0; - while (new_size < fixed_boxes->size()) { - Rectangle& rectangle = (*fixed_boxes)[new_size]; - DCHECK_GT(rectangle.SizeX(), 0); - DCHECK_GT(rectangle.SizeY(), 0); - if (rectangle.x_min < bounding_box.x_min) { - rectangle.x_min = bounding_box.x_min; - changed = true; - } - if (rectangle.x_max > bounding_box.x_max) { - rectangle.x_max = bounding_box.x_max; - changed = true; - } - if (rectangle.y_min < bounding_box.y_min) { - rectangle.y_min = bounding_box.y_min; - changed = true; - } - if (rectangle.y_max > bounding_box.y_max) { - rectangle.y_max = bounding_box.y_max; - changed = true; - } - if (rectangle.SizeX() <= 0 || rectangle.SizeY() <= 0) { - // The whole rectangle was outside of the domain, remove it. - std::swap(rectangle, (*fixed_boxes)[fixed_boxes->size() - 1]); - fixed_boxes->resize(fixed_boxes->size() - 1); - changed = true; - continue; - } else { - new_size++; - } - } - - std::vector optional_boxes = *fixed_boxes; + const Rectangle& bounding_box, IntegerValue min_x_size, + IntegerValue min_y_size) { + std::vector optional_boxes = {fixed_boxes.begin(), + fixed_boxes.end()}; if (bounding_box.x_min > std::numeric_limits::min() && bounding_box.y_min > std::numeric_limits::min() && @@ -228,6 +157,91 @@ bool PresolveFixed2dRectangles( } optional_boxes.erase(optional_boxes.begin(), optional_boxes.begin() + num_optional_boxes_to_remove); + return optional_boxes; +} + +} // namespace + +bool PresolveFixed2dRectangles( + absl::Span non_fixed_boxes, + std::vector* fixed_boxes) { + // This implementation compiles a set of areas that cannot be occupied by any + // item, then calls ReduceNumberofBoxes() to use these areas to minimize + // `fixed_boxes`. + bool changed = false; + + DCHECK(FindPartialRectangleIntersections(*fixed_boxes).empty()); + IntegerValue original_area = 0; + std::vector fixed_boxes_copy; + if (VLOG_IS_ON(1)) { + for (const Rectangle& r : *fixed_boxes) { + original_area += r.Area(); + } + } + if (VLOG_IS_ON(2)) { + fixed_boxes_copy = *fixed_boxes; + } + + const int original_num_boxes = fixed_boxes->size(); + + // The greedy algorithm is really fast. Run it first since it might greatly + // reduce the size of large trivial instances. + std::vector empty_vec; + if (ReduceNumberofBoxesGreedy(fixed_boxes, &empty_vec)) { + changed = true; + } + + IntegerValue min_x_size = std::numeric_limits::max(); + IntegerValue min_y_size = std::numeric_limits::max(); + + CHECK(!non_fixed_boxes.empty()); + Rectangle bounding_box = non_fixed_boxes[0].bounding_area; + + for (const RectangleInRange& box : non_fixed_boxes) { + bounding_box.GrowToInclude(box.bounding_area); + min_x_size = std::min(min_x_size, box.x_size); + min_y_size = std::min(min_y_size, box.y_size); + } + DCHECK_GT(min_x_size, 0); + DCHECK_GT(min_y_size, 0); + + // Fixed items are only useful to constraint where the non-fixed items can be + // placed. This means in particular that any part of a fixed item outside the + // bounding box of the non-fixed items is useless. Clip them. + int new_size = 0; + while (new_size < fixed_boxes->size()) { + Rectangle& rectangle = (*fixed_boxes)[new_size]; + DCHECK_GT(rectangle.SizeX(), 0); + DCHECK_GT(rectangle.SizeY(), 0); + if (rectangle.x_min < bounding_box.x_min) { + rectangle.x_min = bounding_box.x_min; + changed = true; + } + if (rectangle.x_max > bounding_box.x_max) { + rectangle.x_max = bounding_box.x_max; + changed = true; + } + if (rectangle.y_min < bounding_box.y_min) { + rectangle.y_min = bounding_box.y_min; + changed = true; + } + if (rectangle.y_max > bounding_box.y_max) { + rectangle.y_max = bounding_box.y_max; + changed = true; + } + if (rectangle.SizeX() <= 0 || rectangle.SizeY() <= 0) { + // The whole rectangle was outside of the domain, remove it. + std::swap(rectangle, (*fixed_boxes)[fixed_boxes->size() - 1]); + fixed_boxes->resize(fixed_boxes->size() - 1); + changed = true; + continue; + } else { + new_size++; + } + } + + std::vector optional_boxes = FindSpacesThatCannotBeOccupied( + *fixed_boxes, non_fixed_boxes, bounding_box, min_x_size, min_y_size); // TODO(user): instead of doing the greedy algorithm first with optional // boxes, and then the one that is exact for mandatory boxes but weak for @@ -1467,5 +1481,85 @@ bool ReduceNumberOfBoxesExactMandatory( return true; } +Disjoint2dPackingResult DetectDisjointRegionIn2dPacking( + absl::Span non_fixed_boxes, + absl::Span fixed_boxes, int max_num_components) { + if (max_num_components <= 1) return {}; + + IntegerValue min_x_size = std::numeric_limits::max(); + IntegerValue min_y_size = std::numeric_limits::max(); + + CHECK(!non_fixed_boxes.empty()); + Rectangle bounding_box = non_fixed_boxes[0].bounding_area; + + for (const RectangleInRange& box : non_fixed_boxes) { + bounding_box.GrowToInclude(box.bounding_area); + min_x_size = std::min(min_x_size, box.x_size); + min_y_size = std::min(min_y_size, box.y_size); + } + DCHECK_GT(min_x_size, 0); + DCHECK_GT(min_y_size, 0); + + std::vector optional_boxes = FindSpacesThatCannotBeOccupied( + fixed_boxes, non_fixed_boxes, bounding_box, min_x_size, min_y_size); + std::vector unoccupiable_space = {fixed_boxes.begin(), + fixed_boxes.end()}; + unoccupiable_space.insert(unoccupiable_space.end(), optional_boxes.begin(), + optional_boxes.end()); + + std::vector occupiable_space = + FindEmptySpaces(bounding_box, unoccupiable_space); + + std::vector empty; + ReduceNumberofBoxesGreedy(&occupiable_space, &empty); + std::vector> space_components = + SplitInConnectedComponents(BuildNeighboursGraph(occupiable_space)); + + if (space_components.size() == 1 || + space_components.size() > max_num_components) { + return {}; + } + + // If we are here, that means that the space where boxes can be placed is not + // connected. + Disjoint2dPackingResult result; + absl::flat_hash_set component_set; + for (const std::vector& component : space_components) { + Rectangle bin_bounding_box = occupiable_space[component[0]]; + for (int i = 1; i < component.size(); ++i) { + bin_bounding_box.GrowToInclude(occupiable_space[component[i]]); + } + std::optional reachable_area_bounding_box; + Disjoint2dPackingResult::Bin& bin = result.bins.emplace_back(); + for (int idx : component) { + bin.bin_area.push_back(occupiable_space[idx]); + } + for (int i = 0; i < non_fixed_boxes.size(); i++) { + if (!non_fixed_boxes[i].bounding_area.IsDisjoint(bin_bounding_box)) { + if (reachable_area_bounding_box.has_value()) { + reachable_area_bounding_box->GrowToInclude( + non_fixed_boxes[i].bounding_area); + } else { + reachable_area_bounding_box = non_fixed_boxes[i].bounding_area; + } + bin.non_fixed_box_indexes.push_back(i); + } + } + if (bin.non_fixed_box_indexes.empty()) { + result.bins.pop_back(); + continue; + } + bin.fixed_boxes = + FindEmptySpaces(*reachable_area_bounding_box, bin.bin_area); + ReduceNumberofBoxesGreedy(&bin.fixed_boxes, &empty); + } + VLOG_EVERY_N_SEC(1, 1) << "Detected a bin packing problem with " + << result.bins.size() + << " bins. Original problem sizes: " + << non_fixed_boxes.size() << " non-fixed boxes, " + << fixed_boxes.size() << " fixed boxes."; + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/2d_rectangle_presolve.h b/ortools/sat/2d_rectangle_presolve.h index a95fb17761..59fcebf541 100644 --- a/ortools/sat/2d_rectangle_presolve.h +++ b/ortools/sat/2d_rectangle_presolve.h @@ -38,6 +38,70 @@ bool PresolveFixed2dRectangles( absl::Span non_fixed_boxes, std::vector* fixed_boxes); +// Detect whether the fixed boxes of a no_overlap_2d constraint are splitting +// the space into separate components and thus can be replaced by one +// no_overlap_2d constraint per component. If this is not possible, return an +// empty result. Otherwise, return a struct containing what boxes (fixed and +// non-fixed) are needed in each new constraint. +// +// Note that for this to be correct, we need to introduce new boxes to "fill" +// the space occupied by the other components. For example, if we have a +// no_overlap_2d constraint with the fixed boxes and the bounding box of the +// non-fixed boxes as follows: +// +---------------------+ +// | | +// | | +// | | +// | | +// | | +// | ++++++++++++++++++| +// | ++++++++++++++++++| +// | ++++++++++++++++++| +// |***** | +// |***** | +// |***** | +// | | +// | | +// | | +// +---------------------+ +// and we are building the new no_overlap_2d constraint for the space below, we +// would need to add two new fixed boxes (the '.' and the 'o' one): +// +---------------------+ +// |ooooooooooooooooooooo| +// |ooooooooooooooooooooo| +// |ooooooooooooooooooooo| +// |ooooooooooooooooooooo| +// |ooooooooooooooooooooo| +// |...++++++++++++++++++| +// |...++++++++++++++++++| +// |...++++++++++++++++++| +// |***** | +// |***** | +// |***** | +// | | +// | | +// | | +// +---------------------+ +// This ensures that the new no_overlap_2d constraint will impose the box to be +// in that component as it should. Note that in the example above, the number of +// boxes won't really increase after a presolve pass, since the presolve for +// fixed boxes will simplify it to two fixed boxes again. +struct Disjoint2dPackingResult { + struct Bin { + // Fixed boxes that the non-fixed boxes in this bin cannot overlap with. + std::vector fixed_boxes; + // Non-fixed boxes on the original problem to copy to this new constraint. + std::vector non_fixed_box_indexes; + // Area that is covered by the connected component bin represents encoded as + // a non-overlapping set of rectangles. + std::vector bin_area; + }; + std::vector bins; +}; +Disjoint2dPackingResult DetectDisjointRegionIn2dPacking( + absl::Span non_fixed_boxes, + absl::Span fixed_boxes, int max_num_components); + // Given two vectors of non-overlapping rectangles defining two regions of the // space: one mandatory region that must be occupied and one optional region // that can be occupied, try to build a vector of as few non-overlapping diff --git a/ortools/sat/2d_rectangle_presolve_test.cc b/ortools/sat/2d_rectangle_presolve_test.cc index a10bb2d46b..54a08aa8fe 100644 --- a/ortools/sat/2d_rectangle_presolve_test.cc +++ b/ortools/sat/2d_rectangle_presolve_test.cc @@ -981,6 +981,93 @@ TEST(ReduceNumberOfBoxes, Problematic2) { CHECK(input == output); } +TEST(DetectDisjointRegionIn2dPackingTest, Basic) { + // Note that the hole of size 1 is ignored since no box is small enough to + // fit in it. + const std::vector fixed_rectangles = BuildFromAsciiArt(R"( +####################################### +####################################### +####################################### +########### ############## #### +########### ############## #### +########### ####### ### #### +########### ####### ### #### +########### ############## #### +########### ############## #### +########### ####################### +####################################### +####################################### +###################### ################ +#######################################)"); + + Rectangle bb = fixed_rectangles[0]; + for (const Rectangle& r : fixed_rectangles) { + bb.GrowToInclude(r); + } + const std::vector non_fixed_rectangles = { + {.box_index = 0, .bounding_area = bb, .x_size = 2, .y_size = 2}, + }; + + auto res = DetectDisjointRegionIn2dPacking(non_fixed_rectangles, + fixed_rectangles, 100); + EXPECT_THAT(res.bins.size(), 3); + std::sort(res.bins.begin(), res.bins.end(), + [](const Disjoint2dPackingResult::Bin& a, + const Disjoint2dPackingResult::Bin& b) { + return a.bin_area[0].SizeY() < b.bin_area[0].SizeY(); + }); + EXPECT_TRUE( + RectanglesCoverSameArea(res.bins[0].fixed_boxes, BuildFromAsciiArt(R"( +####################################### +####################################### +####################################### +####################################### +####################################### +####################### ############ +####################### ############ +####################################### +####################################### +####################################### +####################################### +####################################### +####################################### +#######################################)"))); + + EXPECT_TRUE( + RectanglesCoverSameArea(res.bins[1].fixed_boxes, BuildFromAsciiArt(R"( +####################################### +####################################### +####################################### +############################## #### +############################## #### +############################## #### +############################## #### +############################## #### +############################## #### +####################################### +####################################### +####################################### +####################################### +#######################################)"))); + + EXPECT_TRUE( + RectanglesCoverSameArea(res.bins[2].fixed_boxes, BuildFromAsciiArt(R"( +####################################### +####################################### +####################################### +########### ####################### +########### ####################### +########### ####################### +########### ####################### +########### ####################### +########### ####################### +########### ####################### +####################################### +####################################### +####################################### +#######################################)"))); +} + } // namespace } // namespace sat } // namespace operations_research diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index e436d225da..f32be59fc4 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -830,7 +830,6 @@ cc_library( ":sat_solver", ":solution_crush", ":util", - "//ortools/algorithms:sparse_permutation", "//ortools/base", "//ortools/port:proto_utils", "//ortools/util:affine_relation", @@ -881,11 +880,14 @@ cc_library( deps = [ ":cp_model_cc_proto", ":cp_model_utils", + ":diffn_util", ":sat_parameters_cc_proto", ":symmetry_util", + ":util", "//ortools/algorithms:sparse_permutation", "//ortools/util:sorted_interval_list", "@com_google_absl//absl/container:flat_hash_map", + "@com_google_absl//absl/container:flat_hash_set", "@com_google_absl//absl/container:inlined_vector", "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", @@ -2056,6 +2058,7 @@ cc_test( "//ortools/base", "//ortools/base:gmock_main", "//ortools/util:sorted_interval_list", + "@com_google_absl//absl/types:span", ], ) @@ -2659,9 +2662,11 @@ cc_library( ":model", ":precedences", ":sat_base", + ":util", "//ortools/base", "//ortools/base:mathutil", "//ortools/base:strong_vector", + "//ortools/graph", "//ortools/graph:max_flow", "//ortools/util:strong_integers", "@com_google_absl//absl/algorithm:container", @@ -2985,6 +2990,7 @@ cc_test( "//ortools/base", "//ortools/base:gmock_main", "//ortools/base:mathutil", + "//ortools/base:stl_util", "//ortools/util:random_engine", "//ortools/util:sorted_interval_list", "@com_google_absl//absl/container:btree", diff --git a/ortools/sat/CMakeLists.txt b/ortools/sat/CMakeLists.txt index a27ad417ce..e49cc257ce 100644 --- a/ortools/sat/CMakeLists.txt +++ b/ortools/sat/CMakeLists.txt @@ -11,7 +11,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -file(GLOB _SRCS "*.h" "*.cc" "python/linear_expr.*") +file(GLOB _SRCS "*.h" "*.cc" "python/linear_expr.*" "c_api/*.h" "c_api/*.cc") list(FILTER _SRCS EXCLUDE REGEX ".*/.*_test.cc") list(FILTER _SRCS EXCLUDE REGEX ".*/.*_fuzz.cc") list(REMOVE_ITEM _SRCS diff --git a/ortools/sat/all_different_test.cc b/ortools/sat/all_different_test.cc index 4ea2114a2f..a69afc5eb7 100644 --- a/ortools/sat/all_different_test.cc +++ b/ortools/sat/all_different_test.cc @@ -20,6 +20,7 @@ #include #include +#include "absl/types/span.h" #include "gtest/gtest.h" #include "ortools/base/logging.h" #include "ortools/sat/integer.h" @@ -36,8 +37,9 @@ namespace { class AllDifferentTest : public ::testing::TestWithParam { public: std::function AllDifferent( - const std::vector& vars) { - return [=](Model* model) { + absl::Span vars) { + return [=, vars = std::vector(vars.begin(), vars.end())]( + Model* model) { if (GetParam() == "binary") { model->Add(AllDifferentBinary(vars)); } else if (GetParam() == "ac") { diff --git a/ortools/sat/c_api/BUILD.bazel b/ortools/sat/c_api/BUILD.bazel new file mode 100644 index 0000000000..ce0b294d2f --- /dev/null +++ b/ortools/sat/c_api/BUILD.bazel @@ -0,0 +1,31 @@ +# Copyright 2010-2025 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. + +load("@rules_cc//cc:defs.bzl", "cc_library") + +package(default_visibility = ["//visibility:public"]) + +cc_library( + name = "cp_solver_c", + srcs = ["cp_solver_c.cc"], + hdrs = ["cp_solver_c.h"], + deps = [ + "//ortools/base:memutil", + "//ortools/sat:cp_model_cc_proto", + "//ortools/sat:cp_model_solver", + "//ortools/sat:model", + "//ortools/sat:sat_parameters_cc_proto", + "//ortools/util:time_limit", + "@com_google_absl//absl/log:check", + ], +) diff --git a/ortools/sat/c_api/cp_solver_c.cc b/ortools/sat/c_api/cp_solver_c.cc new file mode 100644 index 0000000000..17333d1b88 --- /dev/null +++ b/ortools/sat/c_api/cp_solver_c.cc @@ -0,0 +1,83 @@ +// Copyright 2010-2025 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/c_api/cp_solver_c.h" + +#include +#include + +#include "absl/log/check.h" +#include "ortools/base/memutil.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/cp_model_solver.h" +#include "ortools/sat/model.h" +#include "ortools/sat/sat_parameters.pb.h" +#include "ortools/util/time_limit.h" + +namespace operations_research::sat { + +namespace { + +CpSolverResponse solveWithParameters(std::atomic* const limit_reached, + const CpModelProto& proto, + const SatParameters& params) { + Model model; + model.Add(NewSatParameters(params)); + model.GetOrCreate()->RegisterExternalBooleanAsLimit(limit_reached); + return SolveCpModel(proto, &model); +} + +} // namespace + +extern "C" { + +void SolveCpModelWithParameters(const void* creq, int creq_len, + const void* cparams, int cparams_len, + void** cres, int* cres_len) { + return SolveCpInterruptible(nullptr, creq, creq_len, cparams, cparams_len, + cres, cres_len); +} + +void* SolveCpNewAtomicBool() { return new std::atomic(false); } + +void SolveCpDestroyAtomicBool(void* const atomic_bool) { + delete static_cast*>(atomic_bool); +} + +void SolveCpStopSolve(void* const atomic_bool) { + *static_cast*>(atomic_bool) = true; +} + +void SolveCpInterruptible(void* const limit_reached, const void* creq, + int creq_len, const void* cparams, int cparams_len, + void** cres, int* cres_len) { + CpModelProto req; + CHECK(req.ParseFromArray(creq, creq_len)); + + SatParameters params; + CHECK(params.ParseFromArray(cparams, cparams_len)); + + CpSolverResponse res = solveWithParameters( + static_cast*>(limit_reached), req, params); + + std::string res_str; + CHECK(res.SerializeToString(&res_str)); + + *cres_len = static_cast(res_str.size()); + *cres = strings::memdup(res_str.data(), *cres_len); + CHECK(*cres != nullptr); +} + +} // extern "C" + +} // namespace operations_research::sat diff --git a/ortools/sat/c_api/cp_solver_c.h b/ortools/sat/c_api/cp_solver_c.h new file mode 100644 index 0000000000..643cbe8753 --- /dev/null +++ b/ortools/sat/c_api/cp_solver_c.h @@ -0,0 +1,41 @@ +// Copyright 2010-2025 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_C_API_CP_SOLVER_C_H_ +#define OR_TOOLS_SAT_C_API_CP_SOLVER_C_H_ + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +void SolveCpModelWithParameters(const void* creq, int creq_len, + const void* cparams, int cparams_len, + void** cres, int* cres_len); + +void* SolveCpNewAtomicBool(); +void SolveCpDestroyAtomicBool(void* atomic_bool); +void SolveCpStopSolve(void* atomic_bool); +// Allows for interruptible solves. Solves can be interrupted by calling +// `SolveCpStopSolve` with the `limit_reached` atomic Boolean. +void SolveCpInterruptible(void* limit_reached, const void* creq, int creq_len, + const void* cparams, int cparams_len, void** cres, + int* cres_len); + +#ifdef __cplusplus +} // extern "C" +#endif + +#endif // OR_TOOLS_SAT_C_API_CP_SOLVER_C_H_ diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index df68ef5ecb..9ab3045ec1 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -869,8 +869,8 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { // Adds clauses so that literals[i] true <=> encoding[values[i]] true. // This also implicitly use the fact that exactly one alternative is true. -void LinkLiteralsAndValues(const std::vector& literals, - const std::vector& values, +void LinkLiteralsAndValues(absl::Span literals, + absl::Span values, const absl::flat_hash_map& encoding, PresolveContext* context) { CHECK_EQ(literals.size(), values.size()); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index cb1e991332..62cf2d3ffe 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -6005,6 +6005,141 @@ bool CpModelPresolver::PresolveNoOverlap2DFramed( return RemoveConstraint(ct); } +bool CpModelPresolver::ExpandEncoded2DBinPacking( + absl::Span fixed_boxes, + absl::Span non_fixed_boxes, ConstraintProto* ct) { + const Disjoint2dPackingResult disjoint_packing_presolve_result = + DetectDisjointRegionIn2dPacking( + non_fixed_boxes, fixed_boxes, + context_->params() + .maximum_regions_to_split_in_disconnected_no_overlap_2d()); + if (disjoint_packing_presolve_result.bins.empty()) return false; + + const NoOverlap2DConstraintProto& proto = ct->no_overlap_2d(); + std::vector box_in_area_lits; + absl::flat_hash_map> box_to_presence_literal; + // For the boxes that are optional, add a presence literal for each box in a + // fake "absent" bin. + for (int idx = 0; idx < non_fixed_boxes.size(); ++idx) { + const int b = non_fixed_boxes[idx].box_index; + const ConstraintProto& x_interval_ct = + context_->working_model->constraints(proto.x_intervals(b)); + const ConstraintProto& y_interval_ct = + context_->working_model->constraints(proto.y_intervals(b)); + if (x_interval_ct.enforcement_literal().empty() && + y_interval_ct.enforcement_literal().empty()) { + // Mandatory box, cannot be in the "absent" bin -1. + continue; + } + int enforcement_literal = x_interval_ct.enforcement_literal().empty() + ? y_interval_ct.enforcement_literal(0) + : x_interval_ct.enforcement_literal(0); + int potentially_other_enforcement_literal = + y_interval_ct.enforcement_literal().empty() + ? x_interval_ct.enforcement_literal(0) + : y_interval_ct.enforcement_literal(0); + + if (enforcement_literal == potentially_other_enforcement_literal) { + // The box is in the "absent" bin -1. + box_to_presence_literal[idx].push_back(NegatedRef(enforcement_literal)); + } else { + const int interval_is_absent_literal = + context_->NewBoolVarWithConjunction( + {enforcement_literal, potentially_other_enforcement_literal}); + + BoolArgumentProto* bool_or = + context_->working_model->add_constraints()->mutable_bool_or(); + bool_or->add_literals(NegatedRef(interval_is_absent_literal)); + for (const int lit : + {enforcement_literal, potentially_other_enforcement_literal}) { + context_->AddImplication(NegatedRef(interval_is_absent_literal), lit); + bool_or->add_literals(NegatedRef(lit)); + } + box_to_presence_literal[idx].push_back(interval_is_absent_literal); + } + } + // Now create the literals "item i in bin j". + for (int bin_index = 0; + bin_index < disjoint_packing_presolve_result.bins.size(); ++bin_index) { + const Disjoint2dPackingResult::Bin& bin = + disjoint_packing_presolve_result.bins[bin_index]; + NoOverlap2DConstraintProto new_no_overlap_2d; + for (const Rectangle& ret : bin.fixed_boxes) { + new_no_overlap_2d.add_x_intervals( + context_->working_model->constraints_size()); + new_no_overlap_2d.add_y_intervals( + context_->working_model->constraints_size() + 1); + IntervalConstraintProto* new_interval = + context_->working_model->add_constraints()->mutable_interval(); + new_interval->mutable_start()->set_offset(ret.x_min.value()); + new_interval->mutable_size()->set_offset(ret.SizeX().value()); + new_interval->mutable_end()->set_offset(ret.x_max.value()); + + new_interval = + context_->working_model->add_constraints()->mutable_interval(); + new_interval->mutable_start()->set_offset(ret.y_min.value()); + new_interval->mutable_size()->set_offset(ret.SizeY().value()); + new_interval->mutable_end()->set_offset(ret.y_max.value()); + } + for (const int idx : bin.non_fixed_box_indexes) { + int presence_in_box_lit = context_->NewBoolVar("binpacking"); + box_to_presence_literal[idx].push_back(presence_in_box_lit); + const int b = non_fixed_boxes[idx].box_index; + box_in_area_lits.push_back({.box_index = b, + .area_index = bin_index, + .literal = presence_in_box_lit}); + const ConstraintProto& x_interval_ct = + context_->working_model->constraints(proto.x_intervals(b)); + const ConstraintProto& y_interval_ct = + context_->working_model->constraints(proto.y_intervals(b)); + ConstraintProto* new_interval_x = + context_->working_model->add_constraints(); + *new_interval_x = x_interval_ct; + new_interval_x->clear_enforcement_literal(); + new_interval_x->add_enforcement_literal(presence_in_box_lit); + ConstraintProto* new_interval_y = + context_->working_model->add_constraints(); + *new_interval_y = y_interval_ct; + new_interval_y->clear_enforcement_literal(); + new_interval_y->add_enforcement_literal(presence_in_box_lit); + new_no_overlap_2d.add_x_intervals( + context_->working_model->constraints_size() - 2); + new_no_overlap_2d.add_y_intervals( + context_->working_model->constraints_size() - 1); + } + context_->working_model->add_constraints()->mutable_no_overlap_2d()->Swap( + &new_no_overlap_2d); + } + + // Each box is in exactly one bin (including the fake "absent" bin). + for (int box_index = 0; box_index < non_fixed_boxes.size(); ++box_index) { + const std::vector& presence_literals = + box_to_presence_literal[box_index]; + if (presence_literals.empty()) { + return context_->NotifyThatModelIsUnsat( + "A mandatory box cannot be placed in any position"); + } + auto* exactly_one = + context_->working_model->add_constraints()->mutable_exactly_one(); + for (const int presence_literal : presence_literals) { + exactly_one->add_literals(presence_literal); + } + } + CompactVectorVector areas; + for (int bin_index = 0; + bin_index < disjoint_packing_presolve_result.bins.size(); ++bin_index) { + areas.Add(disjoint_packing_presolve_result.bins[bin_index].bin_area); + } + solution_crush_.AssignVariableToPackingArea( + areas, *context_->working_model, proto.x_intervals(), proto.y_intervals(), + box_in_area_lits); + context_->UpdateNewConstraintsVariableUsage(); + context_->UpdateRuleStats( + "no_overlap_2d: fixed boxes partition available space, converted " + "to optional regions"); + return RemoveConstraint(ct); +} + bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) { if (context_->ModelIsUnsat()) { return false; @@ -6250,6 +6385,12 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) { "boxes"); return RemoveConstraint(ct); } + + if (!has_potential_zero_sized_interval) { + if (ExpandEncoded2DBinPacking(fixed_boxes, non_fixed_boxes, ct)) { + return true; + } + } RunPropagatorsForConstraint(*ct); return new_size < initial_num_boxes; } diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h index 04b7f54a49..fff33fe158 100644 --- a/ortools/sat/cp_model_presolve.h +++ b/ortools/sat/cp_model_presolve.h @@ -224,6 +224,19 @@ class CpModelPresolver { absl::Span fixed_boxes, absl::Span non_fixed_boxes, ConstraintProto* ct); + // Detects when the space where items of a no_overlap_2d constraint can placed + // is disjoint (ie., fixed boxes split the domain). When it is the case, we + // can introduce a boolean for each pair encoding whether + // the item is in the component or not. Then we replace the original + // no_overlap_2d constraint by one no_overlap_2d constraint for each + // component, with the new booleans as the enforcement_literal of the + // intervals. This is equivalent to expanding the original no_overlap_2d + // constraint into a bin packing problem with each connected component being a + // bin. + bool ExpandEncoded2DBinPacking( + absl::Span fixed_boxes, + absl::Span non_fixed_boxes, ConstraintProto* ct); + // SetPPC is short for set packing, partitioning and covering constraints. // These are sum of booleans <=, = and >= 1 respectively. // We detect inclusion of these constraint which allows a few simplifications. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 67d3ac6723..f0b124ed9c 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -1005,9 +1005,11 @@ bool SolutionHintIsCompleteAndFeasible( if (model_proto.has_objective()) { absl::StrAppend( &message, " Its objective value is ", - ScaleObjectiveValue( - model_proto.objective(), - ComputeInnerObjective(model_proto.objective(), solution)), + absl::StrFormat( + "%.9g", + ScaleObjectiveValue( + model_proto.objective(), + ComputeInnerObjective(model_proto.objective(), solution))), "."); } SOLVER_LOG(logger, message); diff --git a/ortools/sat/cp_model_solver_helpers.cc b/ortools/sat/cp_model_solver_helpers.cc index 5fd97b89da..1a035ab89f 100644 --- a/ortools/sat/cp_model_solver_helpers.cc +++ b/ortools/sat/cp_model_solver_helpers.cc @@ -1195,6 +1195,8 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { model->GetOrCreate()->ProcessImplicationGraph( model->GetOrCreate()); model->GetOrCreate()->Build(); + + model->GetOrCreate()->Build(); } void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) { @@ -1253,7 +1255,6 @@ void LoadCpModel(const CpModelProto& model_proto, Model* model) { // Note that we do that before we finish loading the problem (objective and // LP relaxation), because propagation will be faster at this point and it // should be enough for the purpose of this auto-detection. - model->GetOrCreate()->Build(); if (parameters.auto_detect_greater_than_at_least_one_of()) { model->GetOrCreate() ->AddGreaterThanAtLeastOneOfConstraints(model); diff --git a/ortools/sat/diffn_cuts.cc b/ortools/sat/diffn_cuts.cc index a2fb810f4f..f372ae3ef4 100644 --- a/ortools/sat/diffn_cuts.cc +++ b/ortools/sat/diffn_cuts.cc @@ -15,7 +15,6 @@ #include #include -#include #include #include #include @@ -33,7 +32,6 @@ #include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_base.h" -#include "ortools/sat/intervals.h" #include "ortools/sat/linear_constraint.h" #include "ortools/sat/linear_constraint_manager.h" #include "ortools/sat/model.h" @@ -393,17 +391,20 @@ std::string DiffnCtEvent::DebugString() const { // The original cut is: // sum(end_min_i * duration_min_i) >= // (sum(duration_min_i^2) + sum(duration_min_i)^2) / 2 -// We strengthen this cuts by noticing that if all tasks starts after S, -// then replacing end_min_i by (end_min_i - S) is still valid. -// -// 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). -// -// TODO(user): merge with Packing cuts. -void GenerateNoOvelap2dCompletionTimeCutsWithEnergy( - absl::string_view cut_name, std::vector events, - bool use_lifting, bool skip_low_sizes, Model* model, - LinearConstraintManager* manager) { +// We apply the following changes (see the code for cumulative constraints): +// - we strengthen this cuts by noticing that if all tasks starts after S, +// then replacing end_min_i by (end_min_i - S) is still valid. +// - we lift rectangles that start before the start of the sequence, but must +// overlap with it. +// - we apply the same transformation that was applied to the cumulative +// constraint to use the no_overlap cut in the no_overlap_2d setting. +// - we use a limited complexity subset-sum to compute reachable capacity +// - we look at a set of intervals starting after a given start_min, sorted by +// relative (end_lp - start_min). +void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name, + std::vector events, + bool use_lifting, Model* model, + LinearConstraintManager* manager) { TopNCuts top_n_cuts(5); // Sort by start min to bucketize by start_min. @@ -413,7 +414,7 @@ void GenerateNoOvelap2dCompletionTimeCutsWithEnergy( std::tie(e2.x_start_min, e2.y_size_min, e2.x_lp_end); }); for (int start = 0; start + 1 < events.size(); ++start) { - // Skip to the next start_min value. + // Skip to the next bucket (of start_min). if (start > 0 && events[start].x_start_min == events[start - 1].x_start_min) { continue; @@ -457,47 +458,54 @@ void GenerateNoOvelap2dCompletionTimeCutsWithEnergy( return e1.x_lp_end < e2.x_lp_end; }); + // Best cut so far for this loop. int best_end = -1; double best_efficacy = 0.01; - IntegerValue best_min_contrib(0); - IntegerValue sum_duration(0); - IntegerValue sum_square_duration(0); - IntegerValue best_capacity(0); - double unscaled_lp_contrib = 0.0; + IntegerValue best_min_contrib = 0; + IntegerValue best_capacity = 0; + IntegerValue best_y_range = 0; + + // Used in the first term of the rhs of the equation. + IntegerValue sum_event_contributions = 0; + // Used in the second term of the rhs of the equation. + IntegerValue sum_energy = 0; + // For normalization. + IntegerValue sum_square_energy = 0; + + double lp_contrib = 0.0; IntegerValue current_start_min(kMaxIntegerValue); - IntegerValue y_min = kMaxIntegerValue; - IntegerValue y_max = kMinIntegerValue; + IntegerValue y_min_of_subset = kMaxIntegerValue; + IntegerValue y_max_of_subset = kMinIntegerValue; + IntegerValue sum_of_y_size_min = 0; bool use_dp = true; MaxBoundedSubsetSum dp(0); + for (int i = 0; i < residual_tasks.size(); ++i) { const DiffnCtEvent& event = residual_tasks[i]; DCHECK_GE(event.x_start_min, sequence_start_min); - const IntegerValue energy = event.energy_min; - sum_duration += energy; - if (!AddProductTo(energy, energy, &sum_square_duration)) break; + // Make sure we do not overflow. + if (!AddTo(event.energy_min, &sum_energy)) break; + if (!AddProductTo(event.energy_min, event.x_size_min, + &sum_event_contributions)) { + break; + } + if (!AddSquareTo(event.energy_min, &sum_square_energy)) break; + if (!AddTo(event.y_size_min, &sum_of_y_size_min)) break; - unscaled_lp_contrib += event.x_lp_end * ToDouble(energy); + lp_contrib += event.x_lp_end * ToDouble(event.energy_min); current_start_min = std::min(current_start_min, event.x_start_min); - // This is competing with the brute force approach. Skip cases covered - // by the other code. - if (skip_low_sizes && i < 7) continue; - // For the capacity, we use the worse |y_max - y_min| and if all the tasks // so far have a fixed demand with a gcd > 1, we can round it down. - // - // TODO(user): Use dynamic programming to compute all possible values for - // the sum of demands as long as the involved numbers are small or the - // number of tasks are small. - y_min = std::min(y_min, event.y_min); - y_max = std::max(y_max, event.y_max); + y_min_of_subset = std::min(y_min_of_subset, event.y_min); + y_max_of_subset = std::max(y_max_of_subset, event.y_max); if (!event.y_size_is_fixed) use_dp = false; if (use_dp) { if (i == 0) { - dp.Reset((y_max - y_min).value()); + dp.Reset((y_max_of_subset - y_min_of_subset).value()); } else { - if (y_max - y_min != dp.Bound()) { + if (y_max_of_subset - y_min_of_subset != dp.Bound()) { use_dp = false; } } @@ -506,36 +514,45 @@ void GenerateNoOvelap2dCompletionTimeCutsWithEnergy( dp.Add(event.y_size_min.value()); } - const IntegerValue capacity = - use_dp ? IntegerValue(dp.CurrentMax()) : y_max - y_min; + const IntegerValue reachable_capacity = + use_dp ? IntegerValue(dp.CurrentMax()) + : y_max_of_subset - y_min_of_subset; - // We compute the cuts like if it was a disjunctive cut with all the - // duration actually equal to energy / capacity. But to keep the - // computation in the integer domain, we multiply by capacity - // everywhere instead. - IntegerValue min_contrib = 0; - if (!AddProductTo(sum_duration, sum_duration, &min_contrib)) break; - if (!AddTo(sum_square_duration, &min_contrib)) break; - min_contrib = min_contrib / 2; // The above is the double of the area. + // If we have not reached capacity, there can be no cuts on ends. + if (sum_of_y_size_min <= reachable_capacity) continue; - const IntegerValue intermediate = CapProdI(sum_duration, capacity); - if (AtMinOrMaxInt64I(intermediate)) break; - const IntegerValue offset = CapProdI(current_start_min, intermediate); - if (AtMinOrMaxInt64I(offset)) break; - if (!AddTo(offset, &min_contrib)) break; + // Do we have a violated cut ? + const IntegerValue large_rectangle_contrib = + CapProdI(sum_energy, sum_energy); + if (AtMinOrMaxInt64I(large_rectangle_contrib)) break; + const IntegerValue mean_large_rectangle_contrib = + CeilRatio(large_rectangle_contrib, reachable_capacity); - // We compute the efficacity in the unscaled domain where the l2 norm of - // the cuts is exactly the sqrt of the sum of squared duration. - const double efficacy = - (ToDouble(min_contrib) / ToDouble(capacity) - unscaled_lp_contrib) / - std::sqrt(ToDouble(sum_square_duration)); + IntegerValue min_contrib = + CapAddI(sum_event_contributions, mean_large_rectangle_contrib); + if (AtMinOrMaxInt64I(min_contrib)) break; + min_contrib = CeilRatio(min_contrib, 2); - // TODO(user): Check overflow and ignore if too big. + // shift contribution by current_start_min. + if (!AddProductTo(sum_energy, current_start_min, &min_contrib)) break; + + // The efficacy of the cut is the normalized violation of the above + // equation. We will normalize by the sqrt of the sum of squared energies. + const double efficacy = (ToDouble(min_contrib) - lp_contrib) / + std::sqrt(ToDouble(sum_square_energy)); + + // For a given start time, we only keep the best cut. + // The reason is that if the cut is strongly violated, we can get a + // sequence of violated cuts as we add more tasks. These new cuts will + // be less violated, but will not bring anything useful to the LP + // relaxation. At the same time, this sequence of cuts can push out + // other cuts from a disjoint set of tasks. if (efficacy > best_efficacy) { best_efficacy = efficacy; best_end = i; best_min_contrib = min_contrib; - best_capacity = capacity; + best_capacity = reachable_capacity; + best_y_range = y_max_of_subset - y_min_of_subset; } } if (best_end != -1) { @@ -546,18 +563,20 @@ void GenerateNoOvelap2dCompletionTimeCutsWithEnergy( const DiffnCtEvent& event = residual_tasks[i]; is_lifted |= event.lifted; add_energy_to_name |= event.use_energy; - cut.AddTerm(event.x_end, event.energy_min * best_capacity); + cut.AddTerm(event.x_end, event.energy_min); } std::string full_name(cut_name); if (is_lifted) full_name.append("_lifted"); if (add_energy_to_name) full_name.append("_energy"); + if (best_capacity < best_y_range) { + full_name.append("_subsetsum"); + } top_n_cuts.AddCut(cut.Build(), full_name, manager->LpValues()); } } top_n_cuts.TransferToManager(manager); } -// TODO(user): Use demands_helper and decomposed energy. CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( NoOverlap2DConstraintHelper* helper, Model* model) { CutGenerator result; @@ -616,10 +635,9 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( events.push_back(event); } - GenerateNoOvelap2dCompletionTimeCutsWithEnergy( - cut_name, std::move(events), - /*use_lifting=*/false, - /*skip_low_sizes=*/false, model, manager); + GenerateNoOvelap2dCompletionTimeCuts(cut_name, std::move(events), + /*use_lifting=*/true, model, + manager); }; if (!helper->SynchronizeAndSetDirection(true, true, false)) { @@ -633,11 +651,11 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( if (!helper->SynchronizeAndSetDirection(false, false, false)) { return false; } - generate_cuts("NoOverlap2dXCompletionTimeMirror"); + generate_cuts("NoOverlap2dXCompletionTime_mirror"); if (!helper->SynchronizeAndSetDirection(false, false, true)) { return false; } - generate_cuts("NoOverlap2dYCompletionTimeMirror"); + generate_cuts("NoOverlap2dYCompletionTime_mirror"); } return true; }; diff --git a/ortools/sat/diffn_cuts.h b/ortools/sat/diffn_cuts.h index ae31c3df3e..967a3e8fd7 100644 --- a/ortools/sat/diffn_cuts.h +++ b/ortools/sat/diffn_cuts.h @@ -17,7 +17,6 @@ #include #include -#include "absl/types/span.h" #include "ortools/sat/cuts.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_base.h" diff --git a/ortools/sat/go/cpmodel/BUILD.bazel b/ortools/sat/go/cpmodel/BUILD.bazel index 69f5a79f32..3b12cff8df 100644 --- a/ortools/sat/go/cpmodel/BUILD.bazel +++ b/ortools/sat/go/cpmodel/BUILD.bazel @@ -20,7 +20,7 @@ go_library( "cp_solver.go", "domain.go", ], - cdeps = [":cp_solver_c"], + cdeps = ["//ortools/sat/c_api:cp_solver_c"], cgo = True, importpath = "github.com/google/or-tools/ortools/sat/go/cpmodel", visibility = ["//visibility:public"], @@ -49,18 +49,3 @@ go_test( "@org_golang_google_protobuf//testing/protocmp", ], ) - -cc_library( - name = "cp_solver_c", - srcs = ["cp_solver_c.cc"], - hdrs = ["cp_solver_c.h"], - deps = [ - "//ortools/base:memutil", - "//ortools/sat:cp_model_cc_proto", - "//ortools/sat:cp_model_solver", - "//ortools/sat:model", - "//ortools/sat:sat_parameters_cc_proto", - "//ortools/util:time_limit", - "@com_google_absl//absl/log:check", - ], -) diff --git a/ortools/sat/go/cpmodel/cp_solver.go b/ortools/sat/go/cpmodel/cp_solver.go index 967827e4c5..2f9b887667 100644 --- a/ortools/sat/go/cpmodel/cp_solver.go +++ b/ortools/sat/go/cpmodel/cp_solver.go @@ -27,7 +27,7 @@ import ( /* #include // for free #include -#include "ortools/sat/go/cpmodel/cp_solver_c.h" +#include "ortools/sat/c_api/cp_solver_c.h" */ import "C" diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index cee9731ca1..5bfaf36fbc 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -1015,12 +1015,12 @@ void AppendNoOverlap2dRelaxationForComponent( if (max_area == kMaxIntegerValue) return; LinearConstraintBuilder lc(model, IntegerValue(0), max_area); - for (int i = 0; i < no_overlap_helper->NumBoxes(); ++i) { - if (no_overlap_helper->IsPresent(i)) { + for (const int b : component) { + if (no_overlap_helper->IsPresent(b)) { const AffineExpression& x_size_affine = - no_overlap_helper->x_helper().Sizes()[i]; + no_overlap_helper->x_helper().Sizes()[b]; const AffineExpression& y_size_affine = - no_overlap_helper->y_helper().Sizes()[i]; + no_overlap_helper->y_helper().Sizes()[b]; const std::vector energy = product_decomposer->TryToDecompose(x_size_affine, y_size_affine); if (!energy.empty()) { @@ -1028,17 +1028,17 @@ void AppendNoOverlap2dRelaxationForComponent( } else { lc.AddQuadraticLowerBound(x_size_affine, y_size_affine, integer_trail); } - } else if (no_overlap_helper->x_helper().IsPresent(i) || - no_overlap_helper->y_helper().IsPresent(i) || - (no_overlap_helper->x_helper().PresenceLiteral(i) == - no_overlap_helper->y_helper().PresenceLiteral(i))) { + } else if (no_overlap_helper->x_helper().IsPresent(b) || + no_overlap_helper->y_helper().IsPresent(b) || + (no_overlap_helper->x_helper().PresenceLiteral(b) == + no_overlap_helper->y_helper().PresenceLiteral(b))) { // We have only one active literal. const Literal presence_literal = - no_overlap_helper->x_helper().IsPresent(i) - ? no_overlap_helper->y_helper().PresenceLiteral(i) - : no_overlap_helper->x_helper().PresenceLiteral(i); + no_overlap_helper->x_helper().IsPresent(b) + ? no_overlap_helper->y_helper().PresenceLiteral(b) + : no_overlap_helper->x_helper().PresenceLiteral(b); const auto& [x_size, y_size] = - no_overlap_helper->GetLevelZeroBoxSizesMin(i); + no_overlap_helper->GetLevelZeroBoxSizesMin(b); const IntegerValue area_min = x_size * y_size; if (area_min > 0) { // Not including the term if we don't have a view is ok. diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index a144fd44a1..aa2533bf6b 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -28,6 +28,7 @@ #include "absl/log/check.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_format.h" +#include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" @@ -752,10 +753,11 @@ SatSolver::Status CoreBasedOptimizer::OptimizeWithSatEncoding( const int num_bools = sat_solver_->NumVariables(); const int num_fixed = sat_solver_->NumFixedVariables(); model_->GetOrCreate()->UpdateInnerObjectiveBounds( - absl::StrFormat( - "bool_core (num_cores=%d [%s] a=%u d=%d fixed=%d/%d clauses=%s)", - iter, previous_core_info, encoder.nodes().size(), max_depth, - num_fixed, num_bools, FormatCounter(clauses_->num_clauses())), + absl::StrFormat("bool_%s (num_cores=%d [%s] a=%u d=%d " + "fixed=%d/%d clauses=%s)", + model_->Name(), iter, previous_core_info, + encoder.nodes().size(), max_depth, num_fixed, + num_bools, FormatCounter(clauses_->num_clauses())), new_obj_lb, integer_trail_->LevelZeroUpperBound(objective_var_)); } diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 25f0d5a60f..0c1fd8ce97 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -11,6 +11,7 @@ // See the License for the specific language governing permissions and // limitations under the License. +// LINT: LEGACY_NAMES syntax = "proto2"; package operations_research.sat; @@ -23,7 +24,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 315 +// NEXT TAG: 316 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -902,6 +903,19 @@ message SatParameters { optional int32 max_pairs_pairwise_reasoning_in_no_overlap_2d = 276 [default = 1250]; + // Detects when the space where items of a no_overlap_2d constraint can placed + // is disjoint (ie., fixed boxes split the domain). When it is the case, we + // can introduce a boolean for each pair encoding whether + // the item is in the component or not. Then we replace the original + // no_overlap_2d constraint by one no_overlap_2d constraint for each + // component, with the new booleans as the enforcement_literal of the + // intervals. This is equivalent to expanding the original no_overlap_2d + // constraint into a bin packing problem with each connected component being a + // bin. This heuristic is only done when the number of regions to split + // is less than this parameter and <= 1 disables it. + optional int32 maximum_regions_to_split_in_disconnected_no_overlap_2d = 315 + [default = 0]; + // When set, it activates a few scheduling parameters to improve the lower // bound of scheduling problems. This is only effective with multiple workers // as it modifies the reduced_cost, lb_tree_search, and probing workers. diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 99bef978e1..20746e0df1 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -1316,22 +1316,6 @@ void AddEventDemandsToCapacitySubsetSum( } } -// In the no_overlap case, we have: -// area = event.x_size_min ^ 2 -// In the simple cumulative case, we split split the task (demand, size) into -// demand tasks in the no_overlap case. -// area = event.y_size_min * event.x_size_min * event.x_size_min -// In the cumulative case, the minimum energy of a task can be greater than -// the product of its size and demand (energy_min > side * demand). -// We use energy_min * size. -IntegerValue SmithRuleIndividualContribution(const CtEvent& event) { - if (event.use_energy) { - return event.energy_min * event.x_size_min; - } else { - return event.x_size_min * event.x_size_min * event.y_size_min; - } -} - } // namespace // We generate the cut from the Smith's rule from: @@ -1345,7 +1329,8 @@ IntegerValue SmithRuleIndividualContribution(const CtEvent& event) { // then replacing end_min_i by (end_min_i - S) is still valid. // // A second difference is that we lift intervals that starts before a given -// value, but are forced to cross it. +// value, but are forced to cross it. This lifting procedure implies trimming +// interval to its part that is after the given value. // // In the case of a cumulative constraint with a capacity of C, we compute a // valid equation by splitting the task (size_min si, demand_min di) into di @@ -1370,8 +1355,8 @@ IntegerValue SmithRuleIndividualContribution(const CtEvent& event) { // is the minimum size of the task, and di its minimum demand). We can // replace the small rectangle area term by ai * si. // -// sum (ai * ei) - sum (ai) * current_start_min >= -// sum(si * ai) / 2 + (sum (ai) ^ 2) / (2 * C) +// sum (ai * ei) - sum (ai) * current_start_min >= +// sum(si * ai) / 2 + (sum (ai) ^ 2) / (2 * C) // // The separation procedure is implemented using two loops: // - first, we loop on potential start times in increasing order. @@ -1451,8 +1436,15 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, // Make sure we do not overflow. if (!AddTo(event.energy_min, &sum_energy)) break; - if (!AddTo(SmithRuleIndividualContribution(event), - &sum_event_contributions)) { + // In the no_overlap case, we have: + // area = event.x_size_min ^ 2 + // In the simple cumulative case, we split split the task + // (demand_min, size_min) into demand_min tasks in the no_overlap case. + // area = event.y_size_min * event.x_size_min * event.x_size_min + // In the cumulative case, we can have energy_min > side_min * demand_min. + // In that case, we use energy_min * size_min. + if (!AddProductTo(event.energy_min, event.x_size_min, + &sum_event_contributions)) { break; } if (!AddSquareTo(event.energy_min, &sum_square_energy)) break; diff --git a/ortools/sat/solution_crush.cc b/ortools/sat/solution_crush.cc index 7681b9e6cb..b3ceae15ef 100644 --- a/ortools/sat/solution_crush.cc +++ b/ortools/sat/solution_crush.cc @@ -25,14 +25,17 @@ #include #include "absl/container/flat_hash_map.h" +#include "absl/container/flat_hash_set.h" #include "absl/log/check.h" #include "absl/types/span.h" #include "ortools/algorithms/sparse_permutation.h" #include "ortools/base/logging.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_utils.h" +#include "ortools/sat/diffn_util.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/symmetry_util.h" +#include "ortools/sat/util.h" #include "ortools/util/sorted_interval_list.h" namespace operations_research { @@ -622,5 +625,85 @@ void SolutionCrush::PermuteVariables(const SparsePermutation& permutation) { permutation.ApplyToDenseCollection(var_values_); } +void SolutionCrush::AssignVariableToPackingArea( + const CompactVectorVector& areas, const CpModelProto& model, + absl::Span x_intervals, absl::Span y_intervals, + absl::Span box_in_area_lits) { + struct RectangleTypeAndIndex { + enum class Type { + kHintedBox, + kArea, + }; + + int index; + Type type; + }; + std::vector rectangles_for_intersections; + std::vector rectangles_index; + + for (int i = 0; i < x_intervals.size(); ++i) { + const ConstraintProto& x_ct = model.constraints(x_intervals[i]); + const ConstraintProto& y_ct = model.constraints(y_intervals[i]); + + const std::optional x_min = + GetExpressionValue(x_ct.interval().start()); + const std::optional x_max = + GetExpressionValue(x_ct.interval().end()); + const std::optional y_min = + GetExpressionValue(y_ct.interval().start()); + const std::optional y_max = + GetExpressionValue(y_ct.interval().end()); + + if (!x_min.has_value() || !x_max.has_value() || !y_min.has_value() || + !y_max.has_value()) { + return; + } + if (*x_min > *x_max || *y_min > *y_max) { + VLOG(2) << "Hinted no_overlap_2d coordinate has max lower than min"; + return; + } + const Rectangle box = {.x_min = x_min.value(), + .x_max = x_max.value(), + .y_min = y_min.value(), + .y_max = y_max.value()}; + rectangles_for_intersections.push_back(box); + rectangles_index.push_back( + {.index = i, .type = RectangleTypeAndIndex::Type::kHintedBox}); + } + + for (int i = 0; i < areas.size(); ++i) { + for (const Rectangle& area : areas[i]) { + rectangles_for_intersections.push_back(area); + rectangles_index.push_back( + {.index = i, .type = RectangleTypeAndIndex::Type::kArea}); + } + } + + const std::vector> intersections = + FindPartialRectangleIntersections(rectangles_for_intersections); + + absl::flat_hash_set> box_to_area_pairs; + + for (const auto& [rec1_index, rec2_index] : intersections) { + RectangleTypeAndIndex rec1 = rectangles_index[rec1_index]; + RectangleTypeAndIndex rec2 = rectangles_index[rec2_index]; + if (rec1.type == rec2.type) { + DCHECK(rec1.type == RectangleTypeAndIndex::Type::kHintedBox); + VLOG(2) << "Hinted position of boxes in no_overlap_2d are overlapping"; + return; + } + if (rec1.type != RectangleTypeAndIndex::Type::kHintedBox) { + std::swap(rec1, rec2); + } + + box_to_area_pairs.insert({rec1.index, rec2.index}); + } + + for (const auto& [box_index, area_index, literal] : box_in_area_lits) { + SetLiteralValue(literal, + box_to_area_pairs.contains({box_index, area_index})); + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/solution_crush.h b/ortools/sat/solution_crush.h index f85511074d..34c10861cf 100644 --- a/ortools/sat/solution_crush.h +++ b/ortools/sat/solution_crush.h @@ -26,6 +26,8 @@ #include "ortools/algorithms/sparse_permutation.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_utils.h" +#include "ortools/sat/diffn_util.h" +#include "ortools/sat/util.h" #include "ortools/util/sorted_interval_list.h" namespace operations_research { @@ -278,6 +280,21 @@ class SolutionCrush { // Stores the solution as a hint in the given model. void StoreSolutionAsHint(CpModelProto& model) const; + // Given a list of N disjoint packing areas (each described by a union of + // rectangles) and a list of M boxes (described by their x and y interval + // constraints in the `model` proto), sets the value of the literals in + // `box_in_area_lits` with whether the box i intersects area j. + struct BoxInAreaLiteral { + int box_index; + int area_index; + int literal; + }; + void AssignVariableToPackingArea( + const CompactVectorVector& areas, + const CpModelProto& model, absl::Span x_intervals, + absl::Span y_intervals, + absl::Span box_in_area_lits); + private: bool HasValue(int var) const { return var_has_value_[var]; } diff --git a/ortools/util/fixed_shape_binary_tree.h b/ortools/util/fixed_shape_binary_tree.h index 86a5587e5a..17008245c5 100644 --- a/ortools/util/fixed_shape_binary_tree.h +++ b/ortools/util/fixed_shape_binary_tree.h @@ -218,6 +218,7 @@ class FixedShapeBinaryTree { template void PartitionIntervalIntoNodes(LeafIndex first_leaf, LeafIndex last_leaf, TypeWithPushBack* result) const { + DCHECK_LE(first_leaf, last_leaf); TreeNodeIndex prev(0); TreeNodeIndex current = GetNodeStartOfRange(first_leaf, last_leaf); if (current == Root()) {