pgRouting  2.2
pgRouting extends the PostGIS / PostgreSQL geospatial database to provide geospatial routing functionality.
 All Classes Functions Variables Pages
tsp.c
1 /*PGR-GNU*****************************************************************
2 
3 Copyright (c) 2015 pgRouting developers
4 Mail: project@pgrouting.org
5 
6 ------
7 
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12 
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
21 
22 ********************************************************************PGR-GNU*/
23 /*
24  * Traveling Salesman Problem solution algorithm for PostgreSQL
25  *
26  * Copyright (c) 2006 Anton A. Patrushev, Orkney, Inc.
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 2 of the License, or
31  * (at your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful,
34  * but WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36  * GNU General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
41  *
42  */
43 
44 #include "tsp.h"
45 
46 #include "postgres.h"
47 #include "executor/spi.h"
48 #include "funcapi.h"
49 #include "catalog/pg_type.h"
50 #if PGSQL_VERSION > 92
51 #include "access/htup_details.h"
52 #endif
53 
54 #include "string.h"
55 #include "math.h"
56 
57 #include "fmgr.h"
58 
59 #ifdef PG_MODULE_MAGIC
60 PG_MODULE_MAGIC;
61 #endif
62 
63 
64 // ------------------------------------------------------------------------
65 
66 /*
67  * Define this to have profiling enabled
68  */
69 //#define PROFILE
70 
71 #ifdef PROFILE
72 #include <sys/time.h>
73 
74 struct timeval prof_tsp, prof_store, prof_extract, prof_total;
75 long proftime[5];
76 long profipts1, profipts2, profopts;
77 #define profstart(x) do { gettimeofday(&x, NULL); } while (0);
78 #define profstop(n, x) do { struct timeval _profstop; \
79  long _proftime; \
80  gettimeofday(&_profstop, NULL); \
81  _proftime = ( _profstop.tv_sec*1000000+_profstop.tv_usec) - \
82  ( x.tv_sec*1000000+x.tv_usec); \
83  elog(NOTICE, \
84  "PRF(%s) %lu (%f ms)", \
85  (n), \
86  _proftime, _proftime / 1000.0); \
87  } while (0);
88 
89 #else
90 
91 #define profstart(x) do { } while (0);
92 #define profstop(n, x) do { } while (0);
93 
94 #endif // PROFILE
95 
96 // ------------------------------------------------------------------------
97 
98 Datum tsp(PG_FUNCTION_ARGS);
99 
100 #undef DEBUG
101 //#define DEBUG 1
102 
103 #ifdef DEBUG
104 #define DBG(format, arg...) \
105  elog(NOTICE, format , ## arg)
106 #else
107 #define DBG(format, arg...) do { ; } while (0)
108 #endif
109 
110 // The number of tuples to fetch from the SPI cursor at each iteration
111 #define TUPLIMIT 1000
112 
113 // macro to store distance values as DISTANCE[MAX_TOWNS][MAX_TOWNS]
114 #define D(i,j) DISTANCE[(i)*MAX_TOWNS + j]
115 
116 DTYPE *DISTANCE;
117 float *x;
118 float *y;
119 int MAX_TOWNS;
120 int total_tuples;
121 
122 
123 static char *
124 text2char(text *in)
125 {
126  char *out = (char*)palloc(VARSIZE(in));
127 
128  memcpy(out, VARDATA(in), VARSIZE(in) - VARHDRSZ);
129  out[VARSIZE(in) - VARHDRSZ] = '\0';
130  return out;
131 }
132 
133 static int
134 finish(int code, int ret)
135 {
136  code = SPI_finish();
137  if (code != SPI_OK_FINISH )
138  {
139  elog(ERROR,"couldn't disconnect from SPI");
140  return -1 ;
141  }
142  return ret;
143 }
144 
145 
146 typedef struct point_columns
147 {
148  int id;
149  float8 x;
150  float8 y;
152 
153 
154 static int
155 fetch_point_columns(SPITupleTable *tuptable, point_columns_t *point_columns)
156 {
157  //DBG("Fetching point");
158 
159  point_columns->id = SPI_fnumber(SPI_tuptable->tupdesc, "source_id");
160  point_columns->x = SPI_fnumber(SPI_tuptable->tupdesc, "x");
161  point_columns->y = SPI_fnumber(SPI_tuptable->tupdesc, "y");
162  if (point_columns->id == SPI_ERROR_NOATTRIBUTE ||
163  point_columns->x == SPI_ERROR_NOATTRIBUTE ||
164  point_columns->y == SPI_ERROR_NOATTRIBUTE)
165  {
166  elog(ERROR, "Error, query must return columns "
167  "'source_id', 'x' and 'y'");
168  return -1;
169  }
170 
171  //DBG("* Point %i [%f, %f]", point_columns->id, point_columns->x, point_columns->y);
172 
173  return 0;
174 }
175 
176 static void
177 fetch_point(HeapTuple *tuple, TupleDesc *tupdesc,
178  point_columns_t *point_columns, point_t *point)
179 {
180  Datum binval;
181  bool isnull;
182 
183  //DBG("inside fetch_point");
184 
185  binval = SPI_getbinval(*tuple, *tupdesc, point_columns->id, &isnull);
186  //DBG("Got id");
187 
188  if (isnull)
189  elog(ERROR, "id contains a null value");
190 
191  point->id = DatumGetInt32(binval);
192 
193  //DBG("id = %i", point->id);
194 
195  binval = SPI_getbinval(*tuple, *tupdesc, point_columns->x, &isnull);
196  if (isnull)
197  elog(ERROR, "x contains a null value");
198 
199  point->x = DatumGetFloat8(binval);
200 
201  //DBG("x = %f", point->x);
202 
203  binval = SPI_getbinval(*tuple, *tupdesc, point_columns->y, &isnull);
204 
205  if (isnull)
206  elog(ERROR, "y contains a null value");
207 
208  point->y = DatumGetFloat8(binval);
209 
210  //DBG("y = %f", point->y);
211 }
212 
213 
214 static int findid(point_t *points, int n, int id)
215 {
216  int i;
217 
218  if (!points || n==0) return 0;
219 
220  for (i=0; i<n; i++) {
221  if (points[i].id == id) return i;
222  }
223 
224  return 0;
225 }
226 
227 
228 static int solve_tsp(char* sql, char* p_ids,
229  int source, path_element_t** path)
230 {
231  int SPIcode;
232  void *SPIplan;
233  Portal SPIportal;
234  bool moredata = TRUE;
235  int ntuples;
236 
237  //todo replace path (vector of path_element_t) with this array
238  int *ids;
239 
240  point_t *points=NULL;
241  point_columns_t point_columns = {.id= -1, .x= -1, .y=-1};
242 
243  char *err_msg = NULL;
244  int ret = -1;
245 
246  char *p;
247  int z = 0;
248  int i;
249 
250  int tt, cc;
251  double dx, dy;
252  float fit=0.0;
253  int got_source = 0;
254 
255  DBG("inside tsp");
256 
257  //int total_tuples = 0;
258  total_tuples = 0;
259 
260  /* count the number of towns and allocate memory */
261  for (i=0, MAX_TOWNS=1; i<strlen(p_ids); i++)
262  if (p_ids[i] == ',') MAX_TOWNS++;
263 
264  DBG("MAX_TOWNS=%d", MAX_TOWNS);
265 
266  if (MAX_TOWNS < 4) {
267  elog(ERROR, "Error: TSP requires 4 or more locations. Only %d were supplied.", MAX_TOWNS);
268  return -1;
269  }
270 
271  *path = (path_element_t *) palloc( MAX_TOWNS * sizeof(path_element_t) );
272  if (! *path) {
273  elog(ERROR, "Failed to alloc memory!");
274  return -1;
275  }
276 
277  ids = (int *)palloc( MAX_TOWNS * sizeof(int) );
278  if (!ids) {
279  elog(ERROR, "Failed to alloc memory!");
280  return -1;
281  }
282 
283  DISTANCE = (DTYPE *) palloc(MAX_TOWNS * MAX_TOWNS * sizeof(DTYPE));
284  if (!DISTANCE) {
285  elog(ERROR, "Failed to alloc memory!");
286  return -1;
287  }
288 
289  x = (float *) palloc(MAX_TOWNS * sizeof(float));
290  y = (float *) palloc(MAX_TOWNS * sizeof(float));
291  if (!x || !y) {
292  elog(ERROR, "Failed to alloc memory!");
293  return -1;
294  }
295 
296  p = strtok(p_ids, ",");
297  while(p != NULL) {
298  ids[z]=atoi(p);
299  if (ids[z] == source) got_source = 1;
300  p = strtok(NULL, ",");
301  z++;
302  if(z > MAX_TOWNS)
303  {
304  elog(ERROR, "Number of points exeeds max number.");
305  break;
306  }
307  }
308 
309  if (!got_source) {
310  elog(ERROR, "tsp: source id not included in list of ids!");
311  return -1;
312  }
313 
314  DBG("start tsp");
315 
316  SPIcode = SPI_connect();
317 
318  if (SPIcode != SPI_OK_CONNECT) {
319  elog(ERROR, "tsp: couldn't open a connection to SPI");
320  return -1;
321  }
322 
323  SPIplan = SPI_prepare(sql, 0, NULL);
324 
325  if (SPIplan == NULL) {
326  elog(ERROR, "tsp: couldn't create query plan via SPI");
327  return finish(SPIcode, -1);
328  }
329 
330  if ((SPIportal = SPI_cursor_open(NULL, SPIplan, NULL, NULL, true)) == NULL) {
331  elog(ERROR, "tsp: SPI_cursor_open('%s') returns NULL", sql);
332  return finish(SPIcode, -1);
333  }
334 
335  DBG("Query: %s",sql);
336 
337  while (moredata == TRUE) {
338  SPI_cursor_fetch(SPIportal, TRUE, TUPLIMIT);
339 
340  if (point_columns.id == -1) {
341  if (fetch_point_columns(SPI_tuptable, &point_columns) == -1)
342  return finish(SPIcode, ret);
343  }
344 
345  ntuples = SPI_processed;
346 
347  total_tuples += ntuples;
348 
349  //DBG("Tuples: %i", total_tuples);
350 
351  if (!points)
352  points = palloc(total_tuples * sizeof(point_t));
353  else
354  points = repalloc(points, total_tuples * sizeof(point_t));
355 
356  if (points == NULL) {
357  elog(ERROR, "Out of memory");
358  return finish(SPIcode, ret);
359  }
360 
361  if (ntuples > 0) {
362  int t;
363  SPITupleTable *tuptable = SPI_tuptable;
364  TupleDesc tupdesc = SPI_tuptable->tupdesc;
365 
366  //DBG("Got tuple desc");
367 
368  for (t = 0; t < ntuples; t++) {
369  HeapTuple tuple = tuptable->vals[t];
370  fetch_point(&tuple, &tupdesc, &point_columns,
371  &points[total_tuples - ntuples + t]);
372  }
373 
374  SPI_freetuptable(tuptable);
375  }
376  else {
377  moredata = FALSE;
378  }
379  }
380 
381 
382  DBG("Calling TSP");
383 
384  profstop("extract", prof_extract);
385  profstart(prof_tsp);
386 
387  DBG("Total tuples: %i", total_tuples);
388 
389  for(tt=0; tt < total_tuples; tt++) {
390  ids[tt] = points[tt].id;
391  x[tt] = points[tt].x;
392  y[tt] = points[tt].y;
393 
394  //DBG("Point at %i: %i [%f, %f]", tt, ids[tt], x[tt], y[tt]);
395 
396  for(cc=0; cc < total_tuples; cc++) {
397  dx = x[tt] - x[cc];
398  dy = y[tt] - y[cc];
399  D(tt, cc) = D(cc, tt) = sqrt(dx*dx + dy*dy);
400  }
401  }
402 
403  DBG("DISTANCE matrix computed");
404 
405  ret = find_tsp_solution(total_tuples, DISTANCE, ids,
406  source, -1, &fit, err_msg);
407 
408  DBG("TSP solved!");
409  DBG("Score: %f", fit);
410 
411  if (ret < 0) {
412  elog(ERROR, "Error computing path: %s", err_msg);
413  return finish(SPIcode, ret);
414  }
415 
416  for(tt=0; tt < total_tuples; tt++) {
417  ((path_element_t*)(*path))[tt].vertex_id = ids[tt];
418  ((path_element_t*)(*path))[tt].edge_id = tt + 1;
419  if (tt == total_tuples-1)
420  ((path_element_t*)(*path))[tt].cost = D(findid(points, total_tuples, ids[tt]), findid(points, total_tuples, ids[0]));
421  else
422  ((path_element_t*)(*path))[tt].cost = D(findid(points, total_tuples, ids[tt]), findid(points, total_tuples, ids[tt+1]));
423  }
424 
425  pfree(points);
426 
427  profstop("tsp", prof_tsp);
428  profstart(prof_store);
429 
430  DBG("Profile changed and ret is %i", ret);
431 
432  if (ret < 0) {
433  //elog(ERROR, "Error computing path: %s", err_msg);
434  ereport(ERROR, (errcode(ERRCODE_E_R_E_CONTAINING_SQL_NOT_PERMITTED), errmsg("Error computing path: %s", err_msg)));
435  }
436 
437  return finish(SPIcode, ret);
438 }
439 
440 PG_FUNCTION_INFO_V1(tsp);
441 Datum
442 tsp(PG_FUNCTION_ARGS)
443 {
444  FuncCallContext *funcctx;
445  int call_cntr;
446  int max_calls;
447  TupleDesc tuple_desc;
448  path_element_t *path;
449 
450  /* stuff done only on the first call of the function */
451  if (SRF_IS_FIRSTCALL())
452  {
453  MemoryContext oldcontext;
454  //int path_count;
455  int ret=-1;
456 
457  // XXX profiling messages are not thread safe
458  profstart(prof_total);
459  profstart(prof_extract);
460 
461  /* create a function context for cross-call persistence */
462  funcctx = SRF_FIRSTCALL_INIT();
463 
464  /* switch to memory context appropriate for multiple function calls */
465  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
466 
467 
468 
469  ret = solve_tsp(text2char(PG_GETARG_TEXT_P(0)),
470  text2char(PG_GETARG_TEXT_P(1)),
471  PG_GETARG_INT32(2),
472  &path);
473 
474  /* total number of tuples to be returned */
475  DBG("Counting tuples number");
476 
477  funcctx->max_calls = total_tuples;
478 
479  funcctx->user_fctx = path;
480 
481  funcctx->tuple_desc = BlessTupleDesc(
482  RelationNameGetTupleDesc("pgr_costResult"));
483  MemoryContextSwitchTo(oldcontext);
484  }
485 
486  /* stuff done on every call of the function */
487  funcctx = SRF_PERCALL_SETUP();
488 
489  call_cntr = funcctx->call_cntr;
490  max_calls = funcctx->max_calls;
491  tuple_desc = funcctx->tuple_desc;
492 
493  path = (path_element_t *)funcctx->user_fctx;
494 
495  DBG("Trying to allocate some memory");
496  DBG("call_cntr = %i, max_calls = %i", call_cntr, max_calls);
497 
498  if (call_cntr < max_calls) /* do when there is more left to send */
499  {
500  HeapTuple tuple;
501  Datum result;
502  Datum *values;
503  char* nulls;
504 
505  values = palloc(4 * sizeof(Datum));
506  nulls = palloc(4 * sizeof(char));
507 
508  values[0] = Int32GetDatum(call_cntr);
509  nulls[0] = ' ';
510  values[1] = Int32GetDatum(path[call_cntr].vertex_id);
511  nulls[1] = ' ';
512  values[2] = Int32GetDatum(path[call_cntr].edge_id);
513  nulls[2] = ' ';
514  values[3] = Float8GetDatum(path[call_cntr].cost);
515  nulls[3] = ' ';
516 
517  DBG("Heap making");
518 
519  tuple = heap_formtuple(tuple_desc, values, nulls);
520 
521  DBG("Datum making");
522 
523  /* make the tuple into a datum */
524  result = HeapTupleGetDatum(tuple);
525 
526  DBG("VAL: %i, %i", values[0], result);
527  DBG("Trying to free some memory");
528 
529  /* clean up (this is not really necessary) */
530  pfree(values);
531  pfree(nulls);
532 
533 
534  SRF_RETURN_NEXT(funcctx, result);
535  }
536  else /* do when there is no more left */
537  {
538  DBG("Ending function");
539  profstop("store", prof_store);
540  profstop("total", prof_total);
541  DBG("Profiles stopped");
542 
543  SRF_RETURN_DONE(funcctx);
544  }
545 }
Definition: tsp.h:34