# HG changeset patch
# User kpeter
# Date 1203305656 0
# Node ID 2c6204d4b0f63a723ac4d796421632c2f30fe6be
# Parent ae092c63d3baba257648106bc8ab725cb2eb8faf
Add a cost scaling min cost flow algorithm.
Add a cost scaling algorithm, which is performing generalized
push-relabel operations. It is almost as efficient as the capacity
scaling algorithm, but slower than network simplex.
diff -r ae092c63d3ba -r 2c6204d4b0f6 lemon/cost_scaling.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/cost_scaling.h Mon Feb 18 03:34:16 2008 +0000
@@ -0,0 +1,561 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_COST_SCALING_H
+#define LEMON_COST_SCALING_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief Cost scaling algorithm for finding a minimum cost flow.
+
+#include <deque>
+#include <lemon/graph_adaptor.h>
+#include <lemon/graph_utils.h>
+#include <lemon/maps.h>
+#include <lemon/math.h>
+
+#include <lemon/circulation.h>
+#include <lemon/bellman_ford.h>
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief Implementation of the cost scaling algorithm for finding a
+ /// minimum cost flow.
+ ///
+ /// \ref CostScaling implements the cost scaling algorithm performing
+ /// generalized push-relabel operations for finding a minimum cost
+ /// flow.
+ ///
+ /// \tparam Graph The directed graph 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.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be \e non-negative \e integers.
+ /// - Supply values should be \e signed \e integers.
+ /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
+ /// - \c CapacityMap::Value and \c SupplyMap::Value must be
+ /// convertible to each other.
+ /// - All value types must be convertible to \c CostMap::Value, which
+ /// must be signed type.
+ ///
+ /// \note Edge 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, <tt>long long int</tt> type is used instead of
+ /// <tt>long int</tt> in the inside computations.
+ ///
+ /// \author Peter Kovacs
+
+ template < typename Graph,
+ typename LowerMap = typename Graph::template EdgeMap<int>,
+ typename CapacityMap = typename Graph::template EdgeMap<int>,
+ typename CostMap = typename Graph::template EdgeMap<int>,
+ typename SupplyMap = typename Graph::template NodeMap<int> >
+ class CostScaling
+ {
+ GRAPH_TYPEDEFS(typename Graph);
+
+ typedef typename CapacityMap::Value Capacity;
+ typedef typename CostMap::Value Cost;
+ typedef typename SupplyMap::Value Supply;
+ typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
+ typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
+
+ typedef ResGraphAdaptor< const Graph, Capacity,
+ CapacityEdgeMap, CapacityEdgeMap > ResGraph;
+ typedef typename ResGraph::Edge ResEdge;
+
+#if defined __GNUC__ && !defined __STRICT_ANSI__
+ typedef long long int LCost;
+#else
+ typedef long int LCost;
+#endif
+ typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
+
+ public:
+
+ /// The type of the flow map.
+ typedef CapacityEdgeMap FlowMap;
+ /// The type of the potential map.
+ typedef typename Graph::template NodeMap<LCost> PotentialMap;
+
+ private:
+
+ /// \brief Map adaptor class for handling residual edge costs.
+ ///
+ /// \ref ResidualCostMap is a map adaptor class for handling
+ /// residual edge costs.
+ class ResidualCostMap : public MapBase<ResEdge, LCost>
+ {
+ private:
+
+ const LargeCostMap &_cost_map;
+
+ public:
+
+ ///\e
+ ResidualCostMap(const LargeCostMap &cost_map) :
+ _cost_map(cost_map) {}
+
+ ///\e
+ LCost operator[](const ResEdge &e) const {
+ return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
+ }
+
+ }; //class ResidualCostMap
+
+ /// \brief Map adaptor class for handling reduced edge costs.
+ ///
+ /// \ref ReducedCostMap is a map adaptor class for handling reduced
+ /// edge costs.
+ class ReducedCostMap : public MapBase<Edge, LCost>
+ {
+ private:
+
+ const Graph &_gr;
+ const LargeCostMap &_cost_map;
+ const PotentialMap &_pot_map;
+
+ public:
+
+ ///\e
+ ReducedCostMap( const Graph &gr,
+ const LargeCostMap &cost_map,
+ const PotentialMap &pot_map ) :
+ _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
+
+ ///\e
+ LCost operator[](const Edge &e) const {
+ return _cost_map[e] + _pot_map[_gr.source(e)]
+ - _pot_map[_gr.target(e)];
+ }
+
+ }; //class ReducedCostMap
+
+ private:
+
+ // Scaling factor
+ static const int ALPHA = 4;
+
+ // Paramters for heuristics
+ static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
+ static const int BF_HEURISTIC_BOUND_FACTOR = 3;
+
+ private:
+
+ // The directed graph the algorithm runs on
+ const Graph &_graph;
+ // The original lower bound map
+ const LowerMap *_lower;
+ // The modified capacity map
+ CapacityEdgeMap _capacity;
+ // The original cost map
+ const CostMap &_orig_cost;
+ // The scaled cost map
+ LargeCostMap _cost;
+ // The modified supply map
+ SupplyNodeMap _supply;
+ bool _valid_supply;
+
+ // Edge map of the current flow
+ FlowMap _flow;
+ // Node map of the current potentials
+ PotentialMap _potential;
+
+ // The residual graph
+ ResGraph _res_graph;
+ // The residual cost map
+ ResidualCostMap _res_cost;
+ // The reduced cost map
+ ReducedCostMap _red_cost;
+ // The excess map
+ SupplyNodeMap _excess;
+ // The epsilon parameter used for cost scaling
+ LCost _epsilon;
+
+ public:
+
+ /// \brief General constructor of the class (with lower bounds).
+ ///
+ /// General constructor of the class (with lower bounds).
+ ///
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param lower The lower bounds of the edges.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \param supply The supply values of the nodes (signed).
+ CostScaling( const Graph &graph,
+ const LowerMap &lower,
+ const CapacityMap &capacity,
+ const CostMap &cost,
+ const SupplyMap &supply ) :
+ _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
+ _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
+ _res_graph(graph, _capacity, _flow), _res_cost(_cost),
+ _red_cost(graph, _cost, _potential), _excess(graph, 0)
+ {
+ // Removing non-zero lower bounds
+ _capacity = subMap(capacity, lower);
+ Supply sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Supply s = supply[n];
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e)
+ s += lower[e];
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
+ s -= lower[e];
+ _supply[n] = s;
+ sum += s;
+ }
+ _valid_supply = sum == 0;
+ }
+
+ /// \brief General constructor of the class (without lower bounds).
+ ///
+ /// General constructor of the class (without lower bounds).
+ ///
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \param supply The supply values of the nodes (signed).
+ CostScaling( const Graph &graph,
+ const CapacityMap &capacity,
+ const CostMap &cost,
+ const SupplyMap &supply ) :
+ _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
+ _cost(graph), _supply(supply), _flow(graph, 0), _potential(graph, 0),
+ _res_graph(graph, _capacity, _flow), _res_cost(_cost),
+ _red_cost(graph, _cost, _potential), _excess(graph, 0)
+ {
+ // Checking the sum of supply values
+ Supply sum = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
+ _valid_supply = sum == 0;
+ }
+
+ /// \brief Simple constructor of the class (with lower bounds).
+ ///
+ /// Simple constructor of the class (with lower bounds).
+ ///
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param lower The lower bounds of the edges.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \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 Graph &graph,
+ const LowerMap &lower,
+ const CapacityMap &capacity,
+ const CostMap &cost,
+ Node s, Node t,
+ Supply flow_value ) :
+ _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
+ _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
+ _res_graph(graph, _capacity, _flow), _res_cost(_cost),
+ _red_cost(graph, _cost, _potential), _excess(graph, 0)
+ {
+ // Removing nonzero lower bounds
+ _capacity = subMap(capacity, lower);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ Supply sum = 0;
+ if (n == s) sum = flow_value;
+ if (n == t) sum = -flow_value;
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e)
+ sum += lower[e];
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
+ sum -= lower[e];
+ _supply[n] = sum;
+ }
+ _valid_supply = true;
+ }
+
+ /// \brief Simple constructor of the class (without lower bounds).
+ ///
+ /// Simple constructor of the class (without lower bounds).
+ ///
+ /// \param graph The directed graph the algorithm runs on.
+ /// \param capacity The capacities (upper bounds) of the edges.
+ /// \param cost The cost (length) values of the edges.
+ /// \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 Graph &graph,
+ const CapacityMap &capacity,
+ const CostMap &cost,
+ Node s, Node t,
+ Supply flow_value ) :
+ _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
+ _cost(graph), _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
+ _res_graph(graph, _capacity, _flow), _res_cost(_cost),
+ _red_cost(graph, _cost, _potential), _excess(graph, 0)
+ {
+ _supply[s] = flow_value;
+ _supply[t] = -flow_value;
+ _valid_supply = true;
+ }
+
+ /// \brief Runs the algorithm.
+ ///
+ /// Runs the algorithm.
+ ///
+ /// \return \c true if a feasible flow can be found.
+ bool run() {
+ init() && start();
+ }
+
+ /// \brief Returns a const reference to the edge map storing the
+ /// found flow.
+ ///
+ /// Returns a const reference to the edge map storing the found flow.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const FlowMap& flowMap() const {
+ return _flow;
+ }
+
+ /// \brief Returns a const reference to the node map storing the
+ /// found potentials (the dual solution).
+ ///
+ /// Returns 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;
+ }
+
+ /// \brief Returns the total cost of the found flow.
+ ///
+ /// Returns the total cost of the found flow. The complexity of the
+ /// function is \f$ O(e) \f$.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ Cost totalCost() const {
+ Cost c = 0;
+ for (EdgeIt e(_graph); e != INVALID; ++e)
+ c += _flow[e] * _orig_cost[e];
+ return c;
+ }
+
+ private:
+
+ /// Initializes the algorithm.
+ bool init() {
+ if (!_valid_supply) return false;
+
+ // Initializing the scaled cost map and the epsilon parameter
+ Cost max_cost = 0;
+ int node_num = countNodes(_graph);
+ for (EdgeIt 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];
+ }
+ _epsilon = max_cost * node_num;
+
+ // Finding a feasible flow using Circulation
+ Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
+ SupplyMap >
+ circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
+ _supply );
+ return circulation.flowMap(_flow).run();
+ }
+
+
+ /// Executes the algorithm.
+ bool start() {
+ std::deque<Node> active_nodes;
+ typename Graph::template NodeMap<bool> hyper(_graph, false);
+
+ int node_num = countNodes(_graph);
+ for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
+ 1 : _epsilon / ALPHA )
+ {
+ // Performing price refinement heuristic using Bellman-Ford
+ // algorithm
+ if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
+ typedef ShiftMap<ResidualCostMap> ShiftCostMap;
+ ShiftCostMap shift_cost(_res_cost, _epsilon);
+ BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
+ bf.init(0);
+ bool done = false;
+ int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
+ for (int i = 0; i < K && !done; ++i)
+ done = bf.processNextWeakRound();
+ if (done) {
+ for (NodeIt n(_graph); n != INVALID; ++n)
+ _potential[n] = bf.dist(n);
+ continue;
+ }
+ }
+
+ // Saturating edges not satisfying the optimality condition
+ Capacity delta;
+ for (EdgeIt 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;
+ }
+ }
+
+ // Finding active nodes (i.e. nodes with positive excess)
+ for (NodeIt n(_graph); n != INVALID; ++n)
+ if (_excess[n] > 0) active_nodes.push_back(n);
+
+ // Performing push and relabel operations
+ while (active_nodes.size() > 0) {
+ Node n = active_nodes[0], t;
+ bool relabel_enabled = true;
+
+ // Performing push operations if there are admissible edges
+ if (_excess[n] > 0) {
+ for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+ if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
+ delta = _capacity[e] - _flow[e] <= _excess[n] ?
+ _capacity[e] - _flow[e] : _excess[n];
+ t = _graph.target(e);
+
+ // Push-look-ahead heuristic
+ Capacity ahead = -_excess[t];
+ for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
+ if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
+ ahead += _capacity[oe] - _flow[oe];
+ }
+ for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
+ if (_flow[ie] > 0 && -_red_cost[ie] < 0)
+ ahead += _flow[ie];
+ }
+ if (ahead < 0) ahead = 0;
+
+ // Pushing flow along the edge
+ 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;
+ }
+ }
+ }
+
+ if (_excess[n] > 0) {
+ for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+ if (_flow[e] > 0 && -_red_cost[e] < 0) {
+ delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
+ t = _graph.source(e);
+
+ // Push-look-ahead heuristic
+ Capacity ahead = -_excess[t];
+ for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
+ if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
+ ahead += _capacity[oe] - _flow[oe];
+ }
+ for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
+ if (_flow[ie] > 0 && -_red_cost[ie] < 0)
+ ahead += _flow[ie];
+ }
+ if (ahead < 0) ahead = 0;
+
+ // Pushing flow along the edge
+ 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;
+ }
+ }
+ }
+
+ if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
+ // Performing relabel operation if the node is still active
+ LCost min_red_cost = std::numeric_limits<LCost>::max();
+ for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
+ if ( _capacity[oe] - _flow[oe] > 0 &&
+ _red_cost[oe] < min_red_cost )
+ min_red_cost = _red_cost[oe];
+ }
+ for (InEdgeIt 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;
+ hyper[n] = false;
+ }
+
+ // Removing active nodes with non-positive excess
+ while ( active_nodes.size() > 0 &&
+ _excess[active_nodes[0]] <= 0 &&
+ !hyper[active_nodes[0]] ) {
+ active_nodes.pop_front();
+ }
+ }
+ }
+
+ // Handling non-zero lower bounds
+ if (_lower) {
+ for (EdgeIt e(_graph); e != INVALID; ++e)
+ _flow[e] += (*_lower)[e];
+ }
+ return true;
+ }
+
+ }; //class CostScaling
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_COST_SCALING_H