diff options
author | Franklin Wei <git@fwei.tk> | 2017-04-29 18:21:56 -0400 |
---|---|---|
committer | Franklin Wei <git@fwei.tk> | 2017-04-29 18:24:42 -0400 |
commit | 881746789a489fad85aae8317555f73dbe261556 (patch) | |
tree | cec2946362c4698c8db3c10f3242ef546c2c22dd /apps/plugins/puzzles/src/maxflow.c | |
parent | 03dd4b92be7dcd5c8ab06da3810887060e06abd5 (diff) | |
download | rockbox-881746789a489fad85aae8317555f73dbe261556.tar.gz rockbox-881746789a489fad85aae8317555f73dbe261556.zip |
puzzles: refactor and resync with upstream
This brings puzzles up-to-date with upstream revision
2d333750272c3967cfd5cd3677572cddeaad5932, though certain changes made
by me, including cursor-only Untangle and some compilation fixes
remain. Upstream code has been moved to its separate subdirectory and
future syncs can be done by simply copying over the new sources.
Change-Id: Ia6506ca5f78c3627165ea6791d38db414ace0804
Diffstat (limited to 'apps/plugins/puzzles/src/maxflow.c')
-rw-r--r-- | apps/plugins/puzzles/src/maxflow.c | 461 |
1 files changed, 461 insertions, 0 deletions
diff --git a/apps/plugins/puzzles/src/maxflow.c b/apps/plugins/puzzles/src/maxflow.c new file mode 100644 index 0000000000..028946b9bd --- /dev/null +++ b/apps/plugins/puzzles/src/maxflow.c | |||
@@ -0,0 +1,461 @@ | |||
1 | /* | ||
2 | * Edmonds-Karp algorithm for finding a maximum flow and minimum | ||
3 | * cut in a network. Almost identical to the Ford-Fulkerson | ||
4 | * algorithm, but apparently using breadth-first search to find the | ||
5 | * _shortest_ augmenting path is a good way to guarantee | ||
6 | * termination and ensure the time complexity is not dependent on | ||
7 | * the actual value of the maximum flow. I don't understand why | ||
8 | * that should be, but it's claimed on the Internet that it's been | ||
9 | * proved, and that's good enough for me. I prefer BFS to DFS | ||
10 | * anyway :-) | ||
11 | */ | ||
12 | |||
13 | #include <assert.h> | ||
14 | #include <stdlib.h> | ||
15 | #include <stdio.h> | ||
16 | |||
17 | #include "maxflow.h" | ||
18 | |||
19 | #include "puzzles.h" /* for snewn/sfree */ | ||
20 | |||
21 | int maxflow_with_scratch(void *scratch, int nv, int source, int sink, | ||
22 | int ne, const int *edges, const int *backedges, | ||
23 | const int *capacity, int *flow, int *cut) | ||
24 | { | ||
25 | int *todo = (int *)scratch; | ||
26 | int *prev = todo + nv; | ||
27 | int *firstedge = todo + 2*nv; | ||
28 | int *firstbackedge = todo + 3*nv; | ||
29 | int i, j, head, tail, from, to; | ||
30 | int totalflow; | ||
31 | |||
32 | /* | ||
33 | * Scan the edges array to find the index of the first edge | ||
34 | * from each node. | ||
35 | */ | ||
36 | j = 0; | ||
37 | for (i = 0; i < ne; i++) | ||
38 | while (j <= edges[2*i]) | ||
39 | firstedge[j++] = i; | ||
40 | while (j < nv) | ||
41 | firstedge[j++] = ne; | ||
42 | assert(j == nv); | ||
43 | |||
44 | /* | ||
45 | * Scan the backedges array to find the index of the first edge | ||
46 | * _to_ each node. | ||
47 | */ | ||
48 | j = 0; | ||
49 | for (i = 0; i < ne; i++) | ||
50 | while (j <= edges[2*backedges[i]+1]) | ||
51 | firstbackedge[j++] = i; | ||
52 | while (j < nv) | ||
53 | firstbackedge[j++] = ne; | ||
54 | assert(j == nv); | ||
55 | |||
56 | /* | ||
57 | * Start the flow off at zero on every edge. | ||
58 | */ | ||
59 | for (i = 0; i < ne; i++) | ||
60 | flow[i] = 0; | ||
61 | totalflow = 0; | ||
62 | |||
63 | /* | ||
64 | * Repeatedly look for an augmenting path, and follow it. | ||
65 | */ | ||
66 | while (1) { | ||
67 | |||
68 | /* | ||
69 | * Set up the prev array. | ||
70 | */ | ||
71 | for (i = 0; i < nv; i++) | ||
72 | prev[i] = -1; | ||
73 | |||
74 | /* | ||
75 | * Initialise the to-do list for BFS. | ||
76 | */ | ||
77 | head = tail = 0; | ||
78 | todo[tail++] = source; | ||
79 | |||
80 | /* | ||
81 | * Now do the BFS loop. | ||
82 | */ | ||
83 | while (head < tail && prev[sink] <= 0) { | ||
84 | from = todo[head++]; | ||
85 | |||
86 | /* | ||
87 | * Try all the forward edges out of node `from'. For a | ||
88 | * forward edge to be valid, it must have flow | ||
89 | * currently less than its capacity. | ||
90 | */ | ||
91 | for (i = firstedge[from]; i < ne && edges[2*i] == from; i++) { | ||
92 | to = edges[2*i+1]; | ||
93 | if (to == source || prev[to] >= 0) | ||
94 | continue; | ||
95 | if (capacity[i] >= 0 && flow[i] >= capacity[i]) | ||
96 | continue; | ||
97 | /* | ||
98 | * This is a valid augmenting edge. Visit node `to'. | ||
99 | */ | ||
100 | prev[to] = 2*i; | ||
101 | todo[tail++] = to; | ||
102 | } | ||
103 | |||
104 | /* | ||
105 | * Try all the backward edges into node `from'. For a | ||
106 | * backward edge to be valid, it must have flow | ||
107 | * currently greater than zero. | ||
108 | */ | ||
109 | for (i = firstbackedge[from]; | ||
110 | j = backedges[i], i < ne && edges[2*j+1]==from; i++) { | ||
111 | to = edges[2*j]; | ||
112 | if (to == source || prev[to] >= 0) | ||
113 | continue; | ||
114 | if (flow[j] <= 0) | ||
115 | continue; | ||
116 | /* | ||
117 | * This is a valid augmenting edge. Visit node `to'. | ||
118 | */ | ||
119 | prev[to] = 2*j+1; | ||
120 | todo[tail++] = to; | ||
121 | } | ||
122 | } | ||
123 | |||
124 | /* | ||
125 | * If prev[sink] is non-null, we have found an augmenting | ||
126 | * path. | ||
127 | */ | ||
128 | if (prev[sink] >= 0) { | ||
129 | int max; | ||
130 | |||
131 | /* | ||
132 | * Work backwards along the path figuring out the | ||
133 | * maximum flow we can add. | ||
134 | */ | ||
135 | to = sink; | ||
136 | max = -1; | ||
137 | while (to != source) { | ||
138 | int spare; | ||
139 | |||
140 | /* | ||
141 | * Find the edge we're currently moving along. | ||
142 | */ | ||
143 | i = prev[to]; | ||
144 | from = edges[i]; | ||
145 | assert(from != to); | ||
146 | |||
147 | /* | ||
148 | * Determine the spare capacity of this edge. | ||
149 | */ | ||
150 | if (i & 1) | ||
151 | spare = flow[i / 2]; /* backward edge */ | ||
152 | else if (capacity[i / 2] >= 0) | ||
153 | spare = capacity[i / 2] - flow[i / 2]; /* forward edge */ | ||
154 | else | ||
155 | spare = -1; /* unlimited forward edge */ | ||
156 | |||
157 | assert(spare != 0); | ||
158 | |||
159 | if (max < 0 || (spare >= 0 && spare < max)) | ||
160 | max = spare; | ||
161 | |||
162 | to = from; | ||
163 | } | ||
164 | /* | ||
165 | * Fail an assertion if max is still < 0, i.e. there is | ||
166 | * an entirely unlimited path from source to sink. Also | ||
167 | * max should not _be_ zero, because by construction | ||
168 | * this _should_ be an augmenting path. | ||
169 | */ | ||
170 | assert(max > 0); | ||
171 | |||
172 | /* | ||
173 | * Now work backwards along the path again, this time | ||
174 | * actually adjusting the flow. | ||
175 | */ | ||
176 | to = sink; | ||
177 | while (to != source) { | ||
178 | /* | ||
179 | * Find the edge we're currently moving along. | ||
180 | */ | ||
181 | i = prev[to]; | ||
182 | from = edges[i]; | ||
183 | assert(from != to); | ||
184 | |||
185 | /* | ||
186 | * Adjust the edge. | ||
187 | */ | ||
188 | if (i & 1) | ||
189 | flow[i / 2] -= max; /* backward edge */ | ||
190 | else | ||
191 | flow[i / 2] += max; /* forward edge */ | ||
192 | |||
193 | to = from; | ||
194 | } | ||
195 | |||
196 | /* | ||
197 | * And adjust the overall flow counter. | ||
198 | */ | ||
199 | totalflow += max; | ||
200 | |||
201 | continue; | ||
202 | } | ||
203 | |||
204 | /* | ||
205 | * If we reach here, we have failed to find an augmenting | ||
206 | * path, which means we're done. Output the `cut' array if | ||
207 | * required, and leave. | ||
208 | */ | ||
209 | if (cut) { | ||
210 | for (i = 0; i < nv; i++) { | ||
211 | if (i == source || prev[i] >= 0) | ||
212 | cut[i] = 0; | ||
213 | else | ||
214 | cut[i] = 1; | ||
215 | } | ||
216 | } | ||
217 | return totalflow; | ||
218 | } | ||
219 | } | ||
220 | |||
221 | int maxflow_scratch_size(int nv) | ||
222 | { | ||
223 | return (nv * 4) * sizeof(int); | ||
224 | } | ||
225 | |||
226 | void maxflow_setup_backedges(int ne, const int *edges, int *backedges) | ||
227 | { | ||
228 | int i, n; | ||
229 | |||
230 | for (i = 0; i < ne; i++) | ||
231 | backedges[i] = i; | ||
232 | |||
233 | /* | ||
234 | * We actually can't use the C qsort() function, because we'd | ||
235 | * need to pass `edges' as a context parameter to its | ||
236 | * comparator function. So instead I'm forced to implement my | ||
237 | * own sorting algorithm internally, which is a pest. I'll use | ||
238 | * heapsort, because I like it. | ||
239 | */ | ||
240 | |||
241 | #define LESS(i,j) ( (edges[2*(i)+1] < edges[2*(j)+1]) || \ | ||
242 | (edges[2*(i)+1] == edges[2*(j)+1] && \ | ||
243 | edges[2*(i)] < edges[2*(j)]) ) | ||
244 | #define PARENT(n) ( ((n)-1)/2 ) | ||
245 | #define LCHILD(n) ( 2*(n)+1 ) | ||
246 | #define RCHILD(n) ( 2*(n)+2 ) | ||
247 | #define SWAP(i,j) do { int swaptmp = (i); (i) = (j); (j) = swaptmp; } while (0) | ||
248 | |||
249 | /* | ||
250 | * Phase 1: build the heap. We want the _largest_ element at | ||
251 | * the top. | ||
252 | */ | ||
253 | n = 0; | ||
254 | while (n < ne) { | ||
255 | n++; | ||
256 | |||
257 | /* | ||
258 | * Swap element n with its parent repeatedly to preserve | ||
259 | * the heap property. | ||
260 | */ | ||
261 | i = n-1; | ||
262 | |||
263 | while (i > 0) { | ||
264 | int p = PARENT(i); | ||
265 | |||
266 | if (LESS(backedges[p], backedges[i])) { | ||
267 | SWAP(backedges[p], backedges[i]); | ||
268 | i = p; | ||
269 | } else | ||
270 | break; | ||
271 | } | ||
272 | } | ||
273 | |||
274 | /* | ||
275 | * Phase 2: repeatedly remove the largest element and stick it | ||
276 | * at the top of the array. | ||
277 | */ | ||
278 | while (n > 0) { | ||
279 | /* | ||
280 | * The largest element is at position 0. Put it at the top, | ||
281 | * and swap the arbitrary element from that position into | ||
282 | * position 0. | ||
283 | */ | ||
284 | n--; | ||
285 | SWAP(backedges[0], backedges[n]); | ||
286 | |||
287 | /* | ||
288 | * Now repeatedly move that arbitrary element down the heap | ||
289 | * by swapping it with the more suitable of its children. | ||
290 | */ | ||
291 | i = 0; | ||
292 | while (1) { | ||
293 | int lc, rc; | ||
294 | |||
295 | lc = LCHILD(i); | ||
296 | rc = RCHILD(i); | ||
297 | |||
298 | if (lc >= n) | ||
299 | break; /* we've hit bottom */ | ||
300 | |||
301 | if (rc >= n) { | ||
302 | /* | ||
303 | * Special case: there is only one child to check. | ||
304 | */ | ||
305 | if (LESS(backedges[i], backedges[lc])) | ||
306 | SWAP(backedges[i], backedges[lc]); | ||
307 | |||
308 | /* _Now_ we've hit bottom. */ | ||
309 | break; | ||
310 | } else { | ||
311 | /* | ||
312 | * The common case: there are two children and we | ||
313 | * must check them both. | ||
314 | */ | ||
315 | if (LESS(backedges[i], backedges[lc]) || | ||
316 | LESS(backedges[i], backedges[rc])) { | ||
317 | /* | ||
318 | * Pick the more appropriate child to swap with | ||
319 | * (i.e. the one which would want to be the | ||
320 | * parent if one were above the other - as one | ||
321 | * is about to be). | ||
322 | */ | ||
323 | if (LESS(backedges[lc], backedges[rc])) { | ||
324 | SWAP(backedges[i], backedges[rc]); | ||
325 | i = rc; | ||
326 | } else { | ||
327 | SWAP(backedges[i], backedges[lc]); | ||
328 | i = lc; | ||
329 | } | ||
330 | } else { | ||
331 | /* This element is in the right place; we're done. */ | ||
332 | break; | ||
333 | } | ||
334 | } | ||
335 | } | ||
336 | } | ||
337 | |||
338 | #undef LESS | ||
339 | #undef PARENT | ||
340 | #undef LCHILD | ||
341 | #undef RCHILD | ||
342 | #undef SWAP | ||
343 | |||
344 | } | ||
345 | |||
346 | int maxflow(int nv, int source, int sink, | ||
347 | int ne, const int *edges, const int *capacity, | ||
348 | int *flow, int *cut) | ||
349 | { | ||
350 | void *scratch; | ||
351 | int *backedges; | ||
352 | int size; | ||
353 | int ret; | ||
354 | |||
355 | /* | ||
356 | * Allocate the space. | ||
357 | */ | ||
358 | size = ne * sizeof(int) + maxflow_scratch_size(nv); | ||
359 | backedges = smalloc(size); | ||
360 | if (!backedges) | ||
361 | return -1; | ||
362 | scratch = backedges + ne; | ||
363 | |||
364 | /* | ||
365 | * Set up the backedges array. | ||
366 | */ | ||
367 | maxflow_setup_backedges(ne, edges, backedges); | ||
368 | |||
369 | /* | ||
370 | * Call the main function. | ||
371 | */ | ||
372 | ret = maxflow_with_scratch(scratch, nv, source, sink, ne, edges, | ||
373 | backedges, capacity, flow, cut); | ||
374 | |||
375 | /* | ||
376 | * Free the scratch space. | ||
377 | */ | ||
378 | sfree(backedges); | ||
379 | |||
380 | /* | ||
381 | * And we're done. | ||
382 | */ | ||
383 | return ret; | ||
384 | } | ||
385 | |||
386 | #ifdef TESTMODE | ||
387 | |||
388 | #define MAXEDGES 256 | ||
389 | #define MAXVERTICES 128 | ||
390 | #define ADDEDGE(i,j) do{edges[ne*2] = (i); edges[ne*2+1] = (j); ne++;}while(0) | ||
391 | |||
392 | int compare_edge(const void *av, const void *bv) | ||
393 | { | ||
394 | const int *a = (const int *)av; | ||
395 | const int *b = (const int *)bv; | ||
396 | |||
397 | if (a[0] < b[0]) | ||
398 | return -1; | ||
399 | else if (a[0] > b[0]) | ||
400 | return +1; | ||
401 | else if (a[1] < b[1]) | ||
402 | return -1; | ||
403 | else if (a[1] > b[1]) | ||
404 | return +1; | ||
405 | else | ||
406 | return 0; | ||
407 | } | ||
408 | |||
409 | int main(void) | ||
410 | { | ||
411 | int edges[MAXEDGES*2], ne, nv; | ||
412 | int capacity[MAXEDGES], flow[MAXEDGES], cut[MAXVERTICES]; | ||
413 | int source, sink, p, q, i, j, ret; | ||
414 | |||
415 | /* | ||
416 | * Use this algorithm to find a maximal complete matching in a | ||
417 | * bipartite graph. | ||
418 | */ | ||
419 | ne = 0; | ||
420 | nv = 0; | ||
421 | source = nv++; | ||
422 | p = nv; | ||
423 | nv += 5; | ||
424 | q = nv; | ||
425 | nv += 5; | ||
426 | sink = nv++; | ||
427 | for (i = 0; i < 5; i++) { | ||
428 | capacity[ne] = 1; | ||
429 | ADDEDGE(source, p+i); | ||
430 | } | ||
431 | for (i = 0; i < 5; i++) { | ||
432 | capacity[ne] = 1; | ||
433 | ADDEDGE(q+i, sink); | ||
434 | } | ||
435 | j = ne; | ||
436 | capacity[ne] = 1; ADDEDGE(p+0,q+0); | ||
437 | capacity[ne] = 1; ADDEDGE(p+1,q+0); | ||
438 | capacity[ne] = 1; ADDEDGE(p+1,q+1); | ||
439 | capacity[ne] = 1; ADDEDGE(p+2,q+1); | ||
440 | capacity[ne] = 1; ADDEDGE(p+2,q+2); | ||
441 | capacity[ne] = 1; ADDEDGE(p+3,q+2); | ||
442 | capacity[ne] = 1; ADDEDGE(p+3,q+3); | ||
443 | capacity[ne] = 1; ADDEDGE(p+4,q+3); | ||
444 | /* capacity[ne] = 1; ADDEDGE(p+2,q+4); */ | ||
445 | qsort(edges, ne, 2*sizeof(int), compare_edge); | ||
446 | |||
447 | ret = maxflow(nv, source, sink, ne, edges, capacity, flow, cut); | ||
448 | |||
449 | printf("ret = %d\n", ret); | ||
450 | |||
451 | for (i = 0; i < ne; i++) | ||
452 | printf("flow %d: %d -> %d\n", flow[i], edges[2*i], edges[2*i+1]); | ||
453 | |||
454 | for (i = 0; i < nv; i++) | ||
455 | if (cut[i] == 0) | ||
456 | printf("difficult set includes %d\n", i); | ||
457 | |||
458 | return 0; | ||
459 | } | ||
460 | |||
461 | #endif | ||