src/work/marci/lp/min_cost_gen_flow.h
changeset 1030 c8a41699e613
parent 1025 3b1ad8bc21da
child 1031 0b7169db694f
equal deleted inserted replaced
1:a62585338e07 2:97b3695fd746
     8 #include <lemon/list_graph.h>
     8 #include <lemon/list_graph.h>
     9 //#include <lemon/dimacs.h>
     9 //#include <lemon/dimacs.h>
    10 //#include <lemon/time_measure.h>
    10 //#include <lemon/time_measure.h>
    11 //#include <graph_wrapper.h>
    11 //#include <graph_wrapper.h>
    12 #include <lemon/preflow.h>
    12 #include <lemon/preflow.h>
       
    13 #include <lemon/min_cost_flow.h>
    13 //#include <augmenting_flow.h>
    14 //#include <augmenting_flow.h>
    14 //#include <preflow_res.h>
    15 //#include <preflow_res.h>
    15 #include <../merge_node_graph_wrapper.h>
    16 #include <work/marci/merge_node_graph_wrapper.h>
    16 #include <lemon/../work/marci/lp/lp_solver_wrapper.h>
    17 #include <work/marci/lp/lp_solver_wrapper.h>
    17 
    18 
    18 namespace lemon {
    19 namespace lemon {
    19 
    20 
    20   template<typename Edge, typename EdgeIndexMap> 
    21   template<typename Edge, typename EdgeIndexMap> 
    21   class PrimalMap {
    22   class PrimalMap {
    28     double operator[](Edge e) const { 
    29     double operator[](Edge e) const { 
    29       return lp->getPrimal((*edge_index_map)[e]);
    30       return lp->getPrimal((*edge_index_map)[e]);
    30     }
    31     }
    31   };
    32   };
    32 
    33 
    33   // excess: rho-delta
    34   // excess: rho-delta egyelore csak =0-ra.
    34   template <typename Graph, typename Num,
    35   template <typename Graph, typename Num,
    35 	    typename Excess=typename Graph::template NodeMap<Num>, 
    36 	    typename Excess=typename Graph::template NodeMap<Num>, 
    36 	    typename LCapMap=typename Graph::template EdgeMap<Num>,
    37 	    typename LCapMap=typename Graph::template EdgeMap<Num>,
    37 	    typename CapMap=typename Graph::template EdgeMap<Num>,
    38 	    typename CapMap=typename Graph::template EdgeMap<Num>,
    38             typename FlowMap=typename Graph::template EdgeMap<Num>,
    39             typename FlowMap=typename Graph::template EdgeMap<Num>,
    72       GWWW gwww(gw, gww);
    73       GWWW gwww(gw, gww);
    73 
    74 
    74       //      std::cout << "making new edges..." << std::endl; 
    75       //      std::cout << "making new edges..." << std::endl; 
    75       typename GWWW::template EdgeMap<Num> translated_cap(gwww);
    76       typename GWWW::template EdgeMap<Num> translated_cap(gwww);
    76 
    77 
    77       for (typename GW::EdgeIt e(gw); e!=INVALID; ++e)
    78       for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
    78       translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
    79 	translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
    79 			 capacity[e]-lcapacity[e]);
    80 			   capacity[e]-lcapacity[e]);
       
    81 	//	cout << "t_cap " << gw.id(e) << " " 
       
    82 	//	     << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
       
    83       }
    80 
    84 
    81       Num expected=0;
    85       Num expected=0;
    82 
    86 
    83       //      std::cout << "making new edges 2..." << std::endl; 
    87       //      std::cout << "making new edges 2..." << std::endl; 
    84       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
    88       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
    90 	if (excess[n]>a) {
    94 	if (excess[n]>a) {
    91 	  typename GWW::Edge e=
    95 	  typename GWW::Edge e=
    92 	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
    96 	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
    93 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
    97 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
    94 			     excess[n]-a);
    98 			     excess[n]-a);
       
    99 	  //	  std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
    95 	}
   100 	}
    96 	if (excess[n]<a) {
   101 	if (excess[n]<a) {
    97 	  typename GWW::Edge e=
   102 	  typename GWW::Edge e=
    98 	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
   103 	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
    99 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
   104 	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
   100 			     a-excess[n]);
   105 			     a-excess[n]);
   101 	  expected+=a-excess[n];
   106 	  expected+=a-excess[n];
       
   107 	  //	  std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
   102 	}
   108 	}
   103       }
   109       }
   104 
   110 
   105       //      std::cout << "preflow..." << std::endl; 
   111       //      std::cout << "preflow..." << std::endl; 
   106       typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
   112       typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
   115 	typename GWWW::Edge ewww(ew, INVALID, false);
   121 	typename GWWW::Edge ewww(ew, INVALID, false);
   116 	flow.set(e, translated_flow[ewww]+lcapacity[e]);
   122 	flow.set(e, translated_flow[ewww]+lcapacity[e]);
   117       }
   123       }
   118       return (expected>=preflow.flowValue());
   124       return (expected>=preflow.flowValue());
   119     }
   125     }
   120     void run() {
   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() {
   121       LPSolverWrapper lp;
   214       LPSolverWrapper lp;
   122       lp.setMinimize();
   215       lp.setMinimize();
   123       typedef LPSolverWrapper::ColIt ColIt;
   216       typedef LPSolverWrapper::ColIt ColIt;
   124       typedef LPSolverWrapper::RowIt RowIt;
   217       typedef LPSolverWrapper::RowIt RowIt;
   125       typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
   218       typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;