1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/lemon/min_cost_flow.h Wed Sep 29 15:30:04 2004 +0000
1.3 @@ -0,0 +1,256 @@
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) 2004 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.8 + * (Egervary Combinatorial Optimization Research Group, 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_wrapper.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
1.43 + /// an algorithm for finding a flow of value \c k
1.44 + /// having minimal total cost
1.45 + /// from a given source node to a given target node in an
1.46 + /// edge-weighted directed graph. To this end,
1.47 + /// the edge-capacities and edge-weitghs have to be nonnegative.
1.48 + /// The edge-capacities should be integers, but the edge-weights can be
1.49 + /// integers, reals or of other comparable numeric type.
1.50 + /// This algorithm is intended to use only for small values of \c k,
1.51 + /// since it is only polynomial in k,
1.52 + /// not in the length of k (which is log k).
1.53 + /// In order to find the minimum cost flow of value \c k it
1.54 + /// finds the minimum cost flow of value \c i for every
1.55 + /// \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::ValueType Length;
1.66 +
1.67 + //Warning: this should be integer type
1.68 + typedef typename CapacityMap::ValueType 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 +
1.77 + typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
1.78 + typedef typename ResGW::Edge ResGraphEdge;
1.79 +
1.80 + class ModLengthMap {
1.81 + typedef typename Graph::template NodeMap<Length> NodeMap;
1.82 + const ResGW& G;
1.83 + const LengthMap &ol;
1.84 + const NodeMap &pot;
1.85 + public :
1.86 + typedef typename LengthMap::KeyType KeyType;
1.87 + typedef typename LengthMap::ValueType ValueType;
1.88 +
1.89 + ValueType operator[](typename ResGW::Edge e) const {
1.90 + if (G.forward(e))
1.91 + return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
1.92 + else
1.93 + return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
1.94 + }
1.95 +
1.96 + ModLengthMap(const ResGW& _G,
1.97 + const LengthMap &o, const NodeMap &p) :
1.98 + G(_G), /*rev(_rev),*/ ol(o), pot(p){};
1.99 + };//ModLengthMap
1.100 +
1.101 +
1.102 + protected:
1.103 +
1.104 + //Input
1.105 + const Graph& G;
1.106 + const LengthMap& length;
1.107 + const CapacityMap& capacity;
1.108 +
1.109 +
1.110 + //auxiliary variables
1.111 +
1.112 + //To store the flow
1.113 + EdgeIntMap flow;
1.114 + //To store the potential (dual variables)
1.115 + typedef typename Graph::template NodeMap<Length> PotentialMap;
1.116 + PotentialMap potential;
1.117 +
1.118 +
1.119 + Length total_length;
1.120 +
1.121 +
1.122 + public :
1.123 +
1.124 + /// The constructor of the class.
1.125 +
1.126 + ///\param _G The directed graph the algorithm runs on.
1.127 + ///\param _length The length (weight or cost) of the edges.
1.128 + ///\param _cap The capacity of the edges.
1.129 + MinCostFlow(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G),
1.130 + length(_length), capacity(_cap), flow(_G), potential(_G){ }
1.131 +
1.132 +
1.133 + ///Runs the algorithm.
1.134 +
1.135 + ///Runs the algorithm.
1.136 + ///Returns k if there is a flow of value at least k edge-disjoint
1.137 + ///from s to t.
1.138 + ///Otherwise it returns the maximum value of a flow from s to t.
1.139 + ///
1.140 + ///\param s The source node.
1.141 + ///\param t The target node.
1.142 + ///\param k The value of the flow we are looking for.
1.143 + ///
1.144 + ///\todo May be it does make sense to be able to start with a nonzero
1.145 + /// feasible primal-dual solution pair as well.
1.146 + int run(Node s, Node t, int k) {
1.147 +
1.148 + //Resetting variables from previous runs
1.149 + total_length = 0;
1.150 +
1.151 + for (typename Graph::EdgeIt e(G); e!=INVALID; ++e) flow.set(e, 0);
1.152 +
1.153 + //Initialize the potential to zero
1.154 + for (typename Graph::NodeIt n(G); n!=INVALID; ++n) potential.set(n, 0);
1.155 +
1.156 +
1.157 + //We need a residual graph
1.158 + ResGW res_graph(G, capacity, flow);
1.159 +
1.160 +
1.161 + ModLengthMap mod_length(res_graph, length, potential);
1.162 +
1.163 + Dijkstra<ResGW, ModLengthMap> dijkstra(res_graph, mod_length);
1.164 +
1.165 + int i;
1.166 + for (i=0; i<k; ++i){
1.167 + dijkstra.run(s);
1.168 + if (!dijkstra.reached(t)){
1.169 + //There are no flow of value k from s to t
1.170 + break;
1.171 + };
1.172 +
1.173 + //We have to change the potential
1.174 + for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
1.175 + potential[n] += dijkstra.distMap()[n];
1.176 +
1.177 +
1.178 + //Augmenting on the sortest path
1.179 + Node n=t;
1.180 + ResGraphEdge e;
1.181 + while (n!=s){
1.182 + e = dijkstra.pred(n);
1.183 + n = dijkstra.predNode(n);
1.184 + res_graph.augment(e,1);
1.185 + //Let's update the total length
1.186 + if (res_graph.forward(e))
1.187 + total_length += length[e];
1.188 + else
1.189 + total_length -= length[e];
1.190 + }
1.191 +
1.192 +
1.193 + }
1.194 +
1.195 +
1.196 + return i;
1.197 + }
1.198 +
1.199 +
1.200 +
1.201 + /// Gives back the total weight of the found flow.
1.202 +
1.203 + ///This function gives back the total weight of the found flow.
1.204 + ///Assumes that \c run() has been run and nothing changed since then.
1.205 + Length totalLength(){
1.206 + return total_length;
1.207 + }
1.208 +
1.209 + ///Returns a const reference to the EdgeMap \c flow.
1.210 +
1.211 + ///Returns a const reference to the EdgeMap \c flow.
1.212 + ///\pre \ref run() must
1.213 + ///be called before using this function.
1.214 + const EdgeIntMap &getFlow() const { return flow;}
1.215 +
1.216 + ///Returns a const reference to the NodeMap \c potential (the dual solution).
1.217 +
1.218 + ///Returns a const reference to the NodeMap \c potential (the dual solution).
1.219 + /// \pre \ref run() must be called before using this function.
1.220 + const PotentialMap &getPotential() const { return potential;}
1.221 +
1.222 + /// Checking the complementary slackness optimality criteria
1.223 +
1.224 + ///This function checks, whether the given solution is optimal
1.225 + ///If executed after the call of \c run() then it should return with true.
1.226 + ///This function only checks optimality, doesn't bother with feasibility.
1.227 + ///It is meant for testing purposes.
1.228 + ///
1.229 + bool checkComplementarySlackness(){
1.230 + Length mod_pot;
1.231 + Length fl_e;
1.232 + for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) {
1.233 + //C^{\Pi}_{i,j}
1.234 + mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
1.235 + fl_e = flow[e];
1.236 + if (0<fl_e && fl_e<capacity[e]) {
1.237 + /// \todo better comparison is needed for real types, moreover,
1.238 + /// this comparison here is superfluous.
1.239 + if (mod_pot != 0)
1.240 + return false;
1.241 + }
1.242 + else {
1.243 + if (mod_pot > 0 && fl_e != 0)
1.244 + return false;
1.245 + if (mod_pot < 0 && fl_e != capacity[e])
1.246 + return false;
1.247 + }
1.248 + }
1.249 + return true;
1.250 + }
1.251 +
1.252 +
1.253 + }; //class MinCostFlow
1.254 +
1.255 + ///@}
1.256 +
1.257 +} //namespace lemon
1.258 +
1.259 +#endif //LEMON_MIN_COST_FLOW_H