1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/marci/lp/max_flow_expression.cc Fri Jan 28 14:33:32 2005 +0000
1.3 @@ -0,0 +1,87 @@
1.4 +// -*- c++ -*-
1.5 +#include <iostream>
1.6 +#include <fstream>
1.7 +
1.8 +#include <lemon/smart_graph.h>
1.9 +#include <lemon/list_graph.h>
1.10 +#include <lemon/dimacs.h>
1.11 +#include <lemon/time_measure.h>
1.12 +#include <lp_solver_wrapper_3.h>
1.13 +
1.14 +using std::cout;
1.15 +using std::endl;
1.16 +using namespace lemon;
1.17 +
1.18 +template<typename Edge, typename EdgeIndexMap>
1.19 +class PrimalMap {
1.20 +protected:
1.21 + LPGLPK* lp;
1.22 + EdgeIndexMap* edge_index_map;
1.23 +public:
1.24 + PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
1.25 + lp(&_lp), edge_index_map(&_edge_index_map) { }
1.26 + double operator[](Edge e) const {
1.27 + return lp->getPrimal((*edge_index_map)[e]);
1.28 + }
1.29 +};
1.30 +
1.31 +// Use a DIMACS max flow file as stdin.
1.32 +// max_flow_expresion < dimacs_max_flow_file
1.33 +
1.34 +int main(int, char **) {
1.35 +
1.36 + typedef ListGraph Graph;
1.37 + typedef Graph::Node Node;
1.38 + typedef Graph::Edge Edge;
1.39 + typedef Graph::EdgeIt EdgeIt;
1.40 +
1.41 + Graph g;
1.42 +
1.43 + Node s, t;
1.44 + Graph::EdgeMap<int> cap(g);
1.45 + readDimacs(std::cin, g, cap, s, t);
1.46 + Timer ts;
1.47 +
1.48 + typedef LPGLPK LPSolver;
1.49 + LPSolver lp;
1.50 + lp.setMaximize();
1.51 + typedef LPSolver::ColIt ColIt;
1.52 + typedef LPSolver::RowIt RowIt;
1.53 + typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
1.54 + EdgeIndexMap edge_index_map(g);
1.55 + PrimalMap<Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
1.56 +
1.57 + // capacity function
1.58 + for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.59 + ColIt col_it=lp.addCol();
1.60 + edge_index_map.set(e, col_it);
1.61 + if (cap[e]==0)
1.62 + lp.setColBounds(col_it, LPSolver::FIXED, 0, cap[e]);
1.63 + else
1.64 + lp.setColBounds(col_it, LPSolver::DOUBLE, 0, cap[e]);
1.65 + }
1.66 +
1.67 + for (Graph::NodeIt n(g); n!=INVALID; ++n) {
1.68 + LPSolver::Expression expr;
1.69 + for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
1.70 + expr+=edge_index_map[e];
1.71 + for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
1.72 + expr-=edge_index_map[e];
1.73 + // cost function
1.74 + if (n==s) {
1.75 + lp.setObjCoeffs(expr);
1.76 + }
1.77 + // flow conservation
1.78 + if ((n!=s) && (n!=t)) {
1.79 + RowIt row_it=lp.addRow();
1.80 + lp.setRowCoeffs(row_it, expr);
1.81 + lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
1.82 + }
1.83 + }
1.84 + lp.solveSimplex();
1.85 + //std::cout << lp.colNum() << std::endl;
1.86 + //std::cout << lp.rowNum() << std::endl;
1.87 + //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
1.88 + for (Graph::EdgeIt e(g); e!=INVALID; ++e)
1.89 + flow.set(e, lp_flow[e]);
1.90 +}