PGROUTING  2.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pgr_tsp.cpp
Go to the documentation of this file.
1 /*PGR-GNU*****************************************************************
2  * File: pgr_tsp.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 
31 #include "./pgr_tsp.hpp"
32 
33 #include <time.h>
34 
35 #include <utility>
36 #include <iomanip>
37 #include <limits>
38 #include <vector>
39 #include <set>
40 #include <algorithm>
41 #include <cmath>
42 
43 #include "../../common/src/pgr_types.h"
44 #include "../../common/src/pgr_assert.h"
45 
46 
47 
48 static
49 size_t
50 rand(size_t n) {
51  pgassert(n > 0);
52  return static_cast< size_t >(std::rand() % n);
53 }
54 
55 
56 static
57 size_t
58 succ(size_t i, size_t n) {
59  pgassert(n > 0);
60  return static_cast<size_t>((i + 1) % n);
61 }
62 
63 static
64 size_t
65 pred(size_t i, size_t n) {
66  return i == 0? n - 1: i -1;
67 }
68 
69 
70 namespace pgrouting {
71 namespace tsp {
72 
73 template < typename MATRIX >
74 void TSP<MATRIX>::invariant() const {
75  /* the calculated value & the actual value are the same */
76  pgassert(std::fabs(tourCost(current_tour) - current_cost) < epsilon);
77  pgassert(std::fabs(tourCost(best_tour) - bestCost) < epsilon);
78  pgassert(n == MATRIX::ids.size());
79  pgassert(n == current_tour.size());
80  pgassert(n == best_tour.size());
81 }
82 
83 template < typename MATRIX >
84 void
86  invariant();
87  ++updatecalls;
88 
89  if (current_cost < bestCost) {
90  ++improve_count;
91  best_tour = current_tour;
92  bestCost = current_cost;
93  }
94 
95  invariant();
96 }
97 
98 
99 
100 template < typename MATRIX >
101 size_t
103  size_t current_city,
104  const std::set<size_t> inserted) const {
105  invariant();
106 
107  auto distance_row(get_row(current_city));
108  pgassert(distance_row.size() == n);
109 
110 #ifndef NDEBUG
111  std::ostringstream err;
112  for (const auto &d : distance_row) {
113  err << d << ", ";
114  }
115 #endif
116 
117  size_t best_city = 0;
118  auto best_distance = (std::numeric_limits<double>::max)();
119 #ifndef NDEBUG
120  bool found(false);
121 #endif
122 
123  for (size_t i = 0; i < distance_row.size(); ++i) {
124  if (i == current_city) continue;
125  if (inserted.find(i) != inserted.end()) continue;
126  if (distance_row[i] < best_distance) {
127  best_city = i;
128  best_distance = distance_row[i];
129 #ifndef NDEBUG
130  found = true;
131 #endif
132  }
133  }
134  pgassertwm(found, err.str());
135 
136  invariant();
137  return best_city;
138 }
139 
140 
141 
142 template < typename MATRIX >
143 void
144 TSP<MATRIX>::greedyInitial(size_t idx_start) {
145  invariant();
146 
147  std::set<size_t> pending(best_tour.cities.begin(), best_tour.cities.end());
148  std::set<size_t> inserted;
149  std::vector<size_t> tour_to_be;
150 
151  auto current_city = idx_start;
152 
153 #ifndef NDEBUG
154  std::ostringstream err;
155  auto ps(pending.size());
156 #endif
157 
158  pending.erase(idx_start);
159 
160 #ifndef NDEBUG
161  pgassert(pending.size() == (ps - 1));
162 #endif
163 
164  tour_to_be.push_back(current_city);
165  inserted.insert(current_city);
166 
167  while (!pending.empty()) {
168  auto next_city = find_closest_city(current_city, inserted);
169  tour_to_be.push_back(next_city);
170  inserted.insert(next_city);
171 
172 #ifndef NDEBUG
173  auto ps(pending.size());
174  err << "before";
175  for (const auto p : pending) {
176  err << p << ",";
177  }
178 #endif
179 
180  pending.erase(next_city);
181 
182 #ifndef NDEBUG
183  err << "\nafter deleting" << next_city << ":\t";
184  for (const auto p : pending) {
185  err << p << ",";
186  }
187  pgassertwm(pending.size() == (ps - 1), err.str());
188 #endif
189 
190  current_city = next_city;
191  }
192 
193  pgassert(tour_to_be.size() == n);
194  current_tour = Tour(tour_to_be);
195  current_cost = tourCost(current_tour);
196  update_if_best();
197  swapClimb();
198 
199  invariant();
200  return;
201 }
202 
203 
204 
205 /*
206  *
207  * 0 1 2 3 4 5 6 7 8 9
208  * p f l
209  * slides [4,5,6] to position p
210  *
211  * 0 1 4 5 6 2 3 7 8 9
212  *
213  *
214  * 0 1 2 3 4 5 6 7 8 9
215  * f l p
216  * slides [2,3,4] to position p
217  *
218  * 0 1 6 7 2 3 4 5 8 9
219  *
220  *
221  */
222 
223 template < typename MATRIX >
224 double
225 TSP<MATRIX>::getDeltaSlide(size_t place, size_t first, size_t last) const {
226  invariant();
227 #ifndef NDEBUG
228  std::ostringstream err;
229  err << "\tplace" << place
230  << "\tfirst" << first
231  << "\tlast" << last
232  << "\tn" << n;
233 #endif
234 
235  pgassertwm(place < first || place > last, err.str());
236  pgassertwm(first < last, err.str());
237  pgassertwm(last < n, err.str());
238  pgassertwm(place < n, err.str());
239  pgassertwm(first < n, err.str());
240 
241  /*
242  * Initial state
243  * [...f] [f+1 ... l] [l+1 ...p] [p+1 ...]
244  *
245  * final state
246  * [...f] [l+1 ... p] [f+1 ...l] [p+1 ...]
247  *
248  *
249  * Initial state
250  * [f+1 ... l]
251  * : :
252  * [...f] [l+1 ...p] [p+1 ...]
253  *
254  * final state
255  * [f+1 ... l]
256  * : :
257  * [...f] [l+1 ...p] [p+1 ...]
258  *
259  */
260 
261  auto cityP = current_tour.cities[place];
262  auto cityF = current_tour.cities[first];
263  auto cityL = current_tour.cities[last];
264  auto cityP1 = current_tour.cities[succ(place, n)];
265  auto cityF1 = current_tour.cities[succ(first, n)];
266  auto cityL1 = current_tour.cities[succ(last, n)];
267 
268  auto delta(
269  distance(cityF, cityL1)
270  + distance(cityP, cityF1)
271  + distance(cityL, cityP1)
272  - distance(cityF, cityF1)
273  - distance(cityL, cityL1)
274  - distance(cityP, cityP1));
275 
276 #ifndef NDEBUG
277  Tour new_tour(current_tour);
278  new_tour.slide(place, first, last);
279 
280  err << "\ncurrent_tour:";
281  for (const auto id : current_tour.cities) {
282  err << id << ", ";
283  }
284 
285  err << "\nnew_tour:";
286  for (const auto id : new_tour.cities) {
287  err << id << ", ";
288  }
289 
290  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
291  err << "\n"
292  << exactDelta
293  << " - " << delta
294  << " = "
295  << exactDelta - delta
296  << " = "
297  << std::fabs(exactDelta - delta);
298  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, err.str());
299 #endif
300 
301  invariant();
302  return delta;
303 }
304 
305 
306 /*
307  * c..d c..d
308  * | | => | |
309  * | | | |
310  * b -- a e --f b -- e a -- f
311  *
312  * a b 1 2 .. n-1 n c d
313  * a c n n-1 .. 2 1 c d
314  */
315 template < typename MATRIX >
316 double
317 TSP<MATRIX>::getDeltaSwap(size_t posA, size_t posE) const {
318  invariant();
319 
320  if (succ(posE, n ) == posA) std::swap(posA, posE);
321  if (succ(posA, n) == posE) {
322  auto b = current_tour.cities[pred(posA, n)];
323  auto a = current_tour.cities[posA];
324 
325  auto e = current_tour.cities[posE];
326  auto f = current_tour.cities[succ(posE, n)];
327  return distance(b, e) + distance(e, a) + distance(a, f)
328  - distance(b, a) - distance(a, e) - distance(e, f);
329  }
330 
331  auto b = current_tour.cities[pred(posA, n)];
332  auto a = current_tour.cities[posA];
333  auto c = current_tour.cities[succ(posA, n)];
334 
335  auto d = current_tour.cities[pred(posE, n)];
336  auto e = current_tour.cities[posE];
337  auto f = current_tour.cities[succ(posE, n)];
338 
339 #ifndef NDEBUG
340  auto delta = distance(b, e)
341  + distance(e, c) + distance(d, a) + distance(a, f)
342  - distance(b, a) - distance(a, c) - distance(d, e) - distance(e, f);
343  auto new_tour(current_tour);
344  new_tour.swap(posA, posE);
345  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
346  std::ostringstream log;
347  log << exactDelta
348  << " - " << delta
349  << " = "
350  << exactDelta - delta
351  << " = "
352  << std::fabs(exactDelta - delta);
353 
354  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, log.str());
355 #endif
356 
357  invariant();
358  return distance(b, e) + distance(e, c) + distance(d, a) + distance(a, f)
359  - distance(b, a) - distance(a, c) - distance(d, e) - distance(e, f);
360 }
361 
362 /*
363  * ..A C
364  * [ )
365  * ..a b 1 2 .. n-1 n c d ..
366  * ..a c n n-1 .. 2 1 b d ..
367  */
368 template < typename MATRIX >
369 double
370 TSP<MATRIX>::getDeltaReverse(size_t posA, size_t posC) const {
371  invariant();
372 
373  if (posA == (posC - 1)) return 0;
374  auto a = current_tour.cities[posA];
375  auto b = current_tour.cities[succ(posA, n)];
376 
377  auto c = current_tour.cities[posC];
378  auto d = current_tour.cities[succ(posC, n)];
379 
380 
381 #ifndef NDEBUG
382  auto delta =
383  distance(a, c) + distance(b, d) - distance(a, b) - distance(c, d);
384  auto new_tour(current_tour);
385  new_tour.reverse(posA, posC);
386  auto exactDelta = tourCost(new_tour) - tourCost(current_tour);
387 
388  std::ostringstream log;
389  log << "exactDelta(" << exactDelta
390  << ") - delta(" << delta
391  << ") = "
392  << exactDelta - delta
393  << " = "
394  << (exactDelta - delta)
395  << " epsilon = " << epsilon;
396  pgassertwm(std::fabs((exactDelta - delta)) < epsilon, log.str());
397 #endif
398 
399  invariant();
400  return distance(a, c) + distance(b, d) - distance(a, b) - distance(c, d);
401 }
402 
403 template < typename MATRIX >
404 void
406  invariant();
407  pgassert(n > 2);
408 
409  // auto first = std::rand() % n;
410  // for (size_t first = std::rand() % n; first < n; first++) {
411  for (size_t first = 0; first < n; first++) {
412  for (size_t last = first + 1; last < n; last++) {
413  pgassert(first < last);
414 
415  auto energyChange = getDeltaSwap(first, last);
416 
417  if (energyChange < 0 && epsilon < std::fabs(energyChange)) {
418  ++swap_count;
419  current_cost += energyChange;
420  current_tour.swap(first, last);
421  update_if_best();
422  }
423  }
424  }
425  invariant();
426 }
427 
428 template < typename MATRIX >
429 void
431  double temperature,
432  double final_temperature,
433  double cooling_factor,
434  int64_t tries_per_temperature,
435  int64_t max_changes_per_temperature,
436  int64_t max_consecutive_non_changes,
437  bool randomize,
438  double time_limit) {
439  invariant();
440  clock_t start_time(clock());
441 
442  if (randomize) {
443  std::srand(static_cast<unsigned int>(time(NULL)));
444  } else {
445  std::srand(1);
446  }
447 
448 
449 
450 
451  /* annealing schedule */
452  for (; final_temperature < temperature; temperature *= cooling_factor) {
453  invariant();
454 
455  log << "\nCycle(" << temperature <<") ";
456 
457  /*
458  how many times the tour changed in current temperature
459  */
460  int64_t pathchg = 0;
461  size_t enchg = 0;
462  int64_t non_change = 0;
463  for (int64_t j = 0; j < tries_per_temperature; j++) {
464  ++non_change;
465 
466  auto which = rand(2);
467  // which = 1;
468  switch (which) {
469  case 0: {
470  /* reverse */
471  pgassert(n > 2);
472 
473  auto c1 = std::rand() % n;
474  auto c2 = std::rand() % n;
475 
476  if (c1 == c2) c2 = succ(c2, n);
477  if (c1 == (c2 - 1)) c2 = succ(c2, n);
478  if (c1 > c2) std::swap(c1, c2);
479 
480  pgassert(c1 != c2);
481  pgassert(c1 < n && c2 < n);
482  pgassert(c1 < c2);
483 
484  auto energyChange = getDeltaReverse(c1, c2);
485 
486  if ( (energyChange < 0
487  && epsilon < std::fabs(energyChange))
488  || (0 < energyChange
489  && (
490  static_cast<double>(std::rand()) /
491  static_cast<double>(RAND_MAX))
492  < exp(-energyChange / temperature))) {
493  if (energyChange < 0) ++enchg;
494  ++reverse_count;
495  ++pathchg;
496  non_change = 0;
497  current_cost += energyChange;
498  current_tour.reverse(c1, c2);
499  update_if_best();
500  }
501  }
502  break;
503  case 1: {
504  /* slide */
505  if (n <= 3) {
506  break;
507  }
508 
509  pgassert(n > 3);
510 
511  auto first = std::rand() % n;
512  auto last = std::rand() % n;
513 
514  if (first == last) last = succ(last, n);
515  if (first > last) std::swap(first, last);
516  if (first == 0 && last == (n - 1)) {
517  first = succ(first, n);
518  }
519 
520  pgassert((n - (last - first) - 1) > 0);
521  auto place = std::rand() % (n - (last - first) - 1);
522  place = place < first? place :
523  last + (place - first) + 1;
524 
525 
526  pgassert((place < first
527  || place > last)
528  && (first < last));
529 
530  auto energyChange = getDeltaSlide(
531  place, first, last);
532 
533  if ((energyChange < 0
534  && epsilon < std::fabs(energyChange))
535  || (0 < energyChange
536  && (static_cast<double>(std::rand())
537  / static_cast<double>(RAND_MAX))
538  < exp(-energyChange / temperature))) {
539  if (energyChange < 0) ++enchg;
540  ++slide_count;
541  ++pathchg;
542  non_change = 0;
543  current_cost += energyChange;
544  current_tour.slide(place, first, last);
545  update_if_best();
546  }
547  }
548  break;
549  } // switch
550 
551 
552  if (max_changes_per_temperature < pathchg
553  && max_consecutive_non_changes < non_change ) {
554  break;
555  }
556  } // for tries per temperature
557 
558  swapClimb();
559  clock_t current_time(clock());
560  double elapsed_time = static_cast<double>(
561  current_time - start_time) / CLOCKS_PER_SEC;
562  if (time_limit < elapsed_time) {
563  break;
564  }
565  log << "\ttotal changes =" << pathchg
566  << "\t" << enchg << " were because delta energy < 0";
567  if (pathchg == 0) break; /* if no change then quit */
568  } // for temperatures
569 }
570 
571 } // namespace tsp
572 } // namespace pgrouting
573 
double getDeltaSlide(size_t posP, size_t posF, size_t posL) const
Definition: pgr_tsp.cpp:225
#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:370
static size_t succ(size_t i, size_t n)
Definition: pgr_tsp.cpp:58
double getDeltaSwap(size_t posA, size_t posC) const
Definition: pgr_tsp.cpp:317
static size_t pred(size_t i, size_t n)
Definition: pgr_tsp.cpp:65
void greedyInitial(size_t idx_start=0)
Definition: pgr_tsp.cpp:144
#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:102
void update_if_best()
Definition: pgr_tsp.cpp:85
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:430
void invariant() const
Definition: pgr_tsp.cpp:74
std::vector< size_t > cities
Definition: tour.h:154
static size_t rand(size_t n)
Definition: pgr_tsp.cpp:50
void slide(size_t place, size_t first, size_t last)
Definition: tour.cpp:56