14 #ifndef OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
15 #define OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
88 #include <type_traits>
92 #include "ortools/base/integral_types.h"
93 #include "ortools/base/logging.h"
94 #include "ortools/util/bitset.h"
95 #include "ortools/util/saturated_arithmetic.h"
96 #include "ortools/util/vector_or_function.h"
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) {
139 static_assert(std::is_integral<Integer>::value,
"Integral type required");
140 static_assert(std::is_unsigned<Integer>::value,
"Unsigned type required");
144 Integer
value()
const {
return value_; }
166 return (value_ & other.value_) == other.value_;
184 DCHECK(
Contains(n)) <<
"n = " << n <<
", value_ = " << value_;
194 return Set(value_ & (singleton.value_ - 1)).Cardinality();
211 return LeastSignificantBitPosition64(value_);
216 return BitCount64(value_);
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)) {
267 DCHECK_LT(0, 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>
358 DCHECK_LT(0, max_card);
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];
398 DCHECK_EQ(1 << n, sum);
400 DCHECK_EQ(0, base_offset_[1]);
401 DCHECK_EQ(max_card_ * (1 << (max_card_ - 1)),
402 base_offset_[max_card_] + max_card_);
406 template <
typename Set,
typename CostType>
411 uint64 local_offset = 0;
413 for (
int node : set) {
416 local_offset += binomial_coefficients_[node][node_rank + 1];
419 DCHECK_EQ(card, node_rank);
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> {
527 static int64 Add(int64 a, int64 b) {
return CapAdd(a, b); }
528 static int64 Sub(int64 a, int64 b) {
return CapSub(a, b); }
531 template <
bool Dummy>
532 struct SaturatedArithmetic<int32, Dummy> {
533 static int32 Add(int32 a, int32 b) {
536 const int64 min_int32 = kint32min;
537 const int64 max_int32 = kint32max;
538 return static_cast<int32
>(
539 std::max(min_int32, std::min(max_int32, a64 + b64)));
541 static int32 Sub(int32 a, int32 b) {
544 const int64 min_int32 = kint32min;
545 const int64 max_int32 = kint32max;
546 return static_cast<int32
>(
547 std::max(min_int32, std::min(max_int32, a64 - b64)));
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),
623 CHECK(cost_.Check());
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);
676 uint64 subset_offset =
677 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());
678 int prev_dest = set.SmallestElement();
680 for (
int dest : set) {
681 CostType min_cost = std::numeric_limits<CostType>::max();
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_);
714 CostType min_hamiltonian_cost = std::numeric_limits<CostType>::max();
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;
759 DCHECK_EQ(0, subset.value());
760 DCHECK(PathIsValid(path, 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_;
787 CostType min_cost = std::numeric_limits<CostType>::max();
788 CostType max_cost = std::numeric_limits<CostType>::min();
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);
915 MatrixOrFunction<CostType, CostFunction, true> cost_;
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) {
1020 #endif // OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_