[1017] | 1 | // -*- c++ -*- |
---|
| 2 | #ifndef LEMON_MIN_COST_GEN_FLOW_H |
---|
| 3 | #define LEMON_MIN_COST_GEN_FLOW_H |
---|
[1025] | 4 | #include <iostream> |
---|
[1017] | 5 | //#include <fstream> |
---|
| 6 | |
---|
[1025] | 7 | #include <lemon/smart_graph.h> |
---|
| 8 | #include <lemon/list_graph.h> |
---|
[1017] | 9 | //#include <lemon/dimacs.h> |
---|
| 10 | //#include <lemon/time_measure.h> |
---|
| 11 | //#include <graph_wrapper.h> |
---|
[1025] | 12 | #include <lemon/preflow.h> |
---|
[1017] | 13 | //#include <augmenting_flow.h> |
---|
| 14 | //#include <preflow_res.h> |
---|
[1025] | 15 | #include <../merge_node_graph_wrapper.h> |
---|
[1017] | 16 | #include <lemon/../work/marci/lp/lp_solver_wrapper.h> |
---|
| 17 | |
---|
| 18 | namespace lemon { |
---|
| 19 | |
---|
| 20 | template<typename Edge, typename EdgeIndexMap> |
---|
| 21 | class PrimalMap { |
---|
| 22 | protected: |
---|
| 23 | LPSolverWrapper* lp; |
---|
| 24 | EdgeIndexMap* edge_index_map; |
---|
| 25 | public: |
---|
| 26 | PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : |
---|
| 27 | lp(&_lp), edge_index_map(&_edge_index_map) { } |
---|
| 28 | double operator[](Edge e) const { |
---|
| 29 | return lp->getPrimal((*edge_index_map)[e]); |
---|
| 30 | } |
---|
| 31 | }; |
---|
| 32 | |
---|
| 33 | // excess: rho-delta |
---|
| 34 | template <typename Graph, typename Num, |
---|
| 35 | typename Excess=typename Graph::template NodeMap<Num>, |
---|
| 36 | typename LCapMap=typename Graph::template EdgeMap<Num>, |
---|
| 37 | typename CapMap=typename Graph::template EdgeMap<Num>, |
---|
| 38 | typename FlowMap=typename Graph::template EdgeMap<Num>, |
---|
| 39 | typename CostMap=typename Graph::template EdgeMap<Num> > |
---|
| 40 | class MinCostGenFlow { |
---|
| 41 | protected: |
---|
| 42 | const Graph& g; |
---|
| 43 | const Excess& excess; |
---|
| 44 | const LCapMap& lcapacity; |
---|
| 45 | const CapMap& capacity; |
---|
| 46 | FlowMap& flow; |
---|
| 47 | const CostMap& cost; |
---|
| 48 | public: |
---|
| 49 | MinCostGenFlow(const Graph& _g, const Excess& _excess, |
---|
| 50 | const LCapMap& _lcapacity, const CapMap& _capacity, |
---|
| 51 | FlowMap& _flow, |
---|
| 52 | const CostMap& _cost) : |
---|
| 53 | g(_g), excess(_excess), lcapacity(_lcapacity), |
---|
| 54 | capacity(_capacity), flow(_flow), cost(_cost) { } |
---|
[1025] | 55 | bool feasible() { |
---|
| 56 | // std::cout << "making new vertices..." << std::endl; |
---|
| 57 | typedef ListGraph Graph2; |
---|
| 58 | Graph2 g2; |
---|
| 59 | typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW; |
---|
| 60 | // std::cout << "merging..." << std::endl; |
---|
| 61 | GW gw(g, g2); |
---|
| 62 | typename GW::Node s(INVALID, g2.addNode(), true); |
---|
| 63 | typename GW::Node t(INVALID, g2.addNode(), true); |
---|
| 64 | typedef SmartGraph Graph3; |
---|
| 65 | // std::cout << "making extender graph..." << std::endl; |
---|
| 66 | typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW; |
---|
| 67 | // { |
---|
| 68 | // checkConcept<StaticGraph, GWW>(); |
---|
| 69 | // } |
---|
| 70 | GWW gww(gw); |
---|
| 71 | typedef AugmentingGraphWrapper<GW, GWW> GWWW; |
---|
| 72 | GWWW gwww(gw, gww); |
---|
| 73 | |
---|
| 74 | // std::cout << "making new edges..." << std::endl; |
---|
| 75 | typename GWWW::template EdgeMap<Num> translated_cap(gwww); |
---|
| 76 | |
---|
| 77 | for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) |
---|
| 78 | translated_cap.set(typename GWWW::Edge(e,INVALID,false), |
---|
| 79 | capacity[e]-lcapacity[e]); |
---|
| 80 | |
---|
| 81 | Num expected=0; |
---|
| 82 | |
---|
| 83 | // std::cout << "making new edges 2..." << std::endl; |
---|
| 84 | for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
---|
| 85 | Num a=0; |
---|
| 86 | for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 87 | a+=lcapacity[e]; |
---|
| 88 | for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 89 | a-=lcapacity[e]; |
---|
| 90 | if (excess[n]>a) { |
---|
| 91 | typename GWW::Edge e= |
---|
| 92 | gww.addEdge(typename GW::Node(n,INVALID,false), t); |
---|
| 93 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 94 | excess[n]-a); |
---|
| 95 | } |
---|
| 96 | if (excess[n]<a) { |
---|
| 97 | typename GWW::Edge e= |
---|
| 98 | gww.addEdge(s, typename GW::Node(n,INVALID,false)); |
---|
| 99 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 100 | a-excess[n]); |
---|
| 101 | expected+=a-excess[n]; |
---|
| 102 | } |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | // std::cout << "preflow..." << std::endl; |
---|
| 106 | typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0); |
---|
| 107 | Preflow<GWWW, Num> preflow(gwww, s, t, |
---|
| 108 | translated_cap, translated_flow); |
---|
| 109 | preflow.run(); |
---|
| 110 | // std::cout << "fv: " << preflow.flowValue() << std::endl; |
---|
| 111 | // std::cout << "expected: " << expected << std::endl; |
---|
| 112 | |
---|
| 113 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 114 | typename GW::Edge ew(e, INVALID, false); |
---|
| 115 | typename GWWW::Edge ewww(ew, INVALID, false); |
---|
| 116 | flow.set(e, translated_flow[ewww]+lcapacity[e]); |
---|
| 117 | } |
---|
| 118 | return (expected>=preflow.flowValue()); |
---|
| 119 | } |
---|
[1017] | 120 | void run() { |
---|
| 121 | LPSolverWrapper lp; |
---|
| 122 | lp.setMinimize(); |
---|
| 123 | typedef LPSolverWrapper::ColIt ColIt; |
---|
| 124 | typedef LPSolverWrapper::RowIt RowIt; |
---|
| 125 | typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap; |
---|
| 126 | EdgeIndexMap edge_index_map(g); |
---|
| 127 | PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map); |
---|
| 128 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 129 | ColIt col_it=lp.addCol(); |
---|
| 130 | edge_index_map.set(e, col_it); |
---|
| 131 | if (lcapacity[e]==capacity[e]) |
---|
| 132 | lp.setColBounds(col_it, LPX_FX, lcapacity[e], capacity[e]); |
---|
| 133 | else |
---|
| 134 | lp.setColBounds(col_it, LPX_DB, lcapacity[e], capacity[e]); |
---|
| 135 | lp.setObjCoef(col_it, cost[e]); |
---|
| 136 | } |
---|
| 137 | for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
---|
| 138 | typename Graph::template EdgeMap<Num> coeffs(g, 0); |
---|
| 139 | for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 140 | coeffs.set(e, coeffs[e]+1); |
---|
| 141 | for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 142 | coeffs.set(e, coeffs[e]-1); |
---|
| 143 | RowIt row_it=lp.addRow(); |
---|
| 144 | typename std::vector< std::pair<ColIt, double> > row; |
---|
| 145 | //std::cout << "node:" <<g.id(n)<<std::endl; |
---|
| 146 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 147 | if (coeffs[e]!=0) { |
---|
| 148 | //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e]; |
---|
| 149 | row.push_back(std::make_pair(edge_index_map[e], coeffs[e])); |
---|
| 150 | } |
---|
| 151 | } |
---|
| 152 | //std::cout << std::endl; |
---|
| 153 | lp.setRowCoeffs(row_it, row.begin(), row.end()); |
---|
| 154 | lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0); |
---|
| 155 | } |
---|
| 156 | lp.solveSimplex(); |
---|
| 157 | //std::cout << lp.colNum() << std::endl; |
---|
| 158 | //std::cout << lp.rowNum() << std::endl; |
---|
| 159 | //std::cout << "flow value: "<< lp.getObjVal() << std::endl; |
---|
| 160 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) |
---|
| 161 | flow.set(e, lp_flow[e]); |
---|
| 162 | } |
---|
| 163 | }; |
---|
| 164 | |
---|
| 165 | } // namespace lemon |
---|
| 166 | |
---|
| 167 | #endif //LEMON_MIN_COST_GEN_FLOW_H |
---|