COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/lp/min_cost_gen_flow.h @ 1027:4ec35d1cd897

Last change on this file since 1027:4ec35d1cd897 was 1025:3b1ad8bc21da, checked in by marci, 19 years ago

MergeGraphWrapper? bug fixes

File size: 5.7 KB
Line 
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 <augmenting_flow.h>
14//#include <preflow_res.h>
15#include <../merge_node_graph_wrapper.h>
16#include <lemon/../work/marci/lp/lp_solver_wrapper.h>
17
18namespace lemon {
19
20  template<typename Edge, typename EdgeIndexMap>
21  class PrimalMap {
22  protected:
23    LPSolverWrapper* lp;
24    EdgeIndexMap* edge_index_map;
25  public:
26    PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) :
27      lp(&_lp), edge_index_map(&_edge_index_map) { }
28    double operator[](Edge e) const {
29      return lp->getPrimal((*edge_index_map)[e]);
30    }
31  };
32
33  // excess: rho-delta
34  template <typename Graph, typename Num,
35            typename Excess=typename Graph::template NodeMap<Num>,
36            typename LCapMap=typename Graph::template EdgeMap<Num>,
37            typename CapMap=typename Graph::template EdgeMap<Num>,
38            typename FlowMap=typename Graph::template EdgeMap<Num>,
39            typename CostMap=typename Graph::template EdgeMap<Num> >
40  class MinCostGenFlow {
41  protected:
42    const Graph& g;
43    const Excess& excess;
44    const LCapMap& lcapacity;
45    const CapMap& capacity;
46    FlowMap& flow;
47    const CostMap& cost;
48  public:
49    MinCostGenFlow(const Graph& _g, const Excess& _excess,
50                   const LCapMap& _lcapacity, const CapMap& _capacity,
51                   FlowMap& _flow,
52                   const CostMap& _cost) :
53      g(_g), excess(_excess), lcapacity(_lcapacity),
54      capacity(_capacity), flow(_flow), cost(_cost) { }
55    bool feasible() {
56      //      std::cout << "making new vertices..." << std::endl;
57      typedef ListGraph Graph2;
58      Graph2 g2;
59      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
60      //      std::cout << "merging..." << std::endl;
61      GW gw(g, g2);
62      typename GW::Node s(INVALID, g2.addNode(), true);
63      typename GW::Node t(INVALID, g2.addNode(), true);
64      typedef SmartGraph Graph3;
65      //      std::cout << "making extender graph..." << std::endl;
66      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
67//       {
68//      checkConcept<StaticGraph, GWW>();   
69//       }
70      GWW gww(gw);
71      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
72      GWWW gwww(gw, gww);
73
74      //      std::cout << "making new edges..." << std::endl;
75      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
76
77      for (typename GW::EdgeIt e(gw); e!=INVALID; ++e)
78      translated_cap.set(typename GWWW::Edge(e,INVALID,false),
79                         capacity[e]-lcapacity[e]);
80
81      Num expected=0;
82
83      //      std::cout << "making new edges 2..." << std::endl;
84      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
85        Num a=0;
86        for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
87          a+=lcapacity[e];
88        for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
89          a-=lcapacity[e];
90        if (excess[n]>a) {
91          typename GWW::Edge e=
92            gww.addEdge(typename GW::Node(n,INVALID,false), t);
93          translated_cap.set(typename GWWW::Edge(INVALID, e, true),
94                             excess[n]-a);
95        }
96        if (excess[n]<a) {
97          typename GWW::Edge e=
98            gww.addEdge(s, typename GW::Node(n,INVALID,false));
99          translated_cap.set(typename GWWW::Edge(INVALID, e, true),
100                             a-excess[n]);
101          expected+=a-excess[n];
102        }
103      }
104
105      //      std::cout << "preflow..." << std::endl;
106      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
107      Preflow<GWWW, Num> preflow(gwww, s, t,
108                                 translated_cap, translated_flow);
109      preflow.run();
110      //      std::cout << "fv: " << preflow.flowValue() << std::endl;
111      //      std::cout << "expected: " << expected << std::endl;
112
113      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
114        typename GW::Edge ew(e, INVALID, false);
115        typename GWWW::Edge ewww(ew, INVALID, false);
116        flow.set(e, translated_flow[ewww]+lcapacity[e]);
117      }
118      return (expected>=preflow.flowValue());
119    }
120    void run() {
121      LPSolverWrapper lp;
122      lp.setMinimize();
123      typedef LPSolverWrapper::ColIt ColIt;
124      typedef LPSolverWrapper::RowIt RowIt;
125      typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
126      EdgeIndexMap edge_index_map(g);
127      PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
128      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
129        ColIt col_it=lp.addCol();
130        edge_index_map.set(e, col_it);
131        if (lcapacity[e]==capacity[e])
132          lp.setColBounds(col_it, LPX_FX, lcapacity[e], capacity[e]);
133        else
134          lp.setColBounds(col_it, LPX_DB, lcapacity[e], capacity[e]);
135        lp.setObjCoef(col_it, cost[e]);
136      }
137      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
138        typename Graph::template EdgeMap<Num> coeffs(g, 0);
139        for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
140        coeffs.set(e, coeffs[e]+1);
141        for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
142        coeffs.set(e, coeffs[e]-1);
143        RowIt row_it=lp.addRow();
144        typename std::vector< std::pair<ColIt, double> > row;
145        //std::cout << "node:" <<g.id(n)<<std::endl;
146        for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
147          if (coeffs[e]!=0) {
148            //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
149            row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
150          }
151        }
152        //std::cout << std::endl;
153        lp.setRowCoeffs(row_it, row.begin(), row.end());
154        lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
155      }
156      lp.solveSimplex();
157      //std::cout << lp.colNum() << std::endl;
158      //std::cout << lp.rowNum() << std::endl;
159      //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
160      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e)
161      flow.set(e, lp_flow[e]);
162    }
163  };
164
165} // namespace lemon
166
167#endif //LEMON_MIN_COST_GEN_FLOW_H
Note: See TracBrowser for help on using the repository browser.