PGROUTING  2.5
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pgr_tsp.hpp
Go to the documentation of this file.
1 /*PGR-GNU*****************************************************************
2  * File: tsp_driver.cpp
3  *
4  * Generated with Template by:
5  * Copyright (c) 2015 pgRouting developers
6  * Mail: project@pgrouting.org
7  *
8  * Function's developer:
9  * Copyright (c) 2015 Celia Virginia Vergara Castillo
10  * Mail: vicky_vergara@hotmail.com
11  *
12  * ------
13  *
14  * This program is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with this program; if not, write to the Free Software
26  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27  *
28  * ********************************************************************PGR-GNU*/
29 
30 #ifndef INCLUDE_TSP_PGR_TSP_HPP_
31 #define INCLUDE_TSP_PGR_TSP_HPP_
32 #pragma once
33 
34 #include <time.h>
35 
36 #include <sstream>
37 #include <vector>
38 #include <set>
39 #include <string>
40 #include <utility>
41 #include <limits>
42 #include <cmath>
43 
44 #include "cpp_common/pgr_assert.h"
45 #include "cpp_common/Dmatrix.h"
46 #include "tsp/eucledianDmatrix.h"
47 #include "tsp/tour.h"
48 
49 
50 
51 static
52 size_t
53 rand(size_t n) {
54  pgassert(n > 0);
55  return static_cast< size_t >(std::rand() % n);
56 }
57 
58 
59 static
60 size_t
61 succ(size_t i, size_t n) {
62  pgassert(n > 0);
63  return static_cast<size_t>((i + 1) % n);
64 }
65 
66 static
67 size_t
68 pred(size_t i, size_t n) {
69  return i == 0? n - 1: i -1;
70 }
71 
72 
73 namespace pgrouting {
74 namespace tsp {
75 
76 template < typename MATRIX >
77 class TSP: public MATRIX {
78  public:
79  using MATRIX::distance;
80  using MATRIX::tourCost;
81  using MATRIX::get_row;
82 
83  /*
84  * function members
85  */
86  explicit TSP(const MATRIX &_costs)
87  : MATRIX(_costs),
88  current_tour(_costs.size()),
89  best_tour(_costs.size()),
90  epsilon(0.000001),
91  n(_costs.size()),
92  updatecalls(0),
93  swap_count(0),
94  slide_count(0),
95  reverse_count(0),
96  improve_count(0) {
97  pgassert(n == MATRIX::size());
98  bestCost = MATRIX::tourCost(best_tour);
99  current_cost = MATRIX::tourCost(current_tour);
101  }
102 
103 
104  Tour get_tour() const {return best_tour;}
105 
106  std::string get_stats() const {
107  std::ostringstream log1;
108  log1
109  << "\nTotal swaps: " << swap_count
110  << "\nTotal slides: " << slide_count
111  << "\nTotal reverses: " << reverse_count
112  << "\nTimes best tour changed: " << improve_count;
113  return log1.str();}
114 
115  std::string get_log() const {
116  return log.str();}
117 
118  void greedyInitial(size_t idx_start = 0);
119  void annealing(
120  double initial_temperature,
121  double final_temperature,
122  double cooling_factor,
123  int64_t tries_per_temperature,
124  int64_t max_changes_per_temperature,
125  int64_t max_consecutive_non_changes,
126  bool randomize,
127  double time_limit);
128 
129 
130  private:
133  double bestCost;
134  double current_cost;
135  double epsilon;
136  size_t n;
137 
139 
140  std::ostringstream log;
141  size_t swap_count;
142  size_t slide_count;
145 
146  private:
147  void invariant() const;
148 
149  size_t find_closest_city(
150  size_t current_city,
151  const std::set<size_t> inserted) const;
152 
153  double getDeltaSlide(
154  size_t posP,
155  size_t posF,
156  size_t posL) const;
157 
158  void swapClimb();
159 
160  double getDeltaSwap(
161  size_t posA,
162  size_t posC) const;
163 
164  double getDeltaReverse(
165  size_t posA,
166  size_t posC) const;
167 
168  void update_if_best();
169 };
170 
171 
172 
173 template < typename MATRIX >
174 void TSP<MATRIX>::invariant() const {
175  /* the calculated value & the actual value are the same */
176  pgassert(std::fabs(tourCost(current_tour) - current_cost) < epsilon);
177  pgassert(std::fabs(tourCost(best_tour) - bestCost) < epsilon);
178  pgassert(n == MATRIX::ids.size());
179  pgassert(n == current_tour.size());
180  pgassert(n == best_tour.size());
181 }
182 
183 template < typename MATRIX >
184 void
186  invariant();
187  ++updatecalls;
188 
189  if (current_cost < bestCost) {
190  ++improve_count;
191  best_tour = current_tour;
192  bestCost = current_cost;
193  }
194 
195  invariant();
196 }
197 
198 
199 
200 template < typename MATRIX >
201 size_t
203  size_t current_city,
204  const std::set<size_t> inserted) const {
205  invariant();
206 
207  auto distance_row(get_row(current_city));
208  pgassert(distance_row.size() == n);
209 
210 #ifndef NDEBUG
211  std::ostringstream err;
212  for (const auto &d : distance_row) {
213  err << d << ", ";
214  }
215 #endif
216 
217  size_t best_city = 0;
218  auto best_distance = (std::numeric_limits<double>::max)();
219 #ifndef NDEBUG
220  bool found(false);
221 #endif
222 
223  for (size_t i = 0; i < distance_row.size(); ++i) {
224  if (i == current_city) continue;
225  if (inserted.find(i) != inserted.end()) continue;
226  if (distance_row[i] < best_distance) {
227  best_city = i;
228  best_distance = distance_row[i];
229 #ifndef NDEBUG
230  found = true;
231 #endif
232  }
233  }
234  pgassertwm(found, err.str());
235 
236  invariant();
237  return best_city;
238 }
239 
240 
241 
242 template < typename MATRIX >
243 void
244 TSP<MATRIX>::greedyInitial(size_t idx_start) {
245  invariant();
246 
247  std::set<size_t> pending(best_tour.cities.begin(), best_tour.cities.end());
248  std::set<size_t> inserted;
249  std::vector<size_t> tour_to_be;
250 
251  auto current_city = idx_start;
252 
253 #ifndef NDEBUG
254  std::ostringstream err;
255  auto ps(pending.size());
256 #endif
257 
258  pending.erase(idx_start);
259 
260 #ifndef NDEBUG
261  pgassert(pending.size() == (ps - 1));
262 #endif
263 
264  tour_to_be.push_back(current_city);
265  inserted.insert(current_city);
266 
267  while (!pending.empty()) {
268  auto next_city = find_closest_city(current_city, inserted);
269  tour_to_be.push_back(next_city);
270  inserted.insert(next_city);
271 
272 #ifndef NDEBUG
273  auto ps(pending.size());
274  err << "before";
275  for (const auto p : pending) {
276  err << p << ",";
277  }
278 #endif
279 
280  pending.erase(next_city);
281 
282 #ifndef NDEBUG
283  err << "\nafter deleting" << next_city << ":\t";
284  for (const auto p : pending) {
285  err << p << ",";
286  }
287  pgassertwm(pending.size() == (ps - 1), err.str());
288 #endif
289 
290  current_city = next_city;
291  }
292 
293  pgassert(tour_to_be.size() == n);
294  current_tour = Tour(tour_to_be);
295  current_cost = tourCost(current_tour);
296  update_if_best();
297  swapClimb();
298 
299  invariant();
300  return;
301 }
302 
303 
304 
305 /*
306  *
307  * 0 1 2 3 4 5 6 7 8 9
308  * p f l
309  * slides [4,5,6] to position p
310  *
311  * 0 1 4 5 6 2 3 7 8 9
312  *
313  *
314  * 0 1 2 3 4 5 6 7 8 9
315  * f l p
316  * slides [2,3,4] to position p
317  *
318  * 0 1 6 7 2 3 4 5 8 9
319  *
320  *
321  */
322 
323 template < typename MATRIX >
324 double
325 TSP<MATRIX>::getDeltaSlide(size_t place, size_t first, size_t last) const {
326  invariant();
327 #ifndef NDEBUG
328  std::ostringstream err;
329  err << "\tplace" << place
330  << "\tfirst" << first
331  << "\tlast" << last
332  << "\tn" << n;
333 #endif
334 
335  pgassertwm(place < first || place > last, err.str());
336  pgassertwm(first < last, err.str());
337  pgassertwm(last < n, err.str());
338  pgassertwm(place < n, err.str());
339  pgassertwm(first < n, err.str());
340 
341  /*
342  * Initial state
343  * [...f] [f+1 ... l] [l+1 ...p] [p+1 ...]
344  *
345  * final state
346  * [...f] [l+1 ... p] [f+1 ...l] [p+1 ...]
347  *
348  *
349  * Initial state
350  * [f+1 ... l]
351  * : :
352  * [...f] [l+1 ...p] [p+1 ...]
353  *
354  * final state
355  * [f+1 ... l]
356  * : :
357  * [...f] [l+1 ...p] [p+1 ...]
358  *
359  */
360 
361  auto cityP = current_tour.cities[place];
362  auto cityF = current_tour.cities[first];
363  auto cityL = current_tour.cities[last];
364  auto cityP1 = current_tour.cities[succ(place, n)];
365  auto cityF1 = current_tour.cities[succ(first, n)];
366  auto cityL1 = current_tour.cities[succ(last, n)];
367 
368  auto delta(
369  distance(cityF, cityL1)
370  + distance(cityP, cityF1)
371  + distance(cityL, cityP1)
372  - distance(cityF, cityF1)
373  - distance(cityL, cityL1)
374  - distance(cityP, cityP1));
375 
376 #ifndef NDEBUG
377  Tour new_tour(current_tour);
378  new_tour.slide(place, first, last);
379 
380  err << "\ncurrent_tour:";
381  for (const auto id : current_tour.cities) {
382  err << id << ", ";
383  }
384 
385  err << "\nnew_tour:";
386  for (const auto id : new_tour.cities) {
387  err << id << ", ";
388  }
389 
390  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
391  err << "\n"
392  << exactDelta
393  << " - " << delta
394  << " = "
395  << exactDelta - delta
396  << " = "
397  << std::fabs(exactDelta - delta);
398  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, err.str());
399 #endif
400 
401  invariant();
402  return delta;
403 }
404 
405 
406 /*
407  * c..d c..d
408  * | | => | |
409  * | | | |
410  * b -- a e --f b -- e a -- f
411  *
412  * a b 1 2 .. n-1 n c d
413  * a c n n-1 .. 2 1 c d
414  */
415 template < typename MATRIX >
416 double
417 TSP<MATRIX>::getDeltaSwap(size_t posA, size_t posE) const {
418  invariant();
419 
420  if (succ(posE, n ) == posA) std::swap(posA, posE);
421  if (succ(posA, n) == posE) {
422  auto b = current_tour.cities[pred(posA, n)];
423  auto a = current_tour.cities[posA];
424 
425  auto e = current_tour.cities[posE];
426  auto f = current_tour.cities[succ(posE, n)];
427  return distance(b, e) + distance(e, a) + distance(a, f)
428  - distance(b, a) - distance(a, e) - distance(e, f);
429  }
430 
431  auto b = current_tour.cities[pred(posA, n)];
432  auto a = current_tour.cities[posA];
433  auto c = current_tour.cities[succ(posA, n)];
434 
435  auto d = current_tour.cities[pred(posE, n)];
436  auto e = current_tour.cities[posE];
437  auto f = current_tour.cities[succ(posE, n)];
438 
439 #ifndef NDEBUG
440  auto delta = distance(b, e)
441  + distance(e, c) + distance(d, a) + distance(a, f)
442  - distance(b, a) - distance(a, c) - distance(d, e) - distance(e, f);
443  auto new_tour(current_tour);
444  new_tour.swap(posA, posE);
445  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
446  std::ostringstream log;
447  log << exactDelta
448  << " - " << delta
449  << " = "
450  << exactDelta - delta
451  << " = "
452  << std::fabs(exactDelta - delta);
453 
454  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, log.str());
455 #endif
456 
457  invariant();
458  return distance(b, e) + distance(e, c) + distance(d, a) + distance(a, f)
459  - distance(b, a) - distance(a, c) - distance(d, e) - distance(e, f);
460 }
461 
462 /*
463  * ..A C
464  * [ )
465  * ..a b 1 2 .. n-1 n c d ..
466  * ..a c n n-1 .. 2 1 b d ..
467  */
468 template < typename MATRIX >
469 double
470 TSP<MATRIX>::getDeltaReverse(size_t posA, size_t posC) const {
471  invariant();
472 
473  if (posA == (posC - 1)) return 0;
474  auto a = current_tour.cities[posA];
475  auto b = current_tour.cities[succ(posA, n)];
476 
477  auto c = current_tour.cities[posC];
478  auto d = current_tour.cities[succ(posC, n)];
479 
480 
481 #ifndef NDEBUG
482  auto delta =
483  distance(a, c) + distance(b, d) - distance(a, b) - distance(c, d);
484  auto new_tour(current_tour);
485  new_tour.reverse(posA, posC);
486  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
487 
488  std::ostringstream log;
489  log << "exactDelta(" << exactDelta
490  << ") - delta(" << delta
491  << ") = "
492  << exactDelta - delta
493  << " = "
494  << (exactDelta - delta)
495  << " epsilon = " << epsilon;
496  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, log.str());
497 #endif
498 
499  invariant();
500  return distance(a, c) + distance(b, d) - distance(a, b) - distance(c, d);
501 }
502 
503 template < typename MATRIX >
504 void
506  invariant();
507  pgassert(n > 2);
508 
509  // auto first = std::rand() % n;
510  // for (size_t first = std::rand() % n; first < n; first++) {
511  for (size_t first = 0; first < n; first++) {
512  for (size_t last = first + 1; last < n; last++) {
513  pgassert(first < last);
514 
515  auto energyChange = getDeltaSwap(first, last);
516 
517  if (energyChange < 0 && epsilon < std::fabs(energyChange)) {
518  ++swap_count;
519  current_cost += energyChange;
520  current_tour.swap(first, last);
521  update_if_best();
522  }
523  }
524  }
525  invariant();
526 }
527 
528 template < typename MATRIX >
529 void
531  double temperature,
532  double final_temperature,
533  double cooling_factor,
534  int64_t tries_per_temperature,
535  int64_t max_changes_per_temperature,
536  int64_t max_consecutive_non_changes,
537  bool randomize,
538  double time_limit) {
539  invariant();
540  clock_t start_time(clock());
541 
542  if (randomize) {
543  std::srand(static_cast<unsigned int>(time(NULL)));
544  } else {
545  std::srand(1);
546  }
547 
548 
549 
550 
551  /* annealing schedule */
552  for (; final_temperature < temperature; temperature *= cooling_factor) {
553  invariant();
554 
555  log << "\nCycle(" << temperature <<") ";
556 
557  /*
558  how many times the tour changed in current temperature
559  */
560  int64_t pathchg = 0;
561  size_t enchg = 0;
562  int64_t non_change = 0;
563  for (int64_t j = 0; j < tries_per_temperature; j++) {
564  ++non_change;
565 
566  auto which = rand(2);
567  // which = 1;
568  switch (which) {
569  case 0: {
570  /* reverse */
571  pgassert(n > 2);
572 
573  auto c1 = std::rand() % n;
574  auto c2 = std::rand() % n;
575 
576  if (c1 == c2) c2 = succ(c2, n);
577  if (c1 == (c2 - 1)) c2 = succ(c2, n);
578  if (c1 > c2) std::swap(c1, c2);
579 
580  pgassert(c1 != c2);
581  pgassert(c1 < n && c2 < n);
582  pgassert(c1 < c2);
583 
584  auto energyChange = getDeltaReverse(c1, c2);
585 
586  if ( (energyChange < 0
587  && epsilon < std::fabs(energyChange))
588  || (0 < energyChange
589  && (
590  static_cast<double>(std::rand()) /
591  static_cast<double>(RAND_MAX))
592  < exp(-energyChange / temperature))) {
593  if (energyChange < 0) ++enchg;
594  ++reverse_count;
595  ++pathchg;
596  non_change = 0;
597  current_cost += energyChange;
598  current_tour.reverse(c1, c2);
599  update_if_best();
600  }
601  }
602  break;
603  case 1: {
604  /* slide */
605  if (n <= 3) {
606  break;
607  }
608 
609  pgassert(n > 3);
610 
611  auto first = std::rand() % n;
612  auto last = std::rand() % n;
613 
614  if (first == last) last = succ(last, n);
615  if (first > last) std::swap(first, last);
616  if (first == 0 && last == (n - 1)) {
617  first = succ(first, n);
618  }
619 
620  pgassert((n - (last - first) - 1) > 0);
621  auto place = std::rand() % (n - (last - first) - 1);
622  place = place < first? place :
623  last + (place - first) + 1;
624 
625 
626  pgassert((place < first
627  || place > last)
628  && (first < last));
629 
630  auto energyChange = getDeltaSlide(
631  place, first, last);
632 
633  if ((energyChange < 0
634  && epsilon < std::fabs(energyChange))
635  || (0 < energyChange
636  && (static_cast<double>(std::rand())
637  / static_cast<double>(RAND_MAX))
638  < exp(-energyChange / temperature))) {
639  if (energyChange < 0) ++enchg;
640  ++slide_count;
641  ++pathchg;
642  non_change = 0;
643  current_cost += energyChange;
644  current_tour.slide(place, first, last);
645  update_if_best();
646  }
647  }
648  break;
649  } // switch
650 
651 
652  if (max_changes_per_temperature < pathchg
653  && max_consecutive_non_changes < non_change ) {
654  break;
655  }
656  } // for tries per temperature
657 
658  swapClimb();
659  clock_t current_time(clock());
660  double elapsed_time = static_cast<double>(
661  current_time - start_time) / CLOCKS_PER_SEC;
662  if (time_limit < elapsed_time) {
663  break;
664  }
665  log << "\ttotal changes =" << pathchg
666  << "\t" << enchg << " were because delta energy < 0";
667  if (pathchg == 0) break; /* if no change then quit */
668  } // for temperatures
669 }
670 
671 } // namespace tsp
672 } // namespace pgrouting
673 
674 #endif // INCLUDE_TSP_PGR_TSP_HPP_
double getDeltaSlide(size_t posP, size_t posF, size_t posL) const
Definition: pgr_tsp.cpp:224
#define pgassertwm(expr, msg)
Adds a message to the assertion.
Definition: pgr_assert.h:104
double getDeltaReverse(size_t posA, size_t posC) const
Definition: pgr_tsp.cpp:369
static size_t rand(size_t n)
Definition: pgr_tsp.hpp:53
double getDeltaSwap(size_t posA, size_t posC) const
Definition: pgr_tsp.cpp:316
void greedyInitial(size_t idx_start=0)
Definition: pgr_tsp.cpp:143
#define pgassert(expr)
Uses the standard assert syntax.
Definition: pgr_assert.h:81
size_t find_closest_city(size_t current_city, const std::set< size_t > inserted) const
Definition: pgr_tsp.cpp:101
void update_if_best()
Definition: pgr_tsp.cpp:84
void annealing(double initial_temperature, double final_temperature, double cooling_factor, int64_t tries_per_temperature, int64_t max_changes_per_temperature, int64_t max_consecutive_non_changes, bool randomize, double time_limit)
Definition: pgr_tsp.cpp:429
void invariant() const
Definition: pgr_tsp.cpp:73
std::string get_stats() const
Definition: pgr_tsp.hpp:106
static size_t succ(size_t i, size_t n)
Definition: pgr_tsp.hpp:61
Assertions Handling.
Tour get_tour() const
Definition: pgr_tsp.hpp:104
static size_t pred(size_t i, size_t n)
Definition: pgr_tsp.hpp:68
std::string get_log() const
Definition: pgr_tsp.hpp:115
std::ostringstream log
Definition: pgr_tsp.hpp:140
TSP(const MATRIX &_costs)
Definition: pgr_tsp.hpp:86