diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h --- a/lemon/cost_scaling.h +++ b/lemon/cost_scaling.h @@ -30,548 +30,912 @@ #include #include #include -#include +#include #include #include namespace lemon { + /// \brief Default traits class of CostScaling algorithm. + /// + /// Default traits class of CostScaling algorithm. + /// \tparam GR Digraph type. + /// \tparam V The value type used for flow amounts, capacity bounds + /// and supply values. By default it is \c int. + /// \tparam C The value type used for costs and potentials. + /// By default it is the same as \c V. +#ifdef DOXYGEN + template +#else + template < typename GR, typename V = int, typename C = V, + bool integer = std::numeric_limits::is_integer > +#endif + struct CostScalingDefaultTraits + { + /// The type of the digraph + typedef GR Digraph; + /// The type of the flow amounts, capacity bounds and supply values + typedef V Value; + /// The type of the arc costs + typedef C Cost; + + /// \brief The large cost type used for internal computations + /// + /// The large cost type used for internal computations. + /// It is \c long \c long if the \c Cost type is integer, + /// otherwise it is \c double. + /// \c Cost must be convertible to \c LargeCost. + typedef double LargeCost; + }; + + // Default traits class for integer cost types + template + struct CostScalingDefaultTraits + { + typedef GR Digraph; + typedef V Value; + typedef C Cost; +#ifdef LEMON_HAVE_LONG_LONG + typedef long long LargeCost; +#else + typedef long LargeCost; +#endif + }; + + /// \addtogroup min_cost_flow_algs /// @{ - /// \brief Implementation of the cost scaling algorithm for finding a - /// minimum cost flow. + /// \brief Implementation of the Cost Scaling algorithm for + /// finding a \ref min_cost_flow "minimum cost flow". /// - /// \ref CostScaling implements the cost scaling algorithm performing - /// augment/push and relabel operations for finding a minimum cost - /// flow. + /// \ref CostScaling implements a cost scaling algorithm that performs + /// push/augment and relabel operations for finding a minimum cost + /// flow. It is an efficient primal-dual solution method, which + /// can be viewed as the generalization of the \ref Preflow + /// "preflow push-relabel" algorithm for the maximum flow problem. /// - /// \tparam Digraph The digraph type the algorithm runs on. - /// \tparam LowerMap The type of the lower bound map. - /// \tparam CapacityMap The type of the capacity (upper bound) map. - /// \tparam CostMap The type of the cost (length) map. - /// \tparam SupplyMap The type of the supply map. + /// Most of the parameters of the problem (except for the digraph) + /// can be given using separate functions, and the algorithm can be + /// executed using the \ref run() function. If some parameters are not + /// specified, then default values will be used. /// - /// \warning - /// - Arc capacities and costs should be \e non-negative \e integers. - /// - Supply values should be \e signed \e integers. - /// - The value types of the maps should be convertible to each other. - /// - \c CostMap::Value must be signed type. + /// \tparam GR The digraph type the algorithm runs on. + /// \tparam V The value type used for flow amounts, capacity bounds + /// and supply values in the algorithm. By default it is \c int. + /// \tparam C The value type used for costs and potentials in the + /// algorithm. By default it is the same as \c V. /// - /// \note Arc costs are multiplied with the number of nodes during - /// the algorithm so overflow problems may arise more easily than with - /// other minimum cost flow algorithms. - /// If it is available, long long int type is used instead of - /// long int in the inside computations. - /// - /// \author Peter Kovacs - template < typename Digraph, - typename LowerMap = typename Digraph::template ArcMap, - typename CapacityMap = typename Digraph::template ArcMap, - typename CostMap = typename Digraph::template ArcMap, - typename SupplyMap = typename Digraph::template NodeMap > + /// \warning Both value types must be signed and all input data must + /// be integer. + /// \warning This algorithm does not support negative costs for such + /// arcs that have infinite upper bound. +#ifdef DOXYGEN + template +#else + template < typename GR, typename V = int, typename C = V, + typename TR = CostScalingDefaultTraits > +#endif class CostScaling { - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); + public: - typedef typename CapacityMap::Value Capacity; - typedef typename CostMap::Value Cost; - typedef typename SupplyMap::Value Supply; - typedef typename Digraph::template ArcMap CapacityArcMap; - typedef typename Digraph::template NodeMap SupplyNodeMap; + /// The type of the digraph + typedef typename TR::Digraph Digraph; + /// The type of the flow amounts, capacity bounds and supply values + typedef typename TR::Value Value; + /// The type of the arc costs + typedef typename TR::Cost Cost; - typedef ResidualDigraph< const Digraph, - CapacityArcMap, CapacityArcMap > ResDigraph; - typedef typename ResDigraph::Arc ResArc; + /// \brief The large cost type + /// + /// The large cost type used for internal computations. + /// Using the \ref CostScalingDefaultTraits "default traits class", + /// it is \c long \c long if the \c Cost type is integer, + /// otherwise it is \c double. + typedef typename TR::LargeCost LargeCost; -#if defined __GNUC__ && !defined __STRICT_ANSI__ - typedef long long int LCost; -#else - typedef long int LCost; -#endif - typedef typename Digraph::template ArcMap LargeCostMap; + /// The \ref CostScalingDefaultTraits "traits class" of the algorithm + typedef TR Traits; public: - /// The type of the flow map. - typedef typename Digraph::template ArcMap FlowMap; - /// The type of the potential map. - typedef typename Digraph::template NodeMap PotentialMap; + /// \brief Problem type constants for the \c run() function. + /// + /// Enum type containing the problem type constants that can be + /// returned by the \ref run() function of the algorithm. + enum ProblemType { + /// The problem has no feasible solution (flow). + INFEASIBLE, + /// The problem has optimal solution (i.e. it is feasible and + /// bounded), and the algorithm has found optimal flow and node + /// potentials (primal and dual solutions). + OPTIMAL, + /// The digraph contains an arc of negative cost and infinite + /// upper bound. It means that the objective function is unbounded + /// on that arc, however note that it could actually be bounded + /// over the feasible flows, but this algroithm cannot handle + /// these cases. + UNBOUNDED + }; private: - /// \brief Map adaptor class for handling residual arc costs. - /// - /// Map adaptor class for handling residual arc costs. - template - class ResidualCostMap : public MapBase - { - private: + TEMPLATE_DIGRAPH_TYPEDEFS(GR); - const Map &_cost_map; + typedef std::vector IntVector; + typedef std::vector BoolVector; + typedef std::vector ValueVector; + typedef std::vector CostVector; + typedef std::vector LargeCostVector; + private: + + template + class VectorMap { public: - - ///\e - ResidualCostMap(const Map &cost_map) : - _cost_map(cost_map) {} - - ///\e - inline typename Map::Value operator[](const ResArc &e) const { - return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e]; + typedef KT Key; + typedef VT Value; + + VectorMap(std::vector& v) : _v(v) {} + + const Value& operator[](const Key& key) const { + return _v[StaticDigraph::id(key)]; } - }; //class ResidualCostMap - - /// \brief Map adaptor class for handling reduced arc costs. - /// - /// Map adaptor class for handling reduced arc costs. - class ReducedCostMap : public MapBase - { - private: - - const Digraph &_gr; - const LargeCostMap &_cost_map; - const PotentialMap &_pot_map; - - public: - - ///\e - ReducedCostMap( const Digraph &gr, - const LargeCostMap &cost_map, - const PotentialMap &pot_map ) : - _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} - - ///\e - inline LCost operator[](const Arc &e) const { - return _cost_map[e] + _pot_map[_gr.source(e)] - - _pot_map[_gr.target(e)]; + Value& operator[](const Key& key) { + return _v[StaticDigraph::id(key)]; + } + + void set(const Key& key, const Value& val) { + _v[StaticDigraph::id(key)] = val; } - }; //class ReducedCostMap + private: + std::vector& _v; + }; + + typedef VectorMap LargeCostNodeMap; + typedef VectorMap LargeCostArcMap; private: - // The digraph the algorithm runs on - const Digraph &_graph; - // The original lower bound map - const LowerMap *_lower; - // The modified capacity map - CapacityArcMap _capacity; - // The original cost map - const CostMap &_orig_cost; - // The scaled cost map - LargeCostMap _cost; - // The modified supply map - SupplyNodeMap _supply; - bool _valid_supply; + // Data related to the underlying digraph + const GR &_graph; + int _node_num; + int _arc_num; + int _res_node_num; + int _res_arc_num; + int _root; - // Arc map of the current flow - FlowMap *_flow; - bool _local_flow; - // Node map of the current potentials - PotentialMap *_potential; - bool _local_potential; + // Parameters of the problem + bool _have_lower; + Value _sum_supply; - // The residual cost map - ResidualCostMap _res_cost; - // The residual digraph - ResDigraph *_res_graph; - // The reduced cost map - ReducedCostMap *_red_cost; - // The excess map - SupplyNodeMap _excess; - // The epsilon parameter used for cost scaling - LCost _epsilon; - // The scaling factor + // Data structures for storing the digraph + IntNodeMap _node_id; + IntArcMap _arc_idf; + IntArcMap _arc_idb; + IntVector _first_out; + BoolVector _forward; + IntVector _source; + IntVector _target; + IntVector _reverse; + + // Node and arc data + ValueVector _lower; + ValueVector _upper; + CostVector _scost; + ValueVector _supply; + + ValueVector _res_cap; + LargeCostVector _cost; + LargeCostVector _pi; + ValueVector _excess; + IntVector _next_out; + std::deque _active_nodes; + + // Data for scaling + LargeCost _epsilon; int _alpha; + // Data for a StaticDigraph structure + typedef std::pair IntPair; + StaticDigraph _sgr; + std::vector _arc_vec; + std::vector _cost_vec; + LargeCostArcMap _cost_map; + LargeCostNodeMap _pi_map; + + public: + + /// \brief Constant for infinite upper bounds (capacities). + /// + /// Constant for infinite upper bounds (capacities). + /// It is \c std::numeric_limits::infinity() if available, + /// \c std::numeric_limits::max() otherwise. + const Value INF; + public: - /// \brief General constructor (with lower bounds). + /// \name Named Template Parameters + /// @{ + + template + struct SetLargeCostTraits : public Traits { + typedef T LargeCost; + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// \c LargeCost type. /// - /// General constructor (with lower bounds). + /// \ref named-templ-param "Named parameter" for setting \c LargeCost + /// type, which is used for internal computations in the algorithm. + /// \c Cost must be convertible to \c LargeCost. + template + struct SetLargeCost + : public CostScaling > { + typedef CostScaling > Create; + }; + + /// @} + + public: + + /// \brief Constructor. /// - /// \param digraph The digraph the algorithm runs on. - /// \param lower The lower bounds of the arcs. - /// \param capacity The capacities (upper bounds) of the arcs. - /// \param cost The cost (length) values of the arcs. - /// \param supply The supply values of the nodes (signed). - CostScaling( const Digraph &digraph, - const LowerMap &lower, - const CapacityMap &capacity, - const CostMap &cost, - const SupplyMap &supply ) : - _graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost), - _cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false), - _potential(NULL), _local_potential(false), _res_cost(_cost), - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) + /// The constructor of the class. + /// + /// \param graph The digraph the algorithm runs on. + CostScaling(const GR& graph) : + _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), + _cost_map(_cost_vec), _pi_map(_pi), + INF(std::numeric_limits::has_infinity ? + std::numeric_limits::infinity() : + std::numeric_limits::max()) { - // Check the sum of supply values - Supply sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; - _valid_supply = sum == 0; + // Check the value types + LEMON_ASSERT(std::numeric_limits::is_signed, + "The flow type of CostScaling must be signed"); + LEMON_ASSERT(std::numeric_limits::is_signed, + "The cost type of CostScaling must be signed"); + + // Resize vectors + _node_num = countNodes(_graph); + _arc_num = countArcs(_graph); + _res_node_num = _node_num + 1; + _res_arc_num = 2 * (_arc_num + _node_num); + _root = _node_num; + + _first_out.resize(_res_node_num + 1); + _forward.resize(_res_arc_num); + _source.resize(_res_arc_num); + _target.resize(_res_arc_num); + _reverse.resize(_res_arc_num); + + _lower.resize(_res_arc_num); + _upper.resize(_res_arc_num); + _scost.resize(_res_arc_num); + _supply.resize(_res_node_num); - for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e]; - for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n]; + _res_cap.resize(_res_arc_num); + _cost.resize(_res_arc_num); + _pi.resize(_res_node_num); + _excess.resize(_res_node_num); + _next_out.resize(_res_node_num); - // Remove non-zero lower bounds - for (ArcIt e(_graph); e != INVALID; ++e) { - if (lower[e] != 0) { - _capacity[e] -= lower[e]; - _supply[_graph.source(e)] -= lower[e]; - _supply[_graph.target(e)] += lower[e]; + _arc_vec.reserve(_res_arc_num); + _cost_vec.reserve(_res_arc_num); + + // Copy the graph + int i = 0, j = 0, k = 2 * _arc_num + _node_num; + for (NodeIt n(_graph); n != INVALID; ++n, ++i) { + _node_id[n] = i; + } + i = 0; + for (NodeIt n(_graph); n != INVALID; ++n, ++i) { + _first_out[i] = j; + for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { + _arc_idf[a] = j; + _forward[j] = true; + _source[j] = i; + _target[j] = _node_id[_graph.runningNode(a)]; } + for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { + _arc_idb[a] = j; + _forward[j] = false; + _source[j] = i; + _target[j] = _node_id[_graph.runningNode(a)]; + } + _forward[j] = false; + _source[j] = i; + _target[j] = _root; + _reverse[j] = k; + _forward[k] = true; + _source[k] = _root; + _target[k] = i; + _reverse[k] = j; + ++j; ++k; } - } -/* - /// \brief General constructor (without lower bounds). - /// - /// General constructor (without lower bounds). - /// - /// \param digraph The digraph the algorithm runs on. - /// \param capacity The capacities (upper bounds) of the arcs. - /// \param cost The cost (length) values of the arcs. - /// \param supply The supply values of the nodes (signed). - CostScaling( const Digraph &digraph, - const CapacityMap &capacity, - const CostMap &cost, - const SupplyMap &supply ) : - _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost), - _cost(digraph), _supply(supply), _flow(NULL), _local_flow(false), - _potential(NULL), _local_potential(false), _res_cost(_cost), - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) - { - // Check the sum of supply values - Supply sum = 0; - for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; - _valid_supply = sum == 0; + _first_out[i] = j; + _first_out[_res_node_num] = k; + for (ArcIt a(_graph); a != INVALID; ++a) { + int fi = _arc_idf[a]; + int bi = _arc_idb[a]; + _reverse[fi] = bi; + _reverse[bi] = fi; + } + + // Reset parameters + reset(); } - /// \brief Simple constructor (with lower bounds). + /// \name Parameters + /// The parameters of the algorithm can be specified using these + /// functions. + + /// @{ + + /// \brief Set the lower bounds on the arcs. /// - /// Simple constructor (with lower bounds). + /// This function sets the lower bounds on the arcs. + /// If it is not used before calling \ref run(), the lower bounds + /// will be set to zero on all arcs. /// - /// \param digraph The digraph the algorithm runs on. - /// \param lower The lower bounds of the arcs. - /// \param capacity The capacities (upper bounds) of the arcs. - /// \param cost The cost (length) values of the arcs. - /// \param s The source node. - /// \param t The target node. - /// \param flow_value The required amount of flow from node \c s - /// to node \c t (i.e. the supply of \c s and the demand of \c t). - CostScaling( const Digraph &digraph, - const LowerMap &lower, - const CapacityMap &capacity, - const CostMap &cost, - Node s, Node t, - Supply flow_value ) : - _graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost), - _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false), - _potential(NULL), _local_potential(false), _res_cost(_cost), - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) - { - // Remove non-zero lower bounds - _supply[s] = flow_value; - _supply[t] = -flow_value; - for (ArcIt e(_graph); e != INVALID; ++e) { - if (lower[e] != 0) { - _capacity[e] -= lower[e]; - _supply[_graph.source(e)] -= lower[e]; - _supply[_graph.target(e)] += lower[e]; - } + /// \param map An arc map storing the lower bounds. + /// Its \c Value type must be convertible to the \c Value type + /// of the algorithm. + /// + /// \return (*this) + template + CostScaling& lowerMap(const LowerMap& map) { + _have_lower = true; + for (ArcIt a(_graph); a != INVALID; ++a) { + _lower[_arc_idf[a]] = map[a]; + _lower[_arc_idb[a]] = map[a]; } - _valid_supply = true; - } - - /// \brief Simple constructor (without lower bounds). - /// - /// Simple constructor (without lower bounds). - /// - /// \param digraph The digraph the algorithm runs on. - /// \param capacity The capacities (upper bounds) of the arcs. - /// \param cost The cost (length) values of the arcs. - /// \param s The source node. - /// \param t The target node. - /// \param flow_value The required amount of flow from node \c s - /// to node \c t (i.e. the supply of \c s and the demand of \c t). - CostScaling( const Digraph &digraph, - const CapacityMap &capacity, - const CostMap &cost, - Node s, Node t, - Supply flow_value ) : - _graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost), - _cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false), - _potential(NULL), _local_potential(false), _res_cost(_cost), - _res_graph(NULL), _red_cost(NULL), _excess(digraph, 0) - { - _supply[s] = flow_value; - _supply[t] = -flow_value; - _valid_supply = true; - } -*/ - /// Destructor. - ~CostScaling() { - if (_local_flow) delete _flow; - if (_local_potential) delete _potential; - delete _res_graph; - delete _red_cost; - } - - /// \brief Set the flow map. - /// - /// Set the flow map. - /// - /// \return \c (*this) - CostScaling& flowMap(FlowMap &map) { - if (_local_flow) { - delete _flow; - _local_flow = false; - } - _flow = ↦ return *this; } - /// \brief Set the potential map. + /// \brief Set the upper bounds (capacities) on the arcs. /// - /// Set the potential map. + /// This function sets the upper bounds (capacities) on the arcs. + /// If it is not used before calling \ref run(), the upper bounds + /// will be set to \ref INF on all arcs (i.e. the flow value will be + /// unbounded from above on each arc). /// - /// \return \c (*this) - CostScaling& potentialMap(PotentialMap &map) { - if (_local_potential) { - delete _potential; - _local_potential = false; + /// \param map An arc map storing the upper bounds. + /// Its \c Value type must be convertible to the \c Value type + /// of the algorithm. + /// + /// \return (*this) + template + CostScaling& upperMap(const UpperMap& map) { + for (ArcIt a(_graph); a != INVALID; ++a) { + _upper[_arc_idf[a]] = map[a]; } - _potential = ↦ return *this; } + /// \brief Set the costs of the arcs. + /// + /// This function sets the costs of the arcs. + /// If it is not used before calling \ref run(), the costs + /// will be set to \c 1 on all arcs. + /// + /// \param map An arc map storing the costs. + /// Its \c Value type must be convertible to the \c Cost type + /// of the algorithm. + /// + /// \return (*this) + template + CostScaling& costMap(const CostMap& map) { + for (ArcIt a(_graph); a != INVALID; ++a) { + _scost[_arc_idf[a]] = map[a]; + _scost[_arc_idb[a]] = -map[a]; + } + return *this; + } + + /// \brief Set the supply values of the nodes. + /// + /// This function sets the supply values of the nodes. + /// If neither this function nor \ref stSupply() is used before + /// calling \ref run(), the supply of each node will be set to zero. + /// + /// \param map A node map storing the supply values. + /// Its \c Value type must be convertible to the \c Value type + /// of the algorithm. + /// + /// \return (*this) + template + CostScaling& supplyMap(const SupplyMap& map) { + for (NodeIt n(_graph); n != INVALID; ++n) { + _supply[_node_id[n]] = map[n]; + } + return *this; + } + + /// \brief Set single source and target nodes and a supply value. + /// + /// This function sets a single source node and a single target node + /// and the required flow value. + /// If neither this function nor \ref supplyMap() is used before + /// calling \ref run(), the supply of each node will be set to zero. + /// + /// Using this function has the same effect as using \ref supplyMap() + /// with such a map in which \c k is assigned to \c s, \c -k is + /// assigned to \c t and all other nodes have zero supply value. + /// + /// \param s The source node. + /// \param t The target node. + /// \param k The required amount of flow from node \c s to node \c t + /// (i.e. the supply of \c s and the demand of \c t). + /// + /// \return (*this) + CostScaling& stSupply(const Node& s, const Node& t, Value k) { + for (int i = 0; i != _res_node_num; ++i) { + _supply[i] = 0; + } + _supply[_node_id[s]] = k; + _supply[_node_id[t]] = -k; + return *this; + } + + /// @} + /// \name Execution control + /// The algorithm can be executed using \ref run(). /// @{ /// \brief Run the algorithm. /// - /// Run the algorithm. + /// This function runs the algorithm. + /// The paramters can be specified using functions \ref lowerMap(), + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). + /// For example, + /// \code + /// CostScaling cs(graph); + /// cs.lowerMap(lower).upperMap(upper).costMap(cost) + /// .supplyMap(sup).run(); + /// \endcode + /// + /// This function can be called more than once. All the parameters + /// that have been given are kept for the next call, unless + /// \ref reset() is called, thus only the modified parameters + /// have to be set again. See \ref reset() for examples. + /// However the underlying digraph must not be modified after this + /// class have been constructed, since it copies the digraph. /// /// \param partial_augment By default the algorithm performs /// partial augment and relabel operations in the cost scaling /// phases. Set this parameter to \c false for using local push and /// relabel operations instead. /// - /// \return \c true if a feasible flow can be found. - bool run(bool partial_augment = true) { - if (partial_augment) { - return init() && startPartialAugment(); - } else { - return init() && startPushRelabel(); + /// \return \c INFEASIBLE if no feasible flow exists, + /// \n \c OPTIMAL if the problem has optimal solution + /// (i.e. it is feasible and bounded), and the algorithm has found + /// optimal flow and node potentials (primal and dual solutions), + /// \n \c UNBOUNDED if the digraph contains an arc of negative cost + /// and infinite upper bound. It means that the objective function + /// is unbounded on that arc, however note that it could actually be + /// bounded over the feasible flows, but this algroithm cannot handle + /// these cases. + /// + /// \see ProblemType + ProblemType run(bool partial_augment = true) { + ProblemType pt = init(); + if (pt != OPTIMAL) return pt; + start(partial_augment); + return OPTIMAL; + } + + /// \brief Reset all the parameters that have been given before. + /// + /// This function resets all the paramaters that have been given + /// before using functions \ref lowerMap(), \ref upperMap(), + /// \ref costMap(), \ref supplyMap(), \ref stSupply(). + /// + /// It is useful for multiple run() calls. If this function is not + /// used, all the parameters given before are kept for the next + /// \ref run() call. + /// However the underlying digraph must not be modified after this + /// class have been constructed, since it copies and extends the graph. + /// + /// For example, + /// \code + /// CostScaling cs(graph); + /// + /// // First run + /// cs.lowerMap(lower).upperMap(upper).costMap(cost) + /// .supplyMap(sup).run(); + /// + /// // Run again with modified cost map (reset() is not called, + /// // so only the cost map have to be set again) + /// cost[e] += 100; + /// cs.costMap(cost).run(); + /// + /// // Run again from scratch using reset() + /// // (the lower bounds will be set to zero on all arcs) + /// cs.reset(); + /// cs.upperMap(capacity).costMap(cost) + /// .supplyMap(sup).run(); + /// \endcode + /// + /// \return (*this) + CostScaling& reset() { + for (int i = 0; i != _res_node_num; ++i) { + _supply[i] = 0; } + int limit = _first_out[_root]; + for (int j = 0; j != limit; ++j) { + _lower[j] = 0; + _upper[j] = INF; + _scost[j] = _forward[j] ? 1 : -1; + } + for (int j = limit; j != _res_arc_num; ++j) { + _lower[j] = 0; + _upper[j] = INF; + _scost[j] = 0; + _scost[_reverse[j]] = 0; + } + _have_lower = false; + return *this; } /// @} /// \name Query Functions - /// The result of the algorithm can be obtained using these + /// The results of the algorithm can be obtained using these /// functions.\n - /// \ref lemon::CostScaling::run() "run()" must be called before - /// using them. + /// The \ref run() function must be called before using them. /// @{ - /// \brief Return a const reference to the arc map storing the - /// found flow. + /// \brief Return the total cost of the found flow. /// - /// Return a const reference to the arc map storing the found flow. + /// This function returns the total cost of the found flow. + /// Its complexity is O(e). + /// + /// \note The return type of the function can be specified as a + /// template parameter. For example, + /// \code + /// cs.totalCost(); + /// \endcode + /// It is useful if the total cost cannot be stored in the \c Cost + /// type of the algorithm, which is the default return type of the + /// function. /// /// \pre \ref run() must be called before using this function. - const FlowMap& flowMap() const { - return *_flow; + template + Number totalCost() const { + Number c = 0; + for (ArcIt a(_graph); a != INVALID; ++a) { + int i = _arc_idb[a]; + c += static_cast(_res_cap[i]) * + (-static_cast(_scost[i])); + } + return c; } - /// \brief Return a const reference to the node map storing the - /// found potentials (the dual solution). - /// - /// Return a const reference to the node map storing the found - /// potentials (the dual solution). - /// - /// \pre \ref run() must be called before using this function. - const PotentialMap& potentialMap() const { - return *_potential; +#ifndef DOXYGEN + Cost totalCost() const { + return totalCost(); } +#endif /// \brief Return the flow on the given arc. /// - /// Return the flow on the given arc. + /// This function returns the flow on the given arc. /// /// \pre \ref run() must be called before using this function. - Capacity flow(const Arc& arc) const { - return (*_flow)[arc]; + Value flow(const Arc& a) const { + return _res_cap[_arc_idb[a]]; } - /// \brief Return the potential of the given node. + /// \brief Return the flow map (the primal solution). /// - /// Return the potential of the given node. + /// This function copies the flow value on each arc into the given + /// map. The \c Value type of the algorithm must be convertible to + /// the \c Value type of the map. /// /// \pre \ref run() must be called before using this function. - Cost potential(const Node& node) const { - return (*_potential)[node]; + template + void flowMap(FlowMap &map) const { + for (ArcIt a(_graph); a != INVALID; ++a) { + map.set(a, _res_cap[_arc_idb[a]]); + } } - /// \brief Return the total cost of the found flow. + /// \brief Return the potential (dual value) of the given node. /// - /// Return the total cost of the found flow. The complexity of the - /// function is \f$ O(e) \f$. + /// This function returns the potential (dual value) of the + /// given node. /// /// \pre \ref run() must be called before using this function. - Cost totalCost() const { - Cost c = 0; - for (ArcIt e(_graph); e != INVALID; ++e) - c += (*_flow)[e] * _orig_cost[e]; - return c; + Cost potential(const Node& n) const { + return static_cast(_pi[_node_id[n]]); + } + + /// \brief Return the potential map (the dual solution). + /// + /// This function copies the potential (dual value) of each node + /// into the given map. + /// The \c Cost type of the algorithm must be convertible to the + /// \c Value type of the map. + /// + /// \pre \ref run() must be called before using this function. + template + void potentialMap(PotentialMap &map) const { + for (NodeIt n(_graph); n != INVALID; ++n) { + map.set(n, static_cast(_pi[_node_id[n]])); + } } /// @} private: - /// Initialize the algorithm. - bool init() { - if (!_valid_supply) return false; - // The scaling factor + // Initialize the algorithm + ProblemType init() { + if (_res_node_num == 0) return INFEASIBLE; + + // Scaling factor _alpha = 8; - // Initialize flow and potential maps - if (!_flow) { - _flow = new FlowMap(_graph); - _local_flow = true; + // Check the sum of supply values + _sum_supply = 0; + for (int i = 0; i != _root; ++i) { + _sum_supply += _supply[i]; } - if (!_potential) { - _potential = new PotentialMap(_graph); - _local_potential = true; + if (_sum_supply > 0) return INFEASIBLE; + + + // Initialize vectors + for (int i = 0; i != _res_node_num; ++i) { + _pi[i] = 0; + _excess[i] = _supply[i]; + } + + // Remove infinite upper bounds and check negative arcs + const Value MAX = std::numeric_limits::max(); + int last_out; + if (_have_lower) { + for (int i = 0; i != _root; ++i) { + last_out = _first_out[i+1]; + for (int j = _first_out[i]; j != last_out; ++j) { + if (_forward[j]) { + Value c = _scost[j] < 0 ? _upper[j] : _lower[j]; + if (c >= MAX) return UNBOUNDED; + _excess[i] -= c; + _excess[_target[j]] += c; + } + } + } + } else { + for (int i = 0; i != _root; ++i) { + last_out = _first_out[i+1]; + for (int j = _first_out[i]; j != last_out; ++j) { + if (_forward[j] && _scost[j] < 0) { + Value c = _upper[j]; + if (c >= MAX) return UNBOUNDED; + _excess[i] -= c; + _excess[_target[j]] += c; + } + } + } + } + Value ex, max_cap = 0; + for (int i = 0; i != _res_node_num; ++i) { + ex = _excess[i]; + _excess[i] = 0; + if (ex < 0) max_cap -= ex; + } + for (int j = 0; j != _res_arc_num; ++j) { + if (_upper[j] >= MAX) _upper[j] = max_cap; } - _red_cost = new ReducedCostMap(_graph, _cost, *_potential); - _res_graph = new ResDigraph(_graph, _capacity, *_flow); + // Initialize the large cost vector and the epsilon parameter + _epsilon = 0; + LargeCost lc; + for (int i = 0; i != _root; ++i) { + last_out = _first_out[i+1]; + for (int j = _first_out[i]; j != last_out; ++j) { + lc = static_cast(_scost[j]) * _res_node_num * _alpha; + _cost[j] = lc; + if (lc > _epsilon) _epsilon = lc; + } + } + _epsilon /= _alpha; - // Initialize the scaled cost map and the epsilon parameter - Cost max_cost = 0; - int node_num = countNodes(_graph); - for (ArcIt e(_graph); e != INVALID; ++e) { - _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha; - if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e]; + // Initialize maps for Circulation and remove non-zero lower bounds + ConstMap low(0); + typedef typename Digraph::template ArcMap ValueArcMap; + typedef typename Digraph::template NodeMap ValueNodeMap; + ValueArcMap cap(_graph), flow(_graph); + ValueNodeMap sup(_graph); + for (NodeIt n(_graph); n != INVALID; ++n) { + sup[n] = _supply[_node_id[n]]; } - _epsilon = max_cost * node_num; + if (_have_lower) { + for (ArcIt a(_graph); a != INVALID; ++a) { + int j = _arc_idf[a]; + Value c = _lower[j]; + cap[a] = _upper[j] - c; + sup[_graph.source(a)] -= c; + sup[_graph.target(a)] += c; + } + } else { + for (ArcIt a(_graph); a != INVALID; ++a) { + cap[a] = _upper[_arc_idf[a]]; + } + } // Find a feasible flow using Circulation - Circulation< Digraph, ConstMap, CapacityArcMap, - SupplyMap > - circulation( _graph, constMap(Capacity(0)), _capacity, - _supply ); - return circulation.flowMap(*_flow).run(); + Circulation, ValueArcMap, ValueNodeMap> + circ(_graph, low, cap, sup); + if (!circ.flowMap(flow).run()) return INFEASIBLE; + + // Set residual capacities and handle GEQ supply type + if (_sum_supply < 0) { + for (ArcIt a(_graph); a != INVALID; ++a) { + Value fa = flow[a]; + _res_cap[_arc_idf[a]] = cap[a] - fa; + _res_cap[_arc_idb[a]] = fa; + sup[_graph.source(a)] -= fa; + sup[_graph.target(a)] += fa; + } + for (NodeIt n(_graph); n != INVALID; ++n) { + _excess[_node_id[n]] = sup[n]; + } + for (int a = _first_out[_root]; a != _res_arc_num; ++a) { + int u = _target[a]; + int ra = _reverse[a]; + _res_cap[a] = -_sum_supply + 1; + _res_cap[ra] = -_excess[u]; + _cost[a] = 0; + _cost[ra] = 0; + _excess[u] = 0; + } + } else { + for (ArcIt a(_graph); a != INVALID; ++a) { + Value fa = flow[a]; + _res_cap[_arc_idf[a]] = cap[a] - fa; + _res_cap[_arc_idb[a]] = fa; + } + for (int a = _first_out[_root]; a != _res_arc_num; ++a) { + int ra = _reverse[a]; + _res_cap[a] = 1; + _res_cap[ra] = 0; + _cost[a] = 0; + _cost[ra] = 0; + } + } + + return OPTIMAL; + } + + // Execute the algorithm and transform the results + void start(bool partial_augment) { + // Execute the algorithm + if (partial_augment) { + startPartialAugment(); + } else { + startPushRelabel(); + } + + // Compute node potentials for the original costs + _arc_vec.clear(); + _cost_vec.clear(); + for (int j = 0; j != _res_arc_num; ++j) { + if (_res_cap[j] > 0) { + _arc_vec.push_back(IntPair(_source[j], _target[j])); + _cost_vec.push_back(_scost[j]); + } + } + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); + + typename BellmanFord + ::template SetDistMap::Create bf(_sgr, _cost_map); + bf.distMap(_pi_map); + bf.init(0); + bf.start(); + + // Handle non-zero lower bounds + if (_have_lower) { + int limit = _first_out[_root]; + for (int j = 0; j != limit; ++j) { + if (!_forward[j]) _res_cap[j] += _lower[j]; + } + } } /// Execute the algorithm performing partial augmentation and - /// relabel operations. - bool startPartialAugment() { + /// relabel operations + void startPartialAugment() { // Paramters for heuristics -// const int BF_HEURISTIC_EPSILON_BOUND = 1000; -// const int BF_HEURISTIC_BOUND_FACTOR = 3; + const int BF_HEURISTIC_EPSILON_BOUND = 1000; + const int BF_HEURISTIC_BOUND_FACTOR = 3; // Maximum augment path length const int MAX_PATH_LENGTH = 4; - // Variables - typename Digraph::template NodeMap pred_arc(_graph); - typename Digraph::template NodeMap forward(_graph); - typename Digraph::template NodeMap next_out(_graph); - typename Digraph::template NodeMap next_in(_graph); - typename Digraph::template NodeMap next_dir(_graph); - std::deque active_nodes; - std::vector path_nodes; - -// int node_num = countNodes(_graph); + // Perform cost scaling phases + IntVector pred_arc(_res_node_num); + std::vector path_nodes; for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1 : _epsilon / _alpha ) { -/* // "Early Termination" heuristic: use Bellman-Ford algorithm // to check if the current flow is optimal if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { - typedef ShiftMap< ResidualCostMap > ShiftCostMap; - ShiftCostMap shift_cost(_res_cost, 1); - BellmanFord bf(*_res_graph, shift_cost); + _arc_vec.clear(); + _cost_vec.clear(); + for (int j = 0; j != _res_arc_num; ++j) { + if (_res_cap[j] > 0) { + _arc_vec.push_back(IntPair(_source[j], _target[j])); + _cost_vec.push_back(_cost[j] + 1); + } + } + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); + + BellmanFord bf(_sgr, _cost_map); bf.init(0); bool done = false; - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); for (int i = 0; i < K && !done; ++i) done = bf.processNextWeakRound(); if (done) break; } -*/ + // Saturate arcs not satisfying the optimality condition - Capacity delta; - for (ArcIt e(_graph); e != INVALID; ++e) { - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { - delta = _capacity[e] - (*_flow)[e]; - _excess[_graph.source(e)] -= delta; - _excess[_graph.target(e)] += delta; - (*_flow)[e] = _capacity[e]; - } - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { - _excess[_graph.target(e)] -= (*_flow)[e]; - _excess[_graph.source(e)] += (*_flow)[e]; - (*_flow)[e] = 0; + for (int a = 0; a != _res_arc_num; ++a) { + if (_res_cap[a] > 0 && + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { + Value delta = _res_cap[a]; + _excess[_source[a]] -= delta; + _excess[_target[a]] += delta; + _res_cap[a] = 0; + _res_cap[_reverse[a]] += delta; } } - + // Find active nodes (i.e. nodes with positive excess) - for (NodeIt n(_graph); n != INVALID; ++n) { - if (_excess[n] > 0) active_nodes.push_back(n); + for (int u = 0; u != _res_node_num; ++u) { + if (_excess[u] > 0) _active_nodes.push_back(u); } - // Initialize the next arc maps - for (NodeIt n(_graph); n != INVALID; ++n) { - next_out[n] = OutArcIt(_graph, n); - next_in[n] = InArcIt(_graph, n); - next_dir[n] = true; + // Initialize the next arcs + for (int u = 0; u != _res_node_num; ++u) { + _next_out[u] = _first_out[u]; } // Perform partial augment and relabel operations - while (active_nodes.size() > 0) { + while (true) { // Select an active node (FIFO selection) - if (_excess[active_nodes[0]] <= 0) { - active_nodes.pop_front(); - continue; + while (_active_nodes.size() > 0 && + _excess[_active_nodes.front()] <= 0) { + _active_nodes.pop_front(); } - Node start = active_nodes[0]; + if (_active_nodes.size() == 0) break; + int start = _active_nodes.front(); path_nodes.clear(); path_nodes.push_back(start); // Find an augmenting path from the start node - Node u, tip = start; - LCost min_red_cost; - while ( _excess[tip] >= 0 && - int(path_nodes.size()) <= MAX_PATH_LENGTH ) - { - if (next_dir[tip]) { - for (OutArcIt e = next_out[tip]; e != INVALID; ++e) { - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { - u = _graph.target(e); - pred_arc[u] = e; - forward[u] = true; - next_out[tip] = e; - tip = u; - path_nodes.push_back(tip); - goto next_step; - } - } - next_dir[tip] = false; - } - for (InArcIt e = next_in[tip]; e != INVALID; ++e) { - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { - u = _graph.source(e); - pred_arc[u] = e; - forward[u] = false; - next_in[tip] = e; + int tip = start; + while (_excess[tip] >= 0 && + int(path_nodes.size()) <= MAX_PATH_LENGTH) { + int u; + LargeCost min_red_cost, rc; + int last_out = _sum_supply < 0 ? + _first_out[tip+1] : _first_out[tip+1] - 1; + for (int a = _next_out[tip]; a != last_out; ++a) { + if (_res_cap[a] > 0 && + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { + u = _target[a]; + pred_arc[u] = a; + _next_out[tip] = a; tip = u; path_nodes.push_back(tip); goto next_step; @@ -579,266 +943,186 @@ } // Relabel tip node - min_red_cost = std::numeric_limits::max() / 2; - for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) { - if ( _capacity[oe] - (*_flow)[oe] > 0 && - (*_red_cost)[oe] < min_red_cost ) - min_red_cost = (*_red_cost)[oe]; + min_red_cost = std::numeric_limits::max() / 2; + for (int a = _first_out[tip]; a != last_out; ++a) { + rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; + if (_res_cap[a] > 0 && rc < min_red_cost) { + min_red_cost = rc; + } } - for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) { - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) - min_red_cost = -(*_red_cost)[ie]; - } - (*_potential)[tip] -= min_red_cost + _epsilon; + _pi[tip] -= min_red_cost + _epsilon; - // Reset the next arc maps - next_out[tip] = OutArcIt(_graph, tip); - next_in[tip] = InArcIt(_graph, tip); - next_dir[tip] = true; + // Reset the next arc of tip + _next_out[tip] = _first_out[tip]; // Step back if (tip != start) { path_nodes.pop_back(); - tip = path_nodes[path_nodes.size()-1]; + tip = path_nodes.back(); } - next_step: - continue; + next_step: ; } // Augment along the found path (as much flow as possible) - Capacity delta; + Value delta; + int u, v = path_nodes.front(), pa; for (int i = 1; i < int(path_nodes.size()); ++i) { - u = path_nodes[i]; - delta = forward[u] ? - _capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] : - (*_flow)[pred_arc[u]]; - delta = std::min(delta, _excess[path_nodes[i-1]]); - (*_flow)[pred_arc[u]] += forward[u] ? delta : -delta; - _excess[path_nodes[i-1]] -= delta; - _excess[u] += delta; - if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u); + u = v; + v = path_nodes[i]; + pa = pred_arc[v]; + delta = std::min(_res_cap[pa], _excess[u]); + _res_cap[pa] -= delta; + _res_cap[_reverse[pa]] += delta; + _excess[u] -= delta; + _excess[v] += delta; + if (_excess[v] > 0 && _excess[v] <= delta) + _active_nodes.push_back(v); } } } - - // Compute node potentials for the original costs - ResidualCostMap res_cost(_orig_cost); - BellmanFord< ResDigraph, ResidualCostMap > - bf(*_res_graph, res_cost); - bf.init(0); bf.start(); - for (NodeIt n(_graph); n != INVALID; ++n) - (*_potential)[n] = bf.dist(n); - - // Handle non-zero lower bounds - if (_lower) { - for (ArcIt e(_graph); e != INVALID; ++e) - (*_flow)[e] += (*_lower)[e]; - } - return true; } - /// Execute the algorithm performing push and relabel operations. - bool startPushRelabel() { + /// Execute the algorithm performing push and relabel operations + void startPushRelabel() { // Paramters for heuristics -// const int BF_HEURISTIC_EPSILON_BOUND = 1000; -// const int BF_HEURISTIC_BOUND_FACTOR = 3; + const int BF_HEURISTIC_EPSILON_BOUND = 1000; + const int BF_HEURISTIC_BOUND_FACTOR = 3; - typename Digraph::template NodeMap hyper(_graph, false); - typename Digraph::template NodeMap pred_arc(_graph); - typename Digraph::template NodeMap forward(_graph); - typename Digraph::template NodeMap next_out(_graph); - typename Digraph::template NodeMap next_in(_graph); - typename Digraph::template NodeMap next_dir(_graph); - std::deque active_nodes; - -// int node_num = countNodes(_graph); + // Perform cost scaling phases + BoolVector hyper(_res_node_num, false); for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 1 : _epsilon / _alpha ) { -/* // "Early Termination" heuristic: use Bellman-Ford algorithm // to check if the current flow is optimal if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { - typedef ShiftMap< ResidualCostMap > ShiftCostMap; - ShiftCostMap shift_cost(_res_cost, 1); - BellmanFord bf(*_res_graph, shift_cost); + _arc_vec.clear(); + _cost_vec.clear(); + for (int j = 0; j != _res_arc_num; ++j) { + if (_res_cap[j] > 0) { + _arc_vec.push_back(IntPair(_source[j], _target[j])); + _cost_vec.push_back(_cost[j] + 1); + } + } + _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); + + BellmanFord bf(_sgr, _cost_map); bf.init(0); bool done = false; - int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); for (int i = 0; i < K && !done; ++i) done = bf.processNextWeakRound(); if (done) break; } -*/ // Saturate arcs not satisfying the optimality condition - Capacity delta; - for (ArcIt e(_graph); e != INVALID; ++e) { - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { - delta = _capacity[e] - (*_flow)[e]; - _excess[_graph.source(e)] -= delta; - _excess[_graph.target(e)] += delta; - (*_flow)[e] = _capacity[e]; - } - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { - _excess[_graph.target(e)] -= (*_flow)[e]; - _excess[_graph.source(e)] += (*_flow)[e]; - (*_flow)[e] = 0; + for (int a = 0; a != _res_arc_num; ++a) { + if (_res_cap[a] > 0 && + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { + Value delta = _res_cap[a]; + _excess[_source[a]] -= delta; + _excess[_target[a]] += delta; + _res_cap[a] = 0; + _res_cap[_reverse[a]] += delta; } } // Find active nodes (i.e. nodes with positive excess) - for (NodeIt n(_graph); n != INVALID; ++n) { - if (_excess[n] > 0) active_nodes.push_back(n); + for (int u = 0; u != _res_node_num; ++u) { + if (_excess[u] > 0) _active_nodes.push_back(u); } - // Initialize the next arc maps - for (NodeIt n(_graph); n != INVALID; ++n) { - next_out[n] = OutArcIt(_graph, n); - next_in[n] = InArcIt(_graph, n); - next_dir[n] = true; + // Initialize the next arcs + for (int u = 0; u != _res_node_num; ++u) { + _next_out[u] = _first_out[u]; } // Perform push and relabel operations - while (active_nodes.size() > 0) { + while (_active_nodes.size() > 0) { + LargeCost min_red_cost, rc; + Value delta; + int n, t, a, last_out = _res_arc_num; + // Select an active node (FIFO selection) - Node n = active_nodes[0], t; - bool relabel_enabled = true; + next_node: + n = _active_nodes.front(); + last_out = _sum_supply < 0 ? + _first_out[n+1] : _first_out[n+1] - 1; // Perform push operations if there are admissible arcs - if (_excess[n] > 0 && next_dir[n]) { - OutArcIt e = next_out[n]; - for ( ; e != INVALID; ++e) { - if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { - delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]); - t = _graph.target(e); + if (_excess[n] > 0) { + for (a = _next_out[n]; a != last_out; ++a) { + if (_res_cap[a] > 0 && + _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { + delta = std::min(_res_cap[a], _excess[n]); + t = _target[a]; // Push-look-ahead heuristic - Capacity ahead = -_excess[t]; - for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) { - if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) - ahead += _capacity[oe] - (*_flow)[oe]; - } - for (InArcIt ie(_graph, t); ie != INVALID; ++ie) { - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) - ahead += (*_flow)[ie]; + Value ahead = -_excess[t]; + int last_out_t = _sum_supply < 0 ? + _first_out[t+1] : _first_out[t+1] - 1; + for (int ta = _next_out[t]; ta != last_out_t; ++ta) { + if (_res_cap[ta] > 0 && + _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0) + ahead += _res_cap[ta]; + if (ahead >= delta) break; } if (ahead < 0) ahead = 0; // Push flow along the arc if (ahead < delta) { - (*_flow)[e] += ahead; + _res_cap[a] -= ahead; + _res_cap[_reverse[a]] += ahead; _excess[n] -= ahead; _excess[t] += ahead; - active_nodes.push_front(t); + _active_nodes.push_front(t); hyper[t] = true; - relabel_enabled = false; - break; + _next_out[n] = a; + goto next_node; } else { - (*_flow)[e] += delta; + _res_cap[a] -= delta; + _res_cap[_reverse[a]] += delta; _excess[n] -= delta; _excess[t] += delta; if (_excess[t] > 0 && _excess[t] <= delta) - active_nodes.push_back(t); + _active_nodes.push_back(t); } - if (_excess[n] == 0) break; + if (_excess[n] == 0) { + _next_out[n] = a; + goto remove_nodes; + } } } - if (e != INVALID) { - next_out[n] = e; - } else { - next_dir[n] = false; - } - } - - if (_excess[n] > 0 && !next_dir[n]) { - InArcIt e = next_in[n]; - for ( ; e != INVALID; ++e) { - if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { - delta = std::min((*_flow)[e], _excess[n]); - t = _graph.source(e); - - // Push-look-ahead heuristic - Capacity ahead = -_excess[t]; - for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) { - if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0) - ahead += _capacity[oe] - (*_flow)[oe]; - } - for (InArcIt ie(_graph, t); ie != INVALID; ++ie) { - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0) - ahead += (*_flow)[ie]; - } - if (ahead < 0) ahead = 0; - - // Push flow along the arc - if (ahead < delta) { - (*_flow)[e] -= ahead; - _excess[n] -= ahead; - _excess[t] += ahead; - active_nodes.push_front(t); - hyper[t] = true; - relabel_enabled = false; - break; - } else { - (*_flow)[e] -= delta; - _excess[n] -= delta; - _excess[t] += delta; - if (_excess[t] > 0 && _excess[t] <= delta) - active_nodes.push_back(t); - } - - if (_excess[n] == 0) break; - } - } - next_in[n] = e; + _next_out[n] = a; } // Relabel the node if it is still active (or hyper) - if (relabel_enabled && (_excess[n] > 0 || hyper[n])) { - LCost min_red_cost = std::numeric_limits::max() / 2; - for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) { - if ( _capacity[oe] - (*_flow)[oe] > 0 && - (*_red_cost)[oe] < min_red_cost ) - min_red_cost = (*_red_cost)[oe]; + if (_excess[n] > 0 || hyper[n]) { + min_red_cost = std::numeric_limits::max() / 2; + for (int a = _first_out[n]; a != last_out; ++a) { + rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; + if (_res_cap[a] > 0 && rc < min_red_cost) { + min_red_cost = rc; + } } - for (InArcIt ie(_graph, n); ie != INVALID; ++ie) { - if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) - min_red_cost = -(*_red_cost)[ie]; - } - (*_potential)[n] -= min_red_cost + _epsilon; + _pi[n] -= min_red_cost + _epsilon; hyper[n] = false; - // Reset the next arc maps - next_out[n] = OutArcIt(_graph, n); - next_in[n] = InArcIt(_graph, n); - next_dir[n] = true; + // Reset the next arc + _next_out[n] = _first_out[n]; } - + // Remove nodes that are not active nor hyper - while ( active_nodes.size() > 0 && - _excess[active_nodes[0]] <= 0 && - !hyper[active_nodes[0]] ) { - active_nodes.pop_front(); + remove_nodes: + while ( _active_nodes.size() > 0 && + _excess[_active_nodes.front()] <= 0 && + !hyper[_active_nodes.front()] ) { + _active_nodes.pop_front(); } } } - - // Compute node potentials for the original costs - ResidualCostMap res_cost(_orig_cost); - BellmanFord< ResDigraph, ResidualCostMap > - bf(*_res_graph, res_cost); - bf.init(0); bf.start(); - for (NodeIt n(_graph); n != INVALID; ++n) - (*_potential)[n] = bf.dist(n); - - // Handle non-zero lower bounds - if (_lower) { - for (ArcIt e(_graph); e != INVALID; ++e) - (*_flow)[e] += (*_lower)[e]; - } - return true; } }; //class CostScaling