1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/ssp_min_cost_flow.h Mon Oct 30 17:22:14 2006 +0000
1.3 @@ -0,0 +1,260 @@
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 +///
1.27 +///\file \brief An algorithm for finding a flow of value \c k (for
1.28 +///small values of \c k) having minimal total cost
1.29 +
1.30 +
1.31 +#include <lemon/dijkstra.h>
1.32 +#include <lemon/graph_adaptor.h>
1.33 +#include <lemon/maps.h>
1.34 +#include <vector>
1.35 +
1.36 +namespace lemon {
1.37 +
1.38 + /// \addtogroup flowalgs
1.39 + /// @{
1.40 +
1.41 + /// \brief Implementation of an algorithm for finding a flow of
1.42 + /// value \c k (for small values of \c k) having minimal total cost
1.43 + /// between two nodes
1.44 + ///
1.45 + ///
1.46 + /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum
1.47 + /// Cost Flow" implements an algorithm for finding a flow of value
1.48 + /// \c k having minimal total cost from a given source node to a
1.49 + /// given target node in a directed graph with a cost function on
1.50 + /// the edges. To this end, the edge-capacities and edge-costs have
1.51 + /// to be nonnegative. The edge-capacities should be integers, but
1.52 + /// the edge-costs can be integers, reals or of other comparable
1.53 + /// numeric type. This algorithm is intended to be used only for
1.54 + /// small values of \c k, since it is only polynomial in k, not in
1.55 + /// the length of k (which is log k): in order to find the minimum
1.56 + /// cost flow of value \c k it finds the minimum cost flow of value
1.57 + /// \c i for every \c i between 0 and \c k.
1.58 + ///
1.59 + ///\param Graph The directed graph type the algorithm runs on.
1.60 + ///\param LengthMap The type of the length map.
1.61 + ///\param CapacityMap The capacity map type.
1.62 + ///
1.63 + ///\author Attila Bernath
1.64 + template <typename Graph, typename LengthMap, typename CapacityMap>
1.65 + class SspMinCostFlow {
1.66 +
1.67 + typedef typename LengthMap::Value Length;
1.68 +
1.69 + //Warning: this should be integer type
1.70 + typedef typename CapacityMap::Value Capacity;
1.71 +
1.72 + typedef typename Graph::Node Node;
1.73 + typedef typename Graph::NodeIt NodeIt;
1.74 + typedef typename Graph::Edge Edge;
1.75 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.76 + typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.77 +
1.78 + typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
1.79 + typedef typename ResGW::Edge ResGraphEdge;
1.80 +
1.81 + protected:
1.82 +
1.83 + const Graph& g;
1.84 + const LengthMap& length;
1.85 + const CapacityMap& capacity;
1.86 +
1.87 + EdgeIntMap flow;
1.88 + typedef typename Graph::template NodeMap<Length> PotentialMap;
1.89 + PotentialMap potential;
1.90 +
1.91 + Node s;
1.92 + Node t;
1.93 +
1.94 + Length total_length;
1.95 +
1.96 + class ModLengthMap {
1.97 + typedef typename Graph::template NodeMap<Length> NodeMap;
1.98 + const ResGW& g;
1.99 + const LengthMap &length;
1.100 + const NodeMap &pot;
1.101 + public :
1.102 + typedef typename LengthMap::Key Key;
1.103 + typedef typename LengthMap::Value Value;
1.104 +
1.105 + ModLengthMap(const ResGW& _g,
1.106 + const LengthMap &_length, const NodeMap &_pot) :
1.107 + g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
1.108 +
1.109 + Value operator[](typename ResGW::Edge e) const {
1.110 + if (g.forward(e))
1.111 + return length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.112 + else
1.113 + return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.114 + }
1.115 +
1.116 + }; //ModLengthMap
1.117 +
1.118 + ResGW res_graph;
1.119 + ModLengthMap mod_length;
1.120 + Dijkstra<ResGW, ModLengthMap> dijkstra;
1.121 +
1.122 + public :
1.123 +
1.124 + /// \brief The constructor of the class.
1.125 + ///
1.126 + /// \param _g The directed graph the algorithm runs on.
1.127 + /// \param _length The length (cost) of the edges.
1.128 + /// \param _cap The capacity of the edges.
1.129 + /// \param _s Source node.
1.130 + /// \param _t Target node.
1.131 + SspMinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
1.132 + Node _s, Node _t) :
1.133 + g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
1.134 + s(_s), t(_t),
1.135 + res_graph(g, capacity, flow),
1.136 + mod_length(res_graph, length, potential),
1.137 + dijkstra(res_graph, mod_length) {
1.138 + reset();
1.139 + }
1.140 +
1.141 + /// \brief Tries to augment the flow between s and t by 1. The
1.142 + /// return value shows if the augmentation is successful.
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
1.182 + /// nonzero feasible primal-dual solution pair as well.
1.183 + ///
1.184 + /// \todo If the actual flow value is bigger than k, then
1.185 + /// everything is cleared and the algorithm starts from zero
1.186 + /// flow. Is it a good approach?
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. The
1.194 + /// class is reset to zero flow and potential.
1.195 + void reset() {
1.196 + total_length=0;
1.197 + for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
1.198 + for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);
1.199 + }
1.200 +
1.201 + /// \brief Returns the value of the actual flow.
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 + /// \brief Total cost of the found flow.
1.210 + ///
1.211 + /// This function gives back the total cost of the found flow.
1.212 + Length totalLength(){
1.213 + return total_length;
1.214 + }
1.215 +
1.216 + /// \brief 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
1.222 + /// (the dual solution).
1.223 + ///
1.224 + /// Returns a const reference to the NodeMap \c potential (the
1.225 + /// dual solution).
1.226 + const PotentialMap &getPotential() const { return potential;}
1.227 +
1.228 + /// \brief Checking the complementary slackness optimality criteria.
1.229 + ///
1.230 + /// This function checks, whether the given flow and potential
1.231 + /// satisfy the complementary slackness conditions (i.e. these are optimal).
1.232 + /// This function only checks optimality, doesn't bother with feasibility.
1.233 + /// For testing purpose.
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 SspMinCostFlow
1.258 +
1.259 + ///@}
1.260 +
1.261 +} //namespace lemon
1.262 +
1.263 +#endif //LEMON_MIN_COST_FLOW_H