1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/athos/min_cost_flow.cc Tue May 25 12:31:18 2004 +0000
1.3 @@ -0,0 +1,110 @@
1.4 +#include <iostream>
1.5 +#include "test_tools.h"
1.6 +#include <hugo/list_graph.h>
1.7 +#include <mincostflow.h>
1.8 +//#include <path.h>
1.9 +//#include <maps.h>
1.10 +
1.11 +using namespace std;
1.12 +using namespace hugo;
1.13 +
1.14 +
1.15 +
1.16 +bool passed = true;
1.17 +/*
1.18 +void check(bool rc, char *msg="") {
1.19 + passed = passed && rc;
1.20 + if(!rc) {
1.21 + std::cerr << "Test failed! ("<< msg << ")" << std::endl; \
1.22 +
1.23 +
1.24 + }
1.25 +}
1.26 +*/
1.27 +
1.28 +
1.29 +int main()
1.30 +{
1.31 +
1.32 + typedef ListGraph::Node Node;
1.33 + typedef ListGraph::Edge Edge;
1.34 +
1.35 + ListGraph graph;
1.36 +
1.37 + //Ahuja könyv példája
1.38 +
1.39 + Node s=graph.addNode();
1.40 + Node v1=graph.addNode();
1.41 + Node v2=graph.addNode();
1.42 + Node v3=graph.addNode();
1.43 + Node v4=graph.addNode();
1.44 + Node v5=graph.addNode();
1.45 + Node t=graph.addNode();
1.46 +
1.47 + Edge s_v1=graph.addEdge(s, v1);
1.48 + Edge v1_v2=graph.addEdge(v1, v2);
1.49 + Edge s_v3=graph.addEdge(s, v3);
1.50 + Edge v2_v4=graph.addEdge(v2, v4);
1.51 + Edge v2_v5=graph.addEdge(v2, v5);
1.52 + Edge v3_v5=graph.addEdge(v3, v5);
1.53 + Edge v4_t=graph.addEdge(v4, t);
1.54 + Edge v5_t=graph.addEdge(v5, t);
1.55 +
1.56 +
1.57 + ListGraph::EdgeMap<int> length(graph);
1.58 +
1.59 + length.set(s_v1, 6);
1.60 + length.set(v1_v2, 4);
1.61 + length.set(s_v3, 10);
1.62 + length.set(v2_v4, 5);
1.63 + length.set(v2_v5, 1);
1.64 + length.set(v3_v5, 4);
1.65 + length.set(v4_t, 8);
1.66 + length.set(v5_t, 8);
1.67 +
1.68 + /*
1.69 + ListGraph::EdgeMap<int> capacity(graph);
1.70 +
1.71 + capacity.set(s_v1, 2);
1.72 + capacity.set(v1_v2, 2);
1.73 + capacity.set(s_v3, 1);
1.74 + capacity.set(v2_v4, 1);
1.75 + capacity.set(v2_v5, 1);
1.76 + capacity.set(v3_v5, 1);
1.77 + capacity.set(v4_t, 1);
1.78 + capacity.set(v5_t, 2);
1.79 + */
1.80 +
1.81 + // ConstMap<Edge, int> const1map(1);
1.82 + std::cout << "Enhanced capacity scaling algorithm test (for the mincostflow problem)..." << std::endl;
1.83 +
1.84 + MinCostFlow< ListGraph, ListGraph::EdgeMap<int>, ListGraph::NodeMap<int> >
1.85 + min_cost_flow_test(graph, length, supply_demand);
1.86 +
1.87 + int k=1;
1.88 +
1.89 + /*
1.90 + check( min_cost_flow_test.run(s,t,k) == 1 && min_cost_flow_test.totalLength() == 19,"One path, total length should be 19");
1.91 +
1.92 + check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
1.93 +
1.94 + k=2;
1.95 +
1.96 + check( min_cost_flow_test.run(s,t,k) == 2 && min_cost_flow_test.totalLength() == 41,"Two paths, total length should be 41");
1.97 +
1.98 + check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
1.99 +
1.100 +
1.101 + k=4;
1.102 +
1.103 + check( min_cost_flow_test.run(s,t,k) == 3 && min_cost_flow_test.totalLength() == 64,"Three paths, total length should be 64");
1.104 +
1.105 + check(min_cost_flow_test.checkComplementarySlackness(), "Is the primal-dual solution pair really optimal?");
1.106 +
1.107 +
1.108 + cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
1.109 + << endl;
1.110 +
1.111 + return passed ? 0 : 1;
1.112 + */
1.113 +}
2.1 --- a/src/work/athos/mincostflow.h Mon May 24 14:13:03 2004 +0000
2.2 +++ b/src/work/athos/mincostflow.h Tue May 25 12:31:18 2004 +0000
2.3 @@ -59,7 +59,7 @@
2.4 class ModLengthMap {
2.5 //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
2.6 typedef typename Graph::template NodeMap<Length> NodeMap;
2.7 - const ResGraphType& G;
2.8 + const ResGraphType& res_graph;
2.9 // const EdgeIntMap& rev;
2.10 const LengthMap &ol;
2.11 const NodeMap &pot;
2.12 @@ -68,22 +68,22 @@
2.13 typedef typename LengthMap::ValueType ValueType;
2.14
2.15 ValueType operator[](typename ResGraphType::Edge e) const {
2.16 - if (G.forward(e))
2.17 - return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
2.18 + if (res_graph.forward(e))
2.19 + return ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
2.20 else
2.21 - return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
2.22 + return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);
2.23 }
2.24
2.25 - ModLengthMap(const ResGraphType& _G,
2.26 + ModLengthMap(const ResGraphType& _res_graph,
2.27 const LengthMap &o, const NodeMap &p) :
2.28 - G(_G), /*rev(_rev),*/ ol(o), pot(p){};
2.29 + res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
2.30 };//ModLengthMap
2.31
2.32
2.33 protected:
2.34
2.35 //Input
2.36 - const Graph& G;
2.37 + const Graph& graph;
2.38 const LengthMap& length;
2.39 const SupplyDemandMap& supply_demand;//supply or demand of nodes
2.40
2.41 @@ -104,8 +104,8 @@
2.42 public :
2.43
2.44
2.45 - MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
2.46 - length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
2.47 + MinCostFlow(Graph& _graph, LengthMap& _length, SupplyDemandMap& _supply_demand) : graph(_graph),
2.48 + length(_length), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
2.49
2.50
2.51 ///Runs the algorithm.
2.52 @@ -114,7 +114,7 @@
2.53
2.54 ///\todo May be it does make sense to be able to start with a nonzero
2.55 /// feasible primal-dual solution pair as well.
2.56 - int run() {
2.57 + void run() {
2.58
2.59 //Resetting variables from previous runs
2.60 //total_length = 0;
2.61 @@ -124,17 +124,18 @@
2.62 std::greater<SupplyDemand> > HeapType;
2.63
2.64 //A heap for the excess nodes
2.65 - HeapMap excess_nodes_map(G,-1);
2.66 + HeapMap excess_nodes_map(graph,-1);
2.67 HeapType excess_nodes(excess_nodes_map);
2.68
2.69 //A heap for the deficit nodes
2.70 - HeapMap deficit_nodes_map(G,-1);
2.71 + HeapMap deficit_nodes_map(graph,-1);
2.72 HeapType deficit_nodes(deficit_nodes_map);
2.73
2.74 //A container to store nonabundant arcs
2.75 list<Edge> nonabundant_arcs;
2.76 -
2.77 - FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
2.78 +
2.79 +
2.80 + FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
2.81 flow.set(e,0);
2.82 nonabundant_arcs.push_back(e);
2.83 }
2.84 @@ -145,12 +146,12 @@
2.85 typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
2.86
2.87 //A union-find structure to store the abundant components
2.88 - UFE::MapType abund_comp_map(G);
2.89 + UFE::MapType abund_comp_map(graph);
2.90 UFE abundant_components(abund_comp_map);
2.91
2.92
2.93
2.94 - FOR_EACH_LOC(typename Graph::NodeIt, n, G){
2.95 + FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
2.96 excess_deficit.set(n,supply_demand[n]);
2.97 //A supply node
2.98 if (excess_deficit[n] > 0){
2.99 @@ -174,10 +175,39 @@
2.100 //can be the maximum deficit here
2.101 SupplyDemand max_excess = delta;
2.102
2.103 + //ConstMap<Edge,SupplyDemand> ConstEdgeMap;
2.104 +
2.105 //We need a residual graph which is uncapacitated
2.106 - ResGraphType res_graph(G, flow);
2.107 + ResGraphType res_graph(graph, flow);
2.108 +
2.109 + //An EdgeMap to tell which arcs are abundant
2.110 + template typename Graph::EdgeMap<bool> abundant_arcs(graph);
2.111
2.112 -
2.113 + //Let's construct the sugraph consisting only of the abundant edges
2.114 + typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
2.115 + ConstNodeMap const_true_map(true);
2.116 + typedef SubGraphWrapper< Graph, ConstNodeMap,
2.117 + template typename Graph::EdgeMap<bool> >
2.118 + AbundantGraph;
2.119 + AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
2.120 +
2.121 + //Let's construct the residual graph for the abundant graph
2.122 + typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap>
2.123 + ResAbGraph;
2.124 + //Again uncapacitated
2.125 + ResAbGraph res_ab_graph(abundant_graph, flow);
2.126 +
2.127 + //We need things for the bfs
2.128 + typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
2.129 + typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>
2.130 + bfs_pred(res_ab_graph);
2.131 + NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
2.132 + //We want to run bfs-es (more) on this graph 'res_ab_graph'
2.133 + Bfs < ResAbGraph ,
2.134 + typename ResAbGraph::NodeMap<bool>,
2.135 + typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
2.136 + NullMap<typename ResAbGraph::Node, int> >
2.137 + bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
2.138
2.139 ModLengthMap mod_length(res_graph, length, potential);
2.140
2.141 @@ -205,13 +235,55 @@
2.142 Node b = abundant_components.find(res_graph.tail(i));
2.143 //Merge
2.144 if (a != b){
2.145 - //Find path and augment
2.146 - //!!!
2.147 - //Find path and augment
2.148 abundant_components.join(a,b);
2.149 + //We want to push the smaller
2.150 + //Which has greater absolut value excess/deficit
2.151 + Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
2.152 + //Which is the other
2.153 + Node non_root = ( a == root ) ? b : a ;
2.154 + abundant_components.makeRep(root);
2.155 + SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
2.156 + //Push the positive value
2.157 + if (excess_deficit[non_root] < 0)
2.158 + swap(root, non_root);
2.159 + //If the non_root node has excess/deficit at all
2.160 + if (qty_to_augment>0){
2.161 + //Find path and augment
2.162 + bfs.run(non_root);
2.163 + //root should be reached
2.164 +
2.165 + //Augmenting on the found path
2.166 + Node n=root;
2.167 + ResGraphEdge e;
2.168 + while (n!=non_root){
2.169 + e = bfs_pred(n);
2.170 + n = res_graph.tail(e);
2.171 + res_graph.augment(e,qty_to_augment);
2.172 + }
2.173 +
2.174 + //We know that non_root had positive excess
2.175 + excess_nodes[non_root] -= qty_to_augment;
2.176 + //But what about root node
2.177 + //It might have been positive and so became larger
2.178 + if (excess_deficit[root]>0){
2.179 + excess_nodes[root] += qty_to_augment;
2.180 + }
2.181 + else{
2.182 + //Or negative but not turned into positive
2.183 + deficit_nodes[root] -= qty_to_augment;
2.184 + }
2.185 +
2.186 + //Update the excess_deficit map
2.187 + excess_deficit[non_root] -= qty_to_augment;
2.188 + excess_deficit[root] += qty_to_augment;
2.189 +
2.190 +
2.191 + }
2.192 }
2.193 //What happens to i?
2.194 - nonabundant_arcs.erase(i);
2.195 + //Marci and Zsolt says I shouldn't do such things
2.196 + nonabundant_arcs.erase(i++);
2.197 + abundant_arcs[i] = true;
2.198 }
2.199 else
2.200 ++i;
2.201 @@ -222,19 +294,19 @@
2.202 Node s = excess_nodes.top();
2.203 SupplyDemand max_excess = excess_nodes[s];
2.204 Node t = deficit_nodes.top();
2.205 - if (max_excess < dificit_nodes[t]){
2.206 - max_excess = dificit_nodes[t];
2.207 + if (max_excess < deficit_nodes[t]){
2.208 + max_excess = deficit_nodes[t];
2.209 }
2.210
2.211
2.212 - while(max_excess > 0){
2.213 -
2.214 + while(max_excess > (n-1)*delta/n){
2.215 +
2.216
2.217 //s es t valasztasa
2.218 -
2.219 +
2.220 //Dijkstra part
2.221 dijkstra.run(s);
2.222 -
2.223 +
2.224 /*We know from theory that t can be reached
2.225 if (!dijkstra.reached(t)){
2.226 //There are no k paths from s to t
2.227 @@ -263,6 +335,11 @@
2.228 total_length -= length[e];
2.229 */
2.230 }
2.231 +
2.232 + //Update the excess_deficit map
2.233 + excess_deficit[s] -= delta;
2.234 + excess_deficit[t] += delta;
2.235 +
2.236
2.237 //Update the excess_nodes heap
2.238 if (delta >= excess_nodes[s]){
2.239 @@ -285,6 +362,15 @@
2.240 deficit_nodes[t] -= delta;
2.241 }
2.242 //Dijkstra part ends here
2.243 +
2.244 + //Choose s and t again
2.245 + s = excess_nodes.top();
2.246 + max_excess = excess_nodes[s];
2.247 + t = deficit_nodes.top();
2.248 + if (max_excess < deficit_nodes[t]){
2.249 + max_excess = deficit_nodes[t];
2.250 + }
2.251 +
2.252 }
2.253
2.254 /*
2.255 @@ -297,7 +383,7 @@
2.256 /*This is not necessary here
2.257 //Update the max_excess
2.258 max_excess = 0;
2.259 - FOR_EACH_LOC(typename Graph::NodeIt, n, G){
2.260 + FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
2.261 if (max_excess < excess_deficit[n]){
2.262 max_excess = excess_deficit[n];
2.263 }
2.264 @@ -336,9 +422,9 @@
2.265 bool checkComplementarySlackness(){
2.266 Length mod_pot;
2.267 Length fl_e;
2.268 - FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
2.269 + FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
2.270 //C^{\Pi}_{i,j}
2.271 - mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
2.272 + mod_pot = length[e]-potential[graph.head(e)]+potential[graph.tail(e)];
2.273 fl_e = flow[e];
2.274 // std::cout << fl_e << std::endl;
2.275 if (0<fl_e && fl_e<capacity[e]){