[Lemon-commits] deba: r3278 - lemon/trunk/lemon
Lemon SVN
svn at lemon.cs.elte.hu
Mon May 7 13:42:19 CEST 2007
Author: deba
Date: Mon May 7 13:42:18 2007
New Revision: 3278
Added:
lemon/trunk/lemon/capacity_scaling.h
lemon/trunk/lemon/cycle_canceling.h
lemon/trunk/lemon/min_cost_flow.h
lemon/trunk/lemon/min_cost_max_flow.h
lemon/trunk/lemon/network_simplex.h
Modified:
lemon/trunk/lemon/Makefile.am
Log:
Various min cost flow solvers
Patch from Peter Kovacs
Modified: lemon/trunk/lemon/Makefile.am
==============================================================================
--- lemon/trunk/lemon/Makefile.am (original)
+++ lemon/trunk/lemon/Makefile.am Mon May 7 13:42:18 2007
@@ -42,12 +42,14 @@
lemon/bp_matching.h \
lemon/bpugraph_adaptor.h \
lemon/bucket_heap.h \
+ lemon/capacity_scaling.h \
lemon/circulation.h \
lemon/color.h \
lemon/config.h \
lemon/concept_check.h \
lemon/counter.h \
lemon/csp.h \
+ lemon/cycle_canceling.h \
lemon/dag_shortest_path.h \
lemon/dfs.h \
lemon/dijkstra.h \
@@ -89,10 +91,13 @@
lemon/matrix_maps.h \
lemon/max_matching.h \
lemon/min_cost_arborescence.h \
+ lemon/min_cost_flow.h \
+ lemon/min_cost_max_flow.h \
lemon/min_mean_cycle.h \
lemon/mip_glpk.h \
lemon/mip_cplex.h \
lemon/nagamochi_ibaraki.h \
+ lemon/network_simplex.h \
lemon/path.h \
lemon/path_utils.h \
lemon/polynomial.h \
Added: lemon/trunk/lemon/capacity_scaling.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/capacity_scaling.h Mon May 7 13:42:18 2007
@@ -0,0 +1,680 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_CAPACITY_SCALING_H
+#define LEMON_CAPACITY_SCALING_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief The capacity scaling algorithm for finding a minimum cost
+/// flow.
+
+#include <vector>
+#include <lemon/dijkstra.h>
+#include <lemon/graph_adaptor.h>
+
+#define WITH_SCALING
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief Implementation of the capacity scaling version of the
+ /// succesive shortest path algorithm for finding a minimum cost flow.
+ ///
+ /// \ref lemon::CapacityScaling "CapacityScaling" implements a
+ /// capacity scaling algorithm for finding a minimum cost flow.
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param LowerMap The type of the lower bound map.
+ /// \param CapacityMap The type of the capacity (upper bound) map.
+ /// \param CostMap The type of the cost (length) map.
+ /// \param SupplyMap The type of the supply map.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be nonnegative integers.
+ /// However \c CostMap::Value should be signed type.
+ /// - Supply values should be integers.
+ /// - \c LowerMap::Value must be convertible to
+ /// \c CapacityMap::Value and \c CapacityMap::Value must be
+ /// convertible to \c SupplyMap::Value.
+ ///
+ /// \author Peter Kovacs
+
+template < typename Graph,
+ typename LowerMap = typename Graph::template EdgeMap<int>,
+ typename CapacityMap = LowerMap,
+ typename CostMap = typename Graph::template EdgeMap<int>,
+ typename SupplyMap = typename Graph::template NodeMap
+ <typename CapacityMap::Value> >
+ class CapacityScaling
+ {
+ typedef typename Graph::Node Node;
+ typedef typename Graph::NodeIt NodeIt;
+ typedef typename Graph::Edge Edge;
+ typedef typename Graph::EdgeIt EdgeIt;
+ typedef typename Graph::InEdgeIt InEdgeIt;
+ typedef typename Graph::OutEdgeIt OutEdgeIt;
+
+ typedef typename LowerMap::Value Lower;
+ typedef typename CapacityMap::Value Capacity;
+ typedef typename CostMap::Value Cost;
+ typedef typename SupplyMap::Value Supply;
+ typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
+ typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
+
+ typedef ResGraphAdaptor< const Graph, Capacity,
+ CapacityRefMap, CapacityRefMap > ResGraph;
+ typedef typename ResGraph::Node ResNode;
+ typedef typename ResGraph::NodeIt ResNodeIt;
+ typedef typename ResGraph::Edge ResEdge;
+ typedef typename ResGraph::EdgeIt ResEdgeIt;
+
+ public:
+
+ /// \brief The type of the flow map.
+ typedef CapacityRefMap FlowMap;
+ /// \brief The type of the potential map.
+ typedef typename Graph::template NodeMap<Cost> PotentialMap;
+
+ protected:
+
+ /// \brief Map adaptor class for handling reduced edge costs.
+ class ReducedCostMap : public MapBase<ResEdge, Cost>
+ {
+ private:
+
+ const ResGraph &gr;
+ const CostMap &cost_map;
+ const PotentialMap &pot_map;
+
+ public:
+
+ typedef typename MapBase<ResEdge, Cost>::Value Value;
+ typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+ ReducedCostMap( const ResGraph &_gr,
+ const CostMap &_cost,
+ const PotentialMap &_pot ) :
+ gr(_gr), cost_map(_cost), pot_map(_pot) {}
+
+ Value operator[](const Key &e) const {
+ return ResGraph::forward(e) ?
+ cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
+ -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+ }
+
+ }; //class ReducedCostMap
+
+ /// \brief Map adaptor for \ref lemon::Dijkstra "Dijkstra" class to
+ /// update node potentials.
+ class PotentialUpdateMap : public MapBase<ResNode, Cost>
+ {
+ private:
+
+ PotentialMap *pot;
+ typedef std::pair<ResNode, Cost> Pair;
+ std::vector<Pair> data;
+
+ public:
+
+ typedef typename MapBase<ResNode, Cost>::Value Value;
+ typedef typename MapBase<ResNode, Cost>::Key Key;
+
+ void potentialMap(PotentialMap &_pot) {
+ pot = &_pot;
+ }
+
+ void init() {
+ data.clear();
+ }
+
+ void set(const Key &n, const Value &v) {
+ data.push_back(Pair(n, v));
+ }
+
+ void update() {
+ Cost val = data[data.size()-1].second;
+ for (int i = 0; i < data.size()-1; ++i)
+ (*pot)[data[i].first] += val - data[i].second;
+ }
+
+ }; //class PotentialUpdateMap
+
+#ifdef WITH_SCALING
+ /// \brief Map adaptor class for identifing deficit nodes.
+ class DeficitBoolMap : public MapBase<ResNode, bool>
+ {
+ private:
+
+ const SupplyRefMap &imb;
+ const Capacity δ
+
+ public:
+
+ DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
+ imb(_imb), delta(_delta) {}
+
+ bool operator[](const ResNode &n) const {
+ return imb[n] <= -delta;
+ }
+
+ }; //class DeficitBoolMap
+
+ /// \brief Map adaptor class for filtering edges with at least
+ /// \c delta residual capacity
+ class DeltaFilterMap : public MapBase<ResEdge, bool>
+ {
+ private:
+
+ const ResGraph &gr;
+ const Capacity δ
+
+ public:
+
+ typedef typename MapBase<ResEdge, Cost>::Value Value;
+ typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+ DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
+ gr(_gr), delta(_delta) {}
+
+ Value operator[](const Key &e) const {
+ return gr.rescap(e) >= delta;
+ }
+
+ }; //class DeltaFilterMap
+
+ typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
+ DeltaResGraph;
+
+ /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
+ class ResDijkstraTraits :
+ public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
+ {
+ public:
+
+ typedef PotentialUpdateMap DistMap;
+
+ static DistMap *createDistMap(const DeltaResGraph&) {
+ return new DistMap();
+ }
+
+ }; //class ResDijkstraTraits
+
+#else //WITHOUT_CAPACITY_SCALING
+ /// \brief Map adaptor class for identifing deficit nodes.
+ class DeficitBoolMap : public MapBase<ResNode, bool>
+ {
+ private:
+
+ const SupplyRefMap &imb;
+
+ public:
+
+ DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
+
+ bool operator[](const ResNode &n) const {
+ return imb[n] < 0;
+ }
+
+ }; //class DeficitBoolMap
+
+ /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
+ class ResDijkstraTraits :
+ public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
+ {
+ public:
+
+ typedef PotentialUpdateMap DistMap;
+
+ static DistMap *createDistMap(const ResGraph&) {
+ return new DistMap();
+ }
+
+ }; //class ResDijkstraTraits
+#endif
+
+ protected:
+
+ /// \brief The directed graph the algorithm runs on.
+ const Graph &graph;
+ /// \brief The original lower bound map.
+ const LowerMap *lower;
+ /// \brief The modified capacity map.
+ CapacityRefMap capacity;
+ /// \brief The cost map.
+ const CostMap &cost;
+ /// \brief The modified supply map.
+ SupplyRefMap supply;
+ /// \brief The sum of supply values equals zero.
+ bool valid_supply;
+
+ /// \brief The edge map of the current flow.
+ FlowMap flow;
+ /// \brief The potential node map.
+ PotentialMap potential;
+ /// \brief The residual graph.
+ ResGraph res_graph;
+ /// \brief The reduced cost map.
+ ReducedCostMap red_cost;
+
+ /// \brief The imbalance map.
+ SupplyRefMap imbalance;
+ /// \brief The excess nodes.
+ std::vector<ResNode> excess_nodes;
+ /// \brief The index of the next excess node.
+ int next_node;
+
+#ifdef WITH_SCALING
+ typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
+ ResDijkstra;
+ /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
+ /// shortest paths in the residual graph with respect to the
+ /// reduced edge costs.
+ ResDijkstra dijkstra;
+
+ /// \brief The delta parameter used for capacity scaling.
+ Capacity delta;
+ /// \brief Edge filter map.
+ DeltaFilterMap delta_filter;
+ /// \brief The delta residual graph.
+ DeltaResGraph dres_graph;
+ /// \brief Map for identifing deficit nodes.
+ DeficitBoolMap delta_deficit;
+
+#else //WITHOUT_CAPACITY_SCALING
+ typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
+ ResDijkstra;
+ /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
+ /// shortest paths in the residual graph with respect to the
+ /// reduced edge costs.
+ ResDijkstra dijkstra;
+ /// \brief Map for identifing deficit nodes.
+ DeficitBoolMap has_deficit;
+#endif
+
+ /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
+ typename ResDijkstra::PredMap pred;
+ /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
+ /// update node potentials.
+ PotentialUpdateMap updater;
+
+ 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).
+ CapacityScaling( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+ supply(_graph), flow(_graph, 0), potential(_graph, 0),
+ res_graph(_graph, capacity, flow),
+ red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+ delta(0), delta_filter(res_graph, delta),
+ dres_graph(res_graph, delta_filter),
+ dijkstra(dres_graph, red_cost), pred(dres_graph),
+ delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+ dijkstra(res_graph, red_cost), pred(res_graph),
+ has_deficit(imbalance)
+#endif
+ {
+ // Removing nonzero 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] = imbalance[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).
+ CapacityScaling( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+ supply(_supply), flow(_graph, 0), potential(_graph, 0),
+ res_graph(_graph, capacity, flow),
+ red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+ delta(0), delta_filter(res_graph, delta),
+ dres_graph(res_graph, delta_filter),
+ dijkstra(dres_graph, red_cost), pred(dres_graph),
+ delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+ dijkstra(res_graph, red_cost), pred(res_graph),
+ has_deficit(imbalance)
+#endif
+ {
+ // 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).
+ CapacityScaling( 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), cost(_cost),
+ supply(_graph), flow(_graph, 0), potential(_graph, 0),
+ res_graph(_graph, capacity, flow),
+ red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+ delta(0), delta_filter(res_graph, delta),
+ dres_graph(res_graph, delta_filter),
+ dijkstra(dres_graph, red_cost), pred(dres_graph),
+ delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+ dijkstra(res_graph, red_cost), pred(res_graph),
+ has_deficit(imbalance)
+#endif
+ {
+ // Removing nonzero lower bounds
+ capacity = subMap(_capacity, _lower);
+ for (NodeIt n(graph); n != INVALID; ++n) {
+ Supply s = 0;
+ if (n == _s) s = _flow_value;
+ if (n == _t) s = -_flow_value;
+ for (InEdgeIt e(graph, n); e != INVALID; ++e)
+ s += _lower[e];
+ for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+ s -= _lower[e];
+ supply[n] = imbalance[n] = s;
+ }
+ 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).
+ CapacityScaling( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ Node _s, Node _t,
+ Supply _flow_value ) :
+ graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+ supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
+ res_graph(_graph, capacity, flow),
+ red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+ delta(0), delta_filter(res_graph, delta),
+ dres_graph(res_graph, delta_filter),
+ dijkstra(dres_graph, red_cost), pred(dres_graph),
+ delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+ dijkstra(res_graph, red_cost), pred(res_graph),
+ has_deficit(imbalance)
+#endif
+ {
+ supply[_s] = _flow_value;
+ supply[_t] = -_flow_value;
+ valid_supply = true;
+ }
+
+ /// \brief Returns a const reference to the flow map.
+ ///
+ /// Returns a const reference to the flow map.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const FlowMap& flowMap() const {
+ return flow;
+ }
+
+ /// \brief Returns a const reference to the potential map (the dual
+ /// solution).
+ ///
+ /// Returns a const reference to the potential map (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] * cost[e];
+ return c;
+ }
+
+ /// \brief Runs the successive shortest path algorithm.
+ ///
+ /// Runs the successive shortest path algorithm.
+ ///
+ /// \return \c true if a feasible flow can be found.
+ bool run() {
+ return init() && start();
+ }
+
+ protected:
+
+ /// \brief Initializes the algorithm.
+ bool init() {
+ if (!valid_supply) return false;
+
+ // Initalizing Dijkstra class
+ updater.potentialMap(potential);
+ dijkstra.distMap(updater).predMap(pred);
+
+#ifdef WITH_SCALING
+ // Initilaizing delta value
+ Capacity max_cap = 0;
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ if (capacity[e] > max_cap) max_cap = capacity[e];
+ }
+ for (delta = 1; 2 * delta < max_cap; delta *= 2) ;
+#endif
+ return true;
+ }
+
+#ifdef WITH_SCALING
+ /// \brief Executes the capacity scaling version of the successive
+ /// shortest path algorithm.
+ bool start() {
+ typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
+ typedef typename DeltaResGraph::Edge DeltaResEdge;
+
+ // Processing capacity scaling phases
+ ResNode s, t;
+ for ( ; delta >= 1; delta = delta < 4 && delta > 1 ?
+ 1 : delta / 4 )
+ {
+ // Saturating edges not satisfying the optimality condition
+ Capacity r;
+ for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
+ if (red_cost[e] < 0) {
+ res_graph.augment(e, r = res_graph.rescap(e));
+ imbalance[dres_graph.target(e)] += r;
+ imbalance[dres_graph.source(e)] -= r;
+ }
+ }
+
+ // Finding excess nodes
+ excess_nodes.clear();
+ for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+ if (imbalance[n] >= delta) excess_nodes.push_back(n);
+ }
+ next_node = 0;
+
+ // Finding successive shortest paths
+ while (next_node < excess_nodes.size()) {
+ // Running Dijkstra
+ s = excess_nodes[next_node];
+ updater.init();
+ dijkstra.init();
+ dijkstra.addSource(s);
+ if ((t = dijkstra.start(delta_deficit)) == INVALID) {
+ if (delta > 1) {
+ ++next_node;
+ continue;
+ }
+ return false;
+ }
+
+ // Updating node potentials
+ updater.update();
+
+ // Augment along a shortest path from s to t
+ Capacity d = imbalance[s] < -imbalance[t] ?
+ imbalance[s] : -imbalance[t];
+ ResNode u = t;
+ ResEdge e;
+ if (d > delta) {
+ while ((e = pred[u]) != INVALID) {
+ if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
+ u = dres_graph.source(e);
+ }
+ }
+ u = t;
+ while ((e = pred[u]) != INVALID) {
+ res_graph.augment(e, d);
+ u = dres_graph.source(e);
+ }
+ imbalance[s] -= d;
+ imbalance[t] += d;
+ if (imbalance[s] < delta) ++next_node;
+ }
+ }
+
+ // Handling nonzero lower bounds
+ if (lower) {
+ for (EdgeIt e(graph); e != INVALID; ++e)
+ flow[e] += (*lower)[e];
+ }
+ return true;
+ }
+
+#else //WITHOUT_CAPACITY_SCALING
+ /// \brief Executes the successive shortest path algorithm without
+ /// capacity scaling.
+ bool start() {
+ // Finding excess nodes
+ for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+ if (imbalance[n] > 0) excess_nodes.push_back(n);
+ }
+ if (excess_nodes.size() == 0) return true;
+ next_node = 0;
+
+ // Finding successive shortest paths
+ ResNode s, t;
+ while ( imbalance[excess_nodes[next_node]] > 0 ||
+ ++next_node < excess_nodes.size() )
+ {
+ // Running Dijkstra
+ s = excess_nodes[next_node];
+ updater.init();
+ dijkstra.init();
+ dijkstra.addSource(s);
+ if ((t = dijkstra.start(has_deficit)) == INVALID)
+ return false;
+
+ // Updating node potentials
+ updater.update();
+
+ // Augmenting along a shortest path from s to t
+ Capacity delta = imbalance[s] < -imbalance[t] ?
+ imbalance[s] : -imbalance[t];
+ ResNode u = t;
+ ResEdge e;
+ while ((e = pred[u]) != INVALID) {
+ if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
+ u = res_graph.source(e);
+ }
+ u = t;
+ while ((e = pred[u]) != INVALID) {
+ res_graph.augment(e, delta);
+ u = res_graph.source(e);
+ }
+ imbalance[s] -= delta;
+ imbalance[t] += delta;
+ }
+
+ // Handling nonzero lower bounds
+ if (lower) {
+ for (EdgeIt e(graph); e != INVALID; ++e)
+ flow[e] += (*lower)[e];
+ }
+ return true;
+ }
+#endif
+
+ }; //class CapacityScaling
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CAPACITY_SCALING_H
Added: lemon/trunk/lemon/cycle_canceling.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/cycle_canceling.h Mon May 7 13:42:18 2007
@@ -0,0 +1,509 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_CYCLE_CANCELING_H
+#define LEMON_CYCLE_CANCELING_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief A cycle canceling algorithm for finding a minimum cost flow.
+
+#include <vector>
+#include <lemon/circulation.h>
+#include <lemon/graph_adaptor.h>
+
+/// \brief The used cycle canceling method.
+#define LIMITED_CYCLE_CANCELING
+//#define MIN_MEAN_CYCLE_CANCELING
+
+#ifdef LIMITED_CYCLE_CANCELING
+ #include <lemon/bellman_ford.h>
+ /// \brief The maximum number of iterations for the first execution
+ /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm.
+ #define STARTING_LIMIT 2
+ /// \brief The iteration limit for the
+ /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
+ /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
+ #define ALPHA_MUL 3
+ /// \brief The iteration limit for the
+ /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
+ /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
+ #define ALPHA_DIV 2
+
+//#define _ONLY_ONE_CYCLE_
+//#define _DEBUG_ITER_
+#endif
+
+#ifdef MIN_MEAN_CYCLE_CANCELING
+ #include <lemon/min_mean_cycle.h>
+ #include <lemon/path.h>
+#endif
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief Implementation of a cycle canceling algorithm for finding
+ /// a minimum cost flow.
+ ///
+ /// \ref lemon::CycleCanceling "CycleCanceling" implements a cycle
+ /// canceling algorithm for finding a minimum cost flow.
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param LowerMap The type of the lower bound map.
+ /// \param CapacityMap The type of the capacity (upper bound) map.
+ /// \param CostMap The type of the cost (length) map.
+ /// \param SupplyMap The type of the supply map.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be nonnegative integers.
+ /// However \c CostMap::Value should be signed type.
+ /// - Supply values should be integers.
+ /// - \c LowerMap::Value must be convertible to
+ /// \c CapacityMap::Value and \c CapacityMap::Value must be
+ /// convertible to \c SupplyMap::Value.
+ ///
+ /// \author Peter Kovacs
+
+template < typename Graph,
+ typename LowerMap = typename Graph::template EdgeMap<int>,
+ typename CapacityMap = LowerMap,
+ typename CostMap = typename Graph::template EdgeMap<int>,
+ typename SupplyMap = typename Graph::template NodeMap
+ <typename CapacityMap::Value> >
+ class CycleCanceling
+ {
+ typedef typename Graph::Node Node;
+ typedef typename Graph::NodeIt NodeIt;
+ typedef typename Graph::Edge Edge;
+ typedef typename Graph::EdgeIt EdgeIt;
+ typedef typename Graph::InEdgeIt InEdgeIt;
+ typedef typename Graph::OutEdgeIt OutEdgeIt;
+
+ typedef typename LowerMap::Value Lower;
+ typedef typename CapacityMap::Value Capacity;
+ typedef typename CostMap::Value Cost;
+ typedef typename SupplyMap::Value Supply;
+ typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
+ typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
+
+ typedef ResGraphAdaptor< const Graph, Capacity,
+ CapacityRefMap, CapacityRefMap > ResGraph;
+ typedef typename ResGraph::Node ResNode;
+ typedef typename ResGraph::NodeIt ResNodeIt;
+ typedef typename ResGraph::Edge ResEdge;
+ typedef typename ResGraph::EdgeIt ResEdgeIt;
+
+ public:
+
+ /// \brief The type of the flow map.
+ typedef CapacityRefMap FlowMap;
+
+ protected:
+
+ /// \brief Map adaptor class for demand map.
+ class DemandMap : public MapBase<Node, Supply>
+ {
+ private:
+
+ const SupplyMap *map;
+
+ public:
+
+ typedef typename MapBase<Node, Supply>::Value Value;
+ typedef typename MapBase<Node, Supply>::Key Key;
+
+ DemandMap(const SupplyMap &_map) : map(&_map) {}
+
+ Value operator[](const Key &e) const {
+ return -(*map)[e];
+ }
+
+ }; //class DemandMap
+
+ /// \brief Map adaptor class for handling residual edge costs.
+ class ResCostMap : public MapBase<ResEdge, Cost>
+ {
+ private:
+
+ const CostMap &cost_map;
+
+ public:
+
+ typedef typename MapBase<ResEdge, Cost>::Value Value;
+ typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+ ResCostMap(const CostMap &_cost) : cost_map(_cost) {}
+
+ Value operator[](const Key &e) const {
+ return ResGraph::forward(e) ? cost_map[e] : -cost_map[e];
+ }
+
+ }; //class ResCostMap
+
+ protected:
+
+ /// \brief The directed graph the algorithm runs on.
+ const Graph &graph;
+ /// \brief The original lower bound map.
+ const LowerMap *lower;
+ /// \brief The modified capacity map.
+ CapacityRefMap capacity;
+ /// \brief The cost map.
+ const CostMap &cost;
+ /// \brief The modified supply map.
+ SupplyRefMap supply;
+ /// \brief The modified demand map.
+ DemandMap demand;
+ /// \brief The sum of supply values equals zero.
+ bool valid_supply;
+
+ /// \brief The current flow.
+ FlowMap flow;
+ /// \brief The residual graph.
+ ResGraph res_graph;
+ /// \brief The residual cost map.
+ ResCostMap res_cost;
+
+ 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).
+ CycleCanceling( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+ supply(_graph), demand(supply), flow(_graph, 0),
+ res_graph(_graph, capacity, flow), res_cost(cost)
+ {
+ // Removing nonzero 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];
+ sum += (supply[n] = 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).
+ CycleCanceling( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+ supply(_supply), demand(supply), flow(_graph, 0),
+ res_graph(_graph, capacity, flow), res_cost(cost)
+ {
+ // 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).
+ CycleCanceling( 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), cost(_cost),
+ supply(_graph), demand(supply), flow(_graph, 0),
+ res_graph(_graph, capacity, flow), res_cost(cost)
+ {
+ // Removing nonzero lower bounds
+ capacity = subMap(_capacity, _lower);
+ for (NodeIt n(graph); n != INVALID; ++n) {
+ Supply s = 0;
+ if (n == _s) s = _flow_value;
+ if (n == _t) s = -_flow_value;
+ 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;
+ }
+ 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).
+ CycleCanceling( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ Node _s, Node _t,
+ Supply _flow_value ) :
+ graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+ supply(_graph, 0), demand(supply), flow(_graph, 0),
+ res_graph(_graph, capacity, flow), res_cost(cost)
+ {
+ supply[_s] = _flow_value;
+ supply[_t] = -_flow_value;
+ valid_supply = true;
+ }
+
+ /// \brief Returns a const reference to the flow map.
+ ///
+ /// Returns a const reference to the flow map.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const FlowMap& flowMap() const {
+ return flow;
+ }
+
+ /// \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] * cost[e];
+ return c;
+ }
+
+ /// \brief Runs the algorithm.
+ ///
+ /// Runs the algorithm.
+ ///
+ /// \return \c true if a feasible flow can be found.
+ bool run() {
+ return init() && start();
+ }
+
+ protected:
+
+ /// \brief Initializes the algorithm.
+ bool init() {
+ // Checking the sum of supply values
+ Supply sum = 0;
+ for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
+ if (sum != 0) return false;
+
+ // Finding a feasible flow
+ Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
+ CapacityRefMap, DemandMap >
+ circulation( graph, constMap<Edge>((Capacity)0),
+ capacity, demand, flow );
+ return circulation.run() == -1;
+ }
+
+#ifdef LIMITED_CYCLE_CANCELING
+ /// \brief Executes a cycle canceling algorithm using
+ /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited
+ /// iteration count.
+ bool start() {
+ typename BellmanFord<ResGraph, ResCostMap>::PredMap pred(res_graph);
+ typename ResGraph::template NodeMap<int> visited(res_graph);
+ std::vector<ResEdge> cycle;
+ int node_num = countNodes(graph);
+
+#ifdef _DEBUG_ITER_
+ int cycle_num = 0;
+#endif
+ int length_bound = STARTING_LIMIT;
+ bool optimal = false;
+ while (!optimal) {
+ BellmanFord<ResGraph, ResCostMap> bf(res_graph, res_cost);
+ bf.predMap(pred);
+ bf.init(0);
+ int iter_num = 0;
+ bool cycle_found = false;
+ while (!cycle_found) {
+ int curr_iter_num = iter_num + length_bound <= node_num ?
+ length_bound : node_num - iter_num;
+ iter_num += curr_iter_num;
+ int real_iter_num = curr_iter_num;
+ for (int i = 0; i < curr_iter_num; ++i) {
+ if (bf.processNextWeakRound()) {
+ real_iter_num = i;
+ break;
+ }
+ }
+ if (real_iter_num < curr_iter_num) {
+ optimal = true;
+ break;
+ } else {
+ // Searching for node disjoint negative cycles
+ for (ResNodeIt n(res_graph); n != INVALID; ++n)
+ visited[n] = 0;
+ int id = 0;
+ for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+ if (visited[n] > 0) continue;
+ visited[n] = ++id;
+ ResNode u = pred[n] == INVALID ?
+ INVALID : res_graph.source(pred[n]);
+ while (u != INVALID && visited[u] == 0) {
+ visited[u] = id;
+ u = pred[u] == INVALID ?
+ INVALID : res_graph.source(pred[u]);
+ }
+ if (u != INVALID && visited[u] == id) {
+ // Finding the negative cycle
+ cycle_found = true;
+ cycle.clear();
+ ResEdge e = pred[u];
+ cycle.push_back(e);
+ Capacity d = res_graph.rescap(e);
+ while (res_graph.source(e) != u) {
+ cycle.push_back(e = pred[res_graph.source(e)]);
+ if (res_graph.rescap(e) < d)
+ d = res_graph.rescap(e);
+ }
+#ifdef _DEBUG_ITER_
+ ++cycle_num;
+#endif
+ // Augmenting along the cycle
+ for (int i = 0; i < cycle.size(); ++i)
+ res_graph.augment(cycle[i], d);
+#ifdef _ONLY_ONE_CYCLE_
+ break;
+#endif
+ }
+ }
+ }
+
+ if (!cycle_found)
+ length_bound = length_bound * ALPHA_MUL / ALPHA_DIV;
+ }
+ }
+
+#ifdef _DEBUG_ITER_
+ std::cout << "Limited cycle canceling algorithm finished. "
+ << "Found " << cycle_num << " negative cycles."
+ << std::endl;
+#endif
+
+ // Handling nonzero lower bounds
+ if (lower) {
+ for (EdgeIt e(graph); e != INVALID; ++e)
+ flow[e] += (*lower)[e];
+ }
+ return true;
+ }
+#endif
+
+#ifdef MIN_MEAN_CYCLE_CANCELING
+ /// \brief Executes the minimum mean cycle canceling algorithm
+ /// using \ref lemon::MinMeanCycle "MinMeanCycle" class.
+ bool start() {
+ typedef Path<ResGraph> ResPath;
+ MinMeanCycle<ResGraph, ResCostMap> mmc(res_graph, res_cost);
+ ResPath cycle;
+
+#ifdef _DEBUG_ITER_
+ int cycle_num = 0;
+#endif
+ mmc.cyclePath(cycle).init();
+ if (mmc.findMinMean()) {
+ while (mmc.cycleLength() < 0) {
+#ifdef _DEBUG_ITER_
+ ++iter;
+#endif
+ // Finding the cycle
+ mmc.findCycle();
+
+ // Finding the largest flow amount that can be augmented
+ // along the cycle
+ Capacity delta = 0;
+ for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
+ if (delta == 0 || res_graph.rescap(e) < delta)
+ delta = res_graph.rescap(e);
+ }
+
+ // Augmenting along the cycle
+ for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
+ res_graph.augment(e, delta);
+
+ // Finding the minimum cycle mean for the modified residual
+ // graph
+ mmc.reset();
+ if (!mmc.findMinMean()) break;
+ }
+ }
+
+#ifdef _DEBUG_ITER_
+ std::cout << "Minimum mean cycle canceling algorithm finished. "
+ << "Found " << cycle_num << " negative cycles."
+ << std::endl;
+#endif
+
+ // Handling nonzero lower bounds
+ if (lower) {
+ for (EdgeIt e(graph); e != INVALID; ++e)
+ flow[e] += (*lower)[e];
+ }
+ return true;
+ }
+#endif
+
+ }; //class CycleCanceling
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CYCLE_CANCELING_H
Added: lemon/trunk/lemon/min_cost_flow.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/min_cost_flow.h Mon May 7 13:42:18 2007
@@ -0,0 +1,127 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_MIN_COST_FLOW_H
+#define LEMON_MIN_COST_FLOW_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief An efficient algorithm for finding a minimum cost flow.
+
+#include <lemon/network_simplex.h>
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief An efficient algorithm for finding a minimum cost flow.
+ ///
+ /// \ref lemon::MinCostFlow "MinCostFlow" implements an efficient
+ /// algorithm for finding a minimum cost flow.
+ ///
+ /// \note \ref lemon::MinCostFlow "MinCostFlow" is just an alias for
+ /// \ref lemon::NetworkSimplex "NetworkSimplex", wich is the most
+ /// efficient implementation for the minimum cost flow problem in the
+ /// Lemon library according to our benchmark tests.
+ ///
+ /// \note There are three implementations for the minimum cost flow
+ /// problem, that can be used in the same way.
+ /// - \ref lemon::CapacityScaling "CapacityScaling"
+ /// - \ref lemon::CycleCanceling "CycleCanceling"
+ /// - \ref lemon::NetworkSimplex "NetworkSimplex"
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param LowerMap The type of the lower bound map.
+ /// \param CapacityMap The type of the capacity (upper bound) map.
+ /// \param CostMap The type of the cost (length) map.
+ /// \param SupplyMap The type of the supply map.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be nonnegative integers.
+ /// However \c CostMap::Value should be signed type.
+ /// - Supply values should be integers.
+ /// - \c LowerMap::Value must be convertible to
+ /// \c CapacityMap::Value and \c CapacityMap::Value must be
+ /// convertible to \c SupplyMap::Value.
+ ///
+ /// \author Peter Kovacs
+
+template < typename Graph,
+ typename LowerMap = typename Graph::template EdgeMap<int>,
+ typename CapacityMap = LowerMap,
+ typename CostMap = typename Graph::template EdgeMap<int>,
+ typename SupplyMap = typename Graph::template NodeMap
+ <typename CapacityMap::Value> >
+ class MinCostFlow :
+ public NetworkSimplex< Graph,
+ LowerMap,
+ CapacityMap,
+ CostMap,
+ SupplyMap >
+ {
+ typedef NetworkSimplex< Graph, LowerMap, CapacityMap,
+ CostMap, SupplyMap >
+ MinCostFlowImpl;
+ typedef typename Graph::Node Node;
+ typedef typename SupplyMap::Value Supply;
+
+ public:
+
+ /// \brief General constructor of the class (with lower bounds).
+ MinCostFlow( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ MinCostFlowImpl(_graph, _lower, _capacity, _cost, _supply) {}
+
+ /// \brief General constructor of the class (without lower bounds).
+ MinCostFlow( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ MinCostFlowImpl(_graph, _capacity, _cost, _supply) {}
+
+ /// \brief Simple constructor of the class (with lower bounds).
+ MinCostFlow( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ Node _s, Node _t,
+ Supply _flow_value ) :
+ MinCostFlowImpl( _graph, _lower, _capacity, _cost,
+ _s, _t, _flow_value ) {}
+
+ /// \brief Simple constructor of the class (without lower bounds).
+ MinCostFlow( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ Node _s, Node _t,
+ Supply _flow_value ) :
+ MinCostFlowImpl( _graph, _capacity, _cost,
+ _s, _t, _flow_value ) {}
+
+ }; //class MinCostFlow
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_MIN_COST_FLOW_H
Added: lemon/trunk/lemon/min_cost_max_flow.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/min_cost_max_flow.h Mon May 7 13:42:18 2007
@@ -0,0 +1,156 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_MIN_COST_MAX_FLOW_H
+#define LEMON_MIN_COST_MAX_FLOW_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief An efficient algorithm for finding a minimum cost maximum
+/// flow.
+
+#include <lemon/preflow.h>
+#include <lemon/network_simplex.h>
+#include <lemon/maps.h>
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief An efficient algorithm for finding a minimum cost maximum
+ /// flow.
+ ///
+ /// \ref lemon::MinCostFlow "MinCostMaxFlow" implements an efficient
+ /// algorithm for finding a maximum flow having minimal total cost
+ /// from a given source node to a given target node in a directed
+ /// graph.
+ ///
+ /// \note \ref lemon::MinCostMaxFlow "MinCostMaxFlow" uses
+ /// \ref lemon::Preflow "Preflow" algorithm for finding the maximum
+ /// flow value and \ref lemon::NetworkSimplex "NetworkSimplex"
+ /// algorithm for finding a minimum cost flow of that value.
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param CapacityMap The type of the capacity (upper bound) map.
+ /// \param CostMap The type of the cost (length) map.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be nonnegative integers.
+ /// However \c CostMap::Value should be signed type.
+ ///
+ /// \author Peter Kovacs
+
+template < typename Graph,
+ typename CapacityMap = typename Graph::template EdgeMap<int>,
+ typename CostMap = typename Graph::template EdgeMap<int> >
+ class MinCostMaxFlow
+ {
+ typedef typename Graph::Node Node;
+ typedef typename Graph::Edge Edge;
+
+ typedef typename CapacityMap::Value Capacity;
+ typedef typename Graph::template NodeMap<Capacity> SupplyMap;
+ typedef NetworkSimplex< Graph, CapacityMap, CapacityMap,
+ CostMap, SupplyMap >
+ MinCostFlowImpl;
+
+ public:
+
+ /// \brief The type of the flow map.
+ typedef typename Graph::template EdgeMap<Capacity> FlowMap;
+
+ private:
+
+ /// \brief The directed graph the algorithm runs on.
+ const Graph &graph;
+ /// \brief The modified capacity map.
+ const CapacityMap &capacity;
+ /// \brief The cost map.
+ const CostMap &cost;
+ /// \brief The source node.
+ Node source;
+ /// \brief The target node.
+ Node target;
+ /// \brief The edge map of the found flow.
+ FlowMap flow;
+
+ typedef Preflow<Graph, Capacity, CapacityMap, FlowMap> PreflowImpl;
+ /// \brief \ref lemon::Preflow "Preflow" class for finding the
+ /// maximum flow value.
+ PreflowImpl preflow;
+
+ public:
+
+ /// \brief The constructor of the class.
+ ///
+ /// The constructor of the class.
+ ///
+ /// \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.
+ MinCostMaxFlow( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ Node _s, Node _t ) :
+ graph(_graph), capacity(_capacity), cost(_cost),
+ source(_s), target(_t), flow(_graph),
+ preflow(_graph, _s, _t, _capacity, flow)
+ {}
+
+ /// \brief Returns a const reference to the flow map.
+ ///
+ /// Returns a const reference to the flow map.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const FlowMap& flowMap() const {
+ return flow;
+ }
+
+ /// \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] * cost[e];
+ return c;
+ }
+
+ /// \brief Runs the algorithm.
+ void run() {
+ preflow.phase1();
+ MinCostFlowImpl mcf_impl( graph, capacity, cost,
+ source, target, preflow.flowValue() );
+ mcf_impl.run();
+ flow = mcf_impl.flowMap();
+ }
+
+ }; //class MinCostMaxFlow
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_MIN_COST_MAX_FLOW_H
Added: lemon/trunk/lemon/network_simplex.h
==============================================================================
--- (empty file)
+++ lemon/trunk/lemon/network_simplex.h Mon May 7 13:42:18 2007
@@ -0,0 +1,945 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_NETWORK_SIMPLEX_H
+#define LEMON_NETWORK_SIMPLEX_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief The network simplex algorithm for finding a minimum cost
+/// flow.
+
+#include <limits>
+#include <lemon/smart_graph.h>
+#include <lemon/graph_utils.h>
+
+/// \brief The pivot rule used in the algorithm.
+#define EDGE_BLOCK_PIVOT
+//#define FIRST_ELIGIBLE_PIVOT
+//#define BEST_ELIGIBLE_PIVOT
+//#define CANDIDATE_LIST_PIVOT
+//#define SORTED_LIST_PIVOT
+
+/// \brief State constant for edges at their lower bounds.
+#define LOWER 1
+/// \brief State constant for edges in the spanning tree.
+#define TREE 0
+/// \brief State constant for edges at their upper bounds.
+#define UPPER -1
+
+#ifdef EDGE_BLOCK_PIVOT
+ /// \brief Number of blocks for the "Edge Block" pivot rule.
+ #define BLOCK_NUM 100
+ /// \brief Lower bound for the number of edges to use "Edge Block"
+ // pivot rule instead of "First Eligible" pivot rule.
+ #define MIN_BOUND 1000
+#endif
+
+#ifdef CANDIDATE_LIST_PIVOT
+ #include <list>
+ /// \brief The maximum length of the edge list for the
+ /// "Candidate List" pivot rule.
+ #define LIST_LENGTH 100
+ /// \brief The maximum number of minor iterations between two major
+ /// itarations.
+ #define MINOR_LIMIT 50
+#endif
+
+#ifdef SORTED_LIST_PIVOT
+ #include <deque>
+ #include <algorithm>
+ /// \brief The maximum length of the edge list for the
+ /// "Sorted List" pivot rule.
+ #define LIST_LENGTH 500
+ #define LOWER_DIV 4
+#endif
+
+//#define _DEBUG_ITER_
+
+namespace lemon {
+
+ /// \addtogroup min_cost_flow
+ /// @{
+
+ /// \brief Implementation of the network simplex algorithm for
+ /// finding a minimum cost flow.
+ ///
+ /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the
+ /// network simplex algorithm for finding a minimum cost flow.
+ ///
+ /// \param Graph The directed graph type the algorithm runs on.
+ /// \param LowerMap The type of the lower bound map.
+ /// \param CapacityMap The type of the capacity (upper bound) map.
+ /// \param CostMap The type of the cost (length) map.
+ /// \param SupplyMap The type of the supply map.
+ ///
+ /// \warning
+ /// - Edge capacities and costs should be nonnegative integers.
+ /// However \c CostMap::Value should be signed type.
+ /// - Supply values should be integers.
+ /// - \c LowerMap::Value must be convertible to
+ /// \c CapacityMap::Value and \c CapacityMap::Value must be
+ /// convertible to \c SupplyMap::Value.
+ ///
+ /// \author Peter Kovacs
+
+template < typename Graph,
+ typename LowerMap = typename Graph::template EdgeMap<int>,
+ typename CapacityMap = LowerMap,
+ typename CostMap = typename Graph::template EdgeMap<int>,
+ typename SupplyMap = typename Graph::template NodeMap
+ <typename CapacityMap::Value> >
+ class NetworkSimplex
+ {
+ typedef typename LowerMap::Value Lower;
+ typedef typename CapacityMap::Value Capacity;
+ typedef typename CostMap::Value Cost;
+ typedef typename SupplyMap::Value Supply;
+
+ typedef SmartGraph SGraph;
+ typedef typename SGraph::Node Node;
+ typedef typename SGraph::NodeIt NodeIt;
+ typedef typename SGraph::Edge Edge;
+ typedef typename SGraph::EdgeIt EdgeIt;
+ typedef typename SGraph::InEdgeIt InEdgeIt;
+ typedef typename SGraph::OutEdgeIt OutEdgeIt;
+
+ typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
+ typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
+ typedef typename SGraph::template EdgeMap<Cost> SCostMap;
+ typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
+ typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
+
+ typedef typename SGraph::template NodeMap<int> IntNodeMap;
+ typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
+ typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
+ typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
+ typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
+
+ typedef typename Graph::template NodeMap<Node> NodeRefMap;
+ typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
+
+ public:
+
+ /// \brief The type of the flow map.
+ typedef typename Graph::template EdgeMap<Capacity> FlowMap;
+ /// \brief The type of the potential map.
+ typedef typename Graph::template NodeMap<Cost> PotentialMap;
+
+ protected:
+
+ /// \brief Map adaptor class for handling reduced edge costs.
+ class ReducedCostMap : public MapBase<Edge, Cost>
+ {
+ private:
+
+ const SGraph &gr;
+ const SCostMap &cost_map;
+ const SPotentialMap &pot_map;
+
+ public:
+
+ typedef typename MapBase<Edge, Cost>::Value Value;
+ typedef typename MapBase<Edge, Cost>::Key Key;
+
+ ReducedCostMap( const SGraph &_gr,
+ const SCostMap &_cm,
+ const SPotentialMap &_pm ) :
+ gr(_gr), cost_map(_cm), pot_map(_pm) {}
+
+ Value operator[](const Key &e) const {
+ return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+ }
+
+ }; //class ReducedCostMap
+
+ protected:
+
+ /// \brief The directed graph the algorithm runs on.
+ SGraph graph;
+ /// \brief The original graph.
+ const Graph &graph_ref;
+ /// \brief The original lower bound map.
+ const LowerMap *lower;
+ /// \brief The capacity map.
+ SCapacityMap capacity;
+ /// \brief The cost map.
+ SCostMap cost;
+ /// \brief The supply map.
+ SSupplyMap supply;
+ /// \brief The reduced cost map.
+ ReducedCostMap red_cost;
+ /// \brief The sum of supply values equals zero.
+ bool valid_supply;
+
+ /// \brief The edge map of the current flow.
+ SCapacityMap flow;
+ /// \brief The edge map of the found flow on the original graph.
+ FlowMap flow_result;
+ /// \brief The potential node map.
+ SPotentialMap potential;
+ /// \brief The potential node map on the original graph.
+ PotentialMap potential_result;
+
+ /// \brief Node reference for the original graph.
+ NodeRefMap node_ref;
+ /// \brief Edge reference for the original graph.
+ EdgeRefMap edge_ref;
+
+ /// \brief The depth node map of the spanning tree structure.
+ IntNodeMap depth;
+ /// \brief The parent node map of the spanning tree structure.
+ NodeNodeMap parent;
+ /// \brief The pred_edge node map of the spanning tree structure.
+ EdgeNodeMap pred_edge;
+ /// \brief The thread node map of the spanning tree structure.
+ NodeNodeMap thread;
+ /// \brief The forward node map of the spanning tree structure.
+ BoolNodeMap forward;
+ /// \brief The state edge map.
+ IntEdgeMap state;
+
+
+#ifdef EDGE_BLOCK_PIVOT
+ /// \brief The size of blocks for the "Edge Block" pivot rule.
+ int block_size;
+#endif
+#ifdef CANDIDATE_LIST_PIVOT
+ /// \brief The list of candidate edges for the "Candidate List"
+ /// pivot rule.
+ std::list<Edge> candidates;
+ /// \brief The number of minor iterations.
+ int minor_count;
+#endif
+#ifdef SORTED_LIST_PIVOT
+ /// \brief The list of candidate edges for the "Sorted List"
+ /// pivot rule.
+ std::deque<Edge> candidates;
+#endif
+
+ // Root node of the starting spanning tree.
+ Node root;
+ // The entering edge of the current pivot iteration.
+ Edge in_edge;
+ // Temporary nodes used in the current pivot iteration.
+ Node join, u_in, v_in, u_out, v_out;
+ Node right, first, second, last;
+ Node stem, par_stem, new_stem;
+ // The maximum augment amount along the cycle in the current pivot
+ // iteration.
+ Capacity delta;
+
+ 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).
+ NetworkSimplex( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
+ supply(graph), flow(graph), flow_result(_graph), potential(graph),
+ potential_result(_graph), depth(graph), parent(graph),
+ pred_edge(graph), thread(graph), forward(graph), state(graph),
+ node_ref(graph_ref), edge_ref(graph_ref),
+ red_cost(graph, cost, potential)
+ {
+ // Checking the sum of supply values
+ Supply sum = 0;
+ for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+ sum += _supply[n];
+ if (!(valid_supply = sum == 0)) return;
+
+ // Copying graph_ref to graph
+ copyGraph(graph, graph_ref)
+ .edgeMap(cost, _cost)
+ .nodeRef(node_ref)
+ .edgeRef(edge_ref)
+ .run();
+
+ // Removing nonzero lower bounds
+ for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
+ capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+ }
+ for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
+ Supply s = _supply[n];
+ for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
+ s += _lower[e];
+ for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
+ s -= _lower[e];
+ supply[node_ref[n]] = s;
+ }
+ }
+
+ /// \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).
+ NetworkSimplex( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ const SupplyMap &_supply ) :
+ graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
+ supply(graph), flow(graph), flow_result(_graph), potential(graph),
+ potential_result(_graph), depth(graph), parent(graph),
+ pred_edge(graph), thread(graph), forward(graph), state(graph),
+ node_ref(graph_ref), edge_ref(graph_ref),
+ red_cost(graph, cost, potential)
+ {
+ // Checking the sum of supply values
+ Supply sum = 0;
+ for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+ sum += _supply[n];
+ if (!(valid_supply = sum == 0)) return;
+
+ // Copying graph_ref to graph
+ copyGraph(graph, graph_ref)
+ .edgeMap(capacity, _capacity)
+ .edgeMap(cost, _cost)
+ .nodeMap(supply, _supply)
+ .nodeRef(node_ref)
+ .edgeRef(edge_ref)
+ .run();
+ }
+
+ /// \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).
+ NetworkSimplex( const Graph &_graph,
+ const LowerMap &_lower,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ typename Graph::Node _s,
+ typename Graph::Node _t,
+ typename SupplyMap::Value _flow_value ) :
+ graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
+ supply(graph), flow(graph), flow_result(_graph), potential(graph),
+ potential_result(_graph), depth(graph), parent(graph),
+ pred_edge(graph), thread(graph), forward(graph), state(graph),
+ node_ref(graph_ref), edge_ref(graph_ref),
+ red_cost(graph, cost, potential)
+ {
+ // Copying graph_ref to graph
+ copyGraph(graph, graph_ref)
+ .edgeMap(cost, _cost)
+ .nodeRef(node_ref)
+ .edgeRef(edge_ref)
+ .run();
+
+ // Removing nonzero lower bounds
+ for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
+ capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+ }
+ for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
+ Supply s = 0;
+ if (n == _s) s = _flow_value;
+ if (n == _t) s = -_flow_value;
+ for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
+ s += _lower[e];
+ for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
+ s -= _lower[e];
+ supply[node_ref[n]] = s;
+ }
+ 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).
+ NetworkSimplex( const Graph &_graph,
+ const CapacityMap &_capacity,
+ const CostMap &_cost,
+ typename Graph::Node _s,
+ typename Graph::Node _t,
+ typename SupplyMap::Value _flow_value ) :
+ graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
+ supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
+ potential_result(_graph), depth(graph), parent(graph),
+ pred_edge(graph), thread(graph), forward(graph), state(graph),
+ node_ref(graph_ref), edge_ref(graph_ref),
+ red_cost(graph, cost, potential)
+ {
+ // Copying graph_ref to graph
+ copyGraph(graph, graph_ref)
+ .edgeMap(capacity, _capacity)
+ .edgeMap(cost, _cost)
+ .nodeRef(node_ref)
+ .edgeRef(edge_ref)
+ .run();
+ supply[node_ref[_s]] = _flow_value;
+ supply[node_ref[_t]] = -_flow_value;
+ valid_supply = true;
+ }
+
+ /// \brief Returns a const reference to the flow map.
+ ///
+ /// Returns a const reference to the flow map.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const FlowMap& flowMap() const {
+ return flow_result;
+ }
+
+ /// \brief Returns a const reference to the potential map (the dual
+ /// solution).
+ ///
+ /// Returns a const reference to the potential map (the dual
+ /// solution).
+ ///
+ /// \pre \ref run() must be called before using this function.
+ const PotentialMap& potentialMap() const {
+ return potential_result;
+ }
+
+ /// \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 (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+ c += flow_result[e] * cost[edge_ref[e]];
+ return c;
+ }
+
+ /// \brief Runs the algorithm.
+ ///
+ /// Runs the algorithm.
+ ///
+ /// \return \c true if a feasible flow can be found.
+ bool run() {
+ return init() && start();
+ }
+
+ protected:
+
+ /// \brief Extends the underlaying graph and initializes all the
+ /// node and edge maps.
+ bool init() {
+ if (!valid_supply) return false;
+
+ // Initializing state and flow maps
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ flow[e] = 0;
+ state[e] = LOWER;
+ }
+
+ // Adding an artificial root node to the graph
+ root = graph.addNode();
+ parent[root] = INVALID;
+ pred_edge[root] = INVALID;
+ depth[root] = supply[root] = potential[root] = 0;
+
+ // Adding artificial edges to the graph and initializing the node
+ // maps of the spanning tree data structure
+ Supply sum = 0;
+ Node last = root;
+ Edge e;
+ Cost max_cost = std::numeric_limits<Cost>::max() / 4;
+ for (NodeIt u(graph); u != INVALID; ++u) {
+ if (u == root) continue;
+ thread[last] = u;
+ last = u;
+ parent[u] = root;
+ depth[u] = 1;
+ sum += supply[u];
+ if (supply[u] >= 0) {
+ e = graph.addEdge(u, root);
+ flow[e] = supply[u];
+ forward[u] = true;
+ potential[u] = max_cost;
+ } else {
+ e = graph.addEdge(root, u);
+ flow[e] = -supply[u];
+ forward[u] = false;
+ potential[u] = -max_cost;
+ }
+ cost[e] = max_cost;
+ capacity[e] = std::numeric_limits<Capacity>::max();
+ state[e] = TREE;
+ pred_edge[u] = e;
+ }
+ thread[last] = root;
+
+#ifdef EDGE_BLOCK_PIVOT
+ // Initializing block_size for the edge block pivot rule
+ int edge_num = countEdges(graph);
+ block_size = edge_num > MIN_BOUND ? edge_num / BLOCK_NUM + 1 : 1;
+#endif
+#ifdef CANDIDATE_LIST_PIVOT
+ minor_count = 0;
+#endif
+
+ return sum == 0;
+ }
+
+#ifdef FIRST_ELIGIBLE_PIVOT
+ /// \brief Finds entering edge according to the "First Eligible"
+ /// pivot rule.
+ bool findEnteringEdge(EdgeIt &next_edge) {
+ for (EdgeIt e = next_edge; e != INVALID; ++e) {
+ if (state[e] * red_cost[e] < 0) {
+ in_edge = e;
+ next_edge = ++e;
+ return true;
+ }
+ }
+ for (EdgeIt e(graph); e != next_edge; ++e) {
+ if (state[e] * red_cost[e] < 0) {
+ in_edge = e;
+ next_edge = ++e;
+ return true;
+ }
+ }
+ return false;
+ }
+#endif
+
+#ifdef BEST_ELIGIBLE_PIVOT
+ /// \brief Finds entering edge according to the "Best Eligible"
+ /// pivot rule.
+ bool findEnteringEdge() {
+ Cost min = 0;
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ if (state[e] * red_cost[e] < min) {
+ min = state[e] * red_cost[e];
+ in_edge = e;
+ }
+ }
+ return min < 0;
+ }
+#endif
+
+#ifdef EDGE_BLOCK_PIVOT
+ /// \brief Finds entering edge according to the "Edge Block"
+ /// pivot rule.
+ bool findEnteringEdge(EdgeIt &next_edge) {
+ if (block_size == 1) {
+ // Performing first eligible selection
+ for (EdgeIt e = next_edge; e != INVALID; ++e) {
+ if (state[e] * red_cost[e] < 0) {
+ in_edge = e;
+ next_edge = ++e;
+ return true;
+ }
+ }
+ for (EdgeIt e(graph); e != next_edge; ++e) {
+ if (state[e] * red_cost[e] < 0) {
+ in_edge = e;
+ next_edge = ++e;
+ return true;
+ }
+ }
+ return false;
+ } else {
+ // Performing edge block selection
+ Cost curr, min = 0;
+ EdgeIt min_edge(graph);
+ int cnt = 0;
+ for (EdgeIt e = next_edge; e != INVALID; ++e) {
+ if ((curr = state[e] * red_cost[e]) < min) {
+ min = curr;
+ min_edge = e;
+ }
+ if (++cnt == block_size) {
+ if (min < 0) break;
+ cnt = 0;
+ }
+ }
+ if (!(min < 0)) {
+ for (EdgeIt e(graph); e != next_edge; ++e) {
+ if ((curr = state[e] * red_cost[e]) < min) {
+ min = curr;
+ min_edge = e;
+ }
+ if (++cnt == block_size) {
+ if (min < 0) break;
+ cnt = 0;
+ }
+ }
+ }
+ in_edge = min_edge;
+ if ((next_edge = ++min_edge) == INVALID)
+ next_edge = EdgeIt(graph);
+ return min < 0;
+ }
+ }
+#endif
+
+#ifdef CANDIDATE_LIST_PIVOT
+ /// \brief Functor class for removing non-eligible edges from the
+ /// candidate list.
+ class RemoveFunc
+ {
+ private:
+ const IntEdgeMap &st;
+ const ReducedCostMap &rc;
+ public:
+ RemoveFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
+ st(_st), rc(_rc) {}
+ bool operator()(const Edge &e) {
+ return st[e] * rc[e] >= 0;
+ }
+ };
+
+ /// \brief Finds entering edge according to the "Candidate List"
+ /// pivot rule.
+ bool findEnteringEdge() {
+ static RemoveFunc remove_func(state, red_cost);
+ typedef typename std::list<Edge>::iterator ListIt;
+
+ candidates.remove_if(remove_func);
+ if (minor_count >= MINOR_LIMIT || candidates.size() == 0) {
+ // Major iteration
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ if (state[e] * red_cost[e] < 0) {
+ candidates.push_back(e);
+ if (candidates.size() == LIST_LENGTH) break;
+ }
+ }
+ if (candidates.size() == 0) return false;
+ }
+
+ // Minor iteration
+ ++minor_count;
+ Cost min = 0;
+ for (ListIt it = candidates.begin(); it != candidates.end(); ++it) {
+ if (state[*it] * red_cost[*it] < min) {
+ min = state[*it] * red_cost[*it];
+ in_edge = *it;
+ }
+ }
+ return true;
+ }
+#endif
+
+#ifdef SORTED_LIST_PIVOT
+ /// \brief Functor class to compare edges during sort of the
+ /// candidate list.
+ class SortFunc
+ {
+ private:
+ const IntEdgeMap &st;
+ const ReducedCostMap &rc;
+ public:
+ SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
+ st(_st), rc(_rc) {}
+ bool operator()(const Edge &e1, const Edge &e2) {
+ return st[e1] * rc[e1] < st[e2] * rc[e2];
+ }
+ };
+
+ /// \brief Finds entering edge according to the "Sorted List"
+ /// pivot rule.
+ bool findEnteringEdge() {
+ static SortFunc sort_func(state, red_cost);
+
+ // Minor iteration
+ while (candidates.size() > 0) {
+ in_edge = candidates.front();
+ candidates.pop_front();
+ if (state[in_edge] * red_cost[in_edge] < 0) return true;
+ }
+
+ // Major iteration
+ Cost curr, min = 0;
+ for (EdgeIt e(graph); e != INVALID; ++e) {
+ if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
+ candidates.push_back(e);
+ if (curr < min) min = curr;
+ if (candidates.size() >= LIST_LENGTH) break;
+ }
+ }
+ if (candidates.size() == 0) return false;
+ sort(candidates.begin(), candidates.end(), sort_func);
+ return true;
+ }
+#endif
+
+ /// \brief Finds the join node.
+ Node findJoinNode() {
+ // Finding the join node
+ Node u = graph.source(in_edge);
+ Node v = graph.target(in_edge);
+ while (u != v) {
+ if (depth[u] == depth[v]) {
+ u = parent[u];
+ v = parent[v];
+ }
+ else if (depth[u] > depth[v]) u = parent[u];
+ else v = parent[v];
+ }
+ return u;
+ }
+
+ /// \brief Finds the leaving edge of the cycle. Returns \c true if
+ /// the leaving edge is not the same as the entering edge.
+ bool findLeavingEdge() {
+ // Initializing first and second nodes according to the direction
+ // of the cycle
+ if (state[in_edge] == LOWER) {
+ first = graph.source(in_edge);
+ second = graph.target(in_edge);
+ } else {
+ first = graph.target(in_edge);
+ second = graph.source(in_edge);
+ }
+ delta = capacity[in_edge];
+ bool result = false;
+ Capacity d;
+ Edge e;
+
+ // Searching the cycle along the path form the first node to the
+ // root node
+ for (Node u = first; u != join; u = parent[u]) {
+ e = pred_edge[u];
+ d = forward[u] ? flow[e] : capacity[e] - flow[e];
+ if (d < delta) {
+ delta = d;
+ u_out = u;
+ u_in = first;
+ v_in = second;
+ result = true;
+ }
+ }
+ // Searching the cycle along the path form the second node to the
+ // root node
+ for (Node u = second; u != join; u = parent[u]) {
+ e = pred_edge[u];
+ d = forward[u] ? capacity[e] - flow[e] : flow[e];
+ if (d <= delta) {
+ delta = d;
+ u_out = u;
+ u_in = second;
+ v_in = first;
+ result = true;
+ }
+ }
+ return result;
+ }
+
+ /// \brief Changes flow and state edge maps.
+ void changeFlows(bool change) {
+ // Augmenting along the cycle
+ if (delta > 0) {
+ Capacity val = state[in_edge] * delta;
+ flow[in_edge] += val;
+ for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
+ flow[pred_edge[u]] += forward[u] ? -val : val;
+ }
+ for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
+ flow[pred_edge[u]] += forward[u] ? val : -val;
+ }
+ }
+ // Updating the state of the entering and leaving edges
+ if (change) {
+ state[in_edge] = TREE;
+ state[pred_edge[u_out]] =
+ (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
+ } else {
+ state[in_edge] = -state[in_edge];
+ }
+ }
+
+ /// \brief Updates thread and parent node maps.
+ void updateThreadParent() {
+ Node u;
+ v_out = parent[u_out];
+
+ // Handling the case when join and v_out coincide
+ bool par_first = false;
+ if (join == v_out) {
+ for (u = join; u != u_in && u != v_in; u = thread[u]) ;
+ if (u == v_in) {
+ par_first = true;
+ while (thread[u] != u_out) u = thread[u];
+ first = u;
+ }
+ }
+
+ // Finding the last successor of u_in (u) and the node after it
+ // (right) according to the thread index
+ for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
+ right = thread[u];
+ if (thread[v_in] == u_out) {
+ for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
+ if (last == u_out) last = thread[last];
+ }
+ else last = thread[v_in];
+
+ // Updating stem nodes
+ thread[v_in] = stem = u_in;
+ par_stem = v_in;
+ while (stem != u_out) {
+ thread[u] = new_stem = parent[stem];
+
+ // Finding the node just before the stem node (u) according to
+ // the original thread index
+ for (u = new_stem; thread[u] != stem; u = thread[u]) ;
+ thread[u] = right;
+
+ // Changing the parent node of stem and shifting stem and
+ // par_stem nodes
+ parent[stem] = par_stem;
+ par_stem = stem;
+ stem = new_stem;
+
+ // Finding the last successor of stem (u) and the node after it
+ // (right) according to the thread index
+ for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
+ right = thread[u];
+ }
+ parent[u_out] = par_stem;
+ thread[u] = last;
+
+ if (join == v_out && par_first) {
+ if (first != v_in) thread[first] = right;
+ } else {
+ for (u = v_out; thread[u] != u_out; u = thread[u]) ;
+ thread[u] = right;
+ }
+ }
+
+ /// \brief Updates pred_edge and forward node maps.
+ void updatePredEdge() {
+ Node u = u_out, v;
+ while (u != u_in) {
+ v = parent[u];
+ pred_edge[u] = pred_edge[v];
+ forward[u] = !forward[v];
+ u = v;
+ }
+ pred_edge[u_in] = in_edge;
+ forward[u_in] = (u_in == graph.source(in_edge));
+ }
+
+ /// \brief Updates depth and potential node maps.
+ void updateDepthPotential() {
+ depth[u_in] = depth[v_in] + 1;
+ potential[u_in] = forward[u_in] ?
+ potential[v_in] + cost[pred_edge[u_in]] :
+ potential[v_in] - cost[pred_edge[u_in]];
+
+ Node u = thread[u_in], v;
+ while (true) {
+ v = parent[u];
+ if (v == INVALID) break;
+ depth[u] = depth[v] + 1;
+ potential[u] = forward[u] ?
+ potential[v] + cost[pred_edge[u]] :
+ potential[v] - cost[pred_edge[u]];
+ if (depth[u] <= depth[v_in]) break;
+ u = thread[u];
+ }
+ }
+
+ /// \brief Executes the algorithm.
+ bool start() {
+ // Processing pivots
+#ifdef _DEBUG_ITER_
+ int iter_num = 0;
+#endif
+#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
+ EdgeIt next_edge(graph);
+ while (findEnteringEdge(next_edge))
+#else
+ while (findEnteringEdge())
+#endif
+ {
+ join = findJoinNode();
+ bool change = findLeavingEdge();
+ changeFlows(change);
+ if (change) {
+ updateThreadParent();
+ updatePredEdge();
+ updateDepthPotential();
+ }
+#ifdef _DEBUG_ITER_
+ ++iter_num;
+#endif
+ }
+
+#ifdef _DEBUG_ITER_
+ t_iter.stop();
+ std::cout << "Network Simplex algorithm finished. " << iter_num
+ << " pivot iterations performed.";
+#endif
+
+ // Checking if the flow amount equals zero on all the
+ // artificial edges
+ for (InEdgeIt e(graph, root); e != INVALID; ++e)
+ if (flow[e] > 0) return false;
+ for (OutEdgeIt e(graph, root); e != INVALID; ++e)
+ if (flow[e] > 0) return false;
+
+ // Copying flow values to flow_result
+ if (lower) {
+ for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+ flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
+ } else {
+ for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+ flow_result[e] = flow[edge_ref[e]];
+ }
+ // Copying potential values to potential_result
+ for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+ potential_result[n] = potential[node_ref[n]];
+
+ return true;
+ }
+
+ }; //class NetworkSimplex
+
+ ///@}
+
+} //namespace lemon
+
+#endif //LEMON_NETWORK_SIMPLEX_H
More information about the Lemon-commits
mailing list