Ticket #219: ns-geq-cf49ce98b6a4.patch
File ns-geq-cf49ce98b6a4.patch, 20.0 KB (added by , 15 years ago) |
---|
-
doc/groups.dox
# HG changeset patch # User Peter Kovacs <kpeter@inf.elte.hu> # Date 1238030755 -3600 # Node ID cf49ce98b6a4e25e6dfb263baf306d76a837324e # Parent c7d160f73d52e2d3f1ab1236cb0b440daf5f7a85 Support <= supply constraints in NetworkSimplex (#219 and #234) The same inequality constraints are supported as by Circulation. The documentation of the modules and the classes NetworkSimplex, Circulation are also improved and extended with important notes and explanations. 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 Edmonds-Karp algorithm. … … 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) = 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 Cycle-canceling algorithms. 371 - \ref CapacityScaling Successive shortest path algorithm with optional 373 Note that \f$\sum_{u\in V} sup(u)\f$ must be zero in order to have a 374 feasible solution. 375 376 NetworkSimplex, which provides the most efficient solution method, 377 supports a more general form of the problem, in which the supply/demand 378 constraints (mass balance constraints) are inequalities 379 (just like the \ref Circulation "network circulation problem"). 380 381 A feasible solution for the minimum cost flow problem can be found 382 using \ref Circulation. 383 384 LEMON contains several algorithms for solving minimum cost flow problems. 385 - \ref CycleCanceling Cycle-Canceling algorithms. 386 - \ref CancelAndTighten The Cancel and Tighten algorithm. 387 - \ref CapacityScaling Successive Shortest %Path algorithm with optional 372 388 capacity scaling. 373 - \ref CostScaling Push- relabel and augment-relabel algorithms based on389 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on 374 390 cost scaling. 375 - \ref NetworkSimplex Primal network simplex algorithm with various391 - \ref NetworkSimplex Primal Network Simplex algorithm with various 376 392 pivot strategies. 393 394 In general NetworkSimplex is the most efficient implementation, 395 but in special cases other algorithms could be faster. 396 For example, if the capacities are very small, CapacityScaling is 397 usually the fastest algorithm (without effective scaling). 377 398 */ 378 399 379 400 /** -
lemon/circulation.h
diff --git a/lemon/circulation.h b/lemon/circulation.h
a b 24 24 25 25 ///\ingroup max_flow 26 26 ///\file 27 ///\brief Push-relabel algorithm for finding a feasible circulation.27 ///\brief Push-relabel algorithm for the network circulation problem. 28 28 /// 29 29 namespace lemon { 30 30 … … 119 119 120 120 The exact formulation of this problem is the following. 121 121 Let \f$G=(V,A)\f$ be a digraph, 122 \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$, 123 \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation 122 \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$ denote the lower and 123 upper bounds for the flow values on the arcs, for which 124 \f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$, and 125 \f$delta: V\rightarrow\mathbf{R}\f$ denotes signed lower bounds for the 126 actual supply of the nodes. Find a feasible circulation 124 127 \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that 125 \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)126 \geq delta(v) \quad \forall v\in V, \f]127 \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]128 \note \f$delta(v)\f$ specifies a lower bound for the supply of node129 \f$v\f$. It can be either positive or negative, however note that130 \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to131 have a feasible solution.132 128 133 \note A special case of this problem is when 134 \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$ 135 will be \e equal \e to \f$delta(v)\f$, if a circulation can be found. 136 Thus a feasible solution for the 137 \ref min_cost_flow "minimum cost flow" problem can be calculated 138 in this way. 129 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) 130 \geq delta(u) \quad \forall u\in V, \f] 131 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f] 132 133 Note that \f$\sum_{u\in V} delta(u)\f$ must be zero or negative in 134 order to have a feasible solution (since the sum of the expressions 135 on the left side of the inequalities are zero). 136 137 A special case of this problem is when the sum of the supply 138 values (\f$\sum_{u\in V} delta(u)\f$) is zero. 139 Then all the inequalities will be equality, i.e. the supply of each 140 node \f$u\f$ will be \e equal \e to \f$delta(u)\f$ if a feasible flow 141 can be found. 142 Thus a feasible solution for the \ref min_cost_flow 143 "minimum cost flow problem" can be calculated with this algorithm. 144 145 \note If you need the opposite inequality in the supply/demand 146 constraints (i.e. the total supply is greater than the total demand), 147 then you could easily transform the problem to this form by reversing 148 the direction of the arcs and taking the negative of the supply values 149 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). 139 150 140 151 \tparam GR The type of the digraph the algorithm runs on. 141 152 \tparam LM The type of the lower bound capacity map. The default -
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/circulation.h> 33 34 34 35 namespace lemon { 35 36 … … 42 43 /// \ref NetworkSimplex implements the primal Network Simplex algorithm 43 44 /// for finding a \ref min_cost_flow "minimum cost flow". 44 45 /// This algorithm is a specialized version of the linear programming 45 /// simplex method directly for the minimum cost flow problem .46 /// It is one of the most efficient solution methods.46 /// simplex method directly for the minimum cost flow problem, 47 /// and it is one of the most efficient solution methods. 47 48 /// 48 49 /// In general this class is the fastest implementation available 49 50 /// in LEMON for the minimum cost flow problem. 50 51 /// 52 /// In fact, this algorithm supports a more general form of the 53 /// \ref min_cost_flow "minimum cost flow problem", in which the 54 /// supply/demand constraints are inequalities, just like the 55 /// \ref Circulation "network circulation problem". 56 /// The exact formulation of this problem is the following 57 /// (using the same notations as in the minimum cost flow problem). 58 /** 59 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 60 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq 61 sup(u) \quad \forall u\in V \f] 62 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 63 */ 64 /// In this form \f$sup(u)\f$ specifies a lower bound for the actual 65 /// supply of node \f$u\f$. 66 /// It can be either positive or negative, however note that 67 /// \f$\sum_{u\in V} sup(u)\f$ must be zero or negative in order to 68 /// have a feasible solution (since the sum of the expressions on the 69 /// left side of the inequalities are zero). 70 /// 71 /// A special case of this problem is when the sum of the supply 72 /// values (\f$\sum_{u\in V} sup(u)\f$) is zero. 73 /// Then all the inequalities will be equality, i.e. the supply of each 74 /// node \f$u\f$ will be \e equal \e to \f$sup(u)\f$ if a feasible flow 75 /// can be found. 76 /// Thus this special case is exactly the same as the \ref min_cost_flow 77 /// "minimum cost flow problem". 78 /// 79 /// \note If you need the opposite inequality in the supply/demand 80 /// constraints (i.e. the total supply is greater than the total demand), 81 /// then you could easily transform the problem to this form by reversing 82 /// the direction of the arcs and taking the negative of the supply values 83 /// (e.g. using \ref ReverseDigraph and \ref NegMap adaptors). 84 /// 51 85 /// \tparam GR The digraph type the algorithm runs on. 52 86 /// \tparam V The value type used in the algorithm. 53 87 /// By default it is \c int. … … 142 176 143 177 // Parameters of the problem 144 178 ValueArcMap *_plower; 145 ValueArcMap *_p upper;179 ValueArcMap *_pcap; 146 180 ValueArcMap *_pcost; 147 181 ValueNodeMap *_psupply; 148 182 bool _pstsup; … … 584 618 /// \param graph The digraph the algorithm runs on. 585 619 NetworkSimplex(const GR& graph) : 586 620 _graph(graph), 587 _plower(NULL), _p upper(NULL), _pcost(NULL),621 _plower(NULL), _pcap(NULL), _pcost(NULL), 588 622 _psupply(NULL), _pstsup(false), 589 623 _flow_map(NULL), _potential_map(NULL), 590 624 _local_flow(false), _local_potential(false), … … 638 672 /// \return <tt>(*this)</tt> 639 673 template<typename UPPER> 640 674 NetworkSimplex& upperMap(const UPPER& map) { 641 delete _p upper;642 _p upper= new ValueArcMap(_graph);675 delete _pcap; 676 _pcap = new ValueArcMap(_graph); 643 677 for (ArcIt a(_graph); a != INVALID; ++a) { 644 (*_p upper)[a] = map[a];678 (*_pcap)[a] = map[a]; 645 679 } 646 680 return *this; 647 681 } … … 851 885 /// \return <tt>(*this)</tt> 852 886 NetworkSimplex& reset() { 853 887 delete _plower; 854 delete _p upper;888 delete _pcap; 855 889 delete _pcost; 856 890 delete _psupply; 857 891 _plower = NULL; 858 _p upper= NULL;892 _pcap = NULL; 859 893 _pcost = NULL; 860 894 _psupply = NULL; 861 895 _pstsup = false; … … 990 1024 991 1025 // Initialize node related data 992 1026 bool valid_supply = true; 993 if (!_pstsup && !_psupply) { 1027 Value total_supply = 0; 1028 if (!_psupply && !_pstsup) { 994 1029 _pstsup = true; 995 1030 _psource = _ptarget = NodeIt(_graph); 996 1031 _pstflow = 0; 997 1032 } 998 1033 if (_psupply) { 999 Value sum = 0;1000 1034 int i = 0; 1001 1035 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1002 1036 _node_id[n] = i; 1003 1037 _supply[i] = (*_psupply)[n]; 1004 sum+= _supply[i];1038 total_supply += _supply[i]; 1005 1039 } 1006 valid_supply = ( sum == 0);1040 valid_supply = (total_supply <= 0); 1007 1041 } else { 1008 1042 int i = 0; 1009 1043 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { … … 1011 1045 _supply[i] = 0; 1012 1046 } 1013 1047 _supply[_node_id[_psource]] = _pstflow; 1014 _supply[_node_id[_ptarget]] 1048 _supply[_node_id[_ptarget]] = -_pstflow; 1015 1049 } 1016 1050 if (!valid_supply) return false; 1017 1051 1052 const Value max_cap = std::numeric_limits<Value>::max(); 1053 const Value max_cost = std::numeric_limits<Value>::max() / 4; 1054 1055 // Run Circulation to check if a feasible solution exists 1056 typedef ConstMap<Arc, Value> ConstArcMap; 1057 ValueNodeMap *csup = NULL; 1058 bool local_csup = false; 1059 if (_psupply) { 1060 csup = _psupply; 1061 } else { 1062 csup = new ValueNodeMap(_graph, 0); 1063 (*csup)[_psource] = _pstflow; 1064 (*csup)[_ptarget] = -_pstflow; 1065 local_csup = true; 1066 } 1067 bool circ_result = false; 1068 if (_plower) { 1069 if (_pcap) { 1070 Circulation<GR, ValueArcMap, ValueArcMap, ValueNodeMap> 1071 circ(_graph, *_plower, *_pcap, *csup); 1072 circ_result = circ.run(); 1073 } else { 1074 Circulation<GR, ValueArcMap, ConstArcMap, ValueNodeMap> 1075 circ(_graph, *_plower, ConstArcMap(max_cap), *csup); 1076 circ_result = circ.run(); 1077 } 1078 } else { 1079 if (_pcap) { 1080 Circulation<GR, ConstArcMap, ValueArcMap, ValueNodeMap> 1081 circ(_graph, ConstArcMap(0), *_pcap, *csup); 1082 circ_result = circ.run(); 1083 } else { 1084 Circulation<GR, ConstArcMap, ConstArcMap, ValueNodeMap> 1085 circ(_graph, ConstArcMap(0), ConstArcMap(max_cap), *csup); 1086 circ_result = circ.run(); 1087 } 1088 } 1089 if (local_csup) delete csup; 1090 if (!circ_result) return false; 1091 1018 1092 // Set data for the artificial root node 1019 1093 _root = _node_num; 1020 1094 _parent[_root] = -1; … … 1023 1097 _rev_thread[0] = _root; 1024 1098 _succ_num[_root] = all_node_num; 1025 1099 _last_succ[_root] = _root - 1; 1026 _supply[_root] = 0;1100 _supply[_root] = -total_supply; 1027 1101 _pi[_root] = 0; 1028 1102 1029 1103 // Store the arcs in a mixed order … … 1035 1109 } 1036 1110 1037 1111 // Initialize arc maps 1038 if (_p upper&& _pcost) {1112 if (_pcap && _pcost) { 1039 1113 for (int i = 0; i != _arc_num; ++i) { 1040 1114 Arc e = _arc_ref[i]; 1041 1115 _source[i] = _node_id[_graph.source(e)]; 1042 1116 _target[i] = _node_id[_graph.target(e)]; 1043 _cap[i] = (*_p upper)[e];1117 _cap[i] = (*_pcap)[e]; 1044 1118 _cost[i] = (*_pcost)[e]; 1045 1119 _flow[i] = 0; 1046 1120 _state[i] = STATE_LOWER; … … 1053 1127 _flow[i] = 0; 1054 1128 _state[i] = STATE_LOWER; 1055 1129 } 1056 if (_p upper) {1130 if (_pcap) { 1057 1131 for (int i = 0; i != _arc_num; ++i) 1058 _cap[i] = (*_p upper)[_arc_ref[i]];1132 _cap[i] = (*_pcap)[_arc_ref[i]]; 1059 1133 } else { 1060 Value val = std::numeric_limits<Value>::max();1134 Value val = max_cap; 1061 1135 for (int i = 0; i != _arc_num; ++i) 1062 1136 _cap[i] = val; 1063 1137 } … … 1083 1157 } 1084 1158 1085 1159 // Add artificial arcs and initialize the spanning tree data structure 1086 Value max_cap = std::numeric_limits<Value>::max();1087 Value max_cost = std::numeric_limits<Value>::max() / 4;1088 1160 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1089 1161 _thread[u] = u + 1; 1090 1162 _rev_thread[u + 1] = u; … … 1095 1167 _cost[e] = max_cost; 1096 1168 _cap[e] = max_cap; 1097 1169 _state[e] = STATE_TREE; 1098 if (_supply[u] > = 0) {1170 if (_supply[u] > 0 || (_supply[u] == 0 && total_supply == 0)) { 1099 1171 _flow[e] = _supply[u]; 1100 1172 _forward[u] = true; 1101 1173 _pi[u] = -max_cost; … … 1371 1443 } 1372 1444 } 1373 1445 1374 // Check if the flow amount equals zero on all the artificial arcs1375 for (int e = _arc_num; e != _arc_num + _node_num; ++e) {1376 if (_flow[e] > 0) return false;1377 }1378 1379 1446 // Copy flow values to _flow_map 1380 1447 if (_plower) { 1381 1448 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\n" 37 " 1 20 27 0 20\n" 38 " 2 -4 0 0 -8\n" 39 " 3 0 0 0 0\n" 40 " 4 0 0 0 0\n" 41 " 5 9 0 0 6\n" 42 " 6 -6 0 0 -5\n" 43 " 7 0 0 0 0\n" 44 " 8 0 0 0 0\n" 45 " 9 3 0 0 0\n" 46 " 10 -2 0 0 -7\n" 47 " 11 0 0 0 -10\n" 48 " 12 -20 -27 0 -30\n" 49 49 "\n" 50 50 "@arcs\n" 51 51 " cost cap low1 low2\n" … … 141 141 template < typename GR, typename LM, typename UM, 142 142 typename SM, typename FM > 143 143 bool checkFlow( const GR& gr, const LM& lower, const UM& upper, 144 const SM& supply, const FM& flow )144 const SM& supply, const FM& flow, bool geq = false ) 145 145 { 146 146 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 147 147 … … 155 155 sum += flow[e]; 156 156 for (InArcIt e(gr, n); e != INVALID; ++e) 157 157 sum -= flow[e]; 158 if (sum != supply[n]) return false; 158 if ( (!geq && sum != supply[n]) || (geq && sum < supply[n]) ) 159 return false; 159 160 } 160 161 161 162 return true; … … 189 190 const GR& gr, const LM& lower, const UM& upper, 190 191 const CM& cost, const SM& supply, 191 192 bool result, typename CM::Value total, 192 const std::string &test_id = "" )193 const std::string &test_id = "", bool geq = false ) 193 194 { 194 195 check(mcf_result == result, "Wrong result " + test_id); 195 196 if (result) { 196 check(checkFlow(gr, lower, upper, supply, mcf.flowMap() ),197 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), geq), 197 198 "The flow is not feasible " + test_id); 198 199 check(mcf.totalCost() == total, "The flow is not optimal " + test_id); 199 200 check(checkPotential(gr, lower, upper, cost, mcf.flowMap(), … … 224 225 // Read the test digraph 225 226 Digraph gr; 226 227 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr); 227 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr) ;228 Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr); 228 229 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); 229 230 Node v, w; 230 231 … … 237 238 .nodeMap("sup1", s1) 238 239 .nodeMap("sup2", s2) 239 240 .nodeMap("sup3", s3) 241 .nodeMap("sup4", s4) 240 242 .node("source", v) 241 243 .node("target", w) 242 244 .run(); … … 245 247 { 246 248 NetworkSimplex<Digraph> mcf(gr); 247 249 250 // Check the equality form 248 251 mcf.upperMap(u).costMap(c); 249 252 checkMcf(mcf, mcf.supplyMap(s1).run(), 250 253 gr, l1, u, c, s1, true, 5240, "#A1"); … … 265 268 gr, l1, cu, cc, s3, true, 0, "#A7"); 266 269 checkMcf(mcf, mcf.boundMaps(l2, u).run(), 267 270 gr, l2, u, cc, s3, false, 0, "#A8"); 271 272 // Check the inequality form 273 mcf.reset().upperMap(u).costMap(c).supplyMap(s4); 274 checkMcf(mcf, mcf.run(), 275 gr, l1, u, c, s4, true, 3530, "#A9", true); 276 checkMcf(mcf, mcf.lowerMap(l2).run(), 277 gr, l2, u, c, s4, true, 4540, "#A10", true); 268 278 } 269 279 270 280 // B. Test NetworkSimplex with each pivot rule