marci@1017: // -*- c++ -*- marci@1017: #ifndef LEMON_MIN_COST_GEN_FLOW_H marci@1017: #define LEMON_MIN_COST_GEN_FLOW_H marci@1025: #include marci@1017: //#include marci@1017: marci@1025: #include marci@1025: #include marci@1017: //#include marci@1017: //#include marci@1017: //#include marci@1025: #include marci@1028: #include marci@1017: //#include marci@1017: //#include marci@1028: #include marci@1074: #include marci@1017: marci@1017: namespace lemon { marci@1017: marci@1017: template marci@1017: class PrimalMap { marci@1017: protected: marci@1081: LPGLPK* lp; marci@1017: EdgeIndexMap* edge_index_map; marci@1017: public: marci@1081: PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : marci@1017: lp(&_lp), edge_index_map(&_edge_index_map) { } marci@1017: double operator[](Edge e) const { marci@1017: return lp->getPrimal((*edge_index_map)[e]); marci@1017: } marci@1017: }; marci@1017: marci@1028: // excess: rho-delta egyelore csak =0-ra. marci@1017: template , marci@1017: typename LCapMap=typename Graph::template EdgeMap, marci@1017: typename CapMap=typename Graph::template EdgeMap, marci@1017: typename FlowMap=typename Graph::template EdgeMap, marci@1017: typename CostMap=typename Graph::template EdgeMap > marci@1017: class MinCostGenFlow { marci@1017: protected: marci@1017: const Graph& g; marci@1017: const Excess& excess; marci@1017: const LCapMap& lcapacity; marci@1017: const CapMap& capacity; marci@1017: FlowMap& flow; marci@1017: const CostMap& cost; marci@1017: public: marci@1017: MinCostGenFlow(const Graph& _g, const Excess& _excess, marci@1017: const LCapMap& _lcapacity, const CapMap& _capacity, marci@1017: FlowMap& _flow, marci@1017: const CostMap& _cost) : marci@1017: g(_g), excess(_excess), lcapacity(_lcapacity), marci@1017: capacity(_capacity), flow(_flow), cost(_cost) { } marci@1025: bool feasible() { marci@1025: // std::cout << "making new vertices..." << std::endl; marci@1025: typedef ListGraph Graph2; marci@1025: Graph2 g2; marci@1025: typedef MergeEdgeGraphWrapper GW; marci@1025: // std::cout << "merging..." << std::endl; marci@1025: GW gw(g, g2); marci@1025: typename GW::Node s(INVALID, g2.addNode(), true); marci@1025: typename GW::Node t(INVALID, g2.addNode(), true); marci@1025: typedef SmartGraph Graph3; marci@1025: // std::cout << "making extender graph..." << std::endl; marci@1025: typedef NewEdgeSetGraphWrapper2 GWW; marci@1025: // { marci@1025: // checkConcept(); marci@1025: // } marci@1025: GWW gww(gw); marci@1025: typedef AugmentingGraphWrapper GWWW; marci@1025: GWWW gwww(gw, gww); marci@1025: marci@1025: // std::cout << "making new edges..." << std::endl; marci@1025: typename GWWW::template EdgeMap translated_cap(gwww); marci@1025: marci@1028: for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) { marci@1028: translated_cap.set(typename GWWW::Edge(e,INVALID,false), marci@1028: capacity[e]-lcapacity[e]); marci@1028: // cout << "t_cap " << gw.id(e) << " " marci@1028: // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl; marci@1028: } marci@1025: marci@1025: Num expected=0; marci@1025: marci@1025: // std::cout << "making new edges 2..." << std::endl; marci@1025: for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { marci@1025: Num a=0; marci@1025: for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) marci@1025: a+=lcapacity[e]; marci@1025: for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) marci@1025: a-=lcapacity[e]; marci@1025: if (excess[n]>a) { marci@1025: typename GWW::Edge e= marci@1025: gww.addEdge(typename GW::Node(n,INVALID,false), t); marci@1025: translated_cap.set(typename GWWW::Edge(INVALID, e, true), marci@1025: excess[n]-a); marci@1028: // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl; marci@1025: } marci@1025: if (excess[n]" << g.id(n) << " "<< a-excess[n] < translated_flow(gwww, 0); marci@1025: Preflow preflow(gwww, s, t, marci@1025: translated_cap, translated_flow); marci@1025: preflow.run(); marci@1025: // std::cout << "fv: " << preflow.flowValue() << std::endl; marci@1025: // std::cout << "expected: " << expected << std::endl; marci@1025: marci@1025: for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { marci@1025: typename GW::Edge ew(e, INVALID, false); marci@1025: typename GWWW::Edge ewww(ew, INVALID, false); marci@1025: flow.set(e, translated_flow[ewww]+lcapacity[e]); marci@1025: } marci@1033: return (preflow.flowValue()>=expected); marci@1025: } marci@1028: // for nonnegative costs marci@1028: bool run() { marci@1028: // std::cout << "making new vertices..." << std::endl; marci@1028: typedef ListGraph Graph2; marci@1028: Graph2 g2; marci@1028: typedef MergeEdgeGraphWrapper GW; marci@1028: // std::cout << "merging..." << std::endl; marci@1028: GW gw(g, g2); marci@1028: typename GW::Node s(INVALID, g2.addNode(), true); marci@1028: typename GW::Node t(INVALID, g2.addNode(), true); marci@1028: typedef SmartGraph Graph3; marci@1028: // std::cout << "making extender graph..." << std::endl; marci@1028: typedef NewEdgeSetGraphWrapper2 GWW; marci@1028: // { marci@1028: // checkConcept(); marci@1028: // } marci@1028: GWW gww(gw); marci@1028: typedef AugmentingGraphWrapper GWWW; marci@1028: GWWW gwww(gw, gww); marci@1028: marci@1028: // std::cout << "making new edges..." << std::endl; marci@1028: typename GWWW::template EdgeMap translated_cap(gwww); marci@1028: marci@1028: for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { marci@1028: typename GW::Edge ew(e, INVALID, false); marci@1028: typename GWWW::Edge ewww(ew, INVALID, false); marci@1028: translated_cap.set(ewww, capacity[e]-lcapacity[e]); marci@1028: // cout << "t_cap " << g.id(e) << " " marci@1028: // << translated_cap[ewww] << endl; marci@1028: } marci@1028: marci@1028: Num expected=0; marci@1028: marci@1028: // std::cout << "making new edges 2..." << std::endl; marci@1028: for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { marci@1028: // std::cout << "node: " << g.id(n) << std::endl; marci@1028: Num a=0; marci@1028: for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) { marci@1028: a+=lcapacity[e]; marci@1028: // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl; marci@1028: } marci@1028: for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) { marci@1028: a-=lcapacity[e]; marci@1028: // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl; marci@1028: } marci@1028: // std::cout << "excess " << g.id(n) << ": " << a << std::endl; marci@1028: if (0>a) { marci@1028: typename GWW::Edge e= marci@1028: gww.addEdge(typename GW::Node(n,INVALID,false), t); marci@1028: translated_cap.set(typename GWWW::Edge(INVALID, e, true), marci@1028: -a); marci@1028: // std::cout << g.id(n) << "->t " << -a << std::endl; marci@1028: } marci@1028: if (0" << g.id(n) << " "<< a < translated_cost(gwww, 0); marci@1028: for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { marci@1028: translated_cost.set(typename GWWW::Edge( marci@1028: typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]); marci@1028: } marci@1028: // typename GWWW::template EdgeMap translated_flow(gwww, 0); marci@1028: MinCostFlow, marci@1028: typename GWWW::template EdgeMap > marci@1028: min_cost_flow(gwww, translated_cost, translated_cap, marci@1028: s, t); marci@1028: while (min_cost_flow.augment()) { } marci@1028: std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; marci@1028: std::cout << "expected: " << expected << std::endl; marci@1028: marci@1028: for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { marci@1028: typename GW::Edge ew(e, INVALID, false); marci@1028: typename GWWW::Edge ewww(ew, INVALID, false); marci@1028: // std::cout << g.id(e) << " " << flow[e] << std::endl; marci@1028: flow.set(e, lcapacity[e]+ marci@1028: min_cost_flow.getFlow()[ewww]); marci@1028: } marci@1033: return (min_cost_flow.flowValue()>=expected); marci@1028: } marci@1028: void runByLP() { marci@1081: typedef LPGLPK LPSolver; marci@1031: LPSolver lp; marci@1017: lp.setMinimize(); marci@1031: typedef LPSolver::ColIt ColIt; marci@1031: typedef LPSolver::RowIt RowIt; marci@1017: typedef typename Graph::template EdgeMap EdgeIndexMap; marci@1017: EdgeIndexMap edge_index_map(g); marci@1017: PrimalMap lp_flow(lp, edge_index_map); marci@1017: for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { marci@1017: ColIt col_it=lp.addCol(); marci@1017: edge_index_map.set(e, col_it); marci@1017: if (lcapacity[e]==capacity[e]) marci@1081: lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]); marci@1017: else marci@1081: lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]); marci@1017: lp.setObjCoef(col_it, cost[e]); marci@1017: } marci@1074: LPSolver::ColIt col_it; marci@1074: for (lp.col_iter_map.first(col_it, lp.VALID_CLASS); marci@1074: lp.col_iter_map.valid(col_it); marci@1074: lp.col_iter_map.next(col_it)) { marci@1081: // std::cout << "ize " << lp.col_iter_map[col_it] << std::endl; marci@1074: } marci@1017: for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { marci@1017: typename Graph::template EdgeMap coeffs(g, 0); marci@1017: for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) marci@1017: coeffs.set(e, coeffs[e]+1); marci@1017: for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) marci@1017: coeffs.set(e, coeffs[e]-1); marci@1017: RowIt row_it=lp.addRow(); marci@1017: typename std::vector< std::pair > row; marci@1017: //std::cout << "node:" <