COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/hugo/min_cost_flow.h @ 906:17f31d280385

Last change on this file since 906:17f31d280385 was 906:17f31d280385, checked in by Alpar Juttner, 20 years ago

Copyright header added.

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