1 // -*- c++ -*- |
|
2 #ifndef LEMON_MIN_COST_GEN_FLOW_H |
|
3 #define LEMON_MIN_COST_GEN_FLOW_H |
|
4 #include <iostream> |
|
5 //#include <fstream> |
|
6 |
|
7 #include <lemon/smart_graph.h> |
|
8 #include <lemon/list_graph.h> |
|
9 //#include <lemon/dimacs.h> |
|
10 //#include <lemon/time_measure.h> |
|
11 //#include <graph_wrapper.h> |
|
12 #include <lemon/preflow.h> |
|
13 #include <lemon/min_cost_flow.h> |
|
14 //#include <augmenting_flow.h> |
|
15 //#include <preflow_res.h> |
|
16 #include <work/marci/merge_node_graph_wrapper.h> |
|
17 #include <work/marci/lp/lp_solver_wrapper_3.h> |
|
18 |
|
19 namespace lemon { |
|
20 |
|
21 template<typename Edge, typename EdgeIndexMap> |
|
22 class PrimalMap { |
|
23 protected: |
|
24 LPGLPK* lp; |
|
25 EdgeIndexMap* edge_index_map; |
|
26 public: |
|
27 PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : |
|
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 |
|
34 // excess: rho-delta egyelore csak =0-ra. |
|
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) { } |
|
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 |
|
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 } |
|
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); |
|
99 // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl; |
|
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]; |
|
107 // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl; |
|
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 } |
|
124 return (preflow.flowValue()>=expected); |
|
125 } |
|
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 } |
|
211 return (min_cost_flow.flowValue()>=expected); |
|
212 } |
|
213 void runByLP() { |
|
214 typedef LPGLPK LPSolver; |
|
215 LPSolver lp; |
|
216 lp.setMinimize(); |
|
217 typedef LPSolver::ColIt ColIt; |
|
218 typedef LPSolver::RowIt RowIt; |
|
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]) |
|
226 lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]); |
|
227 else |
|
228 lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]); |
|
229 lp.setObjCoef(col_it, cost[e]); |
|
230 } |
|
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)) { |
|
235 // std::cout << "ize " << lp.col_iter_map[col_it] << std::endl; |
|
236 } |
|
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; |
|
253 //std::cout << " " << g.id(n) << " " << row.size() << std::endl; |
|
254 lp.setRowCoeffs(row_it, row.begin(), row.end()); |
|
255 lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0); |
|
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 |
|