1.1 --- a/src/work/athos/mincostflow.h Thu May 13 17:33:40 2004 +0000
1.2 +++ b/src/work/athos/mincostflow.h Thu May 13 17:42:23 2004 +0000
1.3 @@ -36,13 +36,13 @@
1.4 /// where \c M is the value of the maximal flow.
1.5 ///
1.6 ///\author Attila Bernath
1.7 - template <typename Graph, typename LengthMap, typename SupplyMap>
1.8 + template <typename Graph, typename LengthMap, typename SupplyDemandMap>
1.9 class MinCostFlow {
1.10
1.11 typedef typename LengthMap::ValueType Length;
1.12
1.13
1.14 - typedef typename SupplyMap::ValueType Supply;
1.15 + typedef typename SupplyDemandMap::ValueType SupplyDemand;
1.16
1.17 typedef typename Graph::Node Node;
1.18 typedef typename Graph::NodeIt NodeIt;
1.19 @@ -84,7 +84,7 @@
1.20 //Input
1.21 const Graph& G;
1.22 const LengthMap& length;
1.23 - const SupplyMap& supply;//supply or demand of nodes
1.24 + const SupplyDemandMap& supply_demand;//supply or demand of nodes
1.25
1.26
1.27 //auxiliary variables
1.28 @@ -94,7 +94,7 @@
1.29 //To store the potentila (dual variables)
1.30 typename Graph::template NodeMap<Length> potential;
1.31 //To store excess-deficit values
1.32 - SupplyMap excess;
1.33 + SupplyDemandMap excess_deficit;
1.34
1.35
1.36 Length total_length;
1.37 @@ -103,39 +103,62 @@
1.38 public :
1.39
1.40
1.41 - MinCostFlow(Graph& _G, LengthMap& _length, SupplyMap& _supply) : G(_G),
1.42 - length(_length), supply(_supply), flow(_G), potential(_G){ }
1.43 + MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
1.44 + length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
1.45
1.46
1.47 ///Runs the algorithm.
1.48
1.49 ///Runs the algorithm.
1.50 - ///Returns k if there are at least k edge-disjoint paths from s to t.
1.51 - ///Otherwise it returns the number of found edge-disjoint paths from s to t.
1.52 +
1.53 ///\todo May be it does make sense to be able to start with a nonzero
1.54 /// feasible primal-dual solution pair as well.
1.55 int run() {
1.56
1.57 //Resetting variables from previous runs
1.58 - total_length = 0;
1.59 + //total_length = 0;
1.60 +
1.61 + typedef typename Graph::template NodeMap<int> HeapMap;
1.62 + typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
1.63 + std::greater<SupplyDemand> > HeapType;
1.64 +
1.65 + //A heap for the excess nodes
1.66 + HeapMap excess_nodes_map(G,-1);
1.67 + HeapType excess_nodes(excess_nodes_map);
1.68 +
1.69 + //A heap for the deficit nodes
1.70 + HeapMap deficit_nodes_map(G,-1);
1.71 + HeapType deficit_nodes(deficit_nodes_map);
1.72 +
1.73
1.74 FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
1.75 flow.set(e,0);
1.76 }
1.77
1.78 //Initial value for delta
1.79 - Supply delta = 0;
1.80 -
1.81 + SupplyDemand delta = 0;
1.82 +
1.83 FOR_EACH_LOC(typename Graph::NodeIt, n, G){
1.84 - if (delta < abs(supply[e])){
1.85 - delta = abs(supply[e]);
1.86 + excess_deficit.set(n,supply_demand[n]);
1.87 + //A supply node
1.88 + if (excess_deficit[n] > 0){
1.89 + excess_nodes.push(n,excess_deficit[n]);
1.90 }
1.91 - excess.set(n,supply[e]);
1.92 + //A demand node
1.93 + if (excess_deficit[n] < 0){
1.94 + deficit_nodes.push(n, - excess_deficit[n]);
1.95 + }
1.96 + //Finding out starting value of delta
1.97 + if (delta < abs(excess_deficit[n])){
1.98 + delta = abs(excess_deficit[n]);
1.99 + }
1.100 //Initialize the copy of the Dijkstra potential to zero
1.101 potential.set(n,0);
1.102 }
1.103 -
1.104
1.105 + //It'll be allright as an initial value, though this value
1.106 + //can be the maximum deficit here
1.107 + SupplyDemand max_excess = delta;
1.108
1.109 //We need a residual graph which is uncapacitated
1.110 ResGraphType res_graph(G, flow);
1.111 @@ -147,45 +170,102 @@
1.112 Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
1.113
1.114
1.115 - int i;
1.116 - for (i=0; i<k; ++i){
1.117 - dijkstra.run(s);
1.118 - if (!dijkstra.reached(t)){
1.119 - //There are no k paths from s to t
1.120 - break;
1.121 - };
1.122 + while (max_excess > 0){
1.123 +
1.124
1.125 - //We have to copy the potential
1.126 - FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
1.127 - potential[n] += dijkstra.distMap()[n];
1.128 + //Merge and stuff
1.129 +
1.130 + Node s = excess_nodes.top();
1.131 + SupplyDemand max_excess = excess_nodes[s];
1.132 + Node t = deficit_nodes.top();
1.133 + if (max_excess < dificit_nodes[t]){
1.134 + max_excess = dificit_nodes[t];
1.135 + }
1.136 +
1.137 +
1.138 + while(max_excess > ){
1.139 +
1.140 +
1.141 + //s es t valasztasa
1.142 +
1.143 + //Dijkstra part
1.144 + dijkstra.run(s);
1.145 +
1.146 + /*We know from theory that t can be reached
1.147 + if (!dijkstra.reached(t)){
1.148 + //There are no k paths from s to t
1.149 + break;
1.150 + };
1.151 + */
1.152 +
1.153 + //We have to change the potential
1.154 + FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
1.155 + potential[n] += dijkstra.distMap()[n];
1.156 + }
1.157 +
1.158 +
1.159 + //Augmenting on the sortest path
1.160 + Node n=t;
1.161 + ResGraphEdge e;
1.162 + while (n!=s){
1.163 + e = dijkstra.pred(n);
1.164 + n = dijkstra.predNode(n);
1.165 + res_graph.augment(e,delta);
1.166 + /*
1.167 + //Let's update the total length
1.168 + if (res_graph.forward(e))
1.169 + total_length += length[e];
1.170 + else
1.171 + total_length -= length[e];
1.172 + */
1.173 + }
1.174 +
1.175 + //Update the excess_nodes heap
1.176 + if (delta >= excess_nodes[s]){
1.177 + if (delta > excess_nodes[s])
1.178 + deficit_nodes.push(s,delta - excess_nodes[s]);
1.179 + excess_nodes.pop();
1.180 +
1.181 + }
1.182 + else{
1.183 + excess_nodes[s] -= delta;
1.184 + }
1.185 + //Update the deficit_nodes heap
1.186 + if (delta >= deficit_nodes[t]){
1.187 + if (delta > deficit_nodes[t])
1.188 + excess_nodes.push(t,delta - deficit_nodes[t]);
1.189 + deficit_nodes.pop();
1.190 +
1.191 + }
1.192 + else{
1.193 + deficit_nodes[t] -= delta;
1.194 + }
1.195 + //Dijkstra part ends here
1.196 }
1.197
1.198 /*
1.199 - {
1.200 + * End of the delta scaling phase
1.201 + */
1.202
1.203 - typename ResGraphType::NodeIt n;
1.204 - for ( res_graph.first(n) ; res_graph.valid(n) ; res_graph.next(n) ) {
1.205 - potential[n] += dijkstra.distMap()[n];
1.206 + //Whatever this means
1.207 + delta = delta / 2;
1.208 +
1.209 + /*This is not necessary here
1.210 + //Update the max_excess
1.211 + max_excess = 0;
1.212 + FOR_EACH_LOC(typename Graph::NodeIt, n, G){
1.213 + if (max_excess < excess_deficit[n]){
1.214 + max_excess = excess_deficit[n];
1.215 }
1.216 }
1.217 */
1.218 -
1.219 - //Augmenting on the sortest path
1.220 - Node n=t;
1.221 - ResGraphEdge e;
1.222 - while (n!=s){
1.223 - e = dijkstra.pred(n);
1.224 - n = dijkstra.predNode(n);
1.225 - res_graph.augment(e,delta);
1.226 - //Let's update the total length
1.227 - if (res_graph.forward(e))
1.228 - total_length += length[e];
1.229 - else
1.230 - total_length -= length[e];
1.231 + //Reset delta if still too big
1.232 + if (8*number_of_nodes*max_excess <= delta){
1.233 + delta = max_excess;
1.234 +
1.235 }
1.236 -
1.237
1.238 - }
1.239 + }//while(max_excess > 0)
1.240
1.241
1.242 return i;