COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/lp/max_flow_by_lp.cc @ 1015:e3bb0e118bb4

Last change on this file since 1015:e3bb0e118bb4 was 1015:e3bb0e118bb4, checked in by marci, 17 years ago

RoadMap? to more general flow algs.

File size: 9.3 KB
Line 
1// -*- c++ -*-
2#include <iostream>
3#include <fstream>
4
5#include <lemon/smart_graph.h>
6#include <lemon/list_graph.h>
7#include <lemon/dimacs.h>
8#include <lemon/time_measure.h>
9//#include <graph_wrapper.h>
10#include <lemon/preflow.h>
11#include <augmenting_flow.h>
12//#include <preflow_res.h>
13#include <lp_solver_wrapper.h>
14
15// Use a DIMACS max flow file as stdin.
16// max_flow_demo < dimacs_max_flow_file
17
18using namespace lemon;
19
20namespace lemon {
21
22  template<typename Edge, typename EdgeIndexMap>
23  class PrimalMap {
24  protected:
25    LPSolverWrapper* lp;
26    EdgeIndexMap* edge_index_map;
27  public:
28    PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) :
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}
103
104int main(int, char **) {
105
106  typedef ListGraph MutableGraph;
107  typedef SmartGraph Graph;
108  typedef Graph::Node Node;
109  typedef Graph::Edge Edge;
110  typedef Graph::EdgeIt EdgeIt;
111
112  Graph g;
113
114  Node s, t;
115  Graph::EdgeMap<int> cap(g);
116  //readDimacsMaxFlow(std::cin, g, s, t, cap);
117  readDimacs(std::cin, g, cap, s, t);
118  Timer ts;
119  Graph::EdgeMap<int> flow(g); //0 flow
120  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
121    max_flow_test(g, s, t, cap, flow);
122  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
123    augmenting_flow_test(g, s, t, cap, flow);
124 
125  Graph::NodeMap<bool> cut(g);
126
127  {
128    std::cout << "preflow ..." << std::endl;
129    ts.reset();
130    max_flow_test.run();
131    std::cout << "elapsed time: " << ts << std::endl;
132    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
133    max_flow_test.minCut(cut);
134
135    for (EdgeIt e(g); e!=INVALID; ++e) {
136      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
137        std::cout << "Slackness does not hold!" << std::endl;
138      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
139        std::cout << "Slackness does not hold!" << std::endl;
140    }
141  }
142
143//   {
144//     std::cout << "preflow ..." << std::endl;
145//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
146//     ts.reset();
147//     max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
148//     std::cout << "elapsed time: " << ts << std::endl;
149//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
150
151//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
152//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
153//      std::cout << "Slackness does not hold!" << std::endl;
154//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
155//      std::cout << "Slackness does not hold!" << std::endl;
156//     }
157//   }
158
159//   {
160//     std::cout << "wrapped preflow ..." << std::endl;
161//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
162//     ts.reset();
163//     pre_flow_res.run();
164//     std::cout << "elapsed time: " << ts << std::endl;
165//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
166//   }
167
168  {
169    std::cout << "physical blocking flow augmentation ..." << std::endl;
170    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
171    ts.reset();
172    int i=0;
173    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
174    std::cout << "elapsed time: " << ts << std::endl;
175    std::cout << "number of augmentation phases: " << i << std::endl;
176    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
177
178    for (EdgeIt e(g); e!=INVALID; ++e) {
179      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
180        std::cout << "Slackness does not hold!" << std::endl;
181      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
182        std::cout << "Slackness does not hold!" << std::endl;
183    }
184  }
185
186//   {
187//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
188//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
189//     ts.reset();
190//     int i=0;
191//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
192//     std::cout << "elapsed time: " << ts << std::endl;
193//     std::cout << "number of augmentation phases: " << i << std::endl;
194//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
195//   }
196
197  {
198    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
199    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
200    ts.reset();
201    int i=0;
202    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
203    std::cout << "elapsed time: " << ts << std::endl;
204    std::cout << "number of augmentation phases: " << i << std::endl;
205    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
206
207    for (EdgeIt e(g); e!=INVALID; ++e) {
208      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
209        std::cout << "Slackness does not hold!" << std::endl;
210      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
211        std::cout << "Slackness does not hold!" << std::endl;
212    }
213  }
214
215//   {
216//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
217//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
218//     ts.reset();
219//     int i=0;
220//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
221//     std::cout << "elapsed time: " << ts << std::endl;
222//     std::cout << "number of augmentation phases: " << i << std::endl;
223//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
224
225//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
226//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
227//      std::cout << "Slackness does not hold!" << std::endl;
228//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
229//      std::cout << "Slackness does not hold!" << std::endl;
230//     }
231//   }
232
233//   {
234//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
235//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
236//     ts.reset();
237//     int i=0;
238//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
239//     std::cout << "elapsed time: " << ts << std::endl;
240//     std::cout << "number of augmentation phases: " << i << std::endl;
241//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
242
243//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
244//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
245//      std::cout << "Slackness does not hold!" << std::endl;
246//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
247//      std::cout << "Slackness does not hold!" << std::endl;
248//     }
249//   }
250
251  ts.reset();
252
253  Edge e=g.addEdge(t, s);
254  Graph::EdgeMap<int> cost(g, 0);
255  cost.set(e, -1);
256  cap.set(e, 10000);
257  typedef ConstMap<Node, int> Excess;
258  Excess excess(0);
259  typedef ConstMap<Edge, int> LCap;
260  LCap lcap(0);
261
262  MinCostGenFlow<Graph, int, Excess, LCap>
263    min_cost(g, excess, lcap, cap, flow, cost);
264  min_cost.run();
265
266  std::cout << "elapsed time: " << ts << std::endl;
267  std::cout << "flow value: "<< flow[e] << std::endl;
268}
Note: See TracBrowser for help on using the repository browser.