Changeset 609:e6927fe719e6 in lemon1.2
 Timestamp:
 04/17/09 18:04:36 (10 years ago)
 Branch:
 default
 Phase:
 public
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

doc/groups.dox
r455 r609 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: … … 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. 363 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 supply(v) \qquad \forall v\in V \f] 367 \f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f] 368 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 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. 367 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] 372 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 
lemon/network_simplex.h
r608 r609 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 { … … 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. … … 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 … … 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: … … 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: … … 156 221 Node _psource, _ptarget; 157 222 Flow _pstflow; 223 ProblemType _ptype; 158 224 159 225 // Result maps … … 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. … … 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), … … 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 /// … … 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. … … 796 882 return *this; 797 883 } 884 885 /// @} 798 886 799 887 /// \name Execution Control … … 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); … … 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 … … 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 } … … 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; … … 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]; 1015 } 1016 valid_supply = (sum == 0); 1112 sum_supply += _supply[i]; 1113 } 1114 valid_supply = (_ptype == GEQ && sum_supply <= 0)  1115 (_ptype == LEQ && sum_supply >= 0); 1017 1116 } else { 1018 1117 int 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; 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; 1027 1211 1028 1212 // Set data for the artificial root node … … 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 … … 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) { … … 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) { … … 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 } … … 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) { 
test/min_cost_flow_test.cc
r607 r609 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" … … 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> … … 98 104 .supplyMap(sup) 99 105 .stSupply(n, n, k) 106 .flowMap(flow) 107 .potentialMap(pot) 100 108 .run(); 101 102 const typename MCF::FlowMap &fm = mcf.flowMap();103 const typename MCF::PotentialMap &pm = mcf.potentialMap(); 104 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);109 110 const MCF& const_mcf = mcf; 111 112 const typename MCF::FlowMap &fm = const_mcf.flowMap(); 113 const typename MCF::PotentialMap &pm = const_mcf.potentialMap(); 114 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); … … 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); … … 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 … … 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); … … 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 } … … 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); … … 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; … … 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) … … 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(), … … 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
Note: See TracChangeset
for help on using the changeset viewer.