1.1 --- a/src/work/athos/mincostflow.h Sun Apr 17 18:57:22 2005 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,513 +0,0 @@
1.4 -// -*- c++ -*-
1.5 -#ifndef LEMON_MINCOSTFLOW_H
1.6 -#define LEMON_MINCOSTFLOW_H
1.7 -
1.8 -///\ingroup galgs
1.9 -///\file
1.10 -///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network
1.11 -
1.12 -#include <lemon/dijkstra.h>
1.13 -#include <lemon/graph_wrapper.h>
1.14 -#include <lemon/maps.h>
1.15 -#include <vector>
1.16 -#include <list>
1.17 -#include <values.h>
1.18 -#include <lemon/for_each_macros.h>
1.19 -#include <lemon/unionfind.h>
1.20 -#include <lemon/bin_heap.h>
1.21 -#include <bfs_dfs.h>
1.22 -
1.23 -namespace lemon {
1.24 -
1.25 -/// \addtogroup galgs
1.26 -/// @{
1.27 -
1.28 - ///\brief Implementation of an algorithm for solving the minimum cost general
1.29 - /// flow problem in an uncapacitated network
1.30 - ///
1.31 - ///
1.32 - /// The class \ref lemon::MinCostFlow "MinCostFlow" implements
1.33 - /// an algorithm for solving the following general minimum cost flow problem>
1.34 - ///
1.35 - ///
1.36 - ///
1.37 - /// \warning It is assumed here that the problem has a feasible solution
1.38 - ///
1.39 - /// The range of the cost (weight) function is nonnegative reals but
1.40 - /// the range of capacity function is the set of nonnegative integers.
1.41 - /// It is not a polinomial time algorithm for counting the minimum cost
1.42 - /// maximal flow, since it counts the minimum cost flow for every value 0..M
1.43 - /// where \c M is the value of the maximal flow.
1.44 - ///
1.45 - ///\author Attila Bernath
1.46 - template <typename Graph, typename CostMap, typename SupplyDemandMap>
1.47 - class MinCostFlow {
1.48 -
1.49 - typedef typename CostMap::Value Cost;
1.50 -
1.51 -
1.52 - typedef typename SupplyDemandMap::Value SupplyDemand;
1.53 -
1.54 - typedef typename Graph::Node Node;
1.55 - typedef typename Graph::NodeIt NodeIt;
1.56 - typedef typename Graph::Edge Edge;
1.57 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.58 - typedef typename Graph::template EdgeMap<SupplyDemand> FlowMap;
1.59 - typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap;
1.60 -
1.61 - // typedef ConstMap<Edge,int> ConstMap;
1.62 -
1.63 - typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
1.64 - typedef typename ResGraph::Edge ResGraphEdge;
1.65 -
1.66 - class ModCostMap {
1.67 - //typedef typename ResGraph::template NodeMap<Cost> NodeMap;
1.68 - typedef typename Graph::template NodeMap<Cost> NodeMap;
1.69 - const ResGraph& res_graph;
1.70 - // const EdgeIntMap& rev;
1.71 - const CostMap &ol;
1.72 - const NodeMap &pot;
1.73 - public :
1.74 - typedef typename CostMap::Key Key;
1.75 - typedef typename CostMap::Value Value;
1.76 -
1.77 - Value operator[](typename ResGraph::Edge e) const {
1.78 - if (res_graph.forward(e))
1.79 - return ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]);
1.80 - else
1.81 - return -ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]);
1.82 - }
1.83 -
1.84 - ModCostMap(const ResGraph& _res_graph,
1.85 - const CostMap &o, const NodeMap &p) :
1.86 - res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
1.87 - };//ModCostMap
1.88 -
1.89 -
1.90 - protected:
1.91 -
1.92 - //Input
1.93 - const Graph& graph;
1.94 - const CostMap& cost;
1.95 - const SupplyDemandMap& supply_demand;//supply or demand of nodes
1.96 -
1.97 -
1.98 - //auxiliary variables
1.99 -
1.100 - //To store the flow
1.101 - FlowMap flow;
1.102 - //To store the potential (dual variables)
1.103 - typedef typename Graph::template NodeMap<Cost> PotentialMap;
1.104 - PotentialMap potential;
1.105 -
1.106 -
1.107 - Cost total_cost;
1.108 -
1.109 -
1.110 - public :
1.111 -
1.112 -
1.113 - MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand):
1.114 - graph(_graph),
1.115 - cost(_cost),
1.116 - supply_demand(_supply_demand),
1.117 - flow(_graph),
1.118 - potential(_graph){ }
1.119 -
1.120 -
1.121 - ///Runs the algorithm.
1.122 -
1.123 - ///Runs the algorithm.
1.124 -
1.125 - ///\todo May be it does make sense to be able to start with a nonzero
1.126 - /// feasible primal-dual solution pair as well.
1.127 - void run() {
1.128 -
1.129 - //To store excess-deficit values
1.130 - SupplyDemandMap excess_deficit(graph);
1.131 -
1.132 - //Resetting variables from previous runs
1.133 - //total_cost = 0;
1.134 -
1.135 -
1.136 - typedef typename Graph::template NodeMap<int> HeapMap;
1.137 - typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
1.138 - std::greater<SupplyDemand> > HeapType;
1.139 -
1.140 - //A heap for the excess nodes
1.141 - HeapMap excess_nodes_map(graph,-1);
1.142 - HeapType excess_nodes(excess_nodes_map);
1.143 -
1.144 - //A heap for the deficit nodes
1.145 - HeapMap deficit_nodes_map(graph,-1);
1.146 - HeapType deficit_nodes(deficit_nodes_map);
1.147 -
1.148 - //A container to store nonabundant arcs
1.149 - std::list<Edge> nonabundant_arcs;
1.150 -
1.151 -
1.152 - FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
1.153 - flow.set(e,0);
1.154 - nonabundant_arcs.push_back(e);
1.155 - }
1.156 -
1.157 - //Initial value for delta
1.158 - SupplyDemand delta = 0;
1.159 -
1.160 - typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
1.161 -
1.162 - //A union-find structure to store the abundant components
1.163 - typename UFE::MapType abund_comp_map(graph);
1.164 - UFE abundant_components(abund_comp_map);
1.165 -
1.166 -
1.167 -
1.168 - FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
1.169 - excess_deficit.set(n,supply_demand[n]);
1.170 - //A supply node
1.171 - if (excess_deficit[n] > 0){
1.172 - excess_nodes.push(n,excess_deficit[n]);
1.173 - }
1.174 - //A demand node
1.175 - if (excess_deficit[n] < 0){
1.176 - deficit_nodes.push(n, - excess_deficit[n]);
1.177 - }
1.178 - //Finding out starting value of delta
1.179 - if (delta < abs(excess_deficit[n])){
1.180 - delta = abs(excess_deficit[n]);
1.181 - }
1.182 - //Initialize the copy of the Dijkstra potential to zero
1.183 - potential.set(n,0);
1.184 - //Every single point is an abundant component initially
1.185 - abundant_components.insert(n);
1.186 - }
1.187 -
1.188 - //It'll be allright as an initial value, though this value
1.189 - //can be the maximum deficit here
1.190 - SupplyDemand max_excess = delta;
1.191 -
1.192 - ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
1.193 - ConstEdgeMap const_inf_map(MAXINT);
1.194 -
1.195 - //We need a residual graph which is uncapacitated
1.196 - ResGraph res_graph(graph, const_inf_map, flow);
1.197 -
1.198 - //An EdgeMap to tell which arcs are abundant
1.199 - typename Graph::template EdgeMap<bool> abundant_arcs(graph);
1.200 -
1.201 - //Let's construct the sugraph consisting only of the abundant edges
1.202 - typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
1.203 -
1.204 - ConstNodeMap const_true_map(true);
1.205 - typedef SubGraphWrapper< const Graph, ConstNodeMap,
1.206 - typename Graph::template EdgeMap<bool> >
1.207 - AbundantGraph;
1.208 - AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
1.209 -
1.210 - //Let's construct the residual graph for the abundant graph
1.211 - typedef ResGraphWrapper<const AbundantGraph,int,ConstEdgeMap,FlowMap>
1.212 - ResAbGraph;
1.213 - //Again uncapacitated
1.214 - ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
1.215 -
1.216 - //We need things for the bfs
1.217 - typename ResAbGraph::template NodeMap<bool> bfs_reached(res_ab_graph);
1.218 - typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>
1.219 - bfs_pred(res_ab_graph);
1.220 - NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy;
1.221 - //Teszt celbol:
1.222 - //BfsIterator<ResAbGraph, typename ResAbGraph::template NodeMap<bool> >
1.223 - //izebize(res_ab_graph, bfs_reached);
1.224 - //We want to run bfs-es (more) on this graph 'res_ab_graph'
1.225 - Bfs < const ResAbGraph ,
1.226 - typename ResAbGraph::template NodeMap<bool>,
1.227 - typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>,
1.228 - NullMap<typename ResAbGraph::Node, int> >
1.229 - bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
1.230 - /*This is what Marci wants for a bfs
1.231 - template <typename Graph,
1.232 - typename ReachedMap=typename Graph::template NodeMap<bool>,
1.233 - typename PredMap
1.234 - =typename Graph::template NodeMap<typename Graph::Edge>,
1.235 - typename DistMap=typename Graph::template NodeMap<int> >
1.236 - class Bfs : public BfsIterator<Graph, ReachedMap> {
1.237 -
1.238 - */
1.239 -
1.240 - ModCostMap mod_cost(res_graph, cost, potential);
1.241 -
1.242 - Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
1.243 -
1.244 - //We will use the number of the nodes of the graph often
1.245 - int number_of_nodes = graph.nodeNum();
1.246 -
1.247 - while (max_excess > 0){
1.248 -
1.249 - //Reset delta if still too big
1.250 - if (8*number_of_nodes*max_excess <= delta){
1.251 - delta = max_excess;
1.252 -
1.253 - }
1.254 -
1.255 - /*
1.256 - * Beginning of the delta scaling phase
1.257 - */
1.258 - //Merge and stuff
1.259 - {
1.260 - SupplyDemand buf=8*number_of_nodes*delta;
1.261 - typename std::list<Edge>::iterator i = nonabundant_arcs.begin();
1.262 - while ( i != nonabundant_arcs.end() ){
1.263 - if (flow[*i]>=buf){
1.264 - Node a = abundant_components.find(res_graph.target(*i));
1.265 - Node b = abundant_components.find(res_graph.source(*i));
1.266 - //Merge
1.267 - if (a != b){
1.268 - abundant_components.join(a,b);
1.269 - //We want to push the smaller
1.270 - //Which has greater absolut value excess/deficit
1.271 - Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
1.272 - //Which is the other
1.273 - Node non_root = ( a == root ) ? b : a ;
1.274 - abundant_components.makeRep(root);
1.275 - SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
1.276 - //Push the positive value
1.277 - if (excess_deficit[non_root] < 0)
1.278 - swap(root, non_root);
1.279 - //If the non_root node has excess/deficit at all
1.280 - if (qty_to_augment>0){
1.281 - //Find path and augment
1.282 - bfs.run(typename AbundantGraph::Node(non_root));
1.283 - //root should be reached
1.284 -
1.285 - //Augmenting on the found path
1.286 - Node n=root;
1.287 - ResGraphEdge e;
1.288 - while (n!=non_root){
1.289 - e = bfs_pred[n];
1.290 - n = res_graph.source(e);
1.291 - res_graph.augment(e,qty_to_augment);
1.292 - }
1.293 -
1.294 - //We know that non_root had positive excess
1.295 - excess_nodes.set(non_root,
1.296 - excess_nodes[non_root] - qty_to_augment);
1.297 - //But what about root node
1.298 - //It might have been positive and so became larger
1.299 - if (excess_deficit[root]>0){
1.300 - excess_nodes.set(root,
1.301 - excess_nodes[root] + qty_to_augment);
1.302 - }
1.303 - else{
1.304 - //Or negative but not turned into positive
1.305 - deficit_nodes.set(root,
1.306 - deficit_nodes[root] - qty_to_augment);
1.307 - }
1.308 -
1.309 - //Update the excess_deficit map
1.310 - excess_deficit[non_root] -= qty_to_augment;
1.311 - excess_deficit[root] += qty_to_augment;
1.312 -
1.313 -
1.314 - }
1.315 - }
1.316 - //What happens to i?
1.317 - //Marci and Zsolt says I shouldn't do such things
1.318 - nonabundant_arcs.erase(i++);
1.319 - abundant_arcs[*i] = true;
1.320 - }
1.321 - else
1.322 - ++i;
1.323 - }
1.324 - }
1.325 -
1.326 -
1.327 - Node s = excess_nodes.top();
1.328 - max_excess = excess_nodes[s];
1.329 - Node t = deficit_nodes.top();
1.330 - if (max_excess < deficit_nodes[t]){
1.331 - max_excess = deficit_nodes[t];
1.332 - }
1.333 -
1.334 -
1.335 - while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){
1.336 -
1.337 -
1.338 - //s es t valasztasa
1.339 -
1.340 - //Dijkstra part
1.341 - dijkstra.run(s);
1.342 -
1.343 - /*We know from theory that t can be reached
1.344 - if (!dijkstra.reached(t)){
1.345 - //There are no k paths from s to t
1.346 - break;
1.347 - };
1.348 - */
1.349 -
1.350 - //We have to change the potential
1.351 - FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
1.352 - potential[n] += dijkstra.distMap()[n];
1.353 - }
1.354 -
1.355 -
1.356 - //Augmenting on the sortest path
1.357 - Node n=t;
1.358 - ResGraphEdge e;
1.359 - while (n!=s){
1.360 - e = dijkstra.pred(n);
1.361 - n = dijkstra.predNode(n);
1.362 - res_graph.augment(e,delta);
1.363 - /*
1.364 - //Let's update the total cost
1.365 - if (res_graph.forward(e))
1.366 - total_cost += cost[e];
1.367 - else
1.368 - total_cost -= cost[e];
1.369 - */
1.370 - }
1.371 -
1.372 - //Update the excess_deficit map
1.373 - excess_deficit[s] -= delta;
1.374 - excess_deficit[t] += delta;
1.375 -
1.376 -
1.377 - //Update the excess_nodes heap
1.378 - if (delta > excess_nodes[s]){
1.379 - if (delta > excess_nodes[s])
1.380 - deficit_nodes.push(s,delta - excess_nodes[s]);
1.381 - excess_nodes.pop();
1.382 -
1.383 - }
1.384 - else{
1.385 - excess_nodes.set(s, excess_nodes[s] - delta);
1.386 - }
1.387 - //Update the deficit_nodes heap
1.388 - if (delta > deficit_nodes[t]){
1.389 - if (delta > deficit_nodes[t])
1.390 - excess_nodes.push(t,delta - deficit_nodes[t]);
1.391 - deficit_nodes.pop();
1.392 -
1.393 - }
1.394 - else{
1.395 - deficit_nodes.set(t, deficit_nodes[t] - delta);
1.396 - }
1.397 - //Dijkstra part ends here
1.398 -
1.399 - //Choose s and t again
1.400 - s = excess_nodes.top();
1.401 - max_excess = excess_nodes[s];
1.402 - t = deficit_nodes.top();
1.403 - if (max_excess < deficit_nodes[t]){
1.404 - max_excess = deficit_nodes[t];
1.405 - }
1.406 -
1.407 - }
1.408 -
1.409 - /*
1.410 - * End of the delta scaling phase
1.411 - */
1.412 -
1.413 - //Whatever this means
1.414 - delta = delta / 2;
1.415 -
1.416 - /*This is not necessary here
1.417 - //Update the max_excess
1.418 - max_excess = 0;
1.419 - FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
1.420 - if (max_excess < excess_deficit[n]){
1.421 - max_excess = excess_deficit[n];
1.422 - }
1.423 - }
1.424 - */
1.425 -
1.426 -
1.427 - }//while(max_excess > 0)
1.428 -
1.429 -
1.430 - //return i;
1.431 - }
1.432 -
1.433 -
1.434 -
1.435 -
1.436 - ///This function gives back the total cost of the found paths.
1.437 - ///Assumes that \c run() has been run and nothing changed since then.
1.438 - Cost totalCost(){
1.439 - return total_cost;
1.440 - }
1.441 -
1.442 - ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
1.443 - ///be called before using this function.
1.444 - const FlowMap &getFlow() const { return flow;}
1.445 -
1.446 - ///Returns a const reference to the NodeMap \c potential (the dual solution).
1.447 - /// \pre \ref run() must be called before using this function.
1.448 - const PotentialMap &getPotential() const { return potential;}
1.449 -
1.450 - ///This function checks, whether the given solution is optimal
1.451 - ///Running after a \c run() should return with true
1.452 - ///In this "state of the art" this only checks optimality, doesn't bother with feasibility
1.453 - ///
1.454 - ///\todo Is this OK here?
1.455 - bool checkComplementarySlackness(){
1.456 - Cost mod_pot;
1.457 - Cost fl_e;
1.458 - FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
1.459 - //C^{\Pi}_{i,j}
1.460 - mod_pot = cost[e]-potential[graph.target(e)]+potential[graph.source(e)];
1.461 - fl_e = flow[e];
1.462 - // std::cout << fl_e << std::endl;
1.463 - if (mod_pot > 0 && fl_e != 0)
1.464 - return false;
1.465 -
1.466 - }
1.467 - return true;
1.468 - }
1.469 -
1.470 - /*
1.471 - //For testing purposes only
1.472 - //Lists the node_properties
1.473 - void write_property_vector(const SupplyDemandMap& a,
1.474 - char* prop_name="property"){
1.475 - FOR_EACH_LOC(typename Graph::NodeIt, i, graph){
1.476 - cout<<"Node id.: "<<graph.id(i)<<", "<<prop_name<<" value: "<<a[i]<<endl;
1.477 - }
1.478 - cout<<endl;
1.479 - }
1.480 - */
1.481 - bool checkFeasibility(){
1.482 - SupplyDemandMap supdem(graph);
1.483 - FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
1.484 -
1.485 - if ( flow[e] < 0){
1.486 -
1.487 - return false;
1.488 - }
1.489 - supdem[graph.source(e)] += flow[e];
1.490 - supdem[graph.target(e)] -= flow[e];
1.491 - }
1.492 - //write_property_vector(supdem, "supdem");
1.493 - //write_property_vector(supply_demand, "supply_demand");
1.494 -
1.495 - FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
1.496 -
1.497 - if ( supdem[n] != supply_demand[n]){
1.498 - //cout<<"Node id.: "<<graph.id(n)<<" : "<<supdem[n]<<", should be: "<<supply_demand[n]<<endl;
1.499 - return false;
1.500 - }
1.501 - }
1.502 -
1.503 - return true;
1.504 - }
1.505 -
1.506 - bool checkOptimality(){
1.507 - return checkFeasibility() && checkComplementarySlackness();
1.508 - }
1.509 -
1.510 - }; //class MinCostFlow
1.511 -
1.512 - ///@}
1.513 -
1.514 -} //namespace lemon
1.515 -
1.516 -#endif //LEMON_MINCOSTFLOW_H