Column generation Bin Packing heuristic based on CFT

ortools/algorithms/bin_packing.h

@@ -0,0 +1,180 @@
+
+// Copyright 2025 Francesco Cavaliere
+// 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_ORTOOLS_ALGORITHMS_BIN_PACKING_H
+#define OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H
+
+#include <absl/algorithm/container.h>
+#include <absl/container/flat_hash_set.h>
+#include <absl/hash/hash.h>
+#include <absl/strings/str_split.h>
+#include <absl/strings/string_view.h>
+
+#include <random>
+#include <vector>
+
+#include "ortools/set_cover/base_types.h"
+#include "ortools/set_cover/set_cover_cft.h"
+
+namespace operations_research {
+
+class BinPackingModel {
+ public:
+  BinPackingModel() = default;
+  BaseInt num_items() const { return weigths_.size(); }
+  Cost bin_capacity() const { return bin_capcaity_; }
+  void set_bin_capacity(Cost capacity);
+  const ElementCostVector& weights() const { return weigths_; }
+  void AddItem(Cost weight);
+  void SortWeights();
+  ElementRange ItemRange() const {
+    return {ElementIndex(), ElementIndex(weigths_.size())};
+  }
+
+ private:
+  bool is_sorted_ = false;
+  Cost bin_capcaity_ = .0;
+  ElementCostVector weigths_ = {};
+};
+
+struct PartialBins {
+  std::vector<SparseColumn> bins;
+  std::vector<Cost> loads;
+};
+
+using SubsetHashVector = util_intops::StrongVector<SubsetIndex, uint64_t>;
+
+class ExpKnap {
+ public:
+  struct Item {
+    Cost profit;  // profit
+    Cost weight;  // weight
+    ElementIndex index;
+  };
+  using ItemIt = std::vector<Item>::const_reverse_iterator;
+
+  void SaveBin();
+  void FindGoodColumns(const ElementCostVector& profits,
+                       const ElementCostVector& weights, Cost capacity,
+                       BaseInt bnb_nodes_limit);
+  void InitSolver(const ElementCostVector& profits,
+                  const ElementCostVector& weights, Cost capacity,
+                  BaseInt bnb_nodes_limit);
+
+  bool EleBranch();
+  bool EleBranch(Cost profit_sum, Cost overweight, ItemIt out_item,
+                 ItemIt in_item);
+
+  void Heuristic();
+
+  const std::vector<SparseColumn>& collected_bins() const {
+    return collected_bins_;
+  }
+  Cost best_cost() const { return break_profit_sum_ + best_delta_; }
+
+ private:
+  Cost capacity_;                                        // capacity
+  util_intops::StrongVector<ElementIndex, Item> items_;  // items
+  ItemIt break_it_;
+  Cost break_profit_sum_;
+  Cost break_weight_sum_;
+  Cost best_delta_;
+  std::vector<ElementIndex> exceptions_;
+  std::vector<SparseColumn> collected_bins_;
+  ElementBoolVector break_selection_;
+  ElementBoolVector inserted_items_;
+  SparseColumn break_solution_;
+  BaseInt bnb_node_countdown_;
+  std::mt19937 rnd_;
+};
+
+class BinPackingSetCoverModel : public scp::FullToCoreModel {
+  using base = scp::FullToCoreModel;
+
+  struct BinPackingModelGlobals {
+    // Dirty hack to avoid invalidation of pointers/references
+    // A pointer to this data structure is used to compute the hash of bins
+    // starting from their indices.
+    scp::Model full_model;
+
+    // External bins do not have a valid index in the model, a temporary pointer
+    // is used insteda (even dirtier hack).
+    const SparseColumn* candidate_bin;
+
+    const SparseColumn& GetBin(SubsetIndex j) const;
+  };
+
+  struct BinHash {
+    const BinPackingModelGlobals* globals;
+    uint64_t operator()(SubsetIndex j) const;
+  };
+
+  struct BinEq {
+    const BinPackingModelGlobals* globals;
+    uint64_t operator()(SubsetIndex j1, SubsetIndex j2) const;
+  };
+
+ public:
+  BinPackingSetCoverModel(const BinPackingModel* bpp_model)
+      : globals_{scp::Model(), nullptr},
+        bpp_model_(bpp_model),
+        knapsack_solver_(),
+        bin_set_({}, 0, BinHash{&globals_}, BinEq{&globals_}),
+        column_gen_countdown_(10),
+        column_gen_period_(10) {}
+  const scp::Model& full_model() const { return globals_.full_model; }
+
+  bool AddBin(const SparseColumn& bin);
+  void CompleteModel() {
+    globals_.full_model.CreateSparseRowView();
+    static_cast<scp::FullToCoreModel&>(*this) =
+        scp::FullToCoreModel(&globals_.full_model);
+  }
+
+  bool UpdateCore(Cost best_lower_bound,
+                  const ElementCostVector& best_multipliers,
+                  const scp::Solution& best_solution, bool force) override;
+
+ private:
+  bool TryInsertBin(const SparseColumn& bin);
+
+  BinPackingModelGlobals globals_;
+  const BinPackingModel* bpp_model_;
+  ExpKnap knapsack_solver_;
+
+  // Contains bin indices, but it really should contains "bins" (aka,
+  // SparseColumn). However to avoid redundant allocations (in scp::Model and
+  // in the set) we cannot store them also here. We cannot also use
+  // iterators/pointers/references because they can be invalidated. So we store
+  // bin indices and do ungodly hacky shenanigans to get the bins from them.
+  absl::flat_hash_set<SubsetIndex, BinHash, BinEq> bin_set_;
+
+  Cost prev_lower_bound_;
+  BaseInt column_gen_countdown_;
+  BaseInt column_gen_period_;
+};
+
+BinPackingModel ReadBpp(absl::string_view filename);
+BinPackingModel ReadCsp(absl::string_view filename);
+
+void BestFit(const BinPackingModel& model,
+             const std::vector<ElementIndex>& items, PartialBins& bins_data);
+
+BinPackingSetCoverModel GenerateInitialBins(const BinPackingModel& model);
+
+void AddRandomizedBins(const BinPackingModel& model, BaseInt num_bins,
+                       BinPackingSetCoverModel& scp_model, std::mt19937& rnd);
+
+}  // namespace operations_research
+#endif /* OR_TOOLS_ORTOOLS_ALGORITHMS_BIN_PACKING_H */

ortools/algorithms/samples/bin_packing_cft.cc

@@ -0,0 +1,144 @@
+
+// Copyright 2025 Francesco Cavaliere
+// 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 <absl/log/initialize.h>
+#include <absl/status/status.h>
+
+#include <cstdlib>
+#include <random>
+
+#include "ortools/algorithms/bin_packing.h"
+#include "ortools/base/init_google.h"
+#include "ortools/set_cover/base_types.h"
+#include "ortools/set_cover/set_cover_cft.h"
+
+using namespace operations_research;
+ABSL_FLAG(std::string, instance, "", "BPP instance int RAIL format.");
+ABSL_FLAG(int, bins, 1000, "Number of bins to generate.");
+
+template <typename Iterable>
+std::string Stringify(const Iterable& col) {
+  std::string result;
+  for (auto i : col) {
+    absl::StrAppend(&result, " ", i);
+  }
+  return result;
+}
+
+bool operator==(const SparseColumn& lhs, const std::vector<BaseInt>& rhs) {
+  if (lhs.size() != rhs.size()) return false;
+  auto lit = lhs.begin();
+  auto rit = rhs.begin();
+  while (lit != lhs.end() && rit != rhs.end()) {
+    if (static_cast<BaseInt>(*lit) != *rit) {
+      return false;
+    }
+    ++lit;
+    ++rit;
+  }
+  return true;
+}
+
+void RunTest(const ElementCostVector& weights, const ElementCostVector& profits,
+             const std::vector<BaseInt>& expected) {
+  ExpKnap knap_solver;
+
+  for (ElementIndex i; i < ElementIndex(weights.size()); ++i) {
+    std::cout << "Item " << i << " -- profit: " << profits[i]
+              << " weight: " << weights[i]
+              << " efficiency: " << profits[i] / weights[i] << "\n";
+  }
+
+  knap_solver.InitSolver(profits, weights, 6, 100000000);
+  knap_solver.Heuristic();
+  std::cout << "Heur solution cost " << knap_solver.best_cost() << " -- "
+            << Stringify(knap_solver.collected_bins()[0]) << "\n";
+
+  knap_solver.EleBranch();
+  std::cout << "B&b solution cost " << knap_solver.best_cost() << " -- "
+            << Stringify(knap_solver.collected_bins()[0]) << "\n";
+
+  const auto& result = knap_solver.collected_bins()[0];
+  if (!(result == expected)) {
+    std::cout << "Error: expected " << Stringify(expected) << " but got "
+              << Stringify(result) << "\n";
+  }
+  std::cout << std::endl;
+}
+
+void KnapsackTest() {
+  std::cout << "Testing knapsack\n";
+  ExpKnap knap_solver;
+  ElementCostVector ws = {1, 2, 3, 4, 5};
+  RunTest(ws, {10, 20, 30, 40, 51}, {0, 4});
+  RunTest(ws, {10, 20, 30, 41, 50}, {1, 3});
+  RunTest(ws, {10, 20, 31, 40, 50}, {0, 1, 2});
+  RunTest(ws, {10, 21, 30, 41, 50}, {1, 3});
+  RunTest(ws, {11, 21, 30, 40, 50}, {0, 1, 2});
+  RunTest(ws, {11, 20, 31, 40, 50}, {0, 1, 2});
+  RunTest(ws, {11, 20, 30, 41, 50}, {0, 4});
+  RunTest(ws, {11, 20, 30, 40, 51}, {0, 4});
+  RunTest(ws, {11, 21, 31, 40, 50}, {0, 1, 2});
+  RunTest({4.1, 2, 2, 2}, {8.5, 3, 3, 3}, {1, 2, 3});
+}
+
+int main(int argc, char** argv) {
+  InitGoogle(argv[0], &argc, &argv, true);
+
+  // KnapsackTest();
+  // return 0;
+
+  BinPackingModel model = ReadBpp(absl::GetFlag(FLAGS_instance));
+
+  // Quick run with a minimal set of bins
+  BinPackingSetCoverModel scp_model = GenerateInitialBins(model);
+  scp::PrimalDualState best_result = scp::RunCftHeuristic(scp_model);
+
+  if (absl::GetFlag(FLAGS_bins) > 0) {
+    // Run the CFT again with more bins to get a better solution
+    std::mt19937 rnd(0);
+    AddRandomizedBins(model, absl::GetFlag(FLAGS_bins), scp_model, rnd);
+    scp::PrimalDualState result =
+        scp::RunCftHeuristic(scp_model, best_result.solution);
+    if (result.solution.cost() < best_result.solution.cost()) {
+      best_result = result;
+    }
+  }
+
+  auto [solution, dual] = best_result;
+  if (solution.subsets().empty()) {
+    std::cerr << "Error: failed to find any solution\n";
+  } else {
+    std::cout << "Solution:         " << solution.cost() << '\n';
+  }
+
+  if (dual.multipliers().empty()) {
+    std::cerr << "Error: failed to find any dual\n";
+  } else {
+    std::cout << "Core Lower bound: " << dual.lower_bound() << '\n';
+  }
+
+  // The lower bound computed on the full model is not a real lower bound unless
+  // the knapsack subproblem failed to fined any negative reduced cost bin to
+  // add to the set cover model.
+  // TODO(anyone): add a flag to indicate if a valid LB has been found or not.
+  if (scp_model.best_dual_state().multipliers().empty()) {
+    std::cerr << "Error: no real dual state has been computed\n";
+  } else {
+    std::cout << "Restricted Lower bound: "
+              << scp_model.best_dual_state().lower_bound() << '\n';
+  }
+
+  return EXIT_SUCCESS;
+}

ortools/set_cover/samples/cft.cc

@@ -0,0 +1,53 @@
+#include <absl/log/initialize.h>
+#include <cstdlib>
+#include <iostream>
+
+#include "ortools/base/init_google.h"
+#include "ortools/set_cover/set_cover_cft.h"
+#include "ortools/set_cover/set_cover_reader.h"
+
+using namespace operations_research;
+ABSL_FLAG(std::string, instance, "", "SCP instance int RAIL format.");
+
+#define DO_PRICING
+int main(int argc, char **argv) {
+  InitGoogle(argv[0], &argc, &argv, true);
+
+  scp::Model original_model = ReadOrlibRail(absl::GetFlag(FLAGS_instance));
+
+#ifdef DO_PRICING
+  scp::FullToCoreModel model(&original_model);
+#else
+  scp::SubModel model(&original_model);
+#endif
+
+  scp::PrimalDualState result = scp::RunCftHeuristic(model);
+  auto [solution, dual] = result;
+  if (solution.subsets().empty()) {
+    std::cerr << "Error: failed to find any solution\n";
+  } else {
+    std::cout << "Solution:         " << solution.cost() << '\n';
+  }
+
+#ifdef DO_PRICING
+  if (dual.multipliers().empty()) {
+    std::cerr << "Error: failed to find any dual\n";
+  } else {
+    std::cout << "Core Lower bound: " << dual.lower_bound() << '\n';
+  }
+  if (model.best_dual_state().multipliers().empty()) {
+    std::cerr << "Error: no real dual state has been computed\n";
+  } else {
+    std::cout << "Full Lower bound: " << model.best_dual_state().lower_bound()
+              << '\n';
+  }
+#else
+  if (dual.multipliers().empty()) {
+    std::cerr << "Error: failed to find any dual\n";
+  } else {
+    std::cout << "Lower bound: " << dual.lower_bound() << '\n';
+  }
+#endif
+
+  return EXIT_SUCCESS;
+}