COIN-OR::LEMON - Graph Library

Changeset 2337:9c3d44ac39fb in lemon-0.x


Ignore:
Timestamp:
01/11/07 22:05:00 (13 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3129
Message:

Adding two heuristics
Based on:
http://www.avglab.com/andrew/pub/neci-tr-96-132.ps

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/nagamochi_ibaraki.h

    r2284 r2337  
    1515 */
    1616
    17 #ifndef LEMON_MIN_CUT_H
    18 #define LEMON_MIN_CUT_H
    19 
    20 
    21 /// \ingroup topology
    22 /// \file
    23 /// \brief Maximum cardinality search and min cut in undirected graphs.
     17#ifndef LEMON_NAGAMOCHI_IBARAKI_H
     18#define LEMON_NAGAMOCHI_IBARAKI_H
     19
     20
     21/// \ingroup topology
     22/// \file
     23/// \brief Maximum cardinality search and minimum cut in undirected
     24/// graphs.
    2425
    2526#include <lemon/list_graph.h>
     
    2728#include <lemon/bucket_heap.h>
    2829
     30#include <lemon/unionfind.h>
     31#include <lemon/topology.h>
     32
    2933#include <lemon/bits/invalid.h>
    3034#include <lemon/error.h>
     
    3236
    3337#include <functional>
     38
     39#include <lemon/graph_writer.h>
     40#include <lemon/time_measure.h>
    3441
    3542namespace lemon {
     
    146153      return new CardinalityMap(graph);
    147154    }
     155
    148156
    149157  };
     
    666674    Value cardinality(Node node) const { return (*_cardinality)[node]; }
    667675
     676    /// \brief The current cardinality of a node.
     677    ///
     678    /// Returns the current cardinality of a node.
     679    /// \pre the given node should be reached but not processed
     680    Value currentCardinality(Node node) const { return (*_heap)[node]; }
     681
    668682    /// \brief Returns a reference to the NodeMap of cardinalities.
    669683    ///
     
    730744    }
    731745
     746    /// \brief The CutValueMap type
     747    ///
     748    /// The type of the map that stores the cut value of a node.
     749    typedef AuxGraph::NodeMap<Value> AuxCutValueMap;
     750
     751    /// \brief Instantiates a AuxCutValueMap.
     752    ///
     753    /// This function instantiates a \ref AuxCutValueMap.
     754    static AuxCutValueMap *createAuxCutValueMap(const AuxGraph& graph) {
     755      return new AuxCutValueMap(graph);
     756    }
     757
    732758    /// \brief The AuxCapacityMap type
    733759    ///
    734     /// The type of the map that stores the auxing edge capacities.
     760    /// The type of the map that stores the auxiliary edge capacities.
    735761    typedef AuxGraph::UEdgeMap<Value> AuxCapacityMap;
    736762
     
    806832
    807833  };
    808 
    809   namespace _min_cut_bits {
    810     template <typename _Key>
    811     class LastTwoMap {
    812     public:
    813       typedef _Key Key;
    814       typedef bool Value;
    815 
    816       LastTwoMap(int _num) : num(_num) {}
    817       void set(const Key& key, bool val) {
    818         if (!val) return;
    819         --num;
    820         if (num > 1) return;
    821         keys[num] = key;
    822       }
    823      
    824       Key operator[](int index) const { return keys[index]; }
    825     private:
    826       Key keys[2];
    827       int num;
    828     };
    829   }
    830834
    831835  /// \ingroup topology
     
    842846  ///
    843847  /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
    844   /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n))
    845   /// \f$. When capacity map is neutral then it uses BucketHeap which
     848  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$.
     849  /// When capacity map is neutral then it uses BucketHeap which
    846850  /// results \f$ O(ne) \f$ time complexity.
     851  ///
     852  /// \warning The value type of the capacity map should be able to hold
     853  /// any cut value of the graph, otherwise the result can overflow.
    847854#ifdef DOXYGEN
    848855  template <typename _Graph, typename _CapacityMap, typename _Traits>
     
    881888    /// The type of the aux capacity map
    882889    typedef typename Traits::AuxCapacityMap AuxCapacityMap;
     890    /// The type of the aux cut value map
     891    typedef typename Traits::AuxCutValueMap AuxCutValueMap;
    883892    /// The cross reference type used for the current heap.
    884893    typedef typename Traits::HeapCrossRef HeapCrossRef;
     
    989998    /// (\c true) or not.
    990999    bool local_aux_capacity;
     1000    /// Pointer to the aux cut value map
     1001    AuxCutValueMap *_aux_cut_value;
     1002    /// \brief Indicates if \ref _aux_cut_value is locally allocated
     1003    /// (\c true) or not.
     1004    bool local_aux_cut_value;
    9911005    /// Pointer to the heap cross references.
    9921006    HeapCrossRef *_heap_cross_ref;
     
    10031017    /// The number of the nodes of the aux graph.
    10041018    int _node_num;
    1005     /// The first and last node of the min cut in the next list;
    1006     typename Graph::Node _first_node, _last_node;
     1019    /// The first and last node of the min cut in the next list.
     1020    std::vector<typename Graph::Node> _cut;
    10071021
    10081022    /// \brief The first and last element in the list associated
     
    10121026    ListRefMap *_next;
    10131027   
    1014     void create_structures() {
     1028    void createStructures() {
    10151029      if (!_capacity) {
    10161030        local_capacity = true;
     
    10251039        _aux_capacity = Traits::createAuxCapacityMap(*_aux_graph);
    10261040      }
     1041      if(!_aux_cut_value) {
     1042        local_aux_cut_value = true;
     1043        _aux_cut_value = Traits::createAuxCutValueMap(*_aux_graph);
     1044      }
    10271045
    10281046      _first = Traits::createNodeRefMap(*_aux_graph);
     
    10301048
    10311049      _next = Traits::createListRefMap(*_graph);
    1032 
    1033       typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
    1034 
    1035       for (typename Graph::NodeIt it(*_graph); it != INVALID; ++it) {
    1036         _next->set(it, INVALID);
    1037         typename AuxGraph::Node node = _aux_graph->addNode();
    1038         _first->set(node, it);
    1039         _last->set(node, it);
    1040         ref.set(it, node);
    1041       }
    1042 
    1043       for (typename Graph::UEdgeIt it(*_graph); it != INVALID; ++it) {
    1044         if (_graph->source(it) == _graph->target(it)) continue;
    1045         typename AuxGraph::UEdge uedge =
    1046           _aux_graph->addEdge(ref[_graph->source(it)],
    1047                                ref[_graph->target(it)]);
    1048         _aux_capacity->set(uedge, (*_capacity)[it]);
    1049        
    1050       }
    10511050
    10521051      if (!_heap_cross_ref) {
     
    10591058      }
    10601059    }
     1060
     1061    void createAuxGraph() {
     1062      typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
     1063
     1064      for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
     1065        _next->set(n, INVALID);
     1066        typename AuxGraph::Node node = _aux_graph->addNode();
     1067        _first->set(node, n);
     1068        _last->set(node, n);
     1069        ref.set(n, node);
     1070      }
     1071
     1072      typename AuxGraph::template NodeMap<typename AuxGraph::UEdge>
     1073      edges(*_aux_graph, INVALID);
     1074
     1075      for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
     1076        for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
     1077          typename Graph::Node tn = _graph->runningNode(e);
     1078          if (n < tn || n == tn) continue;
     1079          if (edges[ref[tn]] != INVALID) {
     1080            Value value =
     1081              (*_aux_capacity)[edges[ref[tn]]] + (*_capacity)[e];
     1082            _aux_capacity->set(edges[ref[tn]], value);
     1083          } else {
     1084            edges.set(ref[tn], _aux_graph->addEdge(ref[n], ref[tn]));
     1085            Value value = (*_capacity)[e];
     1086            _aux_capacity->set(edges[ref[tn]], value);           
     1087          }
     1088        }
     1089        for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
     1090          typename Graph::Node tn = _graph->runningNode(e);
     1091          edges.set(ref[tn], INVALID);
     1092        }
     1093      }
     1094
     1095      _cut.resize(1, INVALID);
     1096      _min_cut = std::numeric_limits<Value>::max();
     1097      for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
     1098        Value value = 0;
     1099        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, n);
     1100             e != INVALID; ++e) {
     1101          value += (*_aux_capacity)[e];
     1102        }
     1103        if (_min_cut > value) {
     1104          _min_cut = value;
     1105          _cut[0] = (*_first)[n];
     1106        }
     1107        (*_aux_cut_value)[n] = value;
     1108      }
     1109    }
     1110   
    10611111
    10621112  public :
     
    10741124        _aux_graph(0), local_aux_graph(false),
    10751125        _aux_capacity(0), local_aux_capacity(false),
     1126        _aux_cut_value(0), local_aux_cut_value(false),
    10761127        _heap_cross_ref(0), local_heap_cross_ref(false),
    10771128        _heap(0), local_heap(false),
     
    10911142        _aux_graph(0), local_aux_graph(false),
    10921143        _aux_capacity(0), local_aux_capacity(false),
     1144        _aux_cut_value(0), local_aux_cut_value(false),
    10931145        _heap_cross_ref(0), local_heap_cross_ref(false),
    10941146        _heap(0), local_heap(false),
     
    11051157      if (_next) delete _next;
    11061158      if (local_aux_capacity) delete _aux_capacity;
     1159      if (local_aux_cut_value) delete _aux_cut_value;
    11071160      if (local_aux_graph) delete _aux_graph;
    11081161      if (local_capacity) delete _capacity;
     
    11791232    /// Initializes the internal data structures.
    11801233    void init() {
    1181       create_structures();
    1182       _first_node = _last_node = INVALID;
    11831234      _node_num = countNodes(*_graph);
    1184     }
     1235      createStructures();
     1236      createAuxGraph();
     1237    }
     1238
     1239  private:
     1240
     1241    struct EdgeInfo {
     1242      typename AuxGraph::Node source, target;
     1243      Value capacity;
     1244    };
     1245   
     1246    struct EdgeInfoLess {
     1247      bool operator()(const EdgeInfo& left, const EdgeInfo& right) const {
     1248        return (left.source < right.source) ||
     1249          (left.source == right.source && left.target < right.target);
     1250      }
     1251    };
     1252
     1253  public:
     1254
    11851255
    11861256    /// \brief Processes the next phase
    11871257    ///
    1188     /// Processes the next phase in the algorithm. The function
    1189     /// should be called countNodes(graph) - 1 times to get
    1190     /// surely the min cut in the graph. The 
     1258    /// Processes the next phase in the algorithm. The function should
     1259    /// be called at most countNodes(graph) - 1 times to get surely
     1260    /// the min cut in the graph.
    11911261    ///
    11921262    ///\return %True when the algorithm finished.
    11931263    bool processNextPhase() {
    11941264      if (_node_num <= 1) return true;
    1195       using namespace _min_cut_bits;
    11961265
    11971266      typedef typename AuxGraph::Node Node;
    11981267      typedef typename AuxGraph::NodeIt NodeIt;
    11991268      typedef typename AuxGraph::UEdge UEdge;
     1269      typedef typename AuxGraph::UEdgeIt UEdgeIt;
    12001270      typedef typename AuxGraph::IncEdgeIt IncEdgeIt;
    12011271     
    1202       typedef typename MaxCardinalitySearch<AuxGraph, AuxCapacityMap>::
    1203       template DefHeap<Heap, HeapCrossRef>::
    1204       template DefCardinalityMap<NullMap<Node, Value> >::
    1205       template DefProcessedMap<LastTwoMap<Node> >::
    1206       Create MaxCardinalitySearch;
     1272      typename AuxGraph::template UEdgeMap<Value> _edge_cut(*_aux_graph);
     1273
     1274
     1275      for (NodeIt n(*_aux_graph); n != INVALID; ++n) {
     1276        _heap_cross_ref->set(n, Heap::PRE_HEAP);
     1277      }
     1278
     1279      std::vector<Node> nodes;
     1280      nodes.reserve(_node_num);
     1281      int sep = 0;
     1282
     1283      Value alpha = 0;
     1284      Value emc = std::numeric_limits<Value>::max();
     1285
     1286      _heap->push(typename AuxGraph::NodeIt(*_aux_graph), 0);
     1287      while (!_heap->empty()) {
     1288        Node node = _heap->top();
     1289        Value value = _heap->prio();
     1290       
     1291        _heap->pop();
     1292        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, node);
     1293             e != INVALID; ++e) {
     1294          Node tn = _aux_graph->runningNode(e);
     1295          switch (_heap->state(tn)) {
     1296          case Heap::PRE_HEAP:
     1297            _heap->push(tn, (*_aux_capacity)[e]);
     1298            _edge_cut[e] = (*_heap)[tn];
     1299            break;
     1300          case Heap::IN_HEAP:
     1301            _heap->decrease(tn, (*_aux_capacity)[e] + (*_heap)[tn]);
     1302            _edge_cut[e] = (*_heap)[tn];
     1303            break;
     1304          case Heap::POST_HEAP:
     1305            break;
     1306          }
     1307        }
     1308
     1309        alpha += (*_aux_cut_value)[node];
     1310        alpha -= 2 * value;
     1311
     1312        nodes.push_back(node);
     1313        if (!_heap->empty()) {
     1314          if (alpha < emc) {
     1315            emc = alpha;
     1316            sep = nodes.size();
     1317          }
     1318        }
     1319      }
     1320
     1321      if ((int)nodes.size() < _node_num) {
     1322        _aux_graph->clear();
     1323        _node_num = 1;
     1324        _cut.clear();
     1325        for (int i = 0; i < (int)nodes.size(); ++i) {
     1326          typename Graph::Node n = (*_first)[nodes[i]];
     1327          while (n != INVALID) {
     1328            _cut.push_back(n);
     1329            n = (*_next)[n];
     1330          }
     1331        }
     1332        _min_cut = 0;
     1333        return true;
     1334      }
     1335
     1336      if (emc < _min_cut) {
     1337        _cut.clear();
     1338        for (int i = 0; i < sep; ++i) {
     1339          typename Graph::Node n = (*_first)[nodes[i]];
     1340          while (n != INVALID) {
     1341            _cut.push_back(n);
     1342            n = (*_next)[n];
     1343          }
     1344        }
     1345        _min_cut = emc;
     1346      }
     1347
     1348      typedef typename AuxGraph::template NodeMap<int> UfeCr;
     1349      UfeCr ufecr(*_aux_graph);
     1350      typedef UnionFindEnum<UfeCr> Ufe;
     1351      Ufe ufe(ufecr);
     1352
     1353      for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
     1354        ufe.insert(n);
     1355      }
     1356
     1357      for (typename AuxGraph::UEdgeIt e(*_aux_graph); e != INVALID; ++e) {
     1358        if (_edge_cut[e] >= emc) {
     1359          ufe.join(_aux_graph->source(e), _aux_graph->target(e));
     1360        }
     1361      }
     1362
     1363      if (ufe.size((typename Ufe::ClassIt)(ufe)) == _node_num) {
     1364        _aux_graph->clear();
     1365        _node_num = 1;
     1366        return true;
     1367      }
    12071368     
    1208       MaxCardinalitySearch mcs(*_aux_graph, *_aux_capacity);
    1209       for (NodeIt it(*_aux_graph); it != INVALID; ++it) {
    1210         _heap_cross_ref->set(it, Heap::PRE_HEAP);
    1211       }
    1212       mcs.heap(*_heap, *_heap_cross_ref);
    1213 
    1214       LastTwoMap<Node> last_two_nodes(_node_num);
    1215       mcs.processedMap(last_two_nodes);
    1216 
    1217       NullMap<Node, Value> cardinality;
    1218       mcs.cardinalityMap(cardinality);
    1219 
    1220       mcs.run();
    1221 
    1222       Node new_node = _aux_graph->addNode();
     1369      std::vector<typename AuxGraph::Node> remnodes;
    12231370
    12241371      typename AuxGraph::template NodeMap<UEdge> edges(*_aux_graph, INVALID);
    1225      
    1226       Node first_node = last_two_nodes[0];
    1227       Node second_node = last_two_nodes[1];
    1228      
    1229       _next->set((*_last)[first_node], (*_first)[second_node]);
    1230       _first->set(new_node, (*_first)[first_node]);
    1231       _last->set(new_node, (*_last)[second_node]);
    1232 
    1233       Value current_cut = 0;     
    1234       for (IncEdgeIt it(*_aux_graph, first_node); it != INVALID; ++it) {
    1235         Node node = _aux_graph->runningNode(it);
    1236         current_cut += (*_aux_capacity)[it];
    1237         if (node == second_node) continue;
    1238         if (edges[node] == INVALID) {
    1239           edges[node] = _aux_graph->addEdge(new_node, node);
    1240           (*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
     1372      for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
     1373        if (ufe.size(c) == 1) continue;
     1374        for (typename Ufe::ItemIt r(ufe, c); r != INVALID; ++r) {
     1375          if ((Node)r == (Node)c) continue;
     1376          _next->set((*_last)[c], (*_first)[r]);
     1377          _last->set(c, (*_last)[r]);
     1378          remnodes.push_back(r);
     1379          --_node_num;
     1380        }
     1381      }
     1382
     1383      std::vector<EdgeInfo> addedges;
     1384      std::vector<UEdge> remedges;
     1385
     1386      for (typename AuxGraph::UEdgeIt e(*_aux_graph);
     1387           e != INVALID; ++e) {
     1388        Node sn = ufe.find(_aux_graph->source(e));
     1389        Node tn = ufe.find(_aux_graph->target(e));
     1390        if ((ufe.size(sn) == 1 && ufe.size(tn) == 1)) {
     1391          continue;
     1392        }
     1393        if (sn == tn) {
     1394          remedges.push_back(e);
     1395          continue;
     1396        }
     1397        EdgeInfo info;
     1398        if (sn < tn) {
     1399          info.source = sn;
     1400          info.target = tn;
    12411401        } else {
    1242           (*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];         
     1402          info.source = tn;
     1403          info.target = sn;
    12431404        }
    1244       }
    1245 
    1246       if (_first_node == INVALID || current_cut < _min_cut) {
    1247         _first_node = (*_first)[first_node];
    1248         _last_node = (*_last)[first_node];
    1249         _min_cut = current_cut;
    1250       }
    1251 
    1252       _aux_graph->erase(first_node);
    1253 
    1254       for (IncEdgeIt it(*_aux_graph, second_node); it != INVALID; ++it) {
    1255         Node node = _aux_graph->runningNode(it);
    1256         if (edges[node] == INVALID) {
    1257           edges[node] = _aux_graph->addEdge(new_node, node);
    1258           (*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
    1259         } else {
    1260           (*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];         
     1405        info.capacity = (*_aux_capacity)[e];
     1406        addedges.push_back(info);
     1407        remedges.push_back(e);
     1408      }
     1409
     1410      for (int i = 0; i < (int)remedges.size(); ++i) {
     1411        _aux_graph->erase(remedges[i]);
     1412      }
     1413
     1414      sort(addedges.begin(), addedges.end(), EdgeInfoLess());
     1415
     1416      {
     1417        int i = 0;
     1418        while (i < (int)addedges.size()) {
     1419          Node sn = addedges[i].source;
     1420          Node tn = addedges[i].target;
     1421          Value ec = addedges[i].capacity;
     1422          ++i;
     1423          while (i < (int)addedges.size() &&
     1424                 sn == addedges[i].source && tn == addedges[i].target) {
     1425            ec += addedges[i].capacity;
     1426            ++i;
     1427          }
     1428          typename AuxGraph::UEdge ne = _aux_graph->addEdge(sn, tn);
     1429          (*_aux_capacity)[ne] = ec;
    12611430        }
    12621431      }
    1263       _aux_graph->erase(second_node);
    1264 
    1265       --_node_num;
     1432
     1433      for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
     1434        if (ufe.size(c) == 1) continue;
     1435        Value cutvalue = 0;
     1436        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, c);
     1437             e != INVALID; ++e) {
     1438          cutvalue += (*_aux_capacity)[e];
     1439        }
     1440       
     1441        (*_aux_cut_value)[c] = cutvalue;
     1442       
     1443      }
     1444
     1445      for (int i = 0; i < (int)remnodes.size(); ++i) {
     1446        _aux_graph->erase(remnodes[i]);
     1447      }
     1448
    12661449      return _node_num == 1;
    12671450    }
     
    13181501    template <typename NodeMap>
    13191502    Value quickMinCut(NodeMap& nodeMap) const {
    1320       for (typename Graph::Node it = _first_node;
    1321            it != _last_node; it = (*_next)[it]) {
    1322              nodeMap.set(it, true);
    1323            }
    1324       nodeMap.set(_last_node, true);
     1503      for (int i = 0; i < (int)_cut.size(); ++i) {
     1504        nodeMap.set(_cut[i], true);
     1505      }
    13251506      return minCut();
    13261507    }
     
    13561537    ///@}
    13571538
     1539  private:
     1540
     1541
    13581542  };
    13591543}
Note: See TracChangeset for help on using the changeset viewer.