1.1 --- a/src/work/marci/lp/min_cost_gen_flow.h Sun Apr 17 18:57:22 2005 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,268 +0,0 @@
1.4 -// -*- c++ -*-
1.5 -#ifndef LEMON_MIN_COST_GEN_FLOW_H
1.6 -#define LEMON_MIN_COST_GEN_FLOW_H
1.7 -#include <iostream>
1.8 -//#include <fstream>
1.9 -
1.10 -#include <lemon/smart_graph.h>
1.11 -#include <lemon/list_graph.h>
1.12 -//#include <lemon/dimacs.h>
1.13 -//#include <lemon/time_measure.h>
1.14 -//#include <graph_wrapper.h>
1.15 -#include <lemon/preflow.h>
1.16 -#include <lemon/min_cost_flow.h>
1.17 -//#include <augmenting_flow.h>
1.18 -//#include <preflow_res.h>
1.19 -#include <work/marci/merge_node_graph_wrapper.h>
1.20 -#include <work/marci/lp/lp_solver_wrapper_3.h>
1.21 -
1.22 -namespace lemon {
1.23 -
1.24 - template<typename Edge, typename EdgeIndexMap>
1.25 - class PrimalMap {
1.26 - protected:
1.27 - LPGLPK* lp;
1.28 - EdgeIndexMap* edge_index_map;
1.29 - public:
1.30 - PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) :
1.31 - lp(&_lp), edge_index_map(&_edge_index_map) { }
1.32 - double operator[](Edge e) const {
1.33 - return lp->getPrimal((*edge_index_map)[e]);
1.34 - }
1.35 - };
1.36 -
1.37 - // excess: rho-delta egyelore csak =0-ra.
1.38 - template <typename Graph, typename Num,
1.39 - typename Excess=typename Graph::template NodeMap<Num>,
1.40 - typename LCapMap=typename Graph::template EdgeMap<Num>,
1.41 - typename CapMap=typename Graph::template EdgeMap<Num>,
1.42 - typename FlowMap=typename Graph::template EdgeMap<Num>,
1.43 - typename CostMap=typename Graph::template EdgeMap<Num> >
1.44 - class MinCostGenFlow {
1.45 - protected:
1.46 - const Graph& g;
1.47 - const Excess& excess;
1.48 - const LCapMap& lcapacity;
1.49 - const CapMap& capacity;
1.50 - FlowMap& flow;
1.51 - const CostMap& cost;
1.52 - public:
1.53 - MinCostGenFlow(const Graph& _g, const Excess& _excess,
1.54 - const LCapMap& _lcapacity, const CapMap& _capacity,
1.55 - FlowMap& _flow,
1.56 - const CostMap& _cost) :
1.57 - g(_g), excess(_excess), lcapacity(_lcapacity),
1.58 - capacity(_capacity), flow(_flow), cost(_cost) { }
1.59 - bool feasible() {
1.60 - // std::cout << "making new vertices..." << std::endl;
1.61 - typedef ListGraph Graph2;
1.62 - Graph2 g2;
1.63 - typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
1.64 - // std::cout << "merging..." << std::endl;
1.65 - GW gw(g, g2);
1.66 - typename GW::Node s(INVALID, g2.addNode(), true);
1.67 - typename GW::Node t(INVALID, g2.addNode(), true);
1.68 - typedef SmartGraph Graph3;
1.69 - // std::cout << "making extender graph..." << std::endl;
1.70 - typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
1.71 -// {
1.72 -// checkConcept<StaticGraph, GWW>();
1.73 -// }
1.74 - GWW gww(gw);
1.75 - typedef AugmentingGraphWrapper<GW, GWW> GWWW;
1.76 - GWWW gwww(gw, gww);
1.77 -
1.78 - // std::cout << "making new edges..." << std::endl;
1.79 - typename GWWW::template EdgeMap<Num> translated_cap(gwww);
1.80 -
1.81 - for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
1.82 - translated_cap.set(typename GWWW::Edge(e,INVALID,false),
1.83 - capacity[e]-lcapacity[e]);
1.84 - // cout << "t_cap " << gw.id(e) << " "
1.85 - // << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
1.86 - }
1.87 -
1.88 - Num expected=0;
1.89 -
1.90 - // std::cout << "making new edges 2..." << std::endl;
1.91 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
1.92 - Num a=0;
1.93 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
1.94 - a+=lcapacity[e];
1.95 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
1.96 - a-=lcapacity[e];
1.97 - if (excess[n]>a) {
1.98 - typename GWW::Edge e=
1.99 - gww.addEdge(typename GW::Node(n,INVALID,false), t);
1.100 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
1.101 - excess[n]-a);
1.102 - // std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
1.103 - }
1.104 - if (excess[n]<a) {
1.105 - typename GWW::Edge e=
1.106 - gww.addEdge(s, typename GW::Node(n,INVALID,false));
1.107 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
1.108 - a-excess[n]);
1.109 - expected+=a-excess[n];
1.110 - // std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
1.111 - }
1.112 - }
1.113 -
1.114 - // std::cout << "preflow..." << std::endl;
1.115 - typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
1.116 - Preflow<GWWW, Num> preflow(gwww, s, t,
1.117 - translated_cap, translated_flow);
1.118 - preflow.run();
1.119 - // std::cout << "fv: " << preflow.flowValue() << std::endl;
1.120 - // std::cout << "expected: " << expected << std::endl;
1.121 -
1.122 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.123 - typename GW::Edge ew(e, INVALID, false);
1.124 - typename GWWW::Edge ewww(ew, INVALID, false);
1.125 - flow.set(e, translated_flow[ewww]+lcapacity[e]);
1.126 - }
1.127 - return (preflow.flowValue()>=expected);
1.128 - }
1.129 - // for nonnegative costs
1.130 - bool run() {
1.131 - // std::cout << "making new vertices..." << std::endl;
1.132 - typedef ListGraph Graph2;
1.133 - Graph2 g2;
1.134 - typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
1.135 - // std::cout << "merging..." << std::endl;
1.136 - GW gw(g, g2);
1.137 - typename GW::Node s(INVALID, g2.addNode(), true);
1.138 - typename GW::Node t(INVALID, g2.addNode(), true);
1.139 - typedef SmartGraph Graph3;
1.140 - // std::cout << "making extender graph..." << std::endl;
1.141 - typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
1.142 -// {
1.143 -// checkConcept<StaticGraph, GWW>();
1.144 -// }
1.145 - GWW gww(gw);
1.146 - typedef AugmentingGraphWrapper<GW, GWW> GWWW;
1.147 - GWWW gwww(gw, gww);
1.148 -
1.149 - // std::cout << "making new edges..." << std::endl;
1.150 - typename GWWW::template EdgeMap<Num> translated_cap(gwww);
1.151 -
1.152 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.153 - typename GW::Edge ew(e, INVALID, false);
1.154 - typename GWWW::Edge ewww(ew, INVALID, false);
1.155 - translated_cap.set(ewww, capacity[e]-lcapacity[e]);
1.156 - // cout << "t_cap " << g.id(e) << " "
1.157 - // << translated_cap[ewww] << endl;
1.158 - }
1.159 -
1.160 - Num expected=0;
1.161 -
1.162 - // std::cout << "making new edges 2..." << std::endl;
1.163 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
1.164 - // std::cout << "node: " << g.id(n) << std::endl;
1.165 - Num a=0;
1.166 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
1.167 - a+=lcapacity[e];
1.168 - // std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
1.169 - }
1.170 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
1.171 - a-=lcapacity[e];
1.172 - // std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
1.173 - }
1.174 - // std::cout << "excess " << g.id(n) << ": " << a << std::endl;
1.175 - if (0>a) {
1.176 - typename GWW::Edge e=
1.177 - gww.addEdge(typename GW::Node(n,INVALID,false), t);
1.178 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
1.179 - -a);
1.180 - // std::cout << g.id(n) << "->t " << -a << std::endl;
1.181 - }
1.182 - if (0<a) {
1.183 - typename GWW::Edge e=
1.184 - gww.addEdge(s, typename GW::Node(n,INVALID,false));
1.185 - translated_cap.set(typename GWWW::Edge(INVALID, e, true),
1.186 - a);
1.187 - expected+=a;
1.188 - // std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
1.189 - }
1.190 - }
1.191 -
1.192 - // std::cout << "preflow..." << std::endl;
1.193 - typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
1.194 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.195 - translated_cost.set(typename GWWW::Edge(
1.196 - typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
1.197 - }
1.198 - // typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
1.199 - MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>,
1.200 - typename GWWW::template EdgeMap<Num> >
1.201 - min_cost_flow(gwww, translated_cost, translated_cap,
1.202 - s, t);
1.203 - while (min_cost_flow.augment()) { }
1.204 - std::cout << "fv: " << min_cost_flow.flowValue() << std::endl;
1.205 - std::cout << "expected: " << expected << std::endl;
1.206 -
1.207 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.208 - typename GW::Edge ew(e, INVALID, false);
1.209 - typename GWWW::Edge ewww(ew, INVALID, false);
1.210 - // std::cout << g.id(e) << " " << flow[e] << std::endl;
1.211 - flow.set(e, lcapacity[e]+
1.212 - min_cost_flow.getFlow()[ewww]);
1.213 - }
1.214 - return (min_cost_flow.flowValue()>=expected);
1.215 - }
1.216 - void runByLP() {
1.217 - typedef LPGLPK LPSolver;
1.218 - LPSolver lp;
1.219 - lp.setMinimize();
1.220 - typedef LPSolver::ColIt ColIt;
1.221 - typedef LPSolver::RowIt RowIt;
1.222 - typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
1.223 - EdgeIndexMap edge_index_map(g);
1.224 - PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
1.225 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.226 - ColIt col_it=lp.addCol();
1.227 - edge_index_map.set(e, col_it);
1.228 - if (lcapacity[e]==capacity[e])
1.229 - lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
1.230 - else
1.231 - lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
1.232 - lp.setObjCoef(col_it, cost[e]);
1.233 - }
1.234 - LPSolver::ColIt col_it;
1.235 - for (lp.col_iter_map.first(col_it, lp.VALID_CLASS);
1.236 - lp.col_iter_map.valid(col_it);
1.237 - lp.col_iter_map.next(col_it)) {
1.238 -// std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
1.239 - }
1.240 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
1.241 - typename Graph::template EdgeMap<Num> coeffs(g, 0);
1.242 - for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
1.243 - coeffs.set(e, coeffs[e]+1);
1.244 - for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
1.245 - coeffs.set(e, coeffs[e]-1);
1.246 - RowIt row_it=lp.addRow();
1.247 - typename std::vector< std::pair<ColIt, double> > row;
1.248 - //std::cout << "node:" <<g.id(n)<<std::endl;
1.249 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.250 - if (coeffs[e]!=0) {
1.251 - //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
1.252 - row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
1.253 - }
1.254 - }
1.255 - //std::cout << std::endl;
1.256 - //std::cout << " " << g.id(n) << " " << row.size() << std::endl;
1.257 - lp.setRowCoeffs(row_it, row.begin(), row.end());
1.258 - lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
1.259 - }
1.260 - lp.solveSimplex();
1.261 - //std::cout << lp.colNum() << std::endl;
1.262 - //std::cout << lp.rowNum() << std::endl;
1.263 - //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
1.264 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)
1.265 - flow.set(e, lp_flow[e]);
1.266 - }
1.267 - };
1.268 -
1.269 -} // namespace lemon
1.270 -
1.271 -#endif //LEMON_MIN_COST_GEN_FLOW_H