# source:lemon-0.x/lemon/ssp_min_cost_flow.h@2418:89cbf0a2ed57

Last change on this file since 2418:89cbf0a2ed57 was 2418:89cbf0a2ed57, checked in by athos, 13 years ago

Slight modifications.

File size: 7.9 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_SSP_MIN_COST_FLOW_H
20#define LEMON_SSP_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 current 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.
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 current 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_SSP_MIN_COST_FLOW_H
Note: See TracBrowser for help on using the repository browser.