athos@610: // -*- c++ -*- athos@633: #ifndef HUGO_MINCOSTFLOW_H athos@633: #define HUGO_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: athos@610: #include <hugo/dijkstra.h> athos@610: #include <hugo/graph_wrapper.h> athos@610: #include <hugo/maps.h> athos@610: #include <vector> athos@657: #include <list> athos@610: #include <for_each_macros.h> athos@657: #include <hugo/union_find.h> athos@610: athos@610: namespace hugo { athos@610: athos@610: /// \addtogroup galgs athos@610: /// @{ athos@610: athos@645: ///\brief Implementation of an algorithm for finding the minimum cost flow athos@645: /// of given value in an uncapacitated network athos@610: /// athos@610: /// athos@633: /// The class \ref hugo::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@610: /// The range of the length (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@635: template <typename Graph, typename LengthMap, typename SupplyDemandMap> athos@633: class MinCostFlow { athos@610: athos@610: typedef typename LengthMap::ValueType Length; athos@610: athos@633: athos@635: typedef typename SupplyDemandMap::ValueType 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@610: typedef typename Graph::template EdgeMap<int> EdgeIntMap; athos@610: athos@610: // typedef ConstMap<Edge,int> ConstMap; athos@610: athos@610: typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType; athos@610: typedef typename ResGraphType::Edge ResGraphEdge; athos@610: athos@610: class ModLengthMap { athos@610: //typedef typename ResGraphType::template NodeMap<Length> NodeMap; athos@610: typedef typename Graph::template NodeMap<Length> NodeMap; athos@610: const ResGraphType& G; athos@610: // const EdgeIntMap& rev; athos@610: const LengthMap &ol; athos@610: const NodeMap &pot; athos@610: public : athos@610: typedef typename LengthMap::KeyType KeyType; athos@610: typedef typename LengthMap::ValueType ValueType; athos@610: athos@610: ValueType operator[](typename ResGraphType::Edge e) const { athos@610: if (G.forward(e)) athos@610: return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]); athos@610: else athos@610: return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]); athos@610: } athos@610: athos@610: ModLengthMap(const ResGraphType& _G, athos@610: const LengthMap &o, const NodeMap &p) : athos@610: G(_G), /*rev(_rev),*/ ol(o), pot(p){}; athos@610: };//ModLengthMap athos@610: athos@610: athos@610: protected: athos@610: athos@610: //Input athos@610: const Graph& G; athos@610: const LengthMap& length; 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@610: EdgeIntMap flow; athos@610: //To store the potentila (dual variables) athos@610: typename Graph::template NodeMap<Length> potential; athos@633: //To store excess-deficit values athos@635: SupplyDemandMap excess_deficit; athos@610: athos@610: athos@610: Length total_length; athos@610: athos@610: athos@610: public : athos@610: athos@610: athos@635: MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G), athos@635: length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ } 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@633: int run() { athos@610: athos@610: //Resetting variables from previous runs athos@635: //total_length = 0; athos@635: athos@635: typedef typename Graph::template NodeMap<int> HeapMap; athos@657: typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>, athos@635: std::greater<SupplyDemand> > HeapType; athos@635: athos@635: //A heap for the excess nodes athos@635: HeapMap excess_nodes_map(G,-1); athos@635: HeapType excess_nodes(excess_nodes_map); athos@635: athos@635: //A heap for the deficit nodes athos@635: HeapMap deficit_nodes_map(G,-1); athos@635: HeapType deficit_nodes(deficit_nodes_map); athos@635: athos@657: //A container to store nonabundant arcs athos@657: list<Edge> nonabundant_arcs; athos@610: athos@610: FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ 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<Node, Graph::template NodeMap> UFE; athos@657: athos@657: //A union-find structure to store the abundant components athos@657: UFE::MapType abund_comp_map(G); athos@657: UFE abundant_components(abund_comp_map); athos@657: athos@657: athos@657: athos@610: FOR_EACH_LOC(typename Graph::NodeIt, n, G){ 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@633: //We need a residual graph which is uncapacitated athos@633: ResGraphType res_graph(G, flow); athos@610: athos@633: athos@610: athos@610: ModLengthMap mod_length(res_graph, length, potential); athos@610: athos@610: Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length); athos@610: 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@657: list<Edge>::iterator i = nonabundant_arcs.begin(); athos@657: while ( i != nonabundant_arcs.end() ){ athos@657: if (flow[i]>=buf){ athos@657: Node a = abundant_components.find(res_graph.head(i)); athos@657: Node b = abundant_components.find(res_graph.tail(i)); athos@657: //Merge athos@657: if (a != b){ athos@657: //Find path and augment athos@657: //!!! athos@657: //Find path and augment athos@657: abundant_components.join(a,b); athos@657: } athos@657: //What happens to i? athos@657: nonabundant_arcs.erase(i); athos@657: } athos@657: else athos@657: ++i; athos@657: } athos@657: } athos@657: athos@635: athos@635: Node s = excess_nodes.top(); athos@635: SupplyDemand max_excess = excess_nodes[s]; athos@635: Node t = deficit_nodes.top(); athos@635: if (max_excess < dificit_nodes[t]){ athos@635: max_excess = dificit_nodes[t]; athos@635: } athos@635: athos@635: athos@645: while(max_excess > 0){ athos@635: athos@635: athos@635: //s es t valasztasa athos@635: athos@635: //Dijkstra part athos@635: dijkstra.run(s); athos@635: 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@635: FOR_EACH_LOC(typename ResGraphType::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@635: //Let's update the total length athos@635: if (res_graph.forward(e)) athos@635: total_length += length[e]; athos@635: else athos@635: total_length -= length[e]; athos@635: */ athos@635: } athos@635: athos@635: //Update the excess_nodes heap athos@635: 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@635: excess_nodes[s] -= delta; athos@635: } athos@635: //Update the deficit_nodes heap athos@635: 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@635: deficit_nodes[t] -= delta; athos@635: } athos@635: //Dijkstra part ends here 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@635: FOR_EACH_LOC(typename Graph::NodeIt, n, G){ 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@610: return i; athos@610: } athos@610: athos@610: athos@610: athos@610: athos@610: ///This function gives back the total length of the found paths. athos@610: ///Assumes that \c run() has been run and nothing changed since then. athos@610: Length totalLength(){ athos@610: return total_length; 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@610: const EdgeIntMap &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@610: const EdgeIntMap &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@610: ///In this "state of the art" this only check optimality, doesn't bother with feasibility athos@610: /// athos@610: ///\todo Is this OK here? athos@610: bool checkComplementarySlackness(){ athos@610: Length mod_pot; athos@610: Length fl_e; athos@610: FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ athos@610: //C^{\Pi}_{i,j} athos@610: mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)]; athos@610: fl_e = flow[e]; athos@610: // std::cout << fl_e << std::endl; athos@610: if (0<fl_e && fl_e<capacity[e]){ athos@610: if (mod_pot != 0) athos@610: return false; athos@610: } athos@610: else{ athos@610: if (mod_pot > 0 && fl_e != 0) athos@610: return false; athos@610: if (mod_pot < 0 && fl_e != capacity[e]) athos@610: return false; athos@610: } athos@610: } athos@610: return true; athos@610: } athos@610: athos@610: athos@633: }; //class MinCostFlow athos@610: athos@610: ///@} athos@610: athos@610: } //namespace hugo athos@610: athos@610: #endif //HUGO_MINCOSTFLOW_H