1.1 --- a/lemon/ssp_min_cost_flow.h Thu Apr 17 21:28:21 2008 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,260 +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-2008
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_SSP_MIN_COST_FLOW_H
1.23 -#define LEMON_SSP_MIN_COST_FLOW_H
1.24 -
1.25 -///\ingroup min_cost_flow
1.26 -///
1.27 -///\file
1.28 -///\brief An algorithm for finding a flow of value \c k (for
1.29 -///small values of \c k) having minimal total cost
1.30 -
1.31 -
1.32 -#include <lemon/dijkstra.h>
1.33 -#include <lemon/graph_adaptor.h>
1.34 -#include <lemon/maps.h>
1.35 -#include <vector>
1.36 -
1.37 -namespace lemon {
1.38 -
1.39 - /// \addtogroup min_cost_flow
1.40 - /// @{
1.41 -
1.42 - /// \brief Implementation of an algorithm for finding a flow of
1.43 - /// value \c k (for small values of \c k) having minimal total cost
1.44 - /// between two nodes
1.45 - ///
1.46 - ///
1.47 - /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum
1.48 - /// Cost Flow" implements an algorithm for finding a flow of value
1.49 - /// \c k having minimal total cost from a given source node to a
1.50 - /// given target node in a directed graph with a cost function on
1.51 - /// the edges. To this end, the edge-capacities and edge-costs have
1.52 - /// to be nonnegative. The edge-capacities should be integers, but
1.53 - /// the edge-costs can be integers, reals or of other comparable
1.54 - /// numeric type. This algorithm is intended to be used only for
1.55 - /// small values of \c k, since it is only polynomial in k, not in
1.56 - /// the length of k (which is log k): in order to find the minimum
1.57 - /// cost flow of value \c k it finds the minimum cost flow of value
1.58 - /// \c i for every \c i between 0 and \c k.
1.59 - ///
1.60 - ///\param Graph The directed graph type the algorithm runs on.
1.61 - ///\param LengthMap The type of the length map.
1.62 - ///\param CapacityMap The capacity map type.
1.63 - ///
1.64 - ///\author Attila Bernath
1.65 - template <typename Graph, typename LengthMap, typename CapacityMap>
1.66 - class SspMinCostFlow {
1.67 -
1.68 - typedef typename LengthMap::Value Length;
1.69 -
1.70 - //Warning: this should be integer type
1.71 - typedef typename CapacityMap::Value Capacity;
1.72 -
1.73 - typedef typename Graph::Node Node;
1.74 - typedef typename Graph::NodeIt NodeIt;
1.75 - typedef typename Graph::Edge Edge;
1.76 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.77 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.78 -
1.79 - typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
1.80 - typedef typename ResGW::Edge ResGraphEdge;
1.81 -
1.82 - protected:
1.83 -
1.84 - const Graph& g;
1.85 - const LengthMap& length;
1.86 - const CapacityMap& capacity;
1.87 -
1.88 - EdgeIntMap flow;
1.89 - typedef typename Graph::template NodeMap<Length> PotentialMap;
1.90 - PotentialMap potential;
1.91 -
1.92 - Node s;
1.93 - Node t;
1.94 -
1.95 - Length total_length;
1.96 -
1.97 - class ModLengthMap {
1.98 - typedef typename Graph::template NodeMap<Length> NodeMap;
1.99 - const ResGW& g;
1.100 - const LengthMap &length;
1.101 - const NodeMap &pot;
1.102 - public :
1.103 - typedef typename LengthMap::Key Key;
1.104 - typedef typename LengthMap::Value Value;
1.105 -
1.106 - ModLengthMap(const ResGW& _g,
1.107 - const LengthMap &_length, const NodeMap &_pot) :
1.108 - g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
1.109 -
1.110 - Value operator[](typename ResGW::Edge e) const {
1.111 - if (g.forward(e))
1.112 - return length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.113 - else
1.114 - return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);
1.115 - }
1.116 -
1.117 - }; //ModLengthMap
1.118 -
1.119 - ResGW res_graph;
1.120 - ModLengthMap mod_length;
1.121 - Dijkstra<ResGW, ModLengthMap> dijkstra;
1.122 -
1.123 - public :
1.124 -
1.125 - /// \brief The constructor of the class.
1.126 - ///
1.127 - /// \param _g The directed graph the algorithm runs on.
1.128 - /// \param _length The length (cost) of the edges.
1.129 - /// \param _cap The capacity of the edges.
1.130 - /// \param _s Source node.
1.131 - /// \param _t Target node.
1.132 - SspMinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap,
1.133 - Node _s, Node _t) :
1.134 - g(_g), length(_length), capacity(_cap), flow(_g), potential(_g),
1.135 - s(_s), t(_t),
1.136 - res_graph(g, capacity, flow),
1.137 - mod_length(res_graph, length, potential),
1.138 - dijkstra(res_graph, mod_length) {
1.139 - reset();
1.140 - }
1.141 -
1.142 - /// \brief Tries to augment the flow between s and t by 1. The
1.143 - /// return value shows if the augmentation is successful.
1.144 - bool augment() {
1.145 - dijkstra.run(s);
1.146 - if (!dijkstra.reached(t)) {
1.147 -
1.148 - //Unsuccessful augmentation.
1.149 - return false;
1.150 - } else {
1.151 -
1.152 - //We have to change the potential
1.153 - for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
1.154 - potential.set(n, potential[n]+dijkstra.distMap()[n]);
1.155 -
1.156 - //Augmenting on the shortest path
1.157 - Node n=t;
1.158 - ResGraphEdge e;
1.159 - while (n!=s){
1.160 - e = dijkstra.predEdge(n);
1.161 - n = dijkstra.predNode(n);
1.162 - res_graph.augment(e,1);
1.163 - //Let's update the total length
1.164 - if (res_graph.forward(e))
1.165 - total_length += length[e];
1.166 - else
1.167 - total_length -= length[e];
1.168 - }
1.169 -
1.170 - return true;
1.171 - }
1.172 - }
1.173 -
1.174 - /// \brief Runs the algorithm.
1.175 - ///
1.176 - /// Runs the algorithm.
1.177 - /// Returns k if there is a flow of value at least k from s to t.
1.178 - /// Otherwise it returns the maximum value of a flow from s to t.
1.179 - ///
1.180 - /// \param k The value of the flow we are looking for.
1.181 - ///
1.182 - /// \todo May be it does make sense to be able to start with a
1.183 - /// nonzero feasible primal-dual solution pair as well.
1.184 - ///
1.185 - /// \todo If the current flow value is bigger than k, then
1.186 - /// everything is cleared and the algorithm starts from zero
1.187 - /// flow. Is it a good approach?
1.188 - int run(int k) {
1.189 - if (flowValue()>k) reset();
1.190 - while (flowValue()<k && augment()) { }
1.191 - return flowValue();
1.192 - }
1.193 -
1.194 - /// \brief The 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 current 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_SSP_MIN_COST_FLOW_H