Changeset 809:22bb98ca0101 in lemon1.2
 Timestamp:
 11/12/09 23:30:45 (10 years ago)
 Branch:
 default
 Phase:
 public
 Rebase:
 35333461653237666330353538396562623465656466303531646265343536366333386130366130
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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