maxflow.cpp 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889
  1. /* maxflow.cpp */
  2. /*
  3. Copyright 2001 Vladimir Kolmogorov (vnk@cs.cornell.edu), Yuri Boykov (yuri@csd.uwo.ca).
  4. This program is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU General Public License as published by
  6. the Free Software Foundation; either version 2 of the License, or
  7. (at your option) any later version.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11. GNU General Public License for more details.
  12. You should have received a copy of the GNU General Public License
  13. along with this program; if not, write to the Free Software
  14. Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  15. */
  16. #include <stdio.h>
  17. #include "graph.h"
  18. using namespace OBJREC;
  19. /*
  20. special constants for node->parent
  21. */
  22. #define TERMINAL ( (arc_forward *) 1 ) /* to terminal */
  23. #define ORPHAN ( (arc_forward *) 2 ) /* orphan */
  24. #define INFINITE_D 1000000000 /* infinite distance to the terminal */
  25. /***********************************************************************/
  26. /*
  27. Functions for processing active list.
  28. i->next points to the next node in the list
  29. (or to i, if i is the last node in the list).
  30. If i->next is NULL iff i is not in the list.
  31. There are two queues. Active nodes are added
  32. to the end of the second queue and read from
  33. the front of the first queue. If the first queue
  34. is empty, it is replaced by the second queue
  35. (and the second queue becomes empty).
  36. */
  37. inline void Graph::set_active(node *i)
  38. {
  39. if (!i->next)
  40. {
  41. /* it's not in the list yet */
  42. if (queue_last[1]) queue_last[1] -> next = i;
  43. else queue_first[1] = i;
  44. queue_last[1] = i;
  45. i -> next = i;
  46. }
  47. }
  48. /*
  49. Returns the next active node.
  50. If it is connected to the sink, it stays in the list,
  51. otherwise it is removed from the list
  52. */
  53. inline Graph::node * Graph::next_active()
  54. {
  55. node *i;
  56. while ( 1 )
  57. {
  58. if (!(i = queue_first[0]))
  59. {
  60. queue_first[0] = i = queue_first[1];
  61. queue_last[0] = queue_last[1];
  62. queue_first[1] = NULL;
  63. queue_last[1] = NULL;
  64. if (!i) return NULL;
  65. }
  66. /* remove it from the active list */
  67. if (i->next == i) queue_first[0] = queue_last[0] = NULL;
  68. else queue_first[0] = i -> next;
  69. i -> next = NULL;
  70. /* a node in the list is active iff it has a parent */
  71. if (i->parent) return i;
  72. }
  73. }
  74. /***********************************************************************/
  75. void Graph::maxflow_init()
  76. {
  77. node *i;
  78. node_block *nb;
  79. queue_first[0] = queue_last[0] = NULL;
  80. queue_first[1] = queue_last[1] = NULL;
  81. orphan_first = NULL;
  82. for (nb = node_block_first; nb; nb = nb->next)
  83. for (i = &nb->nodes[0]; i < nb->current; i++)
  84. {
  85. i -> next = NULL;
  86. i -> TS = 0;
  87. if (i->tr_cap > 0)
  88. {
  89. /* i is connected to the source */
  90. i -> is_sink = 0;
  91. i -> parent = TERMINAL;
  92. set_active(i);
  93. i -> TS = 0;
  94. i -> DIST = 1;
  95. }
  96. else if (i->tr_cap < 0)
  97. {
  98. /* i is connected to the sink */
  99. i -> is_sink = 1;
  100. i -> parent = TERMINAL;
  101. set_active(i);
  102. i -> TS = 0;
  103. i -> DIST = 1;
  104. }
  105. else
  106. {
  107. //seems to be buggy
  108. i -> parent = NULL;
  109. i -> is_sink = 1;
  110. }
  111. }
  112. TIME = 0;
  113. }
  114. /***********************************************************************/
  115. void Graph::augment(node *s_start, node *t_start, captype *cap_middle, captype *rev_cap_middle)
  116. {
  117. node *i;
  118. arc_forward *a;
  119. captype bottleneck;
  120. nodeptr *np;
  121. /* 1. Finding bottleneck capacity */
  122. /* 1a - the source tree */
  123. bottleneck = *cap_middle;
  124. for (i = s_start; ; )
  125. {
  126. a = i -> parent;
  127. if (a == TERMINAL) break;
  128. if (IS_ODD(a))
  129. {
  130. a = MAKE_EVEN(a);
  131. if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
  132. i = NEIGHBOR_NODE_REV(i, a -> shift);
  133. }
  134. else
  135. {
  136. if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
  137. i = NEIGHBOR_NODE(i, a -> shift);
  138. }
  139. }
  140. if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
  141. /* 1b - the sink tree */
  142. for (i = t_start; ; )
  143. {
  144. a = i -> parent;
  145. if (a == TERMINAL) break;
  146. if (IS_ODD(a))
  147. {
  148. a = MAKE_EVEN(a);
  149. if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;
  150. i = NEIGHBOR_NODE_REV(i, a -> shift);
  151. }
  152. else
  153. {
  154. if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
  155. i = NEIGHBOR_NODE(i, a -> shift);
  156. }
  157. }
  158. if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
  159. /* 2. Augmenting */
  160. /* 2a - the source tree */
  161. *rev_cap_middle += bottleneck;
  162. *cap_middle -= bottleneck;
  163. for (i = s_start; ; )
  164. {
  165. a = i -> parent;
  166. if (a == TERMINAL) break;
  167. if (IS_ODD(a))
  168. {
  169. a = MAKE_EVEN(a);
  170. a -> r_rev_cap += bottleneck;
  171. a -> r_cap -= bottleneck;
  172. if (!a->r_cap)
  173. {
  174. /* add i to the adoption list */
  175. i -> parent = ORPHAN;
  176. np = nodeptr_block -> New();
  177. np -> ptr = i;
  178. np -> next = orphan_first;
  179. orphan_first = np;
  180. }
  181. i = NEIGHBOR_NODE_REV(i, a -> shift);
  182. }
  183. else
  184. {
  185. a -> r_cap += bottleneck;
  186. a -> r_rev_cap -= bottleneck;
  187. if (!a->r_rev_cap)
  188. {
  189. /* add i to the adoption list */
  190. i -> parent = ORPHAN;
  191. np = nodeptr_block -> New();
  192. np -> ptr = i;
  193. np -> next = orphan_first;
  194. orphan_first = np;
  195. }
  196. i = NEIGHBOR_NODE(i, a -> shift);
  197. }
  198. }
  199. i -> tr_cap -= bottleneck;
  200. if (!i->tr_cap)
  201. {
  202. /* add i to the adoption list */
  203. i -> parent = ORPHAN;
  204. np = nodeptr_block -> New();
  205. np -> ptr = i;
  206. np -> next = orphan_first;
  207. orphan_first = np;
  208. }
  209. /* 2b - the sink tree */
  210. for (i = t_start; ; )
  211. {
  212. a = i -> parent;
  213. if (a == TERMINAL) break;
  214. if (IS_ODD(a))
  215. {
  216. a = MAKE_EVEN(a);
  217. a -> r_cap += bottleneck;
  218. a -> r_rev_cap -= bottleneck;
  219. if (!a->r_rev_cap)
  220. {
  221. /* add i to the adoption list */
  222. i -> parent = ORPHAN;
  223. np = nodeptr_block -> New();
  224. np -> ptr = i;
  225. np -> next = orphan_first;
  226. orphan_first = np;
  227. }
  228. i = NEIGHBOR_NODE_REV(i, a -> shift);
  229. }
  230. else
  231. {
  232. a -> r_rev_cap += bottleneck;
  233. a -> r_cap -= bottleneck;
  234. if (!a->r_cap)
  235. {
  236. /* add i to the adoption list */
  237. i -> parent = ORPHAN;
  238. np = nodeptr_block -> New();
  239. np -> ptr = i;
  240. np -> next = orphan_first;
  241. orphan_first = np;
  242. }
  243. i = NEIGHBOR_NODE(i, a -> shift);
  244. }
  245. }
  246. i -> tr_cap += bottleneck;
  247. if (!i->tr_cap)
  248. {
  249. /* add i to the adoption list */
  250. i -> parent = ORPHAN;
  251. np = nodeptr_block -> New();
  252. np -> ptr = i;
  253. np -> next = orphan_first;
  254. orphan_first = np;
  255. }
  256. flow += bottleneck;
  257. }
  258. /***********************************************************************/
  259. void Graph::process_source_orphan(node *i)
  260. {
  261. node *j;
  262. arc_forward *a0_for, *a0_for_first, *a0_for_last;
  263. arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
  264. arc_forward *a0_min = NULL, *a;
  265. nodeptr *np;
  266. int d, d_min = INFINITE_D;
  267. /* trying to find a new parent */
  268. a0_for_first = i -> first_out;
  269. if (IS_ODD(a0_for_first))
  270. {
  271. a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
  272. a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
  273. }
  274. else a0_for_last = (i + 1) -> first_out;
  275. a0_rev_first = i -> first_in;
  276. if (IS_ODD(a0_rev_first))
  277. {
  278. a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
  279. a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
  280. }
  281. else a0_rev_last = (i + 1) -> first_in;
  282. for (a0_for = a0_for_first; a0_for < a0_for_last; a0_for++)
  283. if (a0_for->r_rev_cap)
  284. {
  285. j = NEIGHBOR_NODE(i, a0_for -> shift);
  286. if ((a = j->parent) && !j->is_sink)
  287. {
  288. /* checking the origin of j */
  289. d = 0;
  290. while ( 1 )
  291. {
  292. if (j->TS == TIME)
  293. {
  294. d += j -> DIST;
  295. break;
  296. }
  297. a = j -> parent;
  298. d ++;
  299. if (a == TERMINAL)
  300. {
  301. j -> TS = TIME;
  302. j -> DIST = 1;
  303. break;
  304. }
  305. if (a == ORPHAN) {
  306. d = INFINITE_D;
  307. break;
  308. }
  309. if (IS_ODD(a))
  310. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  311. else
  312. j = NEIGHBOR_NODE(j, a -> shift);
  313. }
  314. if (d < INFINITE_D) /* j originates from the source - done */
  315. {
  316. if (d < d_min)
  317. {
  318. a0_min = a0_for;
  319. d_min = d;
  320. }
  321. /* set marks along the path */
  322. for (j = NEIGHBOR_NODE(i, a0_for->shift); j->TS != TIME; )
  323. {
  324. j -> TS = TIME;
  325. j -> DIST = d --;
  326. a = j->parent;
  327. if (IS_ODD(a))
  328. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  329. else
  330. j = NEIGHBOR_NODE(j, a -> shift);
  331. }
  332. }
  333. }
  334. }
  335. for (a0_rev = a0_rev_first; a0_rev < a0_rev_last; a0_rev++)
  336. {
  337. a0_for = a0_rev -> sister;
  338. if (a0_for->r_cap)
  339. {
  340. j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
  341. if (!j->is_sink && (a = j->parent))
  342. {
  343. /* checking the origin of j */
  344. d = 0;
  345. while ( 1 )
  346. {
  347. if (j->TS == TIME)
  348. {
  349. d += j -> DIST;
  350. break;
  351. }
  352. a = j -> parent;
  353. d ++;
  354. if (a == TERMINAL)
  355. {
  356. j -> TS = TIME;
  357. j -> DIST = 1;
  358. break;
  359. }
  360. if (a == ORPHAN) {
  361. d = INFINITE_D;
  362. break;
  363. }
  364. if (IS_ODD(a))
  365. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  366. else
  367. j = NEIGHBOR_NODE(j, a -> shift);
  368. }
  369. if (d < INFINITE_D) /* j originates from the source - done */
  370. {
  371. if (d < d_min)
  372. {
  373. a0_min = MAKE_ODD(a0_for);
  374. d_min = d;
  375. }
  376. /* set marks along the path */
  377. for (j = NEIGHBOR_NODE_REV(i, a0_for->shift); j->TS != TIME; )
  378. {
  379. j -> TS = TIME;
  380. j -> DIST = d --;
  381. a = j->parent;
  382. if (IS_ODD(a))
  383. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  384. else
  385. j = NEIGHBOR_NODE(j, a -> shift);
  386. }
  387. }
  388. }
  389. }
  390. }
  391. if (i->parent = a0_min)
  392. {
  393. i -> TS = TIME;
  394. i -> DIST = d_min + 1;
  395. }
  396. else
  397. {
  398. /* no parent is found */
  399. i -> TS = 0;
  400. /* process neighbors */
  401. for (a0_for = a0_for_first; a0_for < a0_for_last; a0_for++)
  402. {
  403. j = NEIGHBOR_NODE(i, a0_for -> shift);
  404. if (!j->is_sink && (a = j->parent))
  405. {
  406. if (a0_for->r_rev_cap) set_active(j);
  407. if (a != TERMINAL && a != ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift) == i)
  408. {
  409. /* add j to the adoption list */
  410. j -> parent = ORPHAN;
  411. np = nodeptr_block -> New();
  412. np -> ptr = j;
  413. if (orphan_last) orphan_last -> next = np;
  414. else orphan_first = np;
  415. orphan_last = np;
  416. np -> next = NULL;
  417. }
  418. }
  419. }
  420. for (a0_rev = a0_rev_first; a0_rev < a0_rev_last; a0_rev++)
  421. {
  422. a0_for = a0_rev -> sister;
  423. j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
  424. if (!j->is_sink && (a = j->parent))
  425. {
  426. if (a0_for->r_cap) set_active(j);
  427. if (a != TERMINAL && a != ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift) == i)
  428. {
  429. /* add j to the adoption list */
  430. j -> parent = ORPHAN;
  431. np = nodeptr_block -> New();
  432. np -> ptr = j;
  433. if (orphan_last) orphan_last -> next = np;
  434. else orphan_first = np;
  435. orphan_last = np;
  436. np -> next = NULL;
  437. }
  438. }
  439. }
  440. }
  441. }
  442. void Graph::process_sink_orphan(node *i)
  443. {
  444. node *j;
  445. arc_forward *a0_for, *a0_for_first, *a0_for_last;
  446. arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
  447. arc_forward *a0_min = NULL, *a;
  448. nodeptr *np;
  449. int d, d_min = INFINITE_D;
  450. /* trying to find a new parent */
  451. a0_for_first = i -> first_out;
  452. if (IS_ODD(a0_for_first))
  453. {
  454. a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);
  455. a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);
  456. }
  457. else a0_for_last = (i + 1) -> first_out;
  458. a0_rev_first = i -> first_in;
  459. if (IS_ODD(a0_rev_first))
  460. {
  461. a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);
  462. a0_rev_last = (arc_reverse *) ((a0_rev_first ++) -> sister);
  463. }
  464. else a0_rev_last = (i + 1) -> first_in;
  465. for (a0_for = a0_for_first; a0_for < a0_for_last; a0_for++)
  466. if (a0_for->r_cap)
  467. {
  468. j = NEIGHBOR_NODE(i, a0_for -> shift);
  469. if (j->is_sink && (a = j->parent))
  470. {
  471. /* checking the origin of j */
  472. d = 0;
  473. while ( 1 )
  474. {
  475. if (j->TS == TIME)
  476. {
  477. d += j -> DIST;
  478. break;
  479. }
  480. a = j -> parent;
  481. d ++;
  482. if (a == TERMINAL)
  483. {
  484. j -> TS = TIME;
  485. j -> DIST = 1;
  486. break;
  487. }
  488. if (a == ORPHAN) {
  489. d = INFINITE_D;
  490. break;
  491. }
  492. if (IS_ODD(a))
  493. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  494. else
  495. j = NEIGHBOR_NODE(j, a -> shift);
  496. }
  497. if (d < INFINITE_D) /* j originates from the sink - done */
  498. {
  499. if (d < d_min)
  500. {
  501. a0_min = a0_for;
  502. d_min = d;
  503. }
  504. /* set marks along the path */
  505. for (j = NEIGHBOR_NODE(i, a0_for->shift); j->TS != TIME; )
  506. {
  507. j -> TS = TIME;
  508. j -> DIST = d --;
  509. a = j->parent;
  510. if (IS_ODD(a))
  511. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  512. else
  513. j = NEIGHBOR_NODE(j, a -> shift);
  514. }
  515. }
  516. }
  517. }
  518. for (a0_rev = a0_rev_first; a0_rev < a0_rev_last; a0_rev++)
  519. {
  520. a0_for = a0_rev -> sister;
  521. if (a0_for->r_rev_cap)
  522. {
  523. j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
  524. if (j->is_sink && (a = j->parent))
  525. {
  526. /* checking the origin of j */
  527. d = 0;
  528. while ( 1 )
  529. {
  530. if (j->TS == TIME)
  531. {
  532. d += j -> DIST;
  533. break;
  534. }
  535. a = j -> parent;
  536. d ++;
  537. if (a == TERMINAL)
  538. {
  539. j -> TS = TIME;
  540. j -> DIST = 1;
  541. break;
  542. }
  543. if (a == ORPHAN) {
  544. d = INFINITE_D;
  545. break;
  546. }
  547. if (IS_ODD(a))
  548. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  549. else
  550. j = NEIGHBOR_NODE(j, a -> shift);
  551. }
  552. if (d < INFINITE_D) /* j originates from the sink - done */
  553. {
  554. if (d < d_min)
  555. {
  556. a0_min = MAKE_ODD(a0_for);
  557. d_min = d;
  558. }
  559. /* set marks along the path */
  560. for (j = NEIGHBOR_NODE_REV(i, a0_for->shift); j->TS != TIME; )
  561. {
  562. j -> TS = TIME;
  563. j -> DIST = d --;
  564. a = j->parent;
  565. if (IS_ODD(a))
  566. j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
  567. else
  568. j = NEIGHBOR_NODE(j, a -> shift);
  569. }
  570. }
  571. }
  572. }
  573. }
  574. if (i->parent = a0_min)
  575. {
  576. i -> TS = TIME;
  577. i -> DIST = d_min + 1;
  578. }
  579. else
  580. {
  581. /* no parent is found */
  582. i -> TS = 0;
  583. /* process neighbors */
  584. for (a0_for = a0_for_first; a0_for < a0_for_last; a0_for++)
  585. {
  586. j = NEIGHBOR_NODE(i, a0_for -> shift);
  587. if (j->is_sink && (a = j->parent))
  588. {
  589. if (a0_for->r_cap) set_active(j);
  590. if (a != TERMINAL && a != ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift) == i)
  591. {
  592. /* add j to the adoption list */
  593. j -> parent = ORPHAN;
  594. np = nodeptr_block -> New();
  595. np -> ptr = j;
  596. if (orphan_last) orphan_last -> next = np;
  597. else orphan_first = np;
  598. orphan_last = np;
  599. np -> next = NULL;
  600. }
  601. }
  602. }
  603. for (a0_rev = a0_rev_first; a0_rev < a0_rev_last; a0_rev++)
  604. {
  605. a0_for = a0_rev -> sister;
  606. j = NEIGHBOR_NODE_REV(i, a0_for -> shift);
  607. if (j->is_sink && (a = j->parent))
  608. {
  609. if (a0_for->r_rev_cap) set_active(j);
  610. if (a != TERMINAL && a != ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift) == i)
  611. {
  612. /* add j to the adoption list */
  613. j -> parent = ORPHAN;
  614. np = nodeptr_block -> New();
  615. np -> ptr = j;
  616. if (orphan_last) orphan_last -> next = np;
  617. else orphan_first = np;
  618. orphan_last = np;
  619. np -> next = NULL;
  620. }
  621. }
  622. }
  623. }
  624. }
  625. /***********************************************************************/
  626. Graph::flowtype Graph::maxflow()
  627. {
  628. node *i, *j, *current_node = NULL, *s_start, *t_start;
  629. captype *cap_middle, *rev_cap_middle;
  630. arc_forward *a_for, *a_for_first, *a_for_last;
  631. arc_reverse *a_rev, *a_rev_first, *a_rev_last;
  632. nodeptr *np, *np_next;
  633. prepare_graph();
  634. maxflow_init();
  635. nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
  636. while ( 1 )
  637. {
  638. if (i = current_node)
  639. {
  640. i -> next = NULL; /* remove active flag */
  641. if (!i->parent) i = NULL;
  642. }
  643. if (!i)
  644. {
  645. if (!(i = next_active())) break;
  646. }
  647. /* growth */
  648. s_start = NULL;
  649. a_for_first = i -> first_out;
  650. if (IS_ODD(a_for_first))
  651. {
  652. a_for_first = (arc_forward *) (((char *)a_for_first) + 1);
  653. a_for_last = (arc_forward *) ((a_for_first ++) -> shift);
  654. }
  655. else a_for_last = (i + 1) -> first_out;
  656. a_rev_first = i -> first_in;
  657. if (IS_ODD(a_rev_first))
  658. {
  659. a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);
  660. a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);
  661. }
  662. else a_rev_last = (i + 1) -> first_in;
  663. if (!i->is_sink)
  664. {
  665. /* grow source tree */
  666. for (a_for = a_for_first; a_for < a_for_last; a_for++)
  667. if (a_for->r_cap)
  668. {
  669. j = NEIGHBOR_NODE(i, a_for -> shift);
  670. if (!j->parent)
  671. {
  672. j -> is_sink = 0;
  673. j -> parent = MAKE_ODD(a_for);
  674. j -> TS = i -> TS;
  675. j -> DIST = i -> DIST + 1;
  676. set_active(j);
  677. }
  678. else if (j->is_sink)
  679. {
  680. s_start = i;
  681. t_start = j;
  682. cap_middle = & ( a_for -> r_cap );
  683. rev_cap_middle = & ( a_for -> r_rev_cap );
  684. break;
  685. }
  686. else if (j->TS <= i->TS &&
  687. j->DIST > i->DIST)
  688. {
  689. /* heuristic - trying to make the distance from j to the source shorter */
  690. j -> parent = MAKE_ODD(a_for);
  691. j -> TS = i -> TS;
  692. j -> DIST = i -> DIST + 1;
  693. }
  694. }
  695. if (!s_start)
  696. for (a_rev = a_rev_first; a_rev < a_rev_last; a_rev++)
  697. {
  698. a_for = a_rev -> sister;
  699. if (a_for->r_rev_cap)
  700. {
  701. j = NEIGHBOR_NODE_REV(i, a_for -> shift);
  702. if (!j->parent)
  703. {
  704. j -> is_sink = 0;
  705. j -> parent = a_for;
  706. j -> TS = i -> TS;
  707. j -> DIST = i -> DIST + 1;
  708. set_active(j);
  709. }
  710. else if (j->is_sink)
  711. {
  712. s_start = i;
  713. t_start = j;
  714. cap_middle = & ( a_for -> r_rev_cap );
  715. rev_cap_middle = & ( a_for -> r_cap );
  716. break;
  717. }
  718. else if (j->TS <= i->TS &&
  719. j->DIST > i->DIST)
  720. {
  721. /* heuristic - trying to make the distance from j to the source shorter */
  722. j -> parent = a_for;
  723. j -> TS = i -> TS;
  724. j -> DIST = i -> DIST + 1;
  725. }
  726. }
  727. }
  728. }
  729. else
  730. {
  731. /* grow sink tree */
  732. for (a_for = a_for_first; a_for < a_for_last; a_for++)
  733. if (a_for->r_rev_cap)
  734. {
  735. j = NEIGHBOR_NODE(i, a_for -> shift);
  736. if (!j->parent)
  737. {
  738. j -> is_sink = 1;
  739. j -> parent = MAKE_ODD(a_for);
  740. j -> TS = i -> TS;
  741. j -> DIST = i -> DIST + 1;
  742. set_active(j);
  743. }
  744. else if (!j->is_sink)
  745. {
  746. s_start = j;
  747. t_start = i;
  748. cap_middle = & ( a_for -> r_rev_cap );
  749. rev_cap_middle = & ( a_for -> r_cap );
  750. break;
  751. }
  752. else if (j->TS <= i->TS &&
  753. j->DIST > i->DIST)
  754. {
  755. /* heuristic - trying to make the distance from j to the sink shorter */
  756. j -> parent = MAKE_ODD(a_for);
  757. j -> TS = i -> TS;
  758. j -> DIST = i -> DIST + 1;
  759. }
  760. }
  761. for (a_rev = a_rev_first; a_rev < a_rev_last; a_rev++)
  762. {
  763. a_for = a_rev -> sister;
  764. if (a_for->r_cap)
  765. {
  766. j = NEIGHBOR_NODE_REV(i, a_for -> shift);
  767. if (!j->parent)
  768. {
  769. j -> is_sink = 1;
  770. j -> parent = a_for;
  771. j -> TS = i -> TS;
  772. j -> DIST = i -> DIST + 1;
  773. set_active(j);
  774. }
  775. else if (!j->is_sink)
  776. {
  777. s_start = j;
  778. t_start = i;
  779. cap_middle = & ( a_for -> r_cap );
  780. rev_cap_middle = & ( a_for -> r_rev_cap );
  781. break;
  782. }
  783. else if (j->TS <= i->TS &&
  784. j->DIST > i->DIST)
  785. {
  786. /* heuristic - trying to make the distance from j to the sink shorter */
  787. j -> parent = a_for;
  788. j -> TS = i -> TS;
  789. j -> DIST = i -> DIST + 1;
  790. }
  791. }
  792. }
  793. }
  794. TIME ++;
  795. if (s_start)
  796. {
  797. i -> next = i; /* set active flag */
  798. current_node = i;
  799. /* augmentation */
  800. augment(s_start, t_start, cap_middle, rev_cap_middle);
  801. /* augmentation end */
  802. /* adoption */
  803. while (np = orphan_first)
  804. {
  805. np_next = np -> next;
  806. np -> next = NULL;
  807. while (np = orphan_first)
  808. {
  809. orphan_first = np -> next;
  810. i = np -> ptr;
  811. nodeptr_block -> Delete(np);
  812. if (!orphan_first) orphan_last = NULL;
  813. if (i->is_sink) process_sink_orphan(i);
  814. else process_source_orphan(i);
  815. }
  816. orphan_first = np_next;
  817. }
  818. /* adoption end */
  819. }
  820. else current_node = NULL;
  821. }
  822. delete nodeptr_block;
  823. return flow;
  824. }
  825. /***********************************************************************/
  826. Graph::termtype Graph::what_segment(node_id i)
  827. {
  828. if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;
  829. return SINK;
  830. }