1.1 --- a/src/lemon/min_cost_flow.h Sat May 21 21:04:57 2005 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,260 +0,0 @@
1.4 -/* -*- C++ -*-
1.5 - * src/lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library
1.6 - *
1.7 - * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.8 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.9 - *
1.10 - * Permission to use, modify and distribute this software is granted
1.11 - * provided that this copyright notice appears in all copies. For
1.12 - * precise terms see the accompanying LICENSE file.
1.13 - *
1.14 - * This software is provided "AS IS" with no warranty of any kind,
1.15 - * express or implied, and with no claim as to its suitability for any
1.16 - * purpose.
1.17 - *
1.18 - */
1.19 -
1.20 -#ifndef LEMON_MIN_COST_FLOW_H
1.21 -#define LEMON_MIN_COST_FLOW_H
1.22 -
1.23 -///\ingroup flowalgs
1.24 -///\file
1.25 -///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
1.26 -
1.27 -
1.28 -#include <lemon/dijkstra.h>
1.29 -#include <lemon/graph_adaptor.h>
1.30 -#include <lemon/maps.h>
1.31 -#include <vector>
1.32 -
1.33 -namespace lemon {
1.34 -
1.35 -/// \addtogroup flowalgs
1.36 -/// @{
1.37 -
1.38 - ///\brief Implementation of an algorithm for finding a flow of value \c k
1.39 - ///(for small values of \c k) having minimal total cost between 2 nodes
1.40 - ///
1.41 - ///
1.42 - /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an
1.43 - /// algorithm for finding a flow of value \c k having minimal total
1.44 - /// cost from a given source node to a given target node in an
1.45 - /// edge-weighted directed graph. To this end, the edge-capacities
1.46 - /// and edge-weights have to be nonnegative. The edge-capacities
1.47 - /// should be integers, but the edge-weights can be integers, reals
1.48 - /// or of other comparable numeric type. This algorithm is intended
1.49 - /// to be used only for small values of \c k, since it is only
1.50 - /// polynomial in k, not in the length of k (which is log k): in
1.51 - /// order to find the minimum cost flow of value \c k it finds the
1.52 - /// minimum cost flow of value \c i for every \c i between 0 and \c
1.53 - /// k.
1.54 - ///
1.55 - ///\param Graph The directed graph type the algorithm runs on.
1.56 - ///\param LengthMap The type of the length map.
1.57 - ///\param CapacityMap The capacity map type.
1.58 - ///
1.59 - ///\author Attila Bernath
1.60 - template <typename Graph, typename LengthMap, typename CapacityMap>
1.61 - class MinCostFlow {
1.62 -
1.63 - typedef typename LengthMap::Value Length;
1.64 -
1.65 - //Warning: this should be integer type
1.66 - typedef typename CapacityMap::Value Capacity;
1.67 -
1.68 - typedef typename Graph::Node Node;
1.69 - typedef typename Graph::NodeIt NodeIt;
1.70 - typedef typename Graph::Edge Edge;
1.71 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.72 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.73 -
1.74 - typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
1.75 - typedef typename ResGW::Edge ResGraphEdge;
1.76 -
1.77 - protected:
1.78 -
1.79 - const Graph& g;
1.80 - const LengthMap& length;
1.81 - const CapacityMap& capacity;
1.82 -
1.83 - EdgeIntMap flow;
1.84 - typedef typename Graph::template NodeMap<Length> PotentialMap;
1.85 - PotentialMap potential;
1.86 -
1.87 - Node s;
1.88 - Node t;
1.89 -
1.90 - Length total_length;
1.91 -
1.92 - class ModLengthMap {
1.93 - typedef typename Graph::template NodeMap<Length> NodeMap;
1.94 - const ResGW& g;
1.95 - const LengthMap &length;
1.96 - const NodeMap &pot;
1.97 - public :
1.98 - typedef typename LengthMap::Key Key;
1.99 - typedef typename LengthMap::Value Value;
1.100 -
1.101 - ModLengthMap(const ResGW& _g,
1.102 - const LengthMap &_length, const NodeMap &_pot) :
1.103 - g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
1.104 -
1.105 - Value operator[](typename ResGW::Edge e) const {
1.106 - if (g.forward(e))
1.107 - return length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.108 - else
1.109 - return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.110 - }
1.111 -
1.112 - }; //ModLengthMap
1.113 -
1.114 - ResGW res_graph;
1.115 - ModLengthMap mod_length;
1.116 - Dijkstra<ResGW, ModLengthMap> dijkstra;
1.117 -
1.118 - public :
1.119 -
1.120 - /*! \brief The constructor of the class.
1.121 -
1.122 - \param _g The directed graph the algorithm runs on.
1.123 - \param _length The length (weight or cost) of the edges.
1.124 - \param _cap The capacity of the edges.
1.125 - \param _s Source node.
1.126 - \param _t Target node.
1.127 - */
1.128 - MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
1.129 - Node _s, Node _t) :
1.130 - g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
1.131 - s(_s), t(_t),
1.132 - res_graph(g, capacity, flow),
1.133 - mod_length(res_graph, length, potential),
1.134 - dijkstra(res_graph, mod_length) {
1.135 - reset();
1.136 - }
1.137 -
1.138 - /*! Tries to augment the flow between s and t by 1.
1.139 - The return value shows if the augmentation is successful.
1.140 - */
1.141 - bool augment() {
1.142 - dijkstra.run(s);
1.143 - if (!dijkstra.reached(t)) {
1.144 -
1.145 - //Unsuccessful augmentation.
1.146 - return false;
1.147 - } else {
1.148 -
1.149 - //We have to change the potential
1.150 - for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
1.151 - potential.set(n, potential[n]+dijkstra.distMap()[n]);
1.152 -
1.153 - //Augmenting on the shortest path
1.154 - Node n=t;
1.155 - ResGraphEdge e;
1.156 - while (n!=s){
1.157 - e = dijkstra.pred(n);
1.158 - n = dijkstra.predNode(n);
1.159 - res_graph.augment(e,1);
1.160 - //Let's update the total length
1.161 - if (res_graph.forward(e))
1.162 - total_length += length[e];
1.163 - else
1.164 - total_length -= length[e];
1.165 - }
1.166 -
1.167 - return true;
1.168 - }
1.169 - }
1.170 -
1.171 - /*! \brief Runs the algorithm.
1.172 -
1.173 - Runs the algorithm.
1.174 - Returns k if there is a flow of value at least k from s to t.
1.175 - Otherwise it returns the maximum value of a flow from s to t.
1.176 -
1.177 - \param k The value of the flow we are looking for.
1.178 -
1.179 - \todo May be it does make sense to be able to start with a nonzero
1.180 - feasible primal-dual solution pair as well.
1.181 -
1.182 - \todo If the actual flow value is bigger than k, then everything is
1.183 - cleared and the algorithm starts from zero flow. Is it a good approach?
1.184 - */
1.185 - int run(int k) {
1.186 - if (flowValue()>k) reset();
1.187 - while (flowValue()<k && augment()) { }
1.188 - return flowValue();
1.189 - }
1.190 -
1.191 - /*! \brief The class is reset to zero flow and potential.
1.192 - The class is reset to zero flow and potential.
1.193 - */
1.194 - void reset() {
1.195 - total_length=0;
1.196 - for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
1.197 - for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);
1.198 - }
1.199 -
1.200 - /*! Returns the value of the actual flow.
1.201 - */
1.202 - int flowValue() const {
1.203 - int i=0;
1.204 - for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
1.205 - for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
1.206 - return i;
1.207 - }
1.208 -
1.209 - /// Total weight of the found flow.
1.210 -
1.211 - /// This function gives back the total weight of the found flow.
1.212 - Length totalLength(){
1.213 - return total_length;
1.214 - }
1.215 -
1.216 - ///Returns a const reference to the EdgeMap \c flow.
1.217 -
1.218 - ///Returns a const reference to the EdgeMap \c flow.
1.219 - const EdgeIntMap &getFlow() const { return flow;}
1.220 -
1.221 - /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).
1.222 -
1.223 - Returns a const reference to the NodeMap \c potential (the dual solution).
1.224 - */
1.225 - const PotentialMap &getPotential() const { return potential;}
1.226 -
1.227 - /*! \brief Checking the complementary slackness optimality criteria.
1.228 -
1.229 - This function checks, whether the given flow and potential
1.230 - satisfy the complementary slackness conditions (i.e. these are optimal).
1.231 - This function only checks optimality, doesn't bother with feasibility.
1.232 - For testing purpose.
1.233 - */
1.234 - bool checkComplementarySlackness(){
1.235 - Length mod_pot;
1.236 - Length fl_e;
1.237 - for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
1.238 - //C^{\Pi}_{i,j}
1.239 - mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
1.240 - fl_e = flow[e];
1.241 - if (0<fl_e && fl_e<capacity[e]) {
1.242 - /// \todo better comparison is needed for real types, moreover,
1.243 - /// this comparison here is superfluous.
1.244 - if (mod_pot != 0)
1.245 - return false;
1.246 - }
1.247 - else {
1.248 - if (mod_pot > 0 && fl_e != 0)
1.249 - return false;
1.250 - if (mod_pot < 0 && fl_e != capacity[e])
1.251 - return false;
1.252 - }
1.253 - }
1.254 - return true;
1.255 - }
1.256 -
1.257 - }; //class MinCostFlow
1.258 -
1.259 - ///@}
1.260 -
1.261 -} //namespace lemon
1.262 -
1.263 -#endif //LEMON_MIN_COST_FLOW_H