Changeset 881:aef153f430e1 in lemon for lemon/cycle_canceling.h
 Timestamp:
 11/13/09 00:10:33 (10 years ago)
 Branch:
 default
 Phase:
 public
 Rebase:
 62373361373236626665663832376334666366353631373365333236323466343130303133633939
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/cycle_canceling.h
r880 r881 20 20 #define LEMON_CYCLE_CANCELING_H 21 21 22 /// \ingroup min_cost_flow 23 /// 22 /// \ingroup min_cost_flow_algs 24 23 /// \file 25 /// \brief Cyclecanceling algorithm for finding a minimum cost flow.24 /// \brief Cyclecanceling algorithms for finding a minimum cost flow. 26 25 27 26 #include <vector> 27 #include <limits> 28 29 #include <lemon/core.h> 30 #include <lemon/maps.h> 31 #include <lemon/path.h> 32 #include <lemon/math.h> 33 #include <lemon/static_graph.h> 28 34 #include <lemon/adaptors.h> 29 #include <lemon/path.h>30 31 35 #include <lemon/circulation.h> 32 36 #include <lemon/bellman_ford.h> … … 35 39 namespace lemon { 36 40 37 /// \addtogroup min_cost_flow 41 /// \addtogroup min_cost_flow_algs 38 42 /// @{ 39 43 40 /// \brief Implementation of a cyclecanceling algorithmfor41 /// finding a minimum cost flow.44 /// \brief Implementation of cyclecanceling algorithms for 45 /// finding a \ref min_cost_flow "minimum cost flow". 42 46 /// 43 /// \ref CycleCanceling implements a cyclecanceling algorithm for 44 /// finding a minimum cost flow. 47 /// \ref CycleCanceling implements three different cyclecanceling 48 /// algorithms for finding a \ref min_cost_flow "minimum cost flow". 49 /// The most efficent one (both theoretically and practically) 50 /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm, 51 /// thus it is the default method. 52 /// It is strongly polynomial, but in practice, it is typically much 53 /// slower than the scaling algorithms and NetworkSimplex. 45 54 /// 46 /// \tparam Digraph The digraph type the algorithm runs on. 47 /// \tparam LowerMap The type of the lower bound map. 48 /// \tparam CapacityMap The type of the capacity (upper bound) map. 49 /// \tparam CostMap The type of the cost (length) map. 50 /// \tparam SupplyMap The type of the supply map. 55 /// Most of the parameters of the problem (except for the digraph) 56 /// can be given using separate functions, and the algorithm can be 57 /// executed using the \ref run() function. If some parameters are not 58 /// specified, then default values will be used. 51 59 /// 52 /// \ warning53 ///  Arc capacities and costs should be \e nonnegative \e integers.54 ///  Supply values should be \e signed \e integers.55 ///  The value types of the maps should be convertible to each other.56 ///  \c CostMap::Value must be signed type.60 /// \tparam GR The digraph type the algorithm runs on. 61 /// \tparam V The number type used for flow amounts, capacity bounds 62 /// and supply values in the algorithm. By default, it is \c int. 63 /// \tparam C The number type used for costs and potentials in the 64 /// algorithm. By default, it is the same as \c V. 57 65 /// 58 /// \note By default the \ref BellmanFord "BellmanFord" algorithm is 59 /// used for negative cycle detection with limited iteration number. 60 /// However \ref CycleCanceling also provides the "Minimum Mean 61 /// CycleCanceling" algorithm, which is \e strongly \e polynomial, 62 /// but rather slower in practice. 63 /// To use this version of the algorithm, call \ref run() with \c true 64 /// parameter. 66 /// \warning Both number types must be signed and all input data must 67 /// be integer. 68 /// \warning This algorithm does not support negative costs for such 69 /// arcs that have infinite upper bound. 65 70 /// 66 /// \author Peter Kovacs 67 template < typename Digraph, 68 typename LowerMap = typename Digraph::template ArcMap<int>, 69 typename CapacityMap = typename Digraph::template ArcMap<int>, 70 typename CostMap = typename Digraph::template ArcMap<int>, 71 typename SupplyMap = typename Digraph::template NodeMap<int> > 71 /// \note For more information about the three available methods, 72 /// see \ref Method. 73 #ifdef DOXYGEN 74 template <typename GR, typename V, typename C> 75 #else 76 template <typename GR, typename V = int, typename C = V> 77 #endif 72 78 class CycleCanceling 73 79 { 74 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);75 76 typedef typename CapacityMap::Value Capacity;77 typedef typename CostMap::Value Cost;78 typedef typename SupplyMap::Value Supply;79 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;80 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;81 82 typedef ResidualDigraph< const Digraph,83 CapacityArcMap, CapacityArcMap > ResDigraph;84 typedef typename ResDigraph::Node ResNode;85 typedef typename ResDigraph::NodeIt ResNodeIt;86 typedef typename ResDigraph::Arc ResArc;87 typedef typename ResDigraph::ArcIt ResArcIt;88 89 80 public: 90 81 91 /// The type of the flow map. 92 typedef typename Digraph::template ArcMap<Capacity> FlowMap; 93 /// The type of the potential map. 94 typedef typename Digraph::template NodeMap<Cost> PotentialMap; 82 /// The type of the digraph 83 typedef GR Digraph; 84 /// The type of the flow amounts, capacity bounds and supply values 85 typedef V Value; 86 /// The type of the arc costs 87 typedef C Cost; 88 89 public: 90 91 /// \brief Problem type constants for the \c run() function. 92 /// 93 /// Enum type containing the problem type constants that can be 94 /// returned by the \ref run() function of the algorithm. 95 enum ProblemType { 96 /// The problem has no feasible solution (flow). 97 INFEASIBLE, 98 /// The problem has optimal solution (i.e. it is feasible and 99 /// bounded), and the algorithm has found optimal flow and node 100 /// potentials (primal and dual solutions). 101 OPTIMAL, 102 /// The digraph contains an arc of negative cost and infinite 103 /// upper bound. It means that the objective function is unbounded 104 /// on that arc, however, note that it could actually be bounded 105 /// over the feasible flows, but this algroithm cannot handle 106 /// these cases. 107 UNBOUNDED 108 }; 109 110 /// \brief Constants for selecting the used method. 111 /// 112 /// Enum type containing constants for selecting the used method 113 /// for the \ref run() function. 114 /// 115 /// \ref CycleCanceling provides three different cyclecanceling 116 /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" 117 /// is used, which proved to be the most efficient and the most robust 118 /// on various test inputs. 119 /// However, the other methods can be selected using the \ref run() 120 /// function with the proper parameter. 121 enum Method { 122 /// A simple cyclecanceling method, which uses the 123 /// \ref BellmanFord "BellmanFord" algorithm with limited iteration 124 /// number for detecting negative cycles in the residual network. 125 SIMPLE_CYCLE_CANCELING, 126 /// The "Minimum Mean CycleCanceling" algorithm, which is a 127 /// wellknown strongly polynomial method. It improves along a 128 /// \ref min_mean_cycle "minimum mean cycle" in each iteration. 129 /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)). 130 MINIMUM_MEAN_CYCLE_CANCELING, 131 /// The "Cancel And Tighten" algorithm, which can be viewed as an 132 /// improved version of the previous method. 133 /// It is faster both in theory and in practice, its running time 134 /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)). 135 CANCEL_AND_TIGHTEN 136 }; 95 137 96 138 private: 97 139 98 /// \brief Map adaptor class for handling residual arc costs. 99 /// 100 /// Map adaptor class for handling residual arc costs. 101 class ResidualCostMap : public MapBase<ResArc, Cost> 140 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 141 142 typedef std::vector<int> IntVector; 143 typedef std::vector<char> CharVector; 144 typedef std::vector<double> DoubleVector; 145 typedef std::vector<Value> ValueVector; 146 typedef std::vector<Cost> CostVector; 147 148 private: 149 150 template <typename KT, typename VT> 151 class VectorMap { 152 public: 153 typedef KT Key; 154 typedef VT Value; 155 156 VectorMap(std::vector<Value>& v) : _v(v) {} 157 158 const Value& operator[](const Key& key) const { 159 return _v[StaticDigraph::id(key)]; 160 } 161 162 Value& operator[](const Key& key) { 163 return _v[StaticDigraph::id(key)]; 164 } 165 166 void set(const Key& key, const Value& val) { 167 _v[StaticDigraph::id(key)] = val; 168 } 169 170 private: 171 std::vector<Value>& _v; 172 }; 173 174 typedef VectorMap<StaticDigraph::Node, Cost> CostNodeMap; 175 typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap; 176 177 private: 178 179 180 // Data related to the underlying digraph 181 const GR &_graph; 182 int _node_num; 183 int _arc_num; 184 int _res_node_num; 185 int _res_arc_num; 186 int _root; 187 188 // Parameters of the problem 189 bool _have_lower; 190 Value _sum_supply; 191 192 // Data structures for storing the digraph 193 IntNodeMap _node_id; 194 IntArcMap _arc_idf; 195 IntArcMap _arc_idb; 196 IntVector _first_out; 197 CharVector _forward; 198 IntVector _source; 199 IntVector _target; 200 IntVector _reverse; 201 202 // Node and arc data 203 ValueVector _lower; 204 ValueVector _upper; 205 CostVector _cost; 206 ValueVector _supply; 207 208 ValueVector _res_cap; 209 CostVector _pi; 210 211 // Data for a StaticDigraph structure 212 typedef std::pair<int, int> IntPair; 213 StaticDigraph _sgr; 214 std::vector<IntPair> _arc_vec; 215 std::vector<Cost> _cost_vec; 216 IntVector _id_vec; 217 CostArcMap _cost_map; 218 CostNodeMap _pi_map; 219 220 public: 221 222 /// \brief Constant for infinite upper bounds (capacities). 223 /// 224 /// Constant for infinite upper bounds (capacities). 225 /// It is \c std::numeric_limits<Value>::infinity() if available, 226 /// \c std::numeric_limits<Value>::max() otherwise. 227 const Value INF; 228 229 public: 230 231 /// \brief Constructor. 232 /// 233 /// The constructor of the class. 234 /// 235 /// \param graph The digraph the algorithm runs on. 236 CycleCanceling(const GR& graph) : 237 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), 238 _cost_map(_cost_vec), _pi_map(_pi), 239 INF(std::numeric_limits<Value>::has_infinity ? 240 std::numeric_limits<Value>::infinity() : 241 std::numeric_limits<Value>::max()) 102 242 { 103 private: 104 105 const CostMap &_cost_map; 106 107 public: 108 109 ///\e 110 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {} 111 112 ///\e 113 Cost operator[](const ResArc &e) const { 114 return ResDigraph::forward(e) ? _cost_map[e] : _cost_map[e]; 115 } 116 117 }; //class ResidualCostMap 118 119 private: 120 121 // The maximum number of iterations for the first execution of the 122 // BellmanFord algorithm. It should be at least 2. 123 static const int BF_FIRST_LIMIT = 2; 124 // The iteration limit for the BellmanFord algorithm is multiplied 125 // by BF_LIMIT_FACTOR/100 in every round. 126 static const int BF_LIMIT_FACTOR = 150; 127 128 private: 129 130 // The digraph the algorithm runs on 131 const Digraph &_graph; 132 // The original lower bound map 133 const LowerMap *_lower; 134 // The modified capacity map 135 CapacityArcMap _capacity; 136 // The original cost map 137 const CostMap &_cost; 138 // The modified supply map 139 SupplyNodeMap _supply; 140 bool _valid_supply; 141 142 // Arc map of the current flow 143 FlowMap *_flow; 144 bool _local_flow; 145 // Node map of the current potentials 146 PotentialMap *_potential; 147 bool _local_potential; 148 149 // The residual digraph 150 ResDigraph *_res_graph; 151 // The residual cost map 152 ResidualCostMap _res_cost; 153 154 public: 155 156 /// \brief General constructor (with lower bounds). 157 /// 158 /// General constructor (with lower bounds). 159 /// 160 /// \param digraph The digraph the algorithm runs on. 161 /// \param lower The lower bounds of the arcs. 162 /// \param capacity The capacities (upper bounds) of the arcs. 163 /// \param cost The cost (length) values of the arcs. 164 /// \param supply The supply values of the nodes (signed). 165 CycleCanceling( const Digraph &digraph, 166 const LowerMap &lower, 167 const CapacityMap &capacity, 168 const CostMap &cost, 169 const SupplyMap &supply ) : 170 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost), 171 _supply(digraph), _flow(NULL), _local_flow(false), 172 _potential(NULL), _local_potential(false), 173 _res_graph(NULL), _res_cost(_cost) 174 { 175 // Check the sum of supply values 176 Supply sum = 0; 243 // Check the number types 244 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 245 "The flow type of CycleCanceling must be signed"); 246 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 247 "The cost type of CycleCanceling must be signed"); 248 249 // Resize vectors 250 _node_num = countNodes(_graph); 251 _arc_num = countArcs(_graph); 252 _res_node_num = _node_num + 1; 253 _res_arc_num = 2 * (_arc_num + _node_num); 254 _root = _node_num; 255 256 _first_out.resize(_res_node_num + 1); 257 _forward.resize(_res_arc_num); 258 _source.resize(_res_arc_num); 259 _target.resize(_res_arc_num); 260 _reverse.resize(_res_arc_num); 261 262 _lower.resize(_res_arc_num); 263 _upper.resize(_res_arc_num); 264 _cost.resize(_res_arc_num); 265 _supply.resize(_res_node_num); 266 267 _res_cap.resize(_res_arc_num); 268 _pi.resize(_res_node_num); 269 270 _arc_vec.reserve(_res_arc_num); 271 _cost_vec.reserve(_res_arc_num); 272 _id_vec.reserve(_res_arc_num); 273 274 // Copy the graph 275 int i = 0, j = 0, k = 2 * _arc_num + _node_num; 276 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 277 _node_id[n] = i; 278 } 279 i = 0; 280 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 281 _first_out[i] = j; 282 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { 283 _arc_idf[a] = j; 284 _forward[j] = true; 285 _source[j] = i; 286 _target[j] = _node_id[_graph.runningNode(a)]; 287 } 288 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { 289 _arc_idb[a] = j; 290 _forward[j] = false; 291 _source[j] = i; 292 _target[j] = _node_id[_graph.runningNode(a)]; 293 } 294 _forward[j] = false; 295 _source[j] = i; 296 _target[j] = _root; 297 _reverse[j] = k; 298 _forward[k] = true; 299 _source[k] = _root; 300 _target[k] = i; 301 _reverse[k] = j; 302 ++j; ++k; 303 } 304 _first_out[i] = j; 305 _first_out[_res_node_num] = k; 306 for (ArcIt a(_graph); a != INVALID; ++a) { 307 int fi = _arc_idf[a]; 308 int bi = _arc_idb[a]; 309 _reverse[fi] = bi; 310 _reverse[bi] = fi; 311 } 312 313 // Reset parameters 314 reset(); 315 } 316 317 /// \name Parameters 318 /// The parameters of the algorithm can be specified using these 319 /// functions. 320 321 /// @{ 322 323 /// \brief Set the lower bounds on the arcs. 324 /// 325 /// This function sets the lower bounds on the arcs. 326 /// If it is not used before calling \ref run(), the lower bounds 327 /// will be set to zero on all arcs. 328 /// 329 /// \param map An arc map storing the lower bounds. 330 /// Its \c Value type must be convertible to the \c Value type 331 /// of the algorithm. 332 /// 333 /// \return <tt>(*this)</tt> 334 template <typename LowerMap> 335 CycleCanceling& lowerMap(const LowerMap& map) { 336 _have_lower = true; 337 for (ArcIt a(_graph); a != INVALID; ++a) { 338 _lower[_arc_idf[a]] = map[a]; 339 _lower[_arc_idb[a]] = map[a]; 340 } 341 return *this; 342 } 343 344 /// \brief Set the upper bounds (capacities) on the arcs. 345 /// 346 /// This function sets the upper bounds (capacities) on the arcs. 347 /// If it is not used before calling \ref run(), the upper bounds 348 /// will be set to \ref INF on all arcs (i.e. the flow value will be 349 /// unbounded from above). 350 /// 351 /// \param map An arc map storing the upper bounds. 352 /// Its \c Value type must be convertible to the \c Value type 353 /// of the algorithm. 354 /// 355 /// \return <tt>(*this)</tt> 356 template<typename UpperMap> 357 CycleCanceling& upperMap(const UpperMap& map) { 358 for (ArcIt a(_graph); a != INVALID; ++a) { 359 _upper[_arc_idf[a]] = map[a]; 360 } 361 return *this; 362 } 363 364 /// \brief Set the costs of the arcs. 365 /// 366 /// This function sets the costs of the arcs. 367 /// If it is not used before calling \ref run(), the costs 368 /// will be set to \c 1 on all arcs. 369 /// 370 /// \param map An arc map storing the costs. 371 /// Its \c Value type must be convertible to the \c Cost type 372 /// of the algorithm. 373 /// 374 /// \return <tt>(*this)</tt> 375 template<typename CostMap> 376 CycleCanceling& costMap(const CostMap& map) { 377 for (ArcIt a(_graph); a != INVALID; ++a) { 378 _cost[_arc_idf[a]] = map[a]; 379 _cost[_arc_idb[a]] = map[a]; 380 } 381 return *this; 382 } 383 384 /// \brief Set the supply values of the nodes. 385 /// 386 /// This function sets the supply values of the nodes. 387 /// If neither this function nor \ref stSupply() is used before 388 /// calling \ref run(), the supply of each node will be set to zero. 389 /// 390 /// \param map A node map storing the supply values. 391 /// Its \c Value type must be convertible to the \c Value type 392 /// of the algorithm. 393 /// 394 /// \return <tt>(*this)</tt> 395 template<typename SupplyMap> 396 CycleCanceling& supplyMap(const SupplyMap& map) { 177 397 for (NodeIt n(_graph); n != INVALID; ++n) { 178 _supply[n] = supply[n]; 179 sum += _supply[n]; 180 } 181 _valid_supply = sum == 0; 182 183 // Remove nonzero lower bounds 184 for (ArcIt e(_graph); e != INVALID; ++e) { 185 _capacity[e] = capacity[e]; 186 if (lower[e] != 0) { 187 _capacity[e] = lower[e]; 188 _supply[_graph.source(e)] = lower[e]; 189 _supply[_graph.target(e)] += lower[e]; 190 } 191 } 192 } 193 /* 194 /// \brief General constructor (without lower bounds). 195 /// 196 /// General constructor (without lower bounds). 197 /// 198 /// \param digraph The digraph the algorithm runs on. 199 /// \param capacity The capacities (upper bounds) of the arcs. 200 /// \param cost The cost (length) values of the arcs. 201 /// \param supply The supply values of the nodes (signed). 202 CycleCanceling( const Digraph &digraph, 203 const CapacityMap &capacity, 204 const CostMap &cost, 205 const SupplyMap &supply ) : 206 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), 207 _supply(supply), _flow(NULL), _local_flow(false), 208 _potential(NULL), _local_potential(false), _res_graph(NULL), 209 _res_cost(_cost) 210 { 211 // Check the sum of supply values 212 Supply sum = 0; 213 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; 214 _valid_supply = sum == 0; 215 } 216 217 /// \brief Simple constructor (with lower bounds). 218 /// 219 /// Simple constructor (with lower bounds). 220 /// 221 /// \param digraph The digraph the algorithm runs on. 222 /// \param lower The lower bounds of the arcs. 223 /// \param capacity The capacities (upper bounds) of the arcs. 224 /// \param cost The cost (length) values of the arcs. 398 _supply[_node_id[n]] = map[n]; 399 } 400 return *this; 401 } 402 403 /// \brief Set single source and target nodes and a supply value. 404 /// 405 /// This function sets a single source node and a single target node 406 /// and the required flow value. 407 /// If neither this function nor \ref supplyMap() is used before 408 /// calling \ref run(), the supply of each node will be set to zero. 409 /// 410 /// Using this function has the same effect as using \ref supplyMap() 411 /// with such a map in which \c k is assigned to \c s, \c k is 412 /// assigned to \c t and all other nodes have zero supply value. 413 /// 225 414 /// \param s The source node. 226 415 /// \param t The target node. 227 /// \param flow_value The required amount of flow from node \c s 228 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 229 CycleCanceling( const Digraph &digraph, 230 const LowerMap &lower, 231 const CapacityMap &capacity, 232 const CostMap &cost, 233 Node s, Node t, 234 Supply flow_value ) : 235 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost), 236 _supply(digraph, 0), _flow(NULL), _local_flow(false), 237 _potential(NULL), _local_potential(false), _res_graph(NULL), 238 _res_cost(_cost) 239 { 240 // Remove nonzero lower bounds 241 _supply[s] = flow_value; 242 _supply[t] = flow_value; 243 for (ArcIt e(_graph); e != INVALID; ++e) { 244 if (lower[e] != 0) { 245 _capacity[e] = lower[e]; 246 _supply[_graph.source(e)] = lower[e]; 247 _supply[_graph.target(e)] += lower[e]; 248 } 249 } 250 _valid_supply = true; 251 } 252 253 /// \brief Simple constructor (without lower bounds). 254 /// 255 /// Simple constructor (without lower bounds). 256 /// 257 /// \param digraph The digraph the algorithm runs on. 258 /// \param capacity The capacities (upper bounds) of the arcs. 259 /// \param cost The cost (length) values of the arcs. 260 /// \param s The source node. 261 /// \param t The target node. 262 /// \param flow_value The required amount of flow from node \c s 263 /// to node \c t (i.e. the supply of \c s and the demand of \c t). 264 CycleCanceling( const Digraph &digraph, 265 const CapacityMap &capacity, 266 const CostMap &cost, 267 Node s, Node t, 268 Supply flow_value ) : 269 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), 270 _supply(digraph, 0), _flow(NULL), _local_flow(false), 271 _potential(NULL), _local_potential(false), _res_graph(NULL), 272 _res_cost(_cost) 273 { 274 _supply[s] = flow_value; 275 _supply[t] = flow_value; 276 _valid_supply = true; 277 } 278 */ 279 /// Destructor. 280 ~CycleCanceling() { 281 if (_local_flow) delete _flow; 282 if (_local_potential) delete _potential; 283 delete _res_graph; 284 } 285 286 /// \brief Set the flow map. 287 /// 288 /// Set the flow map. 289 /// 290 /// \return \c (*this) 291 CycleCanceling& flowMap(FlowMap &map) { 292 if (_local_flow) { 293 delete _flow; 294 _local_flow = false; 295 } 296 _flow = ↦ 416 /// \param k The required amount of flow from node \c s to node \c t 417 /// (i.e. the supply of \c s and the demand of \c t). 418 /// 419 /// \return <tt>(*this)</tt> 420 CycleCanceling& stSupply(const Node& s, const Node& t, Value k) { 421 for (int i = 0; i != _res_node_num; ++i) { 422 _supply[i] = 0; 423 } 424 _supply[_node_id[s]] = k; 425 _supply[_node_id[t]] = k; 297 426 return *this; 298 427 } 299 300 /// \brief Set the potential map. 301 /// 302 /// Set the potential map. 303 /// 304 /// \return \c (*this) 305 CycleCanceling& potentialMap(PotentialMap &map) { 306 if (_local_potential) { 307 delete _potential; 308 _local_potential = false; 309 } 310 _potential = ↦ 428 429 /// @} 430 431 /// \name Execution control 432 /// The algorithm can be executed using \ref run(). 433 434 /// @{ 435 436 /// \brief Run the algorithm. 437 /// 438 /// This function runs the algorithm. 439 /// The paramters can be specified using functions \ref lowerMap(), 440 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). 441 /// For example, 442 /// \code 443 /// CycleCanceling<ListDigraph> cc(graph); 444 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) 445 /// .supplyMap(sup).run(); 446 /// \endcode 447 /// 448 /// This function can be called more than once. All the parameters 449 /// that have been given are kept for the next call, unless 450 /// \ref reset() is called, thus only the modified parameters 451 /// have to be set again. See \ref reset() for examples. 452 /// However, the underlying digraph must not be modified after this 453 /// class have been constructed, since it copies and extends the graph. 454 /// 455 /// \param method The cyclecanceling method that will be used. 456 /// For more information, see \ref Method. 457 /// 458 /// \return \c INFEASIBLE if no feasible flow exists, 459 /// \n \c OPTIMAL if the problem has optimal solution 460 /// (i.e. it is feasible and bounded), and the algorithm has found 461 /// optimal flow and node potentials (primal and dual solutions), 462 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost 463 /// and infinite upper bound. It means that the objective function 464 /// is unbounded on that arc, however, note that it could actually be 465 /// bounded over the feasible flows, but this algroithm cannot handle 466 /// these cases. 467 /// 468 /// \see ProblemType, Method 469 ProblemType run(Method method = CANCEL_AND_TIGHTEN) { 470 ProblemType pt = init(); 471 if (pt != OPTIMAL) return pt; 472 start(method); 473 return OPTIMAL; 474 } 475 476 /// \brief Reset all the parameters that have been given before. 477 /// 478 /// This function resets all the paramaters that have been given 479 /// before using functions \ref lowerMap(), \ref upperMap(), 480 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). 481 /// 482 /// It is useful for multiple run() calls. If this function is not 483 /// used, all the parameters given before are kept for the next 484 /// \ref run() call. 485 /// However, the underlying digraph must not be modified after this 486 /// class have been constructed, since it copies and extends the graph. 487 /// 488 /// For example, 489 /// \code 490 /// CycleCanceling<ListDigraph> cs(graph); 491 /// 492 /// // First run 493 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) 494 /// .supplyMap(sup).run(); 495 /// 496 /// // Run again with modified cost map (reset() is not called, 497 /// // so only the cost map have to be set again) 498 /// cost[e] += 100; 499 /// cc.costMap(cost).run(); 500 /// 501 /// // Run again from scratch using reset() 502 /// // (the lower bounds will be set to zero on all arcs) 503 /// cc.reset(); 504 /// cc.upperMap(capacity).costMap(cost) 505 /// .supplyMap(sup).run(); 506 /// \endcode 507 /// 508 /// \return <tt>(*this)</tt> 509 CycleCanceling& reset() { 510 for (int i = 0; i != _res_node_num; ++i) { 511 _supply[i] = 0; 512 } 513 int limit = _first_out[_root]; 514 for (int j = 0; j != limit; ++j) { 515 _lower[j] = 0; 516 _upper[j] = INF; 517 _cost[j] = _forward[j] ? 1 : 1; 518 } 519 for (int j = limit; j != _res_arc_num; ++j) { 520 _lower[j] = 0; 521 _upper[j] = INF; 522 _cost[j] = 0; 523 _cost[_reverse[j]] = 0; 524 } 525 _have_lower = false; 311 526 return *this; 312 527 } 313 528 314 /// \name Execution control 529 /// @} 530 531 /// \name Query Functions 532 /// The results of the algorithm can be obtained using these 533 /// functions.\n 534 /// The \ref run() function must be called before using them. 315 535 316 536 /// @{ 317 537 318 /// \brief Run the algorithm. 319 /// 320 /// Run the algorithm. 321 /// 322 /// \param min_mean_cc Set this parameter to \c true to run the 323 /// "Minimum Mean CycleCanceling" algorithm, which is strongly 324 /// polynomial, but rather slower in practice. 325 /// 326 /// \return \c true if a feasible flow can be found. 327 bool run(bool min_mean_cc = false) { 328 return init() && start(min_mean_cc); 538 /// \brief Return the total cost of the found flow. 539 /// 540 /// This function returns the total cost of the found flow. 541 /// Its complexity is O(e). 542 /// 543 /// \note The return type of the function can be specified as a 544 /// template parameter. For example, 545 /// \code 546 /// cc.totalCost<double>(); 547 /// \endcode 548 /// It is useful if the total cost cannot be stored in the \c Cost 549 /// type of the algorithm, which is the default return type of the 550 /// function. 551 /// 552 /// \pre \ref run() must be called before using this function. 553 template <typename Number> 554 Number totalCost() const { 555 Number c = 0; 556 for (ArcIt a(_graph); a != INVALID; ++a) { 557 int i = _arc_idb[a]; 558 c += static_cast<Number>(_res_cap[i]) * 559 (static_cast<Number>(_cost[i])); 560 } 561 return c; 562 } 563 564 #ifndef DOXYGEN 565 Cost totalCost() const { 566 return totalCost<Cost>(); 567 } 568 #endif 569 570 /// \brief Return the flow on the given arc. 571 /// 572 /// This function returns the flow on the given arc. 573 /// 574 /// \pre \ref run() must be called before using this function. 575 Value flow(const Arc& a) const { 576 return _res_cap[_arc_idb[a]]; 577 } 578 579 /// \brief Return the flow map (the primal solution). 580 /// 581 /// This function copies the flow value on each arc into the given 582 /// map. The \c Value type of the algorithm must be convertible to 583 /// the \c Value type of the map. 584 /// 585 /// \pre \ref run() must be called before using this function. 586 template <typename FlowMap> 587 void flowMap(FlowMap &map) const { 588 for (ArcIt a(_graph); a != INVALID; ++a) { 589 map.set(a, _res_cap[_arc_idb[a]]); 590 } 591 } 592 593 /// \brief Return the potential (dual value) of the given node. 594 /// 595 /// This function returns the potential (dual value) of the 596 /// given node. 597 /// 598 /// \pre \ref run() must be called before using this function. 599 Cost potential(const Node& n) const { 600 return static_cast<Cost>(_pi[_node_id[n]]); 601 } 602 603 /// \brief Return the potential map (the dual solution). 604 /// 605 /// This function copies the potential (dual value) of each node 606 /// into the given map. 607 /// The \c Cost type of the algorithm must be convertible to the 608 /// \c Value type of the map. 609 /// 610 /// \pre \ref run() must be called before using this function. 611 template <typename PotentialMap> 612 void potentialMap(PotentialMap &map) const { 613 for (NodeIt n(_graph); n != INVALID; ++n) { 614 map.set(n, static_cast<Cost>(_pi[_node_id[n]])); 615 } 329 616 } 330 617 331 618 /// @} 332 619 333 /// \name Query Functions334 /// The result of the algorithm can be obtained using these335 /// functions.\n336 /// \ref lemon::CycleCanceling::run() "run()" must be called before337 /// using them.338 339 /// @{340 341 /// \brief Return a const reference to the arc map storing the342 /// found flow.343 ///344 /// Return a const reference to the arc map storing the found flow.345 ///346 /// \pre \ref run() must be called before using this function.347 const FlowMap& flowMap() const {348 return *_flow;349 }350 351 /// \brief Return a const reference to the node map storing the352 /// found potentials (the dual solution).353 ///354 /// Return a const reference to the node map storing the found355 /// potentials (the dual solution).356 ///357 /// \pre \ref run() must be called before using this function.358 const PotentialMap& potentialMap() const {359 return *_potential;360 }361 362 /// \brief Return the flow on the given arc.363 ///364 /// Return the flow on the given arc.365 ///366 /// \pre \ref run() must be called before using this function.367 Capacity flow(const Arc& arc) const {368 return (*_flow)[arc];369 }370 371 /// \brief Return the potential of the given node.372 ///373 /// Return the potential of the given node.374 ///375 /// \pre \ref run() must be called before using this function.376 Cost potential(const Node& node) const {377 return (*_potential)[node];378 }379 380 /// \brief Return the total cost of the found flow.381 ///382 /// Return the total cost of the found flow. The complexity of the383 /// function is \f$ O(e) \f$.384 ///385 /// \pre \ref run() must be called before using this function.386 Cost totalCost() const {387 Cost c = 0;388 for (ArcIt e(_graph); e != INVALID; ++e)389 c += (*_flow)[e] * _cost[e];390 return c;391 }392 393 /// @}394 395 620 private: 396 621 397 /// Initialize the algorithm. 398 bool init() { 399 if (!_valid_supply) return false; 400 401 // Initializing flow and potential maps 402 if (!_flow) { 403 _flow = new FlowMap(_graph); 404 _local_flow = true; 405 } 406 if (!_potential) { 407 _potential = new PotentialMap(_graph); 408 _local_potential = true; 409 } 410 411 _res_graph = new ResDigraph(_graph, _capacity, *_flow); 412 413 // Finding a feasible flow using Circulation 414 Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap, 415 SupplyMap > 416 circulation( _graph, constMap<Arc>(Capacity(0)), _capacity, 417 _supply ); 418 return circulation.flowMap(*_flow).run(); 419 } 420 421 bool start(bool min_mean_cc) { 422 if (min_mean_cc) 423 startMinMean(); 424 else 425 start(); 426 427 // Handling nonzero lower bounds 428 if (_lower) { 429 for (ArcIt e(_graph); e != INVALID; ++e) 430 (*_flow)[e] += (*_lower)[e]; 431 } 432 return true; 433 } 434 435 /// \brief Execute the algorithm using \ref BellmanFord. 436 /// 437 /// Execute the algorithm using the \ref BellmanFord 438 /// "BellmanFord" algorithm for negative cycle detection with 439 /// successively larger limit for the number of iterations. 440 void start() { 441 typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph); 442 typename ResDigraph::template NodeMap<int> visited(*_res_graph); 443 std::vector<ResArc> cycle; 444 int node_num = countNodes(_graph); 622 // Initialize the algorithm 623 ProblemType init() { 624 if (_res_node_num <= 1) return INFEASIBLE; 625 626 // Check the sum of supply values 627 _sum_supply = 0; 628 for (int i = 0; i != _root; ++i) { 629 _sum_supply += _supply[i]; 630 } 631 if (_sum_supply > 0) return INFEASIBLE; 632 633 634 // Initialize vectors 635 for (int i = 0; i != _res_node_num; ++i) { 636 _pi[i] = 0; 637 } 638 ValueVector excess(_supply); 639 640 // Remove infinite upper bounds and check negative arcs 641 const Value MAX = std::numeric_limits<Value>::max(); 642 int last_out; 643 if (_have_lower) { 644 for (int i = 0; i != _root; ++i) { 645 last_out = _first_out[i+1]; 646 for (int j = _first_out[i]; j != last_out; ++j) { 647 if (_forward[j]) { 648 Value c = _cost[j] < 0 ? _upper[j] : _lower[j]; 649 if (c >= MAX) return UNBOUNDED; 650 excess[i] = c; 651 excess[_target[j]] += c; 652 } 653 } 654 } 655 } else { 656 for (int i = 0; i != _root; ++i) { 657 last_out = _first_out[i+1]; 658 for (int j = _first_out[i]; j != last_out; ++j) { 659 if (_forward[j] && _cost[j] < 0) { 660 Value c = _upper[j]; 661 if (c >= MAX) return UNBOUNDED; 662 excess[i] = c; 663 excess[_target[j]] += c; 664 } 665 } 666 } 667 } 668 Value ex, max_cap = 0; 669 for (int i = 0; i != _res_node_num; ++i) { 670 ex = excess[i]; 671 if (ex < 0) max_cap = ex; 672 } 673 for (int j = 0; j != _res_arc_num; ++j) { 674 if (_upper[j] >= MAX) _upper[j] = max_cap; 675 } 676 677 // Initialize maps for Circulation and remove nonzero lower bounds 678 ConstMap<Arc, Value> low(0); 679 typedef typename Digraph::template ArcMap<Value> ValueArcMap; 680 typedef typename Digraph::template NodeMap<Value> ValueNodeMap; 681 ValueArcMap cap(_graph), flow(_graph); 682 ValueNodeMap sup(_graph); 683 for (NodeIt n(_graph); n != INVALID; ++n) { 684 sup[n] = _supply[_node_id[n]]; 685 } 686 if (_have_lower) { 687 for (ArcIt a(_graph); a != INVALID; ++a) { 688 int j = _arc_idf[a]; 689 Value c = _lower[j]; 690 cap[a] = _upper[j]  c; 691 sup[_graph.source(a)] = c; 692 sup[_graph.target(a)] += c; 693 } 694 } else { 695 for (ArcIt a(_graph); a != INVALID; ++a) { 696 cap[a] = _upper[_arc_idf[a]]; 697 } 698 } 699 700 // Find a feasible flow using Circulation 701 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap> 702 circ(_graph, low, cap, sup); 703 if (!circ.flowMap(flow).run()) return INFEASIBLE; 704 705 // Set residual capacities and handle GEQ supply type 706 if (_sum_supply < 0) { 707 for (ArcIt a(_graph); a != INVALID; ++a) { 708 Value fa = flow[a]; 709 _res_cap[_arc_idf[a]] = cap[a]  fa; 710 _res_cap[_arc_idb[a]] = fa; 711 sup[_graph.source(a)] = fa; 712 sup[_graph.target(a)] += fa; 713 } 714 for (NodeIt n(_graph); n != INVALID; ++n) { 715 excess[_node_id[n]] = sup[n]; 716 } 717 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 718 int u = _target[a]; 719 int ra = _reverse[a]; 720 _res_cap[a] = _sum_supply + 1; 721 _res_cap[ra] = excess[u]; 722 _cost[a] = 0; 723 _cost[ra] = 0; 724 } 725 } else { 726 for (ArcIt a(_graph); a != INVALID; ++a) { 727 Value fa = flow[a]; 728 _res_cap[_arc_idf[a]] = cap[a]  fa; 729 _res_cap[_arc_idb[a]] = fa; 730 } 731 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 732 int ra = _reverse[a]; 733 _res_cap[a] = 1; 734 _res_cap[ra] = 0; 735 _cost[a] = 0; 736 _cost[ra] = 0; 737 } 738 } 739 740 return OPTIMAL; 741 } 742 743 // Build a StaticDigraph structure containing the current 744 // residual network 745 void buildResidualNetwork() { 746 _arc_vec.clear(); 747 _cost_vec.clear(); 748 _id_vec.clear(); 749 for (int j = 0; j != _res_arc_num; ++j) { 750 if (_res_cap[j] > 0) { 751 _arc_vec.push_back(IntPair(_source[j], _target[j])); 752 _cost_vec.push_back(_cost[j]); 753 _id_vec.push_back(j); 754 } 755 } 756 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 757 } 758 759 // Execute the algorithm and transform the results 760 void start(Method method) { 761 // Execute the algorithm 762 switch (method) { 763 case SIMPLE_CYCLE_CANCELING: 764 startSimpleCycleCanceling(); 765 break; 766 case MINIMUM_MEAN_CYCLE_CANCELING: 767 startMinMeanCycleCanceling(); 768 break; 769 case CANCEL_AND_TIGHTEN: 770 startCancelAndTighten(); 771 break; 772 } 773 774 // Compute node potentials 775 if (method != SIMPLE_CYCLE_CANCELING) { 776 buildResidualNetwork(); 777 typename BellmanFord<StaticDigraph, CostArcMap> 778 ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map); 779 bf.distMap(_pi_map); 780 bf.init(0); 781 bf.start(); 782 } 783 784 // Handle nonzero lower bounds 785 if (_have_lower) { 786 int limit = _first_out[_root]; 787 for (int j = 0; j != limit; ++j) { 788 if (!_forward[j]) _res_cap[j] += _lower[j]; 789 } 790 } 791 } 792 793 // Execute the "Simple Cycle Canceling" method 794 void startSimpleCycleCanceling() { 795 // Constants for computing the iteration limits 796 const int BF_FIRST_LIMIT = 2; 797 const double BF_LIMIT_FACTOR = 1.5; 798 799 typedef VectorMap<StaticDigraph::Arc, Value> FilterMap; 800 typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph; 801 typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap; 802 typedef typename BellmanFord<ResDigraph, CostArcMap> 803 ::template SetDistMap<CostNodeMap> 804 ::template SetPredMap<PredMap>::Create BF; 805 806 // Build the residual network 807 _arc_vec.clear(); 808 _cost_vec.clear(); 809 for (int j = 0; j != _res_arc_num; ++j) { 810 _arc_vec.push_back(IntPair(_source[j], _target[j])); 811 _cost_vec.push_back(_cost[j]); 812 } 813 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 814 815 FilterMap filter_map(_res_cap); 816 ResDigraph rgr(_sgr, filter_map); 817 std::vector<int> cycle; 818 std::vector<StaticDigraph::Arc> pred(_res_arc_num); 819 PredMap pred_map(pred); 820 BF bf(rgr, _cost_map); 821 bf.distMap(_pi_map).predMap(pred_map); 445 822 446 823 int length_bound = BF_FIRST_LIMIT; 447 824 bool optimal = false; 448 825 while (!optimal) { 449 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);450 bf.predMap(pred);451 826 bf.init(0); 452 827 int iter_num = 0; 453 828 bool cycle_found = false; 454 829 while (!cycle_found) { 455 int curr_iter_num = iter_num + length_bound <= node_num ? 456 length_bound : node_num  iter_num; 830 // Perform some iterations of the BellmanFord algorithm 831 int curr_iter_num = iter_num + length_bound <= _node_num ? 832 length_bound : _node_num  iter_num; 457 833 iter_num += curr_iter_num; 458 834 int real_iter_num = curr_iter_num; … … 466 842 // Optimal flow is found 467 843 optimal = true; 468 // Setting node potentials469 for (NodeIt n(_graph); n != INVALID; ++n)470 (*_potential)[n] = bf.dist(n);471 844 break; 472 845 } else { 473 // Searching for node disjoint negative cycles 474 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) 475 visited[n] = 0; 846 // Search for node disjoint negative cycles 847 std::vector<int> state(_res_node_num, 0); 476 848 int id = 0; 477 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) { 478 if (visited[n] > 0) continue; 479 visited[n] = ++id; 480 ResNode u = pred[n] == INVALID ? 481 INVALID : _res_graph>source(pred[n]); 482 while (u != INVALID && visited[u] == 0) { 483 visited[u] = id; 484 u = pred[u] == INVALID ? 485 INVALID : _res_graph>source(pred[u]); 849 for (int u = 0; u != _res_node_num; ++u) { 850 if (state[u] != 0) continue; 851 ++id; 852 int v = u; 853 for (; v != 1 && state[v] == 0; v = pred[v] == INVALID ? 854 1 : rgr.id(rgr.source(pred[v]))) { 855 state[v] = id; 486 856 } 487 if ( u != INVALID && visited[u] == id) {488 // Finding the negative cycle857 if (v != 1 && state[v] == id) { 858 // A negative cycle is found 489 859 cycle_found = true; 490 860 cycle.clear(); 491 ResArc e = pred[u]; 492 cycle.push_back(e); 493 Capacity d = _res_graph>residualCapacity(e); 494 while (_res_graph>source(e) != u) { 495 cycle.push_back(e = pred[_res_graph>source(e)]); 496 if (_res_graph>residualCapacity(e) < d) 497 d = _res_graph>residualCapacity(e); 861 StaticDigraph::Arc a = pred[v]; 862 Value d, delta = _res_cap[rgr.id(a)]; 863 cycle.push_back(rgr.id(a)); 864 while (rgr.id(rgr.source(a)) != v) { 865 a = pred_map[rgr.source(a)]; 866 d = _res_cap[rgr.id(a)]; 867 if (d < delta) delta = d; 868 cycle.push_back(rgr.id(a)); 498 869 } 499 870 500 // Augmenting along the cycle 501 for (int i = 0; i < int(cycle.size()); ++i) 502 _res_graph>augment(cycle[i], d); 871 // Augment along the cycle 872 for (int i = 0; i < int(cycle.size()); ++i) { 873 int j = cycle[i]; 874 _res_cap[j] = delta; 875 _res_cap[_reverse[j]] += delta; 876 } 503 877 } 504 878 } 505 879 } 506 880 507 if (!cycle_found) 508 length_bound = length_bound * BF_LIMIT_FACTOR / 100; 509 } 510 } 511 } 512 513 /// \brief Execute the algorithm using \ref Howard. 514 /// 515 /// Execute the algorithm using \ref Howard for negative 516 /// cycle detection. 517 void startMinMean() { 518 typedef Path<ResDigraph> ResPath; 519 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost); 520 ResPath cycle; 521 881 // Increase iteration limit if no cycle is found 882 if (!cycle_found) { 883 length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR); 884 } 885 } 886 } 887 } 888 889 // Execute the "Minimum Mean Cycle Canceling" method 890 void startMinMeanCycleCanceling() { 891 typedef SimplePath<StaticDigraph> SPath; 892 typedef typename SPath::ArcIt SPathArcIt; 893 typedef typename Howard<StaticDigraph, CostArcMap> 894 ::template SetPath<SPath>::Create MMC; 895 896 SPath cycle; 897 MMC mmc(_sgr, _cost_map); 522 898 mmc.cycle(cycle); 523 if (mmc.findMinMean()) { 524 while (mmc.cycleLength() < 0) { 525 // Finding the cycle 526 mmc.findCycle(); 527 528 // Finding the largest flow amount that can be augmented 529 // along the cycle 530 Capacity delta = 0; 531 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) { 532 if (delta == 0  _res_graph>residualCapacity(e) < delta) 533 delta = _res_graph>residualCapacity(e); 534 } 535 536 // Augmenting along the cycle 537 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) 538 _res_graph>augment(e, delta); 539 540 // Finding the minimum cycle mean for the modified residual 541 // digraph 542 if (!mmc.findMinMean()) break; 543 } 544 } 545 546 // Computing node potentials 547 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost); 548 bf.init(0); bf.start(); 549 for (NodeIt n(_graph); n != INVALID; ++n) 550 (*_potential)[n] = bf.dist(n); 899 buildResidualNetwork(); 900 while (mmc.findMinMean() && mmc.cycleLength() < 0) { 901 // Find the cycle 902 mmc.findCycle(); 903 904 // Compute delta value 905 Value delta = INF; 906 for (SPathArcIt a(cycle); a != INVALID; ++a) { 907 Value d = _res_cap[_id_vec[_sgr.id(a)]]; 908 if (d < delta) delta = d; 909 } 910 911 // Augment along the cycle 912 for (SPathArcIt a(cycle); a != INVALID; ++a) { 913 int j = _id_vec[_sgr.id(a)]; 914 _res_cap[j] = delta; 915 _res_cap[_reverse[j]] += delta; 916 } 917 918 // Rebuild the residual network 919 buildResidualNetwork(); 920 } 921 } 922 923 // Execute the "Cancel And Tighten" method 924 void startCancelAndTighten() { 925 // Constants for the min mean cycle computations 926 const double LIMIT_FACTOR = 1.0; 927 const int MIN_LIMIT = 5; 928 929 // Contruct auxiliary data vectors 930 DoubleVector pi(_res_node_num, 0.0); 931 IntVector level(_res_node_num); 932 CharVector reached(_res_node_num); 933 CharVector processed(_res_node_num); 934 IntVector pred_node(_res_node_num); 935 IntVector pred_arc(_res_node_num); 936 std::vector<int> stack(_res_node_num); 937 std::vector<int> proc_vector(_res_node_num); 938 939 // Initialize epsilon 940 double epsilon = 0; 941 for (int a = 0; a != _res_arc_num; ++a) { 942 if (_res_cap[a] > 0 && _cost[a] > epsilon) 943 epsilon = _cost[a]; 944 } 945 946 // Start phases 947 Tolerance<double> tol; 948 tol.epsilon(1e6); 949 int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num))); 950 if (limit < MIN_LIMIT) limit = MIN_LIMIT; 951 int iter = limit; 952 while (epsilon * _res_node_num >= 1) { 953 // Find and cancel cycles in the admissible network using DFS 954 for (int u = 0; u != _res_node_num; ++u) { 955 reached[u] = false; 956 processed[u] = false; 957 } 958 int stack_head = 1; 959 int proc_head = 1; 960 for (int start = 0; start != _res_node_num; ++start) { 961 if (reached[start]) continue; 962 963 // New start node 964 reached[start] = true; 965 pred_arc[start] = 1; 966 pred_node[start] = 1; 967 968 // Find the first admissible outgoing arc 969 double p = pi[start]; 970 int a = _first_out[start]; 971 int last_out = _first_out[start+1]; 972 for (; a != last_out && (_res_cap[a] == 0  973 !tol.negative(_cost[a] + p  pi[_target[a]])); ++a) ; 974 if (a == last_out) { 975 processed[start] = true; 976 proc_vector[++proc_head] = start; 977 continue; 978 } 979 stack[++stack_head] = a; 980 981 while (stack_head >= 0) { 982 int sa = stack[stack_head]; 983 int u = _source[sa]; 984 int v = _target[sa]; 985 986 if (!reached[v]) { 987 // A new node is reached 988 reached[v] = true; 989 pred_node[v] = u; 990 pred_arc[v] = sa; 991 p = pi[v]; 992 a = _first_out[v]; 993 last_out = _first_out[v+1]; 994 for (; a != last_out && (_res_cap[a] == 0  995 !tol.negative(_cost[a] + p  pi[_target[a]])); ++a) ; 996 stack[++stack_head] = a == last_out ? 1 : a; 997 } else { 998 if (!processed[v]) { 999 // A cycle is found 1000 int n, w = u; 1001 Value d, delta = _res_cap[sa]; 1002 for (n = u; n != v; n = pred_node[n]) { 1003 d = _res_cap[pred_arc[n]]; 1004 if (d <= delta) { 1005 delta = d; 1006 w = pred_node[n]; 1007 } 1008 } 1009 1010 // Augment along the cycle 1011 _res_cap[sa] = delta; 1012 _res_cap[_reverse[sa]] += delta; 1013 for (n = u; n != v; n = pred_node[n]) { 1014 int pa = pred_arc[n]; 1015 _res_cap[pa] = delta; 1016 _res_cap[_reverse[pa]] += delta; 1017 } 1018 for (n = u; stack_head > 0 && n != w; n = pred_node[n]) { 1019 stack_head; 1020 reached[n] = false; 1021 } 1022 u = w; 1023 } 1024 v = u; 1025 1026 // Find the next admissible outgoing arc 1027 p = pi[v]; 1028 a = stack[stack_head] + 1; 1029 last_out = _first_out[v+1]; 1030 for (; a != last_out && (_res_cap[a] == 0  1031 !tol.negative(_cost[a] + p  pi[_target[a]])); ++a) ; 1032 stack[stack_head] = a == last_out ? 1 : a; 1033 } 1034 1035 while (stack_head >= 0 && stack[stack_head] == 1) { 1036 processed[v] = true; 1037 proc_vector[++proc_head] = v; 1038 if (stack_head >= 0) { 1039 // Find the next admissible outgoing arc 1040 v = _source[stack[stack_head]]; 1041 p = pi[v]; 1042 a = stack[stack_head] + 1; 1043 last_out = _first_out[v+1]; 1044 for (; a != last_out && (_res_cap[a] == 0  1045 !tol.negative(_cost[a] + p  pi[_target[a]])); ++a) ; 1046 stack[stack_head] = a == last_out ? 1 : a; 1047 } 1048 } 1049 } 1050 } 1051 1052 // Tighten potentials and epsilon 1053 if (iter > 0) { 1054 for (int u = 0; u != _res_node_num; ++u) { 1055 level[u] = 0; 1056 } 1057 for (int i = proc_head; i > 0; i) { 1058 int u = proc_vector[i]; 1059 double p = pi[u]; 1060 int l = level[u] + 1; 1061 int last_out = _first_out[u+1]; 1062 for (int a = _first_out[u]; a != last_out; ++a) { 1063 int v = _target[a]; 1064 if (_res_cap[a] > 0 && tol.negative(_cost[a] + p  pi[v]) && 1065 l > level[v]) level[v] = l; 1066 } 1067 } 1068 1069 // Modify potentials 1070 double q = std::numeric_limits<double>::max(); 1071 for (int u = 0; u != _res_node_num; ++u) { 1072 int lu = level[u]; 1073 double p, pu = pi[u]; 1074 int last_out = _first_out[u+1]; 1075 for (int a = _first_out[u]; a != last_out; ++a) { 1076 if (_res_cap[a] == 0) continue; 1077 int v = _target[a]; 1078 int ld = lu  level[v]; 1079 if (ld > 0) { 1080 p = (_cost[a] + pu  pi[v] + epsilon) / (ld + 1); 1081 if (p < q) q = p; 1082 } 1083 } 1084 } 1085 for (int u = 0; u != _res_node_num; ++u) { 1086 pi[u] = q * level[u]; 1087 } 1088 1089 // Modify epsilon 1090 epsilon = 0; 1091 for (int u = 0; u != _res_node_num; ++u) { 1092 double curr, pu = pi[u]; 1093 int last_out = _first_out[u+1]; 1094 for (int a = _first_out[u]; a != last_out; ++a) { 1095 if (_res_cap[a] == 0) continue; 1096 curr = _cost[a] + pu  pi[_target[a]]; 1097 if (curr > epsilon) epsilon = curr; 1098 } 1099 } 1100 } else { 1101 typedef Howard<StaticDigraph, CostArcMap> MMC; 1102 typedef typename BellmanFord<StaticDigraph, CostArcMap> 1103 ::template SetDistMap<CostNodeMap>::Create BF; 1104 1105 // Set epsilon to the minimum cycle mean 1106 buildResidualNetwork(); 1107 MMC mmc(_sgr, _cost_map); 1108 mmc.findMinMean(); 1109 epsilon = mmc.cycleMean(); 1110 Cost cycle_cost = mmc.cycleLength(); 1111 int cycle_size = mmc.cycleArcNum(); 1112 1113 // Compute feasible potentials for the current epsilon 1114 for (int i = 0; i != int(_cost_vec.size()); ++i) { 1115 _cost_vec[i] = cycle_size * _cost_vec[i]  cycle_cost; 1116 } 1117 BF bf(_sgr, _cost_map); 1118 bf.distMap(_pi_map); 1119 bf.init(0); 1120 bf.start(); 1121 for (int u = 0; u != _res_node_num; ++u) { 1122 pi[u] = static_cast<double>(_pi[u]) / cycle_size; 1123 } 1124 1125 iter = limit; 1126 } 1127 } 551 1128 } 552 1129
Note: See TracChangeset
for help on using the changeset viewer.