dkh.graph.mincostflow source code

1 module dkh.graph.mincostflow;
2 
3 import dkh.container.deque;
4 
5 /// 最小費用流の情報
6 struct MinCostFlowInfo(C, D, D EPS, T) {
7     T g;
8     int s, t;
9     C nc, capFlow; /// 今の最短路の容量, 今流した量
10     D nd, flow; /// 今の最短路の長さ, 今流したコスト
11     D[] dual; /// 双対問題の答え(=ポテンシャル)
12     private int[] pv, pe;
13     this(T g, int s, int t) {
14         this.g = g;
15         this.s = s;
16         this.t = t;
17         flow = D(0);
18         dual = new D[g.length]; dual[] = D(0);
19         pv = new int[g.length];
20         pe = new int[g.length];
21     }
22 }
23 
24 /// 最小費用流
25 MinCostFlowInfo!(C, D, EPS, T) minCostFlow(C, D, D EPS, T)(T g, int s, int t, bool neg = false) {
26     assert(s != t);
27     auto mcfInfo = MinCostFlowInfo!(C, D, EPS, T)(g, s, t);
28     mcfInfo.dualRef(neg);
29     return mcfInfo;
30 }
31 
32 ///
33 unittest {
34     import std.conv : to;
35     struct Edge {
36         int to, cap, dist, rev;
37     }
38 
39     void addEdge(Edge[][] g, int from, int to, int cap, int dist) {
40         g[from] ~= Edge(to, cap, dist, g[to].length.to!int);
41         g[to] ~= Edge(from, 0, -dist, g[from].length.to!int-1);
42     }
43 
44     auto g = new Edge[][](4);
45 
46     addEdge(g, 0, 1, 10, 3);
47     addEdge(g, 0, 2, 12, 3);
48     addEdge(g, 1, 3, 3, 2);
49     addEdge(g, 2, 3, 20, 4);
50     
51     auto mcfInfo = minCostFlow!(int, int, 0)(g, 0, 3, false);
52     //最初は 0->1->3で容量3, 距離5を流せる
53     assert(mcfInfo.nc == 3 && mcfInfo.nd == 5);
54 
55     //最短経路が変わらない間(=容量3)流す
56     mcfInfo.singleFlow(10^^9);
57 
58     assert(mcfInfo.capFlow == 3 && mcfInfo.flow == 15);
59 
60     //次は 0->2->3で容量12, 距離7を流せる
61     assert(mcfInfo.nc == 12 && mcfInfo.nd == 7);
62     
63     //最短経路が変わらない間(=容量12)流す
64     mcfInfo.singleFlow(10^^9);
65 
66     assert(mcfInfo.capFlow == 3 + 12);
67     assert(mcfInfo.flow == 15 + 12*7);
68 }
69 
70 ///min(nc, c)流す
71 C singleFlow(C, D, alias EPS, T)(ref MinCostFlowInfo!(C, D, EPS, T) mcfInfo, C c) {
72     import std.algorithm;
73     with (mcfInfo) {
74         if (nd == D.max) return nc;
75         c = min(c, nc);
76         for (int v = t; v != s; v = pv[v]) {
77             g[pv[v]][pe[v]].cap -= c;
78             g[v][g[pv[v]][pe[v]].rev].cap += c;
79         }        
80         capFlow += c;
81         flow += c * nd;
82         nc -= c;
83         if (!nc) dualRef(mcfInfo, false);
84     }
85     return c;
86 }
87 
88 ///流量がcになるまで流せるだけ流し続ける
89 void manyFlow(C, D, alias EPS, T)(ref MinCostFlowInfo!(C, D, EPS, T) mcfInfo, C c) {
90     with (mcfInfo) {
91         while (c) {
92             C f = singleFlow(mcfInfo, c);
93             if (!f) break;
94             c -= f;
95         }
96     }
97 }
98 
99 import dkh.container.stackpayload;
100 import dkh.container.radixheap;
101 
102 void dualRef(bool neg, C, D, alias EPS, T)(ref MinCostFlowInfo!(C, D, EPS, T) mcfInfo) {
103     import std.conv : to;
104     import std.traits : isIntegral;
105     import std.typecons;
106     import std.container;
107     import std.algorithm;
108     alias P = Tuple!(int, "to", D, "dist");
109 
110     with(mcfInfo) {
111         int n = g.length.to!int;
112         D[] dist = new D[n]; dist[] = D.max;
113         pv[] = -1; pe[] = -1;
114         StackPayload!int refV;
115         auto que = (){
116             static if (!neg) {
117                 static if (isIntegral!D) {
118                     return RadixHeap!(P, "a.dist")();
119                 } else {
120                     return heapify!"a.dist>b.dist"(make!(Array!P));
121                 }
122             } else {
123                 return Deque!P();
124             }
125         }();
126         void insert(P p) {
127             static if (!neg) {
128                 que.insert(p);
129             } else {
130                 que.insertBack(p);
131             }
132         }
133         P pop() {
134             P p;
135             static if (!neg) {
136                 p = que.front();
137                 que.removeFront();
138             } else {
139                 p = que.back();
140                 que.removeBack();
141             }
142             return p;
143         }
144         insert(P(s, D(0)));
145         dist[s] = D(0);
146         while (!que.empty) {
147             P p = pop();
148             int v = p.to;
149             if (dist[v] < p.dist) continue;
150             if (!neg) {
151                 if (v == t) break;
152                 refV ~= v;
153             }
154             foreach (int i, e; g[v]) {
155                 D ed = e.dist + dual[v] - dual[e.to];
156                 if (e.cap && dist[e.to] > dist[v] + ed + EPS) {
157                     dist[e.to] = dist[v] + ed;
158                     pv[e.to] = v; pe[e.to] = i;
159                     insert(P(e.to, dist[e.to]));
160                 }
161             }
162         }
163         if (dist[t] == D.max) {
164             nd = D.max;
165             nc = 0;
166             return;
167         }
168         static if (!neg) {
169             foreach (v; refV.data) {
170                 if (dist[v] >= dist[t]) continue;
171                 dual[v] += dist[v]-dist[t];
172             }
173         } else {
174             for (int v = 0; v < n; v++) {
175                 if (dist[v] == D.max) dual[v] = D.max;
176                 else dual[v] += dist[v];
177             }
178         }
179         
180         nd = dual[t]-dual[s];
181         nc = C.max;
182         for (int v = t; v != s; v = pv[v]) {
183             nc = min(nc, g[pv[v]][pe[v]].cap);
184         }
185     }
186 }
187 
188 void dualRef(C, D, alias EPS, T)(ref MinCostFlowInfo!(C, D, EPS, T) mcfInfo, bool neg) {
189     if (neg == false) {
190         dualRef!false(mcfInfo);
191     } else {
192         dualRef!true(mcfInfo);
193     }
194 }
195 
196 unittest {
197     import std.algorithm, std.conv, std.stdio, std.range;
198     import std.random;
199     import std.typecons;
200 
201     struct E {
202         int to, cap, dist, rev;
203     }
204     void addEdge(E[][] g, int from, int to, int cap, int dist) {
205         g[from] ~= E(to, cap, dist, g[to].length.to!int);
206         g[to] ~= E(from, 0, -dist, g[from].length.to!int-1);
207     }
208 
209     void f(bool neg)() {
210         int n = uniform(2, 20);
211         int m = uniform(0, 200);
212         int s, t;
213         while (true) {
214             s = uniform(0, n);
215             t = uniform(0, n);
216             if (s != t) break;
217         }
218         auto g = new E[][n];
219         E[][] elist = new E[][n];
220 
221         foreach (i; 0..m) {
222             int x, y;
223             while (true) {
224                 x = uniform(0, n);
225                 y = uniform(0, n);
226                 if (x == y) continue;
227                 break;
228             }
229             int c = uniform(0, 100);
230             int d = uniform(0, 100);
231             addEdge(g, x, y, c, d);
232             elist[x] ~= E(y, c, d, -1);
233         }
234 
235         auto res = minCostFlow!(int, int, 0)(g, s, t, neg);
236         res.manyFlow(10^^9);
237         int sm = (res.dual[t]-res.dual[s]) * res.capFlow;
238         foreach (i, v; elist) {
239             foreach (e; v) {
240                 sm -= (max(0L, (long(res.dual[e.to]) - res.dual[i]) - e.dist) * e.cap).to!long;
241             }
242         }
243         assert(res.flow == sm);
244     }
245     import dkh.stopwatch;
246     auto ti = benchmark!(f!false, f!true)(500);
247     writeln("MinCostFlow int Random500, Neg500: ", ti[].map!(a => a.toMsecs));
248 }
249 
250 unittest {
251     import std.algorithm, std.conv, std.stdio, std.range;
252     import std.random;
253     import std.typecons;
254 
255     struct E {
256         int to, cap;
257         double dist;
258         int rev;
259     }
260     void addEdge(E[][] g, int from, int to, int cap, double dist) {
261         g[from] ~= E(to, cap, dist, g[to].length.to!int);
262         g[to] ~= E(from, 0, -dist, g[from].length.to!int-1);
263     }
264 
265     void f(bool neg)() {
266         int n = uniform(2, 20);
267         int m = uniform(0, 200);
268         int s, t;
269         while (true) {
270             s = uniform(0, n);
271             t = uniform(0, n);
272             if (s != t) break;
273         }
274         auto g = new E[][n];
275         E[][] elist = new E[][n];
276 
277         foreach (i; 0..m) {
278             int x, y;
279             while (true) {
280                 x = uniform(0, n);
281                 y = uniform(0, n);
282                 if (x == y) continue;
283                 break;
284             }
285             int c = uniform(0, 100);
286             double d = uniform(0.0, 100.0);
287             addEdge(g, x, y, c, d);
288             elist[x] ~= E(y, c, d, -1);
289         }
290 
291         auto res = minCostFlow!(int, double, 1e-9)(g, s, t, neg);
292         res.manyFlow(10^^9);
293         double sm = (res.dual[t]-res.dual[s]) * res.capFlow;
294         foreach (i, v; elist) {
295             foreach (e; v) {
296                 sm -= max(0.0, (res.dual[e.to] - res.dual[i]) - e.dist) * e.cap;
297             }
298         }
299         import std.math;
300         assert(abs(res.flow - sm) <= 1e-3);
301     }
302     import dkh.stopwatch;
303     auto ti = benchmark!(f!false, f!true)(500);
304     writeln("MinCostFlow double Random500, Neg500: ", ti[].map!(a => a.toMsecs));
305 }