src/work/edmonds_karp.h
author alpar
Sun, 21 Mar 2004 14:58:48 +0000
changeset 223 02948c4c68e1
parent 198 5cec393baade
child 243 a85fd87460e3
permissions -rw-r--r--
Dijkstra, bin_heap, fib_heap added to the doc.
marci@174
     1
// -*- c++ -*-
marci@174
     2
#ifndef EDMONDS_KARP_H
marci@174
     3
#define EDMONDS_KARP_H
marci@174
     4
marci@174
     5
#include <algorithm>
marci@174
     6
#include <list>
marci@174
     7
#include <iterator>
marci@174
     8
marci@174
     9
#include <bfs_iterator.h>
marci@174
    10
#include <invalid.h>
marci@174
    11
marci@174
    12
namespace hugo {
marci@174
    13
marci@174
    14
  template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
marci@174
    15
  class ResGraph {
marci@174
    16
  public:
marci@174
    17
    typedef typename Graph::Node Node;
marci@174
    18
    typedef typename Graph::NodeIt NodeIt;
marci@174
    19
  private:
marci@174
    20
    typedef typename Graph::SymEdgeIt OldSymEdgeIt;
marci@174
    21
    const Graph& G;
marci@174
    22
    FlowMap& flow;
marci@174
    23
    const CapacityMap& capacity;
marci@174
    24
  public:
marci@174
    25
    ResGraph(const Graph& _G, FlowMap& _flow, 
marci@174
    26
	     const CapacityMap& _capacity) : 
marci@174
    27
      G(_G), flow(_flow), capacity(_capacity) { }
marci@174
    28
marci@174
    29
    class Edge; 
marci@174
    30
    class OutEdgeIt; 
marci@174
    31
    friend class Edge; 
marci@174
    32
    friend class OutEdgeIt; 
marci@174
    33
marci@174
    34
    class Edge {
marci@174
    35
      friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
marci@174
    36
    protected:
marci@174
    37
      const ResGraph<Graph, Number, FlowMap, CapacityMap>* resG;
marci@174
    38
      OldSymEdgeIt sym;
marci@174
    39
    public:
marci@174
    40
      Edge() { } 
marci@174
    41
      //Edge(const Edge& e) : resG(e.resG), sym(e.sym) { }
marci@174
    42
      Number free() const { 
marci@174
    43
	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
marci@174
    44
	  return (resG->capacity.get(sym)-resG->flow.get(sym)); 
marci@174
    45
	} else { 
marci@174
    46
	  return (resG->flow.get(sym)); 
marci@174
    47
	}
marci@174
    48
      }
marci@174
    49
      bool valid() const { return sym.valid(); }
marci@174
    50
      void augment(Number a) const {
marci@174
    51
	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
marci@174
    52
	  resG->flow.set(sym, resG->flow.get(sym)+a);
marci@174
    53
	  //resG->flow[sym]+=a;
marci@174
    54
	} else { 
marci@174
    55
	  resG->flow.set(sym, resG->flow.get(sym)-a);
marci@174
    56
	  //resG->flow[sym]-=a;
marci@174
    57
	}
marci@174
    58
      }
marci@174
    59
    };
marci@174
    60
marci@174
    61
    class OutEdgeIt : public Edge {
marci@174
    62
      friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
marci@174
    63
    public:
marci@174
    64
      OutEdgeIt() { }
marci@174
    65
      //OutEdgeIt(const OutEdgeIt& e) { resG=e.resG; sym=e.sym; }
marci@174
    66
    private:
marci@174
    67
      OutEdgeIt(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _resG, Node v) { 
marci@174
    68
      	resG=&_resG;
marci@174
    69
	sym=resG->G.template first<OldSymEdgeIt>(v);
marci@174
    70
	while( sym.valid() && !(free()>0) ) { ++sym; }
marci@174
    71
      }
marci@174
    72
    public:
marci@174
    73
      OutEdgeIt& operator++() { 
marci@174
    74
	++sym; 
marci@174
    75
	while( sym.valid() && !(free()>0) ) { ++sym; }
marci@174
    76
	return *this; 
marci@174
    77
      }
marci@174
    78
    };
marci@174
    79
marci@174
    80
    void /*getF*/first(OutEdgeIt& e, Node v) const { 
marci@174
    81
      e=OutEdgeIt(*this, v); 
marci@174
    82
    }
marci@174
    83
    void /*getF*/first(NodeIt& v) const { G./*getF*/first(v); }
marci@174
    84
    
marci@174
    85
    template< typename It >
marci@174
    86
    It first() const { 
marci@174
    87
      It e;      
marci@174
    88
      /*getF*/first(e);
marci@174
    89
      return e; 
marci@174
    90
    }
marci@174
    91
marci@174
    92
    template< typename It >
marci@174
    93
    It first(Node v) const { 
marci@174
    94
      It e;
marci@174
    95
      /*getF*/first(e, v);
marci@174
    96
      return e; 
marci@174
    97
    }
marci@174
    98
marci@174
    99
    Node tail(Edge e) const { return G.aNode(e.sym); }
marci@174
   100
    Node head(Edge e) const { return G.bNode(e.sym); }
marci@174
   101
marci@174
   102
    Node aNode(OutEdgeIt e) const { return G.aNode(e.sym); }
marci@174
   103
    Node bNode(OutEdgeIt e) const { return G.bNode(e.sym); }
marci@174
   104
marci@174
   105
    int id(Node v) const { return G.id(v); }
marci@174
   106
marci@174
   107
    template <typename S>
marci@174
   108
    class NodeMap {
marci@174
   109
      typename Graph::NodeMap<S> node_map; 
marci@174
   110
    public:
marci@174
   111
      NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
marci@174
   112
      NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
marci@174
   113
      void set(Node nit, S a) { node_map.set(nit, a); }
marci@174
   114
      S get(Node nit) const { return node_map.get(nit); }
marci@174
   115
      S& operator[](Node nit) { return node_map[nit]; } 
marci@174
   116
      const S& operator[](Node nit) const { return node_map[nit]; } 
marci@174
   117
    };
marci@174
   118
marci@174
   119
  };
marci@174
   120
marci@174
   121
marci@174
   122
  template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
marci@174
   123
  class ResGraph2 {
marci@174
   124
  public:
marci@174
   125
    typedef typename Graph::Node Node;
marci@174
   126
    typedef typename Graph::NodeIt NodeIt;
marci@174
   127
  private:
marci@174
   128
    //typedef typename Graph::SymEdgeIt OldSymEdgeIt;
marci@174
   129
    typedef typename Graph::OutEdgeIt OldOutEdgeIt;
marci@174
   130
    typedef typename Graph::InEdgeIt OldInEdgeIt;
marci@174
   131
    
marci@174
   132
    const Graph& G;
marci@174
   133
    FlowMap& flow;
marci@174
   134
    const CapacityMap& capacity;
marci@174
   135
  public:
marci@174
   136
    ResGraph2(const Graph& _G, FlowMap& _flow, 
marci@174
   137
	     const CapacityMap& _capacity) : 
marci@174
   138
      G(_G), flow(_flow), capacity(_capacity) { }
marci@174
   139
marci@174
   140
    class Edge; 
marci@174
   141
    class OutEdgeIt; 
marci@174
   142
    friend class Edge; 
marci@174
   143
    friend class OutEdgeIt; 
marci@174
   144
marci@174
   145
    class Edge {
marci@174
   146
      friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
marci@174
   147
    protected:
marci@174
   148
      const ResGraph2<Graph, Number, FlowMap, CapacityMap>* resG;
marci@174
   149
      //OldSymEdgeIt sym;
marci@174
   150
      OldOutEdgeIt out;
marci@174
   151
      OldInEdgeIt in;
marci@174
   152
      bool out_or_in; //true, iff out
marci@174
   153
    public:
marci@174
   154
      Edge() : out_or_in(true) { } 
marci@174
   155
      Number free() const { 
marci@174
   156
	if (out_or_in) { 
marci@174
   157
	  return (resG->capacity.get(out)-resG->flow.get(out)); 
marci@174
   158
	} else { 
marci@174
   159
	  return (resG->flow.get(in)); 
marci@174
   160
	}
marci@174
   161
      }
marci@174
   162
      bool valid() const { 
marci@174
   163
	return out_or_in && out.valid() || in.valid(); }
marci@174
   164
      void augment(Number a) const {
marci@174
   165
	if (out_or_in) { 
marci@174
   166
	  resG->flow.set(out, resG->flow.get(out)+a);
marci@174
   167
	} else { 
marci@174
   168
	  resG->flow.set(in, resG->flow.get(in)-a);
marci@174
   169
	}
marci@174
   170
      }
marci@174
   171
    };
marci@174
   172
marci@174
   173
    class OutEdgeIt : public Edge {
marci@174
   174
      friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
marci@174
   175
    public:
marci@174
   176
      OutEdgeIt() { }
marci@174
   177
    private:
marci@174
   178
      OutEdgeIt(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _resG, Node v) { 
marci@174
   179
      	resG=&_resG;
marci@174
   180
	out=resG->G.template first<OldOutEdgeIt>(v);
marci@174
   181
	while( out.valid() && !(free()>0) ) { ++out; }
marci@174
   182
	if (!out.valid()) {
marci@174
   183
	  out_or_in=0;
marci@174
   184
	  in=resG->G.template first<OldInEdgeIt>(v);
marci@174
   185
	  while( in.valid() && !(free()>0) ) { ++in; }
marci@174
   186
	}
marci@174
   187
      }
marci@174
   188
    public:
marci@174
   189
      OutEdgeIt& operator++() { 
marci@174
   190
	if (out_or_in) {
marci@174
   191
	  Node v=resG->G.aNode(out);
marci@174
   192
	  ++out;
marci@174
   193
	  while( out.valid() && !(free()>0) ) { ++out; }
marci@174
   194
	  if (!out.valid()) {
marci@174
   195
	    out_or_in=0;
marci@174
   196
	    in=resG->G.template first<OldInEdgeIt>(v);
marci@174
   197
	    while( in.valid() && !(free()>0) ) { ++in; }
marci@174
   198
	  }
marci@174
   199
	} else {
marci@174
   200
	  ++in;
marci@174
   201
	  while( in.valid() && !(free()>0) ) { ++in; } 
marci@174
   202
	}
marci@174
   203
	return *this; 
marci@174
   204
      }
marci@174
   205
    };
marci@174
   206
marci@174
   207
    void /*getF*/first(OutEdgeIt& e, Node v) const { 
marci@174
   208
      e=OutEdgeIt(*this, v); 
marci@174
   209
    }
marci@174
   210
    void /*getF*/first(NodeIt& v) const { G./*getF*/first(v); }
marci@174
   211
    
marci@174
   212
    template< typename It >
marci@174
   213
    It first() const { 
marci@174
   214
      It e;
marci@174
   215
      /*getF*/first(e);
marci@174
   216
      return e; 
marci@174
   217
    }
marci@174
   218
marci@174
   219
    template< typename It >
marci@174
   220
    It first(Node v) const { 
marci@174
   221
      It e;
marci@174
   222
      /*getF*/first(e, v);
marci@174
   223
      return e; 
marci@174
   224
    }
marci@174
   225
marci@174
   226
    Node tail(Edge e) const { 
marci@174
   227
      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
marci@174
   228
    Node head(Edge e) const { 
marci@174
   229
      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
marci@174
   230
marci@174
   231
    Node aNode(OutEdgeIt e) const { 
marci@174
   232
      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
marci@174
   233
    Node bNode(OutEdgeIt e) const { 
marci@174
   234
      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
marci@174
   235
marci@174
   236
    int id(Node v) const { return G.id(v); }
marci@174
   237
marci@174
   238
    template <typename S>
marci@174
   239
    class NodeMap {
marci@174
   240
      typename Graph::NodeMap<S> node_map; 
marci@174
   241
    public:
marci@174
   242
      NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
marci@174
   243
      NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
marci@174
   244
      void set(Node nit, S a) { node_map.set(nit, a); }
marci@174
   245
      S get(Node nit) const { return node_map.get(nit); }
marci@174
   246
    };
marci@174
   247
  };
marci@174
   248
marci@174
   249
marci@174
   250
  template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
marci@174
   251
  class MaxFlow {
marci@174
   252
  public:
marci@174
   253
    typedef typename Graph::Node Node;
marci@174
   254
    typedef typename Graph::Edge Edge;
marci@174
   255
    typedef typename Graph::EdgeIt EdgeIt;
marci@174
   256
    typedef typename Graph::OutEdgeIt OutEdgeIt;
marci@174
   257
    typedef typename Graph::InEdgeIt InEdgeIt;
marci@174
   258
marci@174
   259
  private:
marci@174
   260
    const Graph* G;
marci@174
   261
    Node s;
marci@174
   262
    Node t;
marci@174
   263
    FlowMap* flow;
marci@174
   264
    const CapacityMap* capacity;
marci@174
   265
    typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
marci@174
   266
    typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
marci@174
   267
    typedef typename AugGraph::Edge AugEdge;
marci@174
   268
marci@174
   269
  public:
marci@174
   270
    MaxFlow(const Graph& _G, Node _s, Node _t, FlowMap& _flow, const CapacityMap& _capacity) : 
marci@191
   271
      G(&_G), s(_s), t(_t), flow(&_flow), capacity(&_capacity) { }
marci@174
   272
    bool augmentOnShortestPath() {
marci@174
   273
      AugGraph res_graph(*G, *flow, *capacity);
marci@174
   274
      bool _augment=false;
marci@174
   275
      
marci@174
   276
      typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@174
   277
      BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
marci@174
   278
      res_bfs.pushAndSetReached(s);
marci@174
   279
	
marci@174
   280
      typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@174
   281
      pred.set(s, AugEdge(INVALID));
marci@174
   282
      
marci@174
   283
      typename AugGraph::NodeMap<Number> free(res_graph);
marci@174
   284
	
marci@174
   285
      //searching for augmenting path
marci@174
   286
      while ( !res_bfs.finished() ) { 
marci@191
   287
	AugOutEdgeIt e=res_bfs;
marci@174
   288
	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
marci@174
   289
	  Node v=res_graph.tail(e);
marci@174
   290
	  Node w=res_graph.head(e);
marci@174
   291
	  pred.set(w, e);
marci@174
   292
	  if (res_graph.valid(pred.get(v))) {
marci@174
   293
	    free.set(w, std::min(free.get(v), res_graph.free(e)));
marci@174
   294
	  } else {
marci@174
   295
	    free.set(w, res_graph.free(e)); 
marci@174
   296
	  }
marci@174
   297
	  if (res_graph.head(e)==t) { _augment=true; break; }
marci@174
   298
	}
marci@174
   299
	
marci@174
   300
	++res_bfs;
marci@174
   301
      } //end of searching augmenting path
marci@174
   302
marci@174
   303
      if (_augment) {
marci@174
   304
	Node n=t;
marci@174
   305
	Number augment_value=free.get(t);
marci@174
   306
	while (res_graph.valid(pred.get(n))) { 
marci@174
   307
	  AugEdge e=pred.get(n);
marci@174
   308
	  res_graph.augment(e, augment_value); 
marci@174
   309
	  n=res_graph.tail(e);
marci@174
   310
	}
marci@174
   311
      }
marci@174
   312
marci@174
   313
      return _augment;
marci@174
   314
    }
marci@174
   315
marci@191
   316
    template<typename MutableGraph> bool augmentOnBlockingFlow() {      
marci@174
   317
      bool _augment=false;
marci@174
   318
marci@174
   319
      AugGraph res_graph(*G, *flow, *capacity);
marci@174
   320
marci@174
   321
      typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@174
   322
      BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
marci@174
   323
marci@174
   324
      bfs.pushAndSetReached(s);
marci@174
   325
      typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
marci@174
   326
      while ( !bfs.finished() ) { 
marci@191
   327
	AugOutEdgeIt e=bfs;
marci@174
   328
	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
marci@174
   329
	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
marci@174
   330
	}
marci@174
   331
	
marci@174
   332
	++bfs;
marci@174
   333
      } //computing distances from s in the residual graph
marci@174
   334
marci@174
   335
      MutableGraph F;
marci@174
   336
      typename AugGraph::NodeMap<typename MutableGraph::Node> 
marci@174
   337
	res_graph_to_F(res_graph);
marci@174
   338
      for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
marci@174
   339
	res_graph_to_F.set(n, F.addNode());
marci@174
   340
      }
marci@174
   341
      
marci@174
   342
      typename MutableGraph::Node sF=res_graph_to_F.get(s);
marci@174
   343
      typename MutableGraph::Node tF=res_graph_to_F.get(t);
marci@174
   344
marci@174
   345
      typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
marci@174
   346
      typename MutableGraph::EdgeMap<Number> residual_capacity(F);
marci@174
   347
marci@174
   348
      //Making F to the graph containing the edges of the residual graph 
marci@174
   349
      //which are in some shortest paths
marci@174
   350
      for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
marci@174
   351
	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
marci@174
   352
	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
marci@174
   353
	  original_edge.update();
marci@174
   354
	  original_edge.set(f, e);
marci@174
   355
	  residual_capacity.update();
marci@174
   356
	  residual_capacity.set(f, res_graph.free(e));
marci@174
   357
	} 
marci@174
   358
      }
marci@174
   359
marci@174
   360
      bool __augment=true;
marci@174
   361
marci@174
   362
      while (__augment) {
marci@174
   363
	__augment=false;
marci@174
   364
	//computing blocking flow with dfs
marci@174
   365
	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
marci@174
   366
	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
marci@174
   367
	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
marci@174
   368
	pred.set(sF, typename MutableGraph::Edge(INVALID));
marci@174
   369
	//invalid iterators for sources
marci@174
   370
marci@174
   371
	typename MutableGraph::NodeMap<Number> free(F);
marci@174
   372
marci@174
   373
	dfs.pushAndSetReached(sF);      
marci@174
   374
	while (!dfs.finished()) {
marci@174
   375
	  ++dfs;
marci@174
   376
	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
marci@174
   377
	    if (dfs.isBNodeNewlyReached()) {
marci@174
   378
	      typename MutableGraph::Node v=F.aNode(dfs);
marci@174
   379
	      typename MutableGraph::Node w=F.bNode(dfs);
marci@174
   380
	      pred.set(w, dfs);
marci@174
   381
	      if (F.valid(pred.get(v))) {
marci@174
   382
		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
marci@174
   383
	      } else {
marci@174
   384
		free.set(w, residual_capacity.get(dfs)); 
marci@174
   385
	      }
marci@174
   386
	      if (w==tF) { 
marci@174
   387
		__augment=true; 
marci@174
   388
		_augment=true;
marci@174
   389
		break; 
marci@174
   390
	      }
marci@174
   391
	      
marci@174
   392
	    } else {
marci@174
   393
	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
marci@174
   394
	    }
marci@174
   395
	  } 
marci@174
   396
	}
marci@174
   397
marci@174
   398
	if (__augment) {
marci@174
   399
	  typename MutableGraph::Node n=tF;
marci@174
   400
	  Number augment_value=free.get(tF);
marci@174
   401
	  while (F.valid(pred.get(n))) { 
marci@174
   402
	    typename MutableGraph::Edge e=pred.get(n);
marci@174
   403
	    res_graph.augment(original_edge.get(e), augment_value); 
marci@174
   404
	    n=F.tail(e);
marci@174
   405
	    if (residual_capacity.get(e)==augment_value) 
marci@174
   406
	      F.erase(e); 
marci@174
   407
	    else 
marci@174
   408
	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
marci@174
   409
	  }
marci@174
   410
	}
marci@174
   411
	
marci@174
   412
      }
marci@174
   413
            
marci@174
   414
      return _augment;
marci@174
   415
    }
marci@206
   416
    template<typename MutableGraph> bool augmentOnBlockingFlow1() {      
marci@206
   417
      bool _augment=false;
marci@206
   418
marci@206
   419
      AugGraph res_graph(*G, *flow, *capacity);
marci@206
   420
marci@206
   421
      //bfs for distances on the residual graph
marci@206
   422
      typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@206
   423
      BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
marci@206
   424
      bfs.pushAndSetReached(s);
marci@206
   425
      typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
marci@206
   426
marci@206
   427
      //F will contain the physical copy of the residual graph
marci@206
   428
      //with the set of edges which areon shortest paths
marci@206
   429
      MutableGraph F;
marci@206
   430
      typename AugGraph::NodeMap<typename MutableGraph::Node> 
marci@206
   431
	res_graph_to_F(res_graph);
marci@206
   432
      for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
marci@206
   433
	res_graph_to_F.set(n, F.addNode());
marci@206
   434
      }
marci@206
   435
      typename MutableGraph::Node sF=res_graph_to_F.get(s);
marci@206
   436
      typename MutableGraph::Node tF=res_graph_to_F.get(t);
marci@206
   437
      typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
marci@206
   438
      typename MutableGraph::EdgeMap<Number> residual_capacity(F);
marci@206
   439
marci@206
   440
      while ( !bfs.finished() ) { 
marci@206
   441
	AugOutEdgeIt e=bfs;
marci@206
   442
	if (res_graph.valid(e)) {
marci@206
   443
	  if (bfs.isBNodeNewlyReached()) {
marci@206
   444
	    dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
marci@206
   445
	    typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
marci@206
   446
	    original_edge.update();
marci@206
   447
	    original_edge.set(f, e);
marci@206
   448
	    residual_capacity.update();
marci@206
   449
	    residual_capacity.set(f, res_graph.free(e));
marci@206
   450
	  } else {
marci@206
   451
	    if (dist.get(res_graph.head(e))==(dist.get(res_graph.tail(e))+1)) {
marci@206
   452
	      typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
marci@206
   453
	      original_edge.update();
marci@206
   454
	      original_edge.set(f, e);
marci@206
   455
	      residual_capacity.update();
marci@206
   456
	      residual_capacity.set(f, res_graph.free(e));
marci@206
   457
	    }
marci@206
   458
	  }
marci@206
   459
	}
marci@206
   460
	++bfs;
marci@206
   461
      } //computing distances from s in the residual graph
marci@206
   462
marci@206
   463
      bool __augment=true;
marci@206
   464
marci@206
   465
      while (__augment) {
marci@206
   466
	__augment=false;
marci@206
   467
	//computing blocking flow with dfs
marci@206
   468
	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
marci@206
   469
	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
marci@206
   470
	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
marci@206
   471
	pred.set(sF, typename MutableGraph::Edge(INVALID));
marci@206
   472
	//invalid iterators for sources
marci@206
   473
marci@206
   474
	typename MutableGraph::NodeMap<Number> free(F);
marci@206
   475
marci@206
   476
	dfs.pushAndSetReached(sF);      
marci@206
   477
	while (!dfs.finished()) {
marci@206
   478
	  ++dfs;
marci@206
   479
	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
marci@206
   480
	    if (dfs.isBNodeNewlyReached()) {
marci@206
   481
	      typename MutableGraph::Node v=F.aNode(dfs);
marci@206
   482
	      typename MutableGraph::Node w=F.bNode(dfs);
marci@206
   483
	      pred.set(w, dfs);
marci@206
   484
	      if (F.valid(pred.get(v))) {
marci@206
   485
		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
marci@206
   486
	      } else {
marci@206
   487
		free.set(w, residual_capacity.get(dfs)); 
marci@206
   488
	      }
marci@206
   489
	      if (w==tF) { 
marci@206
   490
		__augment=true; 
marci@206
   491
		_augment=true;
marci@206
   492
		break; 
marci@206
   493
	      }
marci@206
   494
	      
marci@206
   495
	    } else {
marci@206
   496
	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
marci@206
   497
	    }
marci@206
   498
	  } 
marci@206
   499
	}
marci@206
   500
marci@206
   501
	if (__augment) {
marci@206
   502
	  typename MutableGraph::Node n=tF;
marci@206
   503
	  Number augment_value=free.get(tF);
marci@206
   504
	  while (F.valid(pred.get(n))) { 
marci@206
   505
	    typename MutableGraph::Edge e=pred.get(n);
marci@206
   506
	    res_graph.augment(original_edge.get(e), augment_value); 
marci@206
   507
	    n=F.tail(e);
marci@206
   508
	    if (residual_capacity.get(e)==augment_value) 
marci@206
   509
	      F.erase(e); 
marci@206
   510
	    else 
marci@206
   511
	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
marci@206
   512
	  }
marci@206
   513
	}
marci@206
   514
	
marci@206
   515
      }
marci@206
   516
            
marci@206
   517
      return _augment;
marci@206
   518
    }
marci@174
   519
    bool augmentOnBlockingFlow2() {
marci@174
   520
      bool _augment=false;
marci@174
   521
marci@174
   522
      //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
marci@174
   523
      typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
marci@174
   524
      typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
marci@174
   525
      typedef typename EAugGraph::Edge EAugEdge;
marci@174
   526
marci@174
   527
      EAugGraph res_graph(*G, *flow, *capacity);
marci@174
   528
marci@174
   529
      //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
marci@174
   530
      BfsIterator4< 
marci@174
   531
	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
marci@175
   532
	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt, 
marci@174
   533
	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
marci@174
   534
      
marci@191
   535
      bfs.pushAndSetReached(s);
marci@174
   536
marci@174
   537
      typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
marci@174
   538
	NodeMap<int>& dist=res_graph.dist;
marci@174
   539
marci@174
   540
      while ( !bfs.finished() ) {
marci@175
   541
	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
marci@174
   542
	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
marci@174
   543
	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
marci@174
   544
	}
marci@174
   545
	++bfs;	
marci@174
   546
      } //computing distances from s in the residual graph
marci@174
   547
marci@174
   548
      bool __augment=true;
marci@174
   549
marci@174
   550
      while (__augment) {
marci@174
   551
marci@174
   552
	__augment=false;
marci@174
   553
	//computing blocking flow with dfs
marci@174
   554
	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
marci@174
   555
	DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap > 
marci@174
   556
	  dfs(res_graph);
marci@174
   557
	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph); 
marci@174
   558
	pred.set(s, EAugEdge(INVALID));
marci@174
   559
	//invalid iterators for sources
marci@174
   560
marci@174
   561
	typename EAugGraph::NodeMap<Number> free(res_graph);
marci@174
   562
marci@174
   563
	dfs.pushAndSetReached(s);
marci@174
   564
	while (!dfs.finished()) {
marci@174
   565
	  ++dfs;
marci@174
   566
	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
marci@174
   567
	    if (dfs.isBNodeNewlyReached()) {
marci@174
   568
	  
marci@174
   569
	      typename EAugGraph::Node v=res_graph.aNode(dfs);
marci@174
   570
	      typename EAugGraph::Node w=res_graph.bNode(dfs);
marci@174
   571
marci@174
   572
	      pred.set(w, EAugOutEdgeIt(dfs));
marci@174
   573
	      if (res_graph.valid(pred.get(v))) {
marci@191
   574
		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
marci@174
   575
	      } else {
marci@191
   576
		free.set(w, res_graph.free(dfs)); 
marci@174
   577
	      }
marci@174
   578
	      
marci@174
   579
	      if (w==t) { 
marci@174
   580
		__augment=true; 
marci@174
   581
		_augment=true;
marci@174
   582
		break; 
marci@174
   583
	      }
marci@174
   584
	    } else {
marci@174
   585
	      res_graph.erase(dfs);
marci@174
   586
	    }
marci@174
   587
	  } 
marci@174
   588
marci@174
   589
	}
marci@174
   590
marci@174
   591
	if (__augment) {
marci@174
   592
	  typename EAugGraph::Node n=t;
marci@174
   593
	  Number augment_value=free.get(t);
marci@174
   594
	  while (res_graph.valid(pred.get(n))) { 
marci@174
   595
	    EAugEdge e=pred.get(n);
marci@174
   596
	    res_graph.augment(e, augment_value);
marci@174
   597
	    n=res_graph.tail(e);
marci@174
   598
	    if (res_graph.free(e)==0)
marci@174
   599
	      res_graph.erase(e);
marci@174
   600
	  }
marci@174
   601
	}
marci@174
   602
      
marci@174
   603
      }
marci@174
   604
            
marci@174
   605
      return _augment;
marci@174
   606
    }
marci@174
   607
    void run() {
marci@174
   608
      //int num_of_augmentations=0;
marci@174
   609
      while (augmentOnShortestPath()) { 
marci@174
   610
	//while (augmentOnBlockingFlow<MutableGraph>()) { 
marci@174
   611
	//std::cout << ++num_of_augmentations << " ";
marci@174
   612
	//std::cout<<std::endl;
marci@174
   613
      } 
marci@174
   614
    }
marci@174
   615
    template<typename MutableGraph> void run() {
marci@174
   616
      //int num_of_augmentations=0;
marci@174
   617
      //while (augmentOnShortestPath()) { 
marci@174
   618
	while (augmentOnBlockingFlow<MutableGraph>()) { 
marci@174
   619
	//std::cout << ++num_of_augmentations << " ";
marci@174
   620
	//std::cout<<std::endl;
marci@174
   621
      } 
marci@174
   622
    }
marci@174
   623
    Number flowValue() { 
marci@174
   624
      Number a=0;
marci@174
   625
      OutEdgeIt e;
marci@174
   626
      for(G->/*getF*/first(e, s); G->valid(e); G->next(e)) {
marci@174
   627
	a+=flow->get(e);
marci@174
   628
      }
marci@174
   629
      return a;
marci@174
   630
    }
marci@174
   631
  };
marci@174
   632
marci@193
   633
marci@193
   634
  template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
marci@193
   635
  class MaxMatching {
marci@193
   636
  public:
marci@193
   637
    typedef typename Graph::Node Node;
marci@194
   638
    typedef typename Graph::NodeIt NodeIt;
marci@193
   639
    typedef typename Graph::Edge Edge;
marci@193
   640
    typedef typename Graph::EdgeIt EdgeIt;
marci@193
   641
    typedef typename Graph::OutEdgeIt OutEdgeIt;
marci@193
   642
    typedef typename Graph::InEdgeIt InEdgeIt;
marci@193
   643
marci@193
   644
    typedef typename Graph::NodeMap<bool> SMap;
marci@193
   645
    typedef typename Graph::NodeMap<bool> TMap;
marci@193
   646
  private:
marci@193
   647
    const Graph* G;
marci@193
   648
    SMap* S;
marci@193
   649
    TMap* T;
marci@193
   650
    //Node s;
marci@193
   651
    //Node t;
marci@193
   652
    FlowMap* flow;
marci@193
   653
    const CapacityMap* capacity;
marci@193
   654
    typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
marci@193
   655
    typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
marci@193
   656
    typedef typename AugGraph::Edge AugEdge;
marci@198
   657
    typename Graph::NodeMap<int> used; //0
marci@193
   658
marci@193
   659
  public:
marci@193
   660
    MaxMatching(const Graph& _G, SMap& _S, TMap& _T, FlowMap& _flow, const CapacityMap& _capacity) : 
marci@198
   661
      G(&_G), S(&_S), T(&_T), flow(&_flow), capacity(&_capacity), used(_G) { }
marci@193
   662
    bool augmentOnShortestPath() {
marci@193
   663
      AugGraph res_graph(*G, *flow, *capacity);
marci@193
   664
      bool _augment=false;
marci@193
   665
      
marci@193
   666
      typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@193
   667
      BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
marci@193
   668
      typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@193
   669
      for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
marci@198
   670
	if ((S->get(s)) && (used.get(s)<1) ) {
marci@198
   671
	  //Number u=0;
marci@198
   672
	  //for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
marci@198
   673
	  //u+=flow->get(e);
marci@198
   674
	  //if (u<1) {
marci@196
   675
	    res_bfs.pushAndSetReached(s);
marci@196
   676
	    pred.set(s, AugEdge(INVALID));
marci@198
   677
	    //}
marci@193
   678
	}
marci@193
   679
      }
marci@193
   680
      
marci@193
   681
      typename AugGraph::NodeMap<Number> free(res_graph);
marci@193
   682
	
marci@193
   683
      Node n;
marci@193
   684
      //searching for augmenting path
marci@193
   685
      while ( !res_bfs.finished() ) { 
marci@193
   686
	AugOutEdgeIt e=res_bfs;
marci@193
   687
	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
marci@193
   688
	  Node v=res_graph.tail(e);
marci@193
   689
	  Node w=res_graph.head(e);
marci@193
   690
	  pred.set(w, e);
marci@193
   691
	  if (res_graph.valid(pred.get(v))) {
marci@193
   692
	    free.set(w, std::min(free.get(v), res_graph.free(e)));
marci@193
   693
	  } else {
marci@193
   694
	    free.set(w, res_graph.free(e)); 
marci@193
   695
	  }
marci@197
   696
	  n=res_graph.head(e);
marci@198
   697
	  if (T->get(n) && (used.get(n)<1) ) { 
marci@198
   698
	    //Number u=0;
marci@198
   699
	    //for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
marci@198
   700
	    //u+=flow->get(f);
marci@198
   701
	    //if (u<1) {
marci@197
   702
	      _augment=true; 
marci@197
   703
	      break; 
marci@198
   704
	      //}
marci@193
   705
	  }
marci@193
   706
	}
marci@193
   707
	
marci@193
   708
	++res_bfs;
marci@193
   709
      } //end of searching augmenting path
marci@193
   710
marci@193
   711
      if (_augment) {
marci@193
   712
	//Node n=t;
marci@198
   713
	used.set(n, 1); //mind2 vegen jav
marci@194
   714
	Number augment_value=free.get(n);
marci@193
   715
	while (res_graph.valid(pred.get(n))) { 
marci@193
   716
	  AugEdge e=pred.get(n);
marci@193
   717
	  res_graph.augment(e, augment_value); 
marci@193
   718
	  n=res_graph.tail(e);
marci@193
   719
	}
marci@198
   720
	used.set(n, 1); //mind2 vegen jav
marci@193
   721
      }
marci@193
   722
marci@193
   723
      return _augment;
marci@193
   724
    }
marci@193
   725
marci@193
   726
//     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
marci@193
   727
//       bool _augment=false;
marci@193
   728
marci@193
   729
//       AugGraph res_graph(*G, *flow, *capacity);
marci@193
   730
marci@193
   731
//       typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@193
   732
//       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
marci@193
   733
marci@197
   734
marci@197
   735
marci@197
   736
marci@197
   737
marci@197
   738
//       //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@197
   739
//       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
marci@197
   740
// 	if (S->get(s)) {
marci@197
   741
// 	  Number u=0;
marci@197
   742
// 	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
marci@197
   743
// 	    u+=flow->get(e);
marci@197
   744
// 	  if (u<1) {
marci@197
   745
// 	    res_bfs.pushAndSetReached(s);
marci@197
   746
// 	    //pred.set(s, AugEdge(INVALID));
marci@197
   747
// 	  }
marci@197
   748
// 	}
marci@197
   749
//       }
marci@197
   750
marci@197
   751
marci@197
   752
marci@197
   753
marci@197
   754
//       //bfs.pushAndSetReached(s);
marci@193
   755
//       typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
marci@193
   756
//       while ( !bfs.finished() ) { 
marci@193
   757
// 	AugOutEdgeIt e=bfs;
marci@193
   758
// 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
marci@193
   759
// 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
marci@193
   760
// 	}
marci@193
   761
	
marci@193
   762
// 	++bfs;
marci@193
   763
//       } //computing distances from s in the residual graph
marci@193
   764
marci@193
   765
//       MutableGraph F;
marci@193
   766
//       typename AugGraph::NodeMap<typename MutableGraph::Node> 
marci@193
   767
// 	res_graph_to_F(res_graph);
marci@193
   768
//       for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
marci@193
   769
// 	res_graph_to_F.set(n, F.addNode());
marci@193
   770
//       }
marci@193
   771
      
marci@193
   772
//       typename MutableGraph::Node sF=res_graph_to_F.get(s);
marci@193
   773
//       typename MutableGraph::Node tF=res_graph_to_F.get(t);
marci@193
   774
marci@193
   775
//       typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
marci@193
   776
//       typename MutableGraph::EdgeMap<Number> residual_capacity(F);
marci@193
   777
marci@193
   778
//       //Making F to the graph containing the edges of the residual graph 
marci@193
   779
//       //which are in some shortest paths
marci@193
   780
//       for(typename AugGraph::EdgeIt e=res_graph.template first<typename AugGraph::EdgeIt>(); res_graph.valid(e); res_graph.next(e)) {
marci@193
   781
// 	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
marci@193
   782
// 	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
marci@193
   783
// 	  original_edge.update();
marci@193
   784
// 	  original_edge.set(f, e);
marci@193
   785
// 	  residual_capacity.update();
marci@193
   786
// 	  residual_capacity.set(f, res_graph.free(e));
marci@193
   787
// 	} 
marci@193
   788
//       }
marci@193
   789
marci@193
   790
//       bool __augment=true;
marci@193
   791
marci@193
   792
//       while (__augment) {
marci@193
   793
// 	__augment=false;
marci@193
   794
// 	//computing blocking flow with dfs
marci@193
   795
// 	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
marci@193
   796
// 	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
marci@193
   797
// 	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
marci@193
   798
// 	pred.set(sF, typename MutableGraph::Edge(INVALID));
marci@193
   799
// 	//invalid iterators for sources
marci@193
   800
marci@193
   801
// 	typename MutableGraph::NodeMap<Number> free(F);
marci@193
   802
marci@193
   803
// 	dfs.pushAndSetReached(sF);      
marci@193
   804
// 	while (!dfs.finished()) {
marci@193
   805
// 	  ++dfs;
marci@193
   806
// 	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
marci@193
   807
// 	    if (dfs.isBNodeNewlyReached()) {
marci@193
   808
// 	      typename MutableGraph::Node v=F.aNode(dfs);
marci@193
   809
// 	      typename MutableGraph::Node w=F.bNode(dfs);
marci@193
   810
// 	      pred.set(w, dfs);
marci@193
   811
// 	      if (F.valid(pred.get(v))) {
marci@193
   812
// 		free.set(w, std::min(free.get(v), residual_capacity.get(dfs)));
marci@193
   813
// 	      } else {
marci@193
   814
// 		free.set(w, residual_capacity.get(dfs)); 
marci@193
   815
// 	      }
marci@193
   816
// 	      if (w==tF) { 
marci@193
   817
// 		__augment=true; 
marci@193
   818
// 		_augment=true;
marci@193
   819
// 		break; 
marci@193
   820
// 	      }
marci@193
   821
	      
marci@193
   822
// 	    } else {
marci@193
   823
// 	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
marci@193
   824
// 	    }
marci@193
   825
// 	  } 
marci@193
   826
// 	}
marci@193
   827
marci@193
   828
// 	if (__augment) {
marci@193
   829
// 	  typename MutableGraph::Node n=tF;
marci@193
   830
// 	  Number augment_value=free.get(tF);
marci@193
   831
// 	  while (F.valid(pred.get(n))) { 
marci@193
   832
// 	    typename MutableGraph::Edge e=pred.get(n);
marci@193
   833
// 	    res_graph.augment(original_edge.get(e), augment_value); 
marci@193
   834
// 	    n=F.tail(e);
marci@193
   835
// 	    if (residual_capacity.get(e)==augment_value) 
marci@193
   836
// 	      F.erase(e); 
marci@193
   837
// 	    else 
marci@193
   838
// 	      residual_capacity.set(e, residual_capacity.get(e)-augment_value);
marci@193
   839
// 	  }
marci@193
   840
// 	}
marci@193
   841
	
marci@193
   842
//       }
marci@193
   843
            
marci@193
   844
//       return _augment;
marci@193
   845
//     }
marci@197
   846
    bool augmentOnBlockingFlow2() {
marci@197
   847
      bool _augment=false;
marci@193
   848
marci@197
   849
      //typedef ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> EAugGraph;
marci@197
   850
      typedef FilterGraphWrapper< ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap> > EAugGraph;
marci@197
   851
      typedef typename EAugGraph::OutEdgeIt EAugOutEdgeIt;
marci@197
   852
      typedef typename EAugGraph::Edge EAugEdge;
marci@193
   853
marci@197
   854
      EAugGraph res_graph(*G, *flow, *capacity);
marci@193
   855
marci@197
   856
      //typedef typename EAugGraph::NodeMap<bool> ReachedMap;
marci@197
   857
      BfsIterator4< 
marci@197
   858
	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>, 
marci@197
   859
	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt, 
marci@197
   860
	ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::NodeMap<bool> > bfs(res_graph);
marci@197
   861
marci@197
   862
marci@197
   863
      //typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@197
   864
      for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
marci@197
   865
	if (S->get(s)) {
marci@197
   866
	  Number u=0;
marci@197
   867
	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
marci@197
   868
	    u+=flow->get(e);
marci@197
   869
	  if (u<1) {
marci@197
   870
	    bfs.pushAndSetReached(s);
marci@197
   871
	    //pred.set(s, AugEdge(INVALID));
marci@197
   872
	  }
marci@197
   873
	}
marci@197
   874
      }
marci@197
   875
marci@193
   876
      
marci@197
   877
      //bfs.pushAndSetReached(s);
marci@193
   878
marci@197
   879
      typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::
marci@197
   880
	NodeMap<int>& dist=res_graph.dist;
marci@193
   881
marci@197
   882
      while ( !bfs.finished() ) {
marci@197
   883
	typename ErasingResGraphWrapper<Graph, Number, FlowMap, CapacityMap>::OutEdgeIt e=bfs;
marci@197
   884
	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
marci@197
   885
	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
marci@197
   886
	}
marci@197
   887
	++bfs;	
marci@197
   888
      } //computing distances from s in the residual graph
marci@193
   889
marci@197
   890
      bool __augment=true;
marci@193
   891
marci@197
   892
      while (__augment) {
marci@193
   893
marci@197
   894
	__augment=false;
marci@197
   895
	//computing blocking flow with dfs
marci@197
   896
	typedef typename EAugGraph::NodeMap<bool> BlockingReachedMap;
marci@197
   897
	DfsIterator4< EAugGraph, EAugOutEdgeIt, BlockingReachedMap > 
marci@197
   898
	  dfs(res_graph);
marci@197
   899
	typename EAugGraph::NodeMap<EAugEdge> pred(res_graph, INVALID); 
marci@197
   900
	//pred.set(s, EAugEdge(INVALID));
marci@197
   901
	//invalid iterators for sources
marci@193
   902
marci@197
   903
	typename EAugGraph::NodeMap<Number> free(res_graph);
marci@193
   904
marci@197
   905
marci@197
   906
	//typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@197
   907
      for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
marci@197
   908
	if (S->get(s)) {
marci@197
   909
	  Number u=0;
marci@197
   910
	  for(OutEdgeIt e=G->template first<OutEdgeIt>(s); G->valid(e); G->next(e))
marci@197
   911
	    u+=flow->get(e);
marci@197
   912
	  if (u<1) {
marci@197
   913
	    dfs.pushAndSetReached(s);
marci@197
   914
	    //pred.set(s, AugEdge(INVALID));
marci@197
   915
	  }
marci@197
   916
	}
marci@197
   917
      }
marci@197
   918
marci@197
   919
marci@197
   920
marci@197
   921
      //dfs.pushAndSetReached(s);
marci@197
   922
      typename EAugGraph::Node n;
marci@197
   923
	while (!dfs.finished()) {
marci@197
   924
	  ++dfs;
marci@197
   925
	  if (res_graph.valid(EAugOutEdgeIt(dfs))) { 
marci@197
   926
	    if (dfs.isBNodeNewlyReached()) {
marci@193
   927
	  
marci@197
   928
	      typename EAugGraph::Node v=res_graph.aNode(dfs);
marci@197
   929
	      typename EAugGraph::Node w=res_graph.bNode(dfs);
marci@193
   930
marci@197
   931
	      pred.set(w, EAugOutEdgeIt(dfs));
marci@197
   932
	      if (res_graph.valid(pred.get(v))) {
marci@197
   933
		free.set(w, std::min(free.get(v), res_graph.free(dfs)));
marci@197
   934
	      } else {
marci@197
   935
		free.set(w, res_graph.free(dfs)); 
marci@197
   936
	      }
marci@197
   937
	     
marci@197
   938
	      n=w;
marci@197
   939
	      if (T->get(w)) {
marci@197
   940
		Number u=0;
marci@197
   941
		for(InEdgeIt f=G->template first<InEdgeIt>(n); G->valid(f); G->next(f))
marci@197
   942
		  u+=flow->get(f);
marci@197
   943
		if (u<1) {
marci@197
   944
		  __augment=true; 
marci@197
   945
		  _augment=true;
marci@197
   946
		  break; 
marci@197
   947
		}
marci@197
   948
	      }
marci@197
   949
	    } else {
marci@197
   950
	      res_graph.erase(dfs);
marci@197
   951
	    }
marci@197
   952
	  } 
marci@193
   953
marci@197
   954
	}
marci@193
   955
marci@197
   956
	if (__augment) {
marci@197
   957
	  // typename EAugGraph::Node n=t;
marci@197
   958
	  Number augment_value=free.get(n);
marci@197
   959
	  while (res_graph.valid(pred.get(n))) { 
marci@197
   960
	    EAugEdge e=pred.get(n);
marci@197
   961
	    res_graph.augment(e, augment_value);
marci@197
   962
	    n=res_graph.tail(e);
marci@197
   963
	    if (res_graph.free(e)==0)
marci@197
   964
	      res_graph.erase(e);
marci@197
   965
	  }
marci@197
   966
	}
marci@193
   967
      
marci@197
   968
      }
marci@193
   969
            
marci@197
   970
      return _augment;
marci@197
   971
    }
marci@193
   972
    void run() {
marci@193
   973
      //int num_of_augmentations=0;
marci@193
   974
      while (augmentOnShortestPath()) { 
marci@193
   975
	//while (augmentOnBlockingFlow<MutableGraph>()) { 
marci@193
   976
	//std::cout << ++num_of_augmentations << " ";
marci@193
   977
	//std::cout<<std::endl;
marci@193
   978
      } 
marci@193
   979
    }
marci@194
   980
//     template<typename MutableGraph> void run() {
marci@194
   981
//       //int num_of_augmentations=0;
marci@194
   982
//       //while (augmentOnShortestPath()) { 
marci@194
   983
// 	while (augmentOnBlockingFlow<MutableGraph>()) { 
marci@194
   984
// 	//std::cout << ++num_of_augmentations << " ";
marci@194
   985
// 	//std::cout<<std::endl;
marci@194
   986
//       } 
marci@194
   987
//     } 
marci@193
   988
    Number flowValue() { 
marci@193
   989
      Number a=0;
marci@194
   990
      EdgeIt e;
marci@194
   991
      for(G->/*getF*/first(e); G->valid(e); G->next(e)) {
marci@193
   992
	a+=flow->get(e);
marci@193
   993
      }
marci@193
   994
      return a;
marci@193
   995
    }
marci@193
   996
  };
marci@193
   997
marci@193
   998
marci@193
   999
marci@193
  1000
marci@193
  1001
marci@174
  1002
  
marci@174
  1003
//   template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
marci@174
  1004
//   class MaxFlow2 {
marci@174
  1005
//   public:
marci@174
  1006
//     typedef typename Graph::Node Node;
marci@174
  1007
//     typedef typename Graph::Edge Edge;
marci@174
  1008
//     typedef typename Graph::EdgeIt EdgeIt;
marci@174
  1009
//     typedef typename Graph::OutEdgeIt OutEdgeIt;
marci@174
  1010
//     typedef typename Graph::InEdgeIt InEdgeIt;
marci@174
  1011
//   private:
marci@174
  1012
//     const Graph& G;
marci@174
  1013
//     std::list<Node>& S;
marci@174
  1014
//     std::list<Node>& T;
marci@174
  1015
//     FlowMap& flow;
marci@174
  1016
//     const CapacityMap& capacity;
marci@174
  1017
//     typedef ResGraphWrapper<Graph, Number, FlowMap, CapacityMap > AugGraph;
marci@174
  1018
//     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
marci@174
  1019
//     typedef typename AugGraph::Edge AugEdge;
marci@174
  1020
//     typename Graph::NodeMap<bool> SMap;
marci@174
  1021
//     typename Graph::NodeMap<bool> TMap;
marci@174
  1022
//   public:
marci@174
  1023
//     MaxFlow2(const Graph& _G, std::list<Node>& _S, std::list<Node>& _T, FlowMap& _flow, const CapacityMap& _capacity) : G(_G), S(_S), T(_T), flow(_flow), capacity(_capacity), SMap(_G), TMap(_G) { 
marci@174
  1024
//       for(typename std::list<Node>::const_iterator i=S.begin(); 
marci@174
  1025
// 	  i!=S.end(); ++i) { 
marci@174
  1026
// 	SMap.set(*i, true); 
marci@174
  1027
//       }
marci@174
  1028
//       for (typename std::list<Node>::const_iterator i=T.begin(); 
marci@174
  1029
// 	   i!=T.end(); ++i) { 
marci@174
  1030
// 	TMap.set(*i, true); 
marci@174
  1031
//       }
marci@174
  1032
//     }
marci@174
  1033
//     bool augment() {
marci@174
  1034
//       AugGraph res_graph(G, flow, capacity);
marci@174
  1035
//       bool _augment=false;
marci@174
  1036
//       Node reached_t_node;
marci@174
  1037
      
marci@174
  1038
//       typedef typename AugGraph::NodeMap<bool> ReachedMap;
marci@174
  1039
//       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
marci@174
  1040
//       for(typename std::list<Node>::const_iterator i=S.begin(); 
marci@174
  1041
// 	  i!=S.end(); ++i) {
marci@174
  1042
// 	res_bfs.pushAndSetReached(*i);
marci@174
  1043
//       }
marci@174
  1044
//       //res_bfs.pushAndSetReached(s);
marci@174
  1045
	
marci@174
  1046
//       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
marci@174
  1047
//       //filled up with invalid iterators
marci@174
  1048
      
marci@174
  1049
//       typename AugGraph::NodeMap<Number> free(res_graph);
marci@174
  1050
	
marci@174
  1051
//       //searching for augmenting path
marci@174
  1052
//       while ( !res_bfs.finished() ) { 
marci@174
  1053
// 	AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
marci@174
  1054
// 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
marci@174
  1055
// 	  Node v=res_graph.tail(e);
marci@174
  1056
// 	  Node w=res_graph.head(e);
marci@174
  1057
// 	  pred.set(w, e);
marci@174
  1058
// 	  if (pred.get(v).valid()) {
marci@174
  1059
// 	    free.set(w, std::min(free.get(v), e.free()));
marci@174
  1060
// 	  } else {
marci@174
  1061
// 	    free.set(w, e.free()); 
marci@174
  1062
// 	  }
marci@174
  1063
// 	  if (TMap.get(res_graph.head(e))) { 
marci@174
  1064
// 	    _augment=true; 
marci@174
  1065
// 	    reached_t_node=res_graph.head(e);
marci@174
  1066
// 	    break; 
marci@174
  1067
// 	  }
marci@174
  1068
// 	}
marci@174
  1069
	
marci@174
  1070
// 	++res_bfs;
marci@174
  1071
//       } //end of searching augmenting path
marci@174
  1072
marci@174
  1073
//       if (_augment) {
marci@174
  1074
// 	Node n=reached_t_node;
marci@174
  1075
// 	Number augment_value=free.get(reached_t_node);
marci@174
  1076
// 	while (pred.get(n).valid()) { 
marci@174
  1077
// 	  AugEdge e=pred.get(n);
marci@174
  1078
// 	  e.augment(augment_value); 
marci@174
  1079
// 	  n=res_graph.tail(e);
marci@174
  1080
// 	}
marci@174
  1081
//       }
marci@174
  1082
marci@174
  1083
//       return _augment;
marci@174
  1084
//     }
marci@174
  1085
//     void run() {
marci@174
  1086
//       while (augment()) { } 
marci@174
  1087
//     }
marci@174
  1088
//     Number flowValue() { 
marci@174
  1089
//       Number a=0;
marci@174
  1090
//       for(typename std::list<Node>::const_iterator i=S.begin(); 
marci@174
  1091
// 	  i!=S.end(); ++i) { 
marci@174
  1092
// 	for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
marci@174
  1093
// 	  a+=flow.get(e);
marci@174
  1094
// 	}
marci@174
  1095
// 	for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
marci@174
  1096
// 	  a-=flow.get(e);
marci@174
  1097
// 	}
marci@174
  1098
//       }
marci@174
  1099
//       return a;
marci@174
  1100
//     }
marci@174
  1101
//   };
marci@174
  1102
marci@174
  1103
marci@174
  1104
marci@174
  1105
} // namespace hugo
marci@174
  1106
marci@174
  1107
#endif //EDMONDS_KARP_H