Ticket #291: 291nsfixese239cde7ae57.patch
File 291nsfixese239cde7ae57.patch, 29.0 KB (added by , 15 years ago) 


doc/Makefile.am
# HG changeset patch # User Peter Kovacs <kpeter@inf.elte.hu> # Date 1242122800 7200 # Node ID e239cde7ae572839447f31364ff9e67ae86a059a # Parent ca92c2f936b03d0966865203a090abbd2e703cbc Fix the GEQ/LEQ handling in NetworkSimplex + improve doc (#291)  Fix the optimality conditions for the GEQ/LEQ form.  Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms.  Fix the pivot rules to search all the arcs that have to be allowed to get in the basis.  Better block size for the Block Search pivot rule.  Improve documentation of the problem and move it to a separate page. diff git a/doc/Makefile.am b/doc/Makefile.am
a b 8 8 doc/license.dox \ 9 9 doc/mainpage.dox \ 10 10 doc/migration.dox \ 11 doc/min_cost_flow.dox \ 11 12 doc/namedparam.dox \ 12 13 doc/namespaces.dox \ 13 14 doc/html \ 
doc/groups.dox
diff git a/doc/groups.dox b/doc/groups.dox
a b 345 345 */ 346 346 347 347 /** 348 @defgroup min_cost_flow Minimum Cost Flow Algorithms348 @defgroup min_cost_flow_algs Minimum Cost Flow Algorithms 349 349 @ingroup algs 350 350 351 351 \brief Algorithms for finding minimum cost flows and circulations. 352 352 353 353 This group contains the algorithms for finding minimum cost flows and 354 circulations. 354 circulations. For more information about this problem and its dual 355 solution see \ref min_cost_flow "Minimum Cost Flow Problem". 355 356 356 The \e minimum \e cost \e flow \e problem is to find a feasible flow of 357 minimum total cost from a set of supply nodes to a set of demand nodes 358 in a network with capacity constraints (lower and upper bounds) 359 and arc costs. 360 Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, 361 \f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and 362 upper bounds for the flow values on the arcs, for which 363 \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, 364 \f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow 365 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the 366 signed supply values of the nodes. 367 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ 368 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with 369 \f$sup(u)\f$ demand. 370 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution 371 of the following optimization problem. 372 373 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 374 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 375 sup(u) \quad \forall u\in V \f] 376 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 377 378 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be 379 zero or negative in order to have a feasible solution (since the sum 380 of the expressions on the lefthand side of the inequalities is zero). 381 It means that the total demand must be greater or equal to the total 382 supply and all the supplies have to be carried out from the supply nodes, 383 but there could be demands that are not satisfied. 384 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand 385 constraints have to be satisfied with equality, i.e. all demands 386 have to be satisfied and all supplies have to be used. 387 388 If you need the opposite inequalities in the supply/demand constraints 389 (i.e. the total demand is less than the total supply and all the demands 390 have to be satisfied while there could be supplies that are not used), 391 then you could easily transform the problem to the above form by reversing 392 the direction of the arcs and taking the negative of the supply values 393 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). 394 However \ref NetworkSimplex algorithm also supports this form directly 395 for the sake of convenience. 396 397 A feasible solution for this problem can be found using \ref Circulation. 398 399 Note that the above formulation is actually more general than the usual 400 definition of the minimum cost flow problem, in which strict equalities 401 are required in the supply/demand contraints, i.e. 402 403 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) = 404 sup(u) \quad \forall u\in V. \f] 405 406 However if the sum of the supply values is zero, then these two problems 407 are equivalent. So if you need the equality form, you have to ensure this 408 additional contraint for the algorithms. 409 410 The dual solution of the minimum cost flow problem is represented by node 411 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. 412 An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem 413 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ 414 node potentials the following \e complementary \e slackness optimality 415 conditions hold. 416 417  For all \f$uv\in A\f$ arcs: 418  if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; 419  if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 420  if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 421  For all \f$u\in V\f$ nodes: 422  if \f$\sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \neq sup(u)\f$, 423 then \f$\pi(u)=0\f$. 424 425 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc 426 \f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. 427 \f[ cost^\pi(uv) = cost(uv) + \pi(u)  \pi(v).\f] 428 429 All algorithms provide dual solution (node potentials) as well, 430 if an optimal flow is found. 431 432 LEMON contains several algorithms for solving minimum cost flow problems. 357 LEMON contains several algorithms for this problem. 433 358  \ref NetworkSimplex Primal Network Simplex algorithm with various 434 359 pivot strategies. 435 360  \ref CostScaling PushRelabel and AugmentRelabel algorithms based on … … 439 364  \ref CancelAndTighten The Cancel and Tighten algorithm. 440 365  \ref CycleCanceling CycleCanceling algorithms. 441 366 442 Most of these implementations support the general inequality form of the443 minimum cost flow problem, but CancelAndTighten and CycleCanceling444 only support the equality form due to the primal method they use.445 446 367 In general NetworkSimplex is the most efficient implementation, 447 368 but in special cases other algorithms could be faster. 448 369 For example, if the total supply and/or capacities are rather small, 
new file doc/min_cost_flow.dox
diff git a/doc/min_cost_flow.dox b/doc/min_cost_flow.dox new file mode 100644
 + 1 /* * mode: C++; indenttabsmode: nil; * 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 20032009 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). 8 * 9 * Permission to use, modify and distribute this software is granted 10 * provided that this copyright notice appears in all copies. For 11 * precise terms see the accompanying LICENSE file. 12 * 13 * This software is provided "AS IS" with no warranty of any kind, 14 * express or implied, and with no claim as to its suitability for any 15 * purpose. 16 * 17 */ 18 19 namespace lemon { 20 21 /** 22 \page min_cost_flow Minimum Cost Flow Problem 23 24 \section mcf_def Definition 25 26 The \e minimum \e cost \e flow \e problem is to find a feasible flow of 27 minimum total cost from a set of supply nodes to a set of demand nodes 28 in a network with capacity constraints (lower and upper bounds) 29 and arc costs. 30 31 Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$, 32 \f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and 33 upper bounds for the flow values on the arcs, for which 34 \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, 35 \f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow 36 on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the 37 signed supply values of the nodes. 38 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ 39 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with 40 \f$sup(u)\f$ demand. 41 A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution 42 of the following optimization problem. 43 44 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 45 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 46 sup(u) \quad \forall u\in V \f] 47 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 48 49 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be 50 zero or negative in order to have a feasible solution (since the sum 51 of the expressions on the lefthand side of the inequalities is zero). 52 It means that the total demand must be greater or equal to the total 53 supply and all the supplies have to be carried out from the supply nodes, 54 but there could be demands that are not satisfied. 55 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand 56 constraints have to be satisfied with equality, i.e. all demands 57 have to be satisfied and all supplies have to be used. 58 59 60 \section mcf_algs Algorithms 61 62 LEMON contains several algorithms for solving this problem, for more 63 information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms". 64 65 A feasible solution for this problem can be found using \ref Circulation. 66 67 68 \section mcf_dual Dual Solution 69 70 The dual solution of the minimum cost flow problem is represented by 71 node potentials \f$\pi: V\rightarrow\mathbf{R}\f$. 72 An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal 73 if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials 74 the following \e complementary \e slackness optimality conditions hold. 75 76  For all \f$uv\in A\f$ arcs: 77  if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; 78  if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 79  if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 80  For all \f$u\in V\f$ nodes: 81  \f$\pi(u)<=0\f$; 82  if \f$\sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \neq sup(u)\f$, 83 then \f$\pi(u)=0\f$. 84 85 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc 86 \f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. 87 \f[ cost^\pi(uv) = cost(uv) + \pi(u)  \pi(v).\f] 88 89 All algorithms provide dual solution (node potentials), as well, 90 if an optimal flow is found. 91 92 93 \section mcf_eq Equality Form 94 95 The above \ref mcf_def "definition" is actually more general than the 96 usual formulation of the minimum cost flow problem, in which strict 97 equalities are required in the supply/demand contraints. 98 99 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 100 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) = 101 sup(u) \quad \forall u\in V \f] 102 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 103 104 However if the sum of the supply values is zero, then these two problems 105 are equivalent. 106 The \ref min_cost_flow_algs "algorithms" in LEMON support the general 107 form, so if you need the equality form, you have to ensure this additional 108 contraint manually. 109 110 111 \section mcf_leq Opposite Inequality Form 112 113 Another possible definition of the minimum cost flow problem is 114 when there are <em>"less or equal"</em> (LEQ) supply/demand constraints, 115 instead of the <em>"greater or equal"</em> (GEQ) constraints. 116 117 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 118 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \leq 119 sup(u) \quad \forall u\in V \f] 120 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 121 122 It means that the total demand must be less or equal to the 123 total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 124 positive) and all the demands have to be satisfied, but there 125 could be supplies that are not carried out from the supply 126 nodes. 127 The equality form is also a special case of this form, of course. 128 129 You could easily transform this case to the \ref mcf_def "GEQ form" 130 of the problem by reversing the direction of the arcs and taking the 131 negative of the supply values (e.g. using \ref ReverseDigraph and 132 \ref NegMap adaptors). 133 However \ref NetworkSimplex algorithm also supports this form directly 134 for the sake of convenience. 135 136 Note that the optimality conditions for this supply constraint type are 137 slightly differ from the conditions that are discussed for the GEQ form, 138 namely the potentials have to be nonnegative instead of nonpositive. 139 An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem 140 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ 141 node potentials the following conditions hold. 142 143  For all \f$uv\in A\f$ arcs: 144  if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; 145  if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 146  if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 147  For all \f$u\in V\f$ nodes: 148  \f$\pi(u)>=0\f$; 149  if \f$\sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \neq sup(u)\f$, 150 then \f$\pi(u)=0\f$. 151 152 */ 153 } 
lemon/network_simplex.h
diff git a/lemon/network_simplex.h b/lemon/network_simplex.h
a b 19 19 #ifndef LEMON_NETWORK_SIMPLEX_H 20 20 #define LEMON_NETWORK_SIMPLEX_H 21 21 22 /// \ingroup min_cost_flow 22 /// \ingroup min_cost_flow_algs 23 23 /// 24 24 /// \file 25 25 /// \brief Network Simplex algorithm for finding a minimum cost flow. … … 33 33 34 34 namespace lemon { 35 35 36 /// \addtogroup min_cost_flow 36 /// \addtogroup min_cost_flow_algs 37 37 /// @{ 38 38 39 39 /// \brief Implementation of the primal Network Simplex algorithm … … 102 102 /// i.e. the direction of the inequalities in the supply/demand 103 103 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 104 104 /// 105 /// The default supply type is \c GEQ, since this form is supported 106 /// by other minimum cost flow algorithms and the \ref Circulation 107 /// algorithm, as well. 108 /// The \c LEQ problem type can be selected using the \ref supplyType() 109 /// function. 110 /// 111 /// Note that the equality form is a special case of both supply types. 105 /// The default supply type is \c GEQ, the \c LEQ type can be 106 /// selected using \ref supplyType(). 107 /// The equality form is a special case of both supply types. 112 108 enum SupplyType { 113 114 109 /// This option means that there are <em>"greater or equal"</em> 115 /// supply/demand constraints in the definition, i.e. the exact 116 /// formulation of the problem is the following. 117 /** 118 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 119 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 120 sup(u) \quad \forall u\in V \f] 121 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 122 */ 123 /// It means that the total demand must be greater or equal to the 124 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 125 /// negative) and all the supplies have to be carried out from 126 /// the supply nodes, but there could be demands that are not 127 /// satisfied. 110 /// supply/demand constraints in the definition of the problem. 128 111 GEQ, 129 /// It is just an alias for the \c GEQ option.130 CARRY_SUPPLIES = GEQ,131 132 112 /// This option means that there are <em>"less or equal"</em> 133 /// supply/demand constraints in the definition, i.e. the exact 134 /// formulation of the problem is the following. 135 /** 136 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 137 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \leq 138 sup(u) \quad \forall u\in V \f] 139 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 140 */ 141 /// It means that the total demand must be less or equal to the 142 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 143 /// positive) and all the demands have to be satisfied, but there 144 /// could be supplies that are not carried out from the supply 145 /// nodes. 146 LEQ, 147 /// It is just an alias for the \c LEQ option. 148 SATISFY_DEMANDS = LEQ 113 /// supply/demand constraints in the definition of the problem. 114 LEQ 149 115 }; 150 116 151 117 /// \brief Constants for selecting the pivot rule. … … 215 181 const GR &_graph; 216 182 int _node_num; 217 183 int _arc_num; 184 int _all_arc_num; 185 int _search_arc_num; 218 186 219 187 // Parameters of the problem 220 188 bool _have_lower; … … 277 245 const IntVector &_state; 278 246 const CostVector &_pi; 279 247 int &_in_arc; 280 int _ arc_num;248 int _search_arc_num; 281 249 282 250 // Pivot rule data 283 251 int _next_arc; … … 288 256 FirstEligiblePivotRule(NetworkSimplex &ns) : 289 257 _source(ns._source), _target(ns._target), 290 258 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 291 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 259 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 260 _next_arc(0) 292 261 {} 293 262 294 263 // Find next entering arc 295 264 bool findEnteringArc() { 296 265 Cost c; 297 for (int e = _next_arc; e < _ arc_num; ++e) {266 for (int e = _next_arc; e < _search_arc_num; ++e) { 298 267 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 299 268 if (c < 0) { 300 269 _in_arc = e; … … 328 297 const IntVector &_state; 329 298 const CostVector &_pi; 330 299 int &_in_arc; 331 int _ arc_num;300 int _search_arc_num; 332 301 333 302 public: 334 303 … … 336 305 BestEligiblePivotRule(NetworkSimplex &ns) : 337 306 _source(ns._source), _target(ns._target), 338 307 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 339 _in_arc(ns.in_arc), _ arc_num(ns._arc_num)308 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) 340 309 {} 341 310 342 311 // Find next entering arc 343 312 bool findEnteringArc() { 344 313 Cost c, min = 0; 345 for (int e = 0; e < _ arc_num; ++e) {314 for (int e = 0; e < _search_arc_num; ++e) { 346 315 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 347 316 if (c < min) { 348 317 min = c; … … 367 336 const IntVector &_state; 368 337 const CostVector &_pi; 369 338 int &_in_arc; 370 int _ arc_num;339 int _search_arc_num; 371 340 372 341 // Pivot rule data 373 342 int _block_size; … … 379 348 BlockSearchPivotRule(NetworkSimplex &ns) : 380 349 _source(ns._source), _target(ns._target), 381 350 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 382 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 351 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 352 _next_arc(0) 383 353 { 384 354 // The main parameters of the pivot rule 385 const double BLOCK_SIZE_FACTOR = 2.0;355 const double BLOCK_SIZE_FACTOR = 0.5; 386 356 const int MIN_BLOCK_SIZE = 10; 387 357 388 358 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 389 std::sqrt(double(_ arc_num))),359 std::sqrt(double(_search_arc_num))), 390 360 MIN_BLOCK_SIZE ); 391 361 } 392 362 … … 395 365 Cost c, min = 0; 396 366 int cnt = _block_size; 397 367 int e, min_arc = _next_arc; 398 for (e = _next_arc; e < _ arc_num; ++e) {368 for (e = _next_arc; e < _search_arc_num; ++e) { 399 369 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 400 370 if (c < min) { 401 371 min = c; … … 440 410 const IntVector &_state; 441 411 const CostVector &_pi; 442 412 int &_in_arc; 443 int _ arc_num;413 int _search_arc_num; 444 414 445 415 // Pivot rule data 446 416 IntVector _candidates; … … 454 424 CandidateListPivotRule(NetworkSimplex &ns) : 455 425 _source(ns._source), _target(ns._target), 456 426 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 457 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 427 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 428 _next_arc(0) 458 429 { 459 430 // The main parameters of the pivot rule 460 431 const double LIST_LENGTH_FACTOR = 1.0; … … 463 434 const int MIN_MINOR_LIMIT = 3; 464 435 465 436 _list_length = std::max( int(LIST_LENGTH_FACTOR * 466 std::sqrt(double(_ arc_num))),437 std::sqrt(double(_search_arc_num))), 467 438 MIN_LIST_LENGTH ); 468 439 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), 469 440 MIN_MINOR_LIMIT ); … … 500 471 // Major iteration: build a new candidate list 501 472 min = 0; 502 473 _curr_length = 0; 503 for (e = _next_arc; e < _ arc_num; ++e) {474 for (e = _next_arc; e < _search_arc_num; ++e) { 504 475 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 505 476 if (c < 0) { 506 477 _candidates[_curr_length++] = e; … … 546 517 const IntVector &_state; 547 518 const CostVector &_pi; 548 519 int &_in_arc; 549 int _ arc_num;520 int _search_arc_num; 550 521 551 522 // Pivot rule data 552 523 int _block_size, _head_length, _curr_length; … … 574 545 AlteringListPivotRule(NetworkSimplex &ns) : 575 546 _source(ns._source), _target(ns._target), 576 547 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 577 _in_arc(ns.in_arc), _ arc_num(ns._arc_num),578 _next_arc(0), _cand_cost(ns._ arc_num), _sort_func(_cand_cost)548 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 549 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) 579 550 { 580 551 // The main parameters of the pivot rule 581 552 const double BLOCK_SIZE_FACTOR = 1.5; … … 584 555 const int MIN_HEAD_LENGTH = 3; 585 556 586 557 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 587 std::sqrt(double(_ arc_num))),558 std::sqrt(double(_search_arc_num))), 588 559 MIN_BLOCK_SIZE ); 589 560 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), 590 561 MIN_HEAD_LENGTH ); … … 610 581 int last_arc = 0; 611 582 int limit = _head_length; 612 583 613 for (int e = _next_arc; e < _ arc_num; ++e) {584 for (int e = _next_arc; e < _search_arc_num; ++e) { 614 585 _cand_cost[e] = _state[e] * 615 586 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 616 587 if (_cand_cost[e] < 0) { … … 678 649 _node_num = countNodes(_graph); 679 650 _arc_num = countArcs(_graph); 680 651 int all_node_num = _node_num + 1; 681 int all_arc_num = _arc_num +_node_num;652 int max_arc_num = _arc_num + 2 * _node_num; 682 653 683 _source.resize( all_arc_num);684 _target.resize( all_arc_num);654 _source.resize(max_arc_num); 655 _target.resize(max_arc_num); 685 656 686 _lower.resize( all_arc_num);687 _upper.resize( all_arc_num);688 _cap.resize( all_arc_num);689 _cost.resize( all_arc_num);657 _lower.resize(_arc_num); 658 _upper.resize(_arc_num); 659 _cap.resize(max_arc_num); 660 _cost.resize(max_arc_num); 690 661 _supply.resize(all_node_num); 691 _flow.resize( all_arc_num);662 _flow.resize(max_arc_num); 692 663 _pi.resize(all_node_num); 693 664 694 665 _parent.resize(all_node_num); … … 698 669 _rev_thread.resize(all_node_num); 699 670 _succ_num.resize(all_node_num); 700 671 _last_succ.resize(all_node_num); 701 _state.resize( all_arc_num);672 _state.resize(max_arc_num); 702 673 703 674 // Copy the graph (store the arcs in a mixed order) 704 675 int i = 0; … … 1069 1040 // Initialize artifical cost 1070 1041 Cost ART_COST; 1071 1042 if (std::numeric_limits<Cost>::is_exact) { 1072 ART_COST = std::numeric_limits<Cost>::max() / 4+ 1;1043 ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; 1073 1044 } else { 1074 1045 ART_COST = std::numeric_limits<Cost>::min(); 1075 1046 for (int i = 0; i != _arc_num; ++i) { … … 1093 1064 _succ_num[_root] = _node_num + 1; 1094 1065 _last_succ[_root] = _root  1; 1095 1066 _supply[_root] = _sum_supply; 1096 _pi[_root] = _sum_supply < 0 ? ART_COST : ART_COST;1067 _pi[_root] = 0; 1097 1068 1098 1069 // Add artificial arcs and initialize the spanning tree data structure 1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root; 1101 _pred[u] = e; 1102 _thread[u] = u + 1; 1103 _rev_thread[u + 1] = u; 1104 _succ_num[u] = 1; 1105 _last_succ[u] = u; 1106 _cost[e] = ART_COST; 1107 _cap[e] = INF; 1108 _state[e] = STATE_TREE; 1109 if (_supply[u] > 0  (_supply[u] == 0 && _sum_supply <= 0)) { 1110 _flow[e] = _supply[u]; 1111 _forward[u] = true; 1112 _pi[u] = ART_COST + _pi[_root]; 1113 } else { 1114 _flow[e] = _supply[u]; 1115 _forward[u] = false; 1116 _pi[u] = ART_COST + _pi[_root]; 1070 if (_sum_supply == 0) { 1071 // EQ supply constraints 1072 _search_arc_num = _arc_num; 1073 _all_arc_num = _arc_num + _node_num; 1074 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1075 _parent[u] = _root; 1076 _pred[u] = e; 1077 _thread[u] = u + 1; 1078 _rev_thread[u + 1] = u; 1079 _succ_num[u] = 1; 1080 _last_succ[u] = u; 1081 _cap[e] = INF; 1082 _state[e] = STATE_TREE; 1083 if (_supply[u] >= 0) { 1084 _forward[u] = true; 1085 _pi[u] = 0; 1086 _source[e] = u; 1087 _target[e] = _root; 1088 _flow[e] = _supply[u]; 1089 _cost[e] = 0; 1090 } else { 1091 _forward[u] = false; 1092 _pi[u] = ART_COST; 1093 _source[e] = _root; 1094 _target[e] = u; 1095 _flow[e] = _supply[u]; 1096 _cost[e] = ART_COST; 1097 } 1117 1098 } 1118 1099 } 1100 else if (_sum_supply > 0) { 1101 // LEQ supply constraints 1102 _search_arc_num = _arc_num + _node_num; 1103 int f = _arc_num + _node_num; 1104 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1105 _parent[u] = _root; 1106 _thread[u] = u + 1; 1107 _rev_thread[u + 1] = u; 1108 _succ_num[u] = 1; 1109 _last_succ[u] = u; 1110 if (_supply[u] >= 0) { 1111 _forward[u] = true; 1112 _pi[u] = 0; 1113 _pred[u] = e; 1114 _source[e] = u; 1115 _target[e] = _root; 1116 _cap[e] = INF; 1117 _flow[e] = _supply[u]; 1118 _cost[e] = 0; 1119 _state[e] = STATE_TREE; 1120 } else { 1121 _forward[u] = false; 1122 _pi[u] = ART_COST; 1123 _pred[u] = f; 1124 _source[f] = _root; 1125 _target[f] = u; 1126 _cap[f] = INF; 1127 _flow[f] = _supply[u]; 1128 _cost[f] = ART_COST; 1129 _state[f] = STATE_TREE; 1130 _source[e] = u; 1131 _target[e] = _root; 1132 _cap[e] = INF; 1133 _flow[e] = 0; 1134 _cost[e] = 0; 1135 _state[e] = STATE_LOWER; 1136 ++f; 1137 } 1138 } 1139 _all_arc_num = f; 1140 } 1141 else { 1142 // GEQ supply constraints 1143 _search_arc_num = _arc_num + _node_num; 1144 int f = _arc_num + _node_num; 1145 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1146 _parent[u] = _root; 1147 _thread[u] = u + 1; 1148 _rev_thread[u + 1] = u; 1149 _succ_num[u] = 1; 1150 _last_succ[u] = u; 1151 if (_supply[u] <= 0) { 1152 _forward[u] = false; 1153 _pi[u] = 0; 1154 _pred[u] = e; 1155 _source[e] = _root; 1156 _target[e] = u; 1157 _cap[e] = INF; 1158 _flow[e] = _supply[u]; 1159 _cost[e] = 0; 1160 _state[e] = STATE_TREE; 1161 } else { 1162 _forward[u] = true; 1163 _pi[u] = ART_COST; 1164 _pred[u] = f; 1165 _source[f] = u; 1166 _target[f] = _root; 1167 _cap[f] = INF; 1168 _flow[f] = _supply[u]; 1169 _state[f] = STATE_TREE; 1170 _cost[f] = ART_COST; 1171 _source[e] = _root; 1172 _target[e] = u; 1173 _cap[e] = INF; 1174 _flow[e] = 0; 1175 _cost[e] = 0; 1176 _state[e] = STATE_LOWER; 1177 ++f; 1178 } 1179 } 1180 _all_arc_num = f; 1181 } 1119 1182 1120 1183 return true; 1121 1184 } … … 1374 1437 } 1375 1438 1376 1439 // Check feasibility 1377 if (_sum_supply < 0) { 1378 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1379 if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; 1380 } 1381 } 1382 else if (_sum_supply > 0) { 1383 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1384 if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; 1385 } 1386 } 1387 else { 1388 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1389 if (_flow[e] != 0) return INFEASIBLE; 1390 } 1440 for (int e = _search_arc_num; e != _all_arc_num; ++e) { 1441 if (_flow[e] != 0) return INFEASIBLE; 1391 1442 } 1392 1443 1393 1444 // Transform the solution and the supply map to the original form … … 1401 1452 } 1402 1453 } 1403 1454 } 1455 1456 // Shift potentials to meet the requirements of the GEQ/LEQ type 1457 // optimality conditions 1458 if (_sum_supply == 0) { 1459 if (_stype == GEQ) { 1460 Cost max_pot = std::numeric_limits<Cost>::min(); 1461 for (int i = 0; i != _node_num; ++i) { 1462 if (_pi[i] > max_pot) max_pot = _pi[i]; 1463 } 1464 if (max_pot > 0) { 1465 for (int i = 0; i != _node_num; ++i) 1466 _pi[i] = max_pot; 1467 } 1468 } else { 1469 Cost min_pot = std::numeric_limits<Cost>::max(); 1470 for (int i = 0; i != _node_num; ++i) { 1471 if (_pi[i] < min_pot) min_pot = _pi[i]; 1472 } 1473 if (min_pot < 0) { 1474 for (int i = 0; i != _node_num; ++i) 1475 _pi[i] = min_pot; 1476 } 1477 } 1478 } 1404 1479 1405 1480 return OPTIMAL; 1406 1481 }