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