Ticket #219: nsgeqleqe6927fe719e6.patch
File nsgeqleqe6927fe719e6.patch, 28.1 KB (added by , 10 years ago) 


doc/groups.dox
# HG changeset patch # User Peter Kovacs <kpeter@inf.elte.hu> # Date 1239984276 7200 # Node ID e6927fe719e6e2cc97345c5f89f1cb0d66d2d0c9 # Parent 6ac5d9ae1d3dfdfcc48802e03bc2160d5881a78c Support >= and <= constraints in NetworkSimplex (#219, #234) By default the same inequality constraints are supported as by Circulation (the GEQ form), but the LEQ form can also be selected using the problemType() function. The documentation of the min. cost flow module is reworked and extended with important notes and explanations about the different variants of the problem and about the dual solution and optimality conditions. diff git a/doc/groups.dox b/doc/groups.dox
a b 317 317 318 318 The \e maximum \e flow \e problem is to find a flow of maximum value between 319 319 a single source and a single target. Formally, there is a \f$G=(V,A)\f$ 320 digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and320 digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and 321 321 \f$s, t \in V\f$ source and target nodes. 322 A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the322 A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the 323 323 following optimization problem. 324 324 325 \f[ \max\sum_{ a\in\delta_{out}(s)}f(a)  \sum_{a\in\delta_{in}(s)}f(a) \f]326 \f[ \sum_{ a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)327 \q quad \forall v\in V\setminus\{s,t\} \f]328 \f[ 0 \leq f( a) \leq cap(a) \qquad \forall a\in A \f]325 \f[ \max\sum_{sv\in A} f(sv)  \sum_{vs\in A} f(vs) \f] 326 \f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu) 327 \quad \forall u\in V\setminus\{s,t\} \f] 328 \f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f] 329 329 330 330 LEMON contains several algorithms for solving maximum flow problems: 331 331  \ref EdmondsKarp EdmondsKarp algorithm. … … 345 345 346 346 \brief Algorithms for finding minimum cost flows and circulations. 347 347 348 This group describes the algorithms for finding minimum cost flows and348 This group contains the algorithms for finding minimum cost flows and 349 349 circulations. 350 350 351 351 The \e minimum \e cost \e flow \e problem is to find a feasible flow of 352 352 minimum total cost from a set of supply nodes to a set of demand nodes 353 in a network with capacity constraints and arc costs. 353 in a network with capacity constraints (lower and upper bounds) 354 and arc costs. 354 355 Formally, let \f$G=(V,A)\f$ be a digraph, 355 356 \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and 356 upper bounds for the flow values on the arcs, 357 upper bounds for the flow values on the arcs, for which 358 \f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$. 357 359 \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow 358 on the arcs, and 359 \f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values 360 of the nodes. 361 A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of 362 the following optimization problem. 360 on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the 361 signed supply values of the nodes. 362 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ 363 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with 364 \f$sup(u)\f$ demand. 365 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution 366 of the following optimization problem. 363 367 364 \f[ \min\sum_{ a\in A} f(a) cost(a) \f]365 \f[ \sum_{ a\in\delta_{out}(v)} f(a)  \sum_{a\in\delta_{in}(v)} f(a) =366 sup ply(v) \qquad \forall v\in V \f]367 \f[ lower( a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]368 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 369 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 370 sup(u) \quad \forall u\in V \f] 371 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 368 372 369 LEMON contains several algorithms for solving minimum cost flow problems: 370  \ref CycleCanceling Cyclecanceling algorithms. 371  \ref CapacityScaling Successive shortest path algorithm with optional 373 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be 374 zero or negative in order to have a feasible solution (since the sum 375 of the expressions on the lefthand side of the inequalities is zero). 376 It means that the total demand must be greater or equal to the total 377 supply and all the supplies have to be carried out from the supply nodes, 378 but there could be demands that are not satisfied. 379 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand 380 constraints have to be satisfied with equality, i.e. all demands 381 have to be satisfied and all supplies have to be used. 382 383 If you need the opposite inequalities in the supply/demand constraints 384 (i.e. the total demand is less than the total supply and all the demands 385 have to be satisfied while there could be supplies that are not used), 386 then you could easily transform the problem to the above form by reversing 387 the direction of the arcs and taking the negative of the supply values 388 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). 389 However \ref NetworkSimplex algorithm also supports this form directly 390 for the sake of convenience. 391 392 A feasible solution for this problem can be found using \ref Circulation. 393 394 Note that the above formulation is actually more general than the usual 395 definition of the minimum cost flow problem, in which strict equalities 396 are required in the supply/demand contraints, i.e. 397 398 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) = 399 sup(u) \quad \forall u\in V. \f] 400 401 However if the sum of the supply values is zero, then these two problems 402 are equivalent. So if you need the equality form, you have to ensure this 403 additional contraint for the algorithms. 404 405 The dual solution of the minimum cost flow problem is represented by node 406 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. 407 An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem 408 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ 409 node potentials the following \e complementary \e slackness optimality 410 conditions hold. 411 412  For all \f$uv\in A\f$ arcs: 413  if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; 414  if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; 415  if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. 416  For all \f$u\in V\f$: 417  if \f$\sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \neq sup(u)\f$, 418 then \f$\pi(u)=0\f$. 419 420 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc 421 \f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e. 422 \f[ cost^\pi(uv) = cost(uv) + \pi(u)  \pi(v).\f] 423 424 All algorithms provide dual solution (node potentials) as well 425 if an optimal flow is found. 426 427 LEMON contains several algorithms for solving minimum cost flow problems. 428  \ref NetworkSimplex Primal Network Simplex algorithm with various 429 pivot strategies. 430  \ref CostScaling PushRelabel and AugmentRelabel algorithms based on 431 cost scaling. 432  \ref CapacityScaling Successive Shortest %Path algorithm with optional 372 433 capacity scaling. 373  \ref CostScaling Pushrelabel and augmentrelabel algorithms based on 374 cost scaling. 375  \ref NetworkSimplex Primal network simplex algorithm with various 376 pivot strategies. 434  \ref CancelAndTighten The Cancel and Tighten algorithm. 435  \ref CycleCanceling CycleCanceling algorithms. 436 437 Most of these implementations support the general inequality form of the 438 minimum cost flow problem, but CancelAndTighten and CycleCanceling 439 only support the equality form due to the primal method they use. 440 441 In general NetworkSimplex is the most efficient implementation, 442 but in special cases other algorithms could be faster. 443 For example, if the total supply and/or capacities are rather small, 444 CapacityScaling is usually the fastest algorithm (without effective scaling). 377 445 */ 378 446 379 447 /** 
lemon/network_simplex.h
diff git a/lemon/network_simplex.h b/lemon/network_simplex.h
a b 30 30 31 31 #include <lemon/core.h> 32 32 #include <lemon/math.h> 33 #include <lemon/maps.h> 34 #include <lemon/circulation.h> 35 #include <lemon/adaptors.h> 33 36 34 37 namespace lemon { 35 38 … … 47 50 /// 48 51 /// In general this class is the fastest implementation available 49 52 /// in LEMON for the minimum cost flow problem. 53 /// Moreover it supports both direction of the supply/demand inequality 54 /// constraints. For more information see \ref ProblemType. 50 55 /// 51 56 /// \tparam GR The digraph type the algorithm runs on. 52 57 /// \tparam F The value type used for flow amounts, capacity bounds … … 58 63 /// be integer. 59 64 /// 60 65 /// \note %NetworkSimplex provides five different pivot rule 61 /// implementations. For more information see \ref PivotRule. 66 /// implementations, from which the most efficient one is used 67 /// by default. For more information see \ref PivotRule. 62 68 template <typename GR, typename F = int, typename C = F> 63 69 class NetworkSimplex 64 70 { … … 68 74 typedef F Flow; 69 75 /// The cost type of the algorithm 70 76 typedef C Cost; 77 #ifdef DOXYGEN 78 /// The type of the flow map 79 typedef GR::ArcMap<Flow> FlowMap; 80 /// The type of the potential map 81 typedef GR::NodeMap<Cost> PotentialMap; 82 #else 71 83 /// The type of the flow map 72 84 typedef typename GR::template ArcMap<Flow> FlowMap; 73 85 /// The type of the potential map 74 86 typedef typename GR::template NodeMap<Cost> PotentialMap; 87 #endif 75 88 76 89 public: 77 90 … … 117 130 /// candidate list and extends this list in every iteration. 118 131 ALTERING_LIST 119 132 }; 133 134 /// \brief Enum type for selecting the problem type. 135 /// 136 /// Enum type for selecting the problem type, i.e. the direction of 137 /// the inequalities in the supply/demand constraints of the 138 /// \ref min_cost_flow "minimum cost flow problem". 139 /// 140 /// The default problem type is \c GEQ, since this form is supported 141 /// by other minimum cost flow algorithms and the \ref Circulation 142 /// algorithm as well. 143 /// The \c LEQ problem type can be selected using the \ref problemType() 144 /// function. 145 /// 146 /// Note that the equality form is a special case of both problem type. 147 enum ProblemType { 148 149 /// This option means that there are "<em>greater or equal</em>" 150 /// constraints in the defintion, i.e. the exact formulation of the 151 /// problem is the following. 152 /** 153 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 154 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 155 sup(u) \quad \forall u\in V \f] 156 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 157 */ 158 /// It means that the total demand must be greater or equal to the 159 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 160 /// negative) and all the supplies have to be carried out from 161 /// the supply nodes, but there could be demands that are not 162 /// satisfied. 163 GEQ, 164 /// It is just an alias for the \c GEQ option. 165 CARRY_SUPPLIES = GEQ, 166 167 /// This option means that there are "<em>less or equal</em>" 168 /// constraints in the defintion, i.e. the exact formulation of the 169 /// problem is the following. 170 /** 171 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 172 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \leq 173 sup(u) \quad \forall u\in V \f] 174 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 175 */ 176 /// It means that the total demand must be less or equal to the 177 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 178 /// positive) and all the demands have to be satisfied, but there 179 /// could be supplies that are not carried out from the supply 180 /// nodes. 181 LEQ, 182 /// It is just an alias for the \c LEQ option. 183 SATISFY_DEMANDS = LEQ 184 }; 120 185 121 186 private: 122 187 … … 155 220 bool _pstsup; 156 221 Node _psource, _ptarget; 157 222 Flow _pstflow; 223 ProblemType _ptype; 158 224 159 225 // Result maps 160 226 FlowMap *_flow_map; … … 586 652 587 653 /// \brief Constructor. 588 654 /// 589 /// Constructor.655 /// The constructor of the class. 590 656 /// 591 657 /// \param graph The digraph the algorithm runs on. 592 658 NetworkSimplex(const GR& graph) : 593 659 _graph(graph), 594 660 _plower(NULL), _pupper(NULL), _pcost(NULL), 595 _psupply(NULL), _pstsup(false), 661 _psupply(NULL), _pstsup(false), _ptype(GEQ), 596 662 _flow_map(NULL), _potential_map(NULL), 597 663 _local_flow(false), _local_potential(false), 598 664 _node_id(graph) … … 611 677 if (_local_potential) delete _potential_map; 612 678 } 613 679 680 /// \name Parameters 681 /// The parameters of the algorithm can be specified using these 682 /// functions. 683 684 /// @{ 685 614 686 /// \brief Set the lower bounds on the arcs. 615 687 /// 616 688 /// This function sets the lower bounds on the arcs. … … 760 832 _pstflow = k; 761 833 return *this; 762 834 } 835 836 /// \brief Set the problem type. 837 /// 838 /// This function sets the problem type for the algorithm. 839 /// If it is not used before calling \ref run(), the \ref GEQ problem 840 /// type will be used. 841 /// 842 /// For more information see \ref ProblemType. 843 /// 844 /// \return <tt>(*this)</tt> 845 NetworkSimplex& problemType(ProblemType problem_type) { 846 _ptype = problem_type; 847 return *this; 848 } 763 849 764 850 /// \brief Set the flow map. 765 851 /// … … 795 881 _potential_map = ↦ 796 882 return *this; 797 883 } 884 885 /// @} 798 886 799 887 /// \name Execution Control 800 888 /// The algorithm can be executed using \ref run(). … … 804 892 /// \brief Run the algorithm. 805 893 /// 806 894 /// This function runs the algorithm. 807 /// The paramters can be specified using \ref lowerMap(),895 /// The paramters can be specified using functions \ref lowerMap(), 808 896 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 809 /// \ref costMap(), \ref supplyMap() and \ref stSupply() 810 /// functions. For example, 897 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 898 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). 899 /// For example, 811 900 /// \code 812 901 /// NetworkSimplex<ListDigraph> ns(graph); 813 902 /// ns.boundMaps(lower, upper).costMap(cost) … … 830 919 /// \brief Reset all the parameters that have been given before. 831 920 /// 832 921 /// This function resets all the paramaters that have been given 833 /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(), 834 /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and 835 /// \ref stSupply() functions before. 922 /// before using functions \ref lowerMap(), \ref upperMap(), 923 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), 924 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 925 /// \ref flowMap() and \ref potentialMap(). 836 926 /// 837 927 /// It is useful for multiple run() calls. If this function is not 838 928 /// used, all the parameters given before are kept for the next … … 869 959 _pcost = NULL; 870 960 _psupply = NULL; 871 961 _pstsup = false; 962 _ptype = GEQ; 963 if (_local_flow) delete _flow_map; 964 if (_local_potential) delete _potential_map; 965 _flow_map = NULL; 966 _potential_map = NULL; 967 _local_flow = false; 968 _local_potential = false; 969 872 970 return *this; 873 971 } 874 972 … … 1000 1098 1001 1099 // Initialize node related data 1002 1100 bool valid_supply = true; 1101 Flow sum_supply = 0; 1003 1102 if (!_pstsup && !_psupply) { 1004 1103 _pstsup = true; 1005 1104 _psource = _ptarget = NodeIt(_graph); 1006 1105 _pstflow = 0; 1007 1106 } 1008 1107 if (_psupply) { 1009 Flow sum = 0;1010 1108 int i = 0; 1011 1109 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1012 1110 _node_id[n] = i; 1013 1111 _supply[i] = (*_psupply)[n]; 1014 sum += _supply[i];1112 sum_supply += _supply[i]; 1015 1113 } 1016 valid_supply = (sum == 0); 1114 valid_supply = (_ptype == GEQ && sum_supply <= 0)  1115 (_ptype == LEQ && sum_supply >= 0); 1017 1116 } else { 1018 1117 int i = 0; 1019 1118 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { … … 1021 1120 _supply[i] = 0; 1022 1121 } 1023 1122 _supply[_node_id[_psource]] = _pstflow; 1024 _supply[_node_id[_ptarget]] 1123 _supply[_node_id[_ptarget]] = _pstflow; 1025 1124 } 1026 1125 if (!valid_supply) return false; 1027 1126 1127 // Infinite capacity value 1128 Flow inf_cap = 1129 std::numeric_limits<Flow>::has_infinity ? 1130 std::numeric_limits<Flow>::infinity() : 1131 std::numeric_limits<Flow>::max(); 1132 1133 // Initialize artifical cost 1134 Cost art_cost; 1135 if (std::numeric_limits<Cost>::is_exact) { 1136 art_cost = std::numeric_limits<Cost>::max() / 4 + 1; 1137 } else { 1138 art_cost = std::numeric_limits<Cost>::min(); 1139 for (int i = 0; i != _arc_num; ++i) { 1140 if (_cost[i] > art_cost) art_cost = _cost[i]; 1141 } 1142 art_cost = (art_cost + 1) * _node_num; 1143 } 1144 1145 // Run Circulation to check if a feasible solution exists 1146 typedef ConstMap<Arc, Flow> ConstArcMap; 1147 FlowNodeMap *csup = NULL; 1148 bool local_csup = false; 1149 if (_psupply) { 1150 csup = _psupply; 1151 } else { 1152 csup = new FlowNodeMap(_graph, 0); 1153 (*csup)[_psource] = _pstflow; 1154 (*csup)[_ptarget] = _pstflow; 1155 local_csup = true; 1156 } 1157 bool circ_result = false; 1158 if (_ptype == GEQ  (_ptype == LEQ && sum_supply == 0)) { 1159 // GEQ problem type 1160 if (_plower) { 1161 if (_pupper) { 1162 Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap> 1163 circ(_graph, *_plower, *_pupper, *csup); 1164 circ_result = circ.run(); 1165 } else { 1166 Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap> 1167 circ(_graph, *_plower, ConstArcMap(inf_cap), *csup); 1168 circ_result = circ.run(); 1169 } 1170 } else { 1171 if (_pupper) { 1172 Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap> 1173 circ(_graph, ConstArcMap(0), *_pupper, *csup); 1174 circ_result = circ.run(); 1175 } else { 1176 Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap> 1177 circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup); 1178 circ_result = circ.run(); 1179 } 1180 } 1181 } else { 1182 // LEQ problem type 1183 typedef ReverseDigraph<const GR> RevGraph; 1184 typedef NegMap<FlowNodeMap> NegNodeMap; 1185 RevGraph rgraph(_graph); 1186 NegNodeMap neg_csup(*csup); 1187 if (_plower) { 1188 if (_pupper) { 1189 Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap> 1190 circ(rgraph, *_plower, *_pupper, neg_csup); 1191 circ_result = circ.run(); 1192 } else { 1193 Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap> 1194 circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup); 1195 circ_result = circ.run(); 1196 } 1197 } else { 1198 if (_pupper) { 1199 Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap> 1200 circ(rgraph, ConstArcMap(0), *_pupper, neg_csup); 1201 circ_result = circ.run(); 1202 } else { 1203 Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap> 1204 circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup); 1205 circ_result = circ.run(); 1206 } 1207 } 1208 } 1209 if (local_csup) delete csup; 1210 if (!circ_result) return false; 1211 1028 1212 // Set data for the artificial root node 1029 1213 _root = _node_num; 1030 1214 _parent[_root] = 1; … … 1033 1217 _rev_thread[0] = _root; 1034 1218 _succ_num[_root] = all_node_num; 1035 1219 _last_succ[_root] = _root  1; 1036 _supply[_root] = 0; 1037 _pi[_root] = 0; 1220 _supply[_root] = sum_supply; 1221 if (sum_supply < 0) { 1222 _pi[_root] = art_cost; 1223 } else { 1224 _pi[_root] = art_cost; 1225 } 1038 1226 1039 1227 // Store the arcs in a mixed order 1040 1228 int k = std::max(int(sqrt(_arc_num)), 10); … … 1045 1233 } 1046 1234 1047 1235 // Initialize arc maps 1048 Flow inf_cap =1049 std::numeric_limits<Flow>::has_infinity ?1050 std::numeric_limits<Flow>::infinity() :1051 std::numeric_limits<Flow>::max();1052 1236 if (_pupper && _pcost) { 1053 1237 for (int i = 0; i != _arc_num; ++i) { 1054 1238 Arc e = _arc_ref[i]; … … 1083 1267 } 1084 1268 } 1085 1269 1086 // Initialize artifical cost1087 Cost art_cost;1088 if (std::numeric_limits<Cost>::is_exact) {1089 art_cost = std::numeric_limits<Cost>::max() / 4 + 1;1090 } else {1091 art_cost = std::numeric_limits<Cost>::min();1092 for (int i = 0; i != _arc_num; ++i) {1093 if (_cost[i] > art_cost) art_cost = _cost[i];1094 }1095 art_cost = (art_cost + 1) * _node_num;1096 }1097 1098 1270 // Remove nonzero lower bounds 1099 1271 if (_plower) { 1100 1272 for (int i = 0; i != _arc_num; ++i) { … … 1118 1290 _cost[e] = art_cost; 1119 1291 _cap[e] = inf_cap; 1120 1292 _state[e] = STATE_TREE; 1121 if (_supply[u] > = 0) {1293 if (_supply[u] > 0  (_supply[u] == 0 && sum_supply <= 0)) { 1122 1294 _flow[e] = _supply[u]; 1123 1295 _forward[u] = true; 1124 _pi[u] = art_cost ;1296 _pi[u] = art_cost + _pi[_root]; 1125 1297 } else { 1126 1298 _flow[e] = _supply[u]; 1127 1299 _forward[u] = false; 1128 _pi[u] = art_cost ;1300 _pi[u] = art_cost + _pi[_root]; 1129 1301 } 1130 1302 } 1131 1303 … … 1382 1554 } 1383 1555 } 1384 1556 1385 // Check if the flow amount equals zero on all the artificial arcs1386 for (int e = _arc_num; e != _arc_num + _node_num; ++e) {1387 if (_flow[e] > 0) return false;1388 }1389 1390 1557 // Copy flow values to _flow_map 1391 1558 if (_plower) { 1392 1559 for (int i = 0; i != _arc_num; ++i) { 
test/min_cost_flow_test.cc
diff git a/test/min_cost_flow_test.cc b/test/min_cost_flow_test.cc
a b 33 33 34 34 char test_lgf[] = 35 35 "@nodes\n" 36 "label sup1 sup2 sup3 \n"37 " 1 20 27 0 \n"38 " 2 4 0 0 \n"39 " 3 0 0 0 \n"40 " 4 0 0 0 \n"41 " 5 9 0 0 \n"42 " 6 6 0 0 \n"43 " 7 0 0 0 \n"44 " 8 0 0 0 \n"45 " 9 3 0 0 \n"46 " 10 2 0 0 \n"47 " 11 0 0 0 \n"48 " 12 20 27 0 \n"36 "label sup1 sup2 sup3 sup4 sup5\n" 37 " 1 20 27 0 20 30\n" 38 " 2 4 0 0 8 3\n" 39 " 3 0 0 0 0 0\n" 40 " 4 0 0 0 0 0\n" 41 " 5 9 0 0 6 11\n" 42 " 6 6 0 0 5 6\n" 43 " 7 0 0 0 0 0\n" 44 " 8 0 0 0 0 3\n" 45 " 9 3 0 0 0 0\n" 46 " 10 2 0 0 7 2\n" 47 " 11 0 0 0 10 0\n" 48 " 12 20 27 0 30 20\n" 49 49 "\n" 50 50 "@arcs\n" 51 51 " cost cap low1 low2\n" … … 76 76 "target 12\n"; 77 77 78 78 79 enum ProblemType { 80 EQ, 81 GEQ, 82 LEQ 83 }; 84 79 85 // Check the interface of an MCF algorithm 80 86 template <typename GR, typename Flow, typename Cost> 81 87 class McfClassConcept … … 97 103 .costMap(cost) 98 104 .supplyMap(sup) 99 105 .stSupply(n, n, k) 106 .flowMap(flow) 107 .potentialMap(pot) 100 108 .run(); 109 110 const MCF& const_mcf = mcf; 101 111 102 const typename MCF::FlowMap &fm = mcf.flowMap();103 const typename MCF::PotentialMap &pm = mcf.potentialMap();112 const typename MCF::FlowMap &fm = const_mcf.flowMap(); 113 const typename MCF::PotentialMap &pm = const_mcf.potentialMap(); 104 114 105 v = mcf.totalCost(); 106 double x = mcf.template totalCost<double>(); 107 v = mcf.flow(a); 108 v = mcf.potential(n); 109 mcf.flowMap(flow); 110 mcf.potentialMap(pot); 115 v = const_mcf.totalCost(); 116 double x = const_mcf.template totalCost<double>(); 117 v = const_mcf.flow(a); 118 v = const_mcf.potential(n); 111 119 112 120 ignore_unused_variable_warning(fm); 113 121 ignore_unused_variable_warning(pm); … … 142 150 template < typename GR, typename LM, typename UM, 143 151 typename SM, typename FM > 144 152 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, 145 const SM& supply, const FM& flow ) 153 const SM& supply, const FM& flow, 154 ProblemType type = EQ ) 146 155 { 147 156 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 148 157 … … 156 165 sum += flow[e]; 157 166 for (InArcIt e(gr, n); e != INVALID; ++e) 158 167 sum = flow[e]; 159 if (sum != supply[n]) return false; 168 bool b = (type == EQ && sum == supply[n])  169 (type == GEQ && sum >= supply[n])  170 (type == LEQ && sum <= supply[n]); 171 if (!b) return false; 160 172 } 161 173 162 174 return true; … … 165 177 // Check the feasibility of the given potentials (dual soluiton) 166 178 // using the "Complementary Slackness" optimality condition 167 179 template < typename GR, typename LM, typename UM, 168 typename CM, typename FM, typename PM >180 typename CM, typename SM, typename FM, typename PM > 169 181 bool checkPotential( const GR& gr, const LM& lower, const UM& upper, 170 const CM& cost, const FM& flow, const PM& pi ) 182 const CM& cost, const SM& supply, const FM& flow, 183 const PM& pi ) 171 184 { 172 185 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 173 186 … … 179 192 (red_cost > 0 && flow[e] == lower[e])  180 193 (red_cost < 0 && flow[e] == upper[e]); 181 194 } 195 196 for (NodeIt n(gr); opt && n != INVALID; ++n) { 197 typename SM::Value sum = 0; 198 for (OutArcIt e(gr, n); e != INVALID; ++e) 199 sum += flow[e]; 200 for (InArcIt e(gr, n); e != INVALID; ++e) 201 sum = flow[e]; 202 opt = (sum == supply[n])  (pi[n] == 0); 203 } 204 182 205 return opt; 183 206 } 184 207 … … 190 213 const GR& gr, const LM& lower, const UM& upper, 191 214 const CM& cost, const SM& supply, 192 215 bool result, typename CM::Value total, 193 const std::string &test_id = "" ) 216 const std::string &test_id = "", 217 ProblemType type = EQ ) 194 218 { 195 219 check(mcf_result == result, "Wrong result " + test_id); 196 220 if (result) { 197 check(checkFlow(gr, lower, upper, supply, mcf.flowMap() ),221 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type), 198 222 "The flow is not feasible " + test_id); 199 223 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); 200 check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),224 check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(), 201 225 mcf.potentialMap()), 202 226 "Wrong potentials " + test_id); 203 227 } … … 226 250 // Read the test digraph 227 251 Digraph gr; 228 252 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr); 229 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr) ;253 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr); 230 254 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); 231 255 Node v, w; 232 256 … … 239 263 .nodeMap("sup1", s1) 240 264 .nodeMap("sup2", s2) 241 265 .nodeMap("sup3", s3) 266 .nodeMap("sup4", s4) 267 .nodeMap("sup5", s5) 242 268 .node("source", v) 243 269 .node("target", w) 244 270 .run(); … … 247 273 { 248 274 NetworkSimplex<Digraph> mcf(gr); 249 275 276 // Check the equality form 250 277 mcf.upperMap(u).costMap(c); 251 278 checkMcf(mcf, mcf.supplyMap(s1).run(), 252 279 gr, l1, u, c, s1, true, 5240, "#A1"); … … 267 294 gr, l1, cu, cc, s3, true, 0, "#A7"); 268 295 checkMcf(mcf, mcf.boundMaps(l2, u).run(), 269 296 gr, l2, u, cc, s3, false, 0, "#A8"); 297 298 // Check the GEQ form 299 mcf.reset().upperMap(u).costMap(c).supplyMap(s4); 300 checkMcf(mcf, mcf.run(), 301 gr, l1, u, c, s4, true, 3530, "#A9", GEQ); 302 mcf.problemType(mcf.GEQ); 303 checkMcf(mcf, mcf.lowerMap(l2).run(), 304 gr, l2, u, c, s4, true, 4540, "#A10", GEQ); 305 mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5); 306 checkMcf(mcf, mcf.run(), 307 gr, l2, u, c, s5, false, 0, "#A11", GEQ); 308 309 // Check the LEQ form 310 mcf.reset().problemType(mcf.LEQ); 311 mcf.upperMap(u).costMap(c).supplyMap(s5); 312 checkMcf(mcf, mcf.run(), 313 gr, l1, u, c, s5, true, 5080, "#A12", LEQ); 314 checkMcf(mcf, mcf.lowerMap(l2).run(), 315 gr, l2, u, c, s5, true, 5930, "#A13", LEQ); 316 mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4); 317 checkMcf(mcf, mcf.run(), 318 gr, l2, u, c, s4, false, 0, "#A14", LEQ); 270 319 } 271 320 272 321 // B. Test NetworkSimplex with each pivot rule