Changeset 1015:e3bb0e118bb4 in lemon0.x for src/work/marci/lp/max_flow_by_lp.cc
 Timestamp:
 11/20/04 17:12:47 (20 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@1405
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

src/work/marci/lp/max_flow_by_lp.cc
r1014 r1015 13 13 #include <lp_solver_wrapper.h> 14 14 15 using namespace lemon;16 17 15 // Use a DIMACS max flow file as stdin. 18 16 // max_flow_demo < dimacs_max_flow_file 19 17 20 template<typename Edge, typename EdgeIndexMap> 21 class PrimalMap { 22 protected: 23 LPSolverWrapper* lp; 24 EdgeIndexMap* edge_index_map; 25 public: 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 }; 18 using namespace lemon; 19 20 namespace 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: rhodelta 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 104 int main(int, char **) { … … 36 107 typedef SmartGraph Graph; 37 108 typedef Graph::Node Node; 109 typedef Graph::Edge Edge; 38 110 typedef Graph::EdgeIt EdgeIt; 39 111 … … 61 133 max_flow_test.minCut(cut); 62 134 63 for ( Graph::EdgeIt e(g); e!=INVALID; ++e) {135 for (EdgeIt e(g); e!=INVALID; ++e) { 64 136 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 65 137 std::cout << "Slackness does not hold!" << std::endl; … … 96 168 { 97 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 171 ts.reset(); 100 172 int i=0; … … 104 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 179 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 108 180 std::cout << "Slackness does not hold!" << std::endl; … … 125 197 { 126 198 std::cout << "onthefly 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 200 ts.reset(); 129 201 int i=0; … … 133 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 208 if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 137 209 std::cout << "Slackness does not hold!" << std::endl; … … 178 250 179 251 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); 187 for (Graph::EdgeIt e(g); e!=INVALID; ++e) { 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 } 192 for (Graph::NodeIt n(g); n!=INVALID; ++n) { 193 if (n!=s) { 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; 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 219 266 std::cout << "elapsed time: " << ts << std::endl; 220 221 return 0; 267 std::cout << "flow value: "<< flow[e] << std::endl; 222 268 }
Note: See TracChangeset
for help on using the changeset viewer.