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
tsplib.c
Go to the documentation of this file.
1 /*
2  * THIS CODE IS DEPRECATED
3  */
4 
5 /*PGR-MIT*****************************************************************
6 
7 * $Id: tsplib.c, v 1.1 2006/05/13 23:39:56 woodbri Exp $
8 *
9 * tsplib
10 * Copyright 2005, 2013, Stephen Woodbridge, All rights Reserved
11 * This file is released under MIT-X license as part of pgRouting.
12 
13 ------
14 MIT/X license
15 
16 Permission is hereby granted, free of charge, to any person obtaining a copy
17 of this software and associated documentation files(the "Software"), to deal
18 in the Software without restriction, including without limitation the rights
19 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
20 copies of the Software, and to permit persons to whom the Software is
21 furnished to do so, subject to the following conditions:
22 
23 
24 The above copyright notice and this permission notice shall be included in
25 all copies or substantial portions of the Software.
26 
27 
28 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
29 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
31 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
32 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
33 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
34 THE SOFTWARE.
35 
36 ********************************************************************PGR-MIT*/
37 
38 /*
39 **********************************************************************
40 *
41 * tsplib.c
42 *
43 * Routine to calculate the node order for a Traveling Salesman Problem
44 * from a distance matrix.
45 *
46 * $Log: tsplib.c, v $
47 * Revision 1.1 2006/05/13 23:39:56 woodbri
48 * Initial checkin of IMT::TSP module.
49 *
50 * Initial revision created from router tsp2.c
51 *
52 **********************************************************************
53 *
54 * Simulated annealing and a non symmetric
55 * Euclidian Traveling Salesman Problem.
56 *
57 * Solution based on local search heuristics for
58 * non-crossing paths and nearest neighbors
59 *
60 * Storage Requirements: 4n ints
61 *
62 * Problem: given an nXn matrix of costs for n cities in the plane,
63 * find a permutation pi_1, pi_2, ..., pi_n of 1, 2, ..., n that
64 * minimizes sum for 1<=i<n D(pi_i, pi_i+1), where D(i, j) is the
65 * euclidian distance between cities i and j from the matrix
66 *
67 * This means that the ith row is the "from" city and the jth column
68 * is the cost to the "to" city.
69 *
70 * NOTE: costs are assumed to be ints, so scaling them to meters or
71 * some other int would be advisable.
72 *
73 * NOTE: the internal cost calculations ASSUME a symmetric matrix. There
74 * are notes of what needs to be updated, but all D(i, j)s need to be
75 * reviewed because order is important for asymetric cases.
76 *
77 * Note: with n cities, there is(n-1)!/2 possible tours.
78 * factorial(10)=3628800 factorial(50)=3E+64 factorial(150)=5.7E+262
79 * If we could check one tour per clock cycle on a 100 MHZ computer, we
80 * would still need to wait approximately 10^236 times the age of the
81 * universe to explore all tours for 150 cities.
82 *
83 * Derivied from code 1995 by Maugis Lionel(Sofreavia) on sid1
84 * that was placed in the public domain.
85 *************************************************************************/
86 
87 // #ifdef __MINGW64__
88 // #define ELOG_H
89 // #include <winsock2.h>
90 // #endif
91 #include <postgres.h>
92 #include <string.h> /* memcpy */
93 #include <math.h> /* exp */
94 
95 #include "./tsp.h"
96 
97 #undef DEBUG
98 // #define DEBUG 1
99 #include "../../common/src/debug_macro.h"
100 
101 
102 #define T_INIT 100
103 #define FINAL_T 0.1
104 #define COOLING 0.9 /* to lower down T(< 1) */
105 #define TRIES_PER_T 500*n
106 #define IMPROVED_PATH_PER_T 60*n
107 
108 /*
109  * MACHINE INDEPENDENT RANDOM NUMBER GENERATOR
110  * Written by: DIMACS (modified for TSP)
111 */
112 
113 #define PRANDMAX 1000000000
114 static int a;
115 static int b;
116 static int arr[55];
117 
118 static
119 int Rand(void);
120 
121 static
122 void initRand(int seed) {
123  int i, ii;
124  int last, next;
125 
126  seed %= PRANDMAX;
127  if (seed < 0) seed += PRANDMAX;
128 
129  arr[0] = last = seed;
130  next = 1;
131  for (i = 1; i < 55; i++) {
132  ii =(21 * i) % 55;
133  arr[ii] = next;
134  next = last - next;
135  if (next < 0)
136  next += PRANDMAX;
137  last = arr[ii];
138  }
139  a = 0;
140  b = 24;
141  for (i = 0; i < 165; i++)
142  last = Rand();
143 }
144 
145 static
146 int Rand(void) {
147  int t;
148 
149  if (a-- == 0)
150  a = 54;
151  if (b-- == 0)
152  b = 54;
153 
154  t = arr[a] - arr[b];
155 
156  if (t < 0)
157  t += PRANDMAX;
158 
159  arr[a] = t;
160 
161  return t;
162 }
163 
164 #define RREAL ((double)Rand()/PRANDMAX)
165 #define RANDOM Rand
166 #define unifRand(n) (Rand()%n)
167 
168 
169 /*
170  * Defs
171  */
172 
173 typedef int Path[3]; /* specify how to change path */
174 
175 typedef struct tspstruct {
176  int n;
180  int *iorder;
181  int *jorder;
182  int *border; // best order we find
183  float b[4];
184 } TSP;
185 
186 #define MOD(i, n) ((i) %(n) >= 0 ?(i) %(n) :(i) %(n) +(n))
187 #define D(x, y) dist[(x)*n+y]
188 
189 #define MIN(a, b)((a) < (b)?(a) : (b))
190 #define MAX(a, b)((a) > (b)?(a) : (b))
191 #define sqr(x) ((x)*(x))
192 
193 /*
194  * Prim's approximated TSP tour
195  * See also [Cristophides'92]
196  */
197 static
199  int *mst, *arc;
200  int i, j, k, l, a;
201  int n, *iorder, *jorder;
202  DTYPE d;
203  DTYPE maxd;
204  DTYPE *dist;
205  DTYPE *dis;
206 
207  jorder = tsp->jorder;
208  iorder = tsp->iorder;
209  dist = tsp->dist;
210  maxd = tsp->maxd;
211  n = tsp->n;
212 
213  if (!(mst =(int*) palloc((size_t) n * sizeof(int))) ||
214  !(arc =(int*) palloc((size_t) n * sizeof(int))) ||
215  !(dis =(DTYPE*) palloc((size_t) n * sizeof(DTYPE))) ) {
216  elog(ERROR, "Failed to allocate memory!");
217  return -1;
218  }
219  // PGR_DBG("findEulerianPath: 1");
220 
221  k = -1;
222  j = -1;
223  d = maxd;
224  dis[0] = -1;
225  for (i = 1; i < n; i++) {
226  dis[i] = D(i, 0);
227  arc[i] = 0;
228  if (d > dis[i]) {
229  d = dis[i];
230  j = i;
231  }
232  }
233  // PGR_DBG("findEulerianPath: j=%d", j);
234 
235  if (j == -1)
236  elog(ERROR, "Error TSP fail to findEulerianPath, check your distance matrix is valid.");
237 
238  /*
239  * O(n^2) Minimum Spanning Trees by Prim and Jarnick
240  * for graphs with adjacency matrix.
241  */
242  for (a = 0; a < n - 1; a++) {
243  mst[a] = j * n + arc[j]; /* join fragment j with MST */
244  dis[j] = -1;
245  d = maxd;
246  for (i = 0; i < n; i++) {
247  if (dis[i] >= 0) {
248  /* not connected yet */
249  if (dis[i] > D(i, j)) {
250  dis[i] = D(i, j);
251  arc[i] = j;
252  }
253  if (d > dis[i]) {
254  d = dis[i];
255  k = i;
256  }
257  }
258  }
259  j = k;
260  }
261  // PGR_DBG("findEulerianPath: 3");
262 
263  /*
264  * Preorder Tour of MST
265  */
266 #define VISITED(x) jorder[x]
267 #define NQ(x) arc[l++] = x
268 #define DQ() arc[--l]
269 #define EMPTY (l == 0)
270 
271  for (i = 0; i < n; i++) VISITED(i) = 0;
272  k = 0; l = 0; d = 0; NQ(0);
273  while (!EMPTY) {
274  i = DQ();
275  if (!VISITED(i)) {
276  iorder[k++] = i;
277  VISITED(i) = 1;
278  for (j = 0; j < n - 1; j++) /* push all kids of i */ {
279  if (i == mst[j]%n) NQ(mst[j]/n);
280  }
281  }
282  }
283  // PGR_DBG("findEulerianPath: 4");
284 
285  return 0;
286 }
287 
288 static
290  unsigned int i;
291  DTYPE len = 0;
292 
293  int *iorder = tsp->iorder;
294  DTYPE *dist = tsp->dist;
295  int n = tsp->n;
296 
297  for (i = 0; i < n-1; i++) {
298  len += D(iorder[i], iorder[i+1]);
299  }
300  len += D(iorder[n-1], iorder[0]); /* close path */
301  return(len);
302 }
303 
304 /*
305  * Local Search Heuristics
306  * b-------a b a
307  * . . => .\ /.
308  * . d...e . . e...d .
309  * ./ \. . .
310  * c f c-------f
311  */
312 static
314  int a, b, c, d, e, f;
315  int *iorder = tsp->iorder;
316  DTYPE *dist = tsp->dist;
317  int n = tsp->n;
318 
319  a = iorder[MOD(p[0]-1, n)];
320  b = iorder[p[0]];
321  c = iorder[p[1]];
322  d = iorder[MOD(p[1]+1, n)];
323  e = iorder[p[2]];
324  f = iorder[MOD(p[2]+1, n)];
325 
326  return(D(a, d) + D(e, b) + D(c, f) - D(a, b) - D(c, d) - D(e, f));
327  /* add cost between d and e if non symmetric TSP */
328 }
329 
330 static
331 void doThreeWay(TSP *tsp, Path p) {
332  int i, count, m1, m2, m3, a, b, c, d, e, f;
333  int *iorder = tsp->iorder;
334  int *jorder = tsp->jorder;
335  int n = tsp->n;
336 
337  a = MOD(p[0]-1, n);
338  b = p[0];
339  c = p[1];
340  d = MOD(p[1]+1, n);
341  e = p[2];
342  f = MOD(p[2]+1, n);
343 
344  m1 = MOD(n+c-b, n)+1; /* num cities from b to c */
345  m2 = MOD(n+a-f, n)+1; /* num cities from f to a */
346  m3 = MOD(n+e-d, n)+1; /* num cities from d to e */
347 
348  count = 0;
349  /* [b..c] */
350  for (i = 0; i < m1; i++) jorder[count++] = iorder[MOD(i+b, n)];
351 
352  /* [f..a] */
353  for (i = 0; i < m2; i++) jorder[count++] = iorder[MOD(i+f, n)];
354 
355  /* [d..e] */
356  for (i = 0; i < m3; i++) jorder[count++] = iorder[MOD(i+d, n)];
357 
358  /* copy segment back into iorder */
359  for (i = 0; i < n; i++) iorder[i] = jorder[i];
360 }
361 
362 /*
363  * c..b c..b
364  * \/ => | |
365  * /\ | |
366  * a d a d
367  */
368 static
370  int a, b, c, d;
371  int *iorder = tsp->iorder;
372  DTYPE *dist = tsp->dist;
373  int n = tsp->n;
374 
375  a = iorder[MOD(p[0]-1, n)];
376  b = iorder[p[0]];
377  c = iorder[p[1]];
378  d = iorder[MOD(p[1]+1, n)];
379 
380  return(D(d, b) + D(c, a) - D(a, b) - D(c, d));
381  /* add cost between c and b if non symmetric TSP */
382 }
383 static void doReverse(TSP *tsp, Path p) {
384  int i, nswaps, first, last, tmp;
385  int *iorder = tsp->iorder;
386  int n = tsp->n;
387 
388 
389  /* reverse path b...c */
390  nswaps =(MOD(p[1]-p[0], n)+1)/2;
391  for (i = 0; i < nswaps; i++) {
392  first = MOD(p[0]+i, n);
393  last = MOD(p[1]-i, n);
394  tmp = iorder[first];
395  iorder[first] = iorder[last];
396  iorder[last] = tmp;
397  }
398 }
399 
400 static
401 void annealing(TSP *tsp) {
402  Path p;
403  int i, j, pathchg;
404  int numOnPath, numNotOnPath;
405  DTYPE pathlen;
406  int n = tsp->n;
407  double energyChange, T;
408 
409  pathlen = pathLength(tsp);
410 
411  for (T = T_INIT; T > FINAL_T; T *= COOLING) /* annealing schedule */ {
412  pathchg = 0;
413  for (j = 0; j < TRIES_PER_T; j++) {
414  do {
415  p[0] = unifRand(n);
416  p[1] = unifRand(n);
417  /* non-empty path */
418  if (p[0] == p[1]) p[1] = MOD(p[0]+1, n);
419  numOnPath = MOD(p[1]-p[0], n) + 1;
420  numNotOnPath = n - numOnPath;
421  } while (numOnPath < 2 || numNotOnPath < 2); /* non-empty path */
422 
423  if (RANDOM() % 2) /* threeWay */ {
424  do {
425  p[2] = MOD(unifRand(numNotOnPath)+p[1]+1, n);
426  } while (p[0] == MOD(p[2]+1, n)); /* avoids a non-change */
427 
428  energyChange = getThreeWayCost(tsp, p);
429  if (energyChange < 0 || RREAL < exp(-energyChange/T)) {
430  pathchg++;
431  pathlen += energyChange;
432  doThreeWay(tsp, p);
433  }
434  } else {
435  /* path Reverse */
436  energyChange = getReverseCost(tsp, p);
437  if (energyChange < 0 || RREAL < exp(-energyChange/T)) {
438  pathchg++;
439  pathlen += energyChange;
440  doReverse(tsp, p);
441  }
442  }
443  // if the new length is better than best then save it as best
444  if (pathlen < tsp->bestlen) {
445  tsp->bestlen = pathlen;
446  for (i = 0; i < tsp->n; i++) tsp->border[i] = tsp->iorder[i];
447  }
448  if (pathchg > IMPROVED_PATH_PER_T) break; /* finish early */
449  }
450  PGR_DBG("T:%f L:%f B:%f C:%d", T, pathlen, tsp->bestlen, pathchg);
451  if (pathchg == 0) break; /* if no change then quit */
452  }
453 }
454 
455 
456 static
457 void reverse(int num, int *ids) {
458  int i, j, t;
459  for (i = 0, j = num - 1; i < j; i++, j--) {
460  t = ids[j];
461  ids[j] = ids[i];
462  ids[i] = t;
463  }
464 }
465 
466 
467 int find_tsp_solution(int num, DTYPE *cost, int *ids, int start, int end, DTYPE *total_len, char *err_msg) {
468  int i, j;
469  int istart = 0;
470  int jstart = 0;
471  int iend = -1;
472  int jend = -1;
473  int rev = 0;
474  TSP tsp;
475  int64_t seed = -314159L;
476  DTYPE blength;
477 
478  PGR_DBG("sizeof(int64_t)=%d", (int)sizeof(int64_t));
479 
480  initRand((int) seed);
481 
482 #ifdef DEBUG
483  char bufff[2048];
484  int nnn;
485  PGR_DBG("---------- Matrix[%d][%d] ---------------------\n", num, num);
486  for (i = 0; i < num; i++) {
487  sprintf(bufff, "%d:", i);
488  nnn = 0;
489  for (j = 0; j < num; j++) {
490  nnn += sprintf(bufff + nnn, "\t%.4f", cost[i * num + j]);
491  }
492  PGR_DBG("%s", bufff);
493  }
494 #endif
495 
496  /* initialize tsp struct */
497  tsp.n = num;
498  tsp.dist = NULL;
499  tsp.iorder = NULL;
500  tsp.jorder = NULL;
501  tsp.border = NULL;
502 
503  if (!(tsp.iorder =(int*) palloc((size_t) tsp.n * sizeof(int))) ||
504  !(tsp.jorder =(int*) palloc((size_t) tsp.n * sizeof(int))) ||
505  !(tsp.border =(int*) palloc((size_t) tsp.n * sizeof(int))) ) {
506  elog(FATAL, "Memory allocation failed!");
507  return -1;
508  }
509 
510  tsp.dist = cost;
511  tsp.maxd = 0;
512  for (i=0; i < tsp.n * tsp.n; i++) {
513  tsp.maxd = MAX(tsp.maxd, cost[i]);
514  }
515 
516  /* identity permutation */
517  for (i = 0; i < tsp.n; i++) tsp.iorder[i] = i;
518 
519  tsp.bestlen = pathLength(&tsp);
520  for (i = 0; i < tsp.n; i++) tsp.border[i] = tsp.iorder[i];
521 
522  PGR_DBG("Initial Path Length: %.4f", tsp.bestlen);
523 
524  /*
525  * Set up first eulerian path iorder to be improved by
526  * simulated annealing.
527  */
528  if (findEulerianPath(&tsp))
529  return -1;
530 
531  blength = pathLength(&tsp);
532  if (blength < tsp.bestlen) {
533  tsp.bestlen = blength;
534  for (i = 0; i < tsp.n; i++) tsp.border[i] = tsp.iorder[i];
535  }
536 
537  PGR_DBG("Approximated Path Length: %.4f", blength);
538 
539  annealing(&tsp);
540 
541  *total_len = pathLength(&tsp);
542  PGR_DBG("Final Path Length: %.4f", *total_len);
543 
544  *total_len = tsp.bestlen;
545  for (i = 0; i < tsp.n; i++) tsp.iorder[i] = tsp.border[i];
546  PGR_DBG("Best Path Length: %.4f", *total_len);
547 
548  // reorder ids[] with start as first
549 
550 #ifdef DEBUG
551  for (i = 0; i < tsp.n; i++) {
552  PGR_DBG("i: %d, ids[i]: %d, io[i]: %d, jo[i]: %d, jo[io[i]]: %d",
553  i, ids[i], tsp.iorder[i], tsp.jorder[i], tsp.jorder[tsp.iorder[i]]);
554  }
555 #endif
556 
557  // get index of start node in ids
558  for (i=0; i < tsp.n; i++) {
559  if (ids[i] == start) istart = i;
560  if (ids[i] == end) iend = i;
561  }
562  PGR_DBG("istart: %d, iend: %d", istart, iend);
563 
564  // get the idex of start in iorder
565  for (i=0; i < tsp.n; i++) {
566  if (tsp.iorder[i] == istart) jstart = i;
567  if (tsp.iorder[i] == iend) jend = i;
568  }
569  PGR_DBG("jstart: %d, jend: %d", jstart, jend);
570 
571  /*
572  * If the end is specified and the end point and it follow start
573  * then we swap start and end and extract the list backwards
574  * and later we reverse the list for the desired order.
575  */
576  if ((jend > 0 && jend == jstart+1) ||(jend == 0 && jstart == tsp.n-1)) {
577  int tmp = jend;
578  jend = jstart;
579  jstart = tmp;
580  rev = 1;
581  PGR_DBG("reversed start and end: jstart: %d, jend: %d", jstart, jend);
582  }
583 
584  // copy ids to tsp.jorder so we can rewrite ids
585  memcpy(tsp.jorder, ids, (size_t) tsp.n * sizeof(int));
586 
587  // write reordered ids into ids[]
588  // remember at this point jorder is our list if ids
589  for (i=jstart, j = 0; i < tsp.n; i++, j++)
590  ids[j] = tsp.jorder[tsp.iorder[i]];
591 
592  for (i=0; i < jstart; i++, j++)
593  ids[j] = tsp.jorder[tsp.iorder[i]];
594 
595  // if we reversed the order above, now put it correct.
596  if (rev) {
597  int tmp = jend;
598  jend = jstart;
599  jstart = tmp;
600  reverse(tsp.n, ids);
601  }
602 
603 #ifdef DEBUG
604  PGR_DBG("ids getting returned!");
605  for (i=0; i < tsp.n; i++) {
606  PGR_DBG("i: %d, ids[i]: %d, io[i]: %d, jo[i]: %d",
607  i, ids[i], tsp.iorder[i], tsp.jorder[i]);
608  }
609 #endif
610 
611  PGR_DBG("tsplib: jstart=%d, jend=%d, n=%d, j=%d", jstart, jend, tsp.n, j);
612 
613  return 0;
614 }
615 
616 /* EOF */
struct tspstruct TSP
int n
Definition: tsplib.c:176
int * iorder
Definition: tsplib.c:180
#define COOLING
Definition: tsplib.c:104
#define MOD(i, n)
Definition: tsplib.c:186
#define NQ(x)
static double getThreeWayCost(TSP *tsp, Path p)
Definition: tsplib.c:313
static void doThreeWay(TSP *tsp, Path p)
Definition: tsplib.c:331
#define DQ()
int * jorder
Definition: tsplib.c:181
#define PGR_DBG(...)
Definition: debug_macro.h:33
float b[4]
Definition: tsplib.c:183
#define D(x, y)
Definition: tsplib.c:187
#define RANDOM
Definition: tsplib.c:165
static int a
Definition: tsplib.c:114
#define T_INIT
Definition: tsplib.c:102
double maxd
Definition: tsplib.c:177
#define PRANDMAX
Definition: tsplib.c:113
#define VISITED(x)
double * dist
Definition: tsplib.c:178
static int findEulerianPath(TSP *tsp)
Definition: tsplib.c:198
#define MAX(a, b)
Definition: tsplib.c:190
static void doReverse(TSP *tsp, Path p)
Definition: tsplib.c:383
#define DTYPE
Definition: tsp.h:33
static int b
Definition: tsplib.c:115
#define TRIES_PER_T
Definition: tsplib.c:105
static double pathLength(TSP *tsp)
Definition: tsplib.c:289
int find_tsp_solution(int num, double *cost, int *ids, int start, int end, double *total_len, char *err_msg)
Definition: tsplib.c:467
static void annealing(TSP *tsp)
Definition: tsplib.c:401
#define RREAL
Definition: tsplib.c:164
#define IMPROVED_PATH_PER_T
Definition: tsplib.c:106
#define EMPTY
double bestlen
Definition: tsplib.c:179
#define FINAL_T
Definition: tsplib.c:103
static double getReverseCost(TSP *tsp, Path p)
Definition: tsplib.c:369
static int arr[55]
Definition: tsplib.c:116
char * err_msg
Definition: BDATester.cpp:50
static void reverse(int num, int *ids)
Definition: tsplib.c:457
static void initRand(int seed)
Definition: tsplib.c:122
static int Rand(void)
Definition: tsplib.c:146
int * border
Definition: tsplib.c:182
#define unifRand(n)
Definition: tsplib.c:166