src/work/marci/lp/min_cost_gen_flow.h
author marci
Mon, 06 Dec 2004 11:56:10 +0000
changeset 1031 0b7169db694f
parent 1028 2497336d7e14
child 1033 9fff45a59e92
permissions -rw-r--r--
:-(
     1 // -*- c++ -*-
     2 #ifndef LEMON_MIN_COST_GEN_FLOW_H
     3 #define LEMON_MIN_COST_GEN_FLOW_H
     4 #include <iostream>
     5 //#include <fstream>
     6 
     7 #include <lemon/smart_graph.h>
     8 #include <lemon/list_graph.h>
     9 //#include <lemon/dimacs.h>
    10 //#include <lemon/time_measure.h>
    11 //#include <graph_wrapper.h>
    12 #include <lemon/preflow.h>
    13 #include <lemon/min_cost_flow.h>
    14 //#include <augmenting_flow.h>
    15 //#include <preflow_res.h>
    16 #include <work/marci/merge_node_graph_wrapper.h>
    17 #include <work/marci/lp/lp_solver_wrapper_2.h>
    18 
    19 namespace lemon {
    20 
    21   template<typename Edge, typename EdgeIndexMap> 
    22   class PrimalMap {
    23   protected:
    24     LPSolverWrapper* lp;
    25     EdgeIndexMap* edge_index_map;
    26   public:
    27     PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : 
    28       lp(&_lp), edge_index_map(&_edge_index_map) { }
    29     double operator[](Edge e) const { 
    30       return lp->getPrimal((*edge_index_map)[e]);
    31     }
    32   };
    33 
    34   // excess: rho-delta egyelore csak =0-ra.
    35   template <typename Graph, typename Num,
    36 	    typename Excess=typename Graph::template NodeMap<Num>, 
    37 	    typename LCapMap=typename Graph::template EdgeMap<Num>,
    38 	    typename CapMap=typename Graph::template EdgeMap<Num>,
    39             typename FlowMap=typename Graph::template EdgeMap<Num>,
    40 	    typename CostMap=typename Graph::template EdgeMap<Num> >
    41   class MinCostGenFlow {
    42   protected:
    43     const Graph& g;
    44     const Excess& excess;
    45     const LCapMap& lcapacity;
    46     const CapMap& capacity;
    47     FlowMap& flow;
    48     const CostMap& cost;
    49   public:
    50     MinCostGenFlow(const Graph& _g, const Excess& _excess, 
    51 		   const LCapMap& _lcapacity, const CapMap& _capacity, 
    52 		   FlowMap& _flow, 
    53 		   const CostMap& _cost) :
    54       g(_g), excess(_excess), lcapacity(_lcapacity),
    55       capacity(_capacity), flow(_flow), cost(_cost) { }
    56     bool feasible() {
    57       //      std::cout << "making new vertices..." << std::endl; 
    58       typedef ListGraph Graph2;
    59       Graph2 g2;
    60       typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
    61       //      std::cout << "merging..." << std::endl; 
    62       GW gw(g, g2);
    63       typename GW::Node s(INVALID, g2.addNode(), true);
    64       typename GW::Node t(INVALID, g2.addNode(), true);
    65       typedef SmartGraph Graph3;
    66       //      std::cout << "making extender graph..." << std::endl; 
    67       typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
    68 //       {
    69 // 	checkConcept<StaticGraph, GWW>();   
    70 //       }
    71       GWW gww(gw);
    72       typedef AugmentingGraphWrapper<GW, GWW> GWWW;
    73       GWWW gwww(gw, gww);
    74 
    75       //      std::cout << "making new edges..." << std::endl; 
    76       typename GWWW::template EdgeMap<Num> translated_cap(gwww);
    77 
    78       for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
    79 	translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
    80 			   capacity[e]-lcapacity[e]);
    81 	//	cout << "t_cap " << gw.id(e) << " " 
    82 	//	     << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
    83       }
    84 
    85       Num expected=0;
    86 
    87       //      std::cout << "making new edges 2..." << std::endl; 
    88       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
    89 	Num a=0;
    90 	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
    91 	  a+=lcapacity[e];
    92 	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
    93 	  a-=lcapacity[e];
    94 	if (excess[n]>a) {
    95 	  typename GWW::Edge e=
    96 	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
    97 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
    98 			     excess[n]-a);
    99 	  //	  std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
   100 	}
   101 	if (excess[n]<a) {
   102 	  typename GWW::Edge e=
   103 	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
   104 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
   105 			     a-excess[n]);
   106 	  expected+=a-excess[n];
   107 	  //	  std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
   108 	}
   109       }
   110 
   111       //      std::cout << "preflow..." << std::endl; 
   112       typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
   113       Preflow<GWWW, Num> preflow(gwww, s, t, 
   114 				 translated_cap, translated_flow);
   115       preflow.run();
   116       //      std::cout << "fv: " << preflow.flowValue() << std::endl; 
   117       //      std::cout << "expected: " << expected << std::endl; 
   118 
   119       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   120 	typename GW::Edge ew(e, INVALID, false);
   121 	typename GWWW::Edge ewww(ew, INVALID, false);
   122 	flow.set(e, translated_flow[ewww]+lcapacity[e]);
   123       }
   124       return (expected>=preflow.flowValue());
   125     }
   126     // for nonnegative costs
   127     bool run() {
   128       //      std::cout << "making new vertices..." << std::endl; 
   129       typedef ListGraph Graph2;
   130       Graph2 g2;
   131       typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
   132       //      std::cout << "merging..." << std::endl; 
   133       GW gw(g, g2);
   134       typename GW::Node s(INVALID, g2.addNode(), true);
   135       typename GW::Node t(INVALID, g2.addNode(), true);
   136       typedef SmartGraph Graph3;
   137       //      std::cout << "making extender graph..." << std::endl; 
   138       typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
   139 //       {
   140 // 	checkConcept<StaticGraph, GWW>();   
   141 //       }
   142       GWW gww(gw);
   143       typedef AugmentingGraphWrapper<GW, GWW> GWWW;
   144       GWWW gwww(gw, gww);
   145 
   146       //      std::cout << "making new edges..." << std::endl; 
   147       typename GWWW::template EdgeMap<Num> translated_cap(gwww);
   148 
   149       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   150 	typename GW::Edge ew(e, INVALID, false);
   151 	typename GWWW::Edge ewww(ew, INVALID, false);
   152 	translated_cap.set(ewww, capacity[e]-lcapacity[e]);
   153 	//	cout << "t_cap " << g.id(e) << " " 
   154 	//	     << translated_cap[ewww] << endl;
   155       }
   156 
   157       Num expected=0;
   158 
   159       //      std::cout << "making new edges 2..." << std::endl; 
   160       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
   161 	//	std::cout << "node: " << g.id(n) << std::endl;
   162 	Num a=0;
   163 	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
   164 	  a+=lcapacity[e];
   165 	  //	  std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
   166 	}
   167 	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
   168 	  a-=lcapacity[e];
   169 	  //	  std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
   170 	}
   171 	//	std::cout << "excess " << g.id(n) << ": " << a << std::endl;
   172 	if (0>a) {
   173 	  typename GWW::Edge e=
   174 	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
   175 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
   176 			     -a);
   177 	  //	  std::cout << g.id(n) << "->t " << -a << std::endl;
   178 	}
   179 	if (0<a) {
   180 	  typename GWW::Edge e=
   181 	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
   182 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
   183 			     a);
   184 	  expected+=a;
   185 	  //	  std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
   186 	}
   187       }
   188 
   189       //      std::cout << "preflow..." << std::endl; 
   190       typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
   191       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   192 	translated_cost.set(typename GWWW::Edge(
   193         typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
   194       }
   195       //      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
   196       MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>, 
   197       typename GWWW::template EdgeMap<Num> > 
   198       min_cost_flow(gwww, translated_cost, translated_cap, 
   199 		    s, t);
   200       while (min_cost_flow.augment()) { }
   201       std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; 
   202       std::cout << "expected: " << expected << std::endl; 
   203 
   204       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   205 	typename GW::Edge ew(e, INVALID, false);
   206 	typename GWWW::Edge ewww(ew, INVALID, false);
   207 	//	std::cout << g.id(e) << " " << flow[e] << std::endl;
   208 	flow.set(e, lcapacity[e]+
   209 		 min_cost_flow.getFlow()[ewww]);
   210       }
   211       return (expected>=min_cost_flow.flowValue());
   212     }
   213     void runByLP() {
   214       typedef LPSolverWrapper LPSolver;
   215       LPSolver lp;
   216       lp.setMinimize();
   217       typedef LPSolver::ColIt ColIt;
   218       typedef LPSolver::RowIt RowIt;
   219       typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
   220       EdgeIndexMap edge_index_map(g);
   221       PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
   222       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   223 	ColIt col_it=lp.addCol();
   224 	edge_index_map.set(e, col_it);
   225 	if (lcapacity[e]==capacity[e])
   226 	  lp.setColBounds(col_it, LPX_FX, lcapacity[e], capacity[e]);
   227 	else 
   228 	  lp.setColBounds(col_it, LPX_DB, lcapacity[e], capacity[e]);
   229 	lp.setObjCoef(col_it, cost[e]);
   230       }
   231       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
   232 	typename Graph::template EdgeMap<Num> coeffs(g, 0);
   233 	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
   234 	coeffs.set(e, coeffs[e]+1);
   235 	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
   236 	coeffs.set(e, coeffs[e]-1);
   237 	RowIt row_it=lp.addRow();
   238 	typename std::vector< std::pair<ColIt, double> > row;
   239 	//std::cout << "node:" <<g.id(n)<<std::endl;
   240 	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   241 	  if (coeffs[e]!=0) {
   242 	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
   243 	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
   244 	  }
   245 	}
   246 	//std::cout << std::endl;
   247 	lp.setRowCoeffs(row_it, row.begin(), row.end());
   248 	lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
   249       }
   250       lp.solveSimplex();
   251       //std::cout << lp.colNum() << std::endl;
   252       //std::cout << lp.rowNum() << std::endl;
   253       //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
   254       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) 
   255       flow.set(e, lp_flow[e]);
   256     }
   257   };
   258 
   259 } // namespace lemon
   260 
   261 #endif //LEMON_MIN_COST_GEN_FLOW_H