Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
hamiltonian_path.h
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14#ifndef OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
15#define OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
16
17// Solves the Shortest Hamiltonian Path Problem using a complete algorithm.
18// The algorithm was first described in
19// M. Held, R.M. Karp, A dynamic programming approach to sequencing problems,
20// J. SIAM 10 (1962) 196-210
21//
22// The Shortest Hamiltonian Path Problem (SHPP) is similar to the Traveling
23// Salesperson Problem (TSP).
24// You have to visit all the cities, starting from a given one and you
25// do not need to return to your starting point. With the TSP, you can start
26// anywhere, but you have to return to your start location.
27//
28// By complete we mean that the algorithm guarantees to compute the optimal
29// solution. The algorithm uses dynamic programming. Its time complexity is
30// O(n^2 * 2^(n-1)), where n is the number of nodes to be visited, and '^'
31// denotes exponentiation. Its space complexity is O(n * 2 ^ (n - 1)).
32//
33// Note that the naive implementation of the SHPP
34// exploring all permutations without memorizing intermediate results would
35// have a complexity of (n - 1)! (factorial of (n - 1) ), which is much higher
36// than n^2 * 2^(n-1). To convince oneself of this, just use Stirling's
37// formula: n! ~ sqrt(2 * pi * n)*( n / exp(1)) ^ n.
38// Because of these complexity figures, the algorithm is not practical for
39// problems with more than 20 nodes.
40//
41// Here is how the algorithm works:
42// Let us denote the nodes to be visited by their indices 0 .. n - 1
43// Let us pick 0 as the starting node.
44// Let cost(i,j) denote the cost (or distance) to go from i to j.
45// f(S, j), where S is a set of nodes and j is a node in S, is defined as the
46// total cost of the shortest path from 0 to j going through all nodes of S.
47//
48// We can prove easily that it satisfy the following relation:
49// f(S, j) = min (i in S \ {j}, f(S \ {j}, i) + cost(i, j))
50// (j is an element of S)
51//
52// Note that this formulation, from the original Held-Karp paper is a bit
53// different, but equivalent to the one used in Caseau and Laburthe, Solving
54// Small TSPs with Constraints, 1997, ICLP
55// f(S, j) = min (i in S, f(S \ {i}, i) + cost(i, j))
56// (j is not an element of S)
57//
58// The advantage of the Held and Karp formulation is that it enables:
59// - to build the dynamic programming lattice layer by layer starting from the
60// subsets with cardinality 1, and increasing the cardinality.
61// - to traverse the dynamic programming lattice using sequential memory
62// accesses, making the algorithm cache-friendly, and faster, despite the large
63// amount of computation needed to get the position when f(S, j) is stored.
64// - TODO(user): implement pruning procedures on top of the Held-Karp algorithm.
65//
66// The set S can be represented by an integer where bit i corresponds to
67// element i in the set. In the following S denotes the integer corresponding
68// to set S.
69//
70// The dynamic programming iteration is implemented in the method Solve.
71// The optimal value of the Hamiltonian path starting at 0 is given by
72// min (i in S, f(2 ^ n - 1, i))
73// The optimal value of the Traveling Salesman tour is given by f(2 ^ n, 0).
74// (There is actually no need to duplicate the first node, as all the paths
75// are computed from node 0.)
76//
77// To implement dynamic programming, we store the preceding results of
78// computing f(S,j) in an array M[Offset(S,j)]. See the comments about
79// LatticeMemoryManager::BaseOffset() to see how this is computed.
80// This is really what brings the performance of the algorithm, because memory
81// is accessed in sequential order, without risking to thrash the cache.
82//
83// Keywords: Traveling Salesman, Hamiltonian Path, Dynamic Programming,
84// Held, Karp.
85
86#include <stddef.h>
87
88#include <algorithm>
89#include <cmath>
90#include <cstdint>
91#include <limits>
92#include <stack>
93#include <type_traits>
94#include <utility>
95#include <vector>
96
97#include "absl/log/check.h"
98#include "absl/types/span.h"
99#include "ortools/util/bitset.h"
102
103namespace operations_research {
104
105// TODO(user): Move the Set-related classbelow to util/bitset.h
106// Iterates over the elements of a set represented as an unsigned integer,
107// starting from the smallest element. (See the class Set<Integer> below.)
108template <typename Set>
110 public:
111 explicit ElementIterator(Set set) : current_set_(set) {}
112 bool operator!=(const ElementIterator& other) const {
113 return current_set_ != other.current_set_;
114 }
115
116 // Returns the smallest element in the current_set_.
117 int operator*() const { return current_set_.SmallestElement(); }
118
119 // Advances the iterator by removing its smallest element.
121 current_set_ = current_set_.RemoveSmallestElement();
122 return *this;
123 }
124
125 private:
126 // The current position of the iterator. Stores the set consisting of the
127 // not-yet iterated elements.
128 Set current_set_;
129};
130
131template <typename Integer>
132class Set {
133 public:
134 // Make this visible to classes using this class.
135 typedef Integer IntegerType;
136
137 // Useful constants.
138 static constexpr Integer kOne = Integer{1};
139 static constexpr Integer kZero = Integer{0};
140 static constexpr int kMaxCardinality = std::numeric_limits<Integer>::digits;
141
142 // Construct a set from an Integer.
143 explicit Set(Integer n) : value_(n) {
144 static_assert(std::is_integral<Integer>::value, "Integral type required");
145 static_assert(std::is_unsigned<Integer>::value, "Unsigned type required");
146 }
147
148 // Returns the integer corresponding to the set.
149 Integer value() const { return value_; }
150
151 static Set FullSet(Integer card) {
152 return card == 0 ? Set(0) : Set(~kZero >> (kMaxCardinality - card));
153 }
154
155 // Returns the singleton set with 'n' as its only element.
156 static Set Singleton(Integer n) { return Set(kOne << n); }
157
158 // Returns a set equal to the calling object, with element n added.
159 // If n is already in the set, no operation occurs.
160 Set AddElement(int n) const { return Set(value_ | (kOne << n)); }
161
162 // Returns a set equal to the calling object, with element n removed.
163 // If n is not in the set, no operation occurs.
164 Set RemoveElement(int n) const { return Set(value_ & ~(kOne << n)); }
165
166 // Returns true if the calling set contains element n.
167 bool Contains(int n) const { return ((kOne << n) & value_) != 0; }
168
169 // Returns true if 'other' is included in the calling set.
170 bool Includes(Set other) const {
171 return (value_ & other.value_) == other.value_;
172 }
173
174 // Returns the number of elements in the set. Uses the 32-bit version for
175 // types that have 32-bits or less. Specialized for uint64_t.
176 int Cardinality() const { return BitCount32(value_); }
177
178 // Returns the index of the smallest element in the set. Uses the 32-bit
179 // version for types that have 32-bits or less. Specialized for uint64_t.
180 int SmallestElement() const { return LeastSignificantBitPosition32(value_); }
181
182 // Returns a set equal to the calling object, with its smallest
183 // element removed.
184 Set RemoveSmallestElement() const { return Set(value_ & (value_ - 1)); }
185
186 // Returns the rank of an element in a set. For the set 11100, ElementRank(4)
187 // would return 2. (Ranks start at kZero).
188 int ElementRank(int n) const {
189 DCHECK(Contains(n)) << "n = " << n << ", value_ = " << value_;
190 return SingletonRank(Singleton(n));
191 }
192
193 // Returns the set consisting of the smallest element of the calling object.
194 Set SmallestSingleton() const { return Set(value_ & -value_); }
195
196 // Returns the rank of the singleton's element in the calling Set.
197 int SingletonRank(Set singleton) const {
198 DCHECK_EQ(singleton.value(), singleton.SmallestSingleton().value());
199 return Set(value_ & (singleton.value_ - 1)).Cardinality();
200 }
201
202 // STL iterator-related member functions.
204 return ElementIterator<Set>(Set(value_));
205 }
207 bool operator!=(const Set& other) const { return value_ != other.value_; }
208
209 private:
210 // The Integer representing the set.
211 Integer value_;
212};
213
214template <>
216 return LeastSignificantBitPosition64(value_);
217}
218
219template <>
220inline int Set<uint64_t>::Cardinality() const {
221 return BitCount64(value_);
222}
223
224// An iterator for sets of increasing corresponding values that have the same
225// cardinality. For example, the sets with cardinality 3 will be listed as
226// ...00111, ...01011, ...01101, ...1110, etc...
227template <typename SetRange>
229 public:
230 // Make the parameter types visible to SetRangeWithCardinality.
231 typedef typename SetRange::SetType SetType;
232 typedef typename SetType::IntegerType IntegerType;
233
234 explicit SetRangeIterator(const SetType set) : current_set_(set) {}
235
236 // STL iterator-related methods.
237 SetType operator*() const { return current_set_; }
238 bool operator!=(const SetRangeIterator& other) const {
239 return current_set_ != other.current_set_;
240 }
241
242 // Computes the next set with the same cardinality using Gosper's hack.
243 // ftp://publications.ai.mit.edu/ai-publications/pdf/AIM-239.pdf ITEM 175
244 // Also translated in C https://www.cl.cam.ac.uk/~am21/hakmemc.html
246 const IntegerType c = current_set_.SmallestSingleton().value();
247 const IntegerType a = current_set_.value();
248 const IntegerType r = c + current_set_.value();
249 // Dividing by c as in HAKMEMC can be avoided by taking into account
250 // that c is the smallest singleton of current_set_, and using a shift.
251 const IntegerType shift = current_set_.SmallestElement();
252 current_set_ = r == 0 ? SetType(0) : SetType(((r ^ a) >> (shift + 2)) | r);
253 return *this;
254 }
255
256 private:
257 // The current set of iterator.
258 SetType current_set_;
259};
260
261template <typename Set>
263 public:
264 typedef Set SetType;
265 // The end_ set is the first set with cardinality card, that does not fit
266 // in max_card bits. Thus, its bit at position max_card is set, and the
267 // rightmost (card - 1) bits are set.
268 SetRangeWithCardinality(int card, int max_card)
269 : begin_(Set::FullSet(card)),
270 end_(Set::FullSet(card - 1).AddElement(max_card)) {
271 DCHECK_LT(0, card);
272 DCHECK_LT(0, max_card);
273 DCHECK_EQ(card, begin_.Cardinality());
274 DCHECK_EQ(card, end_.Cardinality());
275 }
276
277 // STL iterator-related methods.
284
285 private:
286 // Keep the beginning and end of the iterator.
287 SetType begin_;
288 SetType end_;
289};
290
291// The Dynamic Programming (DP) algorithm memorizes the values f(set, node) for
292// node in set, for all the subsets of cardinality <= max_card_.
293// LatticeMemoryManager manages the storage of f(set, node) so that the
294// DP iteration access memory in increasing addresses.
295template <typename Set, typename CostType>
297 public:
298 LatticeMemoryManager() : max_card_(0) {}
299
300 // Reserves memory and fills in the data necessary to access memory.
301 void Init(int max_card);
302
303 // Returns the offset in memory for f(s, node), with node contained in s.
304 uint64_t Offset(Set s, int node) const;
305
306 // Returns the base offset in memory for f(s, node), with node contained in s.
307 // This is useful in the Dynamic Programming iterations.
308 // Note(user): inlining this function gains about 5%.
309 // TODO(user): Investigate how to compute BaseOffset(card - 1, s \ { n })
310 // from BaseOffset(card, n) to speed up the DP iteration.
311 inline uint64_t BaseOffset(int card, Set s) const;
312
313 // Returns the offset delta for a set of cardinality 'card', to which
314 // node 'removed_node' is replaced by 'added_node' at 'rank'
315 uint64_t OffsetDelta(int card, int added_node, int removed_node,
316 int rank) const {
317 return card *
318 (binomial_coefficients_[added_node][rank] - // delta for added_node
319 binomial_coefficients_[removed_node][rank]); // for removed_node.
320 }
321
322 // Memorizes the value = f(s, node) at the correct offset.
323 // This is favored in all other uses than the Dynamic Programming iterations.
324 void SetValue(Set s, int node, CostType value);
325
326 // Memorizes 'value' at 'offset'. This is useful in the Dynamic Programming
327 // iterations where we want to avoid compute the offset of a pair (set, node).
328 void SetValueAtOffset(uint64_t offset, CostType value) {
329 memory_[offset] = value;
330 }
331
332 // Returns the memorized value f(s, node) with node in s.
333 // This is favored in all other uses than the Dynamic Programming iterations.
334 CostType Value(Set s, int node) const;
335
336 // Returns the memorized value at 'offset'.
337 // This is useful in the Dynamic Programming iterations.
338 CostType ValueAtOffset(uint64_t offset) const { return memory_[offset]; }
339
340 private:
341 // Returns true if the values used to manage memory are set correctly.
342 // This is intended to only be used in a DCHECK.
343 bool CheckConsistency() const;
344
345 // The maximum cardinality of the set on which the lattice is going to be
346 // used. This is equal to the number of nodes in the TSP.
347 int max_card_;
348
349 // binomial_coefficients_[n][k] contains (n choose k).
350 std::vector<std::vector<uint64_t>> binomial_coefficients_;
351
352 // base_offset_[card] contains the base offset for all f(set, node) with
353 // card(set) == card.
354 std::vector<int64_t> base_offset_;
355
356 // memory_[Offset(set, node)] contains the costs of the partial path
357 // f(set, node).
358 std::vector<CostType> memory_;
359};
360
361template <typename Set, typename CostType>
363 DCHECK_LT(0, max_card);
364 DCHECK_LE(max_card, Set::kMaxCardinality);
365 if (max_card <= max_card_) return;
366 max_card_ = max_card;
367 binomial_coefficients_.resize(max_card_ + 1);
368
369 // Initialize binomial_coefficients_ using Pascal's triangle recurrence
370 for (int n = 0; n <= max_card_; ++n) {
371 binomial_coefficients_[n].resize(n + 2);
372 binomial_coefficients_[n][0] = 1;
373 for (int k = 1; k <= n; ++k) {
374 binomial_coefficients_[n][k] = binomial_coefficients_[n - 1][k - 1] +
375 binomial_coefficients_[n - 1][k];
376 }
377 // Extend to (n, n + 1) to minimize branchings in LatticeMemoryManager().
378 // This also makes the recurrence above work for k = n.
379 binomial_coefficients_[n][n + 1] = 0;
380 }
381 base_offset_.resize(max_card_ + 1);
382 base_offset_[0] = 0;
383 // There are k * binomial_coefficients_[max_card_][k] f(S,j) values to store
384 // for each group of f(S,j), with card(S) = k. Update base_offset[k]
385 // accordingly.
386 for (int k = 0; k < max_card_; ++k) {
387 base_offset_[k + 1] =
388 base_offset_[k] + k * binomial_coefficients_[max_card_][k];
389 }
390 memory_.resize(0);
391 memory_.shrink_to_fit();
392 memory_.resize(max_card_ * (1 << (max_card_ - 1)));
393 DCHECK(CheckConsistency());
394}
395
396template <typename Set, typename CostType>
397bool LatticeMemoryManager<Set, CostType>::CheckConsistency() const {
398 for (int n = 0; n <= max_card_; ++n) {
399 int64_t sum = 0;
400 for (int k = 0; k <= n; ++k) {
401 sum += binomial_coefficients_[n][k];
402 }
403 DCHECK_EQ(1 << n, sum);
404 }
405 DCHECK_EQ(0, base_offset_[1]);
406 DCHECK_EQ(max_card_ * (1 << (max_card_ - 1)),
407 base_offset_[max_card_] + max_card_);
408 return true;
409}
410
411template <typename Set, typename CostType>
413 Set set) const {
414 DCHECK_LT(0, card);
415 DCHECK_EQ(set.Cardinality(), card);
416 uint64_t local_offset = 0;
417 int node_rank = 0;
418 for (int node : set) {
419 // There are binomial_coefficients_[node][node_rank + 1] sets which have
420 // node at node_rank.
421 local_offset += binomial_coefficients_[node][node_rank + 1];
422 ++node_rank;
423 }
424 DCHECK_EQ(card, node_rank);
425 // Note(user): It is possible to get rid of base_offset_[card] by using a 2-D
426 // array. It would also make it possible to free all the memory but the layer
427 // being constructed and the preceding one, if another lattice of paths is
428 // constructed.
429 // TODO(user): Evaluate the interest of the above.
430 // There are 'card' f(set, j) to store. That is why we need to multiply
431 // local_offset by card before adding it to the corresponding base_offset_.
432 return base_offset_[card] + card * local_offset;
433}
434
435template <typename Set, typename CostType>
437 DCHECK(set.Contains(node));
438 return BaseOffset(set.Cardinality(), set) + set.ElementRank(node);
439}
440
441template <typename Set, typename CostType>
442CostType LatticeMemoryManager<Set, CostType>::Value(Set set, int node) const {
443 DCHECK(set.Contains(node));
444 return ValueAtOffset(Offset(set, node));
445}
446
447template <typename Set, typename CostType>
449 CostType value) {
450 DCHECK(set.Contains(node));
451 SetValueAtOffset(Offset(set, node), value);
452}
453
454// Deprecated type.
455typedef int PathNodeIndex;
456
457template <typename CostType, typename CostFunction>
459 // HamiltonianPathSolver computes a minimum Hamiltonian path starting at node
460 // 0 over a graph defined by a cost matrix. The cost function need not be
461 // symmetric.
462 // When the Hamiltonian path is closed, it's a Hamiltonian cycle,
463 // i.e. the algorithm solves the Traveling Salesman Problem.
464 // Example:
465
466 // std::vector<std::vector<int>> cost_mat;
467 // ... fill in cost matrix
468 // HamiltonianPathSolver<int, std::vector<std::vector<int>>>
469 // mhp(cost_mat); // no computation done
470 // printf("%d\n", mhp.TravelingSalesmanCost()); // computation done and
471 // stored
472 public:
473 // In 2010, 26 was the maximum solvable with 24 Gigs of RAM, and it took
474 // several minutes. With this 2014 version of the code, one may go a little
475 // higher, but considering the complexity of the algorithm (n*2^n), and that
476 // there are very good ways to solve TSP with more than 32 cities,
477 // we limit ourselves to 32 cites.
478 // This is why we define the type NodeSet to be 32-bit wide.
479 // TODO(user): remove this limitation by using pruning techniques.
480 using Integer = uint32_t;
482
483 explicit HamiltonianPathSolver(CostFunction cost);
484 HamiltonianPathSolver(int num_nodes, CostFunction cost);
485
486 // Replaces the cost matrix while avoiding re-allocating memory.
487 void ChangeCostMatrix(CostFunction cost);
488 void ChangeCostMatrix(int num_nodes, CostFunction cost);
489
490 // Returns the cost of the Hamiltonian path from 0 to end_node.
491 CostType HamiltonianCost(int end_node);
492
493 // Returns the shortest Hamiltonian path from 0 to end_node.
494 std::vector<int> HamiltonianPath(int end_node);
495
496 // Returns the end-node that yields the shortest Hamiltonian path of
497 // all shortest Hamiltonian path from 0 to end-node (end-node != 0).
499
500 // Deprecated API. Stores HamiltonianPath(BestHamiltonianPathEndNode()) into
501 // *path.
502 void HamiltonianPath(std::vector<PathNodeIndex>* path);
503
504 // Returns the cost of the TSP tour.
506
507 // Returns the TSP tour in the vector pointed to by the argument.
508 std::vector<int> TravelingSalesmanPath();
509
510 // Deprecated API.
511 void TravelingSalesmanPath(std::vector<PathNodeIndex>* path);
512
513 // Returns true if there won't be precision issues.
514 // This is always true for integers, but not for floating-point types.
515 bool IsRobust();
516
517 // Returns true if the cost matrix verifies the triangle inequality.
519
520 private:
521 // Saturated arithmetic helper class.
522 template <typename T,
523 bool = true /* Dummy parameter to allow specialization */>
524 // Returns the saturated addition of a and b. It is specialized below for
525 // int32_t and int64_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; }
529 };
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); }
534 };
535 // TODO(user): implement this natively in saturated_arithmetic.h
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)));
545 }
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)));
553 }
554 };
555
556 template <typename T>
557 using Saturated = SaturatedArithmetic<T>;
558
559 // Returns the cost value between two nodes.
560 CostType Cost(int i, int j) { return cost_(i, j); }
561
562 // Does all the Dynamic Programming iterations.
563 void Solve();
564
565 // Computes a path by looking at the information in mem_.
566 std::vector<int> ComputePath(CostType cost, NodeSet set, int end);
567
568 // Returns true if the path covers all nodes, and its cost is equal to cost.
569 bool PathIsValid(absl::Span<const int> path, CostType cost);
570
571 // Cost function used to build Hamiltonian paths.
572 MatrixOrFunction<CostType, CostFunction, true> cost_;
573
574 // The number of nodes in the problem.
575 int num_nodes_;
576
577 // The cost of the computed TSP path.
578 CostType tsp_cost_;
579
580 // The cost of the computed Hamiltonian path.
581 std::vector<CostType> hamiltonian_costs_;
582
583 bool robust_;
584 bool triangle_inequality_ok_;
585 bool robustness_checked_;
586 bool triangle_inequality_checked_;
587 bool solved_;
588 std::vector<int> tsp_path_;
589
590 // The vector of smallest Hamiltonian paths starting at 0, indexed by their
591 // end nodes.
592 std::vector<std::vector<int>> hamiltonian_paths_;
593
594 // The end node that gives the smallest Hamiltonian path. The smallest
595 // Hamiltonian path starting at 0 of all
596 // is hamiltonian_paths_[best_hamiltonian_path_end_node_].
597 int best_hamiltonian_path_end_node_;
598
599 LatticeMemoryManager<NodeSet, CostType> mem_;
600};
601
602// Utility function to simplify building a HamiltonianPathSolver from a functor.
603template <typename CostType, typename CostFunction>
605 int num_nodes, CostFunction cost) {
607 std::move(cost));
608}
609
610template <typename CostType, typename CostFunction>
612 CostFunction cost)
613 : HamiltonianPathSolver<CostType, CostFunction>(cost.size(), cost) {}
614
615template <typename CostType, typename CostFunction>
617 int num_nodes, CostFunction cost)
618 : cost_(std::move(cost)),
619 num_nodes_(num_nodes),
620 tsp_cost_(0),
621 hamiltonian_costs_(0),
622 robust_(true),
623 triangle_inequality_ok_(true),
624 robustness_checked_(false),
625 triangle_inequality_checked_(false),
626 solved_(false) {
627 CHECK_GE(NodeSet::kMaxCardinality, num_nodes_);
628 CHECK(cost_.Check());
629}
630
631template <typename CostType, typename CostFunction>
633 CostFunction cost) {
634 ChangeCostMatrix(cost.size(), cost);
635}
636
637template <typename CostType, typename CostFunction>
639 int num_nodes, CostFunction cost) {
640 robustness_checked_ = false;
641 triangle_inequality_checked_ = false;
642 solved_ = false;
643 cost_.Reset(cost);
644 num_nodes_ = num_nodes;
645 CHECK_GE(NodeSet::kMaxCardinality, num_nodes_);
646 CHECK(cost_.Check());
647}
648
649template <typename CostType, typename CostFunction>
650void HamiltonianPathSolver<CostType, CostFunction>::Solve() {
651 if (solved_) return;
652 if (num_nodes_ == 0) {
653 tsp_cost_ = 0;
654 tsp_path_ = {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};
660 return;
661 }
662 mem_.Init(num_nodes_);
663 // Initialize the first layer of the search lattice, taking into account
664 // that base_offset_[1] == 0. (This is what the DCHECK_EQ is for).
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));
668 }
669
670 // Populate the dynamic programming lattice layer by layer, by iterating
671 // on cardinality.
672 for (int card = 2; card <= num_nodes_; ++card) {
673 // Iterate on sets of same cardinality.
674 for (NodeSet set :
675 SetRangeWithCardinality<Set<uint32_t>>(card, num_nodes_)) {
676 // Using BaseOffset and maintaining the node ranks, to reduce the
677 // computational effort for accessing the data.
678 const uint64_t set_offset = mem_.BaseOffset(card, set);
679 // The first subset on which we'll iterate is set.RemoveSmallestElement().
680 // Compute its offset. It will be updated incrementaly. This saves about
681 // 30-35% of computation time.
682 uint64_t subset_offset =
683 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());
684 int prev_dest = set.SmallestElement();
685 int dest_rank = 0;
686 for (int dest : set) {
687 CostType min_cost = std::numeric_limits<CostType>::max();
688 const NodeSet subset = set.RemoveElement(dest);
689 // We compute the offset for subset from the preceding iteration
690 // by taking into account that prev_dest is now in subset, and
691 // that dest is now removed from subset.
692 subset_offset += mem_.OffsetDelta(card - 1, prev_dest, dest, dest_rank);
693 int src_rank = 0;
694 for (int src : subset) {
695 min_cost = std::min(
696 min_cost, Saturated<CostType>::Add(
697 Cost(src, dest),
698 mem_.ValueAtOffset(subset_offset + src_rank)));
699 ++src_rank;
700 }
701 prev_dest = dest;
702 mem_.SetValueAtOffset(set_offset + dest_rank, min_cost);
703 ++dest_rank;
704 }
705 }
706 }
707
708 const NodeSet full_set = NodeSet::FullSet(num_nodes_);
709
710 // Get the cost of the TSP from node 0. It is the path that leaves 0 and goes
711 // through all other nodes, and returns at 0, with minimal cost.
712 tsp_cost_ = mem_.Value(full_set, 0);
713 tsp_path_ = ComputePath(tsp_cost_, full_set, 0);
714
715 hamiltonian_paths_.resize(num_nodes_);
716 hamiltonian_costs_.resize(num_nodes_);
717 // Compute the cost of the Hamiltonian paths starting from node 0, going
718 // through all the other nodes, and ending at end_node. Compute the minimum
719 // one along the way.
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;
728 }
729 DCHECK_LE(tsp_cost_, Saturated<CostType>::Add(cost, Cost(end_node, 0)));
730 // Get the Hamiltonian paths.
731 hamiltonian_paths_[end_node] =
732 ComputePath(hamiltonian_costs_[end_node], hamiltonian_set, end_node);
733 }
734
735 solved_ = true;
736}
737
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;
746 int dest = 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));
753 // Take precision into account when CostType is float or double.
754 // There is no visible penalty in the case CostType is an integer type.
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;
759 path[rank] = src;
760 dest = src;
761 break;
762 }
763 }
764 }
765 DCHECK_EQ(0, subset.value());
766 DCHECK(PathIsValid(path, cost));
767 return path;
768}
769
770template <typename CostType, typename CostFunction>
771bool HamiltonianPathSolver<CostType, CostFunction>::PathIsValid(
772 absl::Span<const int> path, CostType cost) {
773 NodeSet coverage(0);
774 for (int node : path) {
775 coverage = coverage.AddElement(node);
776 }
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) {
780 check_cost =
781 Saturated<CostType>::Add(check_cost, Cost(path[i], path[i + 1]));
782 }
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;
786 return true;
787}
788
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();
795
796 // We compute the min and max for the cost matrix.
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));
802 }
803 }
804 // We determine if the range of the cost matrix is going to
805 // make the algorithm not robust because of precision issues.
806 robust_ =
807 min_cost >= 0 && min_cost > num_nodes_ * max_cost *
808 std::numeric_limits<CostType>::epsilon();
809 robustness_checked_ = true;
810 return robust_;
811}
812
813template <typename CostType, typename CostFunction>
814bool HamiltonianPathSolver<CostType,
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_;
827 }
828 }
829 }
830 }
831 return triangle_inequality_ok_;
832}
833
834template <typename CostType, typename CostFunction>
835int HamiltonianPathSolver<CostType,
837 Solve();
838 return best_hamiltonian_path_end_node_;
839}
840
841template <typename CostType, typename CostFunction>
843 int end_node) {
844 Solve();
845 return hamiltonian_costs_[end_node];
846}
847
848template <typename CostType, typename CostFunction>
850 int end_node) {
851 Solve();
852 return hamiltonian_paths_[end_node];
853}
854
855template <typename CostType, typename CostFunction>
857 std::vector<PathNodeIndex>* path) {
858 *path = HamiltonianPath(best_hamiltonian_path_end_node_);
859}
860
861template <typename CostType, typename CostFunction>
862CostType
867
868template <typename CostType, typename CostFunction>
869std::vector<int>
874
875template <typename CostType, typename CostFunction>
877 std::vector<PathNodeIndex>* path) {
878 *path = TravelingSalesmanPath();
879}
880
881template <typename CostType, typename CostFunction>
883 // PruningHamiltonianSolver computes a minimum Hamiltonian path from node 0
884 // over a graph defined by a cost matrix, with pruning. For each search state,
885 // PruningHamiltonianSolver computes the lower bound for the future overall
886 // TSP cost, and stops further search if it exceeds the current best solution.
887
888 // For the heuristics to determine future lower bound over visited nodeset S
889 // and last visited node k, the cost of minimum spanning tree of (V \ S) ∪ {k}
890 // is calculated and added to the current cost(S). The cost of MST is
891 // guaranteed to be smaller than or equal to the cost of Hamiltonian path,
892 // because Hamiltonian path is a spanning tree itself.
893
894 // TODO(user): Use generic map-based cache instead of lattice-based one.
895 // TODO(user): Use SaturatedArithmetic for better precision.
896
897 public:
898 typedef uint32_t Integer;
900
901 explicit PruningHamiltonianSolver(CostFunction cost);
902 PruningHamiltonianSolver(int num_nodes, CostFunction cost);
903
904 // Returns the cost of the Hamiltonian path from 0 to end_node.
905 CostType HamiltonianCost(int end_node);
906
907 // TODO(user): Add function to return an actual path.
908 // TODO(user): Add functions for Hamiltonian cycle.
909
910 private:
911 // Returns the cost value between two nodes.
912 CostType Cost(int i, int j) { return cost_(i, j); }
913
914 // Solve and get TSP cost.
915 void Solve(int end_node);
916
917 // Compute lower bound for remaining subgraph.
918 CostType ComputeFutureLowerBound(NodeSet current_set, int last_visited);
919
920 // Cost function used to build Hamiltonian paths.
922
923 // The number of nodes in the problem.
924 int num_nodes_;
925
926 // The cost of the computed TSP path.
927 CostType tsp_cost_;
928
929 // If already solved.
930 bool solved_;
931
932 // Memoize for dynamic programming.
934};
935
936template <typename CostType, typename CostFunction>
938 CostFunction cost)
939 : PruningHamiltonianSolver<CostType, CostFunction>(cost.size(), cost) {}
940
941template <typename CostType, typename CostFunction>
943 int num_nodes, CostFunction cost)
944 : cost_(std::move(cost)),
945 num_nodes_(num_nodes),
946 tsp_cost_(0),
947 solved_(false) {}
948
949template <typename CostType, typename CostFunction>
950void PruningHamiltonianSolver<CostType, CostFunction>::Solve(int end_node) {
951 if (solved_ || num_nodes_ == 0) return;
952 // TODO(user): Use an approximate solution as a base target before solving.
953
954 // TODO(user): Instead of pure DFS, find out the order of sets to compute
955 // to utilize cache as possible.
956
957 mem_.Init(num_nodes_);
958 NodeSet start_set = NodeSet::Singleton(0);
959 std::stack<std::pair<NodeSet, int>> state_stack;
960 state_stack.push(std::make_pair(start_set, 0));
961
962 while (!state_stack.empty()) {
963 const std::pair<NodeSet, int> current = state_stack.top();
964 state_stack.pop();
965
966 const NodeSet current_set = current.first;
967 const int last_visited = current.second;
968 const CostType current_cost = mem_.Value(current_set, last_visited);
969
970 // TODO(user): Optimize iterating unvisited nodes.
971 for (int next_to_visit = 0; next_to_visit < num_nodes_; next_to_visit++) {
972 // Let's to as much check possible before adding to stack.
973
974 // Skip if this node is already visited.
975 if (current_set.Contains(next_to_visit)) continue;
976
977 // Skip if the end node is prematurely visited.
978 const int next_cardinality = current_set.Cardinality() + 1;
979 if (next_to_visit == end_node && next_cardinality != num_nodes_) continue;
980
981 const NodeSet next_set = current_set.AddElement(next_to_visit);
982 const CostType next_cost =
983 current_cost + Cost(last_visited, next_to_visit);
984
985 // Compare with the best cost found so far, and skip if that is better.
986 const CostType previous_best = mem_.Value(next_set, next_to_visit);
987 if (previous_best != 0 && next_cost >= previous_best) continue;
988
989 // Compute lower bound of Hamiltonian cost, and skip if this is greater
990 // than the best Hamiltonian cost found so far.
991 const CostType lower_bound =
992 ComputeFutureLowerBound(next_set, next_to_visit);
993 if (tsp_cost_ != 0 && next_cost + lower_bound >= tsp_cost_) continue;
994
995 // If next is the last node to visit, update tsp_cost_ and skip.
996 if (next_cardinality == num_nodes_) {
997 tsp_cost_ = next_cost;
998 continue;
999 }
1000
1001 // Add to the stack, finally.
1002 mem_.SetValue(next_set, next_to_visit, next_cost);
1003 state_stack.push(std::make_pair(next_set, next_to_visit));
1004 }
1005 }
1006
1007 solved_ = true;
1008}
1009
1010template <typename CostType, typename CostFunction>
1012 int end_node) {
1013 Solve(end_node);
1014 return tsp_cost_;
1015}
1016
1017template <typename CostType, typename CostFunction>
1018CostType
1019PruningHamiltonianSolver<CostType, CostFunction>::ComputeFutureLowerBound(
1020 NodeSet current_set, int last_visited) {
1021 // TODO(user): Compute MST.
1022 return 0; // For now, return 0 as future lower bound.
1023}
1024} // namespace operations_research
1025
1026#endif // OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
bool operator!=(const ElementIterator &other) const
int operator*() const
Returns the smallest element in the current_set_.
const ElementIterator & operator++()
Advances the iterator by removing its smallest element.
CostType TravelingSalesmanCost()
Returns the cost of the TSP tour.
void ChangeCostMatrix(CostFunction cost)
Replaces the cost matrix while avoiding re-allocating memory.
CostType HamiltonianCost(int end_node)
Returns the cost of the Hamiltonian path from 0 to end_node.
std::vector< int > HamiltonianPath(int end_node)
Returns the shortest Hamiltonian path from 0 to end_node.
std::vector< int > TravelingSalesmanPath()
Returns the TSP tour in the vector pointed to by the argument.
bool VerifiesTriangleInequality()
Returns true if the cost matrix verifies the triangle inequality.
CostType Value(Set s, int node) const
CostType ValueAtOffset(uint64_t offset) const
void SetValue(Set s, int node, CostType value)
uint64_t OffsetDelta(int card, int added_node, int removed_node, int rank) const
uint64_t Offset(Set s, int node) const
Returns the offset in memory for f(s, node), with node contained in s.
void Init(int max_card)
Reserves memory and fills in the data necessary to access memory.
void SetValueAtOffset(uint64_t offset, CostType value)
uint64_t BaseOffset(int card, Set s) const
CostType HamiltonianCost(int end_node)
Returns the cost of the Hamiltonian path from 0 to end_node.
SetType operator*() const
STL iterator-related methods.
SetRange::SetType SetType
Make the parameter types visible to SetRangeWithCardinality.
const SetRangeIterator & operator++()
bool operator!=(const SetRangeIterator &other) const
SetRangeIterator< SetRangeWithCardinality > end() const
SetRangeIterator< SetRangeWithCardinality > begin() const
STL iterator-related methods.
bool operator!=(const Set &other) const
static Set FullSet(Integer card)
bool Includes(Set other) const
Returns true if 'other' is included in the calling set.
Set(Integer n)
Construct a set from an Integer.
bool Contains(int n) const
Returns true if the calling set contains element n.
Integer IntegerType
Make this visible to classes using this class.
static Set Singleton(Integer n)
Returns the singleton set with 'n' as its only element.
ElementIterator< Set > end() const
ElementIterator< Set > begin() const
STL iterator-related member functions.
int SingletonRank(Set singleton) const
Returns the rank of the singleton's element in the calling Set.
Set RemoveElement(int n) const
static constexpr int kMaxCardinality
Set SmallestSingleton() const
Returns the set consisting of the smallest element of the calling object.
static constexpr Integer kZero
Integer value() const
Returns the integer corresponding to the set.
static constexpr Integer kOne
Useful constants.
In SWIG mode, we don't want anything besides these top-level includes.
int64_t CapAdd(int64_t x, int64_t y)
uint32_t BitCount32(uint32_t n)
Definition bitset.h:56
int64_t CapSub(int64_t x, int64_t y)
ClosedInterval::Iterator end(ClosedInterval interval)
uint64_t BitCount64(uint64_t n)
Returns the number of bits set in n.
Definition bitset.h:45
double Cost
Basic non-strict type for cost. The speed penalty for using double is ~2%.
Definition base_types.h:32
int LeastSignificantBitPosition64(uint64_t n)
Definition bitset.h:130
HamiltonianPathSolver< CostType, CostFunction > MakeHamiltonianPathSolver(int num_nodes, CostFunction cost)
Utility function to simplify building a HamiltonianPathSolver from a functor.
int LeastSignificantBitPosition32(uint32_t n)
Definition bitset.h:185
STL namespace.