COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/ssp_min_cost_flow.h @ 2276:1a8a66b6c6ce

Last change on this file since 2276:1a8a66b6c6ce was 2276:1a8a66b6c6ce, checked in by Balazs Dezso, 13 years ago

Min cost flow is renamed to SspMinCostFlow?

File size: 8.0 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_MIN_COST_FLOW_H
20#define LEMON_MIN_COST_FLOW_H
21
22///\ingroup flowalgs
23///
24///\file \brief An algorithm for finding a flow of value \c k (for
25///small values of \c k) having minimal total cost
26
27
28#include <lemon/dijkstra.h>
29#include <lemon/graph_adaptor.h>
30#include <lemon/maps.h>
31#include <vector>
32
33namespace lemon {
34
35  /// \addtogroup flowalgs
36  /// @{
37
38  /// \brief Implementation of an algorithm for finding a flow of
39  /// value \c k (for small values of \c k) having minimal total cost
40  /// between two nodes
41  ///
42  ///
43  /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum
44  /// Cost Flow" implements an algorithm for finding a flow of value
45  /// \c k having minimal total cost from a given source node to a
46  /// given target node in a directed graph with a cost function on
47  /// the edges. To this end, the edge-capacities and edge-costs have
48  /// to be nonnegative.  The edge-capacities should be integers, but
49  /// the edge-costs can be integers, reals or of other comparable
50  /// numeric type.  This algorithm is intended to be used only for
51  /// small values of \c k, since it is only polynomial in k, not in
52  /// the length of k (which is log k): in order to find the minimum
53  /// cost flow of value \c k it finds the minimum cost flow of value
54  /// \c i for every \c i between 0 and \c k.
55  ///
56  ///\param Graph The directed graph type the algorithm runs on.
57  ///\param LengthMap The type of the length map.
58  ///\param CapacityMap The capacity map type.
59  ///
60  ///\author Attila Bernath
61  template <typename Graph, typename LengthMap, typename CapacityMap>
62  class SspMinCostFlow {
63
64    typedef typename LengthMap::Value Length;
65
66    //Warning: this should be integer type
67    typedef typename CapacityMap::Value Capacity;
68   
69    typedef typename Graph::Node Node;
70    typedef typename Graph::NodeIt NodeIt;
71    typedef typename Graph::Edge Edge;
72    typedef typename Graph::OutEdgeIt OutEdgeIt;
73    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
74
75    typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
76    typedef typename ResGW::Edge ResGraphEdge;
77
78  protected:
79
80    const Graph& g;
81    const LengthMap& length;
82    const CapacityMap& capacity;
83
84    EdgeIntMap flow;
85    typedef typename Graph::template NodeMap<Length> PotentialMap;
86    PotentialMap potential;
87
88    Node s;
89    Node t;
90   
91    Length total_length;
92
93    class ModLengthMap {   
94      typedef typename Graph::template NodeMap<Length> NodeMap;
95      const ResGW& g;
96      const LengthMap &length;
97      const NodeMap &pot;
98    public :
99      typedef typename LengthMap::Key Key;
100      typedef typename LengthMap::Value Value;
101
102      ModLengthMap(const ResGW& _g,
103                   const LengthMap &_length, const NodeMap &_pot) :
104        g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
105       
106      Value operator[](typename ResGW::Edge e) const {     
107        if (g.forward(e))
108          return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
109        else
110          return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
111      }     
112       
113    }; //ModLengthMap
114
115    ResGW res_graph;
116    ModLengthMap mod_length;
117    Dijkstra<ResGW, ModLengthMap> dijkstra;
118
119  public :
120
121    /// \brief The constructor of the class.
122    ///
123    /// \param _g The directed graph the algorithm runs on.
124    /// \param _length The length (cost) of the edges.
125    /// \param _cap The capacity of the edges.
126    /// \param _s Source node.
127    /// \param _t Target node.
128    SspMinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
129                   Node _s, Node _t) :
130      g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
131      s(_s), t(_t),
132      res_graph(g, capacity, flow),
133      mod_length(res_graph, length, potential),
134      dijkstra(res_graph, mod_length) {
135      reset();
136    }
137   
138    /// \brief Tries to augment the flow between s and t by 1. The
139    /// return value shows if the augmentation is successful.
140    bool augment() {
141      dijkstra.run(s);
142      if (!dijkstra.reached(t)) {
143
144        //Unsuccessful augmentation.
145        return false;
146      } else {
147
148        //We have to change the potential
149        for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
150          potential.set(n, potential[n]+dijkstra.distMap()[n]);
151       
152        //Augmenting on the shortest path
153        Node n=t;
154        ResGraphEdge e;
155        while (n!=s){
156          e = dijkstra.predEdge(n);
157          n = dijkstra.predNode(n);
158          res_graph.augment(e,1);
159          //Let's update the total length
160          if (res_graph.forward(e))
161            total_length += length[e];
162          else
163            total_length -= length[e];     
164        }
165
166        return true;
167      }
168    }
169   
170    /// \brief Runs the algorithm.
171    ///
172    /// Runs the algorithm.
173    /// Returns k if there is a flow of value at least k from s to t.
174    /// Otherwise it returns the maximum value of a flow from s to t.
175    ///
176    /// \param k The value of the flow we are looking for.
177    ///
178    /// \todo May be it does make sense to be able to start with a
179    /// nonzero feasible primal-dual solution pair as well.
180    ///
181    /// \todo If the actual flow value is bigger than k, then
182    /// everything is cleared and the algorithm starts from zero
183    /// flow. Is it a good approach?
184    int run(int k) {
185      if (flowValue()>k) reset();
186      while (flowValue()<k && augment()) { }
187      return flowValue();
188    }
189
190    /// \brief The class is reset to zero flow and potential. The
191    /// class is reset to zero flow and potential.
192    void reset() {
193      total_length=0;
194      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
195      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0); 
196    }
197
198    /// \brief Returns the value of the actual flow.
199    int flowValue() const {
200      int i=0;
201      for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
202      for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
203      return i;
204    }
205
206    /// \brief Total cost of the found flow.
207    ///
208    /// This function gives back the total cost of the found flow.
209    Length totalLength(){
210      return total_length;
211    }
212
213    /// \brief Returns a const reference to the EdgeMap \c flow.
214    ///
215    /// Returns a const reference to the EdgeMap \c flow.
216    const EdgeIntMap &getFlow() const { return flow;}
217
218    /// \brief Returns a const reference to the NodeMap \c potential
219    /// (the dual solution).
220    ///
221    /// Returns a const reference to the NodeMap \c potential (the
222    /// dual solution).
223    const PotentialMap &getPotential() const { return potential;}
224
225    /// \brief Checking the complementary slackness optimality criteria.
226    ///
227    /// This function checks, whether the given flow and potential
228    /// satisfy the complementary slackness conditions (i.e. these are optimal).
229    /// This function only checks optimality, doesn't bother with feasibility.
230    /// For testing purpose.
231    bool checkComplementarySlackness(){
232      Length mod_pot;
233      Length fl_e;
234        for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
235        //C^{\Pi}_{i,j}
236        mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
237        fl_e = flow[e];
238        if (0<fl_e && fl_e<capacity[e]) {
239          /// \todo better comparison is needed for real types, moreover,
240          /// this comparison here is superfluous.
241          if (mod_pot != 0)
242            return false;
243        }
244        else {
245          if (mod_pot > 0 && fl_e != 0)
246            return false;
247          if (mod_pot < 0 && fl_e != capacity[e])
248            return false;
249        }
250      }
251      return true;
252    }
253   
254  }; //class SspMinCostFlow
255
256  ///@}
257
258} //namespace lemon
259
260#endif //LEMON_MIN_COST_FLOW_H
Note: See TracBrowser for help on using the repository browser.