src/hugo/mincostflows.h
author deba
Wed, 08 Sep 2004 12:06:45 +0000
changeset 822 88226d9fe821
parent 785 a9b0863c2265
child 860 3577b3db6089
permissions -rw-r--r--
The MapFactories have been removed from the code because
if we use macros then they increases only the complexity.

The pair iterators of the maps are separeted from the maps.

Some macros and comments has been changed.
athos@610
     1
// -*- c++ -*-
athos@610
     2
#ifndef HUGO_MINCOSTFLOWS_H
athos@610
     3
#define HUGO_MINCOSTFLOWS_H
athos@610
     4
alpar@758
     5
///\ingroup flowalgs
athos@610
     6
///\file
athos@610
     7
///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost 
athos@610
     8
athos@611
     9
athos@610
    10
#include <hugo/dijkstra.h>
athos@610
    11
#include <hugo/graph_wrapper.h>
athos@610
    12
#include <hugo/maps.h>
athos@610
    13
#include <vector>
athos@610
    14
athos@610
    15
namespace hugo {
athos@610
    16
alpar@758
    17
/// \addtogroup flowalgs
athos@610
    18
/// @{
athos@610
    19
athos@610
    20
  ///\brief Implementation of an algorithm for finding a flow of value \c k 
athos@610
    21
  ///(for small values of \c k) having minimal total cost between 2 nodes 
athos@610
    22
  /// 
athos@610
    23
  ///
athos@610
    24
  /// The class \ref hugo::MinCostFlows "MinCostFlows" implements
athos@610
    25
  /// an algorithm for finding a flow of value \c k 
athos@610
    26
  ///(for small values of \c k) having minimal total cost  
athos@610
    27
  /// from a given source node to a given target node in an
athos@610
    28
  /// edge-weighted directed graph having nonnegative integer capacities.
athos@610
    29
  /// The range of the length (weight) function is nonnegative reals but 
athos@610
    30
  /// the range of capacity function is the set of nonnegative integers. 
athos@610
    31
  /// It is not a polinomial time algorithm for counting the minimum cost
athos@610
    32
  /// maximal flow, since it counts the minimum cost flow for every value 0..M
athos@610
    33
  /// where \c M is the value of the maximal flow.
athos@610
    34
  ///
athos@610
    35
  ///\author Attila Bernath
athos@610
    36
  template <typename Graph, typename LengthMap, typename CapacityMap>
athos@610
    37
  class MinCostFlows {
athos@610
    38
athos@610
    39
    typedef typename LengthMap::ValueType Length;
athos@610
    40
athos@610
    41
    //Warning: this should be integer type
athos@610
    42
    typedef typename CapacityMap::ValueType Capacity;
athos@610
    43
    
athos@610
    44
    typedef typename Graph::Node Node;
athos@610
    45
    typedef typename Graph::NodeIt NodeIt;
athos@610
    46
    typedef typename Graph::Edge Edge;
athos@610
    47
    typedef typename Graph::OutEdgeIt OutEdgeIt;
athos@610
    48
    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
athos@610
    49
athos@610
    50
    //    typedef ConstMap<Edge,int> ConstMap;
athos@610
    51
athos@610
    52
    typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
athos@610
    53
    typedef typename ResGraphType::Edge ResGraphEdge;
athos@610
    54
athos@610
    55
    class ModLengthMap {   
athos@610
    56
      //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
athos@610
    57
      typedef typename Graph::template NodeMap<Length> NodeMap;
athos@610
    58
      const ResGraphType& G;
athos@610
    59
      //      const EdgeIntMap& rev;
athos@610
    60
      const LengthMap &ol;
athos@610
    61
      const NodeMap &pot;
athos@610
    62
    public :
athos@610
    63
      typedef typename LengthMap::KeyType KeyType;
athos@610
    64
      typedef typename LengthMap::ValueType ValueType;
athos@610
    65
	
athos@610
    66
      ValueType operator[](typename ResGraphType::Edge e) const {     
athos@610
    67
	if (G.forward(e))
athos@610
    68
	  return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
athos@610
    69
	else
athos@610
    70
	  return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
athos@610
    71
      }     
athos@610
    72
	
athos@610
    73
      ModLengthMap(const ResGraphType& _G,
athos@610
    74
		   const LengthMap &o,  const NodeMap &p) : 
athos@610
    75
	G(_G), /*rev(_rev),*/ ol(o), pot(p){}; 
athos@610
    76
    };//ModLengthMap
athos@610
    77
athos@610
    78
athos@610
    79
  protected:
athos@610
    80
    
athos@610
    81
    //Input
athos@610
    82
    const Graph& G;
athos@610
    83
    const LengthMap& length;
athos@610
    84
    const CapacityMap& capacity;
athos@610
    85
athos@610
    86
athos@610
    87
    //auxiliary variables
athos@610
    88
athos@610
    89
    //To store the flow
athos@610
    90
    EdgeIntMap flow; 
alpar@785
    91
    //To store the potential (dual variables)
athos@661
    92
    typedef typename Graph::template NodeMap<Length> PotentialMap;
athos@661
    93
    PotentialMap potential;
athos@610
    94
    
athos@610
    95
athos@610
    96
    Length total_length;
athos@610
    97
athos@610
    98
athos@610
    99
  public :
athos@610
   100
athos@610
   101
athos@610
   102
    MinCostFlows(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G), 
athos@610
   103
      length(_length), capacity(_cap), flow(_G), potential(_G){ }
athos@610
   104
athos@610
   105
    
athos@610
   106
    ///Runs the algorithm.
athos@610
   107
athos@610
   108
    ///Runs the algorithm.
athos@610
   109
    ///Returns k if there are at least k edge-disjoint paths from s to t.
athos@610
   110
    ///Otherwise it returns the number of found edge-disjoint paths from s to t.
athos@610
   111
    ///\todo May be it does make sense to be able to start with a nonzero 
athos@610
   112
    /// feasible primal-dual solution pair as well.
athos@610
   113
    int run(Node s, Node t, int k) {
athos@610
   114
athos@610
   115
      //Resetting variables from previous runs
athos@610
   116
      total_length = 0;
athos@610
   117
      
marci@788
   118
      for (typename Graph::EdgeIt e(G); e!=INVALID; ++e) flow.set(e, 0);
athos@634
   119
athos@634
   120
      //Initialize the potential to zero
marci@788
   121
      for (typename Graph::NodeIt n(G); n!=INVALID; ++n) potential.set(n, 0);
athos@610
   122
      
athos@610
   123
      
athos@610
   124
      //We need a residual graph
athos@610
   125
      ResGraphType res_graph(G, capacity, flow);
athos@610
   126
athos@610
   127
athos@610
   128
      ModLengthMap mod_length(res_graph, length, potential);
athos@610
   129
athos@610
   130
      Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
athos@610
   131
athos@610
   132
      int i;
athos@610
   133
      for (i=0; i<k; ++i){
athos@610
   134
	dijkstra.run(s);
athos@610
   135
	if (!dijkstra.reached(t)){
athos@610
   136
	  //There are no k paths from s to t
athos@610
   137
	  break;
athos@610
   138
	};
athos@610
   139
	
athos@634
   140
	//We have to change the potential
marci@788
   141
        for(typename ResGraphType::NodeIt n(res_graph); n!=INVALID; ++n)
athos@633
   142
	  potential[n] += dijkstra.distMap()[n];
athos@634
   143
athos@610
   144
athos@610
   145
	//Augmenting on the sortest path
athos@610
   146
	Node n=t;
athos@610
   147
	ResGraphEdge e;
athos@610
   148
	while (n!=s){
athos@610
   149
	  e = dijkstra.pred(n);
athos@610
   150
	  n = dijkstra.predNode(n);
athos@610
   151
	  res_graph.augment(e,1);
athos@610
   152
	  //Let's update the total length
athos@610
   153
	  if (res_graph.forward(e))
athos@610
   154
	    total_length += length[e];
athos@610
   155
	  else 
athos@610
   156
	    total_length -= length[e];	    
athos@610
   157
	}
athos@610
   158
athos@610
   159
	  
athos@610
   160
      }
athos@610
   161
      
athos@610
   162
athos@610
   163
      return i;
athos@610
   164
    }
athos@610
   165
athos@610
   166
athos@610
   167
athos@610
   168
athos@610
   169
    ///This function gives back the total length of the found paths.
athos@610
   170
    ///Assumes that \c run() has been run and nothing changed since then.
athos@610
   171
    Length totalLength(){
athos@610
   172
      return total_length;
athos@610
   173
    }
athos@610
   174
athos@610
   175
    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
athos@610
   176
    ///be called before using this function.
athos@610
   177
    const EdgeIntMap &getFlow() const { return flow;}
athos@610
   178
athos@610
   179
  ///Returns a const reference to the NodeMap \c potential (the dual solution).
athos@610
   180
    /// \pre \ref run() must be called before using this function.
athos@661
   181
    const PotentialMap &getPotential() const { return potential;}
athos@610
   182
athos@610
   183
    ///This function checks, whether the given solution is optimal
athos@610
   184
    ///Running after a \c run() should return with true
athos@610
   185
    ///In this "state of the art" this only check optimality, doesn't bother with feasibility
athos@610
   186
    ///
athos@610
   187
    ///\todo Is this OK here?
athos@610
   188
    bool checkComplementarySlackness(){
athos@610
   189
      Length mod_pot;
athos@610
   190
      Length fl_e;
marci@788
   191
        for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) {
athos@610
   192
	//C^{\Pi}_{i,j}
athos@610
   193
	mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
athos@610
   194
	fl_e = flow[e];
athos@610
   195
	//	std::cout << fl_e << std::endl;
athos@610
   196
	if (0<fl_e && fl_e<capacity[e]){
athos@610
   197
	  if (mod_pot != 0)
athos@610
   198
	    return false;
athos@610
   199
	}
athos@610
   200
	else{
athos@610
   201
	  if (mod_pot > 0 && fl_e != 0)
athos@610
   202
	    return false;
athos@610
   203
	  if (mod_pot < 0 && fl_e != capacity[e])
athos@610
   204
	    return false;
athos@610
   205
	}
athos@610
   206
      }
athos@610
   207
      return true;
athos@610
   208
    }
athos@610
   209
    
athos@610
   210
athos@610
   211
  }; //class MinCostFlows
athos@610
   212
athos@610
   213
  ///@}
athos@610
   214
athos@610
   215
} //namespace hugo
athos@610
   216
athos@633
   217
#endif //HUGO_MINCOSTFLOWS_H