1.1 --- a/src/lemon/min_cost_flow.h Thu Oct 07 17:21:27 2004 +0000
1.2 +++ b/src/lemon/min_cost_flow.h Fri Oct 08 13:07:51 2004 +0000
1.3 @@ -70,108 +70,85 @@
1.4 typedef typename Graph::OutEdgeIt OutEdgeIt;
1.5 typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.6
1.7 -
1.8 typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
1.9 typedef typename ResGW::Edge ResGraphEdge;
1.10
1.11 + protected:
1.12 +
1.13 + const Graph& g;
1.14 + const LengthMap& length;
1.15 + const CapacityMap& capacity;
1.16 +
1.17 + EdgeIntMap flow;
1.18 + typedef typename Graph::template NodeMap<Length> PotentialMap;
1.19 + PotentialMap potential;
1.20 +
1.21 + Node s;
1.22 + Node t;
1.23 +
1.24 + Length total_length;
1.25 +
1.26 class ModLengthMap {
1.27 typedef typename Graph::template NodeMap<Length> NodeMap;
1.28 - const ResGW& G;
1.29 - const LengthMap &ol;
1.30 + const ResGW& g;
1.31 + const LengthMap &length;
1.32 const NodeMap &pot;
1.33 public :
1.34 typedef typename LengthMap::KeyType KeyType;
1.35 typedef typename LengthMap::ValueType ValueType;
1.36 +
1.37 + ModLengthMap(const ResGW& _g,
1.38 + const LengthMap &_length, const NodeMap &_pot) :
1.39 + g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
1.40
1.41 ValueType operator[](typename ResGW::Edge e) const {
1.42 - if (G.forward(e))
1.43 - return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
1.44 + if (g.forward(e))
1.45 + return length[e]-(pot[g.head(e)]-pot[g.tail(e)]);
1.46 else
1.47 - return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
1.48 + return -length[e]-(pot[g.head(e)]-pot[g.tail(e)]);
1.49 }
1.50
1.51 - ModLengthMap(const ResGW& _G,
1.52 - const LengthMap &o, const NodeMap &p) :
1.53 - G(_G), /*rev(_rev),*/ ol(o), pot(p){};
1.54 - };//ModLengthMap
1.55 + }; //ModLengthMap
1.56
1.57 -
1.58 - protected:
1.59 -
1.60 - //Input
1.61 - const Graph& G;
1.62 - const LengthMap& length;
1.63 - const CapacityMap& capacity;
1.64 -
1.65 -
1.66 - //auxiliary variables
1.67 -
1.68 - //To store the flow
1.69 - EdgeIntMap flow;
1.70 - //To store the potential (dual variables)
1.71 - typedef typename Graph::template NodeMap<Length> PotentialMap;
1.72 - PotentialMap potential;
1.73 -
1.74 -
1.75 - Length total_length;
1.76 -
1.77 + ResGW res_graph;
1.78 + ModLengthMap mod_length;
1.79 + Dijkstra<ResGW, ModLengthMap> dijkstra;
1.80
1.81 public :
1.82
1.83 - /// The constructor of the class.
1.84 + /*! \brief The constructor of the class.
1.85
1.86 - ///\param _G The directed graph the algorithm runs on.
1.87 - ///\param _length The length (weight or cost) of the edges.
1.88 - ///\param _cap The capacity of the edges.
1.89 - MinCostFlow(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G),
1.90 - length(_length), capacity(_cap), flow(_G), potential(_G){ }
1.91 + \param _g The directed graph the algorithm runs on.
1.92 + \param _length The length (weight or cost) of the edges.
1.93 + \param _cap The capacity of the edges.
1.94 + \param _s Source node.
1.95 + \param _t Target node.
1.96 + */
1.97 + MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
1.98 + Node _s, Node _t) :
1.99 + g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
1.100 + s(_s), t(_t),
1.101 + res_graph(g, capacity, flow),
1.102 + mod_length(res_graph, length, potential),
1.103 + dijkstra(res_graph, mod_length) {
1.104 + reset();
1.105 + }
1.106
1.107 -
1.108 - ///Runs the algorithm.
1.109 -
1.110 - ///Runs the algorithm.
1.111 - ///Returns k if there is a flow of value at least k edge-disjoint
1.112 - ///from s to t.
1.113 - ///Otherwise it returns the maximum value of a flow from s to t.
1.114 - ///
1.115 - ///\param s The source node.
1.116 - ///\param t The target node.
1.117 - ///\param k The value of the flow we are looking for.
1.118 - ///
1.119 - ///\todo May be it does make sense to be able to start with a nonzero
1.120 - /// feasible primal-dual solution pair as well.
1.121 - int run(Node s, Node t, int k) {
1.122 + /*! Tries to augment the flow between s and t by 1.
1.123 + The return value shows if the augmentation is successful.
1.124 + */
1.125 + bool augment() {
1.126 + dijkstra.run(s);
1.127 + if (!dijkstra.reached(t)) {
1.128
1.129 - //Resetting variables from previous runs
1.130 - total_length = 0;
1.131 -
1.132 - for (typename Graph::EdgeIt e(G); e!=INVALID; ++e) flow.set(e, 0);
1.133 + //Unsuccessful augmentation.
1.134 + return false;
1.135 + } else {
1.136
1.137 - //Initialize the potential to zero
1.138 - for (typename Graph::NodeIt n(G); n!=INVALID; ++n) potential.set(n, 0);
1.139 -
1.140 -
1.141 - //We need a residual graph
1.142 - ResGW res_graph(G, capacity, flow);
1.143 -
1.144 -
1.145 - ModLengthMap mod_length(res_graph, length, potential);
1.146 -
1.147 - Dijkstra<ResGW, ModLengthMap> dijkstra(res_graph, mod_length);
1.148 -
1.149 - int i;
1.150 - for (i=0; i<k; ++i){
1.151 - dijkstra.run(s);
1.152 - if (!dijkstra.reached(t)){
1.153 - //There are no flow of value k from s to t
1.154 - break;
1.155 - };
1.156 + //We have to change the potential
1.157 + for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
1.158 + potential[n] += dijkstra.distMap()[n];
1.159
1.160 - //We have to change the potential
1.161 - for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
1.162 - potential[n] += dijkstra.distMap()[n];
1.163 -
1.164 -
1.165 //Augmenting on the sortest path
1.166 Node n=t;
1.167 ResGraphEdge e;
1.168 @@ -186,19 +163,51 @@
1.169 total_length -= length[e];
1.170 }
1.171
1.172 -
1.173 + return true;
1.174 }
1.175 -
1.176 + }
1.177 +
1.178 + /*! \brief Runs the algorithm.
1.179 +
1.180 + Runs the algorithm.
1.181 + Returns k if there is a flow of value at least k from s to t.
1.182 + Otherwise it returns the maximum value of a flow from s to t.
1.183 +
1.184 + \param k The value of the flow we are looking for.
1.185 +
1.186 + \todo May be it does make sense to be able to start with a nonzero
1.187 + feasible primal-dual solution pair as well.
1.188 +
1.189 + \todo If the actual flow value is bigger than k, then everything is
1.190 + cleared and the algorithm starts from zero flow. Is it a good approach?
1.191 + */
1.192 + int run(int k) {
1.193 + if (flowValue()>k) reset();
1.194 + while (flowValue()<k && augment()) { }
1.195 + return flowValue();
1.196 + }
1.197
1.198 + /*! \brief The class is reset to zero flow and potential.
1.199 + The class is reset to zero flow and potential.
1.200 + */
1.201 + void reset() {
1.202 + total_length=0;
1.203 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
1.204 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);
1.205 + }
1.206 +
1.207 + /*! Returns the value of the actual flow.
1.208 + */
1.209 + int flowValue() const {
1.210 + int i=0;
1.211 + for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
1.212 + for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
1.213 return i;
1.214 }
1.215
1.216 + /// Total weight of the found flow.
1.217
1.218 -
1.219 - /// Gives back the total weight of the found flow.
1.220 -
1.221 - ///This function gives back the total weight of the found flow.
1.222 - ///Assumes that \c run() has been run and nothing changed since then.
1.223 + /// This function gives back the total weight of the found flow.
1.224 Length totalLength(){
1.225 return total_length;
1.226 }
1.227 @@ -206,29 +215,27 @@
1.228 ///Returns a const reference to the EdgeMap \c flow.
1.229
1.230 ///Returns a const reference to the EdgeMap \c flow.
1.231 - ///\pre \ref run() must
1.232 - ///be called before using this function.
1.233 const EdgeIntMap &getFlow() const { return flow;}
1.234
1.235 - ///Returns a const reference to the NodeMap \c potential (the dual solution).
1.236 + /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).
1.237
1.238 - ///Returns a const reference to the NodeMap \c potential (the dual solution).
1.239 - /// \pre \ref run() must be called before using this function.
1.240 + Returns a const reference to the NodeMap \c potential (the dual solution).
1.241 + */
1.242 const PotentialMap &getPotential() const { return potential;}
1.243
1.244 - /// Checking the complementary slackness optimality criteria
1.245 + /*! \brief Checking the complementary slackness optimality criteria.
1.246
1.247 - ///This function checks, whether the given solution is optimal
1.248 - ///If executed after the call of \c run() then it should return with true.
1.249 - ///This function only checks optimality, doesn't bother with feasibility.
1.250 - ///It is meant for testing purposes.
1.251 - ///
1.252 + This function checks, whether the given flow and potential
1.253 + satisfiy the complementary slackness cnditions (i.e. these are optimal).
1.254 + This function only checks optimality, doesn't bother with feasibility.
1.255 + For testing purpose.
1.256 + */
1.257 bool checkComplementarySlackness(){
1.258 Length mod_pot;
1.259 Length fl_e;
1.260 - for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) {
1.261 + for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.262 //C^{\Pi}_{i,j}
1.263 - mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
1.264 + mod_pot = length[e]-potential[g.head(e)]+potential[g.tail(e)];
1.265 fl_e = flow[e];
1.266 if (0<fl_e && fl_e<capacity[e]) {
1.267 /// \todo better comparison is needed for real types, moreover,
1.268 @@ -246,7 +253,6 @@
1.269 return true;
1.270 }
1.271
1.272 -
1.273 }; //class MinCostFlow
1.274
1.275 ///@}
2.1 --- a/src/lemon/suurballe.h Thu Oct 07 17:21:27 2004 +0000
2.2 +++ b/src/lemon/suurballe.h Fri Oct 08 13:07:51 2004 +0000
2.3 @@ -68,14 +68,16 @@
2.4
2.5 typedef ConstMap<Edge,int> ConstMap;
2.6
2.7 - //Input
2.8 const Graph& G;
2.9
2.10 + Node s;
2.11 + Node t;
2.12 +
2.13 //Auxiliary variables
2.14 //This is the capacity map for the mincostflow problem
2.15 ConstMap const1map;
2.16 //This MinCostFlow instance will actually solve the problem
2.17 - MinCostFlow<Graph, LengthMap, ConstMap> mincost_flow;
2.18 + MinCostFlow<Graph, LengthMap, ConstMap> min_cost_flow;
2.19
2.20 //Container to store found paths
2.21 std::vector< std::vector<Edge> > paths;
2.22 @@ -83,39 +85,40 @@
2.23 public :
2.24
2.25
2.26 - /// The constructor of the class.
2.27 + /*! \brief The constructor of the class.
2.28
2.29 - ///\param _G The directed graph the algorithm runs on.
2.30 - ///\param _length The length (weight or cost) of the edges.
2.31 - Suurballe(Graph& _G, LengthMap& _length) : G(_G),
2.32 - const1map(1), mincost_flow(_G, _length, const1map){}
2.33 + \param _G The directed graph the algorithm runs on.
2.34 + \param _length The length (weight or cost) of the edges.
2.35 + \param _s Source node.
2.36 + \param _t Target node.
2.37 + */
2.38 + Suurballe(Graph& _G, LengthMap& _length, Node _s, Node _t) :
2.39 + G(_G), s(_s), t(_t), const1map(1),
2.40 + min_cost_flow(_G, _length, const1map, _s, _t) { }
2.41
2.42 ///Runs the algorithm.
2.43
2.44 ///Runs the algorithm.
2.45 ///Returns k if there are at least k edge-disjoint paths from s to t.
2.46 - ///Otherwise it returns the number of found edge-disjoint paths from s to t.
2.47 + ///Otherwise it returns the number of edge-disjoint paths found
2.48 + ///from s to t.
2.49 ///
2.50 - ///\param s The source node.
2.51 - ///\param t The target node.
2.52 ///\param k How many paths are we looking for?
2.53 ///
2.54 - int run(Node s, Node t, int k) {
2.55 -
2.56 - int i = mincost_flow.run(s,t,k);
2.57 -
2.58 + int run(int k) {
2.59 + int i = min_cost_flow.run(k);
2.60
2.61 //Let's find the paths
2.62 //We put the paths into stl vectors (as an inner representation).
2.63 //In the meantime we lose the information stored in 'reversed'.
2.64 //We suppose the lengths to be positive now.
2.65
2.66 - //We don't want to change the flow of mincost_flow, so we make a copy
2.67 + //We don't want to change the flow of min_cost_flow, so we make a copy
2.68 //The name here suggests that the flow has only 0/1 values.
2.69 EdgeIntMap reversed(G);
2.70
2.71 for(typename Graph::EdgeIt e(G); e!=INVALID; ++e)
2.72 - reversed[e] = mincost_flow.getFlow()[e];
2.73 + reversed[e] = min_cost_flow.getFlow()[e];
2.74
2.75 paths.clear();
2.76 //total_length=0;
2.77 @@ -143,39 +146,32 @@
2.78 }
2.79
2.80
2.81 - ///Returns the total length of the paths
2.82 + ///Returns the total length of the paths.
2.83
2.84 ///This function gives back the total length of the found paths.
2.85 - ///\pre \ref run() must
2.86 - ///be called before using this function.
2.87 Length totalLength(){
2.88 - return mincost_flow.totalLength();
2.89 + return min_cost_flow.totalLength();
2.90 }
2.91
2.92 ///Returns the found flow.
2.93
2.94 ///This function returns a const reference to the EdgeMap \c flow.
2.95 - ///\pre \ref run() must
2.96 - ///be called before using this function.
2.97 - const EdgeIntMap &getFlow() const { return mincost_flow.flow;}
2.98 + const EdgeIntMap &getFlow() const { return min_cost_flow.flow;}
2.99
2.100 /// Returns the optimal dual solution
2.101
2.102 ///This function returns a const reference to the NodeMap
2.103 ///\c potential (the dual solution).
2.104 - /// \pre \ref run() must be called before using this function.
2.105 - const EdgeIntMap &getPotential() const { return mincost_flow.potential;}
2.106 + const EdgeIntMap &getPotential() const { return min_cost_flow.potential;}
2.107
2.108 ///Checks whether the complementary slackness holds.
2.109
2.110 ///This function checks, whether the given solution is optimal.
2.111 - ///It should return true after calling \ref run()
2.112 ///Currently this function only checks optimality,
2.113 ///doesn't bother with feasibility
2.114 ///It is meant for testing purposes.
2.115 - ///
2.116 bool checkComplementarySlackness(){
2.117 - return mincost_flow.checkComplementarySlackness();
2.118 + return min_cost_flow.checkComplementarySlackness();
2.119 }
2.120
2.121 ///Read the found paths.
3.1 --- a/src/test/min_cost_flow_test.cc Thu Oct 07 17:21:27 2004 +0000
3.2 +++ b/src/test/min_cost_flow_test.cc Fri Oct 08 13:07:51 2004 +0000
3.3 @@ -21,11 +21,9 @@
3.4 //#include <path.h>
3.5 //#include <maps.h>
3.6
3.7 -using namespace std;
3.8 using namespace lemon;
3.9
3.10
3.11 -
3.12 bool passed = true;
3.13 /*
3.14 void check(bool rc, char *msg="") {
3.15 @@ -41,11 +39,11 @@
3.16
3.17 int main()
3.18 {
3.19 + typedef ListGraph Graph;
3.20 + typedef Graph::Node Node;
3.21 + typedef Graph::Edge Edge;
3.22
3.23 - typedef ListGraph::Node Node;
3.24 - typedef ListGraph::Edge Edge;
3.25 -
3.26 - ListGraph graph;
3.27 + Graph graph;
3.28
3.29 //Ahuja könyv példája
3.30
3.31 @@ -67,7 +65,7 @@
3.32 Edge v5_t=graph.addEdge(v5, t);
3.33
3.34
3.35 - ListGraph::EdgeMap<int> length(graph);
3.36 + Graph::EdgeMap<int> length(graph);
3.37
3.38 length.set(s_v1, 6);
3.39 length.set(v1_v2, 4);
3.40 @@ -78,7 +76,7 @@
3.41 length.set(v4_t, 8);
3.42 length.set(v5_t, 8);
3.43
3.44 - ListGraph::EdgeMap<int> capacity(graph);
3.45 + Graph::EdgeMap<int> capacity(graph);
3.46
3.47 capacity.set(s_v1, 2);
3.48 capacity.set(v1_v2, 2);
3.49 @@ -92,31 +90,36 @@
3.50 // ConstMap<Edge, int> const1map(1);
3.51 std::cout << "Mincostflows algorithm test..." << std::endl;
3.52
3.53 - MinCostFlow< ListGraph, ListGraph::EdgeMap<int>, ListGraph::EdgeMap<int> >
3.54 - surb_test(graph, length, capacity);
3.55 + MinCostFlow< Graph, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
3.56 + surb_test(graph, length, capacity, s, t);
3.57
3.58 int k=1;
3.59
3.60 - check( surb_test.run(s,t,k) == 1 && surb_test.totalLength() == 19,"One path, total length should be 19");
3.61 + surb_test.augment();
3.62 + check( surb_test.flowValue() == 1 && surb_test.totalLength() == 19,"One path, total length should be 19");
3.63 +
3.64 + check( surb_test.run(k) == 1 && surb_test.totalLength() == 19,"One path, total length should be 19");
3.65
3.66 check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
3.67
3.68 k=2;
3.69
3.70 - check( surb_test.run(s,t,k) == 2 && surb_test.totalLength() == 41,"Two paths, total length should be 41");
3.71 + check( surb_test.run(k) == 2 && surb_test.totalLength() == 41,"Two paths, total length should be 41");
3.72
3.73 check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
3.74
3.75 -
3.76 + surb_test.augment();
3.77 + surb_test.augment();
3.78 + surb_test.augment();
3.79 k=4;
3.80
3.81 - check( surb_test.run(s,t,k) == 3 && surb_test.totalLength() == 64,"Three paths, total length should be 64");
3.82 + check( surb_test.run(k) == 3 && surb_test.totalLength() == 64,"Three paths, total length should be 64");
3.83
3.84 check(surb_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
3.85
3.86
3.87 - cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
3.88 - << endl;
3.89 + std::cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
3.90 + << std::endl;
3.91
3.92 return passed ? 0 : 1;
3.93
4.1 --- a/src/test/suurballe_test.cc Thu Oct 07 17:21:27 2004 +0000
4.2 +++ b/src/test/suurballe_test.cc Fri Oct 08 13:07:51 2004 +0000
4.3 @@ -20,21 +20,19 @@
4.4 //#include <path.h>
4.5 #include "test_tools.h"
4.6
4.7 -using namespace std;
4.8 using namespace lemon;
4.9
4.10
4.11 -
4.12 bool passed = true;
4.13
4.14
4.15 int main()
4.16 {
4.17 + typedef ListGraph Graph;
4.18 + typedef Graph::Node Node;
4.19 + typedef Graph::Edge Edge;
4.20
4.21 - typedef ListGraph::Node Node;
4.22 - typedef ListGraph::Edge Edge;
4.23 -
4.24 - ListGraph graph;
4.25 + Graph graph;
4.26
4.27 //Ahuja könyv példája
4.28
4.29 @@ -56,7 +54,7 @@
4.30 Edge v5_t=graph.addEdge(v5, t);
4.31
4.32
4.33 - ListGraph::EdgeMap<int> length(graph);
4.34 + Graph::EdgeMap<int> length(graph);
4.35
4.36 length.set(s_v1, 6);
4.37 length.set(v1_v2, 4);
4.38 @@ -71,29 +69,29 @@
4.39
4.40
4.41 int k=3;
4.42 - Suurballe< ListGraph, ListGraph::EdgeMap<int> >
4.43 - surb_test(graph, length);
4.44 + Suurballe< Graph, Graph::EdgeMap<int> >
4.45 + surb_test(graph, length, s, t);
4.46
4.47 - check( surb_test.run(s,t,k) == 2 && surb_test.totalLength() == 46,
4.48 + check( surb_test.run(k) == 2 && surb_test.totalLength() == 46,
4.49 "Two paths, total length should be 46");
4.50
4.51 check( surb_test.checkComplementarySlackness(),
4.52 "Complementary slackness conditions are not met.");
4.53
4.54 - // typedef DirPath<ListGraph> DPath;
4.55 + // typedef DirPath<Graph> DPath;
4.56 // DPath P(graph);
4.57
4.58 /*
4.59 surb_test.getPath(P,0);
4.60 check(P.length() == 4, "First path should contain 4 edges.");
4.61 - cout<<P.length()<<endl;
4.62 + std::cout<<P.length()<<std::endl;
4.63 surb_test.getPath(P,1);
4.64 check(P.length() == 3, "Second path: 3 edges.");
4.65 - cout<<P.length()<<endl;
4.66 + std::cout<<P.length()<<std::endl;
4.67 */
4.68
4.69 k=1;
4.70 - check( surb_test.run(s,t,k) == 1 && surb_test.totalLength() == 19,
4.71 + check( surb_test.run(k) == 1 && surb_test.totalLength() == 19,
4.72 "One path, total length should be 19");
4.73
4.74 check( surb_test.checkComplementarySlackness(),
4.75 @@ -102,8 +100,8 @@
4.76 // surb_test.getPath(P,0);
4.77 // check(P.length() == 4, "First path should contain 4 edges.");
4.78
4.79 - cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
4.80 - << endl;
4.81 + std::cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
4.82 + << std::endl;
4.83
4.84 return passed ? 0 : 1;
4.85