14 #ifndef OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
15 #define OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
88 #include <type_traits>
103 template <
typename Set>
108 return current_set_ != other.current_set_;
126 template <
typename Integer>
133 static constexpr Integer
One =
static_cast<Integer
>(1);
134 static constexpr Integer
Zero =
static_cast<Integer
>(0);
138 explicit Set(Integer n) : value_(n) {
144 Integer
value()
const {
return value_; }
166 return (value_ & other.value_) == other.value_;
194 return Set(value_ & (singleton.value_ - 1)).Cardinality();
222 template <
typename SetRange>
234 return current_set_ != other.current_set_;
241 const IntegerType c = current_set_.SmallestSingleton().value();
246 const IntegerType shift = current_set_.SmallestElement();
247 current_set_ = r == 0 ?
SetType(0) :
SetType(((r ^
a) >> (shift + 2)) | r);
256 template <
typename Set>
264 : begin_(
Set::FullSet(card)),
265 end_(
Set::FullSet(card - 1).AddElement(max_card)) {
290 template <
typename Set,
typename CostType>
313 (binomial_coefficients_[added_node][rank] -
314 binomial_coefficients_[removed_node][rank]);
324 memory_[offset] =
value;
338 bool CheckConsistency()
const;
345 std::vector<std::vector<uint64>> binomial_coefficients_;
349 std::vector<int64> base_offset_;
353 std::vector<CostType> memory_;
356 template <
typename Set,
typename CostType>
360 if (max_card <= max_card_)
return;
361 max_card_ = max_card;
362 binomial_coefficients_.resize(max_card_ + 1);
365 for (
int n = 0; n <= max_card_; ++n) {
366 binomial_coefficients_[n].resize(n + 2);
367 binomial_coefficients_[n][0] = 1;
368 for (
int k = 1; k <= n; ++k) {
369 binomial_coefficients_[n][k] = binomial_coefficients_[n - 1][k - 1] +
370 binomial_coefficients_[n - 1][k];
374 binomial_coefficients_[n][n + 1] = 0;
376 base_offset_.resize(max_card_ + 1);
381 for (
int k = 0; k < max_card_; ++k) {
382 base_offset_[k + 1] =
383 base_offset_[k] + k * binomial_coefficients_[max_card_][k];
386 memory_.shrink_to_fit();
387 memory_.resize(max_card_ * (1 << (max_card_ - 1)));
388 DCHECK(CheckConsistency());
391 template <
typename Set,
typename CostType>
393 for (
int n = 0; n <= max_card_; ++n) {
395 for (
int k = 0; k <= n; ++k) {
396 sum += binomial_coefficients_[n][k];
401 DCHECK_EQ(max_card_ * (1 << (max_card_ - 1)),
402 base_offset_[max_card_] + max_card_);
406 template <
typename Set,
typename CostType>
413 for (
int node : set) {
416 local_offset += binomial_coefficients_[node][node_rank + 1];
427 return base_offset_[card] + card * local_offset;
430 template <
typename Set,
typename CostType>
436 template <
typename Set,
typename CostType>
439 return ValueAtOffset(Offset(set, node));
442 template <
typename Set,
typename CostType>
446 SetValueAtOffset(Offset(set, node),
value);
452 template <
typename CostType,
typename CostFunction>
517 template <
typename T,
521 struct SaturatedArithmetic {
522 static T Add(T
a, T
b) {
return a +
b; }
523 static T Sub(T
a, T
b) {
return a -
b; }
525 template <
bool Dummy>
526 struct SaturatedArithmetic<
int64, Dummy> {
531 template <
bool Dummy>
532 struct SaturatedArithmetic<
int32, Dummy> {
538 return static_cast<int32>(
546 return static_cast<int32>(
551 template <
typename T>
552 using Saturated = SaturatedArithmetic<T>;
555 CostType Cost(
int i,
int j) {
return cost_(i, j); }
561 std::vector<int> ComputePath(CostType
cost,
NodeSet set,
int end);
564 bool PathIsValid(
const std::vector<int>& path, CostType
cost);
567 MatrixOrFunction<CostType, CostFunction, true> cost_;
576 std::vector<CostType> hamiltonian_costs_;
579 bool triangle_inequality_ok_;
580 bool robustness_checked_;
581 bool triangle_inequality_checked_;
583 std::vector<int> tsp_path_;
587 std::vector<std::vector<int>> hamiltonian_paths_;
592 int best_hamiltonian_path_end_node_;
594 LatticeMemoryManager<NodeSet, CostType> mem_;
598 template <
typename CostType,
typename CostFunction>
600 int num_nodes, CostFunction
cost) {
605 template <
typename CostType,
typename CostFunction>
610 template <
typename CostType,
typename CostFunction>
612 int num_nodes, CostFunction
cost)
613 : cost_(std::move(
cost)),
614 num_nodes_(num_nodes),
616 hamiltonian_costs_(0),
618 triangle_inequality_ok_(true),
619 robustness_checked_(false),
620 triangle_inequality_checked_(false),
626 template <
typename CostType,
typename CostFunction>
629 ChangeCostMatrix(
cost.size(),
cost);
632 template <
typename CostType,
typename CostFunction>
634 int num_nodes, CostFunction
cost) {
635 robustness_checked_ =
false;
636 triangle_inequality_checked_ =
false;
639 num_nodes_ = num_nodes;
640 CHECK_GE(NodeSet::MaxCardinality, num_nodes_);
641 CHECK(cost_.Check());
644 template <
typename CostType,
typename CostFunction>
647 if (num_nodes_ == 0) {
650 hamiltonian_paths_.resize(1);
651 hamiltonian_costs_.resize(1);
652 best_hamiltonian_path_end_node_ = 0;
653 hamiltonian_costs_[0] = 0;
654 hamiltonian_paths_[0] = {0};
657 mem_.Init(num_nodes_);
660 for (
int dest = 0; dest < num_nodes_; ++dest) {
661 DCHECK_EQ(dest, mem_.BaseOffset(1, NodeSet::Singleton(dest)));
662 mem_.SetValueAtOffset(dest, Cost(0, dest));
667 for (
int card = 2; card <= num_nodes_; ++card) {
669 for (NodeSet set : SetRangeWithCardinality<Set<uint32>>(card, num_nodes_)) {
672 const uint64 set_offset = mem_.BaseOffset(card, set);
677 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());
678 int prev_dest = set.SmallestElement();
680 for (
int dest : set) {
682 const NodeSet subset = set.RemoveElement(dest);
686 subset_offset += mem_.OffsetDelta(card - 1, prev_dest, dest, dest_rank);
688 for (
int src : subset) {
690 min_cost, Saturated<CostType>::Add(
692 mem_.ValueAtOffset(subset_offset + src_rank)));
696 mem_.SetValueAtOffset(set_offset + dest_rank, min_cost);
702 const NodeSet full_set = NodeSet::FullSet(num_nodes_);
706 tsp_cost_ = mem_.Value(full_set, 0);
707 tsp_path_ = ComputePath(tsp_cost_, full_set, 0);
709 hamiltonian_paths_.resize(num_nodes_);
710 hamiltonian_costs_.resize(num_nodes_);
715 const NodeSet hamiltonian_set = full_set.RemoveElement(0);
716 for (
int end_node : hamiltonian_set) {
717 const CostType
cost = mem_.Value(hamiltonian_set, end_node);
718 hamiltonian_costs_[end_node] =
cost;
719 if (
cost <= min_hamiltonian_cost) {
720 min_hamiltonian_cost =
cost;
721 best_hamiltonian_path_end_node_ = end_node;
723 DCHECK_LE(tsp_cost_, Saturated<CostType>::Add(
cost, Cost(end_node, 0)));
725 hamiltonian_paths_[end_node] =
726 ComputePath(hamiltonian_costs_[end_node], hamiltonian_set, end_node);
732 template <
typename CostType,
typename CostFunction>
733 std::vector<int> HamiltonianPathSolver<CostType, CostFunction>::ComputePath(
734 CostType
cost, NodeSet set,
int end_node) {
735 DCHECK(set.Contains(end_node));
736 const int path_size = set.Cardinality() + 1;
737 std::vector<int> path(path_size, 0);
738 NodeSet subset = set.RemoveElement(end_node);
739 path[path_size - 1] = end_node;
741 CostType current_cost =
cost;
742 for (
int rank = path_size - 2; rank >= 0; --rank) {
743 for (
int src : subset) {
744 const CostType partial_cost = mem_.Value(subset, src);
745 const CostType incumbent_cost =
746 Saturated<CostType>::Add(partial_cost, Cost(src, dest));
749 if (std::abs(Saturated<CostType>::Sub(current_cost, incumbent_cost)) <=
750 std::numeric_limits<CostType>::epsilon() * current_cost) {
751 subset = subset.RemoveElement(src);
752 current_cost = partial_cost;
764 template <
typename CostType,
typename CostFunction>
765 bool HamiltonianPathSolver<CostType, CostFunction>::PathIsValid(
766 const std::vector<int>& path, CostType
cost) {
768 for (
int node : path) {
769 coverage = coverage.AddElement(node);
771 DCHECK_EQ(NodeSet::FullSet(num_nodes_).
value(), coverage.value());
772 CostType check_cost = 0;
773 for (
int i = 0; i < path.size() - 1; ++i) {
775 Saturated<CostType>::Add(check_cost, Cost(path[i], path[i + 1]));
777 DCHECK_LE(std::abs(Saturated<CostType>::Sub(
cost, check_cost)),
778 std::numeric_limits<CostType>::epsilon() *
cost)
779 <<
"cost = " <<
cost <<
" check_cost = " << check_cost;
783 template <
typename CostType,
typename CostFunction>
785 if (std::numeric_limits<CostType>::is_integer)
return true;
786 if (robustness_checked_)
return robust_;
791 for (
int i = 0; i < num_nodes_; ++i) {
792 for (
int j = 0; j < num_nodes_; ++j) {
793 if (i == j)
continue;
794 min_cost =
std::min(min_cost, Cost(i, j));
795 max_cost =
std::max(max_cost, Cost(i, j));
801 min_cost >= 0 && min_cost > num_nodes_ * max_cost *
802 std::numeric_limits<CostType>::epsilon();
803 robustness_checked_ =
true;
807 template <
typename CostType,
typename CostFunction>
809 CostFunction>::VerifiesTriangleInequality() {
810 if (triangle_inequality_checked_)
return triangle_inequality_ok_;
811 triangle_inequality_ok_ =
true;
812 triangle_inequality_checked_ =
true;
813 for (
int k = 0; k < num_nodes_; ++k) {
814 for (
int i = 0; i < num_nodes_; ++i) {
815 for (
int j = 0; j < num_nodes_; ++j) {
816 const CostType detour_cost =
817 Saturated<CostType>::Add(Cost(i, k), Cost(k, j));
818 if (detour_cost < Cost(i, j)) {
819 triangle_inequality_ok_ =
false;
820 return triangle_inequality_ok_;
825 return triangle_inequality_ok_;
828 template <
typename CostType,
typename CostFunction>
830 CostFunction>::BestHamiltonianPathEndNode() {
832 return best_hamiltonian_path_end_node_;
835 template <
typename CostType,
typename CostFunction>
839 return hamiltonian_costs_[end_node];
842 template <
typename CostType,
typename CostFunction>
846 return hamiltonian_paths_[end_node];
849 template <
typename CostType,
typename CostFunction>
851 std::vector<PathNodeIndex>* path) {
852 *path = HamiltonianPath(best_hamiltonian_path_end_node_);
855 template <
typename CostType,
typename CostFunction>
862 template <
typename CostType,
typename CostFunction>
869 template <
typename CostType,
typename CostFunction>
871 std::vector<PathNodeIndex>* path) {
872 *path = TravelingSalesmanPath();
875 template <
typename CostType,
typename CostFunction>
906 CostType Cost(
int i,
int j) {
return cost_(i, j); }
909 void Solve(
int end_node);
912 CostType ComputeFutureLowerBound(
NodeSet current_set,
int last_visited);
930 template <
typename CostType,
typename CostFunction>
935 template <
typename CostType,
typename CostFunction>
937 int num_nodes, CostFunction
cost)
938 : cost_(std::move(
cost)),
939 num_nodes_(num_nodes),
943 template <
typename CostType,
typename CostFunction>
945 if (solved_ || num_nodes_ == 0)
return;
951 mem_.Init(num_nodes_);
952 NodeSet start_set = NodeSet::Singleton(0);
953 std::stack<std::pair<NodeSet, int>> state_stack;
954 state_stack.push(std::make_pair(start_set, 0));
956 while (!state_stack.empty()) {
957 const std::pair<NodeSet, int> current = state_stack.top();
960 const NodeSet current_set = current.first;
961 const int last_visited = current.second;
962 const CostType current_cost = mem_.Value(current_set, last_visited);
965 for (
int next_to_visit = 0; next_to_visit < num_nodes_; next_to_visit++) {
969 if (current_set.Contains(next_to_visit))
continue;
972 const int next_cardinality = current_set.Cardinality() + 1;
973 if (next_to_visit == end_node && next_cardinality != num_nodes_)
continue;
975 const NodeSet next_set = current_set.AddElement(next_to_visit);
976 const CostType next_cost =
977 current_cost + Cost(last_visited, next_to_visit);
980 const CostType previous_best = mem_.Value(next_set, next_to_visit);
981 if (previous_best != 0 && next_cost >= previous_best)
continue;
985 const CostType lower_bound =
986 ComputeFutureLowerBound(next_set, next_to_visit);
987 if (tsp_cost_ != 0 && next_cost + lower_bound >= tsp_cost_)
continue;
990 if (next_cardinality == num_nodes_) {
991 tsp_cost_ = next_cost;
996 mem_.SetValue(next_set, next_to_visit, next_cost);
997 state_stack.push(std::make_pair(next_set, next_to_visit));
1004 template <
typename CostType,
typename CostFunction>
1011 template <
typename CostType,
typename CostFunction>
1014 NodeSet current_set,
int last_visited) {
#define DCHECK_LE(val1, val2)
#define CHECK_GE(val1, val2)
#define DCHECK_GE(val1, val2)
#define DCHECK_LT(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
const ElementIterator & operator++()
bool operator!=(const ElementIterator &other) const
CostType HamiltonianCost(int end_node)
void ChangeCostMatrix(CostFunction cost)
void ChangeCostMatrix(int num_nodes, CostFunction cost)
void TravelingSalesmanPath(std::vector< PathNodeIndex > *path)
std::vector< int > TravelingSalesmanPath()
bool VerifiesTriangleInequality()
std::vector< int > HamiltonianPath(int end_node)
void HamiltonianPath(std::vector< PathNodeIndex > *path)
HamiltonianPathSolver(CostFunction cost)
CostType TravelingSalesmanCost()
HamiltonianPathSolver(int num_nodes, CostFunction cost)
int BestHamiltonianPathEndNode()
void SetValue(Set s, int node, CostType value)
CostType ValueAtOffset(uint64 offset) const
uint64 Offset(Set s, int node) const
void SetValueAtOffset(uint64 offset, CostType value)
uint64 BaseOffset(int card, Set s) const
CostType Value(Set s, int node) const
uint64 OffsetDelta(int card, int added_node, int removed_node, int rank) const
CostType HamiltonianCost(int end_node)
PruningHamiltonianSolver(CostFunction cost)
static const int MaxCardinality
bool Contains(int n) const
Set RemoveElement(int n) const
ElementIterator< Set > begin() const
int SingletonRank(Set singleton) const
int ElementRank(int n) const
bool operator!=(const Set &other) const
Set RemoveSmallestElement() const
Set SmallestSingleton() const
ElementIterator< Set > end() const
static constexpr Integer One
bool Includes(Set other) const
int SmallestElement() const
static Set Singleton(Integer n)
static constexpr Integer Zero
static Set FullSet(Integer card)
Set AddElement(int n) const
SetRange::SetType SetType
SetType::IntegerType IntegerType
bool operator!=(const SetRangeIterator &other) const
const SetRangeIterator & operator++()
SetRangeIterator(const SetType set)
SetType operator*() const
SetRangeIterator< SetRangeWithCardinality > begin() const
SetRangeWithCardinality(int card, int max_card)
SetRangeIterator< SetRangeWithCardinality > end() const
static const int32 kint32max
static const int32 kint32min
CpSolverResponse Solve(const CpModelProto &model_proto)
Solves the given CpModelProto and returns an instance of CpSolverResponse.
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int LeastSignificantBitPosition64(uint64 n)
uint32 BitCount32(uint32 n)
int64 CapSub(int64 x, int64 y)
uint64 BitCount64(uint64 n)
HamiltonianPathSolver< CostType, CostFunction > MakeHamiltonianPathSolver(int num_nodes, CostFunction cost)
int LeastSignificantBitPosition32(uint32 n)
int64 CapAdd(int64 x, int64 y)