GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
sp-template.c
Go to the documentation of this file.
1 /* LIBDGL -- a Directed Graph Library implementation
2  * Copyright (C) 2002 Roberto Micarelli
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 /*
20  * best view with tabstop=4
21  */
22 
23 
24 /*
25  * SHORTEST PATH CACHE
26  *
27  * components:
28  * - start node id
29  * - visited network: a node is marked as visited when its departing
30  * edges have been added to the cache
31  * - predist network: node distances from start node
32  * - NodeHeap: holds unvisited nodes, the next node extracted is the
33  * unvisited node closest to SP start
34  *
35  * not all nodes in the predist network have been visited, SP from start
36  * is known only for visited nodes
37  * unvisited nodes can be reached, but not necessarily on the shortest
38  * possible path
39  * important for DGL_SP_CACHE_DISTANCE_FUNC and DGL_SP_CACHE_REPORT_FUNC
40  */
41 
42 #if !defined(DGL_DEFINE_TREE_PROCS) && !defined(DGL_DEFINE_FLAT_PROCS)
43 
45  dglInt32_t nStart)
46 {
47  pCache->nStartNode = nStart;
48  pCache->pvVisited = NULL;
49  pCache->pvPredist = NULL;
50  dglHeapInit(&pCache->NodeHeap);
51  if ((pCache->pvVisited =
54  return -1;
55  if ((pCache->pvPredist =
58  return -1;
59  return 0;
60 }
61 
63 {
64  if (pCache->pvVisited)
66  if (pCache->pvPredist)
68  dglHeapFree(&pCache->NodeHeap, NULL);
69 }
70 
71 
72 static int DGL_SP_CACHE_DISTANCE_FUNC(dglGraph_s * pgraph,
73  dglSPCache_s * pCache,
74  dglInt32_t * pnDistance,
75  dglInt32_t nStart,
76  dglInt32_t nDestination)
77 {
78  dglTreeTouchI32_s VisitedItem;
79  dglTreePredist_s *pPredistItem, PredistItem;
80 
81  if (pCache->nStartNode != nStart) {
83  return -pgraph->iErrno;
84  }
85 
86  VisitedItem.nKey = nDestination;
87  if (avl_find(pCache->pvVisited, &VisitedItem) == NULL) {
89  return -pgraph->iErrno;
90  }
91 
92  PredistItem.nKey = nDestination;
93  if ((pPredistItem = avl_find(pCache->pvPredist, &PredistItem)) == NULL) {
95  return -pgraph->iErrno;
96  }
97 
98  if (pnDistance)
99  *pnDistance = pPredistItem->nDistance;
100  return 0;
101 }
102 
104  dglSPCache_s * pCache,
105  dglInt32_t nStart,
106  dglInt32_t nDestination)
107 {
108  dglTreeTouchI32_s VisitedItem;
109  dglTreePredist_s *pPredistItem, PredistItem;
110  dglInt32_t *pEdge;
111  dglInt32_t *pDestination;
112  dglSPArc_s arc;
113  long i, istack = 0;
114  unsigned char *pstack = NULL;
115  unsigned char *ppop;
116  dglSPReport_s *pReport = NULL;
117 
118  if (pCache->nStartNode != nStart) {
120  return NULL;
121  }
122 
123  VisitedItem.nKey = nDestination;
124  if (avl_find(pCache->pvVisited, &VisitedItem) == NULL) {
126  return NULL;
127  }
128 
129  PredistItem.nKey = nDestination;
130  if (avl_find(pCache->pvPredist, &PredistItem) == NULL) {
132  return NULL;
133  }
134 
135  for (PredistItem.nKey = nDestination,
136  pPredistItem = avl_find(pCache->pvPredist, &PredistItem);
137  pPredistItem;
138  PredistItem.nKey = pPredistItem->nFrom,
139  pPredistItem = avl_find(pCache->pvPredist, &PredistItem)
140  ) {
141  if (pPredistItem->nFrom < 0) {
142  pgraph->iErrno = DGL_ERR_BadEdge;
143  goto spr_error;
144  }
145 
146  pEdge = (dglInt32_t *) pPredistItem->pnEdge;
147 
148  if (pPredistItem->bFlags == 0) {
149  if (pgraph->Flags & DGL_GS_FLAT) {
150  pDestination =
151  DGL_NODEBUFFER_SHIFT(pgraph,
152  DGL_EDGE_TAILNODE_OFFSET(pEdge));
153  }
154  else {
155  pDestination =
156  DGL_GET_NODE_FUNC(pgraph,
157  DGL_EDGE_TAILNODE_OFFSET(pEdge));
158  }
159  }
160  else {
161  if (pgraph->Flags & DGL_GS_FLAT) {
162  pDestination =
163  DGL_NODEBUFFER_SHIFT(pgraph,
164  DGL_EDGE_HEADNODE_OFFSET(pEdge));
165  }
166  else {
167  pDestination =
168  DGL_GET_NODE_FUNC(pgraph,
169  DGL_EDGE_HEADNODE_OFFSET(pEdge));
170  }
171  }
172 
173  if ((arc.pnEdge = DGL_EDGE_ALLOC(pgraph->EdgeAttrSize)) == NULL)
174  goto spr_error;
175  arc.nFrom = pPredistItem->nFrom;
176  arc.nTo = DGL_NODE_ID(pDestination);
177  arc.nDistance = pPredistItem->nDistance;
178  memcpy(arc.pnEdge, pEdge, DGL_EDGE_SIZEOF(pgraph->EdgeAttrSize));
179  DGL_EDGE_COST(arc.pnEdge) = pPredistItem->nCost;
180 
181  if ((pstack =
182  dgl_mempush(pstack, &istack, sizeof(dglSPArc_s),
183  &arc)) == NULL) {
185  goto spr_error;
186  }
187 
188  if (arc.nFrom == nStart)
189  break;
190  }
191 
192  if (pPredistItem == NULL) {
194  goto spr_error;
195  }
196 
197  if ((pReport = malloc(sizeof(dglSPReport_s))) == NULL) {
199  goto spr_error;
200  }
201  memset(pReport, 0, sizeof(dglSPReport_s));
202 
203  pReport->cArc = istack;
204 
205  if ((pReport->pArc = malloc(sizeof(dglSPArc_s) * pReport->cArc)) == NULL) {
207  goto spr_error;
208  }
209 
210  pReport->nDistance = 0;
211 
212  for (i = 0;
213  (ppop = dgl_mempop(pstack, &istack, sizeof(dglSPArc_s))) != NULL;
214  i++) {
215  memcpy(&pReport->pArc[i], ppop, sizeof(dglSPArc_s));
216  pReport->nDistance += DGL_EDGE_COST(pReport->pArc[i].pnEdge);
217  }
218 
219  pReport->nStartNode = nStart;
220  pReport->nDestinationNode = nDestination;
221 
222  if (pstack)
223  free(pstack);
224 
225  return pReport;
226 
227  spr_error:
228  if (pstack)
229  free(pstack);
230  if (pReport)
231  dglFreeSPReport(pgraph, pReport);
232 
233  return NULL;
234 }
235 #endif
236 
237 #if defined(DGL_DEFINE_TREE_PROCS) || defined(DGL_DEFINE_FLAT_PROCS)
238 
239 #define __EDGELOOP_BODY_1(f) \
240  if ( (f) == 0 ) { \
241  pDestination = _DGL_EDGE_TAILNODE(pgraph, pEdge); \
242  } \
243  else { \
244  pDestination = _DGL_EDGE_HEADNODE(pgraph, pEdge); \
245  } \
246  if ( !(DGL_NODE_STATUS(pDestination) & DGL_NS_TAIL) && pgraph->Version < 3) { \
247  pgraph->iErrno = DGL_ERR_BadEdge; \
248  goto sp_error; \
249  } \
250  clipOutput.nEdgeCost = DGL_EDGE_COST(pEdge); \
251  if ( fnClip ) { \
252  clipInput.pnPrevEdge = NULL; \
253  clipInput.pnNodeFrom = pStart; \
254  clipInput.pnEdge = pEdge; \
255  clipInput.pnNodeTo = pDestination; \
256  clipInput.nFromDistance = 0; \
257  if ( fnClip( pgraph , & clipInput , & clipOutput , pvClipArg ) ) continue; \
258  } \
259  pPredistItem = dglTreePredistAdd( pCache->pvPredist, DGL_NODE_ID(pDestination) ); \
260  if ( pPredistItem == NULL ) goto sp_error; \
261  pPredistItem->nFrom = nStart; \
262  pPredistItem->pnEdge = pEdge; \
263  pPredistItem->nCost = clipOutput.nEdgeCost; \
264  pPredistItem->nDistance = clipOutput.nEdgeCost; \
265  pPredistItem->bFlags = (f); \
266  heapvalue.pv = pEdge; \
267  if ( dglHeapInsertMin( & pCache->NodeHeap, pPredistItem->nDistance , f , heapvalue ) < 0 ) { \
268  pgraph->iErrno = DGL_ERR_HeapError; \
269  goto sp_error; \
270  }
271 
272 #define __EDGELOOP_BODY_2(f) \
273  if ( (f) == 0 ) { \
274  pDestination = _DGL_EDGE_TAILNODE(pgraph, pEdge); \
275  } \
276  else if ( pgraph->Version == 3 ) { \
277  pDestination = _DGL_EDGE_HEADNODE(pgraph, pEdge); \
278  } \
279  if ( !(DGL_NODE_STATUS(pDestination) & DGL_NS_TAIL) && pgraph->Version < 3) { \
280  pgraph->iErrno = DGL_ERR_BadEdge; \
281  goto sp_error; \
282  } \
283  clipOutput.nEdgeCost = DGL_EDGE_COST(pEdge); \
284  if ( fnClip ) { \
285  clipInput.pnPrevEdge = pEdge_prev; \
286  clipInput.pnNodeFrom = pStart; \
287  clipInput.pnEdge = pEdge; \
288  clipInput.pnNodeTo = pDestination; \
289  clipInput.nFromDistance = fromDist; \
290  if ( fnClip( pgraph , & clipInput , & clipOutput , pvClipArg ) ) continue; \
291  } \
292  findPredist.nKey = DGL_NODE_ID(pDestination); \
293  if ( (pPredistItem = avl_find( pCache->pvPredist, &findPredist)) == NULL ) { \
294  if ( (pPredistItem = dglTreePredistAdd( pCache->pvPredist, DGL_NODE_ID(pDestination) )) == NULL ) { \
295  pgraph->iErrno = DGL_ERR_MemoryExhausted; \
296  goto sp_error; \
297  } \
298  } \
299  else { \
300  if ( pPredistItem->nDistance <= fromDist + clipOutput.nEdgeCost ) { \
301  continue; \
302  } \
303  } \
304  pPredistItem->nFrom = DGL_NODE_ID(pStart); \
305  pPredistItem->pnEdge = pEdge; \
306  pPredistItem->nCost = clipOutput.nEdgeCost; \
307  pPredistItem->nDistance = fromDist + clipOutput.nEdgeCost; \
308  pPredistItem->bFlags = (f); \
309  heapvalue.pv = pEdge; \
310  if ( dglHeapInsertMin( & pCache->NodeHeap, pPredistItem->nDistance , f , heapvalue ) < 0 ) { \
311  pgraph->iErrno = DGL_ERR_HeapError; \
312  goto sp_error; \
313  }
314 
315 /*
316  * Dijkstra Shortest Path
317  */
318 int DGL_SP_DIJKSTRA_FUNC(dglGraph_s * pgraph,
319  dglSPReport_s ** ppReport,
320  dglInt32_t * pDistance,
321  dglInt32_t nStart,
322  dglInt32_t nDestination,
323  dglSPClip_fn fnClip,
324  void *pvClipArg, dglSPCache_s * pCache)
325 {
326  dglInt32_t *pStart; /* pointer to the start node (pgraph->pNodeBuffer) */
327  register dglInt32_t *pDestination; /* temporary destination pointer */
328  register dglInt32_t *pEdgeset; /* pointer to the edge (pgraph->pEdgeBuffer) */
329  register dglInt32_t *pEdge; /* pointer to the to-edges in edge */
330  register dglInt32_t *pEdge_prev; /* pointer to the previous edge in path */
331  int nRet;
333 
334  dglSPCache_s spCache;
335  int new_cache = 0;
336 
337  /*
338  * shortest path distance temporary min heap
339  */
340  dglHeapData_u heapvalue;
341  dglHeapNode_s heapnode;
342 
343  /*
344  * shortest path visited network
345  */
346  dglTreeTouchI32_s *pVisitedItem, findVisited;
347 
348  /*
349  * shortest path predecessor and distance network
350  */
351  dglTreePredist_s *pPredistItem, findPredist;
352 
353  /*
354  * args to clip()
355  */
356  dglSPClipInput_s clipInput;
357  dglSPClipOutput_s clipOutput;
358 
359 
360  /*
361  * Initialize the cache: initialize the heap and create temporary networks -
362  * The use of a predist network for predecessor and distance has two important results:
363  * 1) allows us not having to reset the whole graph status at each call;
364  * 2) use of a stack memory area for temporary (and otherwise possibly thread-conflicting) states.
365  * If a cache pointer was supplied, do not initialize it but try to get SP immediately.
366  */
367  if (pCache == NULL) {
368  pCache = &spCache;
369  DGL_SP_CACHE_INITIALIZE_FUNC(pgraph, pCache, nStart);
370  new_cache = 1;
371  }
372  else {
373  if (ppReport) {
374  if ((*ppReport =
375  DGL_SP_CACHE_REPORT_FUNC(pgraph, pCache, nStart,
376  nDestination)) != NULL) {
377  return 1;
378  }
379  }
380  else {
382  (pgraph, pCache, pDistance, nStart, nDestination) >= 0) {
383  return 2;
384  }
385  }
386  if (pgraph->iErrno == DGL_ERR_HeadNodeNotFound) {
387  DGL_SP_CACHE_RELEASE_FUNC(pgraph, pCache);
388  DGL_SP_CACHE_INITIALIZE_FUNC(pgraph, pCache, nStart);
389  new_cache = 1;
390  }
391  else if (pgraph->iErrno != DGL_ERR_TailNodeNotFound) {
392  goto sp_error;
393  }
394  }
395 
396  /*
397  * reset error status after using the cache
398  */
399  pgraph->iErrno = 0;
400 
401  if ((pStart = DGL_GET_NODE_FUNC(pgraph, nStart)) == NULL) {
403  goto sp_error;
404  }
405 
406  if ((pDestination = DGL_GET_NODE_FUNC(pgraph, nDestination)) == NULL) {
408  goto sp_error;
409  }
410 
411  if ((DGL_NODE_STATUS(pStart) & DGL_NS_ALONE) ||
412  (DGL_NODE_STATUS(pDestination) & DGL_NS_ALONE)) {
413  goto sp_error;
414  }
415 
416  if (!(DGL_NODE_STATUS(pStart) & DGL_NS_HEAD) && pgraph->Version < 3) {
417  goto sp_error;
418  }
419 
420  if (!(DGL_NODE_STATUS(pDestination) & DGL_NS_TAIL) && pgraph->Version < 3) {
421  goto sp_error;
422  }
423 
424  /* if we do not need a new cache, we just continue with the unvisited
425  * nodes in the cache */
426  if (new_cache) {
427  /*
428  * now we inspect all edges departing from the start node
429  * - at each loop 'pedge' points to the edge in the edge buffer
430  * - we invoke the caller's clip() and eventually skip the edge (clip() != 0)
431  * - we insert a item in the predist network to set actual predecessor and distance
432  * (there is no precedecessor at this stage) and actual distance from the starting node
433  * (at this stage it equals the edge's cost)
434  * - we insert a item in the node min-heap (sorted on node distance), storing the offset of the
435  * edge in the edge buffer.
436  * In the case of undirected graph (version 3) we inspect input edges as well.
437  */
438  pEdgeset = _DGL_OUTEDGESET(pgraph, pStart);
439  if (DGL_EDGESET_T_INITIALIZE_FUNC(pgraph, &laT, pEdgeset) < 0) {
440  goto sp_error;
441  }
442  for (pEdge = DGL_EDGESET_T_FIRST_FUNC(&laT);
443  pEdge; pEdge = DGL_EDGESET_T_NEXT_FUNC(&laT)
444  ) {
445  __EDGELOOP_BODY_1(0);
446  }
448 
449  if (pgraph->Version == 3) {
450  pEdgeset = _DGL_INEDGESET(pgraph, pStart);
451  if (DGL_EDGESET_T_INITIALIZE_FUNC(pgraph, &laT, pEdgeset) < 0) {
452  goto sp_error;
453  }
454  for (pEdge = DGL_EDGESET_T_FIRST_FUNC(&laT);
455  pEdge; pEdge = DGL_EDGESET_T_NEXT_FUNC(&laT)
456  ) {
457  if (DGL_EDGE_STATUS(pEdge) & DGL_ES_DIRECTED)
458  continue;
459  __EDGELOOP_BODY_1(1);
460  }
462  }
463  }
464 
465  /*
466  * Now we begin extracting nodes from the min-heap. Each node extracted is
467  * the one that is actually closest to the SP start.
468  */
469  while (dglHeapExtractMin(&pCache->NodeHeap, &heapnode) == 1) {
470  dglInt32_t fromDist;
471 
472  /*
473  * recover the stored edge pointer
474  */
475  pEdge = heapnode.value.pv;
476 
477  /*
478  * the new relative head is the tail of the edge
479  * or the head of the edge if the traversal was reversed (undirected edge)
480  */
481  if (heapnode.flags == 0) {
482  pStart = _DGL_EDGE_TAILNODE(pgraph, pEdge); /* continue from previous tail */
483  }
484  else {
485  pStart = _DGL_EDGE_HEADNODE(pgraph, pEdge); /* reversed head/tail */
486  }
487 
488  /*
489  * We do not want to explore twice the same node as a relative starting point,
490  * that's the meaning of 'visited'. We mark actual start node as 'visited' by
491  * inserting it into the visited-network. If we find actual node in the network
492  * we just give up and continue looping. Otherwise we add actual node to the network.
493  */
494  findVisited.nKey = DGL_NODE_ID(pStart);
495  if ((pVisitedItem =
496  avl_find(pCache->pvVisited, &findVisited)) == NULL) {
497  if (dglTreeTouchI32Add(pCache->pvVisited, DGL_NODE_ID(pStart)) ==
498  NULL) {
500  goto sp_error;
501  }
502  }
503 
504  /*
505  * Give up with visited nodes now
506  */
507  if (pVisitedItem) {
508  if (DGL_NODE_ID(pStart) == nDestination) {
509  /* should not happen but does not harm
510  * this case should have been handled above */
511  goto destination_found;
512  }
513  else
514  continue;
515  }
516 
517  /*
518  * If the node is not marked as having departing edges, then we are into a
519  * blind alley. Just give up this direction and continue looping.
520  * This only applies to v1 and v2 (digraphs)
521  */
522  if (!(DGL_NODE_STATUS(pStart) & DGL_NS_HEAD) && pgraph->Version < 3) {
523  if (DGL_NODE_ID(pStart) == nDestination) {
524  goto destination_found;
525  }
526  else
527  continue;
528  }
529 
530  /*
531  * save actual edge for later clip()
532  */
533  pEdge_prev = pEdge;
534 
535  /*
536  * Recover the head node distance from the predist network
537  */
538  findPredist.nKey = DGL_NODE_ID(pStart);
539  if ((pPredistItem =
540  avl_find(pCache->pvPredist, &findPredist)) == NULL) {
542  goto sp_error;
543  }
544 
545  fromDist = pPredistItem->nDistance;
546 
547  /*
548  * Loop on departing edges:
549  * Scan the edgeset and loads pedge at each iteration with next-edge.
550  * iWay == DGL_EDGESET_T_WAY_OUT then pedge is a out arc (departing from node) else ot is a in arc.
551  * V1 has no in-degree support so iWay is always OUT, V2/3 have in-degree support.
552  *
553  * This loop needs to be done also when destination is found, otherwise
554  * the node is marked as visited but its departing edges are not added to the cache
555  * --> loose end, we might need these edges later on
556  */
557  pEdgeset = _DGL_OUTEDGESET(pgraph, pStart);
558  if (DGL_EDGESET_T_INITIALIZE_FUNC(pgraph, &laT, pEdgeset) < 0) {
559  goto sp_error;
560  }
561  for (pEdge = DGL_EDGESET_T_FIRST_FUNC(&laT);
562  pEdge; pEdge = DGL_EDGESET_T_NEXT_FUNC(&laT)
563  ) {
564  __EDGELOOP_BODY_2(0);
565  }
567 
568  if (pgraph->Version == 3) {
569  pEdgeset = _DGL_INEDGESET(pgraph, pStart);
570  if (DGL_EDGESET_T_INITIALIZE_FUNC(pgraph, &laT, pEdgeset) < 0) {
571  goto sp_error;
572  }
573  for (pEdge = DGL_EDGESET_T_FIRST_FUNC(&laT);
574  pEdge; pEdge = DGL_EDGESET_T_NEXT_FUNC(&laT)
575  ) {
576  if (DGL_EDGE_STATUS(pEdge) & DGL_ES_DIRECTED)
577  continue;
578  __EDGELOOP_BODY_2(1);
579  }
581  }
582 
583  /*
584  * Dijkstra algorithm ends when the destination node is extracted from
585  * the min distance heap, that means: no other path exist in the network giving
586  * a shortest output.
587  * If this happens we jump to the epilogue in order to build a path report and return.
588  */
589  if (DGL_NODE_ID(pStart) == nDestination) {
590  goto destination_found;
591  }
592  }
593 
594  sp_error:
595  if (pCache == &spCache) {
596  DGL_SP_CACHE_RELEASE_FUNC(pgraph, pCache);
597  }
598  return -pgraph->iErrno; /* == 0 path not found */
599 
600  destination_found: /* path found - build a shortest path report or report the distance only */
601 
602  if (ppReport) {
603  *ppReport =
604  DGL_SP_CACHE_REPORT_FUNC(pgraph, pCache, nStart, nDestination);
605  if (*ppReport == NULL) {
606  nRet = -pgraph->iErrno;
607  }
608  else {
609  nRet = 1;
610  }
611  }
612  else {
614  (pgraph, pCache, pDistance, nStart, nDestination) < 0) {
615  nRet = -pgraph->iErrno;
616  }
617  else {
618  nRet = 2;
619  }
620  }
621  if (pCache == &spCache) {
622  DGL_SP_CACHE_RELEASE_FUNC(pgraph, pCache);
623  }
624  return nRet;
625 }
626 
627 #endif