marci@764: // -*- c++ -*- marci@764: #include marci@764: #include marci@764: marci@764: #include alpar@921: #include alpar@921: #include alpar@921: #include marci@764: //#include alpar@921: #include marci@764: #include marci@764: //#include marci@764: #include marci@764: #include marci@764: alpar@921: using namespace lemon; marci@764: marci@764: // Use a DIMACS max flow file as stdin. marci@764: // max_flow_demo < dimacs_max_flow_file marci@764: marci@764: template marci@764: class PrimalMap { marci@764: protected: marci@764: LPSolverWrapper* lp; marci@764: EdgeIndexMap* edge_index_map; marci@764: public: marci@764: PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : marci@764: lp(&_lp), edge_index_map(&_edge_index_map) { } marci@764: double operator[](Edge e) const { marci@764: return lp->getPrimal((*edge_index_map)[e]); marci@764: } marci@764: }; marci@764: marci@764: int main(int, char **) { marci@764: marci@764: typedef SageGraph MutableGraph; marci@764: //typedef SmartGraph Graph; marci@764: typedef SageGraph Graph; marci@764: typedef Graph::Node Node; marci@764: typedef Graph::EdgeIt EdgeIt; marci@764: marci@764: Graph g; marci@764: marci@764: Node s, t; marci@764: Graph::EdgeMap cap(g); marci@764: //readDimacsMaxFlow(std::cin, g, s, t, cap); marci@764: readDimacs(std::cin, g, cap, s, t); marci@764: Timer ts; marci@764: Graph::EdgeMap flow(g); //0 flow marci@764: MaxFlow, Graph::EdgeMap > marci@764: max_flow_test(g, s, t, cap, flow); marci@764: AugmentingFlow, Graph::EdgeMap > marci@764: augmenting_flow_test(g, s, t, cap, flow); marci@764: marci@764: Graph::NodeMap cut(g); marci@764: marci@764: { marci@764: std::cout << "preflow ..." << std::endl; marci@764: ts.reset(); marci@764: max_flow_test.run(); marci@764: std::cout << "elapsed time: " << ts << std::endl; marci@764: std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl; marci@764: max_flow_test.actMinCut(cut); marci@764: marci@764: FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: } marci@764: } marci@764: marci@764: // { marci@764: // std::cout << "preflow ..." << std::endl; marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: // ts.reset(); marci@764: // max_flow_test.preflow(MaxFlow, Graph::EdgeMap >::GEN_FLOW); marci@764: // std::cout << "elapsed time: " << ts << std::endl; marci@764: // std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl; marci@764: marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: // if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // } marci@764: // } marci@764: marci@764: // { marci@764: // std::cout << "wrapped preflow ..." << std::endl; marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: // ts.reset(); marci@764: // pre_flow_res.run(); marci@764: // std::cout << "elapsed time: " << ts << std::endl; marci@764: // std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl; marci@764: // } marci@764: marci@764: { marci@764: std::cout << "physical blocking flow augmentation ..." << std::endl; marci@764: FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: ts.reset(); marci@764: int i=0; marci@764: while (augmenting_flow_test.augmentOnBlockingFlow()) { ++i; } marci@764: std::cout << "elapsed time: " << ts << std::endl; marci@764: std::cout << "number of augmentation phases: " << i << std::endl; marci@764: std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; marci@764: marci@764: FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: } marci@764: } marci@764: marci@764: // { marci@764: // std::cout << "faster physical blocking flow augmentation ..." << std::endl; marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: // ts.reset(); marci@764: // int i=0; marci@764: // while (max_flow_test.augmentOnBlockingFlow1()) { ++i; } marci@764: // std::cout << "elapsed time: " << ts << std::endl; marci@764: // std::cout << "number of augmentation phases: " << i << std::endl; marci@764: // std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl; marci@764: // } marci@764: marci@764: { marci@764: std::cout << "on-the-fly blocking flow augmentation ..." << std::endl; marci@764: FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: ts.reset(); marci@764: int i=0; marci@764: while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; } marci@764: std::cout << "elapsed time: " << ts << std::endl; marci@764: std::cout << "number of augmentation phases: " << i << std::endl; marci@764: std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; marci@764: marci@764: FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: std::cout << "Slackness does not hold!" << std::endl; marci@764: } marci@764: } marci@764: marci@764: // { marci@764: // std::cout << "on-the-fly shortest path augmentation ..." << std::endl; marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: // ts.reset(); marci@764: // int i=0; marci@764: // while (augmenting_flow_test.augmentOnShortestPath()) { ++i; } marci@764: // std::cout << "elapsed time: " << ts << std::endl; marci@764: // std::cout << "number of augmentation phases: " << i << std::endl; marci@764: // std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; marci@764: marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: // if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // } marci@764: // } marci@764: marci@764: // { marci@764: // std::cout << "on-the-fly shortest path augmentation ..." << std::endl; marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0); marci@764: // ts.reset(); marci@764: // int i=0; marci@764: // while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; } marci@764: // std::cout << "elapsed time: " << ts << std::endl; marci@764: // std::cout << "number of augmentation phases: " << i << std::endl; marci@764: // std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl; marci@764: marci@764: // FOR_EACH_LOC(Graph::EdgeIt, e, g) { marci@764: // if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) marci@764: // std::cout << "Slackness does not hold!" << std::endl; marci@764: // } marci@764: // } marci@764: marci@764: ts.reset(); marci@764: LPSolverWrapper lp; marci@764: lp.setMaximize(); marci@764: typedef LPSolverWrapper::ColIt ColIt; marci@764: typedef LPSolverWrapper::RowIt RowIt; marci@764: typedef Graph::EdgeMap EdgeIndexMap; marci@764: EdgeIndexMap edge_index_map(g); marci@764: PrimalMap lp_flow(lp, edge_index_map); marci@764: Graph::EdgeIt e; marci@764: for (g.first(e); g.valid(e); g.next(e)) { marci@764: ColIt col_it=lp.addCol(); marci@764: edge_index_map.set(e, col_it); marci@764: lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]); marci@764: } marci@764: Graph::NodeIt n; marci@764: for (g.first(n); g.valid(n); g.next(n)) { marci@764: if (n!=s) { marci@764: //hurokelek miatt marci@764: Graph::EdgeMap coeffs(g, 0); marci@764: { marci@764: Graph::InEdgeIt e; marci@764: for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1); marci@764: } marci@764: { marci@764: Graph::OutEdgeIt e; marci@764: for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1); marci@764: } marci@764: if (n==t) { marci@764: Graph::EdgeIt e; marci@764: //std::vector< std::pair > row; marci@764: for (g.first(e); g.valid(e); g.next(e)) { marci@764: if (coeffs[e]!=0) marci@764: lp.setObjCoef(edge_index_map[e], coeffs[e]); marci@764: } marci@764: } else { marci@764: RowIt row_it=lp.addRow(); marci@764: Graph::EdgeIt e; marci@764: std::vector< std::pair > row; marci@764: for (g.first(e); g.valid(e); g.next(e)) { marci@764: if (coeffs[e]!=0) marci@764: row.push_back(std::make_pair(edge_index_map[e], coeffs[e])); marci@764: } marci@764: lp.setRowCoeffs(row_it, row.begin(), row.end()); marci@764: lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0); marci@764: } marci@764: } marci@764: } marci@764: lp.solveSimplex(); marci@764: std::cout << "flow value: "<< lp.getObjVal() << std::endl; marci@764: std::cout << "elapsed time: " << ts << std::endl; marci@764: marci@764: return 0; marci@764: }