marci@1104: // -*- c++ -*-
marci@1104: #include <iostream>
marci@1104: #include <fstream>
marci@1104: 
marci@1104: #include <lemon/smart_graph.h>
marci@1104: #include <lemon/list_graph.h>
marci@1104: #include <lemon/dimacs.h>
marci@1104: #include <lemon/time_measure.h>
marci@1104: #include <lp_solver_wrapper_3.h>
marci@1104: 
marci@1104: using std::cout;
marci@1104: using std::endl;
marci@1104: using namespace lemon;
marci@1104: 
marci@1104: template<typename Edge, typename EdgeIndexMap> 
marci@1104: class PrimalMap {
marci@1104: protected:
marci@1104:   LPGLPK* lp;
marci@1104:   EdgeIndexMap* edge_index_map;
marci@1104: public:
marci@1104:   PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
marci@1104:     lp(&_lp), edge_index_map(&_edge_index_map) { }
marci@1104:   double operator[](Edge e) const { 
marci@1104:     return lp->getPrimal((*edge_index_map)[e]);
marci@1104:   }
marci@1104: };
marci@1104: 
marci@1104: // Use a DIMACS max flow file as stdin.
marci@1104: // max_flow_expresion < dimacs_max_flow_file
marci@1104: 
marci@1104: int main(int, char **) {
marci@1104: 
marci@1104:   typedef ListGraph Graph;
marci@1104:   typedef Graph::Node Node;
marci@1104:   typedef Graph::Edge Edge;
marci@1104:   typedef Graph::EdgeIt EdgeIt;
marci@1104: 
marci@1104:   Graph g;
marci@1104:   
marci@1104:   Node s, t;
marci@1104:   Graph::EdgeMap<int> cap(g);
marci@1104:   readDimacs(std::cin, g, cap, s, t);
marci@1104:   Timer ts;
marci@1104:   
marci@1104:   typedef LPGLPK LPSolver;
marci@1104:   LPSolver lp;
marci@1104:   lp.setMaximize();
marci@1104:   typedef LPSolver::ColIt ColIt;
marci@1104:   typedef LPSolver::RowIt RowIt;
marci@1104:   typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
marci@1104:   EdgeIndexMap edge_index_map(g);
marci@1104:   PrimalMap<Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
marci@1104: 
marci@1104:   // capacity function
marci@1104:   for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
marci@1104:     ColIt col_it=lp.addCol();
marci@1104:     edge_index_map.set(e, col_it);
marci@1104:     if (cap[e]==0)
marci@1104:       lp.setColBounds(col_it, LPSolver::FIXED, 0, cap[e]);
marci@1104:     else 
marci@1104:       lp.setColBounds(col_it, LPSolver::DOUBLE, 0, cap[e]);
marci@1104:   }
marci@1104:   
marci@1104:   for (Graph::NodeIt n(g); n!=INVALID; ++n) {
marci@1104:     LPSolver::Expression expr;
marci@1104:     for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104:       expr+=edge_index_map[e];
marci@1104:     for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
marci@1104:       expr-=edge_index_map[e];
marci@1104:     // cost function
marci@1104:     if (n==s) {
marci@1104:       lp.setObjCoeffs(expr);      
marci@1104:     }
marci@1104:     // flow conservation
marci@1104:     if ((n!=s) && (n!=t)) {
marci@1104:       RowIt row_it=lp.addRow();
marci@1104:       lp.setRowCoeffs(row_it, expr);
marci@1104:       lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
marci@1104:     }
marci@1104:   }
marci@1104:   lp.solveSimplex();
marci@1104:   //std::cout << lp.colNum() << std::endl;
marci@1104:   //std::cout << lp.rowNum() << std::endl;
marci@1104:   //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
marci@1104:   for (Graph::EdgeIt e(g); e!=INVALID; ++e) 
marci@1104:     flow.set(e, lp_flow[e]);
marci@1104: }