// -*- c++ -*- #ifndef HUGO_MINCOSTFLOW_H #define HUGO_MINCOSTFLOW_H ///\ingroup galgs ///\file ///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network #include #include #include #include #include namespace hugo { /// \addtogroup galgs /// @{ ///\brief Implementation of an algorithm for finding the minimum cost flow /// of given value in an uncapacitated network /// /// /// The class \ref hugo::MinCostFlow "MinCostFlow" implements /// an algorithm for solving the following general minimum cost flow problem> /// /// /// /// \warning It is assumed here that the problem has a feasible solution /// /// The range of the length (weight) function is nonnegative reals but /// the range of capacity function is the set of nonnegative integers. /// It is not a polinomial time algorithm for counting the minimum cost /// maximal flow, since it counts the minimum cost flow for every value 0..M /// where \c M is the value of the maximal flow. /// ///\author Attila Bernath template class MinCostFlow { typedef typename LengthMap::ValueType Length; typedef typename SupplyDemandMap::ValueType SupplyDemand; typedef typename Graph::Node Node; typedef typename Graph::NodeIt NodeIt; typedef typename Graph::Edge Edge; typedef typename Graph::OutEdgeIt OutEdgeIt; typedef typename Graph::template EdgeMap EdgeIntMap; // typedef ConstMap ConstMap; typedef ResGraphWrapper ResGraphType; typedef typename ResGraphType::Edge ResGraphEdge; class ModLengthMap { //typedef typename ResGraphType::template NodeMap NodeMap; typedef typename Graph::template NodeMap NodeMap; const ResGraphType& G; // const EdgeIntMap& rev; const LengthMap &ol; const NodeMap &pot; public : typedef typename LengthMap::KeyType KeyType; typedef typename LengthMap::ValueType ValueType; ValueType operator[](typename ResGraphType::Edge e) const { if (G.forward(e)) return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]); else return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]); } ModLengthMap(const ResGraphType& _G, const LengthMap &o, const NodeMap &p) : G(_G), /*rev(_rev),*/ ol(o), pot(p){}; };//ModLengthMap protected: //Input const Graph& G; const LengthMap& length; const SupplyDemandMap& supply_demand;//supply or demand of nodes //auxiliary variables //To store the flow EdgeIntMap flow; //To store the potentila (dual variables) typename Graph::template NodeMap potential; //To store excess-deficit values SupplyDemandMap excess_deficit; Length total_length; public : MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G), length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ } ///Runs the algorithm. ///Runs the algorithm. ///\todo May be it does make sense to be able to start with a nonzero /// feasible primal-dual solution pair as well. int run() { //Resetting variables from previous runs //total_length = 0; typedef typename Graph::template NodeMap HeapMap; typedef Heap, std::greater > HeapType; //A heap for the excess nodes HeapMap excess_nodes_map(G,-1); HeapType excess_nodes(excess_nodes_map); //A heap for the deficit nodes HeapMap deficit_nodes_map(G,-1); HeapType deficit_nodes(deficit_nodes_map); FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ flow.set(e,0); } //Initial value for delta SupplyDemand delta = 0; FOR_EACH_LOC(typename Graph::NodeIt, n, G){ excess_deficit.set(n,supply_demand[n]); //A supply node if (excess_deficit[n] > 0){ excess_nodes.push(n,excess_deficit[n]); } //A demand node if (excess_deficit[n] < 0){ deficit_nodes.push(n, - excess_deficit[n]); } //Finding out starting value of delta if (delta < abs(excess_deficit[n])){ delta = abs(excess_deficit[n]); } //Initialize the copy of the Dijkstra potential to zero potential.set(n,0); } //It'll be allright as an initial value, though this value //can be the maximum deficit here SupplyDemand max_excess = delta; //We need a residual graph which is uncapacitated ResGraphType res_graph(G, flow); ModLengthMap mod_length(res_graph, length, potential); Dijkstra dijkstra(res_graph, mod_length); while (max_excess > 0){ /* * Beginning of the delta scaling phase */ //Merge and stuff Node s = excess_nodes.top(); SupplyDemand max_excess = excess_nodes[s]; Node t = deficit_nodes.top(); if (max_excess < dificit_nodes[t]){ max_excess = dificit_nodes[t]; } while(max_excess > 0){ //s es t valasztasa //Dijkstra part dijkstra.run(s); /*We know from theory that t can be reached if (!dijkstra.reached(t)){ //There are no k paths from s to t break; }; */ //We have to change the potential FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){ potential[n] += dijkstra.distMap()[n]; } //Augmenting on the sortest path Node n=t; ResGraphEdge e; while (n!=s){ e = dijkstra.pred(n); n = dijkstra.predNode(n); res_graph.augment(e,delta); /* //Let's update the total length if (res_graph.forward(e)) total_length += length[e]; else total_length -= length[e]; */ } //Update the excess_nodes heap if (delta >= excess_nodes[s]){ if (delta > excess_nodes[s]) deficit_nodes.push(s,delta - excess_nodes[s]); excess_nodes.pop(); } else{ excess_nodes[s] -= delta; } //Update the deficit_nodes heap if (delta >= deficit_nodes[t]){ if (delta > deficit_nodes[t]) excess_nodes.push(t,delta - deficit_nodes[t]); deficit_nodes.pop(); } else{ deficit_nodes[t] -= delta; } //Dijkstra part ends here } /* * End of the delta scaling phase */ //Whatever this means delta = delta / 2; /*This is not necessary here //Update the max_excess max_excess = 0; FOR_EACH_LOC(typename Graph::NodeIt, n, G){ if (max_excess < excess_deficit[n]){ max_excess = excess_deficit[n]; } } */ //Reset delta if still too big if (8*number_of_nodes*max_excess <= delta){ delta = max_excess; } }//while(max_excess > 0) return i; } ///This function gives back the total length of the found paths. ///Assumes that \c run() has been run and nothing changed since then. Length totalLength(){ return total_length; } ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must ///be called before using this function. const EdgeIntMap &getFlow() const { return flow;} ///Returns a const reference to the NodeMap \c potential (the dual solution). /// \pre \ref run() must be called before using this function. const EdgeIntMap &getPotential() const { return potential;} ///This function checks, whether the given solution is optimal ///Running after a \c run() should return with true ///In this "state of the art" this only check optimality, doesn't bother with feasibility /// ///\todo Is this OK here? bool checkComplementarySlackness(){ Length mod_pot; Length fl_e; FOR_EACH_LOC(typename Graph::EdgeIt, e, G){ //C^{\Pi}_{i,j} mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)]; fl_e = flow[e]; // std::cout << fl_e << std::endl; if (0 0 && fl_e != 0) return false; if (mod_pot < 0 && fl_e != capacity[e]) return false; } } return true; } }; //class MinCostFlow ///@} } //namespace hugo #endif //HUGO_MINCOSTFLOW_H