COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/ssp_min_cost_flow.h @ 2391:14a343be7a5a

Last change on this file since 2391:14a343be7a5a was 2391:14a343be7a5a, checked in by Alpar Juttner, 12 years ago

Happy New Year to all source files!

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