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