Changeset 872:fa6f37d7a25b in lemon for lemon
 Timestamp:
 11/12/09 23:26:13 (10 years ago)
 Branch:
 default
 Phase:
 public
 Rebase:
 39653339343431333337306539356364643765343037643134313235643439386231623630306163
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/capacity_scaling.h
r871 r872 20 20 #define LEMON_CAPACITY_SCALING_H 21 21 22 /// \ingroup min_cost_flow 22 /// \ingroup min_cost_flow_algs 23 23 /// 24 24 /// \file 25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.25 /// \brief Capacity Scaling algorithm for finding a minimum cost flow. 26 26 27 27 #include <vector> 28 #include <limits> 29 #include <lemon/core.h> 28 30 #include <lemon/bin_heap.h> 29 31 30 32 namespace lemon { 31 33 32 /// \addtogroup min_cost_flow 34 /// \addtogroup min_cost_flow_algs 33 35 /// @{ 34 36 35 /// \brief Implementation of the capacity scaling algorithm for36 /// finding a minimum cost flow.37 /// \brief Implementation of the Capacity Scaling algorithm for 38 /// finding a \ref min_cost_flow "minimum cost flow". 37 39 /// 38 40 /// \ref CapacityScaling implements the capacity scaling version 39 /// of the successive shortest path algorithm for finding a minimum 40 /// cost flow. 41 /// of the successive shortest path algorithm for finding a 42 /// \ref min_cost_flow "minimum cost flow". It is an efficient dual 43 /// solution method. 41 44 /// 42 /// \tparam Digraph The digraph type the algorithm runs on. 43 /// \tparam LowerMap The type of the lower bound map. 44 /// \tparam CapacityMap The type of the capacity (upper bound) map. 45 /// \tparam CostMap The type of the cost (length) map. 46 /// \tparam SupplyMap The type of the supply map. 45 /// Most of the parameters of the problem (except for the digraph) 46 /// can be given using separate functions, and the algorithm can be 47 /// executed using the \ref run() function. If some parameters are not 48 /// specified, then default values will be used. 47 49 /// 48 /// \ warning49 ///  Arc capacities and costs should be \e nonnegative \e integers.50 ///  Supply values should be \e signed \e integers.51 ///  The value types of the maps should be convertible to each other.52 ///  \c CostMap::Value must be signed type.50 /// \tparam GR The digraph type the algorithm runs on. 51 /// \tparam V The value type used for flow amounts, capacity bounds 52 /// and supply values in the algorithm. By default it is \c int. 53 /// \tparam C The value type used for costs and potentials in the 54 /// algorithm. By default it is the same as \c V. 53 55 /// 54 /// \author Peter Kovacs 55 template < typename Digraph, 56 typename LowerMap = typename Digraph::template ArcMap<int>, 57 typename CapacityMap = typename Digraph::template ArcMap<int>, 58 typename CostMap = typename Digraph::template ArcMap<int>, 59 typename SupplyMap = typename Digraph::template NodeMap<int> > 56 /// \warning Both value types must be signed and all input data must 57 /// be integer. 58 /// \warning This algorithm does not support negative costs for such 59 /// arcs that have infinite upper bound. 60 template <typename GR, typename V = int, typename C = V> 60 61 class CapacityScaling 61 62 { 62 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);63 64 typedef typename CapacityMap::Value Capacity;65 typedef typename CostMap::Value Cost;66 typedef typename SupplyMap::Value Supply;67 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;68 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;69 typedef typename Digraph::template NodeMap<Arc> PredMap;70 71 63 public: 72 64 73 /// The type of the flow map. 74 typedef typename Digraph::template ArcMap<Capacity> FlowMap; 75 /// The type of the potential map. 76 typedef typename Digraph::template NodeMap<Cost> PotentialMap; 77 65 /// The type of the flow amounts, capacity bounds and supply values 66 typedef V Value; 67 /// The type of the arc costs 68 typedef C Cost; 69 70 public: 71 72 /// \brief Problem type constants for the \c run() function. 73 /// 74 /// Enum type containing the problem type constants that can be 75 /// returned by the \ref run() function of the algorithm. 76 enum ProblemType { 77 /// The problem has no feasible solution (flow). 78 INFEASIBLE, 79 /// The problem has optimal solution (i.e. it is feasible and 80 /// bounded), and the algorithm has found optimal flow and node 81 /// potentials (primal and dual solutions). 82 OPTIMAL, 83 /// The digraph contains an arc of negative cost and infinite 84 /// upper bound. It means that the objective function is unbounded 85 /// on that arc, however note that it could actually be bounded 86 /// over the feasible flows, but this algroithm cannot handle 87 /// these cases. 88 UNBOUNDED 89 }; 90 78 91 private: 79 92 80 /// \brief Special implementation of the \ref Dijkstra algorithm 81 /// for finding shortest paths in the residual network. 82 /// 83 /// \ref ResidualDijkstra is a special implementation of the 84 /// \ref Dijkstra algorithm for finding shortest paths in the 85 /// residual network of the digraph with respect to the reduced arc 86 /// costs and modifying the node potentials according to the 87 /// distance of the nodes. 93 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 94 95 typedef std::vector<Arc> ArcVector; 96 typedef std::vector<Node> NodeVector; 97 typedef std::vector<int> IntVector; 98 typedef std::vector<bool> BoolVector; 99 typedef std::vector<Value> ValueVector; 100 typedef std::vector<Cost> CostVector; 101 102 private: 103 104 // Data related to the underlying digraph 105 const GR &_graph; 106 int _node_num; 107 int _arc_num; 108 int _res_arc_num; 109 int _root; 110 111 // Parameters of the problem 112 bool _have_lower; 113 Value _sum_supply; 114 115 // Data structures for storing the digraph 116 IntNodeMap _node_id; 117 IntArcMap _arc_idf; 118 IntArcMap _arc_idb; 119 IntVector _first_out; 120 BoolVector _forward; 121 IntVector _source; 122 IntVector _target; 123 IntVector _reverse; 124 125 // Node and arc data 126 ValueVector _lower; 127 ValueVector _upper; 128 CostVector _cost; 129 ValueVector _supply; 130 131 ValueVector _res_cap; 132 CostVector _pi; 133 ValueVector _excess; 134 IntVector _excess_nodes; 135 IntVector _deficit_nodes; 136 137 Value _delta; 138 int _phase_num; 139 IntVector _pred; 140 141 public: 142 143 /// \brief Constant for infinite upper bounds (capacities). 144 /// 145 /// Constant for infinite upper bounds (capacities). 146 /// It is \c std::numeric_limits<Value>::infinity() if available, 147 /// \c std::numeric_limits<Value>::max() otherwise. 148 const Value INF; 149 150 private: 151 152 // Special implementation of the Dijkstra algorithm for finding 153 // shortest paths in the residual network of the digraph with 154 // respect to the reduced arc costs and modifying the node 155 // potentials according to the found distance labels. 88 156 class ResidualDijkstra 89 157 { 90 typedef typename Digraph::template NodeMap<int> HeapCrossRef;158 typedef RangeMap<int> HeapCrossRef; 91 159 typedef BinHeap<Cost, HeapCrossRef> Heap; 92 160 93 161 private: 94 162 95 // The digraph the algorithm runs on 96 const Digraph &_graph; 97 98 // The main maps 99 const FlowMap &_flow; 100 const CapacityArcMap &_res_cap; 101 const CostMap &_cost; 102 const SupplyNodeMap &_excess; 103 PotentialMap &_potential; 104 105 // The distance map 106 PotentialMap _dist; 107 // The pred arc map 108 PredMap &_pred; 109 // The processed (i.e. permanently labeled) nodes 110 std::vector<Node> _proc_nodes; 111 163 int _node_num; 164 const IntVector &_first_out; 165 const IntVector &_target; 166 const CostVector &_cost; 167 const ValueVector &_res_cap; 168 const ValueVector &_excess; 169 CostVector &_pi; 170 IntVector &_pred; 171 172 IntVector _proc_nodes; 173 CostVector _dist; 174 112 175 public: 113 176 114 /// Constructor. 115 ResidualDijkstra( const Digraph &digraph, 116 const FlowMap &flow, 117 const CapacityArcMap &res_cap, 118 const CostMap &cost, 119 const SupplyMap &excess, 120 PotentialMap &potential, 121 PredMap &pred ) : 122 _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost), 123 _excess(excess), _potential(potential), _dist(digraph), 124 _pred(pred) 177 ResidualDijkstra(CapacityScaling& cs) : 178 _node_num(cs._node_num), _first_out(cs._first_out), 179 _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap), 180 _excess(cs._excess), _pi(cs._pi), _pred(cs._pred), 181 _dist(cs._node_num) 125 182 {} 126 183 127 /// Run the algorithm from the given source node. 128 Node run(Node s, Capacity delta = 1) { 129 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); 184 int run(int s, Value delta = 1) { 185 HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP); 130 186 Heap heap(heap_cross_ref); 131 187 heap.push(s, 0); 132 _pred[s] = INVALID;188 _pred[s] = 1; 133 189 _proc_nodes.clear(); 134 190 135 // Process ingnodes191 // Process nodes 136 192 while (!heap.empty() && _excess[heap.top()] > delta) { 137 Nodeu = heap.top(), v;138 Cost d = heap.prio() + _p otential[u], nd;193 int u = heap.top(), v; 194 Cost d = heap.prio() + _pi[u], dn; 139 195 _dist[u] = heap.prio(); 196 _proc_nodes.push_back(u); 140 197 heap.pop(); 141 _proc_nodes.push_back(u); 142 143 // Traversing outgoing arcs 144 for (OutArcIt e(_graph, u); e != INVALID; ++e) { 145 if (_res_cap[e] >= delta) { 146 v = _graph.target(e); 147 switch(heap.state(v)) { 198 199 // Traverse outgoing residual arcs 200 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { 201 if (_res_cap[a] < delta) continue; 202 v = _target[a]; 203 switch (heap.state(v)) { 148 204 case Heap::PRE_HEAP: 149 heap.push(v, d + _cost[ e]  _potential[v]);150 _pred[v] = e;205 heap.push(v, d + _cost[a]  _pi[v]); 206 _pred[v] = a; 151 207 break; 152 208 case Heap::IN_HEAP: 153 nd = d + _cost[e]  _potential[v];154 if ( nd< heap[v]) {155 heap.decrease(v, nd);156 _pred[v] = e;209 dn = d + _cost[a]  _pi[v]; 210 if (dn < heap[v]) { 211 heap.decrease(v, dn); 212 _pred[v] = a; 157 213 } 158 214 break; 159 215 case Heap::POST_HEAP: 160 216 break; 161 }162 217 } 163 218 } 164 165 // Traversing incoming arcs 166 for (InArcIt e(_graph, u); e != INVALID; ++e) { 167 if (_flow[e] >= delta) { 168 v = _graph.source(e); 169 switch(heap.state(v)) { 170 case Heap::PRE_HEAP: 171 heap.push(v, d  _cost[e]  _potential[v]); 172 _pred[v] = e; 173 break; 174 case Heap::IN_HEAP: 175 nd = d  _cost[e]  _potential[v]; 176 if (nd < heap[v]) { 177 heap.decrease(v, nd); 178 _pred[v] = e; 179 } 180 break; 181 case Heap::POST_HEAP: 182 break; 183 } 184 } 185 } 186 } 187 if (heap.empty()) return INVALID; 188 189 // Updating potentials of processed nodes 190 Node t = heap.top(); 191 Cost t_dist = heap.prio(); 192 for (int i = 0; i < int(_proc_nodes.size()); ++i) 193 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]]  t_dist; 219 } 220 if (heap.empty()) return 1; 221 222 // Update potentials of processed nodes 223 int t = heap.top(); 224 Cost dt = heap.prio(); 225 for (int i = 0; i < int(_proc_nodes.size()); ++i) { 226 _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]]  dt; 227 } 194 228 195 229 return t; … … 198 232 }; //class ResidualDijkstra 199 233 200 private:201 202 // The digraph the algorithm runs on203 const Digraph &_graph;204 // The original lower bound map205 const LowerMap *_lower;206 // The modified capacity map207 CapacityArcMap _capacity;208 // The original cost map209 const CostMap &_cost;210 // The modified supply map211 SupplyNodeMap _supply;212 bool _valid_supply;213 214 // Arc map of the current flow215 FlowMap *_flow;216 bool _local_flow;217 // Node map of the current potentials218 PotentialMap *_potential;219 bool _local_potential;220 221 // The residual capacity map222 CapacityArcMap _res_cap;223 // The excess map224 SupplyNodeMap _excess;225 // The excess nodes (i.e. nodes with positive excess)226 std::vector<Node> _excess_nodes;227 // The deficit nodes (i.e. nodes with negative excess)228 std::vector<Node> _deficit_nodes;229 230 // The delta parameter used for capacity scaling231 Capacity _delta;232 // The maximum number of phases233 int _phase_num;234 235 // The pred arc map236 PredMap _pred;237 // Implementation of the Dijkstra algorithm for finding augmenting238 // shortest paths in the residual network239 ResidualDijkstra *_dijkstra;240 241 234 public: 242 235 243 /// \brief General constructor (with lower bounds). 244 /// 245 /// General constructor (with lower bounds). 246 /// 247 /// \param digraph The digraph the algorithm runs on. 248 /// \param lower The lower bounds of the arcs. 249 /// \param capacity The capacities (upper bounds) of the arcs. 250 /// \param cost The cost (length) values of the arcs. 251 /// \param supply The supply values of the nodes (signed). 252 CapacityScaling( const Digraph &digraph, 253 const LowerMap &lower, 254 const CapacityMap &capacity, 255 const CostMap &cost, 256 const SupplyMap &supply ) : 257 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost), 258 _supply(digraph), _flow(NULL), _local_flow(false), 259 _potential(NULL), _local_potential(false), 260 _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL) 236 /// \brief Constructor. 237 /// 238 /// The constructor of the class. 239 /// 240 /// \param graph The digraph the algorithm runs on. 241 CapacityScaling(const GR& graph) : 242 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 243 INF(std::numeric_limits<Value>::has_infinity ? 244 std::numeric_limits<Value>::infinity() : 245 std::numeric_limits<Value>::max()) 261 246 { 262 Supply sum = 0; 247 // Check the value types 248 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 249 "The flow type of CapacityScaling must be signed"); 250 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 251 "The cost type of CapacityScaling must be signed"); 252 253 // Resize vectors 254 _node_num = countNodes(_graph); 255 _arc_num = countArcs(_graph); 256 _res_arc_num = 2 * (_arc_num + _node_num); 257 _root = _node_num; 258 ++_node_num; 259 260 _first_out.resize(_node_num + 1); 261 _forward.resize(_res_arc_num); 262 _source.resize(_res_arc_num); 263 _target.resize(_res_arc_num); 264 _reverse.resize(_res_arc_num); 265 266 _lower.resize(_res_arc_num); 267 _upper.resize(_res_arc_num); 268 _cost.resize(_res_arc_num); 269 _supply.resize(_node_num); 270 271 _res_cap.resize(_res_arc_num); 272 _pi.resize(_node_num); 273 _excess.resize(_node_num); 274 _pred.resize(_node_num); 275 276 // Copy the graph 277 int i = 0, j = 0, k = 2 * _arc_num + _node_num  1; 278 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 279 _node_id[n] = i; 280 } 281 i = 0; 282 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 283 _first_out[i] = j; 284 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { 285 _arc_idf[a] = j; 286 _forward[j] = true; 287 _source[j] = i; 288 _target[j] = _node_id[_graph.runningNode(a)]; 289 } 290 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { 291 _arc_idb[a] = j; 292 _forward[j] = false; 293 _source[j] = i; 294 _target[j] = _node_id[_graph.runningNode(a)]; 295 } 296 _forward[j] = false; 297 _source[j] = i; 298 _target[j] = _root; 299 _reverse[j] = k; 300 _forward[k] = true; 301 _source[k] = _root; 302 _target[k] = i; 303 _reverse[k] = j; 304 ++j; ++k; 305 } 306 _first_out[i] = j; 307 _first_out[_node_num] = k; 308 for (ArcIt a(_graph); a != INVALID; ++a) { 309 int fi = _arc_idf[a]; 310 int bi = _arc_idb[a]; 311 _reverse[fi] = bi; 312 _reverse[bi] = fi; 313 } 314 315 // Reset parameters 316 reset(); 317 } 318 319 /// \name Parameters 320 /// The parameters of the algorithm can be specified using these 321 /// functions. 322 323 /// @{ 324 325 /// \brief Set the lower bounds on the arcs. 326 /// 327 /// This function sets the lower bounds on the arcs. 328 /// If it is not used before calling \ref run(), the lower bounds 329 /// will be set to zero on all arcs. 330 /// 331 /// \param map An arc map storing the lower bounds. 332 /// Its \c Value type must be convertible to the \c Value type 333 /// of the algorithm. 334 /// 335 /// \return <tt>(*this)</tt> 336 template <typename LowerMap> 337 CapacityScaling& lowerMap(const LowerMap& map) { 338 _have_lower = true; 339 for (ArcIt a(_graph); a != INVALID; ++a) { 340 _lower[_arc_idf[a]] = map[a]; 341 _lower[_arc_idb[a]] = map[a]; 342 } 343 return *this; 344 } 345 346 /// \brief Set the upper bounds (capacities) on the arcs. 347 /// 348 /// This function sets the upper bounds (capacities) on the arcs. 349 /// If it is not used before calling \ref run(), the upper bounds 350 /// will be set to \ref INF on all arcs (i.e. the flow value will be 351 /// unbounded from above on each arc). 352 /// 353 /// \param map An arc map storing the upper bounds. 354 /// Its \c Value type must be convertible to the \c Value type 355 /// of the algorithm. 356 /// 357 /// \return <tt>(*this)</tt> 358 template<typename UpperMap> 359 CapacityScaling& upperMap(const UpperMap& map) { 360 for (ArcIt a(_graph); a != INVALID; ++a) { 361 _upper[_arc_idf[a]] = map[a]; 362 } 363 return *this; 364 } 365 366 /// \brief Set the costs of the arcs. 367 /// 368 /// This function sets the costs of the arcs. 369 /// If it is not used before calling \ref run(), the costs 370 /// will be set to \c 1 on all arcs. 371 /// 372 /// \param map An arc map storing the costs. 373 /// Its \c Value type must be convertible to the \c Cost type 374 /// of the algorithm. 375 /// 376 /// \return <tt>(*this)</tt> 377 template<typename CostMap> 378 CapacityScaling& costMap(const CostMap& map) { 379 for (ArcIt a(_graph); a != INVALID; ++a) { 380 _cost[_arc_idf[a]] = map[a]; 381 _cost[_arc_idb[a]] = map[a]; 382 } 383 return *this; 384 } 385 386 /// \brief Set the supply values of the nodes. 387 /// 388 /// This function sets the supply values of the nodes. 389 /// If neither this function nor \ref stSupply() is used before 390 /// calling \ref run(), the supply of each node will be set to zero. 391 /// 392 /// \param map A node map storing the supply values. 393 /// Its \c Value type must be convertible to the \c Value type 394 /// of the algorithm. 395 /// 396 /// \return <tt>(*this)</tt> 397 template<typename SupplyMap> 398 CapacityScaling& supplyMap(const SupplyMap& map) { 263 399 for (NodeIt n(_graph); n != INVALID; ++n) { 264 _supply[n] = supply[n]; 265 _excess[n] = supply[n]; 266 sum += supply[n]; 267 } 268 _valid_supply = sum == 0; 269 for (ArcIt a(_graph); a != INVALID; ++a) { 270 _capacity[a] = capacity[a]; 271 _res_cap[a] = capacity[a]; 272 } 273 274 // Remove nonzero lower bounds 275 typename LowerMap::Value lcap; 276 for (ArcIt e(_graph); e != INVALID; ++e) { 277 if ((lcap = lower[e]) != 0) { 278 _capacity[e] = lcap; 279 _res_cap[e] = lcap; 280 _supply[_graph.source(e)] = lcap; 281 _supply[_graph.target(e)] += lcap; 282 _excess[_graph.source(e)] = lcap; 283 _excess[_graph.target(e)] += lcap; 284 } 285 } 286 } 287 /* 288 /// \brief General constructor (without lower bounds). 289 /// 290 /// General constructor (without lower bounds). 291 /// 292 /// \param digraph The digraph the algorithm runs on. 293 /// \param capacity The capacities (upper bounds) of the arcs. 294 /// \param cost The cost (length) values of the arcs. 295 /// \param supply The supply values of the nodes (signed). 296 CapacityScaling( const Digraph &digraph, 297 const CapacityMap &capacity, 298 const CostMap &cost, 299 const SupplyMap &supply ) : 300 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), 301 _supply(supply), _flow(NULL), _local_flow(false), 302 _potential(NULL), _local_potential(false), 303 _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL) 304 { 305 // Check the sum of supply values 306 Supply sum = 0; 307 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; 308 _valid_supply = sum == 0; 309 } 310 311 /// \brief Simple constructor (with lower bounds). 312 /// 313 /// Simple constructor (with lower bounds). 314 /// 315 /// \param digraph The digraph the algorithm runs on. 316 /// \param lower The lower bounds of the arcs. 317 /// \param capacity The capacities (upper bounds) of the arcs. 318 /// \param cost The cost (length) values of the arcs. 400 _supply[_node_id[n]] = map[n]; 401 } 402 return *this; 403 } 404 405 /// \brief Set single source and target nodes and a supply value. 406 /// 407 /// This function sets a single source node and a single target node 408 /// and the required flow value. 409 /// If neither this function nor \ref supplyMap() is used before 410 /// calling \ref run(), the supply of each node will be set to zero. 411 /// 412 /// Using this function has the same effect as using \ref supplyMap() 413 /// with such a map in which \c k is assigned to \c s, \c k is 414 /// assigned to \c t and all other nodes have zero supply value. 415 /// 319 416 /// \param s The source node. 320 417 /// \param t The target node. 321 /// \param flow_value The required amount of flow from node \c s 322 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 323 CapacityScaling( const Digraph &digraph, 324 const LowerMap &lower, 325 const CapacityMap &capacity, 326 const CostMap &cost, 327 Node s, Node t, 328 Supply flow_value ) : 329 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost), 330 _supply(digraph, 0), _flow(NULL), _local_flow(false), 331 _potential(NULL), _local_potential(false), 332 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) 333 { 334 // Remove nonzero lower bounds 335 _supply[s] = _excess[s] = flow_value; 336 _supply[t] = _excess[t] = flow_value; 337 typename LowerMap::Value lcap; 338 for (ArcIt e(_graph); e != INVALID; ++e) { 339 if ((lcap = lower[e]) != 0) { 340 _capacity[e] = lcap; 341 _res_cap[e] = lcap; 342 _supply[_graph.source(e)] = lcap; 343 _supply[_graph.target(e)] += lcap; 344 _excess[_graph.source(e)] = lcap; 345 _excess[_graph.target(e)] += lcap; 346 } 347 } 348 _valid_supply = true; 349 } 350 351 /// \brief Simple constructor (without lower bounds). 352 /// 353 /// Simple constructor (without lower bounds). 354 /// 355 /// \param digraph The digraph the algorithm runs on. 356 /// \param capacity The capacities (upper bounds) of the arcs. 357 /// \param cost The cost (length) values of the arcs. 358 /// \param s The source node. 359 /// \param t The target node. 360 /// \param flow_value The required amount of flow from node \c s 361 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 362 CapacityScaling( const Digraph &digraph, 363 const CapacityMap &capacity, 364 const CostMap &cost, 365 Node s, Node t, 366 Supply flow_value ) : 367 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), 368 _supply(digraph, 0), _flow(NULL), _local_flow(false), 369 _potential(NULL), _local_potential(false), 370 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) 371 { 372 _supply[s] = _excess[s] = flow_value; 373 _supply[t] = _excess[t] = flow_value; 374 _valid_supply = true; 375 } 376 */ 377 /// Destructor. 378 ~CapacityScaling() { 379 if (_local_flow) delete _flow; 380 if (_local_potential) delete _potential; 381 delete _dijkstra; 382 } 383 384 /// \brief Set the flow map. 385 /// 386 /// Set the flow map. 387 /// 388 /// \return \c (*this) 389 CapacityScaling& flowMap(FlowMap &map) { 390 if (_local_flow) { 391 delete _flow; 392 _local_flow = false; 393 } 394 _flow = ↦ 418 /// \param k The required amount of flow from node \c s to node \c t 419 /// (i.e. the supply of \c s and the demand of \c t). 420 /// 421 /// \return <tt>(*this)</tt> 422 CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { 423 for (int i = 0; i != _node_num; ++i) { 424 _supply[i] = 0; 425 } 426 _supply[_node_id[s]] = k; 427 _supply[_node_id[t]] = k; 395 428 return *this; 396 429 } 397 398 /// \brief Set the potential map. 399 /// 400 /// Set the potential map. 401 /// 402 /// \return \c (*this) 403 CapacityScaling& potentialMap(PotentialMap &map) { 404 if (_local_potential) { 405 delete _potential; 406 _local_potential = false; 407 } 408 _potential = ↦ 409 return *this; 410 } 430 431 /// @} 411 432 412 433 /// \name Execution control … … 417 438 /// 418 439 /// This function runs the algorithm. 440 /// The paramters can be specified using functions \ref lowerMap(), 441 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 442 /// For example, 443 /// \code 444 /// CapacityScaling<ListDigraph> cs(graph); 445 /// cs.lowerMap(lower).upperMap(upper).costMap(cost) 446 /// .supplyMap(sup).run(); 447 /// \endcode 448 /// 449 /// This function can be called more than once. All the parameters 450 /// that have been given are kept for the next call, unless 451 /// \ref reset() is called, thus only the modified parameters 452 /// have to be set again. See \ref reset() for examples. 453 /// However the underlying digraph must not be modified after this 454 /// class have been constructed, since it copies the digraph. 419 455 /// 420 456 /// \param scaling Enable or disable capacity scaling. 421 /// If the maximum arc capacityand/or the amount of total supply457 /// If the maximum upper bound and/or the amount of total supply 422 458 /// is rather small, the algorithm could be slightly faster without 423 459 /// scaling. 424 460 /// 425 /// \return \c true if a feasible flow can be found. 426 bool run(bool scaling = true) { 427 return init(scaling) && start(); 461 /// \return \c INFEASIBLE if no feasible flow exists, 462 /// \n \c OPTIMAL if the problem has optimal solution 463 /// (i.e. it is feasible and bounded), and the algorithm has found 464 /// optimal flow and node potentials (primal and dual solutions), 465 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost 466 /// and infinite upper bound. It means that the objective function 467 /// is unbounded on that arc, however note that it could actually be 468 /// bounded over the feasible flows, but this algroithm cannot handle 469 /// these cases. 470 /// 471 /// \see ProblemType 472 ProblemType run(bool scaling = true) { 473 ProblemType pt = init(scaling); 474 if (pt != OPTIMAL) return pt; 475 return start(); 476 } 477 478 /// \brief Reset all the parameters that have been given before. 479 /// 480 /// This function resets all the paramaters that have been given 481 /// before using functions \ref lowerMap(), \ref upperMap(), 482 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). 483 /// 484 /// It is useful for multiple run() calls. If this function is not 485 /// used, all the parameters given before are kept for the next 486 /// \ref run() call. 487 /// However the underlying digraph must not be modified after this 488 /// class have been constructed, since it copies and extends the graph. 489 /// 490 /// For example, 491 /// \code 492 /// CapacityScaling<ListDigraph> cs(graph); 493 /// 494 /// // First run 495 /// cs.lowerMap(lower).upperMap(upper).costMap(cost) 496 /// .supplyMap(sup).run(); 497 /// 498 /// // Run again with modified cost map (reset() is not called, 499 /// // so only the cost map have to be set again) 500 /// cost[e] += 100; 501 /// cs.costMap(cost).run(); 502 /// 503 /// // Run again from scratch using reset() 504 /// // (the lower bounds will be set to zero on all arcs) 505 /// cs.reset(); 506 /// cs.upperMap(capacity).costMap(cost) 507 /// .supplyMap(sup).run(); 508 /// \endcode 509 /// 510 /// \return <tt>(*this)</tt> 511 CapacityScaling& reset() { 512 for (int i = 0; i != _node_num; ++i) { 513 _supply[i] = 0; 514 } 515 for (int j = 0; j != _res_arc_num; ++j) { 516 _lower[j] = 0; 517 _upper[j] = INF; 518 _cost[j] = _forward[j] ? 1 : 1; 519 } 520 _have_lower = false; 521 return *this; 428 522 } 429 523 … … 433 527 /// The results of the algorithm can be obtained using these 434 528 /// functions.\n 435 /// \ref lemon::CapacityScaling::run() "run()" must be called before 436 /// using them. 529 /// The \ref run() function must be called before using them. 437 530 438 531 /// @{ 439 532 440 /// \brief Return a const reference to the arc map storing the 441 /// found flow. 442 /// 443 /// Return a const reference to the arc map storing the found flow. 533 /// \brief Return the total cost of the found flow. 534 /// 535 /// This function returns the total cost of the found flow. 536 /// Its complexity is O(e). 537 /// 538 /// \note The return type of the function can be specified as a 539 /// template parameter. For example, 540 /// \code 541 /// cs.totalCost<double>(); 542 /// \endcode 543 /// It is useful if the total cost cannot be stored in the \c Cost 544 /// type of the algorithm, which is the default return type of the 545 /// function. 444 546 /// 445 547 /// \pre \ref run() must be called before using this function. 446 const FlowMap& flowMap() const { 447 return *_flow; 448 } 449 450 /// \brief Return a const reference to the node map storing the 451 /// found potentials (the dual solution). 452 /// 453 /// Return a const reference to the node map storing the found 454 /// potentials (the dual solution). 548 template <typename Number> 549 Number totalCost() const { 550 Number c = 0; 551 for (ArcIt a(_graph); a != INVALID; ++a) { 552 int i = _arc_idb[a]; 553 c += static_cast<Number>(_res_cap[i]) * 554 (static_cast<Number>(_cost[i])); 555 } 556 return c; 557 } 558 559 #ifndef DOXYGEN 560 Cost totalCost() const { 561 return totalCost<Cost>(); 562 } 563 #endif 564 565 /// \brief Return the flow on the given arc. 566 /// 567 /// This function returns the flow on the given arc. 455 568 /// 456 569 /// \pre \ref run() must be called before using this function. 457 const PotentialMap& potentialMap() const { 458 return *_potential; 459 } 460 461 /// \brief Return the flow on the given arc. 462 /// 463 /// Return the flow on the given arc. 570 Value flow(const Arc& a) const { 571 return _res_cap[_arc_idb[a]]; 572 } 573 574 /// \brief Return the flow map (the primal solution). 575 /// 576 /// This function copies the flow value on each arc into the given 577 /// map. The \c Value type of the algorithm must be convertible to 578 /// the \c Value type of the map. 464 579 /// 465 580 /// \pre \ref run() must be called before using this function. 466 Capacity flow(const Arc& arc) const { 467 return (*_flow)[arc]; 468 } 469 470 /// \brief Return the potential of the given node. 471 /// 472 /// Return the potential of the given node. 581 template <typename FlowMap> 582 void flowMap(FlowMap &map) const { 583 for (ArcIt a(_graph); a != INVALID; ++a) { 584 map.set(a, _res_cap[_arc_idb[a]]); 585 } 586 } 587 588 /// \brief Return the potential (dual value) of the given node. 589 /// 590 /// This function returns the potential (dual value) of the 591 /// given node. 473 592 /// 474 593 /// \pre \ref run() must be called before using this function. 475 Cost potential(const Node& node) const { 476 return (*_potential)[node]; 477 } 478 479 /// \brief Return the total cost of the found flow. 480 /// 481 /// Return the total cost of the found flow. The complexity of the 482 /// function is \f$ O(e) \f$. 594 Cost potential(const Node& n) const { 595 return _pi[_node_id[n]]; 596 } 597 598 /// \brief Return the potential map (the dual solution). 599 /// 600 /// This function copies the potential (dual value) of each node 601 /// into the given map. 602 /// The \c Cost type of the algorithm must be convertible to the 603 /// \c Value type of the map. 483 604 /// 484 605 /// \pre \ref run() must be called before using this function. 485 Cost totalCost() const {486 Cost c = 0;487 for ( ArcIt e(_graph); e != INVALID; ++e)488 c += (*_flow)[e] * _cost[e];489 return c;606 template <typename PotentialMap> 607 void potentialMap(PotentialMap &map) const { 608 for (NodeIt n(_graph); n != INVALID; ++n) { 609 map.set(n, _pi[_node_id[n]]); 610 } 490 611 } 491 612 … … 494 615 private: 495 616 496 /// Initialize the algorithm. 497 bool init(bool scaling) { 498 if (!_valid_supply) return false; 499 500 // Initializing maps 501 if (!_flow) { 502 _flow = new FlowMap(_graph); 503 _local_flow = true; 504 } 505 if (!_potential) { 506 _potential = new PotentialMap(_graph); 507 _local_potential = true; 508 } 509 for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; 510 for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; 511 512 _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost, 513 _excess, *_potential, _pred ); 514 515 // Initializing delta value 617 // Initialize the algorithm 618 ProblemType init(bool scaling) { 619 if (_node_num == 0) return INFEASIBLE; 620 621 // Check the sum of supply values 622 _sum_supply = 0; 623 for (int i = 0; i != _root; ++i) { 624 _sum_supply += _supply[i]; 625 } 626 if (_sum_supply > 0) return INFEASIBLE; 627 628 // Initialize maps 629 for (int i = 0; i != _root; ++i) { 630 _pi[i] = 0; 631 _excess[i] = _supply[i]; 632 } 633 634 // Remove nonzero lower bounds 635 if (_have_lower) { 636 for (int i = 0; i != _root; ++i) { 637 for (int j = _first_out[i]; j != _first_out[i+1]; ++j) { 638 if (_forward[j]) { 639 Value c = _lower[j]; 640 if (c >= 0) { 641 _res_cap[j] = _upper[j] < INF ? _upper[j]  c : INF; 642 } else { 643 _res_cap[j] = _upper[j] < INF + c ? _upper[j]  c : INF; 644 } 645 _excess[i] = c; 646 _excess[_target[j]] += c; 647 } else { 648 _res_cap[j] = 0; 649 } 650 } 651 } 652 } else { 653 for (int j = 0; j != _res_arc_num; ++j) { 654 _res_cap[j] = _forward[j] ? _upper[j] : 0; 655 } 656 } 657 658 // Handle negative costs 659 for (int u = 0; u != _root; ++u) { 660 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { 661 Value rc = _res_cap[a]; 662 if (_cost[a] < 0 && rc > 0) { 663 if (rc == INF) return UNBOUNDED; 664 _excess[u] = rc; 665 _excess[_target[a]] += rc; 666 _res_cap[a] = 0; 667 _res_cap[_reverse[a]] += rc; 668 } 669 } 670 } 671 672 // Handle GEQ supply type 673 if (_sum_supply < 0) { 674 _pi[_root] = 0; 675 _excess[_root] = _sum_supply; 676 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 677 int u = _target[a]; 678 if (_excess[u] < 0) { 679 _res_cap[a] = _excess[u] + 1; 680 } else { 681 _res_cap[a] = 1; 682 } 683 _res_cap[_reverse[a]] = 0; 684 _cost[a] = 0; 685 _cost[_reverse[a]] = 0; 686 } 687 } else { 688 _pi[_root] = 0; 689 _excess[_root] = 0; 690 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 691 _res_cap[a] = 1; 692 _res_cap[_reverse[a]] = 0; 693 _cost[a] = 0; 694 _cost[_reverse[a]] = 0; 695 } 696 } 697 698 // Initialize delta value 516 699 if (scaling) { 517 700 // With scaling 518 Supplymax_sup = 0, max_dem = 0;519 for ( NodeIt n(_graph); n != INVALID; ++n) {520 if ( _ supply[n] > max_sup) max_sup = _supply[n];521 if (_ supply[n] > max_dem) max_dem = _supply[n];522 } 523 Capacitymax_cap = 0;524 for ( ArcIt e(_graph); e != INVALID; ++e) {525 if (_ capacity[e] > max_cap) max_cap = _capacity[e];701 Value max_sup = 0, max_dem = 0; 702 for (int i = 0; i != _node_num; ++i) { 703 if ( _excess[i] > max_sup) max_sup = _excess[i]; 704 if (_excess[i] > max_dem) max_dem = _excess[i]; 705 } 706 Value max_cap = 0; 707 for (int j = 0; j != _res_arc_num; ++j) { 708 if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; 526 709 } 527 710 max_sup = std::min(std::min(max_sup, max_dem), max_cap); … … 534 717 } 535 718 536 return true; 537 } 538 539 bool start() { 719 return OPTIMAL; 720 } 721 722 ProblemType start() { 723 // Execute the algorithm 724 ProblemType pt; 540 725 if (_delta > 1) 541 returnstartWithScaling();726 pt = startWithScaling(); 542 727 else 543 return startWithoutScaling(); 544 } 545 546 /// Execute the capacity scaling algorithm. 547 bool startWithScaling() { 548 // Processing capacity scaling phases 549 Node s, t; 728 pt = startWithoutScaling(); 729 730 // Handle nonzero lower bounds 731 if (_have_lower) { 732 for (int j = 0; j != _res_arc_num  _node_num + 1; ++j) { 733 if (!_forward[j]) _res_cap[j] += _lower[j]; 734 } 735 } 736 737 // Shift potentials if necessary 738 Cost pr = _pi[_root]; 739 if (_sum_supply < 0  pr > 0) { 740 for (int i = 0; i != _node_num; ++i) { 741 _pi[i] = pr; 742 } 743 } 744 745 return pt; 746 } 747 748 // Execute the capacity scaling algorithm 749 ProblemType startWithScaling() { 750 // Process capacity scaling phases 751 int s, t; 550 752 int phase_cnt = 0; 551 753 int factor = 4; 754 ResidualDijkstra _dijkstra(*this); 552 755 while (true) { 553 // Saturating all arcs not satisfying the optimality condition 554 for (ArcIt e(_graph); e != INVALID; ++e) { 555 Node u = _graph.source(e), v = _graph.target(e); 556 Cost c = _cost[e] + (*_potential)[u]  (*_potential)[v]; 557 if (c < 0 && _res_cap[e] >= _delta) { 558 _excess[u] = _res_cap[e]; 559 _excess[v] += _res_cap[e]; 560 (*_flow)[e] = _capacity[e]; 561 _res_cap[e] = 0; 562 } 563 else if (c > 0 && (*_flow)[e] >= _delta) { 564 _excess[u] += (*_flow)[e]; 565 _excess[v] = (*_flow)[e]; 566 (*_flow)[e] = 0; 567 _res_cap[e] = _capacity[e]; 568 } 569 } 570 571 // Finding excess nodes and deficit nodes 756 // Saturate all arcs not satisfying the optimality condition 757 for (int u = 0; u != _node_num; ++u) { 758 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { 759 int v = _target[a]; 760 Cost c = _cost[a] + _pi[u]  _pi[v]; 761 Value rc = _res_cap[a]; 762 if (c < 0 && rc >= _delta) { 763 _excess[u] = rc; 764 _excess[v] += rc; 765 _res_cap[a] = 0; 766 _res_cap[_reverse[a]] += rc; 767 } 768 } 769 } 770 771 // Find excess nodes and deficit nodes 572 772 _excess_nodes.clear(); 573 773 _deficit_nodes.clear(); 574 for ( NodeIt n(_graph); n != INVALID; ++n) {575 if (_excess[ n] >= _delta) _excess_nodes.push_back(n);576 if (_excess[ n] <= _delta) _deficit_nodes.push_back(n);774 for (int u = 0; u != _node_num; ++u) { 775 if (_excess[u] >= _delta) _excess_nodes.push_back(u); 776 if (_excess[u] <= _delta) _deficit_nodes.push_back(u); 577 777 } 578 778 int next_node = 0, next_def_node = 0; 579 779 580 // Find ingaugmenting shortest paths780 // Find augmenting shortest paths 581 781 while (next_node < int(_excess_nodes.size())) { 582 // Check ingdeficit nodes782 // Check deficit nodes 583 783 if (_delta > 1) { 584 784 bool delta_deficit = false; … … 593 793 } 594 794 595 // Run ning Dijkstra795 // Run Dijkstra in the residual network 596 796 s = _excess_nodes[next_node]; 597 if ((t = _dijkstra >run(s, _delta)) == INVALID) {797 if ((t = _dijkstra.run(s, _delta)) == 1) { 598 798 if (_delta > 1) { 599 799 ++next_node; 600 800 continue; 601 801 } 602 return false;603 } 604 605 // Augment ing along a shortest path from s to t.606 Capacityd = std::min(_excess[s], _excess[t]);607 Nodeu = t;608 Arc e;802 return INFEASIBLE; 803 } 804 805 // Augment along a shortest path from s to t 806 Value d = std::min(_excess[s], _excess[t]); 807 int u = t; 808 int a; 609 809 if (d > _delta) { 610 while ((e = _pred[u]) != INVALID) { 611 Capacity rc; 612 if (u == _graph.target(e)) { 613 rc = _res_cap[e]; 614 u = _graph.source(e); 615 } else { 616 rc = (*_flow)[e]; 617 u = _graph.target(e); 618 } 619 if (rc < d) d = rc; 810 while ((a = _pred[u]) != 1) { 811 if (_res_cap[a] < d) d = _res_cap[a]; 812 u = _source[a]; 620 813 } 621 814 } 622 815 u = t; 623 while ((e = _pred[u]) != INVALID) { 624 if (u == _graph.target(e)) { 625 (*_flow)[e] += d; 626 _res_cap[e] = d; 627 u = _graph.source(e); 628 } else { 629 (*_flow)[e] = d; 630 _res_cap[e] += d; 631 u = _graph.target(e); 632 } 816 while ((a = _pred[u]) != 1) { 817 _res_cap[a] = d; 818 _res_cap[_reverse[a]] += d; 819 u = _source[a]; 633 820 } 634 821 _excess[s] = d; … … 639 826 640 827 if (_delta == 1) break; 641 if (++phase_cnt >_phase_num / 4) factor = 2;828 if (++phase_cnt == _phase_num / 4) factor = 2; 642 829 _delta = _delta <= factor ? 1 : _delta / factor; 643 830 } 644 831 645 // Handling nonzero lower bounds 646 if (_lower) { 647 for (ArcIt e(_graph); e != INVALID; ++e) 648 (*_flow)[e] += (*_lower)[e]; 649 } 650 return true; 651 } 652 653 /// Execute the successive shortest path algorithm. 654 bool startWithoutScaling() { 655 // Finding excess nodes 656 for (NodeIt n(_graph); n != INVALID; ++n) 657 if (_excess[n] > 0) _excess_nodes.push_back(n); 658 if (_excess_nodes.size() == 0) return true; 832 return OPTIMAL; 833 } 834 835 // Execute the successive shortest path algorithm 836 ProblemType startWithoutScaling() { 837 // Find excess nodes 838 _excess_nodes.clear(); 839 for (int i = 0; i != _node_num; ++i) { 840 if (_excess[i] > 0) _excess_nodes.push_back(i); 841 } 842 if (_excess_nodes.size() == 0) return OPTIMAL; 659 843 int next_node = 0; 660 844 661 // Finding shortest paths 662 Node s, t; 845 // Find shortest paths 846 int s, t; 847 ResidualDijkstra _dijkstra(*this); 663 848 while ( _excess[_excess_nodes[next_node]] > 0  664 849 ++next_node < int(_excess_nodes.size()) ) 665 850 { 666 // Run ning Dijkstra851 // Run Dijkstra in the residual network 667 852 s = _excess_nodes[next_node]; 668 if ((t = _dijkstra >run(s)) == INVALID) return false;669 670 // Augment ingalong a shortest path from s to t671 Capacityd = std::min(_excess[s], _excess[t]);672 Nodeu = t;673 Arc e;853 if ((t = _dijkstra.run(s)) == 1) return INFEASIBLE; 854 855 // Augment along a shortest path from s to t 856 Value d = std::min(_excess[s], _excess[t]); 857 int u = t; 858 int a; 674 859 if (d > 1) { 675 while ((e = _pred[u]) != INVALID) { 676 Capacity rc; 677 if (u == _graph.target(e)) { 678 rc = _res_cap[e]; 679 u = _graph.source(e); 680 } else { 681 rc = (*_flow)[e]; 682 u = _graph.target(e); 683 } 684 if (rc < d) d = rc; 860 while ((a = _pred[u]) != 1) { 861 if (_res_cap[a] < d) d = _res_cap[a]; 862 u = _source[a]; 685 863 } 686 864 } 687 865 u = t; 688 while ((e = _pred[u]) != INVALID) { 689 if (u == _graph.target(e)) { 690 (*_flow)[e] += d; 691 _res_cap[e] = d; 692 u = _graph.source(e); 693 } else { 694 (*_flow)[e] = d; 695 _res_cap[e] += d; 696 u = _graph.target(e); 697 } 866 while ((a = _pred[u]) != 1) { 867 _res_cap[a] = d; 868 _res_cap[_reverse[a]] += d; 869 u = _source[a]; 698 870 } 699 871 _excess[s] = d; … … 701 873 } 702 874 703 // Handling nonzero lower bounds 704 if (_lower) { 705 for (ArcIt e(_graph); e != INVALID; ++e) 706 (*_flow)[e] += (*_lower)[e]; 707 } 708 return true; 875 return OPTIMAL; 709 876 } 710 877
Note: See TracChangeset
for help on using the changeset viewer.