athos@610: // -*- c++ -*- alpar@921: #ifndef LEMON_MINCOSTFLOW_H alpar@921: #define LEMON_MINCOSTFLOW_H athos@610: athos@610: ///\ingroup galgs athos@610: ///\file athos@645: ///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network athos@611: alpar@921: #include alpar@921: #include alpar@921: #include athos@610: #include athos@657: #include athos@662: #include alpar@921: #include alpar@921: #include alpar@921: #include athos@662: #include athos@610: alpar@921: namespace lemon { athos@610: athos@610: /// \addtogroup galgs athos@610: /// @{ athos@610: athos@661: ///\brief Implementation of an algorithm for solving the minimum cost general athos@661: /// flow problem in an uncapacitated network athos@610: /// athos@610: /// alpar@921: /// The class \ref lemon::MinCostFlow "MinCostFlow" implements athos@633: /// an algorithm for solving the following general minimum cost flow problem> athos@633: /// athos@633: /// athos@633: /// athos@633: /// \warning It is assumed here that the problem has a feasible solution athos@633: /// athos@661: /// The range of the cost (weight) function is nonnegative reals but athos@610: /// the range of capacity function is the set of nonnegative integers. athos@610: /// It is not a polinomial time algorithm for counting the minimum cost athos@610: /// maximal flow, since it counts the minimum cost flow for every value 0..M athos@610: /// where \c M is the value of the maximal flow. athos@610: /// athos@610: ///\author Attila Bernath athos@661: template athos@633: class MinCostFlow { athos@610: alpar@987: typedef typename CostMap::Value Cost; athos@610: athos@633: alpar@987: typedef typename SupplyDemandMap::Value SupplyDemand; athos@610: athos@610: typedef typename Graph::Node Node; athos@610: typedef typename Graph::NodeIt NodeIt; athos@610: typedef typename Graph::Edge Edge; athos@610: typedef typename Graph::OutEdgeIt OutEdgeIt; athos@661: typedef typename Graph::template EdgeMap FlowMap; athos@661: typedef ConstMap ConstEdgeMap; athos@610: athos@610: // typedef ConstMap ConstMap; athos@610: athos@661: typedef ResGraphWrapper ResGraph; athos@661: typedef typename ResGraph::Edge ResGraphEdge; athos@610: athos@661: class ModCostMap { athos@661: //typedef typename ResGraph::template NodeMap NodeMap; athos@661: typedef typename Graph::template NodeMap NodeMap; athos@661: const ResGraph& res_graph; athos@610: // const EdgeIntMap& rev; athos@661: const CostMap &ol; athos@610: const NodeMap &pot; athos@610: public : alpar@987: typedef typename CostMap::Key Key; alpar@987: typedef typename CostMap::Value Value; athos@610: alpar@987: Value operator[](typename ResGraph::Edge e) const { athos@659: if (res_graph.forward(e)) alpar@986: return ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]); athos@610: else alpar@986: return -ol[e]-(pot[res_graph.target(e)]-pot[res_graph.source(e)]); athos@610: } athos@610: athos@661: ModCostMap(const ResGraph& _res_graph, athos@661: const CostMap &o, const NodeMap &p) : athos@659: res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){}; athos@661: };//ModCostMap athos@610: athos@610: athos@610: protected: athos@610: athos@610: //Input athos@659: const Graph& graph; athos@661: const CostMap& cost; athos@635: const SupplyDemandMap& supply_demand;//supply or demand of nodes athos@610: athos@610: athos@610: //auxiliary variables athos@610: athos@610: //To store the flow athos@661: FlowMap flow; athos@662: //To store the potential (dual variables) athos@662: typedef typename Graph::template NodeMap PotentialMap; athos@662: PotentialMap potential; athos@610: athos@610: athos@661: Cost total_cost; athos@610: athos@610: athos@610: public : athos@610: athos@610: athos@662: MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand): athos@662: graph(_graph), athos@662: cost(_cost), athos@662: supply_demand(_supply_demand), athos@662: flow(_graph), athos@672: potential(_graph){ } athos@610: athos@610: athos@610: ///Runs the algorithm. athos@610: athos@610: ///Runs the algorithm. athos@635: athos@610: ///\todo May be it does make sense to be able to start with a nonzero athos@610: /// feasible primal-dual solution pair as well. athos@659: void run() { athos@610: athos@672: //To store excess-deficit values athos@672: SupplyDemandMap excess_deficit(graph); athos@672: athos@610: //Resetting variables from previous runs athos@661: //total_cost = 0; athos@635: athos@672: athos@635: typedef typename Graph::template NodeMap HeapMap; athos@662: typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap, athos@635: std::greater > HeapType; athos@635: athos@635: //A heap for the excess nodes athos@659: HeapMap excess_nodes_map(graph,-1); athos@635: HeapType excess_nodes(excess_nodes_map); athos@635: athos@635: //A heap for the deficit nodes athos@659: HeapMap deficit_nodes_map(graph,-1); athos@635: HeapType deficit_nodes(deficit_nodes_map); athos@635: athos@657: //A container to store nonabundant arcs athos@662: std::list nonabundant_arcs; athos@659: athos@659: athos@659: FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ athos@610: flow.set(e,0); athos@657: nonabundant_arcs.push_back(e); athos@610: } athos@633: athos@633: //Initial value for delta athos@635: SupplyDemand delta = 0; athos@635: athos@657: typedef UnionFindEnum UFE; athos@657: athos@657: //A union-find structure to store the abundant components athos@662: typename UFE::MapType abund_comp_map(graph); athos@657: UFE abundant_components(abund_comp_map); athos@657: athos@657: athos@657: athos@659: FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ athos@635: excess_deficit.set(n,supply_demand[n]); athos@635: //A supply node athos@635: if (excess_deficit[n] > 0){ athos@635: excess_nodes.push(n,excess_deficit[n]); athos@633: } athos@635: //A demand node athos@635: if (excess_deficit[n] < 0){ athos@635: deficit_nodes.push(n, - excess_deficit[n]); athos@635: } athos@635: //Finding out starting value of delta athos@635: if (delta < abs(excess_deficit[n])){ athos@635: delta = abs(excess_deficit[n]); athos@635: } athos@633: //Initialize the copy of the Dijkstra potential to zero athos@610: potential.set(n,0); athos@657: //Every single point is an abundant component initially athos@657: abundant_components.insert(n); athos@610: } athos@610: athos@635: //It'll be allright as an initial value, though this value athos@635: //can be the maximum deficit here athos@635: SupplyDemand max_excess = delta; athos@610: athos@661: ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph athos@662: ConstEdgeMap const_inf_map(MAXINT); athos@661: athos@633: //We need a residual graph which is uncapacitated athos@661: ResGraph res_graph(graph, const_inf_map, flow); athos@659: athos@659: //An EdgeMap to tell which arcs are abundant athos@662: typename Graph::template EdgeMap abundant_arcs(graph); athos@610: athos@659: //Let's construct the sugraph consisting only of the abundant edges athos@659: typedef ConstMap< typename Graph::Node, bool > ConstNodeMap; athos@672: athos@659: ConstNodeMap const_true_map(true); athos@662: typedef SubGraphWrapper< const Graph, ConstNodeMap, athos@662: typename Graph::template EdgeMap > athos@659: AbundantGraph; athos@659: AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs ); athos@659: athos@659: //Let's construct the residual graph for the abundant graph athos@662: typedef ResGraphWrapper athos@659: ResAbGraph; athos@659: //Again uncapacitated athos@661: ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow); athos@659: athos@659: //We need things for the bfs athos@662: typename ResAbGraph::template NodeMap bfs_reached(res_ab_graph); athos@662: typename ResAbGraph::template NodeMap athos@659: bfs_pred(res_ab_graph); athos@662: NullMap bfs_dist_dummy; athos@671: //Teszt celbol: athos@671: //BfsIterator > athos@671: //izebize(res_ab_graph, bfs_reached); athos@659: //We want to run bfs-es (more) on this graph 'res_ab_graph' athos@671: Bfs < const ResAbGraph , athos@662: typename ResAbGraph::template NodeMap, athos@662: typename ResAbGraph::template NodeMap, athos@659: NullMap > athos@659: bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy); athos@662: /*This is what Marci wants for a bfs athos@662: template , athos@662: typename PredMap athos@662: =typename Graph::template NodeMap, athos@662: typename DistMap=typename Graph::template NodeMap > athos@662: class Bfs : public BfsIterator { athos@662: athos@662: */ athos@610: athos@661: ModCostMap mod_cost(res_graph, cost, potential); athos@610: athos@661: Dijkstra dijkstra(res_graph, mod_cost); athos@610: athos@671: //We will use the number of the nodes of the graph often athos@671: int number_of_nodes = graph.nodeNum(); athos@633: athos@635: while (max_excess > 0){ athos@635: athos@657: //Reset delta if still too big athos@657: if (8*number_of_nodes*max_excess <= delta){ athos@657: delta = max_excess; athos@657: athos@657: } athos@657: athos@645: /* athos@645: * Beginning of the delta scaling phase athos@645: */ athos@635: //Merge and stuff athos@657: { athos@657: SupplyDemand buf=8*number_of_nodes*delta; athos@662: typename std::list::iterator i = nonabundant_arcs.begin(); athos@657: while ( i != nonabundant_arcs.end() ){ athos@671: if (flow[*i]>=buf){ alpar@986: Node a = abundant_components.find(res_graph.target(*i)); alpar@986: Node b = abundant_components.find(res_graph.source(*i)); athos@657: //Merge athos@657: if (a != b){ athos@657: abundant_components.join(a,b); athos@659: //We want to push the smaller athos@659: //Which has greater absolut value excess/deficit athos@659: Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b; athos@659: //Which is the other athos@659: Node non_root = ( a == root ) ? b : a ; athos@659: abundant_components.makeRep(root); athos@659: SupplyDemand qty_to_augment = abs(excess_deficit[non_root]); athos@659: //Push the positive value athos@659: if (excess_deficit[non_root] < 0) athos@659: swap(root, non_root); athos@659: //If the non_root node has excess/deficit at all athos@659: if (qty_to_augment>0){ athos@659: //Find path and augment athos@671: bfs.run(typename AbundantGraph::Node(non_root)); athos@659: //root should be reached athos@659: athos@659: //Augmenting on the found path athos@659: Node n=root; athos@659: ResGraphEdge e; athos@659: while (n!=non_root){ athos@671: e = bfs_pred[n]; alpar@986: n = res_graph.source(e); athos@659: res_graph.augment(e,qty_to_augment); athos@659: } athos@659: athos@659: //We know that non_root had positive excess athos@671: excess_nodes.set(non_root, athos@671: excess_nodes[non_root] - qty_to_augment); athos@659: //But what about root node athos@659: //It might have been positive and so became larger athos@659: if (excess_deficit[root]>0){ athos@671: excess_nodes.set(root, athos@671: excess_nodes[root] + qty_to_augment); athos@659: } athos@659: else{ athos@659: //Or negative but not turned into positive athos@671: deficit_nodes.set(root, athos@671: deficit_nodes[root] - qty_to_augment); athos@659: } athos@659: athos@659: //Update the excess_deficit map athos@659: excess_deficit[non_root] -= qty_to_augment; athos@659: excess_deficit[root] += qty_to_augment; athos@659: athos@659: athos@659: } athos@657: } athos@657: //What happens to i? athos@659: //Marci and Zsolt says I shouldn't do such things athos@659: nonabundant_arcs.erase(i++); athos@671: abundant_arcs[*i] = true; athos@657: } athos@657: else athos@657: ++i; athos@657: } athos@657: } athos@657: athos@635: athos@635: Node s = excess_nodes.top(); athos@672: max_excess = excess_nodes[s]; athos@635: Node t = deficit_nodes.top(); athos@659: if (max_excess < deficit_nodes[t]){ athos@659: max_excess = deficit_nodes[t]; athos@635: } athos@635: athos@635: athos@662: while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){ athos@659: athos@635: athos@635: //s es t valasztasa athos@659: athos@635: //Dijkstra part athos@635: dijkstra.run(s); athos@659: athos@635: /*We know from theory that t can be reached athos@635: if (!dijkstra.reached(t)){ athos@635: //There are no k paths from s to t athos@635: break; athos@635: }; athos@635: */ athos@635: athos@635: //We have to change the potential athos@661: FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){ athos@635: potential[n] += dijkstra.distMap()[n]; athos@635: } athos@635: athos@635: athos@635: //Augmenting on the sortest path athos@635: Node n=t; athos@635: ResGraphEdge e; athos@635: while (n!=s){ athos@635: e = dijkstra.pred(n); athos@635: n = dijkstra.predNode(n); athos@635: res_graph.augment(e,delta); athos@635: /* athos@661: //Let's update the total cost athos@635: if (res_graph.forward(e)) athos@661: total_cost += cost[e]; athos@635: else athos@661: total_cost -= cost[e]; athos@635: */ athos@635: } athos@659: athos@659: //Update the excess_deficit map athos@659: excess_deficit[s] -= delta; athos@659: excess_deficit[t] += delta; athos@659: athos@635: athos@635: //Update the excess_nodes heap athos@672: if (delta > excess_nodes[s]){ athos@635: if (delta > excess_nodes[s]) athos@635: deficit_nodes.push(s,delta - excess_nodes[s]); athos@635: excess_nodes.pop(); athos@635: athos@635: } athos@635: else{ athos@671: excess_nodes.set(s, excess_nodes[s] - delta); athos@635: } athos@635: //Update the deficit_nodes heap athos@672: if (delta > deficit_nodes[t]){ athos@635: if (delta > deficit_nodes[t]) athos@635: excess_nodes.push(t,delta - deficit_nodes[t]); athos@635: deficit_nodes.pop(); athos@635: athos@635: } athos@635: else{ athos@671: deficit_nodes.set(t, deficit_nodes[t] - delta); athos@635: } athos@635: //Dijkstra part ends here athos@659: athos@659: //Choose s and t again athos@659: s = excess_nodes.top(); athos@659: max_excess = excess_nodes[s]; athos@659: t = deficit_nodes.top(); athos@659: if (max_excess < deficit_nodes[t]){ athos@659: max_excess = deficit_nodes[t]; athos@659: } athos@659: athos@633: } athos@633: athos@633: /* athos@635: * End of the delta scaling phase athos@635: */ athos@633: athos@635: //Whatever this means athos@635: delta = delta / 2; athos@635: athos@635: /*This is not necessary here athos@635: //Update the max_excess athos@635: max_excess = 0; athos@659: FOR_EACH_LOC(typename Graph::NodeIt, n, graph){ athos@635: if (max_excess < excess_deficit[n]){ athos@635: max_excess = excess_deficit[n]; athos@610: } athos@610: } athos@633: */ athos@657: athos@610: athos@635: }//while(max_excess > 0) athos@610: athos@610: athos@671: //return i; athos@610: } athos@610: athos@610: athos@610: athos@610: athos@661: ///This function gives back the total cost of the found paths. athos@610: ///Assumes that \c run() has been run and nothing changed since then. athos@661: Cost totalCost(){ athos@661: return total_cost; athos@610: } athos@610: athos@610: ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must athos@610: ///be called before using this function. athos@662: const FlowMap &getFlow() const { return flow;} athos@610: athos@610: ///Returns a const reference to the NodeMap \c potential (the dual solution). athos@610: /// \pre \ref run() must be called before using this function. athos@662: const PotentialMap &getPotential() const { return potential;} athos@610: athos@610: ///This function checks, whether the given solution is optimal athos@610: ///Running after a \c run() should return with true athos@672: ///In this "state of the art" this only checks optimality, doesn't bother with feasibility athos@610: /// athos@610: ///\todo Is this OK here? athos@610: bool checkComplementarySlackness(){ athos@661: Cost mod_pot; athos@661: Cost fl_e; athos@659: FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){ athos@610: //C^{\Pi}_{i,j} alpar@986: mod_pot = cost[e]-potential[graph.target(e)]+potential[graph.source(e)]; athos@610: fl_e = flow[e]; athos@610: // std::cout << fl_e << std::endl; athos@672: if (mod_pot > 0 && fl_e != 0) athos@672: return false; athos@672: athos@610: } athos@610: return true; athos@610: } athos@672: athos@672: /* athos@672: //For testing purposes only athos@672: //Lists the node_properties athos@672: void write_property_vector(const SupplyDemandMap& a, athos@672: char* prop_name="property"){ athos@672: FOR_EACH_LOC(typename Graph::NodeIt, i, graph){ athos@672: cout<<"Node id.: "<