lemon-0.x: src/work/marci/lp/min_cost_gen

     1 // -*- c++ -*-

     2 #ifndef LEMON_MIN_COST_GEN_FLOW_H

     3 #define LEMON_MIN_COST_GEN_FLOW_H

     4 #include <iostream>

     5 //#include <fstream>

     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_3.h>

    19 namespace lemon {

    21   template<typename Edge, typename EdgeIndexMap>

    22   class PrimalMap {

    23   protected:

    24     LPGLPK* lp;

    25     EdgeIndexMap* edge_index_map;

    26   public:

    27     PrimalMap(LPGLPK& _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   };

    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);

    75       //      std::cout << "making new edges..." << std::endl;

    76       typename GWWW::template EdgeMap<Num> translated_cap(gwww);

    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       }

    85       Num expected=0;

    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       }

   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;

   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 (preflow.flowValue()>=expected);

   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);

   146       //      std::cout << "making new edges..." << std::endl;

   147       typename GWWW::template EdgeMap<Num> translated_cap(gwww);

   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       }

   157       Num expected=0;

   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       }

   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;

   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 (min_cost_flow.flowValue()>=expected);

   212     }

   213     void runByLP() {

   214       typedef LPGLPK 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, LPSolver::FIXED, lcapacity[e], capacity[e]);

   227 	else

   228 	  lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);

   229 	lp.setObjCoef(col_it, cost[e]);

   230       }

   231       LPSolver::ColIt col_it;

   232       for (lp.col_iter_map.first(col_it, lp.VALID_CLASS);

   233 	   lp.col_iter_map.valid(col_it);

   234 	   lp.col_iter_map.next(col_it)) {

   235 //	std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;

   236       }

   237       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {

   238 	typename Graph::template EdgeMap<Num> coeffs(g, 0);

   239 	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)

   240 	coeffs.set(e, coeffs[e]+1);

   241 	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)

   242 	coeffs.set(e, coeffs[e]-1);

   243 	RowIt row_it=lp.addRow();

   244 	typename std::vector< std::pair<ColIt, double> > row;

   245 	//std::cout << "node:" <<g.id(n)<<std::endl;

   246 	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {

   247 	  if (coeffs[e]!=0) {

   248 	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];

   249 	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));

   250 	  }

   251 	}

   252 	//std::cout << std::endl;

   253 	//std::cout << " " << g.id(n) << " " << row.size() << std::endl;

   254 	lp.setRowCoeffs(row_it, row.begin(), row.end());

   255 	lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);

   256       }

   257       lp.solveSimplex();

   258       //std::cout << lp.colNum() << std::endl;

   259       //std::cout << lp.rowNum() << std::endl;

   260       //std::cout << "flow value: "<< lp.getObjVal() << std::endl;

   261       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)

   262       flow.set(e, lp_flow[e]);

   263     }

   264   };

   266 } // namespace lemon

   268 #endif //LEMON_MIN_COST_GEN_FLOW_H

author	marci
	Thu, 17 Feb 2005 15:14:13 +0000
changeset 1153	4b0468de3a31
parent 1074	4a24a46407db
permissions	-rw-r--r--