COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/lp/max_flow_by_lp.cc @ 888:cc3590763f7f

Last change on this file since 888:cc3590763f7f was 764:615aca7091d2, checked in by marci, 17 years ago

An experimental LPSolverWrapper class which uses glpk. For a short
demo, max flow problems are solved with it. This demo does not
demonstrates, but the main aims of this class are row and column
generation capabilities, i.e. to be a core for easily
implementable branch-and-cut a column generetion algorithms.

File size: 7.8 KB
Line 
1// -*- c++ -*-
2#include <iostream>
3#include <fstream>
4
5#include <sage_graph.h>
6#include <hugo/smart_graph.h>
7#include <hugo/dimacs.h>
8#include <hugo/time_measure.h>
9//#include <graph_wrapper.h>
10#include <hugo/max_flow.h>
11#include <augmenting_flow.h>
12//#include <preflow_res.h>
13#include <for_each_macros.h>
14#include <lp_solver_wrapper.h>
15
16using namespace hugo;
17
18// Use a DIMACS max flow file as stdin.
19// max_flow_demo < dimacs_max_flow_file
20
21template<typename Edge, typename EdgeIndexMap>
22class PrimalMap {
23protected:
24  LPSolverWrapper* lp;
25  EdgeIndexMap* edge_index_map;
26public:
27  PrimalMap(LPSolverWrapper& _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
34int main(int, char **) {
35
36  typedef SageGraph MutableGraph;
37  //typedef SmartGraph Graph;
38  typedef SageGraph Graph;
39  typedef Graph::Node Node;
40  typedef Graph::EdgeIt EdgeIt;
41
42  Graph g;
43
44  Node s, t;
45  Graph::EdgeMap<int> cap(g);
46  //readDimacsMaxFlow(std::cin, g, s, t, cap);
47  readDimacs(std::cin, g, cap, s, t);
48  Timer ts;
49  Graph::EdgeMap<int> flow(g); //0 flow
50  MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
51    max_flow_test(g, s, t, cap, flow);
52  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
53    augmenting_flow_test(g, s, t, cap, flow);
54 
55  Graph::NodeMap<bool> cut(g);
56
57  {
58    std::cout << "preflow ..." << std::endl;
59    ts.reset();
60    max_flow_test.run();
61    std::cout << "elapsed time: " << ts << std::endl;
62    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
63    max_flow_test.actMinCut(cut);
64
65    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
66      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
67        std::cout << "Slackness does not hold!" << std::endl;
68      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
69        std::cout << "Slackness does not hold!" << std::endl;
70    }
71  }
72
73//   {
74//     std::cout << "preflow ..." << std::endl;
75//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
76//     ts.reset();
77//     max_flow_test.preflow(MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
78//     std::cout << "elapsed time: " << ts << std::endl;
79//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
80
81//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
82//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
83//      std::cout << "Slackness does not hold!" << std::endl;
84//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
85//      std::cout << "Slackness does not hold!" << std::endl;
86//     }
87//   }
88
89//   {
90//     std::cout << "wrapped preflow ..." << std::endl;
91//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
92//     ts.reset();
93//     pre_flow_res.run();
94//     std::cout << "elapsed time: " << ts << std::endl;
95//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
96//   }
97
98  {
99    std::cout << "physical blocking flow augmentation ..." << std::endl;
100    FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
101    ts.reset();
102    int i=0;
103    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
104    std::cout << "elapsed time: " << ts << std::endl;
105    std::cout << "number of augmentation phases: " << i << std::endl;
106    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
107
108    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
109      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
110        std::cout << "Slackness does not hold!" << std::endl;
111      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
112        std::cout << "Slackness does not hold!" << std::endl;
113    }
114  }
115
116//   {
117//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
118//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
119//     ts.reset();
120//     int i=0;
121//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
122//     std::cout << "elapsed time: " << ts << std::endl;
123//     std::cout << "number of augmentation phases: " << i << std::endl;
124//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
125//   }
126
127  {
128    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
129    FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
130    ts.reset();
131    int i=0;
132    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
133    std::cout << "elapsed time: " << ts << std::endl;
134    std::cout << "number of augmentation phases: " << i << std::endl;
135    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
136
137    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
138      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
139        std::cout << "Slackness does not hold!" << std::endl;
140      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
141        std::cout << "Slackness does not hold!" << std::endl;
142    }
143  }
144
145//   {
146//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
147//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
148//     ts.reset();
149//     int i=0;
150//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
151//     std::cout << "elapsed time: " << ts << std::endl;
152//     std::cout << "number of augmentation phases: " << i << std::endl;
153//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
154
155//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
156//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
157//      std::cout << "Slackness does not hold!" << std::endl;
158//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
159//      std::cout << "Slackness does not hold!" << std::endl;
160//     }
161//   }
162
163//   {
164//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
165//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
166//     ts.reset();
167//     int i=0;
168//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
169//     std::cout << "elapsed time: " << ts << std::endl;
170//     std::cout << "number of augmentation phases: " << i << std::endl;
171//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
172
173//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
174//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
175//      std::cout << "Slackness does not hold!" << std::endl;
176//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
177//      std::cout << "Slackness does not hold!" << std::endl;
178//     }
179//   }
180
181  ts.reset();
182  LPSolverWrapper lp;
183  lp.setMaximize();
184  typedef LPSolverWrapper::ColIt ColIt;
185  typedef LPSolverWrapper::RowIt RowIt;
186  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
187  EdgeIndexMap edge_index_map(g);
188  PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
189  Graph::EdgeIt e;
190  for (g.first(e); g.valid(e); g.next(e)) {
191    ColIt col_it=lp.addCol();
192    edge_index_map.set(e, col_it);
193    lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
194  }
195  Graph::NodeIt n;
196  for (g.first(n); g.valid(n); g.next(n)) {
197    if (n!=s) {
198      //hurokelek miatt
199      Graph::EdgeMap<int> coeffs(g, 0);
200      {
201        Graph::InEdgeIt e;
202        for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1);
203      }
204      {
205        Graph::OutEdgeIt e;
206        for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1);
207      }
208      if (n==t) {
209        Graph::EdgeIt e;
210        //std::vector< std::pair<ColIt, double> > row;
211        for (g.first(e); g.valid(e); g.next(e)) {
212          if (coeffs[e]!=0)
213            lp.setObjCoef(edge_index_map[e], coeffs[e]);
214        }
215      } else  {
216        RowIt row_it=lp.addRow();
217        Graph::EdgeIt e;
218        std::vector< std::pair<ColIt, double> > row;
219        for (g.first(e); g.valid(e); g.next(e)) {
220          if (coeffs[e]!=0)
221            row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
222        }       
223        lp.setRowCoeffs(row_it, row.begin(), row.end());
224        lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
225      }
226    }
227  }
228  lp.solveSimplex();
229  std::cout << "flow value: "<< lp.getObjVal() << std::endl;
230  std::cout << "elapsed time: " << ts << std::endl;
231
232  return 0;
233}
Note: See TracBrowser for help on using the repository browser.