diff --git a/ortools/set_cover/python/BUILD.bazel b/ortools/set_cover/python/BUILD.bazel index faad17f4672..722bad9e493 100644 --- a/ortools/set_cover/python/BUILD.bazel +++ b/ortools/set_cover/python/BUILD.bazel @@ -17,8 +17,6 @@ load("@pip_deps//:requirements.bzl", "requirement") load("@pybind11_bazel//:build_defs.bzl", "pybind_extension") load("@rules_python//python:py_test.bzl", "py_test") -package(default_visibility = ["//visibility:public"]) - pybind_extension( name = "set_cover", srcs = ["set_cover.cc"], @@ -27,6 +25,7 @@ pybind_extension( "//ortools/set_cover:base_types", "//ortools/set_cover:set_cover_heuristics", "//ortools/set_cover:set_cover_invariant", + "//ortools/set_cover:set_cover_lagrangian", "//ortools/set_cover:set_cover_model", "//ortools/set_cover:set_cover_reader", "@abseil-cpp//absl/types:span", diff --git a/ortools/set_cover/python/set_cover.cc b/ortools/set_cover/python/set_cover.cc index cdc34787615..cb9885cb5b1 100644 --- a/ortools/set_cover/python/set_cover.cc +++ b/ortools/set_cover/python/set_cover.cc @@ -23,6 +23,7 @@ #include "ortools/set_cover/base_types.h" #include "ortools/set_cover/set_cover_heuristics.h" #include "ortools/set_cover/set_cover_invariant.h" +#include "ortools/set_cover/set_cover_lagrangian.h" #include "ortools/set_cover/set_cover_model.h" #include "ortools/set_cover/set_cover_reader.h" #include "pybind11/numpy.h" @@ -32,7 +33,10 @@ #include "pybind11_protobuf/native_proto_caster.h" using ::operations_research::BaseInt; +using ::operations_research::ClearMostCoveredElements; using ::operations_research::ClearRandomSubsets; +using ::operations_research::Cost; +using ::operations_research::ElementCostVector; using ::operations_research::ElementDegreeSolutionGenerator; using ::operations_research::ElementIndex; using ::operations_research::GreedySolutionGenerator; @@ -54,6 +58,7 @@ using ::operations_research::WriteSetCoverSolutionText; using ::operations_research::SetCoverDecision; using ::operations_research::SetCoverInvariant; +using ::operations_research::SetCoverLagrangian; using ::operations_research::SetCoverModel; using ::operations_research::SparseColumn; using ::operations_research::SparseRow; @@ -71,7 +76,8 @@ using ::py::make_iterator; std::vector VectorIntToVectorSubsetIndex( absl::Span ints) { std::vector subs; - std::transform(ints.begin(), ints.end(), subs.begin(), + subs.reserve(ints.size()); + std::transform(ints.begin(), ints.end(), std::back_inserter(subs), [](int subset) -> SubsetIndex { return SubsetIndex(subset); }); return subs; } @@ -82,6 +88,12 @@ SubsetCostVector VectorDoubleToSubsetCostVector( return costs; } +ElementCostVector VectorDoubleToElementCostVector( + absl::Span doubles) { + ElementCostVector costs(doubles.begin(), doubles.end()); + return costs; +} + SubsetBoolVector BoolVectorToSubsetBoolVector(const std::vector& bools) { SubsetBoolVector bools_copy(bools.begin(), bools.end()); return bools_copy; @@ -199,9 +211,11 @@ PYBIND11_MODULE(set_cover, m) { .def_property_readonly("all_subsets", [](SetCoverModel& model) -> std::vector { std::vector subsets; + subsets.reserve(model.all_subsets().size()); std::transform( model.all_subsets().begin(), - model.all_subsets().end(), subsets.begin(), + model.all_subsets().end(), + std::back_inserter(subsets), [](const SubsetIndex element) -> BaseInt { return element.value(); }); @@ -555,9 +569,9 @@ PYBIND11_MODULE(set_cover, m) { return heuristic.NextSolution( BoolVectorToSubsetBoolVector(in_focus)); }) - .def("get_lagrangian_factor", &GuidedTabuSearch::SetLagrangianFactor, + .def("set_lagrangian_factor", &GuidedTabuSearch::SetLagrangianFactor, arg("factor")) - .def("set_lagrangian_factor", &GuidedTabuSearch::GetLagrangianFactor) + .def("get_lagrangian_factor", &GuidedTabuSearch::GetLagrangianFactor) .def("set_epsilon", &GuidedTabuSearch::SetEpsilon, arg("r")) .def("get_epsilon", &GuidedTabuSearch::GetEpsilon) .def("set_penalty_factor", &GuidedTabuSearch::SetPenaltyFactor, @@ -611,5 +625,131 @@ PYBIND11_MODULE(set_cover, m) { m.def("read_set_cover_solution_proto", &ReadSetCoverSolutionProto); // set_cover_lagrangian.h - // TODO(user): add support for SetCoverLagrangian. + py::class_(m, "SetCoverLagrangian") + .def(py::init()) + .def("use_num_threads", &SetCoverLagrangian::UseNumThreads, + arg("num_threads"), py::return_value_policy::reference_internal) + .def( + "initialize_lagrange_multipliers", + [](SetCoverLagrangian& lagrangian) -> std::vector { + return lagrangian.InitializeLagrangeMultipliers().get(); + }) + .def( + "compute_reduced_costs", + [](SetCoverLagrangian& lagrangian, + absl::Span costs, + absl::Span multipliers) -> std::vector { + return lagrangian + .ComputeReducedCosts(VectorDoubleToSubsetCostVector(costs), + VectorDoubleToElementCostVector(multipliers)) + .get(); + }, + arg("costs"), arg("multipliers")) + .def( + "parallel_compute_reduced_costs", + [](SetCoverLagrangian& lagrangian, + absl::Span costs, + absl::Span multipliers) -> std::vector { + return lagrangian + .ParallelComputeReducedCosts( + VectorDoubleToSubsetCostVector(costs), + VectorDoubleToElementCostVector(multipliers)) + .get(); + }, + arg("costs"), arg("multipliers")) + .def( + "compute_subgradient", + [](SetCoverLagrangian& lagrangian, + absl::Span reduced_costs) -> std::vector { + return lagrangian + .ComputeSubgradient( + VectorDoubleToSubsetCostVector(reduced_costs)) + .get(); + }, + arg("reduced_costs")) + .def( + "parallel_compute_subgradient", + [](SetCoverLagrangian& lagrangian, + absl::Span reduced_costs) -> std::vector { + return lagrangian + .ParallelComputeSubgradient( + VectorDoubleToSubsetCostVector(reduced_costs)) + .get(); + }, + arg("reduced_costs")) + .def( + "compute_lagrangian_value", + [](SetCoverLagrangian& lagrangian, + absl::Span reduced_costs, + absl::Span multipliers) -> double { + return lagrangian.ComputeLagrangianValue( + VectorDoubleToSubsetCostVector(reduced_costs), + VectorDoubleToElementCostVector(multipliers)); + }, + arg("reduced_costs"), arg("multipliers")) + .def( + "parallel_compute_lagrangian_value", + [](SetCoverLagrangian& lagrangian, + absl::Span reduced_costs, + absl::Span multipliers) -> double { + return lagrangian.ParallelComputeLagrangianValue( + VectorDoubleToSubsetCostVector(reduced_costs), + VectorDoubleToElementCostVector(multipliers)); + }, + arg("reduced_costs"), arg("multipliers")) + .def( + "update_multipliers", + [](SetCoverLagrangian& lagrangian, double step_size, + double lagrangian_value, double upper_bound, + absl::Span reduced_costs, + std::vector multipliers) -> std::vector { + ElementCostVector mult = + VectorDoubleToElementCostVector(multipliers); + lagrangian.UpdateMultipliers( + step_size, lagrangian_value, upper_bound, + VectorDoubleToSubsetCostVector(reduced_costs), &mult); + return mult.get(); + }, + arg("step_size"), arg("lagrangian_value"), arg("upper_bound"), + arg("reduced_costs"), arg("multipliers")) + .def( + "parallel_update_multipliers", + [](SetCoverLagrangian& lagrangian, double step_size, + double lagrangian_value, double upper_bound, + absl::Span reduced_costs, + std::vector multipliers) -> std::vector { + ElementCostVector mult = + VectorDoubleToElementCostVector(multipliers); + lagrangian.ParallelUpdateMultipliers( + step_size, lagrangian_value, upper_bound, + VectorDoubleToSubsetCostVector(reduced_costs), &mult); + return mult.get(); + }, + arg("step_size"), arg("lagrangian_value"), arg("upper_bound"), + arg("reduced_costs"), arg("multipliers")) + .def( + "compute_gap", + [](SetCoverLagrangian& lagrangian, + absl::Span reduced_costs, + const std::vector& solution, + absl::Span multipliers) -> double { + return lagrangian.ComputeGap( + VectorDoubleToSubsetCostVector(reduced_costs), + BoolVectorToSubsetBoolVector(solution), + VectorDoubleToElementCostVector(multipliers)); + }, + arg("reduced_costs"), arg("solution"), arg("multipliers")) + .def("three_phase", &SetCoverLagrangian::ThreePhase, arg("upper_bound")) + .def( + "compute_lower_bound", + [](SetCoverLagrangian& lagrangian, + absl::Span costs, + double upper_bound) + -> std::tuple, std::vector> { + auto [lb, reduced_costs, multipliers] = + lagrangian.ComputeLowerBound( + VectorDoubleToSubsetCostVector(costs), upper_bound); + return {lb, reduced_costs.get(), multipliers.get()}; + }, + arg("costs"), arg("upper_bound")); } diff --git a/ortools/set_cover/python/set_cover_test.py b/ortools/set_cover/python/set_cover_test.py index ee1ed4b9f54..7279e485d97 100644 --- a/ortools/set_cover/python/set_cover_test.py +++ b/ortools/set_cover/python/set_cover_test.py @@ -222,6 +222,125 @@ def test_knights_cover_trivial(self): inv.check_consistency(set_cover.consistency_level.FREE_AND_UNCOVERED) ) + def test_all_subsets_property(self): + model = create_knights_cover_model(4, 4) + all_subs = model.all_subsets + self.assertLen(all_subs, model.num_subsets) + self.assertEqual(all_subs, list(range(model.num_subsets))) + + def test_focus_based_next_solution(self): + model = create_knights_cover_model(8, 8) + self.assertTrue(model.compute_feasibility()) + inv = set_cover.SetCoverInvariant(model) + + focus = list(range(model.num_subsets)) + greedy = set_cover.GreedySolutionGenerator(inv) + self.assertTrue(greedy.next_solution(focus)) + self.assertEqual(inv.num_uncovered_elements(), 0) + + def test_guided_tabu_search_lagrangian_factor(self): + model = create_knights_cover_model(8, 8) + self.assertTrue(model.compute_feasibility()) + inv = set_cover.SetCoverInvariant(model) + + greedy = set_cover.GreedySolutionGenerator(inv) + self.assertTrue(greedy.next_solution()) + + gts = set_cover.GuidedTabuSearch(inv) + gts.initialize() + gts.set_lagrangian_factor(0.5) + self.assertAlmostEqual(gts.get_lagrangian_factor(), 0.5) + + def test_lagrangian_initialize_multipliers(self): + model = create_knights_cover_model(8, 8) + inv = set_cover.SetCoverInvariant(model) + + lagrangian = set_cover.SetCoverLagrangian(inv) + multipliers = lagrangian.initialize_lagrange_multipliers() + self.assertLen(multipliers, model.num_elements) + # All multipliers should be non-negative. + for m in multipliers: + self.assertGreaterEqual(m, 0.0) + + def test_lagrangian_compute_reduced_costs(self): + model = create_knights_cover_model(4, 4) + inv = set_cover.SetCoverInvariant(model) + + lagrangian = set_cover.SetCoverLagrangian(inv) + multipliers = lagrangian.initialize_lagrange_multipliers() + costs = list(model.subset_costs) + reduced = lagrangian.compute_reduced_costs(costs, multipliers) + self.assertLen(reduced, model.num_subsets) + + def test_lagrangian_compute_lower_bound(self): + model = create_knights_cover_model(8, 8) + inv = set_cover.SetCoverInvariant(model) + + # Get an upper bound from greedy. + greedy = set_cover.GreedySolutionGenerator(inv) + self.assertTrue(greedy.next_solution()) + upper_bound = inv.cost() + + lagrangian = set_cover.SetCoverLagrangian(inv) + costs = list(model.subset_costs) + lb, reduced_costs, multipliers = lagrangian.compute_lower_bound( + costs, upper_bound + ) + self.assertLen(reduced_costs, model.num_subsets) + self.assertLen(multipliers, model.num_elements) + self.assertLessEqual(lb, upper_bound) + + def test_lagrangian_subgradient_and_value(self): + model = create_knights_cover_model(4, 4) + inv = set_cover.SetCoverInvariant(model) + + lagrangian = set_cover.SetCoverLagrangian(inv) + multipliers = lagrangian.initialize_lagrange_multipliers() + costs = list(model.subset_costs) + reduced = lagrangian.compute_reduced_costs(costs, multipliers) + + subgradient = lagrangian.compute_subgradient(reduced) + self.assertLen(subgradient, model.num_elements) + + value = lagrangian.compute_lagrangian_value(reduced, multipliers) + self.assertIsInstance(value, float) + + def test_lagrangian_update_multipliers(self): + model = create_knights_cover_model(4, 4) + inv = set_cover.SetCoverInvariant(model) + + lagrangian = set_cover.SetCoverLagrangian(inv) + multipliers = lagrangian.initialize_lagrange_multipliers() + costs = list(model.subset_costs) + reduced = lagrangian.compute_reduced_costs(costs, multipliers) + value = lagrangian.compute_lagrangian_value(reduced, multipliers) + + updated = lagrangian.update_multipliers( + 1.0, value, 100.0, reduced, multipliers + ) + self.assertLen(updated, model.num_elements) + # Updated multipliers should still be non-negative. + for m in updated: + self.assertGreaterEqual(m, 0.0) + + def test_lagrangian_use_num_threads(self): + model = create_knights_cover_model(4, 4) + inv = set_cover.SetCoverInvariant(model) + + lagrangian = set_cover.SetCoverLagrangian(inv) + # Should return self for chaining. + result = lagrangian.use_num_threads(2) + self.assertIs(result, lagrangian) + + # Parallel methods should produce the same results as serial. + multipliers = lagrangian.initialize_lagrange_multipliers() + costs = list(model.subset_costs) + serial = lagrangian.compute_reduced_costs(costs, multipliers) + parallel = lagrangian.parallel_compute_reduced_costs(costs, multipliers) + self.assertEqual(len(serial), len(parallel)) + for s, p in zip(serial, parallel): + self.assertAlmostEqual(s, p) + # TODO(user): KnightsCoverGreedyAndTabu, KnightsCoverGreedyRandomClear, # KnightsCoverElementDegreeRandomClear, KnightsCoverRandomClearMip, # KnightsCoverMip