1.1 --- a/src/work/jacint/max_flow.h Wed May 19 16:06:57 2004 +0000
1.2 +++ b/src/work/jacint/max_flow.h Wed May 19 16:09:38 2004 +0000
1.3 @@ -25,22 +25,22 @@
1.4 ///This class provides various algorithms for finding a flow of
1.5 ///maximum value in a directed graph. The \e source node, the \e
1.6 ///target node, the \e capacity of the edges and the \e starting \e
1.7 - ///flow value of the edges can be passed to the algorithm through the
1.8 + ///flow value of the edges should be passed to the algorithm through the
1.9 ///constructor. It is possible to change these quantities using the
1.10 ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
1.11 ///\ref resetFlow. Before any subsequent runs of any algorithm of
1.12 - ///the class \ref resetFlow should be called, otherwise it will
1.13 - ///start from a maximum flow.
1.14 - ///After running an algorithm of the class, the maximum value of a
1.15 - ///value can be obtained by calling \ref flowValue(). The minimum
1.16 + ///the class \ref resetFlow should be called.
1.17 +
1.18 + ///After running an algorithm of the class, the actual flow value
1.19 + ///can be obtained by calling \ref flowValue(). The minimum
1.20 ///value cut can be written into a \c node map of \c bools by
1.21 ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
1.22 ///the inclusionwise minimum and maximum of the minimum value
1.23 ///cuts, resp.)
1.24 ///\param Graph The directed graph type the algorithm runs on.
1.25 ///\param Num The number type of the capacities and the flow values.
1.26 - ///\param CapMap The type of the capacity map.
1.27 - ///\param FlowMap The type of the flow map.
1.28 + ///\param CapMap The capacity map type.
1.29 + ///\param FlowMap The flow map type.
1.30 ///\author Marton Makai, Jacint Szabo
1.31 template <typename Graph, typename Num,
1.32 typename CapMap=typename Graph::template EdgeMap<Num>,
1.33 @@ -111,17 +111,46 @@
1.34 ///least the sum of the out-flows in every node except the \e source.
1.35 ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
1.36 ///set to the constant zero flow in the beginning of the algorithm in this case.
1.37 - enum flowEnum{
1.38 + enum FlowEnum{
1.39 ZERO_FLOW,
1.40 GEN_FLOW,
1.41 PRE_FLOW,
1.42 NO_FLOW
1.43 };
1.44
1.45 + enum StatusEnum {
1.46 + AFTER_NOTHING,
1.47 + AFTER_AUGMENTING,
1.48 + AFTER_PRE_FLOW_PHASE_1,
1.49 + AFTER_PRE_FLOW_PHASE_2
1.50 + };
1.51 +
1.52 + /// Don not needle this flag only if necessary.
1.53 + StatusEnum status;
1.54 + int number_of_augmentations;
1.55 +
1.56 +
1.57 + template<typename IntMap>
1.58 + class TrickyReachedMap {
1.59 + protected:
1.60 + IntMap* map;
1.61 + int* number_of_augmentations;
1.62 + public:
1.63 + TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
1.64 + map(&_map), number_of_augmentations(&_number_of_augmentations) { }
1.65 + void set(const Node& n, bool b) {
1.66 + map->set(n, *number_of_augmentations);
1.67 + }
1.68 + bool operator[](const Node& n) const {
1.69 + return (*map)[n]==*number_of_augmentations;
1.70 + }
1.71 + };
1.72 +
1.73 MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
1.74 FlowMap& _flow) :
1.75 g(&_G), s(_s), t(_t), capacity(&_capacity),
1.76 - flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
1.77 + flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
1.78 + status(AFTER_NOTHING), number_of_augmentations(0) { }
1.79
1.80 ///Runs a maximum flow algorithm.
1.81
1.82 @@ -132,11 +161,11 @@
1.83 /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.84 /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.85 /// - any map if \c fe is NO_FLOW.
1.86 - void run(flowEnum fe=ZERO_FLOW) {
1.87 + void run(FlowEnum fe=ZERO_FLOW) {
1.88 preflow(fe);
1.89 }
1.90
1.91 -
1.92 +
1.93 ///Runs a preflow algorithm.
1.94
1.95 ///Runs a preflow algorithm. The preflow algorithms provide the
1.96 @@ -146,7 +175,7 @@
1.97 /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.98 /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.99 /// - any map if \c fe is NO_FLOW.
1.100 - void preflow(flowEnum fe) {
1.101 + void preflow(FlowEnum fe) {
1.102 preflowPhase1(fe);
1.103 preflowPhase2();
1.104 }
1.105 @@ -173,7 +202,7 @@
1.106 /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.107 /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.108 /// - any map if \c fe is NO_FLOW.
1.109 - void preflowPhase1( flowEnum fe );
1.110 + void preflowPhase1(FlowEnum fe);
1.111
1.112 ///Runs the second phase of the preflow algorithm.
1.113
1.114 @@ -189,6 +218,7 @@
1.115 /// and augments the flow on if any.
1.116 /// The return value shows if the augmentation was succesful.
1.117 bool augmentOnShortestPath();
1.118 + bool augmentOnShortestPath2();
1.119
1.120 /// Starting from a flow, this method searches for an augmenting blocking
1.121 /// flow according to Dinits' algorithm and augments the flow on if any.
1.122 @@ -207,7 +237,7 @@
1.123 /// Returns the maximum value of a flow, by counting the
1.124 /// over-flow of the target node \ref t.
1.125 /// It can be called already after running \ref preflowPhase1.
1.126 - Num flowValue() {
1.127 + Num flowValue() const {
1.128 Num a=0;
1.129 FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
1.130 FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
1.131 @@ -228,14 +258,31 @@
1.132 /// of the class. This enables us to determine which methods are valid
1.133 /// for MinCut computation
1.134 template<typename _CutMap>
1.135 - void actMinCut(_CutMap& M) {
1.136 + void actMinCut(_CutMap& M) const {
1.137 NodeIt v;
1.138 - for(g->first(v); g->valid(v); g->next(v)) {
1.139 - if ( level[v] < n ) {
1.140 - M.set(v,false);
1.141 - } else {
1.142 - M.set(v,true);
1.143 + switch (status) {
1.144 + case AFTER_PRE_FLOW_PHASE_1:
1.145 + for(g->first(v); g->valid(v); g->next(v)) {
1.146 + if (level[v] < n) {
1.147 + M.set(v, false);
1.148 + } else {
1.149 + M.set(v, true);
1.150 + }
1.151 }
1.152 + break;
1.153 + case AFTER_PRE_FLOW_PHASE_2:
1.154 + case AFTER_NOTHING:
1.155 + minMinCut(M);
1.156 + break;
1.157 + case AFTER_AUGMENTING:
1.158 + for(g->first(v); g->valid(v); g->next(v)) {
1.159 + if (level[v]) {
1.160 + M.set(v, true);
1.161 + } else {
1.162 + M.set(v, false);
1.163 + }
1.164 + }
1.165 + break;
1.166 }
1.167 }
1.168
1.169 @@ -247,8 +294,7 @@
1.170 ///\pre M should be a node map of bools initialized to false.
1.171 ///\pre \c flow must be a maximum flow.
1.172 template<typename _CutMap>
1.173 - void minMinCut(_CutMap& M) {
1.174 -
1.175 + void minMinCut(_CutMap& M) const {
1.176 std::queue<Node> queue;
1.177
1.178 M.set(s,true);
1.179 @@ -286,7 +332,7 @@
1.180 ///\pre M should be a node map of bools initialized to false.
1.181 ///\pre \c flow must be a maximum flow.
1.182 template<typename _CutMap>
1.183 - void maxMinCut(_CutMap& M) {
1.184 + void maxMinCut(_CutMap& M) const {
1.185
1.186 NodeIt v;
1.187 for(g->first(v) ; g->valid(v); g->next(v)) {
1.188 @@ -328,31 +374,31 @@
1.189 ///\pre M should be a node map of bools initialized to false.
1.190 ///\pre \c flow must be a maximum flow.
1.191 template<typename CutMap>
1.192 - void minCut(CutMap& M) { minMinCut(M); }
1.193 + void minCut(CutMap& M) const { minMinCut(M); }
1.194
1.195 ///Resets the source node to \c _s.
1.196
1.197 ///Resets the source node to \c _s.
1.198 ///
1.199 - void resetSource(Node _s) { s=_s; }
1.200 + void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
1.201
1.202 ///Resets the target node to \c _t.
1.203
1.204 ///Resets the target node to \c _t.
1.205 ///
1.206 - void resetTarget(Node _t) { t=_t; }
1.207 + void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
1.208
1.209 /// Resets the edge map of the capacities to _cap.
1.210
1.211 /// Resets the edge map of the capacities to _cap.
1.212 ///
1.213 - void resetCap(const CapMap& _cap) { capacity=&_cap; }
1.214 + void resetCap(const CapMap& _cap) { capacity=&_cap; status=AFTER_NOTHING; }
1.215
1.216 /// Resets the edge map of the flows to _flow.
1.217
1.218 /// Resets the edge map of the flows to _flow.
1.219 ///
1.220 - void resetFlow(FlowMap& _flow) { flow=&_flow; }
1.221 + void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
1.222
1.223
1.224 private:
1.225 @@ -434,7 +480,7 @@
1.226 }
1.227
1.228
1.229 - void preflowPreproc(flowEnum fe, VecStack& active,
1.230 + void preflowPreproc(FlowEnum fe, VecStack& active,
1.231 VecNode& level_list, NNMap& left, NNMap& right)
1.232 {
1.233 std::queue<Node> bfs_queue;
1.234 @@ -638,8 +684,9 @@
1.235 void set(const typename MapGraphWrapper::Node& n, int a) {
1.236 dist.set(n, a);
1.237 }
1.238 - int operator[](const typename MapGraphWrapper::Node& n)
1.239 - { return dist[n]; }
1.240 + int operator[](const typename MapGraphWrapper::Node& n) const {
1.241 + return dist[n];
1.242 + }
1.243 // int get(const typename MapGraphWrapper::Node& n) const {
1.244 // return dist[n]; }
1.245 // bool get(const typename MapGraphWrapper::Edge& e) const {
1.246 @@ -653,7 +700,7 @@
1.247
1.248
1.249 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.250 - void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1( flowEnum fe )
1.251 + void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1(FlowEnum fe)
1.252 {
1.253
1.254 int heur0=(int)(H0*n); //time while running 'bound decrease'
1.255 @@ -766,6 +813,8 @@
1.256 }
1.257 }
1.258 }
1.259 +
1.260 + status=AFTER_PRE_FLOW_PHASE_1;
1.261 }
1.262
1.263
1.264 @@ -830,6 +879,8 @@
1.265 }
1.266 } // if stack[b] is nonempty
1.267 } // while(true)
1.268 +
1.269 + status=AFTER_PRE_FLOW_PHASE_2;
1.270 }
1.271
1.272
1.273 @@ -878,13 +929,67 @@
1.274 }
1.275 }
1.276
1.277 + status=AFTER_AUGMENTING;
1.278 return _augment;
1.279 }
1.280
1.281
1.282 + template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.283 + bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
1.284 + {
1.285 + ResGW res_graph(*g, *capacity, *flow);
1.286 + bool _augment=false;
1.287
1.288 + if (status!=AFTER_AUGMENTING) {
1.289 + FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, -1);
1.290 + number_of_augmentations=0;
1.291 + } else {
1.292 + ++number_of_augmentations;
1.293 + }
1.294 + TrickyReachedMap<ReachedMap>
1.295 + tricky_reached_map(level, number_of_augmentations);
1.296 + //ReachedMap level(res_graph);
1.297 +// FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.298 + BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >
1.299 + bfs(res_graph, tricky_reached_map);
1.300 + bfs.pushAndSetReached(s);
1.301
1.302 + typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
1.303 + pred.set(s, INVALID);
1.304
1.305 + typename ResGW::template NodeMap<Num> free(res_graph);
1.306 +
1.307 + //searching for augmenting path
1.308 + while ( !bfs.finished() ) {
1.309 + ResGWOutEdgeIt e=bfs;
1.310 + if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1.311 + Node v=res_graph.tail(e);
1.312 + Node w=res_graph.head(e);
1.313 + pred.set(w, e);
1.314 + if (res_graph.valid(pred[v])) {
1.315 + free.set(w, std::min(free[v], res_graph.resCap(e)));
1.316 + } else {
1.317 + free.set(w, res_graph.resCap(e));
1.318 + }
1.319 + if (res_graph.head(e)==t) { _augment=true; break; }
1.320 + }
1.321 +
1.322 + ++bfs;
1.323 + } //end of searching augmenting path
1.324 +
1.325 + if (_augment) {
1.326 + Node n=t;
1.327 + Num augment_value=free[t];
1.328 + while (res_graph.valid(pred[n])) {
1.329 + ResGWEdge e=pred[n];
1.330 + res_graph.augment(e, augment_value);
1.331 + n=res_graph.tail(e);
1.332 + }
1.333 + }
1.334 +
1.335 + status=AFTER_AUGMENTING;
1.336 + return _augment;
1.337 + }
1.338
1.339
1.340 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.341 @@ -999,6 +1104,7 @@
1.342
1.343 }
1.344
1.345 + status=AFTER_AUGMENTING;
1.346 return _augment;
1.347 }
1.348
1.349 @@ -1129,6 +1235,7 @@
1.350
1.351 } //while (__augment)
1.352
1.353 + status=AFTER_AUGMENTING;
1.354 return _augment;
1.355 }
1.356