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