lemon-0.x: src/work/marci/lp/max_flow_by

     1 // -*- c++ -*-

     2 #include <iostream>

     3 #include <fstream>

     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>

    16 using namespace hugo;

    18 // Use a DIMACS max flow file as stdin.

    19 // max_flow_demo < dimacs_max_flow_file

    21 template<typename Edge, typename EdgeIndexMap>

    22 class PrimalMap {

    23 protected:

    24   LPSolverWrapper* lp;

    25   EdgeIndexMap* edge_index_map;

    26 public:

    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 };

    34 int main(int, char **) {

    36   typedef SageGraph MutableGraph;

    37   //typedef SmartGraph Graph;

    38   typedef SageGraph Graph;

    39   typedef Graph::Node Node;

    40   typedef Graph::EdgeIt EdgeIt;

    42   Graph g;

    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);

    55   Graph::NodeMap<bool> cut(g);

    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);

    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   }

    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;

    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 //   }

    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 //   }

    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;

   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   }

   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 //   }

   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;

   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   }

   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;

   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 //   }

   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;

   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 //   }

   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;

   232   return 0;

   233 }

author	marci
	Wed, 22 Sep 2004 10:47:59 +0000
changeset 901	69a8e672acb1
child 921	818510fa3d99
permissions	-rw-r--r--