1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the program and library             */
4    	/*         SCIP --- Solving Constraint Integer Programs                      */
5    	/*                                                                           */
6    	/*  Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SCIP; see the file LICENSE. If not visit scipopt.org.         */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file   dijkstra.c
26   	 * @ingroup OTHER_CFILES
27   	 * @brief  C implementation of Dijkstra's algorithm
28   	 * @author Thorsten Koch
29   	 * @author Marc Pfetsch
30   	 */
31   	
32   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
33   	
34   	#include <stdio.h>
35   	#include <stdlib.h>
36   	#include <assert.h>
37   	
38   	#include "dijkstra.h"
39   	
40   	
41   	/** Check whether the data structures of the graph are valid. */
42   	DIJKSTRA_Bool dijkstraGraphIsValid(
43   	   const DIJKSTRA_GRAPH* G                   /**< directed graph to be checked */
44   	   )
45   	{
46   	   unsigned int count = 0;
47   	   unsigned int i;
48   	   unsigned int k;
49   	
50   	   if ( G == NULL || G->outbeg == NULL || G->outcnt == NULL || G->weight == NULL || G->head == NULL )
51   	      abort();
52   	
53   	   for (i = 0; i < G->nodes; ++i)
54   	   {
55   	      for (k = G->outbeg[i]; k < G->outbeg[i] + G->outcnt[i]; ++k)
56   	      {
57   	         if ( G->head[k] >= G->nodes )
58   	            abort();
59   	
60   	         if ( G->weight[k] > G->maxweight || G->weight[k] < G->minweight )
61   	            abort();
62   	
63   	         ++count;
64   	      }
65   	      if ( G->head[k] != DIJKSTRA_UNUSED )
66   	         abort();
67   	
68   	      ++count;
69   	   }
70   	   if ( count > G->arcs )
71   	      abort();
72   	
73   	   return TRUE;
74   	}
75   	
76   	
77   	#ifndef NDEBUG
78   	/** Check whether heap is valid.
79   	 *
80   	 *  @note Sift up/down do not use order, only for the last the changed one is entered.
81   	 */
82   	static
83   	DIJKSTRA_Bool dijkstraHeapIsValid(
84   	   const unsigned int*   entry,              /**< entries of heap */
85   	   const unsigned long long* value,          /**< values in heap */
86   	   const unsigned int*   order,              /**< order of entries */
87   	   const unsigned int    used,               /**< number of used entries */
88   	   const unsigned int    size                /**< size of entry array */
89   	   )
90   	{
91   	   unsigned int i;
92   	
93   	   if ( entry == NULL || value == NULL || order == NULL || used  >  size )
94   	      return FALSE;
95   	
96   	   /* check heap property */
97   	   for (i = 0; i < used / 2; ++i)
98   	   {
99   	      if ( value[entry[i]] > value[entry[i + i]] )
100  	         return FALSE;
101  	      if ( i + i + 1 < used && value[entry[i]] > value[entry[i + i + 1]] )
102  	         return FALSE;
103  	   }
104  	
105  	   return TRUE;
106  	}
107  	#endif
108  	
109  	
110  	/** Moves an entry down in the vector until the sorting is valid again. */
111  	static
112  	void dijkstraSiftDown(
113  	   unsigned int*         entry,              /**< entries of heap */
114  	   const unsigned long long* value,          /**< values in heap */
115  	   unsigned int*         order,              /**< order of entries */
116  	   unsigned int          used,               /**< number of used entries */
117  	   unsigned int          current             /**< current entry to be sifted */
118  	   )
119  	{
120  	   unsigned long long val;
121  	   unsigned int child;
122  	   unsigned int ent;
123  	   unsigned int e;
124  	
125  	   child = current + current;
126  	   ent = entry[current];
127  	   val = value[ent];
128  	
129  	   while ( child < used )
130  	   {
131  	      e = entry[child];
132  	
133  	      if ( child + 1 < used )
134  	      {
135  	         if ( value[entry[child + 1]] < value[e] )
136  	         {
137  	            ++child;
138  	            e = entry[child];
139  	         }
140  	      }
141  	      if ( value[e] >= val )
142  	         break;
143  	
144  	      entry[current] = e;
145  	      order[e] = current;
146  	
147  	      current = child;
148  	      child += child;
149  	   }
150  	   entry[current] = ent;
151  	   order[ent] = current;
152  	}
153  	
154  	
155  	/** Moves an entry up in the vector until the sorting is valid again. */
156  	static
157  	void dijkstraSiftUp(
158  	   unsigned int*         entry,              /**< entries of heap */
159  	   const unsigned long long* value,          /**< values in heap */
160  	   unsigned int*         order,              /**< order of entries */
161  	   unsigned int          current             /**< current entry to be sifted */
162  	   )
163  	{
164  	   unsigned long long val;
165  	   unsigned int parent;
166  	   unsigned int ent;
167  	   unsigned int e;
168  	
169  	   ent = entry[current];
170  	   val = value[ent];
171  	
172  	   while ( current > 0 )
173  	   {
174  	      parent = current / 2;
175  	      e = entry[parent];
176  	
177  	      if ( value[e] <= val )
178  	         break;
179  	
180  	      entry[current] = e;
181  	      order[e] = current;
182  	      current = parent;
183  	   }
184  	   entry[current] = ent;
185  	   order[ent] = current;
186  	}
187  	
188  	
189  	/** Dijkstra's algorithm for shortest paths from a single source using binary heaps */
190  	unsigned int dijkstra(
191  	   const DIJKSTRA_GRAPH* G,                  /**< directed graph */
192  	   unsigned int          source,             /**< source node */
193  	   unsigned long long*   dist,               /**< node distances (allocated by user) */
194  	   unsigned int*         pred,               /**< node predecessors in final shortest path tree (allocated by user) */
195  	   unsigned int*         entry,              /**< temporary storage (for each node - must be allocated by user) */
196  	   unsigned int*         order               /**< temporary storage (for each node - must be allocated by user) */
197  	   )
198  	{
199  	   unsigned long long weight;
200  	   unsigned int iters = 0;
201  	   unsigned int used = 0;
202  	   unsigned int head;
203  	   unsigned int tail;
204  	   unsigned int i;
205  	   unsigned int e;
206  	
207  	   assert( dijkstraGraphIsValid(G) );
208  	   assert( source < G->nodes );
209  	   assert( dist != NULL );
210  	   assert( pred != NULL );
211  	
212  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
213  	
214  	   /* initialize nodes */
215  	   for (i = 0; i < G->nodes; ++i)
216  	   {
217  	      dist[i] = DIJKSTRA_FARAWAY;
218  	      order[i] = DIJKSTRA_UNUSED;
219  	      pred[i] = DIJKSTRA_UNUSED;
220  	   }
221  	
222  	   /* enter source node into heap */
223  	   entry[0] = source;
224  	   order[source] = 0;
225  	   pred[source] = DIJKSTRA_UNUSED;
226  	   dist[source] = 0;
227  	
228  	   ++used;
229  	
230  	   /* loop while heap is not empty */
231  	   while ( used > 0 )
232  	   {
233  	      /* get next node */
234  	      tail = entry[0];
235  	
236  	      /* remove node from heap */
237  	      --used;
238  	      entry[0] = entry[used];
239  	      order[entry[0]] = 0;
240  	      order[tail] = DIJKSTRA_UNUSED;
241  	
242  	      dijkstraSiftDown(entry, dist, order, used, 0);
243  	
244  	      assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
245  	      assert( entry[used] < G->nodes );
246  	
247  	      /* check adjacent nodes */
248  	      for (e = G->outbeg[tail]; G->head[e] != DIJKSTRA_UNUSED; ++e)
249  	      {
250  	         head = G->head[e];
251  	         weight = G->weight[e] + dist[tail];
252  	
253  		 /* Can we improve the current shortest path? */
254  	         if ( dist[head] > weight )
255  	         {
256  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
257  	            assert( used < G->nodes );
258  	            assert( head <= G->nodes );
259  	
260  	            pred[head] = tail;
261  	            dist[head] = weight;
262  	
263  	            if ( order[head] == DIJKSTRA_UNUSED )
264  	            {
265  	               assert( head < G->nodes );
266  	
267  	               entry[used] = head;
268  	               order[head] = used;
269  	
270  	               dijkstraSiftUp(entry, dist, order, used);
271  	               ++used;
272  	            }
273  	            else
274  	            {
275  	               dijkstraSiftUp(entry, dist, order, order[head]);
276  	            }
277  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
278  	
279  	            ++iters;
280  	         }
281  	      }
282  	   }
283  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
284  	
285  	   return iters;
286  	}
287  	
288  	
289  	/** Dijkstra's algorithm for shortest paths between a pair of nodes using binary heaps */
290  	unsigned int dijkstraPair(
291  	   const DIJKSTRA_GRAPH* G,                  /**< directed graph */
292  	   unsigned int          source,             /**< source node */
293  	   unsigned int          target,             /**< target node */
294  	   unsigned long long*   dist,               /**< node distances (allocated by user) */
295  	   unsigned int*         pred,               /**< node predecessors in final shortest path tree (allocated by user) */
296  	   unsigned int*         entry,              /**< temporary storage (for each node - must be allocated by user) */
297  	   unsigned int*         order               /**< temporary storage (for each node - must be allocated by user) */
298  	   )
299  	{
300  	   unsigned long long weight;
301  	   unsigned int iters = 0;
302  	   unsigned int used = 0;
303  	   unsigned int head;
304  	   unsigned int tail;
305  	   unsigned int i;
306  	   unsigned int e;
307  	
308  	   assert( dijkstraGraphIsValid(G) );
309  	   assert( source < G->nodes );
310  	   assert( target < G->nodes );
311  	   assert( dist != NULL );
312  	   assert( pred != NULL );
313  	
314  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
315  	
316  	   /* initialize nodes */
317  	   for (i = 0; i < G->nodes; ++i)
318  	   {
319  	      dist[i] = DIJKSTRA_FARAWAY;
320  	      order[i] = DIJKSTRA_UNUSED;
321  	      pred[i] = DIJKSTRA_UNUSED;
322  	   }
323  	
324  	   /* enter source node into heap */
325  	   entry[0] = source;
326  	   order[source] = 0;
327  	   pred[source] = DIJKSTRA_UNUSED;
328  	   dist[source] = 0;
329  	
330  	   ++used;
331  	
332  	   /* loop while heap is not empty */
333  	   while ( used > 0 )
334  	   {
335  	      /* get next node */
336  	      tail = entry[0];
337  	
338  	      /* stop if we have found the target node */
339  	      if ( tail == target )
340  		 break;
341  	
342  	      /* remove node from heap */
343  	      --used;
344  	      entry[0] = entry[used];
345  	      order[entry[0]] = 0;
346  	      order[tail] = DIJKSTRA_UNUSED;
347  	
348  	      dijkstraSiftDown(entry, dist, order, used, 0);
349  	
350  	      assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
351  	      assert( entry[used] < G->nodes );
352  	
353  	      /* check adjacent nodes */
354  	      for (e = G->outbeg[tail]; G->head[e] != DIJKSTRA_UNUSED; ++e)
355  	      {
356  	         head = G->head[e];
357  	         weight = G->weight[e] + dist[tail];
358  	
359  		 /* Can we improve the current shortest path? */
360  	         if ( dist[head] > weight )
361  	         {
362  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
363  	            assert( used < G->nodes );
364  	            assert( head <= G->nodes );
365  	
366  	            pred[head] = tail;
367  	            dist[head] = weight;
368  	
369  	            if ( order[head] == DIJKSTRA_UNUSED )
370  	            {
371  	               assert( head < G->nodes );
372  	
373  	               entry[used] = head;
374  	               order[head] = used;
375  	
376  	               dijkstraSiftUp(entry, dist, order, used);
377  	               ++used;
378  	            }
379  	            else
380  	            {
381  	               dijkstraSiftUp(entry, dist, order, order[head]);
382  	            }
383  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
384  	
385  	            ++iters;
386  	         }
387  	      }
388  	   }
389  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
390  	
391  	   return iters;
392  	}
393  	
394  	
395  	/** Dijkstra's algorithm for shortest paths between a pair of nodes using binary heaps and truncated at cutoff */
396  	unsigned int dijkstraPairCutoff(
397  	   const DIJKSTRA_GRAPH* G,                  /**< directed graph */
398  	   unsigned int          source,             /**< source node */
399  	   unsigned int          target,             /**< target node */
400  	   unsigned long long    cutoff,             /**< if the distance of a node reached this value, we truncate the search */
401  	   unsigned long long*   dist,               /**< node distances (allocated by user) */
402  	   unsigned int*         pred,               /**< node predecessors in final shortest path tree (allocated by user) */
403  	   unsigned int*         entry,              /**< temporary storage (for each node - must be allocated by user) */
404  	   unsigned int*         order               /**< temporary storage (for each node - must be allocated by user) */
405  	   )
406  	{
407  	   unsigned long long weight;
408  	   unsigned int iters = 0;
409  	   unsigned int used = 0;
410  	   unsigned int head;
411  	   unsigned int tail;
412  	   unsigned int i;
413  	   unsigned int e;
414  	
415  	   assert( dijkstraGraphIsValid(G) );
416  	   assert( source < G->nodes );
417  	   assert( target < G->nodes );
418  	   assert( dist != NULL );
419  	   assert( pred != NULL );
420  	
421  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
422  	
423  	   /* initialize nodes */
424  	   for (i = 0; i < G->nodes; ++i)
425  	   {
426  	      dist[i] = DIJKSTRA_FARAWAY;
427  	      order[i] = DIJKSTRA_UNUSED;
428  	      pred[i] = DIJKSTRA_UNUSED;
429  	   }
430  	
431  	   /* enter source node into heap */
432  	   entry[0] = source;
433  	   order[source] = 0;
434  	   pred[source] = DIJKSTRA_UNUSED;
435  	   dist[source] = 0;
436  	
437  	   ++used;
438  	
439  	   /* loop while heap is not empty */
440  	   while ( used > 0 )
441  	   {
442  	      /* get next node */
443  	      tail = entry[0];
444  	
445  	      /* stop if we have found the target node */
446  	      if ( tail == target )
447  		 break;
448  	
449  	      /* remove node from heap */
450  	      --used;
451  	      entry[0] = entry[used];
452  	      order[entry[0]] = 0;
453  	      order[tail] = DIJKSTRA_UNUSED;
454  	
455  	      dijkstraSiftDown(entry, dist, order, used, 0);
456  	
457  	      assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
458  	      assert( entry[used] < G->nodes );
459  	
460  	      /* only work on nodes if their distance is less than the cutoff */
461  	      if ( dist[tail] >= cutoff )
462  	         continue;
463  	
464  	      /* check adjacent nodes */
465  	      for (e = G->outbeg[tail]; G->head[e] != DIJKSTRA_UNUSED; ++e)
466  	      {
467  	         head = G->head[e];
468  	         weight = G->weight[e] + dist[tail];
469  	
470  		 /* Can we improve the current shortest path? */
471  	         if ( dist[head] > weight )
472  	         {
473  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
474  	            assert( used < G->nodes );
475  	            assert( head <= G->nodes );
476  	
477  	            pred[head] = tail;
478  	            dist[head] = weight;
479  	
480  	            if ( order[head] == DIJKSTRA_UNUSED )
481  	            {
482  	               assert( head < G->nodes );
483  	
484  	               entry[used] = head;
485  	               order[head] = used;
486  	
487  	               dijkstraSiftUp(entry, dist, order, used);
488  	               ++used;
489  	            }
490  	            else
491  	            {
492  	               dijkstraSiftUp(entry, dist, order, order[head]);
493  	            }
494  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
495  	
496  	            ++iters;
497  	         }
498  	      }
499  	   }
500  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
501  	
502  	   return iters;
503  	}
504  	
505  	
506  	/** Dijkstra's algorithm for shortest paths between a pair of nodes ignoring nodes, using binary heaps, and truncated at cutoff */
507  	unsigned int dijkstraPairCutoffIgnore(
508  	   const DIJKSTRA_GRAPH* G,                  /**< directed graph */
509  	   unsigned int          source,             /**< source node */
510  	   unsigned int          target,             /**< target node */
511  	   unsigned int*         ignore,             /**< marking nodes to be ignored (if value is nonzero) */
512  	   unsigned long long    cutoff,             /**< if the distance of a node reached this value, we truncate the search */
513  	   unsigned long long*   dist,               /**< node distances (allocated by user) */
514  	   unsigned int*         pred,               /**< node predecessors in final shortest path tree (allocated by user) */
515  	   unsigned int*         entry,              /**< temporary storage (for each node - must be allocated by user) */
516  	   unsigned int*         order               /**< temporary storage (for each node - must be allocated by user) */
517  	   )
518  	{
519  	   unsigned long long weight;
520  	   unsigned int iters = 0;
521  	   unsigned int used = 0;
522  	   unsigned int head;
523  	   unsigned int tail;
524  	   unsigned int i;
525  	   unsigned int e;
526  	
527  	   assert( dijkstraGraphIsValid(G) );
528  	   assert( source < G->nodes );
529  	   assert( target < G->nodes );
530  	   assert( dist != NULL );
531  	   assert( pred != NULL );
532  	   assert( ignore[source] == 0 );
533  	   assert( ignore[target] == 0 );
534  	
535  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
536  	
537  	   /* initialize nodes */
538  	   for (i = 0; i < G->nodes; ++i)
539  	   {
540  	      dist[i] = DIJKSTRA_FARAWAY;
541  	      order[i] = DIJKSTRA_UNUSED;
542  	      pred[i] = DIJKSTRA_UNUSED;
543  	   }
544  	
545  	   /* enter source node into heap */
546  	   entry[0] = source;
547  	   order[source] = 0;
548  	   pred[source] = DIJKSTRA_UNUSED;
549  	   dist[source] = 0;
550  	
551  	   ++used;
552  	
553  	   /* loop while heap is not empty */
554  	   while ( used > 0 )
555  	   {
556  	      /* get next node */
557  	      tail = entry[0];
558  	
559  	      /* stop if we have found the target node */
560  	      if ( tail == target )
561  		 break;
562  	
563  	      /* remove node from heap */
564  	      --used;
565  	      entry[0] = entry[used];
566  	      order[entry[0]] = 0;
567  	      order[tail] = DIJKSTRA_UNUSED;
568  	
569  	      dijkstraSiftDown(entry, dist, order, used, 0);
570  	
571  	      assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
572  	      assert( entry[used] < G->nodes );
573  	
574  	      /* only work on nodes if their distance is less than the cutoff */
575  	      if ( dist[tail] >= cutoff )
576  	         continue;
577  	
578  	      /* check adjacent nodes */
579  	      for (e = G->outbeg[tail]; G->head[e] != DIJKSTRA_UNUSED; ++e)
580  	      {
581  	         head = G->head[e];
582  	
583  	         /* skip ignored nodes */
584  	         if ( ignore[head] != 0 )
585  	            continue;
586  	
587  	         weight = G->weight[e] + dist[tail];
588  	
589  		 /* Can we improve the current shortest path? */
590  	         if ( dist[head] > weight )
591  	         {
592  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
593  	            assert( used < G->nodes );
594  	            assert( head <= G->nodes );
595  	
596  	            pred[head] = tail;
597  	            dist[head] = weight;
598  	
599  	            if ( order[head] == DIJKSTRA_UNUSED )
600  	            {
601  	               assert( head < G->nodes );
602  	
603  	               entry[used] = head;
604  	               order[head] = used;
605  	
606  	               dijkstraSiftUp(entry, dist, order, used);
607  	               ++used;
608  	            }
609  	            else
610  	            {
611  	               dijkstraSiftUp(entry, dist, order, order[head]);
612  	            }
613  	            assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
614  	
615  	            ++iters;
616  	         }
617  	      }
618  	   }
619  	   assert( dijkstraHeapIsValid(entry, dist, order, used, G->nodes) );
620  	
621  	   return iters;
622  	}
623