522 template <
typename T,
526 struct SaturatedArithmetic {
527 static T Add(T a, T b) {
return a + b; }
528 static T Sub(T a, T b) {
return a - b; }
530 template <
bool Dummy>
531 struct SaturatedArithmetic<int64_t, Dummy> {
532 static int64_t Add(int64_t a, int64_t b) {
return CapAdd(a, b); }
533 static int64_t Sub(int64_t a, int64_t b) {
return CapSub(a, b); }
536 template <
bool Dummy>
537 struct SaturatedArithmetic<int32_t, Dummy> {
538 static int32_t Add(int32_t a, int32_t b) {
539 const int64_t a64 = a;
540 const int64_t b64 =
b;
541 const int64_t min_int32 = std::numeric_limits<int32_t>::min();
542 const int64_t max_int32 = std::numeric_limits<int32_t>::max();
543 return static_cast<int32_t
>(
544 std::max(min_int32, std::min(max_int32, a64 + b64)));
546 static int32_t Sub(int32_t a, int32_t b) {
547 const int64_t a64 = a;
548 const int64_t b64 =
b;
549 const int64_t min_int32 = std::numeric_limits<int32_t>::min();
550 const int64_t max_int32 = std::numeric_limits<int32_t>::max();
551 return static_cast<int32_t
>(
552 std::max(min_int32, std::min(max_int32, a64 - b64)));
556 template <
typename T>
557 using Saturated = SaturatedArithmetic<T>;
560 CostType
Cost(
int i,
int j) {
return cost_(i, j); }
566 std::vector<int> ComputePath(CostType cost,
NodeSet set,
int end);
569 bool PathIsValid(absl::Span<const int> path, CostType cost);
572 MatrixOrFunction<CostType, CostFunction, true> cost_;
581 std::vector<CostType> hamiltonian_costs_;
584 bool triangle_inequality_ok_;
585 bool robustness_checked_;
586 bool triangle_inequality_checked_;
588 std::vector<int> tsp_path_;
592 std::vector<std::vector<int>> hamiltonian_paths_;
597 int best_hamiltonian_path_end_node_;
599 LatticeMemoryManager<NodeSet, CostType> mem_;
603template <
typename CostType,
typename CostFunction>
605 int num_nodes, CostFunction cost) {
610template <
typename CostType,
typename CostFunction>
615template <
typename CostType,
typename CostFunction>
617 int num_nodes, CostFunction cost)
618 : cost_(
std::move(cost)),
619 num_nodes_(num_nodes),
621 hamiltonian_costs_(0),
623 triangle_inequality_ok_(true),
624 robustness_checked_(false),
625 triangle_inequality_checked_(false),
628 CHECK(cost_.Check());
631template <
typename CostType,
typename CostFunction>
637template <
typename CostType,
typename CostFunction>
639 int num_nodes, CostFunction cost) {
640 robustness_checked_ =
false;
641 triangle_inequality_checked_ =
false;
644 num_nodes_ = num_nodes;
646 CHECK(cost_.Check());
649template <
typename CostType,
typename CostFunction>
650void HamiltonianPathSolver<CostType, CostFunction>::Solve() {
652 if (num_nodes_ == 0) {
655 hamiltonian_paths_.resize(1);
656 hamiltonian_costs_.resize(1);
657 best_hamiltonian_path_end_node_ = 0;
658 hamiltonian_costs_[0] = 0;
659 hamiltonian_paths_[0] = {0};
662 mem_.Init(num_nodes_);
665 for (
int dest = 0; dest < num_nodes_; ++dest) {
666 DCHECK_EQ(dest, mem_.BaseOffset(1, NodeSet::Singleton(dest)));
667 mem_.SetValueAtOffset(dest,
Cost(0, dest));
672 for (
int card = 2; card <= num_nodes_; ++card) {
678 const uint64_t set_offset = mem_.BaseOffset(card, set);
682 uint64_t subset_offset =
683 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());
684 int prev_dest = set.SmallestElement();
686 for (
int dest : set) {
687 CostType min_cost = std::numeric_limits<CostType>::max();
688 const NodeSet subset = set.RemoveElement(dest);
692 subset_offset += mem_.OffsetDelta(card - 1, prev_dest, dest, dest_rank);
694 for (
int src : subset) {
696 min_cost, Saturated<CostType>::Add(
698 mem_.ValueAtOffset(subset_offset + src_rank)));
702 mem_.SetValueAtOffset(set_offset + dest_rank, min_cost);
708 const NodeSet full_set = NodeSet::FullSet(num_nodes_);
712 tsp_cost_ = mem_.Value(full_set, 0);
713 tsp_path_ = ComputePath(tsp_cost_, full_set, 0);
715 hamiltonian_paths_.resize(num_nodes_);
716 hamiltonian_costs_.resize(num_nodes_);
720 CostType min_hamiltonian_cost = std::numeric_limits<CostType>::max();
721 const NodeSet hamiltonian_set = full_set.RemoveElement(0);
722 for (
int end_node : hamiltonian_set) {
723 const CostType cost = mem_.Value(hamiltonian_set, end_node);
724 hamiltonian_costs_[end_node] = cost;
725 if (cost <= min_hamiltonian_cost) {
726 min_hamiltonian_cost = cost;
727 best_hamiltonian_path_end_node_ = end_node;
729 DCHECK_LE(tsp_cost_, Saturated<CostType>::Add(cost,
Cost(end_node, 0)));
731 hamiltonian_paths_[end_node] =
732 ComputePath(hamiltonian_costs_[end_node], hamiltonian_set, end_node);
738template <
typename CostType,
typename CostFunction>
739std::vector<int> HamiltonianPathSolver<CostType, CostFunction>::ComputePath(
740 CostType cost, NodeSet set,
int end_node) {
741 DCHECK(set.Contains(end_node));
742 const int path_size = set.Cardinality() + 1;
743 std::vector<int> path(path_size, 0);
744 NodeSet subset = set.RemoveElement(end_node);
745 path[path_size - 1] = end_node;
747 CostType current_cost = cost;
748 for (
int rank = path_size - 2; rank >= 0; --rank) {
749 for (
int src : subset) {
750 const CostType partial_cost = mem_.Value(subset, src);
751 const CostType incumbent_cost =
752 Saturated<CostType>::Add(partial_cost,
Cost(src, dest));
755 if (std::abs(Saturated<CostType>::Sub(current_cost, incumbent_cost)) <=
756 std::numeric_limits<CostType>::epsilon() * current_cost) {
757 subset = subset.RemoveElement(src);
758 current_cost = partial_cost;
765 DCHECK_EQ(0, subset.value());
766 DCHECK(PathIsValid(path, cost));
770template <
typename CostType,
typename CostFunction>
771bool HamiltonianPathSolver<CostType, CostFunction>::PathIsValid(
772 absl::Span<const int> path, CostType cost) {
774 for (
int node : path) {
775 coverage = coverage.AddElement(node);
777 DCHECK_EQ(NodeSet::FullSet(num_nodes_).value(), coverage.value());
778 CostType check_cost = 0;
779 for (
int i = 0;
i < path.size() - 1; ++
i) {
781 Saturated<CostType>::Add(check_cost,
Cost(path[i], path[i + 1]));
783 DCHECK_LE(std::abs(Saturated<CostType>::Sub(cost, check_cost)),
784 std::numeric_limits<CostType>::epsilon() * cost)
785 <<
"cost = " << cost <<
" check_cost = " << check_cost;
789template <
typename CostType,
typename CostFunction>
791 if (std::numeric_limits<CostType>::is_integer)
return true;
792 if (robustness_checked_)
return robust_;
793 CostType min_cost = std::numeric_limits<CostType>::max();
794 CostType max_cost = std::numeric_limits<CostType>::min();
797 for (
int i = 0; i < num_nodes_; ++i) {
798 for (
int j = 0; j < num_nodes_; ++j) {
799 if (i == j)
continue;
800 min_cost = std::min(min_cost,
Cost(i, j));
801 max_cost = std::max(max_cost,
Cost(i, j));
807 min_cost >= 0 && min_cost > num_nodes_ * max_cost *
808 std::numeric_limits<CostType>::epsilon();
809 robustness_checked_ =
true;
813template <
typename CostType,
typename CostFunction>
816 if (triangle_inequality_checked_)
return triangle_inequality_ok_;
817 triangle_inequality_ok_ =
true;
818 triangle_inequality_checked_ =
true;
819 for (
int k = 0; k < num_nodes_; ++k) {
820 for (
int i = 0; i < num_nodes_; ++i) {
821 for (
int j = 0; j < num_nodes_; ++j) {
822 const CostType detour_cost =
823 Saturated<CostType>::Add(
Cost(i, k),
Cost(k, j));
824 if (detour_cost <
Cost(i, j)) {
825 triangle_inequality_ok_ =
false;
826 return triangle_inequality_ok_;
831 return triangle_inequality_ok_;