[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> |
---|
[1028] | 13 | #include <lemon/min_cost_flow.h> |
---|
[1017] | 14 | //#include <augmenting_flow.h> |
---|
| 15 | //#include <preflow_res.h> |
---|
[1028] | 16 | #include <work/marci/merge_node_graph_wrapper.h> |
---|
[1074] | 17 | #include <work/marci/lp/lp_solver_wrapper_3.h> |
---|
[1017] | 18 | |
---|
| 19 | namespace lemon { |
---|
| 20 | |
---|
| 21 | template<typename Edge, typename EdgeIndexMap> |
---|
| 22 | class PrimalMap { |
---|
| 23 | protected: |
---|
[1081] | 24 | LPGLPK* lp; |
---|
[1017] | 25 | EdgeIndexMap* edge_index_map; |
---|
| 26 | public: |
---|
[1081] | 27 | PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : |
---|
[1017] | 28 | lp(&_lp), edge_index_map(&_edge_index_map) { } |
---|
| 29 | double operator[](Edge e) const { |
---|
| 30 | return lp->getPrimal((*edge_index_map)[e]); |
---|
| 31 | } |
---|
| 32 | }; |
---|
| 33 | |
---|
[1028] | 34 | // excess: rho-delta egyelore csak =0-ra. |
---|
[1017] | 35 | template <typename Graph, typename Num, |
---|
| 36 | typename Excess=typename Graph::template NodeMap<Num>, |
---|
| 37 | typename LCapMap=typename Graph::template EdgeMap<Num>, |
---|
| 38 | typename CapMap=typename Graph::template EdgeMap<Num>, |
---|
| 39 | typename FlowMap=typename Graph::template EdgeMap<Num>, |
---|
| 40 | typename CostMap=typename Graph::template EdgeMap<Num> > |
---|
| 41 | class MinCostGenFlow { |
---|
| 42 | protected: |
---|
| 43 | const Graph& g; |
---|
| 44 | const Excess& excess; |
---|
| 45 | const LCapMap& lcapacity; |
---|
| 46 | const CapMap& capacity; |
---|
| 47 | FlowMap& flow; |
---|
| 48 | const CostMap& cost; |
---|
| 49 | public: |
---|
| 50 | MinCostGenFlow(const Graph& _g, const Excess& _excess, |
---|
| 51 | const LCapMap& _lcapacity, const CapMap& _capacity, |
---|
| 52 | FlowMap& _flow, |
---|
| 53 | const CostMap& _cost) : |
---|
| 54 | g(_g), excess(_excess), lcapacity(_lcapacity), |
---|
| 55 | capacity(_capacity), flow(_flow), cost(_cost) { } |
---|
[1025] | 56 | bool feasible() { |
---|
| 57 | // std::cout << "making new vertices..." << std::endl; |
---|
| 58 | typedef ListGraph Graph2; |
---|
| 59 | Graph2 g2; |
---|
| 60 | typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW; |
---|
| 61 | // std::cout << "merging..." << std::endl; |
---|
| 62 | GW gw(g, g2); |
---|
| 63 | typename GW::Node s(INVALID, g2.addNode(), true); |
---|
| 64 | typename GW::Node t(INVALID, g2.addNode(), true); |
---|
| 65 | typedef SmartGraph Graph3; |
---|
| 66 | // std::cout << "making extender graph..." << std::endl; |
---|
| 67 | typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW; |
---|
| 68 | // { |
---|
| 69 | // checkConcept<StaticGraph, GWW>(); |
---|
| 70 | // } |
---|
| 71 | GWW gww(gw); |
---|
| 72 | typedef AugmentingGraphWrapper<GW, GWW> GWWW; |
---|
| 73 | GWWW gwww(gw, gww); |
---|
| 74 | |
---|
| 75 | // std::cout << "making new edges..." << std::endl; |
---|
| 76 | typename GWWW::template EdgeMap<Num> translated_cap(gwww); |
---|
| 77 | |
---|
[1028] | 78 | for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) { |
---|
| 79 | translated_cap.set(typename GWWW::Edge(e,INVALID,false), |
---|
| 80 | capacity[e]-lcapacity[e]); |
---|
| 81 | // cout << "t_cap " << gw.id(e) << " " |
---|
| 82 | // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl; |
---|
| 83 | } |
---|
[1025] | 84 | |
---|
| 85 | Num expected=0; |
---|
| 86 | |
---|
| 87 | // std::cout << "making new edges 2..." << std::endl; |
---|
| 88 | for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
---|
| 89 | Num a=0; |
---|
| 90 | for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 91 | a+=lcapacity[e]; |
---|
| 92 | for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 93 | a-=lcapacity[e]; |
---|
| 94 | if (excess[n]>a) { |
---|
| 95 | typename GWW::Edge e= |
---|
| 96 | gww.addEdge(typename GW::Node(n,INVALID,false), t); |
---|
| 97 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 98 | excess[n]-a); |
---|
[1028] | 99 | // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl; |
---|
[1025] | 100 | } |
---|
| 101 | if (excess[n]<a) { |
---|
| 102 | typename GWW::Edge e= |
---|
| 103 | gww.addEdge(s, typename GW::Node(n,INVALID,false)); |
---|
| 104 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 105 | a-excess[n]); |
---|
| 106 | expected+=a-excess[n]; |
---|
[1028] | 107 | // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl; |
---|
[1025] | 108 | } |
---|
| 109 | } |
---|
| 110 | |
---|
| 111 | // std::cout << "preflow..." << std::endl; |
---|
| 112 | typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0); |
---|
| 113 | Preflow<GWWW, Num> preflow(gwww, s, t, |
---|
| 114 | translated_cap, translated_flow); |
---|
| 115 | preflow.run(); |
---|
| 116 | // std::cout << "fv: " << preflow.flowValue() << std::endl; |
---|
| 117 | // std::cout << "expected: " << expected << std::endl; |
---|
| 118 | |
---|
| 119 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 120 | typename GW::Edge ew(e, INVALID, false); |
---|
| 121 | typename GWWW::Edge ewww(ew, INVALID, false); |
---|
| 122 | flow.set(e, translated_flow[ewww]+lcapacity[e]); |
---|
| 123 | } |
---|
[1033] | 124 | return (preflow.flowValue()>=expected); |
---|
[1025] | 125 | } |
---|
[1028] | 126 | // for nonnegative costs |
---|
| 127 | bool run() { |
---|
| 128 | // std::cout << "making new vertices..." << std::endl; |
---|
| 129 | typedef ListGraph Graph2; |
---|
| 130 | Graph2 g2; |
---|
| 131 | typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW; |
---|
| 132 | // std::cout << "merging..." << std::endl; |
---|
| 133 | GW gw(g, g2); |
---|
| 134 | typename GW::Node s(INVALID, g2.addNode(), true); |
---|
| 135 | typename GW::Node t(INVALID, g2.addNode(), true); |
---|
| 136 | typedef SmartGraph Graph3; |
---|
| 137 | // std::cout << "making extender graph..." << std::endl; |
---|
| 138 | typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW; |
---|
| 139 | // { |
---|
| 140 | // checkConcept<StaticGraph, GWW>(); |
---|
| 141 | // } |
---|
| 142 | GWW gww(gw); |
---|
| 143 | typedef AugmentingGraphWrapper<GW, GWW> GWWW; |
---|
| 144 | GWWW gwww(gw, gww); |
---|
| 145 | |
---|
| 146 | // std::cout << "making new edges..." << std::endl; |
---|
| 147 | typename GWWW::template EdgeMap<Num> translated_cap(gwww); |
---|
| 148 | |
---|
| 149 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 150 | typename GW::Edge ew(e, INVALID, false); |
---|
| 151 | typename GWWW::Edge ewww(ew, INVALID, false); |
---|
| 152 | translated_cap.set(ewww, capacity[e]-lcapacity[e]); |
---|
| 153 | // cout << "t_cap " << g.id(e) << " " |
---|
| 154 | // << translated_cap[ewww] << endl; |
---|
| 155 | } |
---|
| 156 | |
---|
| 157 | Num expected=0; |
---|
| 158 | |
---|
| 159 | // std::cout << "making new edges 2..." << std::endl; |
---|
| 160 | for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
---|
| 161 | // std::cout << "node: " << g.id(n) << std::endl; |
---|
| 162 | Num a=0; |
---|
| 163 | for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) { |
---|
| 164 | a+=lcapacity[e]; |
---|
| 165 | // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl; |
---|
| 166 | } |
---|
| 167 | for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) { |
---|
| 168 | a-=lcapacity[e]; |
---|
| 169 | // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl; |
---|
| 170 | } |
---|
| 171 | // std::cout << "excess " << g.id(n) << ": " << a << std::endl; |
---|
| 172 | if (0>a) { |
---|
| 173 | typename GWW::Edge e= |
---|
| 174 | gww.addEdge(typename GW::Node(n,INVALID,false), t); |
---|
| 175 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 176 | -a); |
---|
| 177 | // std::cout << g.id(n) << "->t " << -a << std::endl; |
---|
| 178 | } |
---|
| 179 | if (0<a) { |
---|
| 180 | typename GWW::Edge e= |
---|
| 181 | gww.addEdge(s, typename GW::Node(n,INVALID,false)); |
---|
| 182 | translated_cap.set(typename GWWW::Edge(INVALID, e, true), |
---|
| 183 | a); |
---|
| 184 | expected+=a; |
---|
| 185 | // std::cout << "s->" << g.id(n) << " "<< a <<std:: endl; |
---|
| 186 | } |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | // std::cout << "preflow..." << std::endl; |
---|
| 190 | typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0); |
---|
| 191 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 192 | translated_cost.set(typename GWWW::Edge( |
---|
| 193 | typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]); |
---|
| 194 | } |
---|
| 195 | // typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0); |
---|
| 196 | MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>, |
---|
| 197 | typename GWWW::template EdgeMap<Num> > |
---|
| 198 | min_cost_flow(gwww, translated_cost, translated_cap, |
---|
| 199 | s, t); |
---|
| 200 | while (min_cost_flow.augment()) { } |
---|
| 201 | std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; |
---|
| 202 | std::cout << "expected: " << expected << std::endl; |
---|
| 203 | |
---|
| 204 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 205 | typename GW::Edge ew(e, INVALID, false); |
---|
| 206 | typename GWWW::Edge ewww(ew, INVALID, false); |
---|
| 207 | // std::cout << g.id(e) << " " << flow[e] << std::endl; |
---|
| 208 | flow.set(e, lcapacity[e]+ |
---|
| 209 | min_cost_flow.getFlow()[ewww]); |
---|
| 210 | } |
---|
[1033] | 211 | return (min_cost_flow.flowValue()>=expected); |
---|
[1028] | 212 | } |
---|
| 213 | void runByLP() { |
---|
[1081] | 214 | typedef LPGLPK LPSolver; |
---|
[1031] | 215 | LPSolver lp; |
---|
[1017] | 216 | lp.setMinimize(); |
---|
[1031] | 217 | typedef LPSolver::ColIt ColIt; |
---|
| 218 | typedef LPSolver::RowIt RowIt; |
---|
[1017] | 219 | typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap; |
---|
| 220 | EdgeIndexMap edge_index_map(g); |
---|
| 221 | PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map); |
---|
| 222 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 223 | ColIt col_it=lp.addCol(); |
---|
| 224 | edge_index_map.set(e, col_it); |
---|
| 225 | if (lcapacity[e]==capacity[e]) |
---|
[1081] | 226 | lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]); |
---|
[1017] | 227 | else |
---|
[1081] | 228 | lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]); |
---|
[1017] | 229 | lp.setObjCoef(col_it, cost[e]); |
---|
| 230 | } |
---|
[1074] | 231 | LPSolver::ColIt col_it; |
---|
| 232 | for (lp.col_iter_map.first(col_it, lp.VALID_CLASS); |
---|
| 233 | lp.col_iter_map.valid(col_it); |
---|
| 234 | lp.col_iter_map.next(col_it)) { |
---|
[1081] | 235 | // std::cout << "ize " << lp.col_iter_map[col_it] << std::endl; |
---|
[1074] | 236 | } |
---|
[1017] | 237 | for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
---|
| 238 | typename Graph::template EdgeMap<Num> coeffs(g, 0); |
---|
| 239 | for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 240 | coeffs.set(e, coeffs[e]+1); |
---|
| 241 | for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
---|
| 242 | coeffs.set(e, coeffs[e]-1); |
---|
| 243 | RowIt row_it=lp.addRow(); |
---|
| 244 | typename std::vector< std::pair<ColIt, double> > row; |
---|
| 245 | //std::cout << "node:" <<g.id(n)<<std::endl; |
---|
| 246 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
---|
| 247 | if (coeffs[e]!=0) { |
---|
| 248 | //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e]; |
---|
| 249 | row.push_back(std::make_pair(edge_index_map[e], coeffs[e])); |
---|
| 250 | } |
---|
| 251 | } |
---|
| 252 | //std::cout << std::endl; |
---|
[1081] | 253 | //std::cout << " " << g.id(n) << " " << row.size() << std::endl; |
---|
[1017] | 254 | lp.setRowCoeffs(row_it, row.begin(), row.end()); |
---|
[1081] | 255 | lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0); |
---|
[1017] | 256 | } |
---|
| 257 | lp.solveSimplex(); |
---|
| 258 | //std::cout << lp.colNum() << std::endl; |
---|
| 259 | //std::cout << lp.rowNum() << std::endl; |
---|
| 260 | //std::cout << "flow value: "<< lp.getObjVal() << std::endl; |
---|
| 261 | for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) |
---|
| 262 | flow.set(e, lp_flow[e]); |
---|
| 263 | } |
---|
| 264 | }; |
---|
| 265 | |
---|
| 266 | } // namespace lemon |
---|
| 267 | |
---|
| 268 | #endif //LEMON_MIN_COST_GEN_FLOW_H |
---|