Unionfind changes induced some bugs here. Also some augmentations made.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_MIN_COST_FLOW_H
20 #define LEMON_MIN_COST_FLOW_H
24 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
27 #include <lemon/dijkstra.h>
28 #include <lemon/graph_adaptor.h>
29 #include <lemon/maps.h>
34 /// \addtogroup flowalgs
37 ///\brief Implementation of an algorithm for finding a flow of value \c k
38 ///(for small values of \c k) having minimal total cost between 2 nodes
41 /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an
42 /// algorithm for finding a flow of value \c k having minimal total
43 /// cost from a given source node to a given target node in a
44 /// directed graph with a cost function on the edges. To
45 /// this end, the edge-capacities and edge-costs have to be
46 /// nonnegative. The edge-capacities should be integers, but the
47 /// edge-costs can be integers, reals or of other comparable
48 /// numeric type. This algorithm is intended to be used only for
49 /// small values of \c k, since it is only polynomial in k, not in
50 /// the length of k (which is log k): in order to find the minimum
51 /// cost flow of value \c k it finds the minimum cost flow of value
52 /// \c i for every \c i between 0 and \c k.
54 ///\param Graph The directed graph type the algorithm runs on.
55 ///\param LengthMap The type of the length map.
56 ///\param CapacityMap The capacity map type.
58 ///\author Attila Bernath
59 template <typename Graph, typename LengthMap, typename CapacityMap>
62 typedef typename LengthMap::Value Length;
64 //Warning: this should be integer type
65 typedef typename CapacityMap::Value Capacity;
67 typedef typename Graph::Node Node;
68 typedef typename Graph::NodeIt NodeIt;
69 typedef typename Graph::Edge Edge;
70 typedef typename Graph::OutEdgeIt OutEdgeIt;
71 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
73 typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
74 typedef typename ResGW::Edge ResGraphEdge;
79 const LengthMap& length;
80 const CapacityMap& capacity;
83 typedef typename Graph::template NodeMap<Length> PotentialMap;
84 PotentialMap potential;
92 typedef typename Graph::template NodeMap<Length> NodeMap;
94 const LengthMap &length;
97 typedef typename LengthMap::Key Key;
98 typedef typename LengthMap::Value Value;
100 ModLengthMap(const ResGW& _g,
101 const LengthMap &_length, const NodeMap &_pot) :
102 g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
104 Value operator[](typename ResGW::Edge e) const {
106 return length[e]-(pot[g.target(e)]-pot[g.source(e)]);
108 return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);
114 ModLengthMap mod_length;
115 Dijkstra<ResGW, ModLengthMap> dijkstra;
119 /*! \brief The constructor of the class.
121 \param _g The directed graph the algorithm runs on.
122 \param _length The length (cost) of the edges.
123 \param _cap The capacity of the edges.
124 \param _s Source node.
125 \param _t Target node.
127 MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
129 g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
131 res_graph(g, capacity, flow),
132 mod_length(res_graph, length, potential),
133 dijkstra(res_graph, mod_length) {
137 /*! Tries to augment the flow between s and t by 1.
138 The return value shows if the augmentation is successful.
142 if (!dijkstra.reached(t)) {
144 //Unsuccessful augmentation.
148 //We have to change the potential
149 for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
150 potential.set(n, potential[n]+dijkstra.distMap()[n]);
152 //Augmenting on the shortest path
156 e = dijkstra.predEdge(n);
157 n = dijkstra.predNode(n);
158 res_graph.augment(e,1);
159 //Let's update the total length
160 if (res_graph.forward(e))
161 total_length += length[e];
163 total_length -= length[e];
170 /*! \brief Runs the algorithm.
173 Returns k if there is a flow of value at least k from s to t.
174 Otherwise it returns the maximum value of a flow from s to t.
176 \param k The value of the flow we are looking for.
178 \todo May be it does make sense to be able to start with a nonzero
179 feasible primal-dual solution pair as well.
181 \todo If the actual flow value is bigger than k, then everything is
182 cleared and the algorithm starts from zero flow. Is it a good approach?
185 if (flowValue()>k) reset();
186 while (flowValue()<k && augment()) { }
190 /*! \brief The class is reset to zero flow and potential.
191 The class is reset to zero flow and potential.
195 for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
196 for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);
199 /*! Returns the value of the actual flow.
201 int flowValue() const {
203 for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
204 for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
208 /// Total cost of the found flow.
210 /// This function gives back the total cost of the found flow.
211 Length totalLength(){
215 ///Returns a const reference to the EdgeMap \c flow.
217 ///Returns a const reference to the EdgeMap \c flow.
218 const EdgeIntMap &getFlow() const { return flow;}
220 /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).
222 Returns a const reference to the NodeMap \c potential (the dual solution).
224 const PotentialMap &getPotential() const { return potential;}
226 /*! \brief Checking the complementary slackness optimality criteria.
228 This function checks, whether the given flow and potential
229 satisfy the complementary slackness conditions (i.e. these are optimal).
230 This function only checks optimality, doesn't bother with feasibility.
233 bool checkComplementarySlackness(){
236 for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
238 mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
240 if (0<fl_e && fl_e<capacity[e]) {
241 /// \todo better comparison is needed for real types, moreover,
242 /// this comparison here is superfluous.
247 if (mod_pot > 0 && fl_e != 0)
249 if (mod_pot < 0 && fl_e != capacity[e])
256 }; //class MinCostFlow
262 #endif //LEMON_MIN_COST_FLOW_H