COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/lp/max_flow_by_lp.cc @ 1014:aae850a2394d

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

Modifications for hugo 0.2

File size: 7.6 KB
RevLine 
[764]1// -*- c++ -*-
2#include <iostream>
3#include <fstream>
4
[921]5#include <lemon/smart_graph.h>
[1014]6#include <lemon/list_graph.h>
[921]7#include <lemon/dimacs.h>
8#include <lemon/time_measure.h>
[764]9//#include <graph_wrapper.h>
[1014]10#include <lemon/preflow.h>
[764]11#include <augmenting_flow.h>
12//#include <preflow_res.h>
13#include <lp_solver_wrapper.h>
14
[921]15using namespace lemon;
[764]16
17// Use a DIMACS max flow file as stdin.
18// max_flow_demo < dimacs_max_flow_file
19
20template<typename Edge, typename EdgeIndexMap>
21class PrimalMap {
22protected:
23  LPSolverWrapper* lp;
24  EdgeIndexMap* edge_index_map;
25public:
26  PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) :
27    lp(&_lp), edge_index_map(&_edge_index_map) { }
28  double operator[](Edge e) const {
29    return lp->getPrimal((*edge_index_map)[e]);
30  }
31};
32
33int main(int, char **) {
34
[1014]35  typedef ListGraph MutableGraph;
36  typedef SmartGraph Graph;
[764]37  typedef Graph::Node Node;
38  typedef Graph::EdgeIt EdgeIt;
39
40  Graph g;
41
42  Node s, t;
43  Graph::EdgeMap<int> cap(g);
44  //readDimacsMaxFlow(std::cin, g, s, t, cap);
45  readDimacs(std::cin, g, cap, s, t);
46  Timer ts;
47  Graph::EdgeMap<int> flow(g); //0 flow
[1014]48  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
[764]49    max_flow_test(g, s, t, cap, flow);
50  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
51    augmenting_flow_test(g, s, t, cap, flow);
52 
53  Graph::NodeMap<bool> cut(g);
54
55  {
56    std::cout << "preflow ..." << std::endl;
57    ts.reset();
58    max_flow_test.run();
59    std::cout << "elapsed time: " << ts << std::endl;
60    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
[1014]61    max_flow_test.minCut(cut);
[764]62
[1014]63    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[986]64      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]65        std::cout << "Slackness does not hold!" << std::endl;
[986]66      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]67        std::cout << "Slackness does not hold!" << std::endl;
68    }
69  }
70
71//   {
72//     std::cout << "preflow ..." << std::endl;
73//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
74//     ts.reset();
[1014]75//     max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
[764]76//     std::cout << "elapsed time: " << ts << std::endl;
77//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
78
79//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
[986]80//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]81//      std::cout << "Slackness does not hold!" << std::endl;
[986]82//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]83//      std::cout << "Slackness does not hold!" << std::endl;
84//     }
85//   }
86
87//   {
88//     std::cout << "wrapped preflow ..." << std::endl;
89//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
90//     ts.reset();
91//     pre_flow_res.run();
92//     std::cout << "elapsed time: " << ts << std::endl;
93//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
94//   }
95
96  {
97    std::cout << "physical blocking flow augmentation ..." << std::endl;
[1014]98    for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
[764]99    ts.reset();
100    int i=0;
101    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
102    std::cout << "elapsed time: " << ts << std::endl;
103    std::cout << "number of augmentation phases: " << i << std::endl;
104    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
105
[1014]106    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[986]107      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]108        std::cout << "Slackness does not hold!" << std::endl;
[986]109      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]110        std::cout << "Slackness does not hold!" << std::endl;
111    }
112  }
113
114//   {
115//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
116//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
117//     ts.reset();
118//     int i=0;
119//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
120//     std::cout << "elapsed time: " << ts << std::endl;
121//     std::cout << "number of augmentation phases: " << i << std::endl;
122//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
123//   }
124
125  {
126    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
[1014]127    for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
[764]128    ts.reset();
129    int i=0;
130    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
131    std::cout << "elapsed time: " << ts << std::endl;
132    std::cout << "number of augmentation phases: " << i << std::endl;
133    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
134
[1014]135    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[986]136      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]137        std::cout << "Slackness does not hold!" << std::endl;
[986]138      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]139        std::cout << "Slackness does not hold!" << std::endl;
140    }
141  }
142
143//   {
144//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
145//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
146//     ts.reset();
147//     int i=0;
148//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
149//     std::cout << "elapsed time: " << ts << std::endl;
150//     std::cout << "number of augmentation phases: " << i << std::endl;
151//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
152
153//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
[986]154//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]155//      std::cout << "Slackness does not hold!" << std::endl;
[986]156//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]157//      std::cout << "Slackness does not hold!" << std::endl;
158//     }
159//   }
160
161//   {
162//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
163//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
164//     ts.reset();
165//     int i=0;
166//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
167//     std::cout << "elapsed time: " << ts << std::endl;
168//     std::cout << "number of augmentation phases: " << i << std::endl;
169//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
170
171//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
[986]172//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e])
[764]173//      std::cout << "Slackness does not hold!" << std::endl;
[986]174//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0)
[764]175//      std::cout << "Slackness does not hold!" << std::endl;
176//     }
177//   }
178
179  ts.reset();
180  LPSolverWrapper lp;
181  lp.setMaximize();
182  typedef LPSolverWrapper::ColIt ColIt;
183  typedef LPSolverWrapper::RowIt RowIt;
184  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
185  EdgeIndexMap edge_index_map(g);
186  PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
[1014]187  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[764]188    ColIt col_it=lp.addCol();
189    edge_index_map.set(e, col_it);
190    lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
191  }
[1014]192  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
[764]193    if (n!=s) {
194      //hurokelek miatt
195      Graph::EdgeMap<int> coeffs(g, 0);
[1014]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);
[764]200      if (n==t) {
[1014]201        for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[764]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;
[1014]208        for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
[764]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;
220
221  return 0;
222}
Note: See TracBrowser for help on using the repository browser.