Document the file itself.
2 #ifndef HUGO_MINCOSTFLOWS_H
3 #define HUGO_MINCOSTFLOWS_H
7 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
10 #include <hugo/dijkstra.h>
11 #include <hugo/graph_wrapper.h>
12 #include <hugo/maps.h>
17 /// \addtogroup flowalgs
20 ///\brief Implementation of an algorithm for finding a flow of value \c k
21 ///(for small values of \c k) having minimal total cost between 2 nodes
24 /// The class \ref hugo::MinCostFlows "MinCostFlows" implements
25 /// an algorithm for finding a flow of value \c k
26 /// having minimal total cost
27 /// from a given source node to a given target node in an
28 /// edge-weighted directed graph having nonnegative integer capacities.
29 /// The range of the length (weight or cost) function can be nonnegative reals but
30 /// the range of the capacity function has to be the set of nonnegative integers.
31 /// This algorithm is intended to use only for for small values of \c k, since /// it is not a polinomial time algorithm for finding the minimum cost
32 /// maximal flow (in order to find the minimum cost flow of value \c k it
33 /// finds the minimum cost flow of value \c i for every
34 /// \c i between 0 and \c k).
36 ///\param Graph The directed graph type the algorithm runs on.
37 ///\param LengthMap The type of the length map.
38 ///\param CapacityMap The capacity map type.
40 ///\author Attila Bernath
41 template <typename Graph, typename LengthMap, typename CapacityMap>
46 typedef typename LengthMap::ValueType Length;
48 //Warning: this should be integer type
49 typedef typename CapacityMap::ValueType Capacity;
51 typedef typename Graph::Node Node;
52 typedef typename Graph::NodeIt NodeIt;
53 typedef typename Graph::Edge Edge;
54 typedef typename Graph::OutEdgeIt OutEdgeIt;
55 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
58 typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
59 typedef typename ResGraphType::Edge ResGraphEdge;
62 typedef typename Graph::template NodeMap<Length> NodeMap;
63 const ResGraphType& G;
67 typedef typename LengthMap::KeyType KeyType;
68 typedef typename LengthMap::ValueType ValueType;
70 ValueType operator[](typename ResGraphType::Edge e) const {
72 return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
74 return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
77 ModLengthMap(const ResGraphType& _G,
78 const LengthMap &o, const NodeMap &p) :
79 G(_G), /*rev(_rev),*/ ol(o), pot(p){};
87 const LengthMap& length;
88 const CapacityMap& capacity;
95 //To store the potential (dual variables)
96 typedef typename Graph::template NodeMap<Length> PotentialMap;
97 PotentialMap potential;
105 /// The constructor of the class.
107 ///\param _G The directed graph the algorithm runs on.
108 ///\param _length The length (weight or cost) of the edges.
109 ///\param _cap The capacity of the edges.
110 MinCostFlows(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G),
111 length(_length), capacity(_cap), flow(_G), potential(_G){ }
114 ///Runs the algorithm.
116 ///Runs the algorithm.
117 ///Returns k if there is a flow of value at least k edge-disjoint
119 ///Otherwise it returns the maximum value of a flow from s to t.
121 ///\param s The source node.
122 ///\param t The target node.
123 ///\param k The value of the flow we are looking for.
125 ///\todo May be it does make sense to be able to start with a nonzero
126 /// feasible primal-dual solution pair as well.
127 int run(Node s, Node t, int k) {
129 //Resetting variables from previous runs
132 for (typename Graph::EdgeIt e(G); e!=INVALID; ++e) flow.set(e, 0);
134 //Initialize the potential to zero
135 for (typename Graph::NodeIt n(G); n!=INVALID; ++n) potential.set(n, 0);
138 //We need a residual graph
139 ResGraphType res_graph(G, capacity, flow);
142 ModLengthMap mod_length(res_graph, length, potential);
144 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
149 if (!dijkstra.reached(t)){
150 //There are no flow of value k from s to t
154 //We have to change the potential
155 for(typename ResGraphType::NodeIt n(res_graph); n!=INVALID; ++n)
156 potential[n] += dijkstra.distMap()[n];
159 //Augmenting on the sortest path
163 e = dijkstra.pred(n);
164 n = dijkstra.predNode(n);
165 res_graph.augment(e,1);
166 //Let's update the total length
167 if (res_graph.forward(e))
168 total_length += length[e];
170 total_length -= length[e];
182 /// Gives back the total weight of the found flow.
184 ///This function gives back the total weight of the found flow.
185 ///Assumes that \c run() has been run and nothing changed since then.
186 Length totalLength(){
190 ///Returns a const reference to the EdgeMap \c flow.
192 ///Returns a const reference to the EdgeMap \c flow.
193 ///\pre \ref run() must
194 ///be called before using this function.
195 const EdgeIntMap &getFlow() const { return flow;}
197 ///Returns a const reference to the NodeMap \c potential (the dual solution).
199 ///Returns a const reference to the NodeMap \c potential (the dual solution).
200 /// \pre \ref run() must be called before using this function.
201 const PotentialMap &getPotential() const { return potential;}
203 /// Checking the complementary slackness optimality criteria
205 ///This function checks, whether the given solution is optimal
206 ///If executed after the call of \c run() then it should return with true.
207 ///This function only checks optimality, doesn't bother with feasibility.
208 ///It is meant for testing purposes.
210 bool checkComplementarySlackness(){
213 for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) {
215 mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
217 // std::cout << fl_e << std::endl;
218 if (0<fl_e && fl_e<capacity[e]){
223 if (mod_pot > 0 && fl_e != 0)
225 if (mod_pot < 0 && fl_e != capacity[e])
233 }; //class MinCostFlows
239 #endif //HUGO_MINCOSTFLOWS_H