COIN-OR::LEMON - Graph Library

Changeset 862:732f2acb7239 in lemon-0.x


Ignore:
Timestamp:
09/16/04 12:59:52 (20 years ago)
Author:
marci
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@1162
Message:

bug correction

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/work/marci/augmenting_flow.h

    r854 r862  
    44
    55#include <vector>
    6 #include <queue>
    7 #include <stack>
     6//#include <queue>
     7//#include <stack>
    88#include <iostream>
    99
     
    1212#include <hugo/invalid.h>
    1313#include <hugo/maps.h>
    14 //#include <for_each_macros.h>
    1514
    1615/// \file
     
    2019namespace hugo {
    2120
     21  /// \brief A map for filtering the edge-set to those edges
     22  /// which are tight w.r.t. some node_potential map and
     23  /// edge_distance map.
     24  ///
     25  /// A node-map node_potential is said to be a potential w.r.t.
     26  /// an edge-map edge_distance
     27  /// if and only if for each edge e, node_potential[g.head(e)]
     28  /// <= edge_distance[e]+node_potential[g.tail(e)]
     29  /// (or the reverse inequality holds for each edge).
     30  /// An edge is said to be tight if this inequality holds with equality,
     31  /// and the map returns true exactly for those edges.
     32  /// To avoid rounding errors, it is recommended to use this class with exact
     33  /// types, e.g. with int.
     34  template<typename Graph,
     35           typename NodePotentialMap, typename EdgeDistanceMap>
     36  class TightEdgeFilterMap {
     37  protected:
     38    const Graph* g;
     39    NodePotentialMap* node_potential;
     40    EdgeDistanceMap* edge_distance;
     41  public:
     42    TightEdgeFilterMap(Graph& _g, NodePotentialMap& _node_potential,
     43                       EdgeDistanceMap& _edge_distance) :
     44      g(&_g), node_potential(&_node_potential),
     45      edge_distance(&_edge_distance) { }
     46//     void set(const typename Graph::Node& n, int a) {
     47//       pot->set(n, a);
     48//     }
     49//     int operator[](const typename Graph::Node& n) const {
     50//       return (*node_potential)[n];
     51//     }
     52    bool operator[](const typename Graph::Edge& e) const {
     53      return ((*node_potential)[g->head(e)] ==
     54              (*edge_distance)[e]+(*node_potential)[g->tail(e)]);
     55    }
     56  };
     57
    2258  /// \addtogroup galgs
    2359  /// @{                                                                                                                                       
    24   ///Maximum flow algorithms class.
    25 
    26   ///This class provides various algorithms for finding a flow of
    27   ///maximum value in a directed graph. The \e source node, the \e
    28   ///target node, the \e capacity of the edges and the \e starting \e
    29   ///flow value of the edges should be passed to the algorithm through the
    30   ///constructor. It is possible to change these quantities using the
    31   ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
    32   ///\ref resetFlow. Before any subsequent runs of any algorithm of
    33   ///the class \ref resetFlow should be called.
    34 
    35   ///After running an algorithm of the class, the actual flow value
    36   ///can be obtained by calling \ref flowValue(). The minimum
    37   ///value cut can be written into a \c node map of \c bools by
    38   ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
    39   ///the inclusionwise minimum and maximum of the minimum value
    40   ///cuts, resp.)                                                                                                                               
     60  /// Class for augmenting path flow algorithms.
     61
     62  /// This class provides various algorithms for finding a flow of
     63  /// maximum value in a directed graph. The \e source node, the \e
     64  /// target node, the \e capacity of the edges and the \e starting \e
     65  /// flow value of the edges should be passed to the algorithm through the
     66  /// constructor.
     67//   /// It is possible to change these quantities using the
     68//   /// functions \ref resetSource, \ref resetTarget, \ref resetCap and
     69//   /// \ref resetFlow. Before any subsequent runs of any algorithm of
     70//   /// the class \ref resetFlow should be called.
     71
     72  /// After running an algorithm of the class, the actual flow value
     73  /// can be obtained by calling \ref flowValue(). The minimum
     74  /// value cut can be written into a \c node map of \c bools by
     75  /// calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
     76  /// the inclusionwise minimum and maximum of the minimum value
     77  /// cuts, resp.)                                                                                                                               
    4178  ///\param Graph The directed graph type the algorithm runs on.
    4279  ///\param Num The number type of the capacities and the flow values.
    4380  ///\param CapMap The capacity map type.
    4481  ///\param FlowMap The flow map type.                                                                                                           
    45   ///\author Marton Makai, Jacint Szabo
    46 //   template <typename Graph, typename Num,
    47 //          typename CapMap=typename Graph::template EdgeMap<Num>,
    48 //             typename FlowMap=typename Graph::template EdgeMap<Num> >
    49 //   class MaxFlow {
    50 //   protected:
    51 //     typedef typename Graph::Node Node;
    52 //     typedef typename Graph::NodeIt NodeIt;
    53 //     typedef typename Graph::EdgeIt EdgeIt;
    54 //     typedef typename Graph::OutEdgeIt OutEdgeIt;
    55 //     typedef typename Graph::InEdgeIt InEdgeIt;
    56 
    57 //     typedef typename std::vector<std::stack<Node> > VecStack;
    58 //     typedef typename Graph::template NodeMap<Node> NNMap;
    59 //     typedef typename std::vector<Node> VecNode;
    60 
    61 //     const Graph* g;
    62 //     Node s;
    63 //     Node t;
    64 //     const CapMap* capacity;
    65 //     FlowMap* flow;
    66 //     int n;      //the number of nodes of G
    67 //     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    68 //     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    69 //     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    70 //     typedef typename ResGW::Edge ResGWEdge;
    71 //     //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    72 //     typedef typename Graph::template NodeMap<int> ReachedMap;
    73 
    74 
    75 //     //level works as a bool map in augmenting path algorithms and is
    76 //     //used by bfs for storing reached information.  In preflow, it
    77 //     //shows the levels of nodes.     
    78 //     ReachedMap level;
    79 
    80 //     //excess is needed only in preflow
    81 //     typename Graph::template NodeMap<Num> excess;
    82 
    83 //     //fixme   
    84 // //   protected:
    85 //     //     MaxFlow() { }
    86 //     //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
    87 //     //            FlowMap& _flow)
    88 //     //       {
    89 //     //       g=&_G;
    90 //     //       s=_s;
    91 //     //       t=_t;
    92 //     //       capacity=&_capacity;
    93 //     //       flow=&_flow;
    94 //     //       n=_G.nodeNum;
    95 //     //       level.set (_G); //kellene vmi ilyesmi fv
    96 //     //       excess(_G,0); //itt is
    97 //     //       }
    98 
    99 //     // constants used for heuristics
    100 //     static const int H0=20;
    101 //     static const int H1=1;
    102 
    103 //   public:
    104 
    105 //     ///Indicates the property of the starting flow.
    106 
    107 //     ///Indicates the property of the starting flow. The meanings are as follows:
    108 //     ///- \c ZERO_FLOW: constant zero flow
    109 //     ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
    110 //     ///the sum of the out-flows in every node except the \e source and
    111 //     ///the \e target.
    112 //     ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
    113 //     ///least the sum of the out-flows in every node except the \e source.
    114 //     ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
    115 //     ///set to the constant zero flow in the beginning of the algorithm in this case.
    116 //     enum FlowEnum{
    117 //       ZERO_FLOW,
    118 //       GEN_FLOW,
    119 //       PRE_FLOW,
    120 //       NO_FLOW
    121 //     };
    122 
    123 //     enum StatusEnum {
    124 //       AFTER_NOTHING,
    125 //       AFTER_AUGMENTING,
    126 //       AFTER_FAST_AUGMENTING,
    127 //       AFTER_PRE_FLOW_PHASE_1,     
    128 //       AFTER_PRE_FLOW_PHASE_2
    129 //     };
    130 
    131 //     /// Don not needle this flag only if necessary.
    132 //     StatusEnum status;
    133 // //     int number_of_augmentations;
    134 
    135 
    136 // //     template<typename IntMap>
    137 // //     class TrickyReachedMap {
    138 // //     protected:
    139 // //       IntMap* map;
    140 // //       int* number_of_augmentations;
    141 // //     public:
    142 // //       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
    143 // //   map(&_map), number_of_augmentations(&_number_of_augmentations) { }
    144 // //       void set(const Node& n, bool b) {
    145 // //   if (b)
    146 // //     map->set(n, *number_of_augmentations);
    147 // //   else
    148 // //     map->set(n, *number_of_augmentations-1);
    149 // //       }
    150 // //       bool operator[](const Node& n) const {
    151 // //   return (*map)[n]==*number_of_augmentations;
    152 // //       }
    153 // //     };
    154    
    155 //     MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
    156 //          FlowMap& _flow) :
    157 //       g(&_G), s(_s), t(_t), capacity(&_capacity),
    158 //       flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
    159 //       status(AFTER_NOTHING) { }
    160 
    161 //     ///Runs a maximum flow algorithm.
    162 
    163 //     ///Runs a preflow algorithm, which is the fastest maximum flow
    164 //     ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
    165 //     ///\pre The starting flow must be
    166 //     /// - a constant zero flow if \c fe is \c ZERO_FLOW,
    167 //     /// - an arbitary flow if \c fe is \c GEN_FLOW,
    168 //     /// - an arbitary preflow if \c fe is \c PRE_FLOW,
    169 //     /// - any map if \c fe is NO_FLOW.
    170 //     void run(FlowEnum fe=ZERO_FLOW) {
    171 //       preflow(fe);
    172 //     }
    173 
    174                                                                              
    175 //     ///Runs a preflow algorithm. 
    176 
    177 //     ///Runs a preflow algorithm. The preflow algorithms provide the
    178 //     ///fastest way to compute a maximum flow in a directed graph.
    179 //     ///\pre The starting flow must be
    180 //     /// - a constant zero flow if \c fe is \c ZERO_FLOW,
    181 //     /// - an arbitary flow if \c fe is \c GEN_FLOW,
    182 //     /// - an arbitary preflow if \c fe is \c PRE_FLOW,
    183 //     /// - any map if \c fe is NO_FLOW.
    184 //     void preflow(FlowEnum fe) {
    185 //       preflowPhase1(fe);
    186 //       preflowPhase2();
    187 //     }
    188 //     // Heuristics:
    189 //     //   2 phase
    190 //     //   gap
    191 //     //   list 'level_list' on the nodes on level i implemented by hand
    192 //     //   stack 'active' on the active nodes on level i                                                                                   
    193 //     //   runs heuristic 'highest label' for H1*n relabels
    194 //     //   runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    195 //     //   Parameters H0 and H1 are initialized to 20 and 1.
    196 
    197 //     ///Runs the first phase of the preflow algorithm.
    198 
    199 //     ///The preflow algorithm consists of two phases, this method runs the
    200 //     ///first phase. After the first phase the maximum flow value and a
    201 //     ///minimum value cut can already be computed, though a maximum flow
    202 //     ///is net yet obtained. So after calling this method \ref flowValue
    203 //     ///and \ref actMinCut gives proper results.
    204 //     ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
    205 //     ///give minimum value cuts unless calling \ref preflowPhase2.
    206 //     ///\pre The starting flow must be
    207 //     /// - a constant zero flow if \c fe is \c ZERO_FLOW,
    208 //     /// - an arbitary flow if \c fe is \c GEN_FLOW,
    209 //     /// - an arbitary preflow if \c fe is \c PRE_FLOW,
    210 //     /// - any map if \c fe is NO_FLOW.
    211 //     void preflowPhase1(FlowEnum fe);
    212 
    213 //     ///Runs the second phase of the preflow algorithm.
    214 
    215 //     ///The preflow algorithm consists of two phases, this method runs
    216 //     ///the second phase. After calling \ref preflowPhase1 and then
    217 //     ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
    218 //     ///\ref minMinCut and \ref maxMinCut give proper results.
    219 //     ///\pre \ref preflowPhase1 must be called before.
    220 //     void preflowPhase2();
    221 
    222 //     /// Returns the maximum value of a flow.
    223 
    224 //     /// Returns the maximum value of a flow, by counting the
    225 //     /// over-flow of the target node \ref t.
    226 //     /// It can be called already after running \ref preflowPhase1.
    227 //     Num flowValue() const {
    228 //       Num a=0;
    229 //       FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
    230 //       FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
    231 //       return a;
    232 //       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
    233 //     }
    234 
    235 //     ///Returns a minimum value cut after calling \ref preflowPhase1.
    236 
    237 //     ///After the first phase of the preflow algorithm the maximum flow
    238 //     ///value and a minimum value cut can already be computed. This
    239 //     ///method can be called after running \ref preflowPhase1 for
    240 //     ///obtaining a minimum value cut.
    241 //     /// \warning Gives proper result only right after calling \ref
    242 //     /// preflowPhase1.
    243 //     /// \todo We have to make some status variable which shows the
    244 //     /// actual state
    245 //     /// of the class. This enables us to determine which methods are valid
    246 //     /// for MinCut computation
    247 //     template<typename _CutMap>
    248 //     void actMinCut(_CutMap& M) const {
    249 //       NodeIt v;
    250 //       switch (status) {
    251 //       case AFTER_PRE_FLOW_PHASE_1:
    252 //      for(g->first(v); g->valid(v); g->next(v)) {
    253 //        if (level[v] < n) {
    254 //          M.set(v, false);
    255 //        } else {
    256 //          M.set(v, true);
    257 //        }
    258 //      }
    259 //      break;
    260 //       case AFTER_PRE_FLOW_PHASE_2:
    261 //       case AFTER_NOTHING:
    262 //       case AFTER_AUGMENTING:
    263 //       case AFTER_FAST_AUGMENTING:
    264 //      minMinCut(M);
    265 //      break;
    266 // //       case AFTER_AUGMENTING:
    267 // //   for(g->first(v); g->valid(v); g->next(v)) {
    268 // //     if (level[v]) {
    269 // //       M.set(v, true);
    270 // //     } else {
    271 // //       M.set(v, false);
    272 // //     }
    273 // //   }
    274 // //   break;
    275 // //       case AFTER_FAST_AUGMENTING:
    276 // //   for(g->first(v); g->valid(v); g->next(v)) {
    277 // //     if (level[v]==number_of_augmentations) {
    278 // //       M.set(v, true);
    279 // //     } else {
    280 // //       M.set(v, false);
    281 // //     }
    282 // //   }
    283 // //   break;
    284 //       }
    285 //     }
    286 
    287 //     ///Returns the inclusionwise minimum of the minimum value cuts.
    288 
    289 //     ///Sets \c M to the characteristic vector of the minimum value cut
    290 //     ///which is inclusionwise minimum. It is computed by processing
    291 //     ///a bfs from the source node \c s in the residual graph.
    292 //     ///\pre M should be a node map of bools initialized to false.
    293 //     ///\pre \c flow must be a maximum flow.
    294 //     template<typename _CutMap>
    295 //     void minMinCut(_CutMap& M) const {
    296 //       std::queue<Node> queue;
    297 
    298 //       M.set(s,true);
    299 //       queue.push(s);
    300 
    301 //       while (!queue.empty()) {
    302 //         Node w=queue.front();
    303 //      queue.pop();
    304 
    305 //      OutEdgeIt e;
    306 //      for(g->first(e,w) ; g->valid(e); g->next(e)) {
    307 //        Node v=g->head(e);
    308 //        if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
    309 //          queue.push(v);
    310 //          M.set(v, true);
    311 //        }
    312 //      }
    313 
    314 //      InEdgeIt f;
    315 //      for(g->first(f,w) ; g->valid(f); g->next(f)) {
    316 //        Node v=g->tail(f);
    317 //        if (!M[v] && (*flow)[f] > 0 ) {
    318 //          queue.push(v);
    319 //          M.set(v, true);
    320 //        }
    321 //      }
    322 //       }
    323 //     }
    324 
    325 //     ///Returns the inclusionwise maximum of the minimum value cuts.
    326 
    327 //     ///Sets \c M to the characteristic vector of the minimum value cut
    328 //     ///which is inclusionwise maximum. It is computed by processing a
    329 //     ///backward bfs from the target node \c t in the residual graph.
    330 //     ///\pre M should be a node map of bools initialized to false.
    331 //     ///\pre \c flow must be a maximum flow.
    332 //     template<typename _CutMap>
    333 //     void maxMinCut(_CutMap& M) const {
    334 
    335 //       NodeIt v;
    336 //       for(g->first(v) ; g->valid(v); g->next(v)) {
    337 //      M.set(v, true);
    338 //       }
    339 
    340 //       std::queue<Node> queue;
    341 
    342 //       M.set(t,false);
    343 //       queue.push(t);
    344 
    345 //       while (!queue.empty()) {
    346 //         Node w=queue.front();
    347 //      queue.pop();
    348 
    349 //      InEdgeIt e;
    350 //      for(g->first(e,w) ; g->valid(e); g->next(e)) {
    351 //        Node v=g->tail(e);
    352 //        if (M[v] && (*flow)[e] < (*capacity)[e] ) {
    353 //          queue.push(v);
    354 //          M.set(v, false);
    355 //        }
    356 //      }
    357 
    358 //      OutEdgeIt f;
    359 //      for(g->first(f,w) ; g->valid(f); g->next(f)) {
    360 //        Node v=g->head(f);
    361 //        if (M[v] && (*flow)[f] > 0 ) {
    362 //          queue.push(v);
    363 //          M.set(v, false);
    364 //        }
    365 //      }
    366 //       }
    367 //     }
    368 
    369 //     ///Returns a minimum value cut.
    370 
    371 //     ///Sets \c M to the characteristic vector of a minimum value cut.
    372 //     ///\pre M should be a node map of bools initialized to false.
    373 //     ///\pre \c flow must be a maximum flow.   
    374 //     template<typename CutMap>
    375 //     void minCut(CutMap& M) const { minMinCut(M); }
    376 
    377 //     ///Resets the source node to \c _s.
    378 
    379 //     ///Resets the source node to \c _s.
    380 //     ///
    381 //     void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
    382 
    383 //     ///Resets the target node to \c _t.
    384 
    385 //     ///Resets the target node to \c _t.
    386 //     ///
    387 //     void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
    388 
    389 //     /// Resets the edge map of the capacities to _cap.
    390 
    391 //     /// Resets the edge map of the capacities to _cap.
    392 //     ///
    393 //     void resetCap(const CapMap& _cap) { capacity=&_cap; status=AFTER_NOTHING; }
    394 
    395 //     /// Resets the edge map of the flows to _flow.
    396 
    397 //     /// Resets the edge map of the flows to _flow.
    398 //     ///
    399 //     void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
    400 
    401 
    402 //   private:
    403 
    404 //     int push(Node w, VecStack& active) {
    405 
    406 //       int lev=level[w];
    407 //       Num exc=excess[w];
    408 //       int newlevel=n;       //bound on the next level of w
    409 
    410 //       OutEdgeIt e;
    411 //       for(g->first(e,w); g->valid(e); g->next(e)) {
    412 
    413 //      if ( (*flow)[e] >= (*capacity)[e] ) continue;
    414 //      Node v=g->head(e);
    415 
    416 //      if( lev > level[v] ) { //Push is allowed now
    417 
    418 //        if ( excess[v]<=0 && v!=t && v!=s ) {
    419 //          int lev_v=level[v];
    420 //          active[lev_v].push(v);
    421 //        }
    422 
    423 //        Num cap=(*capacity)[e];
    424 //        Num flo=(*flow)[e];
    425 //        Num remcap=cap-flo;
    426 
    427 //        if ( remcap >= exc ) { //A nonsaturating push.
    428 
    429 //          flow->set(e, flo+exc);
    430 //          excess.set(v, excess[v]+exc);
    431 //          exc=0;
    432 //          break;
    433 
    434 //        } else { //A saturating push.
    435 //          flow->set(e, cap);
    436 //          excess.set(v, excess[v]+remcap);
    437 //          exc-=remcap;
    438 //        }
    439 //      } else if ( newlevel > level[v] ) newlevel = level[v];
    440 //       } //for out edges wv
    441 
    442 //       if ( exc > 0 ) {
    443 //      InEdgeIt e;
    444 //      for(g->first(e,w); g->valid(e); g->next(e)) {
    445 
    446 //        if( (*flow)[e] <= 0 ) continue;
    447 //        Node v=g->tail(e);
    448 
    449 //        if( lev > level[v] ) { //Push is allowed now
    450 
    451 //          if ( excess[v]<=0 && v!=t && v!=s ) {
    452 //            int lev_v=level[v];
    453 //            active[lev_v].push(v);
    454 //          }
    455 
    456 //          Num flo=(*flow)[e];
    457 
    458 //          if ( flo >= exc ) { //A nonsaturating push.
    459 
    460 //            flow->set(e, flo-exc);
    461 //            excess.set(v, excess[v]+exc);
    462 //            exc=0;
    463 //            break;
    464 //          } else {  //A saturating push.
    465 
    466 //            excess.set(v, excess[v]+flo);
    467 //            exc-=flo;
    468 //            flow->set(e,0);
    469 //          }
    470 //        } else if ( newlevel > level[v] ) newlevel = level[v];
    471 //      } //for in edges vw
    472 
    473 //       } // if w still has excess after the out edge for cycle
    474 
    475 //       excess.set(w, exc);
    476 
    477 //       return newlevel;
    478 //     }
    479 
    480 
    481 //     void preflowPreproc(FlowEnum fe, VecStack& active,
    482 //                      VecNode& level_list, NNMap& left, NNMap& right)
    483 //     {
    484 //       std::queue<Node> bfs_queue;
    485 
    486 //       switch (fe) {
    487 //       case NO_FLOW:   //flow is already set to const zero in this case
    488 //       case ZERO_FLOW:
    489 //      {
    490 //        //Reverse_bfs from t, to find the starting level.
    491 //        level.set(t,0);
    492 //        bfs_queue.push(t);
    493 
    494 //        while (!bfs_queue.empty()) {
    495 
    496 //          Node v=bfs_queue.front();
    497 //          bfs_queue.pop();
    498 //          int l=level[v]+1;
    499 
    500 //          InEdgeIt e;
    501 //          for(g->first(e,v); g->valid(e); g->next(e)) {
    502 //            Node w=g->tail(e);
    503 //            if ( level[w] == n && w != s ) {
    504 //              bfs_queue.push(w);
    505 //              Node first=level_list[l];
    506 //              if ( g->valid(first) ) left.set(first,w);
    507 //              right.set(w,first);
    508 //              level_list[l]=w;
    509 //              level.set(w, l);
    510 //            }
    511 //          }
    512 //        }
    513 
    514 //        //the starting flow
    515 //        OutEdgeIt e;
    516 //        for(g->first(e,s); g->valid(e); g->next(e))
    517 //          {
    518 //            Num c=(*capacity)[e];
    519 //            if ( c <= 0 ) continue;
    520 //            Node w=g->head(e);
    521 //            if ( level[w] < n ) {
    522 //              if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    523 //              flow->set(e, c);
    524 //              excess.set(w, excess[w]+c);
    525 //            }
    526 //          }
    527 //        break;
    528 //      }
    529 
    530 //       case GEN_FLOW:
    531 //       case PRE_FLOW:
    532 //      {
    533 //        //Reverse_bfs from t in the residual graph,
    534 //        //to find the starting level.
    535 //        level.set(t,0);
    536 //        bfs_queue.push(t);
    537 
    538 //        while (!bfs_queue.empty()) {
    539 
    540 //          Node v=bfs_queue.front();
    541 //          bfs_queue.pop();
    542 //          int l=level[v]+1;
    543 
    544 //          InEdgeIt e;
    545 //          for(g->first(e,v); g->valid(e); g->next(e)) {
    546 //            if ( (*capacity)[e] <= (*flow)[e] ) continue;
    547 //            Node w=g->tail(e);
    548 //            if ( level[w] == n && w != s ) {
    549 //              bfs_queue.push(w);
    550 //              Node first=level_list[l];
    551 //              if ( g->valid(first) ) left.set(first,w);
    552 //              right.set(w,first);
    553 //              level_list[l]=w;
    554 //              level.set(w, l);
    555 //            }
    556 //          }
    557 
    558 //          OutEdgeIt f;
    559 //          for(g->first(f,v); g->valid(f); g->next(f)) {
    560 //            if ( 0 >= (*flow)[f] ) continue;
    561 //            Node w=g->head(f);
    562 //            if ( level[w] == n && w != s ) {
    563 //              bfs_queue.push(w);
    564 //              Node first=level_list[l];
    565 //              if ( g->valid(first) ) left.set(first,w);
    566 //              right.set(w,first);
    567 //              level_list[l]=w;
    568 //              level.set(w, l);
    569 //            }
    570 //          }
    571 //        }
    572 
    573 
    574 //        //the starting flow
    575 //        OutEdgeIt e;
    576 //        for(g->first(e,s); g->valid(e); g->next(e))
    577 //          {
    578 //            Num rem=(*capacity)[e]-(*flow)[e];
    579 //            if ( rem <= 0 ) continue;
    580 //            Node w=g->head(e);
    581 //            if ( level[w] < n ) {
    582 //              if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    583 //              flow->set(e, (*capacity)[e]);
    584 //              excess.set(w, excess[w]+rem);
    585 //            }
    586 //          }
    587 
    588 //        InEdgeIt f;
    589 //        for(g->first(f,s); g->valid(f); g->next(f))
    590 //          {
    591 //            if ( (*flow)[f] <= 0 ) continue;
    592 //            Node w=g->tail(f);
    593 //            if ( level[w] < n ) {
    594 //              if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
    595 //              excess.set(w, excess[w]+(*flow)[f]);
    596 //              flow->set(f, 0);
    597 //            }
    598 //          }
    599 //        break;
    600 //      } //case PRE_FLOW
    601 //       }
    602 //     } //preflowPreproc
    603 
    604 
    605 
    606 //     void relabel(Node w, int newlevel, VecStack& active,
    607 //               VecNode& level_list, NNMap& left,
    608 //               NNMap& right, int& b, int& k, bool what_heur )
    609 //     {
    610 
    611 //       //FIXME jacint: ez mitol num
    612 // //      Num lev=level[w];
    613 //       int lev=level[w];
    614 
    615 //       Node right_n=right[w];
    616 //       Node left_n=left[w];
    617 
    618 //       //unlacing starts
    619 //       if ( g->valid(right_n) ) {
    620 //      if ( g->valid(left_n) ) {
    621 //        right.set(left_n, right_n);
    622 //        left.set(right_n, left_n);
    623 //      } else {
    624 //        level_list[lev]=right_n;
    625 //        left.set(right_n, INVALID);
    626 //      }
    627 //       } else {
    628 //      if ( g->valid(left_n) ) {
    629 //        right.set(left_n, INVALID);
    630 //      } else {
    631 //        level_list[lev]=INVALID;
    632 //      }
    633 //       }
    634 //       //unlacing ends
    635 
    636 //       if ( !g->valid(level_list[lev]) ) {
    637 
    638 //      //gapping starts
    639 //      for (int i=lev; i!=k ; ) {
    640 //        Node v=level_list[++i];
    641 //        while ( g->valid(v) ) {
    642 //          level.set(v,n);
    643 //          v=right[v];
    644 //        }
    645 //        level_list[i]=INVALID;
    646 //        if ( !what_heur ) {
    647 //          while ( !active[i].empty() ) {
    648 //            active[i].pop();    //FIXME: ezt szebben kene
    649 //          }
    650 //        }
    651 //      }
    652 
    653 //      level.set(w,n);
    654 //      b=lev-1;
    655 //      k=b;
    656 //      //gapping ends
    657 
    658 //       } else {
    659 
    660 //      if ( newlevel == n ) level.set(w,n);
    661 //      else {
    662 //        level.set(w,++newlevel);
    663 //        active[newlevel].push(w);
    664 //        if ( what_heur ) b=newlevel;
    665 //        if ( k < newlevel ) ++k;      //now k=newlevel
    666 //        Node first=level_list[newlevel];
    667 //        if ( g->valid(first) ) left.set(first,w);
    668 //        right.set(w,first);
    669 //        left.set(w,INVALID);
    670 //        level_list[newlevel]=w;
    671 //      }
    672 //       }
    673 
    674 //     } //relabel
    675 
    676 //   };
    677 
    678 
    679 
    680 //   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
    681 //   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1(FlowEnum fe)
    682 //   {
    683 
    684 //     int heur0=(int)(H0*n);  //time while running 'bound decrease'
    685 //     int heur1=(int)(H1*n);  //time while running 'highest label'
    686 //     int heur=heur1;         //starting time interval (#of relabels)
    687 //     int numrelabel=0;
    688 
    689 //     bool what_heur=1;
    690 //     //It is 0 in case 'bound decrease' and 1 in case 'highest label'
    691 
    692 //     bool end=false;
    693 //     //Needed for 'bound decrease', true means no active nodes are above bound
    694 //     //b.
    695 
    696 //     int k=n-2;  //bound on the highest level under n containing a node
    697 //     int b=k;    //bound on the highest level under n of an active node
    698 
    699 //     VecStack active(n);
    700 
    701 //     NNMap left(*g, INVALID);
    702 //     NNMap right(*g, INVALID);
    703 //     VecNode level_list(n,INVALID);
    704 //     //List of the nodes in level i<n, set to n.
    705 
    706 //     NodeIt v;
    707 //     for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
    708 //     //setting each node to level n
    709 
    710 //     if ( fe == NO_FLOW ) {
    711 //       EdgeIt e;
    712 //       for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
    713 //     }
    714 
    715 //     switch (fe) { //computing the excess
    716 //     case PRE_FLOW:
    717 //       {
    718 //      NodeIt v;
    719 //      for(g->first(v); g->valid(v); g->next(v)) {
    720 //        Num exc=0;
    721 
    722 //        InEdgeIt e;
    723 //        for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
    724 //        OutEdgeIt f;
    725 //        for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
    726 
    727 //        excess.set(v,exc);
    728 
    729 //        //putting the active nodes into the stack
    730 //        int lev=level[v];
    731 //        if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
    732 //      }
    733 //      break;
    734 //       }
    735 //     case GEN_FLOW:
    736 //       {
    737 //      NodeIt v;
    738 //      for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
    739 
    740 //      Num exc=0;
    741 //      InEdgeIt e;
    742 //      for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
    743 //      OutEdgeIt f;
    744 //      for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
    745 //      excess.set(t,exc);
    746 //      break;
    747 //       }
    748 //     case ZERO_FLOW:
    749 //     case NO_FLOW:
    750 //       {
    751 //      NodeIt v;
    752 //         for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
    753 //      break;
    754 //       }
    755 //     }
    756 
    757 //     preflowPreproc(fe, active, level_list, left, right);
    758 //     //End of preprocessing
    759 
    760 
    761 //     //Push/relabel on the highest level active nodes.
    762 //     while ( true ) {
    763 //       if ( b == 0 ) {
    764 //      if ( !what_heur && !end && k > 0 ) {
    765 //        b=k;
    766 //        end=true;
    767 //      } else break;
    768 //       }
    769 
    770 //       if ( active[b].empty() ) --b;
    771 //       else {
    772 //      end=false;
    773 //      Node w=active[b].top();
    774 //      active[b].pop();
    775 //      int newlevel=push(w,active);
    776 //      if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
    777 //                                   left, right, b, k, what_heur);
    778 
    779 //      ++numrelabel;
    780 //      if ( numrelabel >= heur ) {
    781 //        numrelabel=0;
    782 //        if ( what_heur ) {
    783 //          what_heur=0;
    784 //          heur=heur0;
    785 //          end=false;
    786 //        } else {
    787 //          what_heur=1;
    788 //          heur=heur1;
    789 //          b=k;
    790 //        }
    791 //      }
    792 //       }
    793 //     }
    794 
    795 //     status=AFTER_PRE_FLOW_PHASE_1;
    796 //   }
    797 
    798 
    799 
    800 //   template <typename Graph, typename Num, typename CapMap, typename FlowMap>
    801 //   void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase2()
    802 //   {
    803 
    804 //     int k=n-2;  //bound on the highest level under n containing a node
    805 //     int b=k;    //bound on the highest level under n of an active node
    806 
    807 //     VecStack active(n);
    808 //     level.set(s,0);
    809 //     std::queue<Node> bfs_queue;
    810 //     bfs_queue.push(s);
    811 
    812 //     while (!bfs_queue.empty()) {
    813 
    814 //       Node v=bfs_queue.front();
    815 //       bfs_queue.pop();
    816 //       int l=level[v]+1;
    817 
    818 //       InEdgeIt e;
    819 //       for(g->first(e,v); g->valid(e); g->next(e)) {
    820 //      if ( (*capacity)[e] <= (*flow)[e] ) continue;
    821 //      Node u=g->tail(e);
    822 //      if ( level[u] >= n ) {
    823 //        bfs_queue.push(u);
    824 //        level.set(u, l);
    825 //        if ( excess[u] > 0 ) active[l].push(u);
    826 //      }
    827 //       }
    828 
    829 //       OutEdgeIt f;
    830 //       for(g->first(f,v); g->valid(f); g->next(f)) {
    831 //      if ( 0 >= (*flow)[f] ) continue;
    832 //      Node u=g->head(f);
    833 //      if ( level[u] >= n ) {
    834 //        bfs_queue.push(u);
    835 //        level.set(u, l);
    836 //        if ( excess[u] > 0 ) active[l].push(u);
    837 //      }
    838 //       }
    839 //     }
    840 //     b=n-2;
    841 
    842 //     while ( true ) {
    843 
    844 //       if ( b == 0 ) break;
    845 
    846 //       if ( active[b].empty() ) --b;
    847 //       else {
    848 //      Node w=active[b].top();
    849 //      active[b].pop();
    850 //      int newlevel=push(w,active);
    851 
    852 //      //relabel
    853 //      if ( excess[w] > 0 ) {
    854 //        level.set(w,++newlevel);
    855 //        active[newlevel].push(w);
    856 //        b=newlevel;
    857 //      }
    858 //       }  // if stack[b] is nonempty
    859 //     } // while(true)
    860 
    861 //     status=AFTER_PRE_FLOW_PHASE_2;
    862 //   }
    863 
    864 
     82  ///\author Marton Makai
    86583  template <typename Graph, typename Num,
    86684            typename CapMap=typename Graph::template EdgeMap<Num>,
     
    87391    typedef typename Graph::OutEdgeIt OutEdgeIt;
    87492    typedef typename Graph::InEdgeIt InEdgeIt;
    875 
    876 //    typedef typename std::vector<std::stack<Node> > VecStack;
    877 //    typedef typename Graph::template NodeMap<Node> NNMap;
    878 //    typedef typename std::vector<Node> VecNode;
    87993
    88094    const Graph* g;
     
    891105    typedef typename Graph::template NodeMap<int> ReachedMap;
    892106
    893 
    894107    //level works as a bool map in augmenting path algorithms and is
    895108    //used by bfs for storing reached information.  In preflow, it
     
    897110    ReachedMap level;
    898111
    899     //excess is needed only in preflow
    900 //    typename Graph::template NodeMap<Num> excess;
    901 
    902     //fixme   
    903 //   protected:
    904     //     MaxFlow() { }
    905     //     void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
    906     //       FlowMap& _flow)
    907     //       {
    908     //  g=&_G;
    909     //  s=_s;
    910     //  t=_t;
    911     //  capacity=&_capacity;
    912     //  flow=&_flow;
    913     //  n=_G.nodeNum;
    914     //  level.set (_G); //kellene vmi ilyesmi fv
    915     //  excess(_G,0); //itt is
    916     //       }
    917 
    918     // constants used for heuristics
    919 //    static const int H0=20;
    920 //    static const int H1=1;
    921 
    922112  public:
    923 
    924113    ///Indicates the property of the starting flow.
    925114
     
    1089278    }
    1090279
    1091     template<typename MapGraphWrapper>
    1092     class DistanceMap {
    1093     protected:
    1094       const MapGraphWrapper* g;
    1095       typename MapGraphWrapper::template NodeMap<int> dist;
    1096     public:
    1097       DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
    1098       void set(const typename MapGraphWrapper::Node& n, int a) {
    1099         dist.set(n, a);
    1100       }
    1101       int operator[](const typename MapGraphWrapper::Node& n) const {
    1102         return dist[n];
    1103       }
    1104       //       int get(const typename MapGraphWrapper::Node& n) const {
    1105       //        return dist[n]; }
    1106       //       bool get(const typename MapGraphWrapper::Edge& e) const {
    1107       //        return (dist.get(g->tail(e))<dist.get(g->head(e))); }
    1108       bool operator[](const typename MapGraphWrapper::Edge& e) const {
    1109         return (dist[g->tail(e)]<dist[g->head(e)]);
    1110       }
    1111     };
    1112 
    1113280  };
    1114281
     
    1245412    {
    1246413      typename ResGW::NodeIt n;
    1247       for(res_graph.first(n); n!=INVALID; ++n) {
     414      for(res_graph.first(n); n!=INVALID; ++n)
    1248415        res_graph_to_F.set(n, F.addNode());
    1249       }
    1250416    }
    1251417
     
    1337503  }
    1338504
    1339 
     505  /// Blocking flow augmentation without constructing the layered
     506  /// graph physically in which the blocking flow is computed.
    1340507  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
    1341508  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
     
    1345512    ResGW res_graph(*g, *capacity, *flow);
    1346513
    1347     //ReachedMap level(res_graph);
     514    //Potential map, for distances from s
     515    typename ResGW::template NodeMap<int> potential(res_graph, 0);
     516    typedef ConstMap<typename ResGW::Edge, int> Const1Map;
     517    Const1Map const_1_map(1);
     518    TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,
     519      Const1Map> tight_edge_filter(res_graph, potential, const_1_map);
     520
    1348521    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
    1349522    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    1350 
    1351523    bfs.pushAndSetReached(s);
    1352     DistanceMap<ResGW> dist(res_graph);
     524
     525    //computing distances from s in the residual graph
    1353526    while ( !bfs.finished() ) {
    1354527      ResGWEdge e=bfs;
    1355       if (e!=INVALID && bfs.isBNodeNewlyReached()) {
    1356         dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
    1357       }
     528      if (e!=INVALID && bfs.isBNodeNewlyReached())
     529        potential.set(res_graph.head(e), potential[res_graph.tail(e)]+1);
    1358530      ++bfs;
    1359     } //computing distances from s in the residual graph
    1360 
    1361     //Subgraph containing the edges on some shortest paths
     531    }
     532
     533    //Subgraph containing the edges on some shortest paths
     534    //(i.e. tight edges)
    1362535    ConstMap<typename ResGW::Node, bool> true_map(true);
    1363536    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
    1364       DistanceMap<ResGW> > FilterResGW;
    1365     FilterResGW filter_res_graph(res_graph, true_map, dist);
     537      TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,
     538      Const1Map> > FilterResGW;
     539    FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);
    1366540
    1367541    //Subgraph, which is able to delete edges which are already
     
    1369543    typename FilterResGW::template NodeMap<typename FilterResGW::Edge>
    1370544      first_out_edges(filter_res_graph);
    1371     typename FilterResGW::NodeIt v;
    1372     for(filter_res_graph.first(v); v!=INVALID; ++v)
    1373       {
    1374         typename FilterResGW::OutEdgeIt e;
    1375         filter_res_graph.first(e, v);
    1376         first_out_edges.set(v, e);
    1377       }
     545    for (typename FilterResGW::NodeIt v(filter_res_graph); v!=INVALID; ++v)
     546      first_out_edges.set
     547        (v, typename FilterResGW::OutEdgeIt(filter_res_graph, v));
     548
    1378549    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
    1379550      template NodeMap<typename FilterResGW::Edge> > ErasingResGW;
     
    1408579      while (!dfs.finished()) {
    1409580        ++dfs;
    1410         if (typename ErasingResGW::Edge(dfs)!=INVALID)
    1411           {
    1412             if (dfs.isBNodeNewlyReached()) {
    1413 
    1414               typename ErasingResGW::Node v=erasing_res_graph.tail(dfs);
    1415               typename ErasingResGW::Node w=erasing_res_graph.head(dfs);
    1416 
    1417               pred.set(w, typename ErasingResGW::Edge(dfs));
    1418               if (pred[v]!=INVALID) {
    1419                 free1.set
    1420                   (w, std::min(free1[v], res_graph.resCap
    1421                                (typename ErasingResGW::Edge(dfs))));
    1422               } else {
    1423                 free1.set
    1424                   (w, res_graph.resCap
    1425                    (typename ErasingResGW::Edge(dfs)));
    1426               }
    1427 
    1428               if (w==t) {
    1429                 __augment=true;
    1430                 _augment=true;
    1431                 break;
    1432               }
    1433             } else {
    1434               erasing_res_graph.erase(dfs);
     581        if (typename ErasingResGW::Edge(dfs)!=INVALID) {
     582          if (dfs.isBNodeNewlyReached()) {
     583           
     584            typename ErasingResGW::Node v=erasing_res_graph.tail(dfs);
     585            typename ErasingResGW::Node w=erasing_res_graph.head(dfs);
     586
     587            pred.set(w, typename ErasingResGW::Edge(dfs));
     588            if (pred[v]!=INVALID) {
     589              free1.set
     590                (w, std::min(free1[v], res_graph.resCap
     591                             (typename ErasingResGW::Edge(dfs))));
     592            } else {
     593              free1.set
     594                (w, res_graph.resCap
     595                 (typename ErasingResGW::Edge(dfs)));
    1435596            }
     597
     598            if (w==t) {
     599              __augment=true;
     600              _augment=true;
     601              break;
     602            }
     603          } else {
     604            erasing_res_graph.erase(dfs);
    1436605          }
     606        }
    1437607      }
    1438608
     
    1440610        typename ErasingResGW::Node
    1441611          n=typename FilterResGW::Node(typename ResGW::Node(t));
    1442         //        typename ResGW::NodeMap<Num> a(res_graph);
    1443         //        typename ResGW::Node b;
    1444         //        Num j=a[b];
    1445         //        typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
    1446         //        typename FilterResGW::Node b1;
    1447         //        Num j1=a1[b1];
    1448         //        typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
    1449         //        typename ErasingResGW::Node b2;
    1450         //        Num j2=a2[b2];
    1451612        Num augment_value=free1[n];
    1452613        while (pred[n]!=INVALID) {
     
    1471632
    1472633
    1473 
    1474 
Note: See TracChangeset for help on using the changeset viewer.