10 #include <lemon/preflow.h> |
10 #include <lemon/preflow.h> |
11 #include <augmenting_flow.h> |
11 #include <augmenting_flow.h> |
12 //#include <preflow_res.h> |
12 //#include <preflow_res.h> |
13 #include <lp_solver_wrapper.h> |
13 #include <lp_solver_wrapper.h> |
14 |
14 |
15 using namespace lemon; |
|
16 |
|
17 // Use a DIMACS max flow file as stdin. |
15 // Use a DIMACS max flow file as stdin. |
18 // max_flow_demo < dimacs_max_flow_file |
16 // max_flow_demo < dimacs_max_flow_file |
19 |
17 |
20 template<typename Edge, typename EdgeIndexMap> |
18 using namespace lemon; |
21 class PrimalMap { |
19 |
22 protected: |
20 namespace lemon { |
23 LPSolverWrapper* lp; |
21 |
24 EdgeIndexMap* edge_index_map; |
22 template<typename Edge, typename EdgeIndexMap> |
25 public: |
23 class PrimalMap { |
26 PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : |
24 protected: |
27 lp(&_lp), edge_index_map(&_edge_index_map) { } |
25 LPSolverWrapper* lp; |
28 double operator[](Edge e) const { |
26 EdgeIndexMap* edge_index_map; |
29 return lp->getPrimal((*edge_index_map)[e]); |
27 public: |
30 } |
28 PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : |
31 }; |
29 lp(&_lp), edge_index_map(&_edge_index_map) { } |
|
30 double operator[](Edge e) const { |
|
31 return lp->getPrimal((*edge_index_map)[e]); |
|
32 } |
|
33 }; |
|
34 |
|
35 // excess: rho-delta |
|
36 template <typename Graph, typename Num, |
|
37 typename Excess=typename Graph::template NodeMap<Num>, |
|
38 typename LCapMap=typename Graph::template EdgeMap<Num>, |
|
39 typename CapMap=typename Graph::template EdgeMap<Num>, |
|
40 typename FlowMap=typename Graph::template EdgeMap<Num>, |
|
41 typename CostMap=typename Graph::template EdgeMap<Num> > |
|
42 class MinCostGenFlow { |
|
43 protected: |
|
44 const Graph& g; |
|
45 const Excess& excess; |
|
46 const LCapMap& lcapacity; |
|
47 const CapMap& capacity; |
|
48 FlowMap& flow; |
|
49 const CostMap& cost; |
|
50 public: |
|
51 MinCostGenFlow(const Graph& _g, const Excess& _excess, |
|
52 const LCapMap& _lcapacity, const CapMap& _capacity, |
|
53 FlowMap& _flow, |
|
54 const CostMap& _cost) : |
|
55 g(_g), excess(_excess), lcapacity(_lcapacity), |
|
56 capacity(_capacity), flow(_flow), cost(_cost) { } |
|
57 void run() { |
|
58 LPSolverWrapper lp; |
|
59 lp.setMinimize(); |
|
60 typedef LPSolverWrapper::ColIt ColIt; |
|
61 typedef LPSolverWrapper::RowIt RowIt; |
|
62 typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap; |
|
63 EdgeIndexMap edge_index_map(g); |
|
64 PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map); |
|
65 for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
|
66 ColIt col_it=lp.addCol(); |
|
67 edge_index_map.set(e, col_it); |
|
68 if (lcapacity[e]==capacity[e]) |
|
69 lp.setColBounds(col_it, LPX_FX, lcapacity[e], capacity[e]); |
|
70 else |
|
71 lp.setColBounds(col_it, LPX_DB, lcapacity[e], capacity[e]); |
|
72 lp.setObjCoef(col_it, cost[e]); |
|
73 } |
|
74 for (typename Graph::NodeIt n(g); n!=INVALID; ++n) { |
|
75 typename Graph::template EdgeMap<Num> coeffs(g, 0); |
|
76 for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
|
77 coeffs.set(e, coeffs[e]+1); |
|
78 for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
|
79 coeffs.set(e, coeffs[e]-1); |
|
80 RowIt row_it=lp.addRow(); |
|
81 typename std::vector< std::pair<ColIt, double> > row; |
|
82 //std::cout << "node:" <<g.id(n)<<std::endl; |
|
83 for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) { |
|
84 if (coeffs[e]!=0) { |
|
85 //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e]; |
|
86 row.push_back(std::make_pair(edge_index_map[e], coeffs[e])); |
|
87 } |
|
88 } |
|
89 //std::cout << std::endl; |
|
90 lp.setRowCoeffs(row_it, row.begin(), row.end()); |
|
91 lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0); |
|
92 } |
|
93 lp.solveSimplex(); |
|
94 //std::cout << lp.colNum() << std::endl; |
|
95 //std::cout << lp.rowNum() << std::endl; |
|
96 //std::cout << "flow value: "<< lp.getObjVal() << std::endl; |
|
97 for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) |
|
98 flow.set(e, lp_flow[e]); |
|
99 } |
|
100 }; |
|
101 |
|
102 } |
32 |
103 |
33 int main(int, char **) { |
104 int main(int, char **) { |
34 |
105 |
35 typedef ListGraph MutableGraph; |
106 typedef ListGraph MutableGraph; |
36 typedef SmartGraph Graph; |
107 typedef SmartGraph Graph; |
37 typedef Graph::Node Node; |
108 typedef Graph::Node Node; |
|
109 typedef Graph::Edge Edge; |
38 typedef Graph::EdgeIt EdgeIt; |
110 typedef Graph::EdgeIt EdgeIt; |
39 |
111 |
40 Graph g; |
112 Graph g; |
41 |
113 |
42 Node s, t; |
114 Node s, t; |
93 // std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl; |
165 // std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl; |
94 // } |
166 // } |
95 |
167 |
96 { |
168 { |
97 std::cout << "physical blocking flow augmentation ..." << std::endl; |
169 std::cout << "physical blocking flow augmentation ..." << std::endl; |
98 for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0); |
170 for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0); |
99 ts.reset(); |
171 ts.reset(); |
100 int i=0; |
172 int i=0; |
101 while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; } |
173 while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; } |
102 std::cout << "elapsed time: " << ts << std::endl; |
174 std::cout << "elapsed time: " << ts << std::endl; |
103 std::cout << "number of augmentation phases: " << i << std::endl; |
175 std::cout << "number of augmentation phases: " << i << std::endl; |
104 std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; |
176 std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; |
105 |
177 |
106 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { |
178 for (EdgeIt e(g); e!=INVALID; ++e) { |
107 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) |
179 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) |
108 std::cout << "Slackness does not hold!" << std::endl; |
180 std::cout << "Slackness does not hold!" << std::endl; |
109 if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) |
181 if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) |
110 std::cout << "Slackness does not hold!" << std::endl; |
182 std::cout << "Slackness does not hold!" << std::endl; |
111 } |
183 } |
122 // std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl; |
194 // std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl; |
123 // } |
195 // } |
124 |
196 |
125 { |
197 { |
126 std::cout << "on-the-fly blocking flow augmentation ..." << std::endl; |
198 std::cout << "on-the-fly blocking flow augmentation ..." << std::endl; |
127 for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0); |
199 for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0); |
128 ts.reset(); |
200 ts.reset(); |
129 int i=0; |
201 int i=0; |
130 while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; } |
202 while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; } |
131 std::cout << "elapsed time: " << ts << std::endl; |
203 std::cout << "elapsed time: " << ts << std::endl; |
132 std::cout << "number of augmentation phases: " << i << std::endl; |
204 std::cout << "number of augmentation phases: " << i << std::endl; |
133 std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; |
205 std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; |
134 |
206 |
135 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { |
207 for (EdgeIt e(g); e!=INVALID; ++e) { |
136 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) |
208 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) |
137 std::cout << "Slackness does not hold!" << std::endl; |
209 std::cout << "Slackness does not hold!" << std::endl; |
138 if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) |
210 if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) |
139 std::cout << "Slackness does not hold!" << std::endl; |
211 std::cout << "Slackness does not hold!" << std::endl; |
140 } |
212 } |
175 // std::cout << "Slackness does not hold!" << std::endl; |
247 // std::cout << "Slackness does not hold!" << std::endl; |
176 // } |
248 // } |
177 // } |
249 // } |
178 |
250 |
179 ts.reset(); |
251 ts.reset(); |
180 LPSolverWrapper lp; |
252 |
181 lp.setMaximize(); |
253 Edge e=g.addEdge(t, s); |
182 typedef LPSolverWrapper::ColIt ColIt; |
254 Graph::EdgeMap<int> cost(g, 0); |
183 typedef LPSolverWrapper::RowIt RowIt; |
255 cost.set(e, -1); |
184 typedef Graph::EdgeMap<ColIt> EdgeIndexMap; |
256 cap.set(e, 10000); |
185 EdgeIndexMap edge_index_map(g); |
257 typedef ConstMap<Node, int> Excess; |
186 PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map); |
258 Excess excess(0); |
187 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { |
259 typedef ConstMap<Edge, int> LCap; |
188 ColIt col_it=lp.addCol(); |
260 LCap lcap(0); |
189 edge_index_map.set(e, col_it); |
261 |
190 lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]); |
262 MinCostGenFlow<Graph, int, Excess, LCap> |
191 } |
263 min_cost(g, excess, lcap, cap, flow, cost); |
192 for (Graph::NodeIt n(g); n!=INVALID; ++n) { |
264 min_cost.run(); |
193 if (n!=s) { |
265 |
194 //hurokelek miatt |
|
195 Graph::EdgeMap<int> coeffs(g, 0); |
|
196 for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e) |
|
197 coeffs.set(e, coeffs[e]+1); |
|
198 for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) |
|
199 coeffs.set(e, coeffs[e]-1); |
|
200 if (n==t) { |
|
201 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { |
|
202 if (coeffs[e]!=0) |
|
203 lp.setObjCoef(edge_index_map[e], coeffs[e]); |
|
204 } |
|
205 } else { |
|
206 RowIt row_it=lp.addRow(); |
|
207 std::vector< std::pair<ColIt, double> > row; |
|
208 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { |
|
209 if (coeffs[e]!=0) |
|
210 row.push_back(std::make_pair(edge_index_map[e], coeffs[e])); |
|
211 } |
|
212 lp.setRowCoeffs(row_it, row.begin(), row.end()); |
|
213 lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0); |
|
214 } |
|
215 } |
|
216 } |
|
217 lp.solveSimplex(); |
|
218 std::cout << "flow value: "<< lp.getObjVal() << std::endl; |
|
219 std::cout << "elapsed time: " << ts << std::endl; |
266 std::cout << "elapsed time: " << ts << std::endl; |
220 |
267 std::cout << "flow value: "<< flow[e] << std::endl; |
221 return 0; |
|
222 } |
268 } |