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