1.1 --- a/src/work/marci/lp/lp_solver_wrapper.h	Sat Nov 20 14:23:27 2004 +0000
     1.2 +++ b/src/work/marci/lp/lp_solver_wrapper.h	Sat Nov 20 16:12:47 2004 +0000
     1.3 @@ -1,6 +1,6 @@
     1.4  // -*- c++ -*-
     1.5  #ifndef LEMON_LP_SOLVER_WRAPPER_H
     1.6 -#define LEMON_LP_SOLVER_WRAPPER
     1.7 +#define LEMON_LP_SOLVER_WRAPPER_H
     1.8  
     1.9  ///\ingroup misc
    1.10  ///\file

     2.1 --- a/src/work/marci/lp/max_flow_by_lp.cc	Sat Nov 20 14:23:27 2004 +0000
     2.2 +++ b/src/work/marci/lp/max_flow_by_lp.cc	Sat Nov 20 16:12:47 2004 +0000
     2.3 @@ -12,29 +12,101 @@
     2.4  //#include <preflow_res.h>
     2.5  #include <lp_solver_wrapper.h>
     2.6  
     2.7 -using namespace lemon;
     2.8 -
     2.9  // Use a DIMACS max flow file as stdin.
    2.10  // max_flow_demo < dimacs_max_flow_file
    2.11  
    2.12 -template<typename Edge, typename EdgeIndexMap> 
    2.13 -class PrimalMap {
    2.14 -protected:
    2.15 -  LPSolverWrapper* lp;
    2.16 -  EdgeIndexMap* edge_index_map;
    2.17 -public:
    2.18 -  PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : 
    2.19 -    lp(&_lp), edge_index_map(&_edge_index_map) { }
    2.20 -  double operator[](Edge e) const { 
    2.21 -    return lp->getPrimal((*edge_index_map)[e]);
    2.22 -  }
    2.23 -};
    2.24 +using namespace lemon;
    2.25 +
    2.26 +namespace lemon {
    2.27 +
    2.28 +  template<typename Edge, typename EdgeIndexMap> 
    2.29 +  class PrimalMap {
    2.30 +  protected:
    2.31 +    LPSolverWrapper* lp;
    2.32 +    EdgeIndexMap* edge_index_map;
    2.33 +  public:
    2.34 +    PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : 
    2.35 +      lp(&_lp), edge_index_map(&_edge_index_map) { }
    2.36 +    double operator[](Edge e) const { 
    2.37 +      return lp->getPrimal((*edge_index_map)[e]);
    2.38 +    }
    2.39 +  };
    2.40 +
    2.41 +  // excess: rho-delta
    2.42 +  template <typename Graph, typename Num,
    2.43 +	    typename Excess=typename Graph::template NodeMap<Num>, 
    2.44 +	    typename LCapMap=typename Graph::template EdgeMap<Num>,
    2.45 +	    typename CapMap=typename Graph::template EdgeMap<Num>,
    2.46 +            typename FlowMap=typename Graph::template EdgeMap<Num>,
    2.47 +	    typename CostMap=typename Graph::template EdgeMap<Num> >
    2.48 +  class MinCostGenFlow {
    2.49 +  protected:
    2.50 +    const Graph& g;
    2.51 +    const Excess& excess;
    2.52 +    const LCapMap& lcapacity;
    2.53 +    const CapMap& capacity;
    2.54 +    FlowMap& flow;
    2.55 +    const CostMap& cost;
    2.56 +  public:
    2.57 +    MinCostGenFlow(const Graph& _g, const Excess& _excess, 
    2.58 +		   const LCapMap& _lcapacity, const CapMap& _capacity, 
    2.59 +		   FlowMap& _flow, 
    2.60 +		   const CostMap& _cost) :
    2.61 +      g(_g), excess(_excess), lcapacity(_lcapacity),
    2.62 +      capacity(_capacity), flow(_flow), cost(_cost) { }
    2.63 +    void run() {
    2.64 +      LPSolverWrapper lp;
    2.65 +      lp.setMinimize();
    2.66 +      typedef LPSolverWrapper::ColIt ColIt;
    2.67 +      typedef LPSolverWrapper::RowIt RowIt;
    2.68 +      typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
    2.69 +      EdgeIndexMap edge_index_map(g);
    2.70 +      PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
    2.71 +      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
    2.72 +	ColIt col_it=lp.addCol();
    2.73 +	edge_index_map.set(e, col_it);
    2.74 +	if (lcapacity[e]==capacity[e])
    2.75 +	  lp.setColBounds(col_it, LPX_FX, lcapacity[e], capacity[e]);
    2.76 +	else 
    2.77 +	  lp.setColBounds(col_it, LPX_DB, lcapacity[e], capacity[e]);
    2.78 +	lp.setObjCoef(col_it, cost[e]);
    2.79 +      }
    2.80 +      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
    2.81 +	typename Graph::template EdgeMap<Num> coeffs(g, 0);
    2.82 +	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
    2.83 +	coeffs.set(e, coeffs[e]+1);
    2.84 +	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
    2.85 +	coeffs.set(e, coeffs[e]-1);
    2.86 +	RowIt row_it=lp.addRow();
    2.87 +	typename std::vector< std::pair<ColIt, double> > row;
    2.88 +	//std::cout << "node:" <<g.id(n)<<std::endl;
    2.89 +	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
    2.90 +	  if (coeffs[e]!=0) {
    2.91 +	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
    2.92 +	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
    2.93 +	  }
    2.94 +	}
    2.95 +	//std::cout << std::endl;
    2.96 +	lp.setRowCoeffs(row_it, row.begin(), row.end());
    2.97 +	lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
    2.98 +      }
    2.99 +      lp.solveSimplex();
   2.100 +      //std::cout << lp.colNum() << std::endl;
   2.101 +      //std::cout << lp.rowNum() << std::endl;
   2.102 +      //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
   2.103 +      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) 
   2.104 +      flow.set(e, lp_flow[e]);
   2.105 +    }
   2.106 +  };
   2.107 +
   2.108 +}
   2.109  
   2.110  int main(int, char **) {
   2.111  
   2.112    typedef ListGraph MutableGraph;
   2.113    typedef SmartGraph Graph;
   2.114    typedef Graph::Node Node;
   2.115 +  typedef Graph::Edge Edge;
   2.116    typedef Graph::EdgeIt EdgeIt;
   2.117  
   2.118    Graph g;
   2.119 @@ -60,7 +132,7 @@
   2.120      std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
   2.121      max_flow_test.minCut(cut);
   2.122  
   2.123 -    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.124 +    for (EdgeIt e(g); e!=INVALID; ++e) {
   2.125        if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
   2.126  	std::cout << "Slackness does not hold!" << std::endl;
   2.127        if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
   2.128 @@ -95,7 +167,7 @@
   2.129  
   2.130    {
   2.131      std::cout << "physical blocking flow augmentation ..." << std::endl;
   2.132 -    for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
   2.133 +    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
   2.134      ts.reset();
   2.135      int i=0;
   2.136      while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
   2.137 @@ -103,7 +175,7 @@
   2.138      std::cout << "number of augmentation phases: " << i << std::endl; 
   2.139      std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
   2.140  
   2.141 -    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.142 +    for (EdgeIt e(g); e!=INVALID; ++e) {
   2.143        if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
   2.144  	std::cout << "Slackness does not hold!" << std::endl;
   2.145        if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
   2.146 @@ -124,7 +196,7 @@
   2.147  
   2.148    {
   2.149      std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
   2.150 -    for (Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
   2.151 +    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
   2.152      ts.reset();
   2.153      int i=0;
   2.154      while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
   2.155 @@ -132,7 +204,7 @@
   2.156      std::cout << "number of augmentation phases: " << i << std::endl; 
   2.157      std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
   2.158  
   2.159 -    for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.160 +    for (EdgeIt e(g); e!=INVALID; ++e) {
   2.161        if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
   2.162  	std::cout << "Slackness does not hold!" << std::endl;
   2.163        if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
   2.164 @@ -177,46 +249,20 @@
   2.165  //   }
   2.166  
   2.167    ts.reset();
   2.168 -  LPSolverWrapper lp;
   2.169 -  lp.setMaximize();
   2.170 -  typedef LPSolverWrapper::ColIt ColIt;
   2.171 -  typedef LPSolverWrapper::RowIt RowIt;
   2.172 -  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
   2.173 -  EdgeIndexMap edge_index_map(g);
   2.174 -  PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
   2.175 -  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.176 -    ColIt col_it=lp.addCol();
   2.177 -    edge_index_map.set(e, col_it);
   2.178 -    lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
   2.179 -  }
   2.180 -  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
   2.181 -    if (n!=s) {
   2.182 -      //hurokelek miatt
   2.183 -      Graph::EdgeMap<int> coeffs(g, 0);
   2.184 -      for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e) 
   2.185 -	coeffs.set(e, coeffs[e]+1);
   2.186 -      for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
   2.187 -	coeffs.set(e, coeffs[e]-1);
   2.188 -      if (n==t) {
   2.189 -	for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.190 -	  if (coeffs[e]!=0) 
   2.191 -	    lp.setObjCoef(edge_index_map[e], coeffs[e]);
   2.192 -	}
   2.193 -      } else  {
   2.194 -	RowIt row_it=lp.addRow();
   2.195 -	std::vector< std::pair<ColIt, double> > row;
   2.196 -	for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
   2.197 -	  if (coeffs[e]!=0) 
   2.198 -	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
   2.199 -	}	
   2.200 -	lp.setRowCoeffs(row_it, row.begin(), row.end());
   2.201 -	lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
   2.202 -      }
   2.203 -    }
   2.204 -  }
   2.205 -  lp.solveSimplex();
   2.206 -  std::cout << "flow value: "<< lp.getObjVal() << std::endl;
   2.207 +
   2.208 +  Edge e=g.addEdge(t, s);
   2.209 +  Graph::EdgeMap<int> cost(g, 0);
   2.210 +  cost.set(e, -1);
   2.211 +  cap.set(e, 10000);
   2.212 +  typedef ConstMap<Node, int> Excess;
   2.213 +  Excess excess(0);
   2.214 +  typedef ConstMap<Edge, int> LCap;
   2.215 +  LCap lcap(0);
   2.216 +
   2.217 +  MinCostGenFlow<Graph, int, Excess, LCap> 
   2.218 +    min_cost(g, excess, lcap, cap, flow, cost); 
   2.219 +  min_cost.run();
   2.220 +
   2.221    std::cout << "elapsed time: " << ts << std::endl;
   2.222 -
   2.223 -  return 0;
   2.224 +  std::cout << "flow value: "<< flow[e] << std::endl;
   2.225  }