// -*- c++ -*-

 #include <iostream>

 #include <fstream>

 #include <lemon/graph_utils.h>

 #include <lemon/smart_graph.h>

 #include <lemon/list_graph.h>

 #include <lemon/dimacs.h>

 #include <lemon/time_measure.h>

 #include <lp_solver_glpk.h>

 using std::cout;

 using std::endl;

 using namespace lemon;

 template<typename Edge, typename EdgeIndexMap> 

 class PrimalMap {

 protected:

   LpGlpk* lp;

   EdgeIndexMap* edge_index_map;

 public:

   PrimalMap(LpGlpk& _lp, EdgeIndexMap& _edge_index_map) : 

     lp(&_lp), edge_index_map(&_edge_index_map) { }

   double operator[](Edge e) const { 

     return lp->getPrimal((*edge_index_map)[e]);

   }

 };

 // Use a DIMACS max flow file as stdin.

 // max_flow_expresion < dimacs_max_flow_file

 int main(int, char **) {

   typedef ListGraph Graph;

   typedef Graph::Node Node;

   typedef Graph::Edge Edge;

   typedef Graph::EdgeIt EdgeIt;

   Graph g;

   Node s, t;

   Graph::EdgeMap<int> cap(g);

   readDimacs(std::cin, g, cap, s, t);

   Timer ts;

   typedef LpGlpk LPSolver;

   LPSolver lp;

   lp.setMaximize();

   typedef LPSolver::Col Col;

   typedef LPSolver::Row Row;

   typedef Graph::EdgeMap<Col> EdgeIndexMap;

   typedef Graph::NodeMap<Row> NodeIndexMap;

   EdgeIndexMap edge_index_map(g);

   NodeIndexMap node_index_map(g);

   PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);

   // nonnegativity of flow and capacity function

   for (Graph::EdgeIt e(g); e!=INVALID; ++e) {

     Col col=lp.addCol();

     edge_index_map.set(e, col);

     // interesting property in GLPK:

     // if you change the order of the following two lines, the 

     // two runs of GLPK are extremely different

       lp.setColLowerBound(col, 0);

       lp.setColUpperBound(col, cap[e]);

   }

   for (Graph::NodeIt n(g); n!=INVALID; ++n) {

     LPSolver::Expression expr;

     for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)

       expr+=edge_index_map[e];

     for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)

       expr-=edge_index_map[e];

     // cost function

     if (n==s) {

       lp.setObjCoeffs(expr);      

     }

     // flow conservation constraints

     if ((n!=s) && (n!=t)) {

       node_index_map.set(n, lp.addRow(expr == 0.0));

     }

   }

   lp.solveSimplex();

   cout << "elapsed time: " << ts << endl;

 //   cout << "rows:" << endl;

 //   for (

 //        LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);

 //        i!=INVALID;

 //        ++i) { 

 //     cout << i << " ";

 //   }

 //   cout << endl;

 //   cout << "cols:" << endl;

 //   for (

 //        LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);

 //        i!=INVALID;

 //        ++i) { 

 //     cout << i << " ";

 //   }

 //   cout << endl;

   lp.setMIP();

   cout << "elapsed time: " << ts << endl;

   for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {

     lp.setColInt(it);

   }

   cout << "elapsed time: " << ts << endl;

   lp.solveBandB();

   cout << "elapsed time: " << ts << endl;

 }