HCoop / hcoop / debian / mlton.git / blame
| Commit | Line | Data |
|---|---|---|
| 7f918cf1 CE |
1 | (* |
| 2 | * Translated from Jeff Siskind's Scheme code by Stephen Weeks | |
| 3 | * (sweeks@sweeks.com). | |
| 4 | * Here is the description from Jeff Siskind (qobi@research.nj.nec.com) | |
| 5 | * | |
| 6 | * It is an implementation of Ratio | |
| 7 | * Regions, an image segmentation/contour finding technique due to Ingemar Cox, | |
| 8 | * Satish Rao, and Yu Zhong. The algorithm is a reduction to max flow, an | |
| 9 | * unpublished technique that Satish described to me. Peter Blicher originally | |
| 10 | * implemented this via a translation to Andrew Goldberg's generic max-flow code. | |
| 11 | * I've reimplemented it, specializing the max-flow algorithm to the particular | |
| 12 | * graphs that are produced by Satish's reduction instead of using Andrew's code. | |
| 13 | * The max-flow algorithm is preflow-push with periodic relabeling and a | |
| 14 | * wave-based heuristic for scheduling pushes and lifts due to Sebastien Roy. | |
| 15 | *) | |
| 16 | ||
| 17 | fun print _ = () | |
| 18 | ||
| 19 | local | |
| 20 | ||
| 21 | fun doo(max: int, f: int -> unit): unit = | |
| 22 | let fun loop i = if i >= max then () else (f i; loop(i + 1)) | |
| 23 | in loop 0 | |
| 24 | end | |
| 25 | ||
| 26 | fun zero x = x = 0 | |
| 27 | val cons = op :: | |
| 28 | val make_vector = Array.array | |
| 29 | val vector_length = Array.length | |
| 30 | val vector_ref = Array.sub | |
| 31 | val vector_set = Array.update | |
| 32 | val map_n_vector = Array.tabulate | |
| 33 | val string_length = String.size | |
| 34 | val string_ref = String.sub | |
| 35 | fun write_char c = () (* TextIO.output1(TextIO.stdOut, c) *) | |
| 36 | val modulo = Int.mod | |
| 37 | val quotient = Int.quot | |
| 38 | fun for_each(l, f) = List.app f l | |
| 39 | fun negative x = x < 0 | |
| 40 | fun positive x = x > 0 | |
| 41 | ||
| 42 | fun min l = | |
| 43 | case l of | |
| 44 | x :: l => let | |
| 45 | fun loop(l, min) = | |
| 46 | case l of | |
| 47 | [] => min | |
| 48 | | x :: l => loop(l, Int.min(min, x)) | |
| 49 | in loop(l, x) | |
| 50 | end | |
| 51 | | _ => raise Fail "min" | |
| 52 | ||
| 53 | fun every_n(n, p) = | |
| 54 | let fun loop i = i >= n orelse (p i andalso loop(i + 1)) | |
| 55 | in loop 0 | |
| 56 | end | |
| 57 | ||
| 58 | fun some(l, p) = List.exists p l | |
| 59 | ||
| 60 | fun some_n(n, p) = | |
| 61 | let fun loop i = i < n andalso (p i orelse loop(i + 1)) | |
| 62 | in loop 0 | |
| 63 | end | |
| 64 | ||
| 65 | fun some_vector(v, p) = | |
| 66 | let | |
| 67 | fun loop i = | |
| 68 | i < vector_length v | |
| 69 | andalso (p(vector_ref(v, i)) | |
| 70 | orelse loop(i + 1)) | |
| 71 | in loop 0 | |
| 72 | end | |
| 73 | ||
| 74 | fun x(x, _) = x | |
| 75 | fun y(_, y) = y | |
| 76 | ||
| 77 | datatype 'a matrix = Matrix of 'a array array | |
| 78 | ||
| 79 | fun make_matrix(m: int, n: int, a: 'a): 'a matrix = | |
| 80 | Matrix(map_n_vector(m, fn i => make_vector(n, a))) | |
| 81 | ||
| 82 | fun matrix_rows(Matrix a) = vector_length a | |
| 83 | fun matrix_columns(Matrix a) = vector_length(vector_ref(a, 0)) | |
| 84 | fun matrix_ref(Matrix a, i, j) = vector_ref(vector_ref(a, i), j) | |
| 85 | fun matrix_set(Matrix a, i, j, x) = vector_set(vector_ref(a, i), j, x) | |
| 86 | ||
| 87 | datatype pormatValue = | |
| 88 | Int of int | |
| 89 | | String of string | |
| 90 | ||
| 91 | fun pormat(control_string: string, values: pormatValue list): unit = | |
| 92 | let | |
| 93 | fun loop(i: int, values: pormatValue list): unit = | |
| 94 | if not(i = string_length control_string) | |
| 95 | then | |
| 96 | let val c = string_ref(control_string, i) | |
| 97 | in if c = #"~" | |
| 98 | then let val c2 = string_ref(control_string, i + 1) | |
| 99 | in case (c2, values) of | |
| 100 | (#"s", Int n :: values) => | |
| 101 | (print(Int.toString n) ; loop(i + 2, values)) | |
| 102 | | (#"a", String s :: values) => | |
| 103 | (print s ; loop(i + 2, values)) | |
| 104 | | (#"%", _) => | |
| 105 | (print "\n"; loop(i + 2, values)) | |
| 106 | | _ => (write_char c; loop(i + 1, values)) | |
| 107 | end | |
| 108 | else (write_char c ; loop(i + 1, values)) | |
| 109 | end | |
| 110 | else () | |
| 111 | in loop(0, values) | |
| 112 | end | |
| 113 | ||
| 114 | (* The vertices are s, t, and (y,x). | |
| 115 | * C_RIGHT[y,x] is the capacity from (y,x) to (y,x+1) which is the same as the | |
| 116 | * capacity from (y,x+1) to (y,x). | |
| 117 | * C_DOWN[y,x] is the capacity from (y,x) to (y+1,x) which is the same as the | |
| 118 | * capacity from (y+1,x) to (y,x). | |
| 119 | * The capacity from s to (y,0), (0,x), (y,Y_1), (0,X_1) is implicitly | |
| 120 | * infinite. | |
| 121 | * The capacity from (x,y) to t is V*W[y,x]. | |
| 122 | * F_RIGHT[y,x] is the preflow from (y,x) to (y,x+1) which is the negation of | |
| 123 | * the preflow from (y,x+1) to (y,x). | |
| 124 | * F_DOWN[y,x] is the preflow from (y,x) to (y+1,x) which is the negation of | |
| 125 | * the preflow from (y+1,x) to (y,x). | |
| 126 | * We do not record the preflow from s to (y,X_1), (y,0), (Y_1,x), and (0,x) | |
| 127 | * and from (y,X_1), (y,0), (Y_1,x), and (0,x) to s. | |
| 128 | * F_T[y,x] is the preflow from (y,x) to t. | |
| 129 | * We do not record the preflow from t to (y,x). | |
| 130 | * {C,F}_RIGHT[0:Y_1,0:X_2]. | |
| 131 | * {C,F}_DOWN[0:Y_2,0:X_1]. | |
| 132 | * F_T[0:Y_1,0:X_1] | |
| 133 | * For now, we will keep all capacities (and thus all preflows) as integers. | |
| 134 | * (CF_RIGHT y x) is the residual capacity from (y,x) to (y,x+1). | |
| 135 | * (CF_LEFT y x) is the residual capacity from (y,x) to (y,x_1). | |
| 136 | * (CF_DOWN y x) is the residual capacity from (y,x) to (y+1,x). | |
| 137 | * (CF_UP y x) is the residual capacity from (y,x) to (y_1,x). | |
| 138 | * We do not compute the residual capacities from s to (y,X_1), (y,0), | |
| 139 | * (Y_1,x), and (0,x) because they are all infinite. | |
| 140 | * We do not compute the residual capacities from (y,X_1), (y,0), (Y_1,x), | |
| 141 | * and (0,x) to s because they will never be used. | |
| 142 | * (CF_T y x) is the residual capacity from (y,x) to t. | |
| 143 | * We do not compute the residual capacity from t to (y,x) because it will | |
| 144 | * be used. | |
| 145 | * (EF_RIGHT? y x) is true if there is an edge from (y,x) to (y,x+1) in the | |
| 146 | * residual network. | |
| 147 | * (EF_LEFT? y x) is true if there is an edge from (y,x) to (y,x_1) in the | |
| 148 | * residual network. | |
| 149 | * (EF_DOWN? y x) is true if there is an edge from (y,x) to (y+1,x) in the | |
| 150 | * residual network. | |
| 151 | * (EF_UP? y x) is true if there is an edge from (y,x) to (y_1,x) in the | |
| 152 | * residual network. | |
| 153 | * (EF_T? y x) is true if there is an edge from (y,x) to t in the | |
| 154 | * residual network. | |
| 155 | * There are always edges in the residual network from s to (y,X_1), (y,0), | |
| 156 | * (Y_1,x), and (0,x). | |
| 157 | * We don't care whether there are edges in the residual network from | |
| 158 | * (y,X_1), (y,0), (Y_1,x), and (0,x) to s because they will never be used. | |
| 159 | * We don't care whether there are edges in the residual network from t to | |
| 160 | * (y,x) because they will never be used. | |
| 161 | *) | |
| 162 | ||
| 163 | fun positive_min(x, y) = if negative x then y else Int.min(x, y) | |
| 164 | ||
| 165 | fun positive_minus(x, y) = if negative x then x else x - y | |
| 166 | ||
| 167 | fun positive_plus(x, y) = if negative x then x else x + y | |
| 168 | ||
| 169 | fun rao_ratio_region(c_right, c_down, w, lg_max_v) = | |
| 170 | let val height = matrix_rows w | |
| 171 | val width = matrix_columns w | |
| 172 | val f_right = make_matrix(height, width - 1, 0) | |
| 173 | val f_down = make_matrix(height - 1, width, 0) | |
| 174 | val f_t = make_matrix(height, width, 0) | |
| 175 | val h = make_matrix(height, width, 0) | |
| 176 | val e = make_matrix(height, width, 0) | |
| 177 | val marked = make_matrix(height, width, false) | |
| 178 | val m1 = height * width + 2 | |
| 179 | val m2 = 2 * height * width + 2 | |
| 180 | val q = make_vector(2 * height * width + 3, []) | |
| 181 | fun cf_right(y, x) = | |
| 182 | matrix_ref(c_right, y, x) - matrix_ref(f_right, y, x) | |
| 183 | fun cf_left(y, x) = | |
| 184 | matrix_ref(c_right, y, x - 1) + matrix_ref(f_right, y, x - 1) | |
| 185 | fun cf_down(y, x) = | |
| 186 | matrix_ref(c_down, y, x) - matrix_ref(f_down, y, x) | |
| 187 | fun cf_up(y, x) = | |
| 188 | matrix_ref(c_down, y - 1, x) + matrix_ref(f_down, y - 1, x) | |
| 189 | fun ef_right(y, x) = positive(cf_right(y, x)) | |
| 190 | fun ef_left(y, x) = positive(cf_left(y, x)) | |
| 191 | fun ef_down(y, x) = positive(cf_down(y, x)) | |
| 192 | fun ef_up(y, x) = positive(cf_up(y, x)) | |
| 193 | fun preflow_push v = | |
| 194 | let | |
| 195 | fun enqueue(y, x) = | |
| 196 | if not(matrix_ref(marked, y, x)) | |
| 197 | then | |
| 198 | (vector_set(q, | |
| 199 | matrix_ref(h, y, x), | |
| 200 | (cons((x, y), | |
| 201 | vector_ref(q, matrix_ref(h, y, x))))) | |
| 202 | ; matrix_set(marked, y, x, true)) | |
| 203 | else () | |
| 204 | fun cf_t(y, x) = v * matrix_ref(w, y, x) - matrix_ref(f_t, y, x) | |
| 205 | fun ef_t(y, x) = positive(cf_t(y, x)) | |
| 206 | fun can_push_right(y, x) = | |
| 207 | x < width - 1 | |
| 208 | andalso not(zero(matrix_ref(e, y, x))) | |
| 209 | andalso ef_right(y, x) | |
| 210 | andalso matrix_ref(h, y, x) = matrix_ref(h, y, x + 1) + 1 | |
| 211 | fun can_push_left(y, x) = | |
| 212 | x > 0 | |
| 213 | andalso not(zero(matrix_ref(e, y, x))) | |
| 214 | andalso ef_left(y, x) | |
| 215 | andalso matrix_ref(h, y, x) = matrix_ref(h, y, x - 1) + 1 | |
| 216 | fun can_push_down(y, x) = | |
| 217 | y < height - 1 | |
| 218 | andalso not(zero(matrix_ref(e, y, x))) | |
| 219 | andalso ef_down(y, x) | |
| 220 | andalso matrix_ref(h, y, x) = matrix_ref(h, y + 1, x) + 1 | |
| 221 | fun can_push_up(y, x) = | |
| 222 | y > 0 | |
| 223 | andalso not(zero(matrix_ref(e, y, x))) | |
| 224 | andalso ef_up(y, x) | |
| 225 | andalso matrix_ref(h, y, x) = matrix_ref(h, y - 1, x) + 1 | |
| 226 | fun can_push_t(y, x) = | |
| 227 | not(zero(matrix_ref(e, y, x))) | |
| 228 | andalso ef_t(y, x) | |
| 229 | andalso matrix_ref(h, y, x) = 1 | |
| 230 | fun can_lift(y, x) = | |
| 231 | not(zero(matrix_ref(e, y, x))) | |
| 232 | andalso (if x = width - 1 | |
| 233 | then matrix_ref(h, y, x) <= m1 | |
| 234 | else (not(ef_right(y, x)) | |
| 235 | orelse | |
| 236 | matrix_ref(h, y, x) <= matrix_ref(h, y, x + 1))) | |
| 237 | andalso (if x = 0 | |
| 238 | then matrix_ref(h, y, x) <= m1 | |
| 239 | else (not(ef_left(y, x)) | |
| 240 | orelse | |
| 241 | matrix_ref(h, y, x) <= matrix_ref(h, y, x - 1))) | |
| 242 | andalso (if y = height - 1 | |
| 243 | then matrix_ref(h, y, x) <= m1 | |
| 244 | else (not(ef_down(y, x)) | |
| 245 | orelse | |
| 246 | matrix_ref(h, y, x) <= matrix_ref(h, y + 1, x))) | |
| 247 | andalso (if y = 0 | |
| 248 | then matrix_ref(h, y, x) <= m1 | |
| 249 | else (not(ef_up(y, x)) | |
| 250 | orelse | |
| 251 | matrix_ref(h, y, x) <= matrix_ref(h, y - 1, x))) | |
| 252 | andalso (not(ef_t(y, x)) orelse matrix_ref(h, y, x) = 0) | |
| 253 | fun push_right(y, x) = | |
| 254 | (* (pormat "Push right ~s ~s~%" y x) *) | |
| 255 | let val df_u_v = positive_min(matrix_ref(e, y, x), cf_right(y, x)) | |
| 256 | in matrix_set(f_right, y, x, matrix_ref(f_right, y, x) + df_u_v) | |
| 257 | ; matrix_set(e, y, x, | |
| 258 | positive_minus(matrix_ref(e, y, x), df_u_v)) | |
| 259 | ; matrix_set(e, y, x + 1, | |
| 260 | positive_plus(matrix_ref(e, y, x + 1), df_u_v)) | |
| 261 | ; enqueue(y, x + 1) | |
| 262 | end | |
| 263 | fun push_left(y, x) = | |
| 264 | (* (pormat "Push left ~s ~s~%" y x) *) | |
| 265 | let val df_u_v = positive_min(matrix_ref(e, y, x), cf_left(y, x)) | |
| 266 | in matrix_set(f_right, y, x - 1, | |
| 267 | matrix_ref(f_right, y, x - 1) - df_u_v) | |
| 268 | ; matrix_set(e, y, x, | |
| 269 | positive_minus(matrix_ref(e, y, x), df_u_v)) | |
| 270 | ; matrix_set(e, y, x - 1, | |
| 271 | positive_plus(matrix_ref(e, y, x - 1), df_u_v)) | |
| 272 | ; enqueue(y, x - 1) | |
| 273 | end | |
| 274 | ||
| 275 | fun push_down(y, x) = | |
| 276 | (* (pormat "Push down ~s ~s~%" y x) *) | |
| 277 | let val df_u_v = positive_min(matrix_ref(e, y, x), cf_down(y, x)) | |
| 278 | in matrix_set(f_down, y, x, matrix_ref(f_down, y, x) + df_u_v) | |
| 279 | ; matrix_set(e, y, x, | |
| 280 | positive_minus(matrix_ref(e, y, x), df_u_v)) | |
| 281 | ; matrix_set(e, y + 1, x, | |
| 282 | positive_plus(matrix_ref(e, y + 1, x), df_u_v)) | |
| 283 | ; enqueue(y + 1, x) | |
| 284 | end | |
| 285 | fun push_up(y, x) = | |
| 286 | (* ;;(pormat "Push up ~s ~s~%" y x) *) | |
| 287 | let val df_u_v = positive_min(matrix_ref(e, y, x), cf_up(y, x)) | |
| 288 | in matrix_set(f_down, y - 1, x, | |
| 289 | matrix_ref(f_down, y - 1, x) - df_u_v) | |
| 290 | ; matrix_set(e, y, x, | |
| 291 | positive_minus(matrix_ref(e, y, x), df_u_v)) | |
| 292 | ; matrix_set(e, y - 1, x, | |
| 293 | positive_plus(matrix_ref(e, y - 1, x), df_u_v)) | |
| 294 | ; enqueue(y - 1, x) | |
| 295 | end | |
| 296 | fun push_t(y, x) = | |
| 297 | (* ;;(pormat "Push t ~s ~s~%" y x) *) | |
| 298 | let val df_u_v = positive_min(matrix_ref(e, y, x), cf_t(y, x)) | |
| 299 | in matrix_set(f_t, y, x, matrix_ref(f_t, y, x) + df_u_v) | |
| 300 | ; matrix_set(e, y, x, | |
| 301 | positive_minus(matrix_ref(e, y, x), df_u_v)) | |
| 302 | end | |
| 303 | fun lift(y, x) = | |
| 304 | (* ;;(pormat "Lift ~s ~s~%" y x) *) | |
| 305 | matrix_set | |
| 306 | (h, y, x, | |
| 307 | 1 + min[if x = width - 1 | |
| 308 | then m1 | |
| 309 | else if ef_right(y, x) | |
| 310 | then matrix_ref(h, y, x + 1) | |
| 311 | else m2, | |
| 312 | if x = 0 | |
| 313 | then m1 | |
| 314 | else if ef_left(y, x) | |
| 315 | then matrix_ref(h, y, x - 1) | |
| 316 | else m2, | |
| 317 | if y = height - 1 | |
| 318 | then m1 | |
| 319 | else if ef_down(y, x) | |
| 320 | then matrix_ref(h, y + 1, x) | |
| 321 | else m2, | |
| 322 | if y = 0 | |
| 323 | then m1 | |
| 324 | else if ef_up(y, x) | |
| 325 | then matrix_ref(h, y - 1, x) | |
| 326 | else m2, | |
| 327 | if ef_t(y, x) then 0 else m2]) | |
| 328 | fun relabel() = | |
| 329 | (* ;;(pormat "Relabel~%") *) | |
| 330 | let | |
| 331 | datatype 'a queue = | |
| 332 | Nil | |
| 333 | | Cons of 'a * 'a queue ref | |
| 334 | fun null(q: 'q queue ref) = | |
| 335 | case !q of | |
| 336 | Nil => true | |
| 337 | | _ => false | |
| 338 | val q: (int * int) queue ref = ref Nil | |
| 339 | val tail: (int * int) queue ref = ref Nil | |
| 340 | fun enqueue(y, x, value) = | |
| 341 | if value < matrix_ref(h, y, x) | |
| 342 | then (matrix_set(h, y, x, value) | |
| 343 | ; if not(matrix_ref(marked, y, x)) | |
| 344 | then (matrix_set(marked, y, x, true) | |
| 345 | ; (case !tail of | |
| 346 | Nil => | |
| 347 | (tail := Cons((x, y), ref Nil) | |
| 348 | ; q := !tail) | |
| 349 | | Cons(_, cdr) => | |
| 350 | (cdr := Cons((x, y), ref Nil) | |
| 351 | ; tail := !cdr))) | |
| 352 | else ()) | |
| 353 | else () | |
| 354 | fun dequeue() = | |
| 355 | case !q of | |
| 356 | Nil => raise Fail "dequeue" | |
| 357 | | Cons(p, rest) => | |
| 358 | (matrix_set(marked, y p, x p, false) | |
| 359 | ; q := !rest | |
| 360 | ; if null q then tail := Nil else () | |
| 361 | ; p) | |
| 362 | in doo(height, fn y => | |
| 363 | doo(width, fn x => | |
| 364 | (matrix_set(h, y, x, m1) | |
| 365 | ; matrix_set(marked, y, x, false)))) | |
| 366 | ; doo(height, fn y => | |
| 367 | doo(width, fn x => | |
| 368 | if ef_t(y, x) | |
| 369 | andalso matrix_ref(h, y, x) > 1 | |
| 370 | then enqueue(y, x, 1) | |
| 371 | else ())) | |
| 372 | ; let | |
| 373 | fun loop() = | |
| 374 | if not(null q) | |
| 375 | then | |
| 376 | (let val p = dequeue() | |
| 377 | val x = x p | |
| 378 | val y = y p | |
| 379 | val value = matrix_ref(h, y, x) + 1 | |
| 380 | in if x > 0 andalso ef_right(y, x - 1) | |
| 381 | then enqueue(y, x - 1, value) | |
| 382 | else () | |
| 383 | ; if x < width - 1 andalso ef_left(y, x + 1) | |
| 384 | then enqueue(y, x + 1, value) | |
| 385 | else () | |
| 386 | ; if y > 0 andalso ef_down(y - 1, x) | |
| 387 | then enqueue(y - 1, x, value) | |
| 388 | else () | |
| 389 | ; if y < height - 1 andalso ef_up(y + 1, x) | |
| 390 | then enqueue(y + 1, x, value) | |
| 391 | else () | |
| 392 | end | |
| 393 | ; loop()) | |
| 394 | else () | |
| 395 | in loop() | |
| 396 | end | |
| 397 | end (* relabel *) | |
| 398 | in doo(height, fn y => | |
| 399 | doo(width, fn x => | |
| 400 | (matrix_set(e, y, x, 0) | |
| 401 | ; matrix_set(f_t, y, x, 0)))) | |
| 402 | ; doo(height, fn y => | |
| 403 | doo(width - 1, fn x => | |
| 404 | matrix_set(f_right, y, x, 0))) | |
| 405 | ; doo(height - 1, fn y => | |
| 406 | doo(width, fn x => | |
| 407 | matrix_set(f_down, y, x, 0))) | |
| 408 | ; doo(height, fn y => | |
| 409 | (matrix_set(e, y, width - 1, ~1) | |
| 410 | ; matrix_set(e, y, 0, ~1))) | |
| 411 | ; doo(width - 1, fn x => | |
| 412 | (matrix_set(e, height - 1, x, ~1) | |
| 413 | ; matrix_set(e, 0, x, ~1))) | |
| 414 | ; let val pushes = ref 0 | |
| 415 | val lifts = ref 0 | |
| 416 | val relabels = ref 0 | |
| 417 | fun loop(i, p) = | |
| 418 | if zero(modulo(i, 6)) andalso not p | |
| 419 | then (relabel() | |
| 420 | ; relabels := !relabels + 1 | |
| 421 | ; if every_n(height, fn y => | |
| 422 | every_n(width, fn x => | |
| 423 | zero(matrix_ref(e, y, x)) | |
| 424 | orelse | |
| 425 | matrix_ref(h, y, x) = m1)) | |
| 426 | then | |
| 427 | (* Every vertex with excess capacity is not reachable from the sink in | |
| 428 | * the inverse residual network. So terminate early because we have | |
| 429 | * already found a min cut. In this case, the preflows and excess | |
| 430 | * capacities will not be correct. But the cut is indicated by the | |
| 431 | * heights. Vertices reachable from the source have height | |
| 432 | * HEIGHT * WIDTH + 2 while vertices reachable from the sink have | |
| 433 | * smaller height. Early termination is necessary with relabeling to | |
| 434 | * prevent an infinite loop. The loop arises because vertices that are | |
| 435 | * not reachable from the sink in the inverse residual network have | |
| 436 | * their height reset to HEIGHT * WIDTH + 2 by the relabeling | |
| 437 | * process. If there are such vertices with excess capacity, this is | |
| 438 | * not high enough for the excess capacity to be pushed back to the | |
| 439 | * perimeter. So after relabeling, vertices get lifted to try to push | |
| 440 | * excess capacity back to the perimeter but then a relabeling happens | |
| 441 | * to soon and foils this lifting. Terminating when all vertices with | |
| 442 | * excess capacity are not reachable from the sink in the inverse | |
| 443 | * residual network eliminates this problem. | |
| 444 | *) | |
| 445 | (pormat | |
| 446 | ("~s push~a, ~s lift~a, ~s relabel~a, ~s wave~a, terminated early~%", | |
| 447 | [Int(! pushes), | |
| 448 | String(if !pushes = 1 then "" else "es"), | |
| 449 | Int(! lifts), | |
| 450 | String(if !lifts = 1 then "" else "s"), | |
| 451 | Int(! relabels), | |
| 452 | String(if !relabels = 1 then "" else "s"), | |
| 453 | Int i, | |
| 454 | String(if i = 1 then "" else "s")])) | |
| 455 | else | |
| 456 | (* We need to rebuild the priority queue after relabeling since the | |
| 457 | * heights might have changed and the priority queue is indexed by | |
| 458 | * height. This also assumes that a relabel is done before any pushes | |
| 459 | * or lifts. | |
| 460 | *) | |
| 461 | (doo(vector_length q, fn k => | |
| 462 | vector_set(q, k, [])) | |
| 463 | ; doo(height, fn y => | |
| 464 | doo(width, fn x => | |
| 465 | matrix_set(marked, y, x, false))) | |
| 466 | ; doo(height, fn y => | |
| 467 | doo(width, fn x => | |
| 468 | if not(zero(matrix_ref(e, y, x))) | |
| 469 | then enqueue(y, x) | |
| 470 | else ())) | |
| 471 | ; loop(i, true))) | |
| 472 | else if some_vector(q, fn ps => | |
| 473 | some(ps, fn p => | |
| 474 | let val x = x p | |
| 475 | val y = y p | |
| 476 | in can_push_right(y, x) | |
| 477 | orelse can_push_left(y, x) | |
| 478 | orelse can_push_down(y, x) | |
| 479 | orelse can_push_up(y, x) | |
| 480 | orelse can_push_t(y, x) | |
| 481 | orelse can_lift(y, x) | |
| 482 | end)) | |
| 483 | then | |
| 484 | ( | |
| 485 | let fun loop k = | |
| 486 | if not(negative k) | |
| 487 | then | |
| 488 | ( | |
| 489 | let val ps = vector_ref(q, k) | |
| 490 | in vector_set(q, k, []) | |
| 491 | ; (for_each | |
| 492 | (ps, fn p => | |
| 493 | matrix_set(marked, y p, x p, | |
| 494 | false))) | |
| 495 | ; (for_each | |
| 496 | (ps, fn p => | |
| 497 | let val x = x p | |
| 498 | val y = y p | |
| 499 | in if can_push_right(y, x) | |
| 500 | then (pushes := !pushes + 1 | |
| 501 | ; push_right(y, x)) | |
| 502 | else () | |
| 503 | ; if can_push_left(y, x) | |
| 504 | then (pushes := !pushes + 1 | |
| 505 | ; push_left(y, x)) | |
| 506 | else () | |
| 507 | ; if can_push_down(y, x) | |
| 508 | then (pushes := !pushes + 1 | |
| 509 | ; push_down(y, x)) | |
| 510 | else () | |
| 511 | ; if can_push_up(y, x) | |
| 512 | then (pushes := !pushes + 1 | |
| 513 | ; push_up(y, x)) | |
| 514 | else () | |
| 515 | ; if can_push_t(y, x) | |
| 516 | then (pushes := !pushes + 1 | |
| 517 | ; push_t(y, x)) | |
| 518 | else () | |
| 519 | ; if can_lift(y, x) | |
| 520 | then (lifts := !lifts + 1 | |
| 521 | ; lift(y, x)) | |
| 522 | else () | |
| 523 | ; if not(zero(matrix_ref(e, y, x))) | |
| 524 | then enqueue(y, x) | |
| 525 | else () | |
| 526 | end)) | |
| 527 | end | |
| 528 | ; loop(k - 1)) | |
| 529 | else () | |
| 530 | in loop(vector_length q - 1) | |
| 531 | end | |
| 532 | ; loop(i + 1, false)) | |
| 533 | else | |
| 534 | (* This is so MIN_CUT and MIN_CUT_INCLUDES_EVERY_EDGE_TO_T work. *) | |
| 535 | (relabel() | |
| 536 | ; relabels := !relabels + 1 | |
| 537 | ; (pormat("~s push~a, ~s lift~a, ~s relabel~a, ~s wave~a~%", | |
| 538 | [Int(! pushes), | |
| 539 | String(if !pushes = 1 then "" else "es"), | |
| 540 | Int(! lifts), | |
| 541 | String(if !lifts = 1 then "" else "s"), | |
| 542 | Int(! relabels), | |
| 543 | String(if !relabels = 1 then "" else "s"), | |
| 544 | Int i, | |
| 545 | String(if i = 1 then "" else "s")]))) | |
| 546 | in loop(0, false) | |
| 547 | end | |
| 548 | end | |
| 549 | fun min_cut_includes_every_edge_to_t() = | |
| 550 | (* This requires that a relabel was done immediately before returning from | |
| 551 | * PREFLOW_PUSH. | |
| 552 | *) | |
| 553 | every_n(height, fn y => | |
| 554 | every_n(width, fn x => | |
| 555 | matrix_ref(h, y, x) = m1)) | |
| 556 | fun min_cut() = | |
| 557 | (* This requires that a relabel was done immediately before returning from | |
| 558 | * PREFLOW_PUSH | |
| 559 | *) | |
| 560 | map_n_vector | |
| 561 | (height, fn y => | |
| 562 | map_n_vector(width, fn x => | |
| 563 | not(matrix_ref(h, y, x) = m1))) | |
| 564 | fun loop(lg_v, v_max) = | |
| 565 | if negative lg_v | |
| 566 | then (pormat("V-MAX=~s~%",[Int v_max]) | |
| 567 | ; preflow_push(v_max + 1) | |
| 568 | ; min_cut()) | |
| 569 | else let val v = v_max + let | |
| 570 | fun loop(i, c) = | |
| 571 | if (zero i) | |
| 572 | then c | |
| 573 | else loop(i - 1, c + c) | |
| 574 | in loop(lg_v, 1) | |
| 575 | end | |
| 576 | in pormat("LG-V=~s, V-MAX=~s, V=~s~%", | |
| 577 | [Int lg_v, Int v_max, Int v]) | |
| 578 | ; preflow_push v | |
| 579 | ; loop(lg_v - 1, | |
| 580 | if min_cut_includes_every_edge_to_t() | |
| 581 | then v | |
| 582 | else v_max) | |
| 583 | end | |
| 584 | in loop(lg_max_v, 0) | |
| 585 | end | |
| 586 | ||
| 587 | in | |
| 588 | ||
| 589 | fun doit n = | |
| 590 | let val height = n | |
| 591 | val width = n | |
| 592 | val lg_max_v = 15 | |
| 593 | val c_right = make_matrix(height, width - 1, ~1) | |
| 594 | val c_down = make_matrix(height - 1, width, ~1) | |
| 595 | in doo(height, fn y => | |
| 596 | doo(width - 1, fn x => | |
| 597 | matrix_set | |
| 598 | (c_right, y, x, | |
| 599 | if (y >= quotient(height, 4) | |
| 600 | andalso y < quotient(3 * height, 4) | |
| 601 | andalso (x = quotient(width, 4) - 1 | |
| 602 | orelse x = quotient(3 * width, 4) - 1)) | |
| 603 | then 1 | |
| 604 | else 128))) | |
| 605 | ; doo(height - 1, fn y => | |
| 606 | doo(width, fn x => | |
| 607 | matrix_set | |
| 608 | (c_down, y, x, | |
| 609 | if (x >= quotient(width, 4) | |
| 610 | andalso x < quotient(3 * width, 4) | |
| 611 | andalso (y = quotient(height, 4) - 1 | |
| 612 | orelse y = quotient(3 * height, 4) - 1)) | |
| 613 | then 1 | |
| 614 | else 128))) | |
| 615 | ; rao_ratio_region(c_right, c_down, | |
| 616 | make_matrix(height, width, 1), | |
| 617 | lg_max_v) | |
| 618 | end | |
| 619 | ||
| 620 | end | |
| 621 | ||
| 622 | structure Main = | |
| 623 | struct | |
| 624 | val doit = doit | |
| 625 | end |