athos@276
|
1 |
// -*- c++ -*-
|
athos@523
|
2 |
#ifndef HUGO_MINCOSTFLOWS_H
|
athos@523
|
3 |
#define HUGO_MINCOSTFLOWS_H
|
athos@276
|
4 |
|
klao@491
|
5 |
///\ingroup galgs
|
alpar@294
|
6 |
///\file
|
athos@523
|
7 |
///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
|
alpar@294
|
8 |
|
athos@276
|
9 |
#include <iostream>
|
athos@551
|
10 |
#include <hugo/dijkstra.h>
|
athos@607
|
11 |
#include <hugo/graph_wrapper.h>
|
athos@551
|
12 |
#include <hugo/maps.h>
|
athos@551
|
13 |
#include <vector>
|
athos@530
|
14 |
#include <for_each_macros.h>
|
athos@306
|
15 |
|
athos@276
|
16 |
namespace hugo {
|
athos@276
|
17 |
|
alpar@430
|
18 |
/// \addtogroup galgs
|
alpar@430
|
19 |
/// @{
|
athos@322
|
20 |
|
athos@523
|
21 |
///\brief Implementation of an algorithm for finding a flow of value \c k
|
athos@523
|
22 |
///(for small values of \c k) having minimal total cost between 2 nodes
|
athos@523
|
23 |
///
|
klao@310
|
24 |
///
|
athos@523
|
25 |
/// The class \ref hugo::MinCostFlows "MinCostFlows" implements
|
athos@523
|
26 |
/// an algorithm for finding a flow of value \c k
|
athos@523
|
27 |
///(for small values of \c k) having minimal total cost
|
klao@310
|
28 |
/// from a given source node to a given target node in an
|
athos@523
|
29 |
/// edge-weighted directed graph having nonnegative integer capacities.
|
athos@523
|
30 |
/// The range of the length (weight) function is nonnegative reals but
|
athos@523
|
31 |
/// the range of capacity function is the set of nonnegative integers.
|
athos@523
|
32 |
/// It is not a polinomial time algorithm for counting the minimum cost
|
athos@523
|
33 |
/// maximal flow, since it counts the minimum cost flow for every value 0..M
|
athos@523
|
34 |
/// where \c M is the value of the maximal flow.
|
alpar@456
|
35 |
///
|
alpar@456
|
36 |
///\author Attila Bernath
|
athos@530
|
37 |
template <typename Graph, typename LengthMap, typename CapacityMap>
|
athos@523
|
38 |
class MinCostFlows {
|
athos@276
|
39 |
|
klao@310
|
40 |
typedef typename LengthMap::ValueType Length;
|
athos@527
|
41 |
|
athos@530
|
42 |
//Warning: this should be integer type
|
athos@530
|
43 |
typedef typename CapacityMap::ValueType Capacity;
|
athos@511
|
44 |
|
athos@276
|
45 |
typedef typename Graph::Node Node;
|
athos@276
|
46 |
typedef typename Graph::NodeIt NodeIt;
|
athos@276
|
47 |
typedef typename Graph::Edge Edge;
|
athos@276
|
48 |
typedef typename Graph::OutEdgeIt OutEdgeIt;
|
athos@511
|
49 |
typedef typename Graph::template EdgeMap<int> EdgeIntMap;
|
athos@306
|
50 |
|
athos@527
|
51 |
// typedef ConstMap<Edge,int> ConstMap;
|
athos@306
|
52 |
|
athos@530
|
53 |
typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
|
athos@530
|
54 |
typedef typename ResGraphType::Edge ResGraphEdge;
|
athos@547
|
55 |
|
athos@306
|
56 |
class ModLengthMap {
|
athos@547
|
57 |
//typedef typename ResGraphType::template NodeMap<Length> NodeMap;
|
athos@547
|
58 |
typedef typename Graph::template NodeMap<Length> NodeMap;
|
athos@306
|
59 |
const ResGraphType& G;
|
athos@527
|
60 |
// const EdgeIntMap& rev;
|
klao@310
|
61 |
const LengthMap &ol;
|
klao@310
|
62 |
const NodeMap &pot;
|
athos@306
|
63 |
public :
|
athos@306
|
64 |
typedef typename LengthMap::KeyType KeyType;
|
athos@306
|
65 |
typedef typename LengthMap::ValueType ValueType;
|
athos@511
|
66 |
|
athos@306
|
67 |
ValueType operator[](typename ResGraphType::Edge e) const {
|
athos@527
|
68 |
if (G.forward(e))
|
athos@527
|
69 |
return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
|
athos@527
|
70 |
else
|
athos@527
|
71 |
return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
|
athos@306
|
72 |
}
|
athos@511
|
73 |
|
athos@530
|
74 |
ModLengthMap(const ResGraphType& _G,
|
klao@310
|
75 |
const LengthMap &o, const NodeMap &p) :
|
athos@527
|
76 |
G(_G), /*rev(_rev),*/ ol(o), pot(p){};
|
athos@511
|
77 |
};//ModLengthMap
|
athos@511
|
78 |
|
athos@511
|
79 |
|
athos@607
|
80 |
protected:
|
athos@306
|
81 |
|
athos@527
|
82 |
//Input
|
athos@276
|
83 |
const Graph& G;
|
athos@276
|
84 |
const LengthMap& length;
|
athos@530
|
85 |
const CapacityMap& capacity;
|
athos@276
|
86 |
|
athos@607
|
87 |
|
alpar@328
|
88 |
//auxiliary variables
|
athos@322
|
89 |
|
athos@551
|
90 |
//To store the flow
|
athos@527
|
91 |
EdgeIntMap flow;
|
athos@551
|
92 |
//To store the potentila (dual variables)
|
athos@547
|
93 |
typename Graph::template NodeMap<Length> potential;
|
athos@276
|
94 |
|
athos@322
|
95 |
//Container to store found paths
|
athos@551
|
96 |
//std::vector< std::vector<Edge> > paths;
|
athos@511
|
97 |
//typedef DirPath<Graph> DPath;
|
athos@511
|
98 |
//DPath paths;
|
athos@511
|
99 |
|
athos@511
|
100 |
|
athos@511
|
101 |
Length total_length;
|
athos@322
|
102 |
|
athos@607
|
103 |
|
athos@276
|
104 |
public :
|
klao@310
|
105 |
|
athos@276
|
106 |
|
athos@530
|
107 |
MinCostFlows(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G),
|
athos@547
|
108 |
length(_length), capacity(_cap), flow(_G), potential(_G){ }
|
athos@276
|
109 |
|
alpar@294
|
110 |
|
alpar@329
|
111 |
///Runs the algorithm.
|
alpar@329
|
112 |
|
alpar@329
|
113 |
///Runs the algorithm.
|
athos@306
|
114 |
///Returns k if there are at least k edge-disjoint paths from s to t.
|
alpar@329
|
115 |
///Otherwise it returns the number of found edge-disjoint paths from s to t.
|
athos@607
|
116 |
///\todo May be it does make sense to be able to start with a nonzero
|
athos@607
|
117 |
/// feasible primal-dual solution pair as well.
|
athos@306
|
118 |
int run(Node s, Node t, int k) {
|
athos@276
|
119 |
|
athos@530
|
120 |
//Resetting variables from previous runs
|
athos@530
|
121 |
total_length = 0;
|
athos@547
|
122 |
|
athos@530
|
123 |
FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
|
athos@530
|
124 |
flow.set(e,0);
|
athos@530
|
125 |
}
|
athos@547
|
126 |
|
athos@547
|
127 |
FOR_EACH_LOC(typename Graph::NodeIt, n, G){
|
athos@547
|
128 |
//cout << potential[n]<<endl;
|
athos@547
|
129 |
potential.set(n,0);
|
athos@547
|
130 |
}
|
athos@547
|
131 |
|
athos@511
|
132 |
|
athos@530
|
133 |
|
athos@527
|
134 |
//We need a residual graph
|
athos@527
|
135 |
ResGraphType res_graph(G, capacity, flow);
|
athos@306
|
136 |
|
athos@306
|
137 |
//Initialize the copy of the Dijkstra potential to zero
|
athos@547
|
138 |
|
athos@547
|
139 |
//typename ResGraphType::template NodeMap<Length> potential(res_graph);
|
athos@547
|
140 |
|
athos@547
|
141 |
|
athos@547
|
142 |
ModLengthMap mod_length(res_graph, length, potential);
|
athos@306
|
143 |
|
athos@306
|
144 |
Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
|
athos@322
|
145 |
|
athos@322
|
146 |
int i;
|
athos@322
|
147 |
for (i=0; i<k; ++i){
|
athos@276
|
148 |
dijkstra.run(s);
|
athos@276
|
149 |
if (!dijkstra.reached(t)){
|
athos@314
|
150 |
//There are no k paths from s to t
|
athos@322
|
151 |
break;
|
athos@276
|
152 |
};
|
athos@306
|
153 |
|
athos@306
|
154 |
{
|
athos@306
|
155 |
//We have to copy the potential
|
athos@306
|
156 |
typename ResGraphType::NodeIt n;
|
athos@306
|
157 |
for ( res_graph.first(n) ; res_graph.valid(n) ; res_graph.next(n) ) {
|
athos@547
|
158 |
potential[n] += dijkstra.distMap()[n];
|
athos@306
|
159 |
}
|
athos@306
|
160 |
}
|
athos@306
|
161 |
|
athos@306
|
162 |
|
athos@527
|
163 |
//Augmenting on the sortest path
|
athos@276
|
164 |
Node n=t;
|
athos@530
|
165 |
ResGraphEdge e;
|
athos@276
|
166 |
while (n!=s){
|
athos@291
|
167 |
e = dijkstra.pred(n);
|
athos@291
|
168 |
n = dijkstra.predNode(n);
|
athos@530
|
169 |
res_graph.augment(e,1);
|
athos@530
|
170 |
//Let's update the total length
|
athos@530
|
171 |
if (res_graph.forward(e))
|
athos@530
|
172 |
total_length += length[e];
|
athos@530
|
173 |
else
|
athos@530
|
174 |
total_length -= length[e];
|
athos@276
|
175 |
}
|
athos@276
|
176 |
|
athos@276
|
177 |
|
athos@276
|
178 |
}
|
athos@322
|
179 |
|
athos@322
|
180 |
|
athos@322
|
181 |
return i;
|
athos@276
|
182 |
}
|
athos@276
|
183 |
|
athos@530
|
184 |
|
athos@530
|
185 |
|
athos@547
|
186 |
|
athos@511
|
187 |
///This function gives back the total length of the found paths.
|
athos@511
|
188 |
///Assumes that \c run() has been run and nothing changed since then.
|
athos@511
|
189 |
Length totalLength(){
|
athos@511
|
190 |
return total_length;
|
athos@511
|
191 |
}
|
athos@511
|
192 |
|
athos@607
|
193 |
///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
|
athos@607
|
194 |
///be called before using this function.
|
athos@607
|
195 |
const EdgeIntMap &getFlow() const { return flow;}
|
athos@607
|
196 |
|
athos@607
|
197 |
///Returns a const reference to the NodeMap \c potential (the dual solution).
|
athos@607
|
198 |
/// \pre \ref run() must be called before using this function.
|
athos@607
|
199 |
const EdgeIntMap &getPotential() const { return potential;}
|
athos@607
|
200 |
|
athos@607
|
201 |
///This function checks, whether the given solution is optimal
|
athos@607
|
202 |
///Running after a \c run() should return with true
|
athos@607
|
203 |
///In this "state of the art" this only check optimality, doesn't bother with feasibility
|
athos@607
|
204 |
///
|
athos@607
|
205 |
///\todo Is this OK here?
|
athos@607
|
206 |
bool checkComplementarySlackness(){
|
athos@551
|
207 |
Length mod_pot;
|
athos@551
|
208 |
Length fl_e;
|
athos@551
|
209 |
FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
|
athos@551
|
210 |
//C^{\Pi}_{i,j}
|
athos@554
|
211 |
mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
|
athos@551
|
212 |
fl_e = flow[e];
|
athos@554
|
213 |
// std::cout << fl_e << std::endl;
|
athos@551
|
214 |
if (0<fl_e && fl_e<capacity[e]){
|
athos@551
|
215 |
if (mod_pot != 0)
|
athos@551
|
216 |
return false;
|
athos@551
|
217 |
}
|
athos@551
|
218 |
else{
|
athos@551
|
219 |
if (mod_pot > 0 && fl_e != 0)
|
athos@551
|
220 |
return false;
|
athos@551
|
221 |
if (mod_pot < 0 && fl_e != capacity[e])
|
athos@551
|
222 |
return false;
|
athos@551
|
223 |
}
|
athos@551
|
224 |
}
|
athos@551
|
225 |
return true;
|
athos@551
|
226 |
}
|
athos@551
|
227 |
|
athos@530
|
228 |
/*
|
athos@530
|
229 |
///\todo To be implemented later
|
athos@530
|
230 |
|
athos@511
|
231 |
///This function gives back the \c j-th path in argument p.
|
athos@511
|
232 |
///Assumes that \c run() has been run and nothing changed since then.
|
athos@519
|
233 |
/// \warning It is assumed that \c p is constructed to be a path of graph \c G. If \c j is greater than the result of previous \c run, then the result here will be an empty path.
|
athos@511
|
234 |
template<typename DirPath>
|
athos@511
|
235 |
void getPath(DirPath& p, int j){
|
athos@511
|
236 |
p.clear();
|
athos@511
|
237 |
typename DirPath::Builder B(p);
|
athos@511
|
238 |
for(typename std::vector<Edge>::iterator i=paths[j].begin();
|
athos@511
|
239 |
i!=paths[j].end(); ++i ){
|
athos@520
|
240 |
B.pushBack(*i);
|
athos@511
|
241 |
}
|
athos@511
|
242 |
|
athos@511
|
243 |
B.commit();
|
athos@511
|
244 |
}
|
athos@276
|
245 |
|
athos@530
|
246 |
*/
|
athos@530
|
247 |
|
athos@530
|
248 |
}; //class MinCostFlows
|
athos@276
|
249 |
|
alpar@430
|
250 |
///@}
|
athos@276
|
251 |
|
athos@276
|
252 |
} //namespace hugo
|
athos@276
|
253 |
|
athos@527
|
254 |
#endif //HUGO_MINCOSTFLOW_H
|