dkh.graph.dijkstra source code

1 module dkh.graph.dijkstra;
2 
3 import dkh.algorithm, dkh.container.radixheap;
4 
5 struct DijkstraInfo(T) {
6     T[] dist;
7     int[] from;
8     this(int n, T inf) {
9         dist = new T[n]; dist[] = inf;
10         from = new int[n];
11     }
12 }
13 
14 DijkstraInfo!D dijkstra(D, T)(T g, int s, D inf = D.max) {
15     import std.conv : to;
16     import std.typecons : Tuple;
17     import std.container : make, Array, heapify;
18     import std.container.binaryheap : BinaryHeap;
19     import std.traits : isIntegral;
20 
21     int V = g.length.to!int;
22     auto dijk = DijkstraInfo!D(V, inf);
23     with (dijk) {        
24         alias P = Tuple!(int, "to", D, "dist");
25         auto q = (){
26             static if (isIntegral!D) {
27                 return RadixHeap!(P, "a.dist")();
28             } else {
29                 return heapify!"a.dist>b.dist"(make!(Array!P));
30             }
31         }();
32         q.insert(P(s, D(0)));
33         dist[s] = D(0);
34         from[s] = -1;
35         while (!q.empty) {
36             P p = q.front; q.removeFront();
37             if (dist[p.to] < p.dist) continue;
38             foreach (e; g[p.to]) {
39                 if (p.dist+e.dist < dist[e.to]) {
40                     dist[e.to] = p.dist+e.dist;
41                     from[e.to] = p.to;
42                     q.insert(P(e.to, dist[e.to]));
43                 }
44             }
45         }
46     }
47     return dijk;
48 }
49 
50 DijkstraInfo!D dijkstraDense(D, T)(T g, int s, D inf = D.max) {
51     import std.conv : to;
52     import std.typecons : Tuple;
53     import std.container : make, Array, heapify;
54     import std.range : enumerate;
55     import std.algorithm : filter;
56 
57     int V = g.length.to!int;
58     auto dijk = DijkstraInfo!D(V, inf);
59     with (dijk) {
60         alias P = Tuple!(int, "to", D, "dist");
61 
62         bool[] used = new bool[](V);
63         dist[s] = D(0);
64         from[s] = -1;
65         while (true) {
66             //todo can optimize
67             auto rng = dist.enumerate.filter!(a => !used[a.index]);
68             if (rng.empty) break;
69             auto nx = rng.minimum!"a.value < b.value";
70             used[nx.index] = true;
71             P p = P(nx.index.to!int, nx.value); 
72             if (p.dist == inf) continue;
73             foreach (e; g[p.to]) {
74                 if (p.dist+e.dist < dist[e.to]) {
75                     dist[e.to] = p.dist+e.dist;
76                     from[e.to] = p.to;
77                 }
78             }
79         }
80     }
81     return dijk;
82 }
83 
84 unittest {
85     import std.algorithm, std.conv, std.stdio, std.range;
86     import std.random;
87     import std.typecons;
88 
89     alias E = Tuple!(int, "to", int, "dist");
90 
91     void f(alias pred)() {
92         int n = uniform(1, 50);
93         int m = uniform(1, 500);
94         E[][] g = new E[][n];
95         int[][] dist = new int[][](n, n);
96         foreach (i, ref v; dist) {
97             v[] = 10^^9;
98         }
99         iota(n).each!(i => dist[i][i] = 0);
100         foreach (i; 0..m) {
101             int a = uniform(0, n);
102             int b = uniform(0, n);
103             int c = uniform(0, 1000);
104             g[a] ~= E(b, c);
105             dist[a][b] = min(dist[a][b], c);
106         }
107         foreach (k; 0..n) {
108             foreach (i; 0..n) {
109                 foreach (j; 0..n) {
110                     dist[i][j] = min(dist[i][j], dist[i][k]+dist[k][j]);
111                 }
112             }
113         }
114         foreach (i; 0..n) {
115             auto dijk = pred!int(g, i, 10^^9);
116             foreach (j; 0..n) {
117                 assert(dist[i][j] == dijk.dist[j]);
118             }
119             assert(dijk.from[i] == -1);
120             foreach (j; 0..n) {
121                 if (i == j) continue;
122                 if (dist[i][j] == 10^^9) continue;
123                 int b = dijk.from[j];
124                 assert(dijk.dist[j] == dist[i][b]+dist[b][j]);
125             }
126         }
127     }
128     import dkh.stopwatch;
129     auto ti = benchmark!(f!dijkstra, f!dijkstraDense)(100);
130     writeln("Dijkstra Random100: ", ti[].map!(a => a.toMsecs()));
131 }
132 
133 unittest {
134     import std.algorithm, std.conv, std.stdio, std.range;
135     import std.random;
136     import std.typecons;
137 
138     alias E = Tuple!(int, "to", int, "dist");
139 
140     void f(alias pred)() {
141         int n = uniform(1, 50);
142         int m = uniform(1, 500);
143         E[][] g = new E[][n];
144         int[][] dist = new int[][](n, n);
145         foreach (i, ref v; dist) {
146             v[] = int.max;
147         }
148         iota(n).each!(i => dist[i][i] = 0);
149         foreach (i; 0..m) {
150             int a = uniform(0, n);
151             int b = uniform(0, n);
152             int c = uniform(0, 1000);
153             g[a] ~= E(b, c);
154             dist[a][b] = min(dist[a][b], c);
155         }
156         foreach (k; 0..n) {
157             foreach (i; 0..n) {
158                 foreach (j; 0..n) {
159                     if (dist[i][k] == int.max) continue;
160                     if (dist[k][j] == int.max) continue;
161                     dist[i][j] = min(dist[i][j], dist[i][k]+dist[k][j]);
162                 }
163             }
164         }
165         foreach (i; 0..n) {
166             auto dijk = pred!int(g, i);
167             foreach (j; 0..n) {
168                 assert(dist[i][j] == dijk.dist[j]);
169             }
170             assert(dijk.from[i] == -1);
171             foreach (j; 0..n) {
172                 if (i == j) continue;
173                 if (dist[i][j] == int.max) continue;
174                 int b = dijk.from[j];
175                 assert(dijk.dist[j] == dist[i][b]+dist[b][j]);
176             }
177         }
178     }
179     import dkh.stopwatch;
180     auto ti = benchmark!(f!dijkstra, f!dijkstraDense)(100);
181     writeln("Dijkstra INF Random100: ", ti[].map!(a => a.toMsecs()));
182 }