COIN-OR::LEMON - Graph Library

Changeset 2530:f86f7e4eb2ba in lemon-0.x for lemon


Ignore:
Timestamp:
12/04/07 11:55:27 (12 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3407
Message:

Reimplementation of Hao-Orlin algorithm
Little modifictaion in NagamochiIbaraki?
More docs for minimum cut algorithms

Location:
lemon
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/hao_orlin.h

    r2391 r2530  
    2020#define LEMON_HAO_ORLIN_H
    2121
    22 #include <cassert>
    23  
    24 
    25 
    2622#include <vector>
    27 #include <queue>
    2823#include <list>
     24#include <ext/hash_set>
    2925#include <limits>
    3026
     
    3834/// \brief Implementation of the Hao-Orlin algorithm.
    3935///
    40 /// Implementation of the HaoOrlin algorithms class for testing network
     36/// Implementation of the Hao-Orlin algorithm class for testing network
    4137/// reliability.
    4238
     
    4743  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    4844  ///
    49   /// Hao-Orlin calculates a minimum cut in a directed graph
    50   /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
    51   /// of two phases: in the first phase it determines a minimum cut
    52   /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
    53   /// with \f$ source \in X \f$ and minimal out-degree) and in the
    54   /// second phase it determines a minimum cut with \f$ source \f$ on the
    55   /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
    56   /// and minimal out-degree). Obviously, the smaller of these two
    57   /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
    58   /// modified push-relabel preflow algorithm and our implementation
    59   /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
    60   /// highest-label rule). The purpose of such an algorithm is testing
    61   /// network reliability. For an undirected graph with \f$ n \f$
    62   /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
    63   /// and Ibaraki which solves the undirected problem in
    64   /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut
    65   /// algorithm
    66   /// class.
     45  /// Hao-Orlin calculates a minimum cut in a directed graph
     46  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
     47  /// consists of two phases: in the first phase it determines a
     48  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
     49  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
     50  /// out-degree) and in the second phase it determines a minimum cut
     51  /// with \f$ source \f$ on the sink-side (i.e. a set
     52  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
     53  /// out-degree). Obviously, the smaller of these two cuts will be a
     54  /// minimum cut of \f$ D \f$. The algorithm is a modified
     55  /// push-relabel preflow algorithm and our implementation calculates
     56  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
     57  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
     58  /// purpose of such algorithm is testing network reliability. For an
     59  /// undirected graph you can run just the first phase of the
     60  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
     61  /// which solves the undirected problem in
     62  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
     63  /// NagamochiIbaraki algorithm class.
    6764  ///
    6865  /// \param _Graph is the graph type of the algorithm.
     
    8178#endif
    8279  class HaoOrlin {
    83   protected:
     80  private:
    8481
    8582    typedef _Graph Graph;
     
    8986    typedef typename CapacityMap::Value Value;
    9087
     88    GRAPH_TYPEDEFS(typename Graph);
    9189   
    92     typedef typename Graph::Node Node;
    93     typedef typename Graph::NodeIt NodeIt;
    94     typedef typename Graph::EdgeIt EdgeIt;
    95     typedef typename Graph::OutEdgeIt OutEdgeIt;
    96     typedef typename Graph::InEdgeIt InEdgeIt;
    97 
    98     const Graph* _graph;
    99 
     90    const Graph& _graph;
    10091    const CapacityMap* _capacity;
    10192
    10293    typedef typename Graph::template EdgeMap<Value> FlowMap;
    103 
    104     FlowMap* _preflow;
    105 
    106     Node _source, _target;
     94    FlowMap* _flow;
     95
     96    Node _source;
     97
    10798    int _node_num;
    10899
    109     typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
    110                             FlowMap, Tolerance> OutResGraph;
    111     typedef typename OutResGraph::Edge OutResEdge;
     100    // Bucketing structure
     101    std::vector<Node> _first, _last;
     102    typename Graph::template NodeMap<Node>* _next;
     103    typename Graph::template NodeMap<Node>* _prev;   
     104    typename Graph::template NodeMap<bool>* _active;
     105    typename Graph::template NodeMap<int>* _bucket;
    112106   
    113     OutResGraph* _out_res_graph;
    114 
    115     typedef RevGraphAdaptor<const Graph> RevGraph;
    116     RevGraph* _rev_graph;
    117 
    118     typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap,
    119                             FlowMap, Tolerance> InResGraph;
    120     typedef typename InResGraph::Edge InResEdge;
    121    
    122     InResGraph* _in_res_graph;
    123 
    124     typedef IterableBoolMap<Graph, Node> WakeMap;
    125     WakeMap* _wake;
    126 
    127     typedef typename Graph::template NodeMap<int> DistMap;
    128     DistMap* _dist; 
     107    std::vector<bool> _dormant;
     108
     109    std::list<std::list<int> > _sets;
     110    std::list<int>::iterator _highest;
    129111   
    130112    typedef typename Graph::template NodeMap<Value> ExcessMap;
     
    133115    typedef typename Graph::template NodeMap<bool> SourceSetMap;
    134116    SourceSetMap* _source_set;
    135 
    136     std::vector<int> _level_size;
    137 
    138     int _highest_active;
    139     std::vector<std::list<Node> > _active_nodes;
    140 
    141     int _dormant_max;
    142     std::vector<std::list<Node> > _dormant;
    143 
    144117
    145118    Value _min_cut;
     
    157130    HaoOrlin(const Graph& graph, const CapacityMap& capacity,
    158131             const Tolerance& tolerance = Tolerance()) :
    159       _graph(&graph), _capacity(&capacity),
    160       _preflow(0), _source(), _target(),
    161       _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
    162       _wake(0),_dist(0), _excess(0), _source_set(0),
    163       _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
    164       _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
     132      _graph(graph), _capacity(&capacity), _flow(0), _source(),
     133      _node_num(), _first(), _last(), _next(0), _prev(0),
     134      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
     135      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
     136      _tolerance(tolerance) {}
    165137
    166138    ~HaoOrlin() {
     
    168140        delete _min_cut_map;
    169141      }
    170       if (_in_res_graph) {
    171         delete _in_res_graph;
    172       }
    173       if (_rev_graph) {
    174         delete _rev_graph;
    175       }
    176       if (_out_res_graph) {
    177         delete _out_res_graph;
    178       }
    179142      if (_source_set) {
    180143        delete _source_set;
     
    183146        delete _excess;
    184147      }
    185       if (_dist) {
    186         delete _dist;
    187       }
    188       if (_wake) {
    189         delete _wake;
    190       }
    191       if (_preflow) {
    192         delete _preflow;
     148      if (_next) {
     149        delete _next;
     150      }
     151      if (_prev) {
     152        delete _prev;
     153      }
     154      if (_active) {
     155        delete _active;
     156      }
     157      if (_bucket) {
     158        delete _bucket;
     159      }
     160      if (_flow) {
     161        delete _flow;
    193162      }
    194163    }
    195164   
    196165  private:
     166
     167    void activate(const Node& i) {
     168      _active->set(i, true);
     169
     170      int bucket = (*_bucket)[i];
     171
     172      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;       
     173      //unlace
     174      _next->set((*_prev)[i], (*_next)[i]);
     175      if ((*_next)[i] != INVALID) {
     176        _prev->set((*_next)[i], (*_prev)[i]);
     177      } else {
     178        _last[bucket] = (*_prev)[i];
     179      }
     180      //lace
     181      _next->set(i, _first[bucket]);
     182      _prev->set(_first[bucket], i);
     183      _prev->set(i, INVALID);
     184      _first[bucket] = i;
     185    }
     186
     187    void deactivate(const Node& i) {
     188      _active->set(i, false);
     189      int bucket = (*_bucket)[i];
     190
     191      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
     192     
     193      //unlace
     194      _prev->set((*_next)[i], (*_prev)[i]);
     195      if ((*_prev)[i] != INVALID) {
     196        _next->set((*_prev)[i], (*_next)[i]);
     197      } else {
     198        _first[bucket] = (*_next)[i];
     199      }
     200      //lace
     201      _prev->set(i, _last[bucket]);
     202      _next->set(_last[bucket], i);
     203      _next->set(i, INVALID);
     204      _last[bucket] = i;
     205    }
     206
     207    void addItem(const Node& i, int bucket) {
     208      (*_bucket)[i] = bucket;
     209      if (_last[bucket] != INVALID) {
     210        _prev->set(i, _last[bucket]);
     211        _next->set(_last[bucket], i);
     212        _next->set(i, INVALID);
     213        _last[bucket] = i;
     214      } else {
     215        _prev->set(i, INVALID);
     216        _first[bucket] = i;
     217        _next->set(i, INVALID);
     218        _last[bucket] = i;
     219      }
     220    }
    197221   
    198     template <typename ResGraph>
    199     void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
    200       typedef typename ResGraph::Edge ResEdge;
    201       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
    202 
    203       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
    204         (*_preflow)[it] = 0;     
    205       }
    206       for (NodeIt it(*_graph); it != INVALID; ++it) {
    207         (*_wake)[it] = true;
    208         (*_dist)[it] = 1;
    209         (*_excess)[it] = 0;
    210         (*_source_set)[it] = false;
    211       }
    212 
    213       _dormant[0].push_front(_source);
    214       (*_source_set)[_source] = true;
    215       _dormant_max = 0;
    216       (*_wake)[_source] = false;
    217 
    218       _level_size[0] = 1;
    219       _level_size[1] = _node_num - 1;
    220 
    221       _target = target;
    222       (*_dist)[target] = 0;
    223 
    224       for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
    225         Value delta = res_graph.rescap(it);
    226         (*_excess)[_source] -= delta;
    227         res_graph.augment(it, delta);
    228         Node a = res_graph.target(it);
    229         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
    230           _active_nodes[(*_dist)[a]].push_front(a);
    231           if (_highest_active < (*_dist)[a]) {
    232             _highest_active = (*_dist)[a];
    233           }
    234         }
    235         (*_excess)[a] += delta;
    236       }
    237 
    238 
    239       do {
    240         Node n;
    241         while ((n = findActiveNode()) != INVALID) {
    242           for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
    243             Node a = res_graph.target(e);
    244             if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
    245             Value delta = res_graph.rescap(e);
    246             if (_tolerance.positive((*_excess)[n] - delta)) {
    247               (*_excess)[n] -= delta;
     222    void findMinCutOut() {
     223
     224      for (NodeIt n(_graph); n != INVALID; ++n) {
     225        _excess->set(n, 0);
     226      }
     227
     228      for (EdgeIt e(_graph); e != INVALID; ++e) {
     229        _flow->set(e, 0);
     230      }
     231
     232      int bucket_num = 1;
     233     
     234      {
     235        typename Graph::template NodeMap<bool> reached(_graph, false);
     236       
     237        reached.set(_source, true);
     238
     239        bool first_set = true;
     240
     241        for (NodeIt t(_graph); t != INVALID; ++t) {
     242          if (reached[t]) continue;
     243          _sets.push_front(std::list<int>());
     244          _sets.front().push_front(bucket_num);
     245          _dormant[bucket_num] = !first_set;
     246
     247          _bucket->set(t, bucket_num);
     248          _first[bucket_num] = _last[bucket_num] = t;
     249          _next->set(t, INVALID);
     250          _prev->set(t, INVALID);
     251
     252          ++bucket_num;
     253         
     254          std::vector<Node> queue;
     255          queue.push_back(t);
     256          reached.set(t, true);
     257         
     258          while (!queue.empty()) {
     259            _sets.front().push_front(bucket_num);
     260            _dormant[bucket_num] = !first_set;
     261            _first[bucket_num] = _last[bucket_num] = INVALID;
     262           
     263            std::vector<Node> nqueue;
     264            for (int i = 0; i < int(queue.size()); ++i) {
     265              Node n = queue[i];
     266              for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     267                Node u = _graph.source(e);
     268                if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
     269                  reached.set(u, true);
     270                  addItem(u, bucket_num);
     271                  nqueue.push_back(u);
     272                }
     273              }
     274            }
     275            queue.swap(nqueue);
     276            ++bucket_num;
     277          }
     278          _sets.front().pop_front();
     279          --bucket_num;
     280          first_set = false;
     281        }
     282
     283        _bucket->set(_source, 0);
     284        _dormant[0] = true;
     285      }
     286      _source_set->set(_source, true);   
     287         
     288      Node target = _last[_sets.back().back()];
     289      {
     290        for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
     291          if (_tolerance.positive((*_capacity)[e])) {
     292            Node u = _graph.target(e);
     293            _flow->set(e, (*_capacity)[e]);
     294            _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
     295            if (!(*_active)[u] && u != _source) {
     296              activate(u);
     297            }
     298          }
     299        }
     300        if ((*_active)[target]) {
     301          deactivate(target);
     302        }
     303       
     304        _highest = _sets.back().begin();
     305        while (_highest != _sets.back().end() &&
     306               !(*_active)[_first[*_highest]]) {
     307          ++_highest;
     308        }
     309      }
     310
     311
     312      while (true) {
     313        while (_highest != _sets.back().end()) {
     314          Node n = _first[*_highest];
     315          Value excess = (*_excess)[n];
     316          int next_bucket = _node_num;
     317
     318          int under_bucket;
     319          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
     320            under_bucket = -1;
     321          } else {
     322            under_bucket = *(++std::list<int>::iterator(_highest));
     323          }
     324
     325          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     326            Node v = _graph.target(e);
     327            if (_dormant[(*_bucket)[v]]) continue;
     328            Value rem = (*_capacity)[e] - (*_flow)[e];
     329            if (!_tolerance.positive(rem)) continue;
     330            if ((*_bucket)[v] == under_bucket) {
     331              if (!(*_active)[v] && v != target) {
     332                activate(v);
     333              }
     334              if (!_tolerance.less(rem, excess)) {
     335                _flow->set(e, (*_flow)[e] + excess);
     336                _excess->set(v, (*_excess)[v] + excess);
     337                excess = 0;
     338                goto no_more_push;
     339              } else {
     340                excess -= rem;
     341                _excess->set(v, (*_excess)[v] + rem);
     342                _flow->set(e, (*_capacity)[e]);
     343              }
     344            } else if (next_bucket > (*_bucket)[v]) {
     345              next_bucket = (*_bucket)[v];
     346            }
     347          }
     348
     349          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     350            Node v = _graph.source(e);
     351            if (_dormant[(*_bucket)[v]]) continue;
     352            Value rem = (*_flow)[e];
     353            if (!_tolerance.positive(rem)) continue;
     354            if ((*_bucket)[v] == under_bucket) {
     355              if (!(*_active)[v] && v != target) {
     356                activate(v);
     357              }
     358              if (!_tolerance.less(rem, excess)) {
     359                _flow->set(e, (*_flow)[e] - excess);
     360                _excess->set(v, (*_excess)[v] + excess);
     361                excess = 0;
     362                goto no_more_push;
     363              } else {
     364                excess -= rem;
     365                _excess->set(v, (*_excess)[v] + rem);
     366                _flow->set(e, 0);
     367              }
     368            } else if (next_bucket > (*_bucket)[v]) {
     369              next_bucket = (*_bucket)[v];
     370            }
     371          }
     372         
     373        no_more_push:
     374         
     375          _excess->set(n, excess);
     376         
     377          if (excess != 0) {
     378            if ((*_next)[n] == INVALID) {
     379              typename std::list<std::list<int> >::iterator new_set =
     380                _sets.insert(--_sets.end(), std::list<int>());
     381              new_set->splice(new_set->end(), _sets.back(),
     382                              _sets.back().begin(), ++_highest);
     383              for (std::list<int>::iterator it = new_set->begin();
     384                   it != new_set->end(); ++it) {
     385                _dormant[*it] = true;
     386              }
     387              while (_highest != _sets.back().end() &&
     388                     !(*_active)[_first[*_highest]]) {
     389                ++_highest;
     390              }
     391            } else if (next_bucket == _node_num) {
     392              _first[(*_bucket)[n]] = (*_next)[n];
     393              _prev->set((*_next)[n], INVALID);
     394             
     395              std::list<std::list<int> >::iterator new_set =
     396                _sets.insert(--_sets.end(), std::list<int>());
     397             
     398              new_set->push_front(bucket_num);
     399              _bucket->set(n, bucket_num);
     400              _first[bucket_num] = _last[bucket_num] = n;
     401              _next->set(n, INVALID);
     402              _prev->set(n, INVALID);
     403              _dormant[bucket_num] = true;
     404              ++bucket_num;
     405
     406              while (_highest != _sets.back().end() &&
     407                     !(*_active)[_first[*_highest]]) {
     408                ++_highest;
     409              }
    248410            } else {
    249               delta = (*_excess)[n];
    250               (*_excess)[n] = 0;
    251             }
    252             res_graph.augment(e, delta);
    253             if ((*_excess)[a] == 0 && a != _target) {
    254               _active_nodes[(*_dist)[a]].push_front(a);
    255             }
    256             (*_excess)[a] += delta;
    257             if ((*_excess)[n] == 0) break;
    258           }
    259           if ((*_excess)[n] != 0) {
    260             relabel(n, res_graph);
    261           }
    262         }
    263 
    264         Value current_value = cutValue(out);
    265         if (_min_cut > current_value){
    266           if (out) {
    267             for (NodeIt it(*_graph); it != INVALID; ++it) {
    268               _min_cut_map->set(it, !(*_wake)[it]);
    269             }
    270           } else {
    271             for (NodeIt it(*_graph); it != INVALID; ++it) {
    272               _min_cut_map->set(it, (*_wake)[it]);
    273             }
    274           }
    275 
    276           _min_cut = current_value;
    277         }
    278 
    279       } while (selectNewSink(res_graph));
    280     }
    281 
    282     template <typename ResGraph>
    283     void relabel(const Node& n, ResGraph& res_graph) {
    284       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
    285 
    286       int k = (*_dist)[n];
    287       if (_level_size[k] == 1) {
    288         ++_dormant_max;
    289         for (NodeIt it(*_graph); it != INVALID; ++it) {
    290           if ((*_wake)[it] && (*_dist)[it] >= k) {
    291             (*_wake)[it] = false;
    292             _dormant[_dormant_max].push_front(it);
    293             --_level_size[(*_dist)[it]];
    294           }
    295         }
    296         --_highest_active;
    297       } else { 
    298         int new_dist = _node_num;
    299         for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
    300           Node t = res_graph.target(e);
    301           if ((*_wake)[t] && new_dist > (*_dist)[t]) {
    302             new_dist = (*_dist)[t];
    303           }
    304         }
    305         if (new_dist == _node_num) {
    306           ++_dormant_max;
    307           (*_wake)[n] = false;
    308           _dormant[_dormant_max].push_front(n);
    309           --_level_size[(*_dist)[n]];
    310         } else {           
    311           --_level_size[(*_dist)[n]];
    312           (*_dist)[n] = new_dist + 1;
    313           _highest_active = (*_dist)[n];
    314           _active_nodes[_highest_active].push_front(n);
    315           ++_level_size[(*_dist)[n]];
    316         }
    317       }
    318     }
    319 
    320     template <typename ResGraph>
    321     bool selectNewSink(ResGraph& res_graph) {
    322       typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
    323 
    324       Node old_target = _target;
    325       (*_wake)[_target] = false;
    326       --_level_size[(*_dist)[_target]];
    327       _dormant[0].push_front(_target);
    328       (*_source_set)[_target] = true;
    329       if (int(_dormant[0].size()) == _node_num){
    330         _dormant[0].clear();
    331         return false;
    332       }
    333 
    334       bool wake_was_empty = false;
    335 
    336       if(_wake->trueNum() == 0) {
    337         while (!_dormant[_dormant_max].empty()){
    338           (*_wake)[_dormant[_dormant_max].front()] = true;
    339           ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
    340           _dormant[_dormant_max].pop_front();
    341         }
    342         --_dormant_max;
    343         wake_was_empty = true;
    344       }
    345 
    346       int min_dist = std::numeric_limits<int>::max();
    347       for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
    348         if (min_dist > (*_dist)[it]){
    349           _target = it;
    350           min_dist = (*_dist)[it];
    351         }
    352       }
    353 
    354       if (wake_was_empty){
    355         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
    356           if ((*_excess)[it] != 0 && it != _target) {
    357             _active_nodes[(*_dist)[it]].push_front(it);
    358             if (_highest_active < (*_dist)[it]) {
    359               _highest_active = (*_dist)[it];               
    360             }
    361           }
    362         }
    363       }
    364 
    365       Node n = old_target;
    366       for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
    367         Node a = res_graph.target(e);
    368         if (!(*_wake)[a]) continue;
    369         Value delta = res_graph.rescap(e);
    370         res_graph.augment(e, delta);
    371         (*_excess)[n] -= delta;
    372         if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
    373           _active_nodes[(*_dist)[a]].push_front(a);
    374           if (_highest_active < (*_dist)[a]) {
    375             _highest_active = (*_dist)[a];
    376           }
    377         }
    378         (*_excess)[a] += delta;
    379       }
     411              _first[*_highest] = (*_next)[n];
     412              _prev->set((*_next)[n], INVALID);
     413             
     414              while (next_bucket != *_highest) {
     415                --_highest;
     416              }
     417
     418              if (_highest == _sets.back().begin()) {
     419                _sets.back().push_front(bucket_num);
     420                _dormant[bucket_num] = false;
     421                _first[bucket_num] = _last[bucket_num] = INVALID;
     422                ++bucket_num;
     423              }
     424              --_highest;
     425
     426              _bucket->set(n, *_highest);
     427              _next->set(n, _first[*_highest]);
     428              if (_first[*_highest] != INVALID) {
     429                _prev->set(_first[*_highest], n);
     430              } else {
     431                _last[*_highest] = n;
     432              }
     433              _first[*_highest] = n;         
     434            }
     435          } else {
     436
     437            deactivate(n);
     438            if (!(*_active)[_first[*_highest]]) {
     439              ++_highest;
     440              if (_highest != _sets.back().end() &&
     441                  !(*_active)[_first[*_highest]]) {
     442                _highest = _sets.back().end();
     443              }
     444            }
     445          }
     446        }
     447
     448        if ((*_excess)[target] < _min_cut) {
     449          _min_cut = (*_excess)[target];
     450          for (NodeIt i(_graph); i != INVALID; ++i) {
     451            _min_cut_map->set(i, true);
     452          }
     453          for (std::list<int>::iterator it = _sets.back().begin();
     454               it != _sets.back().end(); ++it) {
     455            Node n = _first[*it];
     456            while (n != INVALID) {
     457              _min_cut_map->set(n, false);
     458              n = (*_next)[n];
     459            }
     460          }
     461        }
     462
     463        {
     464          Node new_target;
     465          if ((*_prev)[target] != INVALID) {
     466            _last[(*_bucket)[target]] = (*_prev)[target];
     467            _next->set((*_prev)[target], INVALID);
     468            new_target = (*_prev)[target];
     469          } else {
     470            _sets.back().pop_back();
     471            if (_sets.back().empty()) {
     472              _sets.pop_back();
     473              if (_sets.empty())
     474                break;
     475              for (std::list<int>::iterator it = _sets.back().begin();
     476                   it != _sets.back().end(); ++it) {
     477                _dormant[*it] = false;
     478              }
     479            }
     480            new_target = _last[_sets.back().back()];
     481          }
     482
     483          _bucket->set(target, 0);
     484
     485          _source_set->set(target, true);         
     486          for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
     487            Value rem = (*_capacity)[e] - (*_flow)[e];
     488            if (!_tolerance.positive(rem)) continue;
     489            Node v = _graph.target(e);
     490            if (!(*_active)[v] && !(*_source_set)[v]) {
     491              activate(v);
     492            }
     493            _excess->set(v, (*_excess)[v] + rem);
     494            _flow->set(e, (*_capacity)[e]);
     495          }
     496         
     497          for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
     498            Value rem = (*_flow)[e];
     499            if (!_tolerance.positive(rem)) continue;
     500            Node v = _graph.source(e);
     501            if (!(*_active)[v] && !(*_source_set)[v]) {
     502              activate(v);
     503            }
     504            _excess->set(v, (*_excess)[v] + rem);
     505            _flow->set(e, 0);
     506          }
     507         
     508          target = new_target;
     509          if ((*_active)[target]) {
     510            deactivate(target);
     511          }
     512
     513          _highest = _sets.back().begin();
     514          while (_highest != _sets.back().end() &&
     515                 !(*_active)[_first[*_highest]]) {
     516            ++_highest;
     517          }
     518        }
     519      }
     520    }   
     521
     522    void findMinCutIn() {
     523
     524      for (NodeIt n(_graph); n != INVALID; ++n) {
     525        _excess->set(n, 0);
     526      }
     527
     528      for (EdgeIt e(_graph); e != INVALID; ++e) {
     529        _flow->set(e, 0);
     530      }
     531
     532      int bucket_num = 1;
    380533     
    381       return true;
    382     }
    383 
    384     Node findActiveNode() {
    385       while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
    386         --_highest_active;
    387       }
    388       if( _highest_active > 0) {
    389         Node n = _active_nodes[_highest_active].front();
    390         _active_nodes[_highest_active].pop_front();
    391         return n;
    392       } else {
    393         return INVALID;
    394       }
    395     }
    396 
    397     Value cutValue(bool out) {
    398       Value value = 0;
    399       if (out) {
    400         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
    401           for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
    402             if (!(*_wake)[_graph->source(e)]){
    403               value += (*_capacity)[e];
    404             }   
    405           }
    406         }
    407       } else {
    408         for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
    409           for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
    410             if (!(*_wake)[_graph->target(e)]){
    411               value += (*_capacity)[e];
    412             }   
    413           }
    414         }
    415       }
    416       return value;
    417     }
    418 
     534      {
     535        typename Graph::template NodeMap<bool> reached(_graph, false);
     536       
     537        reached.set(_source, true);
     538
     539        bool first_set = true;
     540
     541        for (NodeIt t(_graph); t != INVALID; ++t) {
     542          if (reached[t]) continue;
     543          _sets.push_front(std::list<int>());
     544          _sets.front().push_front(bucket_num);
     545          _dormant[bucket_num] = !first_set;
     546
     547          _bucket->set(t, bucket_num);
     548          _first[bucket_num] = _last[bucket_num] = t;
     549          _next->set(t, INVALID);
     550          _prev->set(t, INVALID);
     551
     552          ++bucket_num;
     553         
     554          std::vector<Node> queue;
     555          queue.push_back(t);
     556          reached.set(t, true);
     557         
     558          while (!queue.empty()) {
     559            _sets.front().push_front(bucket_num);
     560            _dormant[bucket_num] = !first_set;
     561            _first[bucket_num] = _last[bucket_num] = INVALID;
     562           
     563            std::vector<Node> nqueue;
     564            for (int i = 0; i < int(queue.size()); ++i) {
     565              Node n = queue[i];
     566              for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     567                Node u = _graph.target(e);
     568                if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
     569                  reached.set(u, true);
     570                  addItem(u, bucket_num);
     571                  nqueue.push_back(u);
     572                }
     573              }
     574            }
     575            queue.swap(nqueue);
     576            ++bucket_num;
     577          }
     578          _sets.front().pop_front();
     579          --bucket_num;
     580          first_set = false;
     581        }
     582
     583        _bucket->set(_source, 0);
     584        _dormant[0] = true;
     585      }
     586      _source_set->set(_source, true);   
     587         
     588      Node target = _last[_sets.back().back()];
     589      {
     590        for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
     591          if (_tolerance.positive((*_capacity)[e])) {
     592            Node u = _graph.source(e);
     593            _flow->set(e, (*_capacity)[e]);
     594            _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
     595            if (!(*_active)[u] && u != _source) {
     596              activate(u);
     597            }
     598          }
     599        }
     600        if ((*_active)[target]) {
     601          deactivate(target);
     602        }
     603       
     604        _highest = _sets.back().begin();
     605        while (_highest != _sets.back().end() &&
     606               !(*_active)[_first[*_highest]]) {
     607          ++_highest;
     608        }
     609      }
     610
     611
     612      while (true) {
     613        while (_highest != _sets.back().end()) {
     614          Node n = _first[*_highest];
     615          Value excess = (*_excess)[n];
     616          int next_bucket = _node_num;
     617
     618          int under_bucket;
     619          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
     620            under_bucket = -1;
     621          } else {
     622            under_bucket = *(++std::list<int>::iterator(_highest));
     623          }
     624
     625          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     626            Node v = _graph.source(e);
     627            if (_dormant[(*_bucket)[v]]) continue;
     628            Value rem = (*_capacity)[e] - (*_flow)[e];
     629            if (!_tolerance.positive(rem)) continue;
     630            if ((*_bucket)[v] == under_bucket) {
     631              if (!(*_active)[v] && v != target) {
     632                activate(v);
     633              }
     634              if (!_tolerance.less(rem, excess)) {
     635                _flow->set(e, (*_flow)[e] + excess);
     636                _excess->set(v, (*_excess)[v] + excess);
     637                excess = 0;
     638                goto no_more_push;
     639              } else {
     640                excess -= rem;
     641                _excess->set(v, (*_excess)[v] + rem);
     642                _flow->set(e, (*_capacity)[e]);
     643              }
     644            } else if (next_bucket > (*_bucket)[v]) {
     645              next_bucket = (*_bucket)[v];
     646            }
     647          }
     648
     649          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     650            Node v = _graph.target(e);
     651            if (_dormant[(*_bucket)[v]]) continue;
     652            Value rem = (*_flow)[e];
     653            if (!_tolerance.positive(rem)) continue;
     654            if ((*_bucket)[v] == under_bucket) {
     655              if (!(*_active)[v] && v != target) {
     656                activate(v);
     657              }
     658              if (!_tolerance.less(rem, excess)) {
     659                _flow->set(e, (*_flow)[e] - excess);
     660                _excess->set(v, (*_excess)[v] + excess);
     661                excess = 0;
     662                goto no_more_push;
     663              } else {
     664                excess -= rem;
     665                _excess->set(v, (*_excess)[v] + rem);
     666                _flow->set(e, 0);
     667              }
     668            } else if (next_bucket > (*_bucket)[v]) {
     669              next_bucket = (*_bucket)[v];
     670            }
     671          }
     672         
     673        no_more_push:
     674         
     675          _excess->set(n, excess);
     676         
     677          if (excess != 0) {
     678            if ((*_next)[n] == INVALID) {
     679              typename std::list<std::list<int> >::iterator new_set =
     680                _sets.insert(--_sets.end(), std::list<int>());
     681              new_set->splice(new_set->end(), _sets.back(),
     682                              _sets.back().begin(), ++_highest);
     683              for (std::list<int>::iterator it = new_set->begin();
     684                   it != new_set->end(); ++it) {
     685                _dormant[*it] = true;
     686              }
     687              while (_highest != _sets.back().end() &&
     688                     !(*_active)[_first[*_highest]]) {
     689                ++_highest;
     690              }
     691            } else if (next_bucket == _node_num) {
     692              _first[(*_bucket)[n]] = (*_next)[n];
     693              _prev->set((*_next)[n], INVALID);
     694             
     695              std::list<std::list<int> >::iterator new_set =
     696                _sets.insert(--_sets.end(), std::list<int>());
     697             
     698              new_set->push_front(bucket_num);
     699              _bucket->set(n, bucket_num);
     700              _first[bucket_num] = _last[bucket_num] = n;
     701              _next->set(n, INVALID);
     702              _prev->set(n, INVALID);
     703              _dormant[bucket_num] = true;
     704              ++bucket_num;
     705
     706              while (_highest != _sets.back().end() &&
     707                     !(*_active)[_first[*_highest]]) {
     708                ++_highest;
     709              }
     710            } else {
     711              _first[*_highest] = (*_next)[n];
     712              _prev->set((*_next)[n], INVALID);
     713
     714              while (next_bucket != *_highest) {
     715                --_highest;
     716              }
     717              if (_highest == _sets.back().begin()) {
     718                _sets.back().push_front(bucket_num);
     719                _dormant[bucket_num] = false;
     720                _first[bucket_num] = _last[bucket_num] = INVALID;
     721                ++bucket_num;
     722              }
     723              --_highest;
     724
     725              _bucket->set(n, *_highest);
     726              _next->set(n, _first[*_highest]);
     727              if (_first[*_highest] != INVALID) {
     728                _prev->set(_first[*_highest], n);
     729              } else {
     730                _last[*_highest] = n;
     731              }
     732              _first[*_highest] = n;         
     733            }
     734          } else {
     735
     736            deactivate(n);
     737            if (!(*_active)[_first[*_highest]]) {
     738              ++_highest;
     739              if (_highest != _sets.back().end() &&
     740                  !(*_active)[_first[*_highest]]) {
     741                _highest = _sets.back().end();
     742              }
     743            }
     744          }
     745        }
     746
     747        if ((*_excess)[target] < _min_cut) {
     748          _min_cut = (*_excess)[target];
     749          for (NodeIt i(_graph); i != INVALID; ++i) {
     750            _min_cut_map->set(i, false);
     751          }
     752          for (std::list<int>::iterator it = _sets.back().begin();
     753               it != _sets.back().end(); ++it) {
     754            Node n = _first[*it];
     755            while (n != INVALID) {
     756              _min_cut_map->set(n, true);
     757              n = (*_next)[n];
     758            }
     759          }
     760        }
     761
     762        {
     763          Node new_target;
     764          if ((*_prev)[target] != INVALID) {
     765            _last[(*_bucket)[target]] = (*_prev)[target];
     766            _next->set((*_prev)[target], INVALID);
     767            new_target = (*_prev)[target];
     768          } else {
     769            _sets.back().pop_back();
     770            if (_sets.back().empty()) {
     771              _sets.pop_back();
     772              if (_sets.empty())
     773                break;
     774              for (std::list<int>::iterator it = _sets.back().begin();
     775                   it != _sets.back().end(); ++it) {
     776                _dormant[*it] = false;
     777              }
     778            }
     779            new_target = _last[_sets.back().back()];
     780          }
     781
     782          _bucket->set(target, 0);
     783
     784          _source_set->set(target, true);         
     785          for (InEdgeIt e(_graph, target); e != INVALID; ++e) {
     786            Value rem = (*_capacity)[e] - (*_flow)[e];
     787            if (!_tolerance.positive(rem)) continue;
     788            Node v = _graph.source(e);
     789            if (!(*_active)[v] && !(*_source_set)[v]) {
     790              activate(v);
     791            }
     792            _excess->set(v, (*_excess)[v] + rem);
     793            _flow->set(e, (*_capacity)[e]);
     794          }
     795         
     796          for (OutEdgeIt e(_graph, target); e != INVALID; ++e) {
     797            Value rem = (*_flow)[e];
     798            if (!_tolerance.positive(rem)) continue;
     799            Node v = _graph.target(e);
     800            if (!(*_active)[v] && !(*_source_set)[v]) {
     801              activate(v);
     802            }
     803            _excess->set(v, (*_excess)[v] + rem);
     804            _flow->set(e, 0);
     805          }
     806         
     807          target = new_target;
     808          if ((*_active)[target]) {
     809            deactivate(target);
     810          }
     811
     812          _highest = _sets.back().begin();
     813          while (_highest != _sets.back().end() &&
     814                 !(*_active)[_first[*_highest]]) {
     815            ++_highest;
     816          }
     817        }
     818      }
     819    }
    419820
    420821  public:
     
    436837    /// for the algorithm.
    437838    void init() {
    438       init(NodeIt(*_graph));
     839      init(NodeIt(_graph));
    439840    }
    440841
     
    447848    void init(const Node& source) {
    448849      _source = source;
    449       _node_num = countNodes(*_graph);
     850     
     851      _node_num = countNodes(_graph);
     852     
     853      _first.resize(_node_num);
     854      _last.resize(_node_num);
    450855
    451856      _dormant.resize(_node_num);
    452       _level_size.resize(_node_num, 0);
    453       _active_nodes.resize(_node_num);
    454 
    455       if (!_preflow) {
    456         _preflow = new FlowMap(*_graph);
    457       }
    458       if (!_wake) {
    459         _wake = new WakeMap(*_graph);
    460       }
    461       if (!_dist) {
    462         _dist = new DistMap(*_graph);
     857
     858      if (!_flow) {
     859        _flow = new FlowMap(_graph);
     860      }
     861      if (!_next) {
     862        _next = new typename Graph::template NodeMap<Node>(_graph);
     863      }
     864      if (!_prev) {
     865        _prev = new typename Graph::template NodeMap<Node>(_graph);
     866      }
     867      if (!_active) {
     868        _active = new typename Graph::template NodeMap<bool>(_graph);
     869      }
     870      if (!_bucket) {
     871        _bucket = new typename Graph::template NodeMap<int>(_graph);
    463872      }
    464873      if (!_excess) {
    465         _excess = new ExcessMap(*_graph);
     874        _excess = new ExcessMap(_graph);
    466875      }
    467876      if (!_source_set) {
    468         _source_set = new SourceSetMap(*_graph);
    469       }
    470       if (!_out_res_graph) {
    471         _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
    472       }
    473       if (!_rev_graph) {
    474         _rev_graph = new RevGraph(*_graph);
    475       }
    476       if (!_in_res_graph) {
    477         _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
     877        _source_set = new SourceSetMap(_graph);
    478878      }
    479879      if (!_min_cut_map) {
    480         _min_cut_map = new MinCutMap(*_graph);
     880        _min_cut_map = new MinCutMap(_graph);
    481881      }
    482882
     
    488888    /// source-side.
    489889    ///
     890    /// Calculates a minimum cut with \f$ source \f$ on the
     891    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
     892    /// \in X \f$ and minimal out-degree).
     893    void calculateOut() {
     894      findMinCutOut();
     895    }
     896
    490897    /// \brief Calculates a minimum cut with \f$ source \f$ on the
    491     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
    492     ///  and minimal out-degree).
    493     void calculateOut() {
    494       for (NodeIt it(*_graph); it != INVALID; ++it) {
    495         if (it != _source) {
    496           calculateOut(it);
    497           return;
    498         }
    499       }
    500     }
    501 
    502     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    503     /// source-side.
     898    /// target-side.
    504899    ///
    505     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    506     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
    507     ///  and minimal out-degree). The \c target is the initial target
    508     /// for the push-relabel algorithm.
    509     void calculateOut(const Node& target) {
    510       findMinCut(target, true, *_out_res_graph);
    511     }
    512 
    513     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    514     /// sink-side.
    515     ///
    516     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    517     /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
    518     /// \f$ source \notin X \f$
    519     /// and minimal out-degree).
     900    /// Calculates a minimum cut with \f$ source \f$ on the
     901    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
     902    /// \in X \f$ and minimal out-degree).
    520903    void calculateIn() {
    521       for (NodeIt it(*_graph); it != INVALID; ++it) {
    522         if (it != _source) {
    523           calculateIn(it);
    524           return;
    525         }
    526       }
    527     }
    528 
    529     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    530     /// sink-side.
    531     ///
    532     /// \brief Calculates a minimum cut with \f$ source \f$ on the
    533     /// sink-side (i.e. a set \f$ X\subsetneq V
    534     /// \f$ with \f$ source \notin X \f$ and minimal out-degree). 
    535     /// The \c target is the initial
    536     /// target for the push-relabel algorithm.
    537     void calculateIn(const Node& target) {
    538       findMinCut(target, false, *_in_res_graph);
    539     }
     904      findMinCutIn();
     905    }
     906
    540907
    541908    /// \brief Runs the algorithm.
     
    546913    void run() {
    547914      init();
    548       for (NodeIt it(*_graph); it != INVALID; ++it) {
    549         if (it != _source) {
    550           calculateOut(it);
    551           calculateIn(it);
    552           return;
    553         }
    554       }
     915      calculateOut();
     916      calculateIn();
    555917    }
    556918
     
    562924    void run(const Node& s) {
    563925      init(s);
    564       for (NodeIt it(*_graph); it != INVALID; ++it) {
    565         if (it != _source) {
    566           calculateOut(it);
    567           calculateIn(it);
    568           return;
    569         }
    570       }
    571     }
    572 
    573     /// \brief Runs the algorithm.
    574     ///
    575     /// Runs the algorithm. It just calls the \ref init() and then
    576     /// \ref calculateOut() and \ref calculateIn().
    577     void run(const Node& s, const Node& t) {
    578       init(s);
    579       calculateOut(t);
    580       calculateIn(t);
     926      calculateOut();
     927      calculateIn();
    581928    }
    582929
     
    609956    template <typename NodeMap>
    610957    Value minCut(NodeMap& nodeMap) const {
    611       for (NodeIt it(*_graph); it != INVALID; ++it) {
     958      for (NodeIt it(_graph); it != INVALID; ++it) {
    612959        nodeMap.set(it, (*_min_cut_map)[it]);
    613960      }
  • lemon/nagamochi_ibaraki.h

    r2506 r2530  
    848848  ///
    849849  /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
    850   /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$. 
    851   /// When capacity map is neutral then it uses BucketHeap which
    852   /// results \f$ O(ne) \f$ time complexity.
     850  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$.
     851  /// When unit capacity minimum cut is computed then it uses
     852  /// BucketHeap which results \f$ O(ne) \f$ time complexity.
    853853  ///
    854854  /// \warning The value type of the capacity map should be able to hold
     
    907907    ///@{
    908908
    909     struct DefNeutralCapacityTraits : public Traits {
     909    struct DefUnitCapacityTraits : public Traits {
    910910      typedef ConstMap<typename Graph::UEdge, Const<int, 1> > CapacityMap;
    911911      static CapacityMap *createCapacityMap(const Graph&) {
     
    918918    /// \ref named-templ-param "Named parameter" for setting
    919919    /// the capacity type to constMap<UEdge, int, 1>()
    920     struct DefNeutralCapacity
     920    struct DefUnitCapacity
    921921      : public NagamochiIbaraki<Graph, CapacityMap,
    922                                 DefNeutralCapacityTraits> {
     922                                DefUnitCapacityTraits> {
    923923      typedef NagamochiIbaraki<Graph, CapacityMap,
    924                                DefNeutralCapacityTraits> Create;
     924                               DefUnitCapacityTraits> Create;
    925925    };
    926926
     
    11351135    /// This constructor can be used only when the Traits class
    11361136    /// defines how can we instantiate a local capacity map.
    1137     /// If the DefNeutralCapacity used the algorithm automatically
     1137    /// If the DefUnitCapacity used the algorithm automatically
    11381138    /// construct the capacity map.
    11391139    ///
Note: See TracChangeset for help on using the changeset viewer.