[ENH] set_cover Python: add SetCoverLagrangian bindings

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",

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<SubsetIndex> VectorIntToVectorSubsetIndex(
     absl::Span<const BaseInt> ints) {
   std::vector<SubsetIndex> 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<const double> doubles) {
+  ElementCostVector costs(doubles.begin(), doubles.end());
+  return costs;
+}
+
 SubsetBoolVector BoolVectorToSubsetBoolVector(const std::vector<bool>& 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<BaseInt> {
                                std::vector<BaseInt> 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_<SetCoverLagrangian>(m, "SetCoverLagrangian")
+      .def(py::init<SetCoverInvariant*>())
+      .def("use_num_threads", &SetCoverLagrangian::UseNumThreads,
+           arg("num_threads"), py::return_value_policy::reference_internal)
+      .def(
+          "initialize_lagrange_multipliers",
+          [](SetCoverLagrangian& lagrangian) -> std::vector<double> {
+            return lagrangian.InitializeLagrangeMultipliers().get();
+          })
+      .def(
+          "compute_reduced_costs",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> costs,
+             absl::Span<const double> multipliers) -> std::vector<double> {
+            return lagrangian
+                .ComputeReducedCosts(VectorDoubleToSubsetCostVector(costs),
+                                     VectorDoubleToElementCostVector(multipliers))
+                .get();
+          },
+          arg("costs"), arg("multipliers"))
+      .def(
+          "parallel_compute_reduced_costs",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> costs,
+             absl::Span<const double> multipliers) -> std::vector<double> {
+            return lagrangian
+                .ParallelComputeReducedCosts(
+                    VectorDoubleToSubsetCostVector(costs),
+                    VectorDoubleToElementCostVector(multipliers))
+                .get();
+          },
+          arg("costs"), arg("multipliers"))
+      .def(
+          "compute_subgradient",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> reduced_costs) -> std::vector<double> {
+            return lagrangian
+                .ComputeSubgradient(
+                    VectorDoubleToSubsetCostVector(reduced_costs))
+                .get();
+          },
+          arg("reduced_costs"))
+      .def(
+          "parallel_compute_subgradient",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> reduced_costs) -> std::vector<double> {
+            return lagrangian
+                .ParallelComputeSubgradient(
+                    VectorDoubleToSubsetCostVector(reduced_costs))
+                .get();
+          },
+          arg("reduced_costs"))
+      .def(
+          "compute_lagrangian_value",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> reduced_costs,
+             absl::Span<const double> multipliers) -> double {
+            return lagrangian.ComputeLagrangianValue(
+                VectorDoubleToSubsetCostVector(reduced_costs),
+                VectorDoubleToElementCostVector(multipliers));
+          },
+          arg("reduced_costs"), arg("multipliers"))
+      .def(
+          "parallel_compute_lagrangian_value",
+          [](SetCoverLagrangian& lagrangian,
+             absl::Span<const double> reduced_costs,
+             absl::Span<const double> 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<const double> reduced_costs,
+             std::vector<double> multipliers) -> std::vector<double> {
+            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<const double> reduced_costs,
+             std::vector<double> multipliers) -> std::vector<double> {
+            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<const double> reduced_costs,
+             const std::vector<bool>& solution,
+             absl::Span<const double> 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<const double> costs,
+             double upper_bound)
+              -> std::tuple<double, std::vector<double>, std::vector<double>> {
+            auto [lb, reduced_costs, multipliers] =
+                lagrangian.ComputeLowerBound(
+                    VectorDoubleToSubsetCostVector(costs), upper_bound);
+            return {lb, reduced_costs.get(), multipliers.get()};
+          },
+          arg("costs"), arg("upper_bound"));
 }

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