C++ Reference

C++ Reference: Graph

hamiltonian_path.h
Go to the documentation of this file.
1// Copyright 2010-2018 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 d(i,j) denote the distance (or cost) 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 follows:
46// f(S, j) = min (i in S \ {j}, f(S \ {j}, i) + cost(i, j))
47// (j is an element of S)
48// Note that this formulation, from the original Held-Karp paper is a bit
49// different, but equivalent to the one used in Caseau and Laburthe, Solving
50// Small TSPs with Constraints, 1997, ICLP
51// f(S, j) = min (i in S, f(S \ {i}, i) + cost(i, j))
52// (j is not an element of S)
53//
54// The advantage of the Held and Karp formulation is that it enables:
55// - to build the dynamic programming lattice layer by layer starting from the
56// subsets with cardinality 1, and increasing the cardinality.
57// - to traverse the dynamic programming lattice using sequential memory
58// accesses, making the algorithm cache-friendly, and faster, despite the large
59// amount of computation needed to get the position when f(S, j) is stored.
60// - TODO(user): implement pruning procedures on top of the Held-Karp algorithm.
61//
62// The set S can be represented by an integer where bit i corresponds to
63// element i in the set. In the following S denotes the integer corresponding
64// to set S.
65//
66// The dynamic programming iteration is implemented in the method Solve.
67// The optimal value of the Hamiltonian path starting at 0 is given by
68// min (i in S, f(2 ^ n - 1, i))
69// The optimal value of the Traveling Salesman tour is given by f(2 ^ n, 0).
70// (There is actually no need to duplicate the first node, as all the paths
71// are computed from node 0.)
72//
73// To implement dynamic programming, we store the preceding results of
74// computing f(S,j) in an array M[Offset(S,j)]. See the comments about
75// LatticeMemoryManager::BaseOffset() to see how this is computed.
76//
77// Keywords: Traveling Salesman, Hamiltonian Path, Dynamic Programming,
78// Held, Karp.
79
80#include <math.h>
81#include <stddef.h>
82
83#include <algorithm>
84#include <cmath>
85#include <limits>
86#include <memory>
87#include <stack>
88#include <type_traits>
89#include <utility>
90#include <vector>
91
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"
97
98namespace operations_research {
99
100// TODO(user): Move the Set-related classbelow to util/bitset.h
101// Iterates over the elements of a set represented as an unsigned integer,
102// starting from the smallest element. (See the class Set<Integer> below.)
103template <typename Set>
105 public:
106 explicit ElementIterator(Set set) : current_set_(set) {}
107 bool operator!=(const ElementIterator& other) const {
108 return current_set_ != other.current_set_;
109 }
110
111 // Returns the smallest element in the current_set_.
112 int operator*() const { return current_set_.SmallestElement(); }
113
114 // Advances the iterator by removing its smallest element.
116 current_set_ = current_set_.RemoveSmallestElement();
117 return *this;
118 }
119
120 private:
121 // The current position of the iterator. Stores the set consisting of the
122 // not-yet iterated elements.
123 Set current_set_;
124};
125
126template <typename Integer>
127class Set {
128 public:
129 // Make this visible to classes using this class.
130 typedef Integer IntegerType;
131
132 // Useful constants.
133 static constexpr Integer One = static_cast<Integer>(1);
134 static constexpr Integer Zero = static_cast<Integer>(0);
135 static const int MaxCardinality = 8 * sizeof(Integer); // NOLINT
136
137 // Construct a set from an Integer.
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");
141 }
142
143 // Returns the integer corresponding to the set.
144 Integer value() const { return value_; }
145
146 static Set FullSet(Integer card) {
147 return card == 0 ? Set(0) : Set(~Zero >> (MaxCardinality - card));
148 }
149
150 // Returns the singleton set with 'n' as its only element.
151 static Set Singleton(Integer n) { return Set(One << n); }
152
153 // Returns a set equal to the calling object, with element n added.
154 // If n is already in the set, no operation occurs.
155 Set AddElement(int n) const { return Set(value_ | (One << n)); }
156
157 // Returns a set equal to the calling object, with element n removed.
158 // If n is not in the set, no operation occurs.
159 Set RemoveElement(int n) const { return Set(value_ & ~(One << n)); }
160
161 // Returns true if the calling set contains element n.
162 bool Contains(int n) const { return ((One << n) & value_) != 0; }
163
164 // Returns true if 'other' is included in the calling set.
165 bool Includes(Set other) const {
166 return (value_ & other.value_) == other.value_;
167 }
168
169 // Returns the number of elements in the set. Uses the 32-bit version for
170 // types that have 32-bits or less. Specialized for uint64.
171 int Cardinality() const { return BitCount32(value_); }
172
173 // Returns the index of the smallest element in the set. Uses the 32-bit
174 // version for types that have 32-bits or less. Specialized for uint64.
175 int SmallestElement() const { return LeastSignificantBitPosition32(value_); }
176
177 // Returns a set equal to the calling object, with its smallest
178 // element removed.
179 Set RemoveSmallestElement() const { return Set(value_ & (value_ - 1)); }
180
181 // Returns the rank of an element in a set. For the set 11100, ElementRank(4)
182 // would return 2. (Ranks start at zero).
183 int ElementRank(int n) const {
184 DCHECK(Contains(n)) << "n = " << n << ", value_ = " << value_;
185 return SingletonRank(Singleton(n));
186 }
187
188 // Returns the set consisting of the smallest element of the calling object.
189 Set SmallestSingleton() const { return Set(value_ & -value_); }
190
191 // Returns the rank of the singleton's element in the calling Set.
192 int SingletonRank(Set singleton) const {
193 DCHECK_EQ(singleton.value(), singleton.SmallestSingleton().value());
194 return Set(value_ & (singleton.value_ - 1)).Cardinality();
195 }
196
197 // STL iterator-related member functions.
199 return ElementIterator<Set>(Set(value_));
200 }
202 bool operator!=(const Set& other) const { return value_ != other.value_; }
203
204 private:
205 // The Integer representing the set.
206 Integer value_;
207};
208
209template <>
211 return LeastSignificantBitPosition64(value_);
212}
213
214template <>
215inline int Set<uint64>::Cardinality() const {
216 return BitCount64(value_);
217}
218
219// An iterator for sets of increasing corresponding values that have the same
220// cardinality. For example, the sets with cardinality 3 will be listed as
221// ...00111, ...01011, ...01101, ...1110, etc...
222template <typename SetRange>
224 public:
225 // Make the parameter types visible to SetRangeWithCardinality.
226 typedef typename SetRange::SetType SetType;
227 typedef typename SetType::IntegerType IntegerType;
228
229 explicit SetRangeIterator(const SetType set) : current_set_(set) {}
230
231 // STL iterator-related methods.
232 SetType operator*() const { return current_set_; }
233 bool operator!=(const SetRangeIterator& other) const {
234 return current_set_ != other.current_set_;
235 }
236
237 // Computes the next set with the same cardinality using Gosper's hack.
238 // ftp://publications.ai.mit.edu/ai-publications/pdf/AIM-239.pdf ITEM 175
239 // Also translated in C https://www.cl.cam.ac.uk/~am21/hakmemc.html
241 const IntegerType c = current_set_.SmallestSingleton().value();
242 const IntegerType a = current_set_.value();
243 const IntegerType r = c + current_set_.value();
244 // Dividing by c as in HAKMEMC can be avoided by taking into account
245 // that c is the smallest singleton of current_set_, and using a shift.
246 const IntegerType shift = current_set_.SmallestElement();
247 current_set_ = r == 0 ? SetType(0) : SetType(((r ^ a) >> (shift + 2)) | r);
248 return *this;
249 }
250
251 private:
252 // The current set of iterator.
253 SetType current_set_;
254};
255
256template <typename Set>
258 public:
259 typedef Set SetType;
260 // The end_ set is the first set with cardinality card, that does not fit
261 // in max_card bits. Thus, its bit at position max_card is set, and the
262 // rightmost (card - 1) bits are set.
263 SetRangeWithCardinality(int card, int max_card)
264 : begin_(Set::FullSet(card)),
265 end_(Set::FullSet(card - 1).AddElement(max_card)) {
266 DCHECK_LT(0, card);
267 DCHECK_LT(0, max_card);
268 DCHECK_EQ(card, begin_.Cardinality());
269 DCHECK_EQ(card, end_.Cardinality());
270 }
271
272 // STL iterator-related methods.
275 }
278 }
279
280 private:
281 // Keep the beginning and end of the iterator.
282 SetType begin_;
283 SetType end_;
284};
285
286// The Dynamic Programming (DP) algorithm memorizes the values f(set, node) for
287// node in set, for all the subsets of cardinality <= max_card_.
288// LatticeMemoryManager manages the storage of f(set, node) so that the
289// DP iteration access memory in increasing addresses.
290template <typename Set, typename CostType>
292 public:
293 LatticeMemoryManager() : max_card_(0) {}
294
295 // Reserves memory and fills in the data necessary to access memory.
296 void Init(int max_card);
297
298 // Returns the offset in memory for f(s, node), with node contained in s.
299 uint64 Offset(Set s, int node) const;
300
301 // Returns the base offset in memory for f(s, node), with node contained in s.
302 // This is useful in the Dynamic Programming iterations.
303 // Note(user): inlining this function gains about 5%.
304 // TODO(user): Investigate how to compute BaseOffset(card - 1, s \ { n })
305 // from BaseOffset(card, n) to speed up the DP iteration.
306 inline uint64 BaseOffset(int card, Set s) const;
307
308 // Returns the offset delta for a set of cardinality 'card', to which
309 // node 'removed_node' is replaced by 'added_node' at 'rank'
310 uint64 OffsetDelta(int card, int added_node, int removed_node,
311 int rank) const {
312 return card *
313 (binomial_coefficients_[added_node][rank] - // delta for added_node
314 binomial_coefficients_[removed_node][rank]); // for removed_node.
315 }
316
317 // Memorizes the value = f(s, node) at the correct offset.
318 // This is favored in all other uses than the Dynamic Programming iterations.
319 void SetValue(Set s, int node, CostType value);
320
321 // Memorizes 'value' at 'offset'. This is useful in the Dynamic Programming
322 // iterations where we want to avoid compute the offset of a pair (set, node).
323 void SetValueAtOffset(uint64 offset, CostType value) {
324 memory_[offset] = value;
325 }
326
327 // Returns the memorized value f(s, node) with node in s.
328 // This is favored in all other uses than the Dynamic Programming iterations.
329 CostType Value(Set s, int node) const;
330
331 // Returns the memorized value at 'offset'.
332 // This is useful in the Dynamic Programming iterations.
333 CostType ValueAtOffset(uint64 offset) const { return memory_[offset]; }
334
335 private:
336 // Returns true if the values used to manage memory are set correctly.
337 // This is intended to only be used in a DCHECK.
338 bool CheckConsistency() const;
339
340 // The maximum cardinality of the set on which the lattice is going to be
341 // used. This is equal to the number of nodes in the TSP.
342 int max_card_;
343
344 // binomial_coefficients_[n][k] contains (n choose k).
345 std::vector<std::vector<uint64>> binomial_coefficients_;
346
347 // base_offset_[card] contains the base offset for all f(set, node) with
348 // card(set) == card.
349 std::vector<int64> base_offset_;
350
351 // memory_[Offset(set, node)] contains the costs of the partial path
352 // f(set, node).
353 std::vector<CostType> memory_;
354};
355
356template <typename Set, typename CostType>
358 DCHECK_LT(0, max_card);
359 DCHECK_GE(Set::MaxCardinality, max_card);
360 if (max_card <= max_card_) return;
361 max_card_ = max_card;
362 binomial_coefficients_.resize(max_card_ + 1);
363
364 // Initialize binomial_coefficients_ using Pascal's triangle recursion.
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];
371 }
372 // Extend to (n, n + 1) to minimize branchings in LatticeMemoryManager().
373 // This also makes the recurrence above work for k = n.
374 binomial_coefficients_[n][n + 1] = 0;
375 }
376 base_offset_.resize(max_card_ + 1);
377 base_offset_[0] = 0;
378 // There are k * binomial_coefficients_[max_card_][k] f(S,j) values to store
379 // for each group of f(S,j), with card(S) = k. Update base_offset[k]
380 // accordingly.
381 for (int k = 0; k < max_card_; ++k) {
382 base_offset_[k + 1] =
383 base_offset_[k] + k * binomial_coefficients_[max_card_][k];
384 }
385 memory_.resize(0);
386 memory_.shrink_to_fit();
387 memory_.resize(max_card_ * (1 << (max_card_ - 1)));
388 DCHECK(CheckConsistency());
389}
390
391template <typename Set, typename CostType>
393 for (int n = 0; n <= max_card_; ++n) {
394 int64 sum = 0;
395 for (int k = 0; k <= n; ++k) {
396 sum += binomial_coefficients_[n][k];
397 }
398 DCHECK_EQ(1 << n, sum);
399 }
400 DCHECK_EQ(0, base_offset_[1]);
401 DCHECK_EQ(max_card_ * (1 << (max_card_ - 1)),
402 base_offset_[max_card_] + max_card_);
403 return true;
404}
405
406template <typename Set, typename CostType>
408 Set set) const {
409 DCHECK_LT(0, card);
410 DCHECK_EQ(set.Cardinality(), card);
411 uint64 local_offset = 0;
412 int node_rank = 0;
413 for (int node : set) {
414 // There are binomial_coefficients_[node][node_rank + 1] sets which have
415 // node at node_rank.
416 local_offset += binomial_coefficients_[node][node_rank + 1];
417 ++node_rank;
418 }
419 DCHECK_EQ(card, node_rank);
420 // Note(user): It is possible to get rid of base_offset_[card] by using a 2-D
421 // array. It would also make it possible to free all the memory but the layer
422 // being constructed and the preceding one, if another lattice of paths is
423 // constructed.
424 // TODO(user): Evaluate the interest of the above.
425 // There are 'card' f(set, j) to store. That is why we need to multiply
426 // local_offset by card before adding it to the corresponding base_offset_.
427 return base_offset_[card] + card * local_offset;
428}
429
430template <typename Set, typename CostType>
432 DCHECK(set.Contains(node));
433 return BaseOffset(set.Cardinality(), set) + set.ElementRank(node);
434}
435
436template <typename Set, typename CostType>
437CostType LatticeMemoryManager<Set, CostType>::Value(Set set, int node) const {
438 DCHECK(set.Contains(node));
439 return ValueAtOffset(Offset(set, node));
440}
441
442template <typename Set, typename CostType>
444 CostType value) {
445 DCHECK(set.Contains(node));
446 SetValueAtOffset(Offset(set, node), value);
447}
448
449// Deprecated type.
450typedef int PathNodeIndex;
451
452template <typename CostType, typename CostFunction>
454 // HamiltonianPathSolver computes a minimum Hamiltonian path starting at node
455 // 0 over a graph defined by a cost matrix. The cost function need not be
456 // symmetric.
457 // When the Hamiltonian path is closed, it's a Hamiltonian cycle,
458 // i.e. the algorithm solves the Traveling Salesman Problem.
459 // Example:
460
461 // std::vector<std::vector<int>> cost_mat;
462 // ... fill in cost matrix
463 // HamiltonianPathSolver<int, std::vector<std::vector<int>>>
464 // mhp(cost_mat); // no computation done
465 // printf("%d\n", mhp.TravelingSalesmanCost()); // computation done and
466 // stored
467 public:
468 // In 2010, 26 was the maximum solvable with 24 Gigs of RAM, and it took
469 // several minutes. With this 2014 version of the code, one may go a little
470 // higher, but considering the complexity of the algorithm (n*2^n), and that
471 // there are very good ways to solve TSP with more than 32 cities,
472 // we limit ourselves to 32 cites.
473 // This is why we define the type NodeSet to be 32-bit wide.
474 // TODO(user): remove this limitation by using pruning techniques.
475 typedef uint32 Integer;
477
478 explicit HamiltonianPathSolver(CostFunction cost);
479 HamiltonianPathSolver(int num_nodes, CostFunction cost);
480
481 // Replaces the cost matrix while avoiding re-allocating memory.
482 void ChangeCostMatrix(CostFunction cost);
483 void ChangeCostMatrix(int num_nodes, CostFunction cost);
484
485 // Returns the cost of the Hamiltonian path from 0 to end_node.
486 CostType HamiltonianCost(int end_node);
487
488 // Returns the shortest Hamiltonian path from 0 to end_node.
489 std::vector<int> HamiltonianPath(int end_node);
490
491 // Returns the end-node that yields the shortest Hamiltonian path of
492 // all shortest Hamiltonian path from 0 to end-node (end-node != 0).
494
495 // Deprecated API. Stores HamiltonianPath(BestHamiltonianPathEndNode()) into
496 // *path.
497 void HamiltonianPath(std::vector<PathNodeIndex>* path);
498
499 // Returns the cost of the TSP tour.
500 CostType TravelingSalesmanCost();
501
502 // Returns the TSP tour in the vector pointed to by the argument.
503 std::vector<int> TravelingSalesmanPath();
504
505 // Deprecated API.
506 void TravelingSalesmanPath(std::vector<PathNodeIndex>* path);
507
508 // Returns true if there won't be precision issues.
509 // This is always true for integers, but not for floating-point types.
510 bool IsRobust();
511
512 // Returns true if the cost matrix verifies the triangle inequality.
514
515 private:
516 // Saturated arithmetic helper class.
517 template <typename T,
518 bool = true /* Dummy parameter to allow specialization */>
519 // Returns the saturated addition of a and b. It is specialized below for
520 // int32 and int64.
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; }
524 };
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); }
529 };
530 // TODO(user): implement this natively in saturated_arithmetic.h
531 template <bool Dummy>
532 struct SaturatedArithmetic<int32, Dummy> {
533 static int32 Add(int32 a, int32 b) {
534 const int64 a64 = a;
535 const int64 b64 = 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)));
540 }
541 static int32 Sub(int32 a, int32 b) {
542 const int64 a64 = a;
543 const int64 b64 = 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)));
548 }
549 };
550
551 template <typename T>
552 using Saturated = SaturatedArithmetic<T>;
553
554 // Returns the cost value between two nodes.
555 CostType Cost(int i, int j) { return cost_(i, j); }
556
557 // Does all the Dynamic Progamming iterations.
558 void Solve();
559
560 // Computes a path by looking at the information in mem_.
561 std::vector<int> ComputePath(CostType cost, NodeSet set, int end);
562
563 // Returns true if the path covers all nodes, and its cost is equal to cost.
564 bool PathIsValid(const std::vector<int>& path, CostType cost);
565
566 // Cost function used to build Hamiltonian paths.
567 MatrixOrFunction<CostType, CostFunction, true> cost_;
568
569 // The number of nodes in the problem.
570 int num_nodes_;
571
572 // The cost of the computed TSP path.
573 CostType tsp_cost_;
574
575 // The cost of the computed Hamiltonian path.
576 std::vector<CostType> hamiltonian_costs_;
577
578 bool robust_;
579 bool triangle_inequality_ok_;
580 bool robustness_checked_;
581 bool triangle_inequality_checked_;
582 bool solved_;
583 std::vector<int> tsp_path_;
584
585 // The vector of smallest Hamiltonian paths starting at 0, indexed by their
586 // end nodes.
587 std::vector<std::vector<int>> hamiltonian_paths_;
588
589 // The end node that gives the smallest Hamiltonian path. The smallest
590 // Hamiltonian path starting at 0 of all
591 // is hamiltonian_paths_[best_hamiltonian_path_end_node_].
592 int best_hamiltonian_path_end_node_;
593
594 LatticeMemoryManager<NodeSet, CostType> mem_;
595};
596
597// Utility function to simplify building a HamiltonianPathSolver from a functor.
598template <typename CostType, typename CostFunction>
600 int num_nodes, CostFunction cost) {
602 std::move(cost));
603}
604
605template <typename CostType, typename CostFunction>
607 CostFunction cost)
608 : HamiltonianPathSolver<CostType, CostFunction>(cost.size(), cost) {}
609
610template <typename CostType, typename CostFunction>
612 int num_nodes, CostFunction cost)
613 : cost_(std::move(cost)),
614 num_nodes_(num_nodes),
615 tsp_cost_(0),
616 hamiltonian_costs_(0),
617 robust_(true),
618 triangle_inequality_ok_(true),
619 robustness_checked_(false),
620 triangle_inequality_checked_(false),
621 solved_(false) {
622 CHECK_GE(NodeSet::MaxCardinality, num_nodes_);
623 CHECK(cost_.Check());
624}
625
626template <typename CostType, typename CostFunction>
628 CostFunction cost) {
629 ChangeCostMatrix(cost.size(), cost);
630}
631
632template <typename CostType, typename CostFunction>
634 int num_nodes, CostFunction cost) {
635 robustness_checked_ = false;
636 triangle_inequality_checked_ = false;
637 solved_ = false;
638 cost_.Reset(cost);
639 num_nodes_ = num_nodes;
640 CHECK_GE(NodeSet::MaxCardinality, num_nodes_);
641 CHECK(cost_.Check());
642}
643
644template <typename CostType, typename CostFunction>
646 if (solved_) return;
647 if (num_nodes_ == 0) {
648 tsp_cost_ = 0;
649 tsp_path_ = {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};
655 return;
656 }
657 mem_.Init(num_nodes_);
658 // Initialize the first layer of the search lattice, taking into account
659 // that base_offset_[1] == 0. (This is what the DCHECK_EQ is for).
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));
663 }
664
665 // Populate the dynamic programming lattice layer by layer, by iterating
666 // on cardinality.
667 for (int card = 2; card <= num_nodes_; ++card) {
668 // Iterate on sets of same cardinality.
669 for (NodeSet set : SetRangeWithCardinality<Set<uint32>>(card, num_nodes_)) {
670 // Using BaseOffset and maintaining the node ranks, to reduce the
671 // computational effort for accessing the data.
672 const uint64 set_offset = mem_.BaseOffset(card, set);
673 // The first subset on which we'll iterate is set.RemoveSmallestElement().
674 // Compute its offset. It will be updated incrementaly. This saves about
675 // 30-35% of computation time.
676 uint64 subset_offset =
677 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());
678 int prev_dest = set.SmallestElement();
679 int dest_rank = 0;
680 for (int dest : set) {
681 CostType min_cost = std::numeric_limits<CostType>::max();
682 const NodeSet subset = set.RemoveElement(dest);
683 // We compute the offset for subset from the preceding iteration
684 // by taking into account that prev_dest is now in subset, and
685 // that dest is now removed from subset.
686 subset_offset += mem_.OffsetDelta(card - 1, prev_dest, dest, dest_rank);
687 int src_rank = 0;
688 for (int src : subset) {
689 min_cost = std::min(
690 min_cost, Saturated<CostType>::Add(
691 Cost(src, dest),
692 mem_.ValueAtOffset(subset_offset + src_rank)));
693 ++src_rank;
694 }
695 prev_dest = dest;
696 mem_.SetValueAtOffset(set_offset + dest_rank, min_cost);
697 ++dest_rank;
698 }
699 }
700 }
701
702 const NodeSet full_set = NodeSet::FullSet(num_nodes_);
703
704 // Get the cost of the tsp from node 0. It is the path that leaves 0 and goes
705 // through all other nodes, and returns at 0, with minimal cost.
706 tsp_cost_ = mem_.Value(full_set, 0);
707 tsp_path_ = ComputePath(tsp_cost_, full_set, 0);
708
709 hamiltonian_paths_.resize(num_nodes_);
710 hamiltonian_costs_.resize(num_nodes_);
711 // Compute the cost of the Hamiltonian paths starting from node 0, going
712 // through all the other nodes, and ending at end_node. Compute the minimum
713 // one along the way.
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;
722 }
723 DCHECK_LE(tsp_cost_, Saturated<CostType>::Add(cost, Cost(end_node, 0)));
724 // Get the Hamiltonian paths.
725 hamiltonian_paths_[end_node] =
726 ComputePath(hamiltonian_costs_[end_node], hamiltonian_set, end_node);
727 }
728
729 solved_ = true;
730}
731
732template <typename CostType, typename CostFunction>
733std::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;
740 int dest = 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));
747 // Take precision into account when CosttType is float or double.
748 // There is no visible penalty in the case CostType is an integer type.
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;
753 path[rank] = src;
754 dest = src;
755 break;
756 }
757 }
758 }
759 DCHECK_EQ(0, subset.value());
760 DCHECK(PathIsValid(path, cost));
761 return path;
762}
763
764template <typename CostType, typename CostFunction>
765bool HamiltonianPathSolver<CostType, CostFunction>::PathIsValid(
766 const std::vector<int>& path, CostType cost) {
767 NodeSet coverage(0);
768 for (int node : path) {
769 coverage = coverage.AddElement(node);
770 }
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) {
774 check_cost =
775 Saturated<CostType>::Add(check_cost, Cost(path[i], path[i + 1]));
776 }
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;
780 return true;
781}
782
783template <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();
789
790 // We compute the min and max for the cost matrix.
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));
796 }
797 }
798 // We determine if the range of the cost matrix is going to
799 // make the algorithm not robust because of precision issues.
800 robust_ =
801 min_cost >= 0 && min_cost > num_nodes_ * max_cost *
802 std::numeric_limits<CostType>::epsilon();
803 robustness_checked_ = true;
804 return robust_;
805}
806
807template <typename CostType, typename CostFunction>
808bool HamiltonianPathSolver<CostType,
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_;
821 }
822 }
823 }
824 }
825 return triangle_inequality_ok_;
826}
827
828template <typename CostType, typename CostFunction>
829int HamiltonianPathSolver<CostType,
830 CostFunction>::BestHamiltonianPathEndNode() {
831 Solve();
832 return best_hamiltonian_path_end_node_;
833}
834
835template <typename CostType, typename CostFunction>
837 int end_node) {
838 Solve();
839 return hamiltonian_costs_[end_node];
840}
841
842template <typename CostType, typename CostFunction>
844 int end_node) {
845 Solve();
846 return hamiltonian_paths_[end_node];
847}
848
849template <typename CostType, typename CostFunction>
851 std::vector<PathNodeIndex>* path) {
852 *path = HamiltonianPath(best_hamiltonian_path_end_node_);
853}
854
855template <typename CostType, typename CostFunction>
856CostType
858 Solve();
859 return tsp_cost_;
860}
861
862template <typename CostType, typename CostFunction>
863std::vector<int>
865 Solve();
866 return tsp_path_;
867}
868
869template <typename CostType, typename CostFunction>
871 std::vector<PathNodeIndex>* path) {
872 *path = TravelingSalesmanPath();
873}
874
875template <typename CostType, typename CostFunction>
877 // PruningHamiltonianSolver computes a minimum Hamiltonian path from node 0
878 // over a graph defined by a cost matrix, with pruning. For each search state,
879 // PruningHamiltonianSolver computes the lower bound for the future overall
880 // TSP cost, and stops further search if it exceeds the current best solution.
881
882 // For the heuristics to determine future lower bound over visited nodeset S
883 // and last visited node k, the cost of minimum spanning tree of (V \ S) ∪ {k}
884 // is calculated and added to the current cost(S). The cost of MST is
885 // guaranteed to be smaller than or equal to the cost of Hamiltonian path,
886 // because Hamiltonian path is a spanning tree itself.
887
888 // TODO(user): Use generic map-based cache instead of lattice-based one.
889 // TODO(user): Use SaturatedArithmetic for better precision.
890
891 public:
892 typedef uint32 Integer;
894
895 explicit PruningHamiltonianSolver(CostFunction cost);
896 PruningHamiltonianSolver(int num_nodes, CostFunction cost);
897
898 // Returns the cost of the Hamiltonian path from 0 to end_node.
899 CostType HamiltonianCost(int end_node);
900
901 // TODO(user): Add function to return an actual path.
902 // TODO(user): Add functions for Hamiltonian cycle.
903
904 private:
905 // Returns the cost value between two nodes.
906 CostType Cost(int i, int j) { return cost_(i, j); }
907
908 // Solve and get TSP cost.
909 void Solve(int end_node);
910
911 // Compute lower bound for remaining subgraph.
912 CostType ComputeFutureLowerBound(NodeSet current_set, int last_visited);
913
914 // Cost function used to build Hamiltonian paths.
915 MatrixOrFunction<CostType, CostFunction, true> cost_;
916
917 // The number of nodes in the problem.
918 int num_nodes_;
919
920 // The cost of the computed TSP path.
921 CostType tsp_cost_;
922
923 // If already solved.
924 bool solved_;
925
926 // Memoize for dynamic programming.
928};
929
930template <typename CostType, typename CostFunction>
932 CostFunction cost)
933 : PruningHamiltonianSolver<CostType, CostFunction>(cost.size(), cost) {}
934
935template <typename CostType, typename CostFunction>
937 int num_nodes, CostFunction cost)
938 : cost_(std::move(cost)),
939 num_nodes_(num_nodes),
940 tsp_cost_(0),
941 solved_(false) {}
942
943template <typename CostType, typename CostFunction>
945 if (solved_ || num_nodes_ == 0) return;
946 // TODO(user): Use an approximate solution as a base target before solving.
947
948 // TODO(user): Instead of pure DFS, find out the order of sets to compute
949 // to utilize cache as possible.
950
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));
955
956 while (!state_stack.empty()) {
957 const std::pair<NodeSet, int> current = state_stack.top();
958 state_stack.pop();
959
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);
963
964 // TODO(user): Optimize iterating unvisited nodes.
965 for (int next_to_visit = 0; next_to_visit < num_nodes_; next_to_visit++) {
966 // Let's to as much check possible before adding to stack.
967
968 // Skip if this node is already visited.
969 if (current_set.Contains(next_to_visit)) continue;
970
971 // Skip if the end node is prematurely visited.
972 const int next_cardinality = current_set.Cardinality() + 1;
973 if (next_to_visit == end_node && next_cardinality != num_nodes_) continue;
974
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);
978
979 // Compare with the best cost found so far, and skip if that is better.
980 const CostType previous_best = mem_.Value(next_set, next_to_visit);
981 if (previous_best != 0 && next_cost >= previous_best) continue;
982
983 // Compute lower bound of Hamiltonian cost, and skip if this is greater
984 // than the best Hamiltonian cost found so far.
985 const CostType lower_bound =
986 ComputeFutureLowerBound(next_set, next_to_visit);
987 if (tsp_cost_ != 0 && next_cost + lower_bound >= tsp_cost_) continue;
988
989 // If next is the last node to visit, update tsp_cost_ and skip.
990 if (next_cardinality == num_nodes_) {
991 tsp_cost_ = next_cost;
992 continue;
993 }
994
995 // Add to the stack, finally.
996 mem_.SetValue(next_set, next_to_visit, next_cost);
997 state_stack.push(std::make_pair(next_set, next_to_visit));
998 }
999 }
1000
1001 solved_ = true;
1002}
1003
1004template <typename CostType, typename CostFunction>
1006 int end_node) {
1007 Solve(end_node);
1008 return tsp_cost_;
1009}
1010
1011template <typename CostType, typename CostFunction>
1012CostType
1014 NodeSet current_set, int last_visited) {
1015 // TODO(user): Compute MST.
1016 return 0; // For now, return 0 as future lower bound.
1017}
1018} // namespace operations_research
1019
1020#endif // OR_TOOLS_GRAPH_HAMILTONIAN_PATH_H_
const ElementIterator & operator++()
bool operator!=(const ElementIterator &other) const
std::vector< int > HamiltonianPath(int end_node)
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
static const int MaxCardinality
bool Contains(int n) const
ElementIterator< Set > end() const
Set RemoveElement(int n) const
int SingletonRank(Set singleton) const
int ElementRank(int n) const
bool operator!=(const Set &other) const
ElementIterator< Set > begin() const
static constexpr Integer One
bool Includes(Set other) const
static Set Singleton(Integer n)
static constexpr Integer Zero
static Set FullSet(Integer card)
Set AddElement(int n) const
bool operator!=(const SetRangeIterator &other) const
const SetRangeIterator & operator++()
SetRangeIterator< SetRangeWithCardinality > begin() const
SetRangeIterator< SetRangeWithCardinality > end() const
HamiltonianPathSolver< CostType, CostFunction > MakeHamiltonianPathSolver(int num_nodes, CostFunction cost)