1.1 --- a/src/hugo/mincostflows.h Thu May 13 11:25:52 2004 +0000
1.2 +++ b/src/hugo/mincostflows.h Thu May 13 16:00:18 2004 +0000
1.3 @@ -92,11 +92,6 @@
1.4 //To store the potentila (dual variables)
1.5 typename Graph::template NodeMap<Length> potential;
1.6
1.7 - //Container to store found paths
1.8 - //std::vector< std::vector<Edge> > paths;
1.9 - //typedef DirPath<Graph> DPath;
1.10 - //DPath paths;
1.11 -
1.12
1.13 Length total_length;
1.14
1.15 @@ -151,6 +146,11 @@
1.16 break;
1.17 };
1.18
1.19 + //We have to copy the potential
1.20 + FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
1.21 + potential[n] += dijkstra.distMap()[n];
1.22 + }
1.23 + /*
1.24 {
1.25 //We have to copy the potential
1.26 typename ResGraphType::NodeIt n;
1.27 @@ -158,7 +158,7 @@
1.28 potential[n] += dijkstra.distMap()[n];
1.29 }
1.30 }
1.31 -
1.32 + */
1.33
1.34 //Augmenting on the sortest path
1.35 Node n=t;
1.36 @@ -225,25 +225,6 @@
1.37 return true;
1.38 }
1.39
1.40 - /*
1.41 - ///\todo To be implemented later
1.42 -
1.43 - ///This function gives back the \c j-th path in argument p.
1.44 - ///Assumes that \c run() has been run and nothing changed since then.
1.45 - /// \warning It is assumed that \c p is constructed to be a path of graph \c G. If \c j is greater than the result of previous \c run, then the result here will be an empty path.
1.46 - template<typename DirPath>
1.47 - void getPath(DirPath& p, int j){
1.48 - p.clear();
1.49 - typename DirPath::Builder B(p);
1.50 - for(typename std::vector<Edge>::iterator i=paths[j].begin();
1.51 - i!=paths[j].end(); ++i ){
1.52 - B.pushBack(*i);
1.53 - }
1.54 -
1.55 - B.commit();
1.56 - }
1.57 -
1.58 - */
1.59
1.60 }; //class MinCostFlows
1.61
1.62 @@ -251,4 +232,4 @@
1.63
1.64 } //namespace hugo
1.65
1.66 -#endif //HUGO_MINCOSTFLOW_H
1.67 +#endif //HUGO_MINCOSTFLOWS_H
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/work/athos/mincostflow.h Thu May 13 16:00:18 2004 +0000
2.3 @@ -0,0 +1,245 @@
2.4 +// -*- c++ -*-
2.5 +#ifndef HUGO_MINCOSTFLOW_H
2.6 +#define HUGO_MINCOSTFLOW_H
2.7 +
2.8 +///\ingroup galgs
2.9 +///\file
2.10 +///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
2.11 +
2.12 +
2.13 +#include <hugo/dijkstra.h>
2.14 +#include <hugo/graph_wrapper.h>
2.15 +#include <hugo/maps.h>
2.16 +#include <vector>
2.17 +#include <for_each_macros.h>
2.18 +
2.19 +namespace hugo {
2.20 +
2.21 +/// \addtogroup galgs
2.22 +/// @{
2.23 +
2.24 + ///\brief Implementation of an algorithm for finding a flow of value \c k
2.25 + ///(for small values of \c k) having minimal total cost between 2 nodes
2.26 + ///
2.27 + ///
2.28 + /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
2.29 + /// an algorithm for solving the following general minimum cost flow problem>
2.30 + ///
2.31 + ///
2.32 + ///
2.33 + /// \warning It is assumed here that the problem has a feasible solution
2.34 + ///
2.35 + /// The range of the length (weight) function is nonnegative reals but
2.36 + /// the range of capacity function is the set of nonnegative integers.
2.37 + /// It is not a polinomial time algorithm for counting the minimum cost
2.38 + /// maximal flow, since it counts the minimum cost flow for every value 0..M
2.39 + /// where \c M is the value of the maximal flow.
2.40 + ///
2.41 + ///\author Attila Bernath
2.42 + template <typename Graph, typename LengthMap, typename SupplyMap>
2.43 + class MinCostFlow {
2.44 +
2.45 + typedef typename LengthMap::ValueType Length;
2.46 +
2.47 +
2.48 + typedef typename SupplyMap::ValueType Supply;
2.49 +
2.50 + typedef typename Graph::Node Node;
2.51 + typedef typename Graph::NodeIt NodeIt;
2.52 + typedef typename Graph::Edge Edge;
2.53 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.54 + typedef typename Graph::template EdgeMap<int> EdgeIntMap;
2.55 +
2.56 + // typedef ConstMap<Edge,int> ConstMap;
2.57 +
2.58 + typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
2.59 + typedef typename ResGraphType::Edge ResGraphEdge;
2.60 +
2.61 + class ModLengthMap {
2.62 + //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
2.63 + typedef typename Graph::template NodeMap<Length> NodeMap;
2.64 + const ResGraphType& G;
2.65 + // const EdgeIntMap& rev;
2.66 + const LengthMap &ol;
2.67 + const NodeMap &pot;
2.68 + public :
2.69 + typedef typename LengthMap::KeyType KeyType;
2.70 + typedef typename LengthMap::ValueType ValueType;
2.71 +
2.72 + ValueType operator[](typename ResGraphType::Edge e) const {
2.73 + if (G.forward(e))
2.74 + return ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
2.75 + else
2.76 + return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);
2.77 + }
2.78 +
2.79 + ModLengthMap(const ResGraphType& _G,
2.80 + const LengthMap &o, const NodeMap &p) :
2.81 + G(_G), /*rev(_rev),*/ ol(o), pot(p){};
2.82 + };//ModLengthMap
2.83 +
2.84 +
2.85 + protected:
2.86 +
2.87 + //Input
2.88 + const Graph& G;
2.89 + const LengthMap& length;
2.90 + const SupplyMap& supply;//supply or demand of nodes
2.91 +
2.92 +
2.93 + //auxiliary variables
2.94 +
2.95 + //To store the flow
2.96 + EdgeIntMap flow;
2.97 + //To store the potentila (dual variables)
2.98 + typename Graph::template NodeMap<Length> potential;
2.99 + //To store excess-deficit values
2.100 + SupplyMap excess;
2.101 +
2.102 +
2.103 + Length total_length;
2.104 +
2.105 +
2.106 + public :
2.107 +
2.108 +
2.109 + MinCostFlow(Graph& _G, LengthMap& _length, SupplyMap& _supply) : G(_G),
2.110 + length(_length), supply(_supply), flow(_G), potential(_G){ }
2.111 +
2.112 +
2.113 + ///Runs the algorithm.
2.114 +
2.115 + ///Runs the algorithm.
2.116 + ///Returns k if there are at least k edge-disjoint paths from s to t.
2.117 + ///Otherwise it returns the number of found edge-disjoint paths from s to t.
2.118 + ///\todo May be it does make sense to be able to start with a nonzero
2.119 + /// feasible primal-dual solution pair as well.
2.120 + int run() {
2.121 +
2.122 + //Resetting variables from previous runs
2.123 + total_length = 0;
2.124 +
2.125 + FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
2.126 + flow.set(e,0);
2.127 + }
2.128 +
2.129 + //Initial value for delta
2.130 + Supply delta = 0;
2.131 +
2.132 + FOR_EACH_LOC(typename Graph::NodeIt, n, G){
2.133 + if (delta < abs(supply[e])){
2.134 + delta = abs(supply[e]);
2.135 + }
2.136 + excess.set(n,supply[e]);
2.137 + //Initialize the copy of the Dijkstra potential to zero
2.138 + potential.set(n,0);
2.139 + }
2.140 +
2.141 +
2.142 +
2.143 + //We need a residual graph which is uncapacitated
2.144 + ResGraphType res_graph(G, flow);
2.145 +
2.146 +
2.147 +
2.148 + ModLengthMap mod_length(res_graph, length, potential);
2.149 +
2.150 + Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
2.151 +
2.152 +
2.153 + int i;
2.154 + for (i=0; i<k; ++i){
2.155 + dijkstra.run(s);
2.156 + if (!dijkstra.reached(t)){
2.157 + //There are no k paths from s to t
2.158 + break;
2.159 + };
2.160 +
2.161 + //We have to copy the potential
2.162 + FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
2.163 + potential[n] += dijkstra.distMap()[n];
2.164 + }
2.165 +
2.166 + /*
2.167 + {
2.168 +
2.169 + typename ResGraphType::NodeIt n;
2.170 + for ( res_graph.first(n) ; res_graph.valid(n) ; res_graph.next(n) ) {
2.171 + potential[n] += dijkstra.distMap()[n];
2.172 + }
2.173 + }
2.174 + */
2.175 +
2.176 + //Augmenting on the sortest path
2.177 + Node n=t;
2.178 + ResGraphEdge e;
2.179 + while (n!=s){
2.180 + e = dijkstra.pred(n);
2.181 + n = dijkstra.predNode(n);
2.182 + res_graph.augment(e,delta);
2.183 + //Let's update the total length
2.184 + if (res_graph.forward(e))
2.185 + total_length += length[e];
2.186 + else
2.187 + total_length -= length[e];
2.188 + }
2.189 +
2.190 +
2.191 + }
2.192 +
2.193 +
2.194 + return i;
2.195 + }
2.196 +
2.197 +
2.198 +
2.199 +
2.200 + ///This function gives back the total length of the found paths.
2.201 + ///Assumes that \c run() has been run and nothing changed since then.
2.202 + Length totalLength(){
2.203 + return total_length;
2.204 + }
2.205 +
2.206 + ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
2.207 + ///be called before using this function.
2.208 + const EdgeIntMap &getFlow() const { return flow;}
2.209 +
2.210 + ///Returns a const reference to the NodeMap \c potential (the dual solution).
2.211 + /// \pre \ref run() must be called before using this function.
2.212 + const EdgeIntMap &getPotential() const { return potential;}
2.213 +
2.214 + ///This function checks, whether the given solution is optimal
2.215 + ///Running after a \c run() should return with true
2.216 + ///In this "state of the art" this only check optimality, doesn't bother with feasibility
2.217 + ///
2.218 + ///\todo Is this OK here?
2.219 + bool checkComplementarySlackness(){
2.220 + Length mod_pot;
2.221 + Length fl_e;
2.222 + FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
2.223 + //C^{\Pi}_{i,j}
2.224 + mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
2.225 + fl_e = flow[e];
2.226 + // std::cout << fl_e << std::endl;
2.227 + if (0<fl_e && fl_e<capacity[e]){
2.228 + if (mod_pot != 0)
2.229 + return false;
2.230 + }
2.231 + else{
2.232 + if (mod_pot > 0 && fl_e != 0)
2.233 + return false;
2.234 + if (mod_pot < 0 && fl_e != capacity[e])
2.235 + return false;
2.236 + }
2.237 + }
2.238 + return true;
2.239 + }
2.240 +
2.241 +
2.242 + }; //class MinCostFlow
2.243 +
2.244 + ///@}
2.245 +
2.246 +} //namespace hugo
2.247 +
2.248 +#endif //HUGO_MINCOSTFLOW_H