COIN-OR::LEMON - Graph Library

Changeset 327:91d63b8b1a4c in lemon-1.2


Ignore:
Timestamp:
10/13/08 14:00:11 (11 years ago)
Author:
Balazs Dezso <deba@…>
Branch:
default
Phase:
public
Message:

Several improvements in maximum matching algorithms

  • The interface of MaxMatching? is changed to be similar to the weighted algorithms
  • The internal data structure (the queue implementation and the matching map) is changed in the MaxMatching? algorithm, which provides better runtime properties
  • The Blossom iterators are changed slightly in the weighted matching algorithms
  • Several documentation improvments
  • The test files are merged
Files:
1 deleted
4 edited

Legend:

Unmodified
Added
Removed
  • lemon/max_matching.h

    r326 r327  
    3232///\ingroup matching
    3333///\file
    34 ///\brief Maximum matching algorithms in graph.
     34///\brief Maximum matching algorithms in general graphs.
    3535
    3636namespace lemon {
    3737
    38   ///\ingroup matching
     38  /// \ingroup matching
    3939  ///
    40   ///\brief Edmonds' alternating forest maximum matching algorithm.
     40  /// \brief Edmonds' alternating forest maximum matching algorithm.
    4141  ///
    42   ///This class provides Edmonds' alternating forest matching
    43   ///algorithm. The starting matching (if any) can be passed to the
    44   ///algorithm using some of init functions.
     42  /// This class provides Edmonds' alternating forest matching
     43  /// algorithm. The starting matching (if any) can be passed to the
     44  /// algorithm using some of init functions.
    4545  ///
    46   ///The dual side of a matching is a map of the nodes to
    47   ///MaxMatching::DecompType, having values \c D, \c A and \c C
    48   ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
    49   ///in \c D induce a digraph with factor-critical components, the nodes
    50   ///in \c A form the barrier, and the nodes in \c C induce a digraph
    51   ///having a perfect matching. This decomposition can be attained by
    52   ///calling \c decomposition() after running the algorithm.
     46  /// The dual side of a matching is a map of the nodes to
     47  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
     48  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
     49  /// graph. The nodes in \c EVEN/D induce a graph with
     50  /// factor-critical components, the nodes in \c ODD/A form the
     51  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
     52  /// perfect matching. The number of the fractor critical components
     53  /// minus the number of barrier nodes is a lower bound on the
     54  /// unmatched nodes, and if the matching is optimal this bound is
     55  /// tight. This decomposition can be attained by calling \c
     56  /// decomposition() after running the algorithm.
    5357  ///
    54   ///\param Digraph The graph type the algorithm runs on.
    55   template <typename Graph>
     58  /// \param _Graph The graph type the algorithm runs on.
     59  template <typename _Graph>
    5660  class MaxMatching {
    57 
    58   protected:
     61  public:
     62
     63    typedef _Graph Graph;
     64    typedef typename Graph::template NodeMap<typename Graph::Arc>
     65    MatchingMap;
     66
     67    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
     68    ///
     69    ///Indicates the Gallai-Edmonds decomposition of the graph, which
     70    ///shows an upper bound on the size of a maximum matching. The
     71    ///nodes with Status \c EVEN/D induce a graph with factor-critical
     72    ///components, the nodes in \c ODD/A form the canonical barrier,
     73    ///and the nodes in \c MATCHED/C induce a graph having a perfect
     74    ///matching.
     75    enum Status {
     76      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
     77    };
     78
     79    typedef typename Graph::template NodeMap<Status> StatusMap;
     80
     81  private:
    5982
    6083    TEMPLATE_GRAPH_TYPEDEFS(Graph);
    6184
    62     typedef typename Graph::template NodeMap<int> UFECrossRef;
    63     typedef UnionFindEnum<UFECrossRef> UFE;
    64     typedef std::vector<Node> NV;
    65 
    66     typedef typename Graph::template NodeMap<int> EFECrossRef;
    67     typedef ExtendFindEnum<EFECrossRef> EFE;
     85    typedef UnionFindEnum<IntNodeMap> BlossomSet;
     86    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     87    typedef RangeMap<Node> NodeIntMap;
     88    typedef MatchingMap EarMap;
     89    typedef std::vector<Node> NodeQueue;
     90
     91    const Graph& _graph;
     92    MatchingMap* _matching;
     93    StatusMap* _status;
     94
     95    EarMap* _ear;
     96
     97    IntNodeMap* _blossom_set_index;
     98    BlossomSet* _blossom_set;
     99    NodeIntMap* _blossom_rep;
     100
     101    IntNodeMap* _tree_set_index;
     102    TreeSet* _tree_set;
     103
     104    NodeQueue _node_queue;
     105    int _process, _postpone, _last;
     106
     107    int _node_num;
     108
     109  private:
     110
     111    void createStructures() {
     112      _node_num = countNodes(_graph);
     113      if (!_matching) {
     114        _matching = new MatchingMap(_graph);
     115      }
     116      if (!_status) {
     117        _status = new StatusMap(_graph);
     118      }
     119      if (!_ear) {
     120        _ear = new EarMap(_graph);
     121      }
     122      if (!_blossom_set) {
     123        _blossom_set_index = new IntNodeMap(_graph);
     124        _blossom_set = new BlossomSet(*_blossom_set_index);
     125      }
     126      if (!_blossom_rep) {
     127        _blossom_rep = new NodeIntMap(_node_num);
     128      }
     129      if (!_tree_set) {
     130        _tree_set_index = new IntNodeMap(_graph);
     131        _tree_set = new TreeSet(*_tree_set_index);
     132      }
     133      _node_queue.resize(_node_num);
     134    }
     135
     136    void destroyStructures() {
     137      if (_matching) {
     138        delete _matching;
     139      }
     140      if (_status) {
     141        delete _status;
     142      }
     143      if (_ear) {
     144        delete _ear;
     145      }
     146      if (_blossom_set) {
     147        delete _blossom_set;
     148        delete _blossom_set_index;
     149      }
     150      if (_blossom_rep) {
     151        delete _blossom_rep;
     152      }
     153      if (_tree_set) {
     154        delete _tree_set_index;
     155        delete _tree_set;
     156      }
     157    }
     158
     159    void processDense(const Node& n) {
     160      _process = _postpone = _last = 0;
     161      _node_queue[_last++] = n;
     162
     163      while (_process != _last) {
     164        Node u = _node_queue[_process++];
     165        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     166          Node v = _graph.target(a);
     167          if ((*_status)[v] == MATCHED) {
     168            extendOnArc(a);
     169          } else if ((*_status)[v] == UNMATCHED) {
     170            augmentOnArc(a);
     171            return;
     172          }
     173        }
     174      }
     175
     176      while (_postpone != _last) {
     177        Node u = _node_queue[_postpone++];
     178
     179        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
     180          Node v = _graph.target(a);
     181
     182          if ((*_status)[v] == EVEN) {
     183            if (_blossom_set->find(u) != _blossom_set->find(v)) {
     184              shrinkOnEdge(a);
     185            }
     186          }
     187
     188          while (_process != _last) {
     189            Node w = _node_queue[_process++];
     190            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
     191              Node x = _graph.target(b);
     192              if ((*_status)[x] == MATCHED) {
     193                extendOnArc(b);
     194              } else if ((*_status)[x] == UNMATCHED) {
     195                augmentOnArc(b);
     196                return;
     197              }
     198            }
     199          }
     200        }
     201      }
     202    }
     203
     204    void processSparse(const Node& n) {
     205      _process = _last = 0;
     206      _node_queue[_last++] = n;
     207      while (_process != _last) {
     208        Node u = _node_queue[_process++];
     209        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     210          Node v = _graph.target(a);
     211
     212          if ((*_status)[v] == EVEN) {
     213            if (_blossom_set->find(u) != _blossom_set->find(v)) {
     214              shrinkOnEdge(a);
     215            }
     216          } else if ((*_status)[v] == MATCHED) {
     217            extendOnArc(a);
     218          } else if ((*_status)[v] == UNMATCHED) {
     219            augmentOnArc(a);
     220            return;
     221          }
     222        }
     223      }
     224    }
     225
     226    void shrinkOnEdge(const Edge& e) {
     227      Node nca = INVALID;
     228
     229      {
     230        std::set<Node> left_set, right_set;
     231
     232        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
     233        left_set.insert(left);
     234
     235        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
     236        right_set.insert(right);
     237
     238        while (true) {
     239          if ((*_matching)[left] == INVALID) break;
     240          left = _graph.target((*_matching)[left]);
     241          left = (*_blossom_rep)[_blossom_set->
     242                                 find(_graph.target((*_ear)[left]))];
     243          if (right_set.find(left) != right_set.end()) {
     244            nca = left;
     245            break;
     246          }
     247          left_set.insert(left);
     248
     249          if ((*_matching)[right] == INVALID) break;
     250          right = _graph.target((*_matching)[right]);
     251          right = (*_blossom_rep)[_blossom_set->
     252                                  find(_graph.target((*_ear)[right]))];
     253          if (left_set.find(right) != left_set.end()) {
     254            nca = right;
     255            break;
     256          }
     257          right_set.insert(right);
     258        }
     259
     260        if (nca == INVALID) {
     261          if ((*_matching)[left] == INVALID) {
     262            nca = right;
     263            while (left_set.find(nca) == left_set.end()) {
     264              nca = _graph.target((*_matching)[nca]);
     265              nca =(*_blossom_rep)[_blossom_set->
     266                                   find(_graph.target((*_ear)[nca]))];
     267            }
     268          } else {
     269            nca = left;
     270            while (right_set.find(nca) == right_set.end()) {
     271              nca = _graph.target((*_matching)[nca]);
     272              nca = (*_blossom_rep)[_blossom_set->
     273                                   find(_graph.target((*_ear)[nca]))];
     274            }
     275          }
     276        }
     277      }
     278
     279      {
     280
     281        Node node = _graph.u(e);
     282        Arc arc = _graph.direct(e, true);
     283        Node base = (*_blossom_rep)[_blossom_set->find(node)];
     284
     285        while (base != nca) {
     286          _ear->set(node, arc);
     287
     288          Node n = node;
     289          while (n != base) {
     290            n = _graph.target((*_matching)[n]);
     291            Arc a = (*_ear)[n];
     292            n = _graph.target(a);
     293            _ear->set(n, _graph.oppositeArc(a));
     294          }
     295          node = _graph.target((*_matching)[base]);
     296          _tree_set->erase(base);
     297          _tree_set->erase(node);
     298          _blossom_set->insert(node, _blossom_set->find(base));
     299          _status->set(node, EVEN);
     300          _node_queue[_last++] = node;
     301          arc = _graph.oppositeArc((*_ear)[node]);
     302          node = _graph.target((*_ear)[node]);
     303          base = (*_blossom_rep)[_blossom_set->find(node)];
     304          _blossom_set->join(_graph.target(arc), base);
     305        }
     306      }
     307
     308      _blossom_rep->set(_blossom_set->find(nca), nca);
     309
     310      {
     311
     312        Node node = _graph.v(e);
     313        Arc arc = _graph.direct(e, false);
     314        Node base = (*_blossom_rep)[_blossom_set->find(node)];
     315
     316        while (base != nca) {
     317          _ear->set(node, arc);
     318
     319          Node n = node;
     320          while (n != base) {
     321            n = _graph.target((*_matching)[n]);
     322            Arc a = (*_ear)[n];
     323            n = _graph.target(a);
     324            _ear->set(n, _graph.oppositeArc(a));
     325          }
     326          node = _graph.target((*_matching)[base]);
     327          _tree_set->erase(base);
     328          _tree_set->erase(node);
     329          _blossom_set->insert(node, _blossom_set->find(base));
     330          _status->set(node, EVEN);
     331          _node_queue[_last++] = node;
     332          arc = _graph.oppositeArc((*_ear)[node]);
     333          node = _graph.target((*_ear)[node]);
     334          base = (*_blossom_rep)[_blossom_set->find(node)];
     335          _blossom_set->join(_graph.target(arc), base);
     336        }
     337      }
     338
     339      _blossom_rep->set(_blossom_set->find(nca), nca);
     340    }
     341
     342
     343
     344    void extendOnArc(const Arc& a) {
     345      Node base = _graph.source(a);
     346      Node odd = _graph.target(a);
     347
     348      _ear->set(odd, _graph.oppositeArc(a));
     349      Node even = _graph.target((*_matching)[odd]);
     350      _blossom_rep->set(_blossom_set->insert(even), even);
     351      _status->set(odd, ODD);
     352      _status->set(even, EVEN);
     353      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
     354      _tree_set->insert(odd, tree);
     355      _tree_set->insert(even, tree);
     356      _node_queue[_last++] = even;
     357
     358    }
     359
     360    void augmentOnArc(const Arc& a) {
     361      Node even = _graph.source(a);
     362      Node odd = _graph.target(a);
     363
     364      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
     365
     366      _matching->set(odd, _graph.oppositeArc(a));
     367      _status->set(odd, MATCHED);
     368
     369      Arc arc = (*_matching)[even];
     370      _matching->set(even, a);
     371
     372      while (arc != INVALID) {
     373        odd = _graph.target(arc);
     374        arc = (*_ear)[odd];
     375        even = _graph.target(arc);
     376        _matching->set(odd, arc);
     377        arc = (*_matching)[even];
     378        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
     379      }
     380
     381      for (typename TreeSet::ItemIt it(*_tree_set, tree);
     382           it != INVALID; ++it) {
     383        if ((*_status)[it] == ODD) {
     384          _status->set(it, MATCHED);
     385        } else {
     386          int blossom = _blossom_set->find(it);
     387          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
     388               jt != INVALID; ++jt) {
     389            _status->set(jt, MATCHED);
     390          }
     391          _blossom_set->eraseClass(blossom);
     392        }
     393      }
     394      _tree_set->eraseClass(tree);
     395
     396    }
    68397
    69398  public:
    70399
    71     ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
    72     ///
    73     ///Indicates the Gallai-Edmonds decomposition of the digraph, which
    74     ///shows an upper bound on the size of a maximum matching. The
    75     ///nodes with DecompType \c D induce a digraph with factor-critical
    76     ///components, the nodes in \c A form the canonical barrier, and the
    77     ///nodes in \c C induce a digraph having a perfect matching.
    78     enum DecompType {
    79       D=0,
    80       A=1,
    81       C=2
    82     };
    83 
    84   protected:
    85 
    86     static const int HEUR_density=2;
    87     const Graph& g;
    88     typename Graph::template NodeMap<Node> _mate;
    89     typename Graph::template NodeMap<DecompType> position;
    90 
    91   public:
    92 
    93     MaxMatching(const Graph& _g)
    94       : g(_g), _mate(_g), position(_g) {}
    95 
    96     ///\brief Sets the actual matching to the empty matching.
    97     ///
    98     ///Sets the actual matching to the empty matching.
     400    /// \brief Constructor
     401    ///
     402    /// Constructor.
     403    MaxMatching(const Graph& graph)
     404      : _graph(graph), _matching(0), _status(0), _ear(0),
     405        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
     406        _tree_set_index(0), _tree_set(0) {}
     407
     408    ~MaxMatching() {
     409      destroyStructures();
     410    }
     411
     412    /// \name Execution control
     413    /// The simplest way to execute the algorithm is to use the member
     414    /// \c run() member function.
     415    /// \n
     416
     417    /// If you need more control on the execution, first you must call
     418    /// \ref init(), \ref greedyInit() or \ref matchingInit()
     419    /// functions, then you can start the algorithm with the \ref
     420    /// startParse() or startDense() functions.
     421
     422    ///@{
     423
     424    /// \brief Sets the actual matching to the empty matching.
     425    ///
     426    /// Sets the actual matching to the empty matching.
    99427    ///
    100428    void init() {
    101       for(NodeIt v(g); v!=INVALID; ++v) {
    102         _mate.set(v,INVALID);
    103         position.set(v,C);
     429      createStructures();
     430      for(NodeIt n(_graph); n != INVALID; ++n) {
     431        _matching->set(n, INVALID);
     432        _status->set(n, UNMATCHED);
    104433      }
    105434    }
     
    109438    ///For initial matchig it finds a maximal greedy matching.
    110439    void greedyInit() {
    111       for(NodeIt v(g); v!=INVALID; ++v) {
    112         _mate.set(v,INVALID);
    113         position.set(v,C);
    114       }
    115       for(NodeIt v(g); v!=INVALID; ++v)
    116         if ( _mate[v]==INVALID ) {
    117           for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
    118             Node y=g.runningNode(e);
    119             if ( _mate[y]==INVALID && y!=v ) {
    120               _mate.set(v,y);
    121               _mate.set(y,v);
     440      createStructures();
     441      for (NodeIt n(_graph); n != INVALID; ++n) {
     442        _matching->set(n, INVALID);
     443        _status->set(n, UNMATCHED);
     444      }
     445      for (NodeIt n(_graph); n != INVALID; ++n) {
     446        if ((*_matching)[n] == INVALID) {
     447          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
     448            Node v = _graph.target(a);
     449            if ((*_matching)[v] == INVALID && v != n) {
     450              _matching->set(n, a);
     451              _status->set(n, MATCHED);
     452              _matching->set(v, _graph.oppositeArc(a));
     453              _status->set(v, MATCHED);
    122454              break;
    123455            }
    124456          }
    125457        }
    126     }
    127 
    128     ///\brief Initialize the matching from each nodes' mate.
    129     ///
    130     ///Initialize the matching from a \c Node valued \c Node map. This
    131     ///map must be \e symmetric, i.e. if \c map[u]==v then \c
    132     ///map[v]==u must hold, and \c uv will be an arc of the initial
    133     ///matching.
    134     template <typename MateMap>
    135     void mateMapInit(MateMap& map) {
    136       for(NodeIt v(g); v!=INVALID; ++v) {
    137         _mate.set(v,map[v]);
    138         position.set(v,C);
    139       }
    140     }
    141 
    142     ///\brief Initialize the matching from a node map with the
    143     ///incident matching arcs.
    144     ///
    145     ///Initialize the matching from an \c Edge valued \c Node map. \c
    146     ///map[v] must be an \c Edge incident to \c v. This map must have
    147     ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
    148     ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
    149     ///u to \c v will be an arc of the matching.
    150     template<typename MatchingMap>
    151     void matchingMapInit(MatchingMap& map) {
    152       for(NodeIt v(g); v!=INVALID; ++v) {
    153         position.set(v,C);
    154         Edge e=map[v];
    155         if ( e!=INVALID )
    156           _mate.set(v,g.oppositeNode(v,e));
    157         else
    158           _mate.set(v,INVALID);
    159       }
    160     }
    161 
    162     ///\brief Initialize the matching from the map containing the
    163     ///undirected matching arcs.
    164     ///
    165     ///Initialize the matching from a \c bool valued \c Edge map. This
    166     ///map must have the property that there are no two incident arcs
    167     ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
    168     ///map[e]==true form the matching.
     458      }
     459    }
     460
     461
     462    /// \brief Initialize the matching from the map containing a matching.
     463    ///
     464    /// Initialize the matching from a \c bool valued \c Edge map. This
     465    /// map must have the property that there are no two incident edges
     466    /// with true value, ie. it contains a matching.
     467    /// \return %True if the map contains a matching.
    169468    template <typename MatchingMap>
    170     void matchingInit(MatchingMap& map) {
    171       for(NodeIt v(g); v!=INVALID; ++v) {
    172         _mate.set(v,INVALID);
    173         position.set(v,C);
    174       }
    175       for(EdgeIt e(g); e!=INVALID; ++e) {
    176         if ( map[e] ) {
    177           Node u=g.u(e);
    178           Node v=g.v(e);
    179           _mate.set(u,v);
    180           _mate.set(v,u);
    181         }
    182       }
    183     }
    184 
    185 
    186     ///\brief Runs Edmonds' algorithm.
    187     ///
    188     ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
    189     ///2*number of nodes), and a heuristical Edmonds' algorithm with a
    190     ///heuristic of postponing shrinks for dense digraphs.
     469    bool matchingInit(const MatchingMap& matching) {
     470      createStructures();
     471
     472      for (NodeIt n(_graph); n != INVALID; ++n) {
     473        _matching->set(n, INVALID);
     474        _status->set(n, UNMATCHED);
     475      }
     476      for(EdgeIt e(_graph); e!=INVALID; ++e) {
     477        if (matching[e]) {
     478
     479          Node u = _graph.u(e);
     480          if ((*_matching)[u] != INVALID) return false;
     481          _matching->set(u, _graph.direct(e, true));
     482          _status->set(u, MATCHED);
     483
     484          Node v = _graph.v(e);
     485          if ((*_matching)[v] != INVALID) return false;
     486          _matching->set(v, _graph.direct(e, false));
     487          _status->set(v, MATCHED);
     488        }
     489      }
     490      return true;
     491    }
     492
     493    /// \brief Starts Edmonds' algorithm
     494    ///
     495    /// If runs the original Edmonds' algorithm.
     496    void startSparse() {
     497      for(NodeIt n(_graph); n != INVALID; ++n) {
     498        if ((*_status)[n] == UNMATCHED) {
     499          (*_blossom_rep)[_blossom_set->insert(n)] = n;
     500          _tree_set->insert(n);
     501          _status->set(n, EVEN);
     502          processSparse(n);
     503        }
     504      }
     505    }
     506
     507    /// \brief Starts Edmonds' algorithm.
     508    ///
     509    /// It runs Edmonds' algorithm with a heuristic of postponing
     510    /// shrinks, giving a faster algorithm for dense graphs.
     511    void startDense() {
     512      for(NodeIt n(_graph); n != INVALID; ++n) {
     513        if ((*_status)[n] == UNMATCHED) {
     514          (*_blossom_rep)[_blossom_set->insert(n)] = n;
     515          _tree_set->insert(n);
     516          _status->set(n, EVEN);
     517          processDense(n);
     518        }
     519      }
     520    }
     521
     522
     523    /// \brief Runs Edmonds' algorithm
     524    ///
     525    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
     526    /// or Edmonds' algorithm with a heuristic of
     527    /// postponing shrinks for dense graphs.
    191528    void run() {
    192       if (countEdges(g) < HEUR_density * countNodes(g)) {
     529      if (countEdges(_graph) < 2 * countNodes(_graph)) {
    193530        greedyInit();
    194531        startSparse();
     
    199536    }
    200537
    201 
    202     ///\brief Starts Edmonds' algorithm.
    203     ///
    204     ///If runs the original Edmonds' algorithm.
    205     void startSparse() {
    206 
    207       typename Graph::template NodeMap<Node> ear(g,INVALID);
    208       //undefined for the base nodes of the blossoms (i.e. for the
    209       //representative elements of UFE blossom) and for the nodes in C
    210 
    211       UFECrossRef blossom_base(g);
    212       UFE blossom(blossom_base);
    213       NV rep(countNodes(g));
    214 
    215       EFECrossRef tree_base(g);
    216       EFE tree(tree_base);
    217 
    218       //If these UFE's would be members of the class then also
    219       //blossom_base and tree_base should be a member.
    220 
    221       //We build only one tree and the other vertices uncovered by the
    222       //matching belong to C. (They can be considered as singleton
    223       //trees.) If this tree can be augmented or no more
    224       //grow/augmentation/shrink is possible then we return to this
    225       //"for" cycle.
    226       for(NodeIt v(g); v!=INVALID; ++v) {
    227         if (position[v]==C && _mate[v]==INVALID) {
    228           rep[blossom.insert(v)] = v;
    229           tree.insert(v);
    230           position.set(v,D);
    231           normShrink(v, ear, blossom, rep, tree);
    232         }
    233       }
    234     }
    235 
    236     ///\brief Starts Edmonds' algorithm.
    237     ///
    238     ///It runs Edmonds' algorithm with a heuristic of postponing
    239     ///shrinks, giving a faster algorithm for dense digraphs.
    240     void startDense() {
    241 
    242       typename Graph::template NodeMap<Node> ear(g,INVALID);
    243       //undefined for the base nodes of the blossoms (i.e. for the
    244       //representative elements of UFE blossom) and for the nodes in C
    245 
    246       UFECrossRef blossom_base(g);
    247       UFE blossom(blossom_base);
    248       NV rep(countNodes(g));
    249 
    250       EFECrossRef tree_base(g);
    251       EFE tree(tree_base);
    252 
    253       //If these UFE's would be members of the class then also
    254       //blossom_base and tree_base should be a member.
    255 
    256       //We build only one tree and the other vertices uncovered by the
    257       //matching belong to C. (They can be considered as singleton
    258       //trees.) If this tree can be augmented or no more
    259       //grow/augmentation/shrink is possible then we return to this
    260       //"for" cycle.
    261       for(NodeIt v(g); v!=INVALID; ++v) {
    262         if ( position[v]==C && _mate[v]==INVALID ) {
    263           rep[blossom.insert(v)] = v;
    264           tree.insert(v);
    265           position.set(v,D);
    266           lateShrink(v, ear, blossom, rep, tree);
    267         }
    268       }
    269     }
    270 
    271 
     538    /// @}
     539
     540    /// \name Primal solution
     541    /// Functions for get the primal solution, ie. the matching.
     542
     543    /// @{
    272544
    273545    ///\brief Returns the size of the actual matching stored.
    274546    ///
    275547    ///Returns the size of the actual matching stored. After \ref
    276     ///run() it returns the size of a maximum matching in the digraph.
    277     int size() const {
    278       int s=0;
    279       for(NodeIt v(g); v!=INVALID; ++v) {
    280         if ( _mate[v]!=INVALID ) {
    281           ++s;
    282         }
    283       }
    284       return s/2;
    285     }
    286 
     548    ///run() it returns the size of the maximum matching in the graph.
     549    int matchingSize() const {
     550      int size = 0;
     551      for (NodeIt n(_graph); n != INVALID; ++n) {
     552        if ((*_matching)[n] != INVALID) {
     553          ++size;
     554        }
     555      }
     556      return size / 2;
     557    }
     558
     559    /// \brief Returns true when the edge is in the matching.
     560    ///
     561    /// Returns true when the edge is in the matching.
     562    bool matching(const Edge& edge) const {
     563      return edge == (*_matching)[_graph.u(edge)];
     564    }
     565
     566    /// \brief Returns the matching edge incident to the given node.
     567    ///
     568    /// Returns the matching edge of a \c node in the actual matching or
     569    /// INVALID if the \c node is not covered by the actual matching.
     570    Arc matching(const Node& n) const {
     571      return (*_matching)[n];
     572    }
    287573
    288574    ///\brief Returns the mate of a node in the actual matching.
    289575    ///
    290     ///Returns the mate of a \c node in the actual matching.
    291     ///Returns INVALID if the \c node is not covered by the actual matching.
    292     Node mate(const Node& node) const {
    293       return _mate[node];
    294     }
    295 
    296     ///\brief Returns the matching arc incident to the given node.
    297     ///
    298     ///Returns the matching arc of a \c node in the actual matching.
    299     ///Returns INVALID if the \c node is not covered by the actual matching.
    300     Edge matchingArc(const Node& node) const {
    301       if (_mate[node] == INVALID) return INVALID;
    302       Node n = node < _mate[node] ? node : _mate[node];
    303       for (IncEdgeIt e(g, n); e != INVALID; ++e) {
    304         if (g.oppositeNode(n, e) == _mate[n]) {
    305           return e;
    306         }
    307       }
    308       return INVALID;
    309     }
     576    ///Returns the mate of a \c node in the actual matching or
     577    ///INVALID if the \c node is not covered by the actual matching.
     578    Node mate(const Node& n) const {
     579      return (*_matching)[n] != INVALID ?
     580        _graph.target((*_matching)[n]) : INVALID;
     581    }
     582
     583    /// @}
     584
     585    /// \name Dual solution
     586    /// Functions for get the dual solution, ie. the decomposition.
     587
     588    /// @{
    310589
    311590    /// \brief Returns the class of the node in the Edmonds-Gallai
     
    314593    /// Returns the class of the node in the Edmonds-Gallai
    315594    /// decomposition.
    316     DecompType decomposition(const Node& n) {
    317       return position[n] == A;
     595    Status decomposition(const Node& n) const {
     596      return (*_status)[n];
    318597    }
    319598
     
    321600    ///
    322601    /// Returns true when the node is in the barrier.
    323     bool barrier(const Node& n) {
    324       return position[n] == A;
    325     }
    326 
    327     ///\brief Gives back the matching in a \c Node of mates.
    328     ///
    329     ///Writes the stored matching to a \c Node valued \c Node map. The
    330     ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
    331     ///map[v]==u will hold, and now \c uv is an arc of the matching.
    332     template <typename MateMap>
    333     void mateMap(MateMap& map) const {
    334       for(NodeIt v(g); v!=INVALID; ++v) {
    335         map.set(v,_mate[v]);
    336       }
    337     }
    338 
    339     ///\brief Gives back the matching in an \c Edge valued \c Node
    340     ///map.
    341     ///
    342     ///Writes the stored matching to an \c Edge valued \c Node
    343     ///map. \c map[v] will be an \c Edge incident to \c v. This
    344     ///map will have the property that if \c g.oppositeNode(u,map[u])
    345     ///== v then \c map[u]==map[v] holds, and now this arc is an arc
    346     ///of the matching.
    347     template <typename MatchingMap>
    348     void matchingMap(MatchingMap& map)  const {
    349       typename Graph::template NodeMap<bool> todo(g,true);
    350       for(NodeIt v(g); v!=INVALID; ++v) {
    351         if (_mate[v]!=INVALID && v < _mate[v]) {
    352           Node u=_mate[v];
    353           for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
    354             if ( g.runningNode(e) == u ) {
    355               map.set(u,e);
    356               map.set(v,e);
    357               todo.set(u,false);
    358               todo.set(v,false);
    359               break;
    360             }
    361           }
    362         }
    363       }
    364     }
    365 
    366 
    367     ///\brief Gives back the matching in a \c bool valued \c Edge
    368     ///map.
    369     ///
    370     ///Writes the matching stored to a \c bool valued \c Arc
    371     ///map. This map will have the property that there are no two
    372     ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
    373     ///arcs \c e with \c map[e]==true form the matching.
    374     template<typename MatchingMap>
    375     void matching(MatchingMap& map) const {
    376       for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
    377 
    378       typename Graph::template NodeMap<bool> todo(g,true);
    379       for(NodeIt v(g); v!=INVALID; ++v) {
    380         if ( todo[v] && _mate[v]!=INVALID ) {
    381           Node u=_mate[v];
    382           for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
    383             if ( g.runningNode(e) == u ) {
    384               map.set(e,true);
    385               todo.set(u,false);
    386               todo.set(v,false);
    387               break;
    388             }
    389           }
    390         }
    391       }
    392     }
    393 
    394 
    395     ///\brief Returns the canonical decomposition of the digraph after running
    396     ///the algorithm.
    397     ///
    398     ///After calling any run methods of the class, it writes the
    399     ///Gallai-Edmonds canonical decomposition of the digraph. \c map
    400     ///must be a node map of \ref DecompType 's.
    401     template <typename DecompositionMap>
    402     void decomposition(DecompositionMap& map) const {
    403       for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
    404     }
    405 
    406     ///\brief Returns a barrier on the nodes.
    407     ///
    408     ///After calling any run methods of the class, it writes a
    409     ///canonical barrier on the nodes. The odd component number of the
    410     ///remaining digraph minus the barrier size is a lower bound for the
    411     ///uncovered nodes in the digraph. The \c map must be a node map of
    412     ///bools.
    413     template <typename BarrierMap>
    414     void barrier(BarrierMap& barrier) {
    415       for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
    416     }
    417 
    418   private:
    419 
    420 
    421     void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
    422                     UFE& blossom, NV& rep, EFE& tree) {
    423       //We have one tree which we grow, and also shrink but only if it
    424       //cannot be postponed. If we augment then we return to the "for"
    425       //cycle of runEdmonds().
    426 
    427       std::queue<Node> Q;   //queue of the totally unscanned nodes
    428       Q.push(v);
    429       std::queue<Node> R;
    430       //queue of the nodes which must be scanned for a possible shrink
    431 
    432       while ( !Q.empty() ) {
    433         Node x=Q.front();
    434         Q.pop();
    435         for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
    436           Node y=g.runningNode(e);
    437           //growOrAugment grows if y is covered by the matching and
    438           //augments if not. In this latter case it returns 1.
    439           if (position[y]==C &&
    440               growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
    441         }
    442         R.push(x);
    443       }
    444 
    445       while ( !R.empty() ) {
    446         Node x=R.front();
    447         R.pop();
    448 
    449         for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
    450           Node y=g.runningNode(e);
    451 
    452           if ( position[y] == D && blossom.find(x) != blossom.find(y) )
    453             //Recall that we have only one tree.
    454             shrink( x, y, ear, blossom, rep, tree, Q);
    455 
    456           while ( !Q.empty() ) {
    457             Node z=Q.front();
    458             Q.pop();
    459             for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
    460               Node w=g.runningNode(f);
    461               //growOrAugment grows if y is covered by the matching and
    462               //augments if not. In this latter case it returns 1.
    463               if (position[w]==C &&
    464                   growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
    465             }
    466             R.push(z);
    467           }
    468         } //for e
    469       } // while ( !R.empty() )
    470     }
    471 
    472     void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
    473                     UFE& blossom, NV& rep, EFE& tree) {
    474       //We have one tree, which we grow and shrink. If we augment then we
    475       //return to the "for" cycle of runEdmonds().
    476 
    477       std::queue<Node> Q;   //queue of the unscanned nodes
    478       Q.push(v);
    479       while ( !Q.empty() ) {
    480 
    481         Node x=Q.front();
    482         Q.pop();
    483 
    484         for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
    485           Node y=g.runningNode(e);
    486 
    487           switch ( position[y] ) {
    488           case D:          //x and y must be in the same tree
    489             if ( blossom.find(x) != blossom.find(y))
    490               //x and y are in the same tree
    491               shrink(x, y, ear, blossom, rep, tree, Q);
    492             break;
    493           case C:
    494             //growOrAugment grows if y is covered by the matching and
    495             //augments if not. In this latter case it returns 1.
    496             if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
    497             break;
    498           default: break;
    499           }
    500         }
    501       }
    502     }
    503 
    504     void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
    505                 UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
    506       //x and y are the two adjacent vertices in two blossoms.
    507 
    508       typename Graph::template NodeMap<bool> path(g,false);
    509 
    510       Node b=rep[blossom.find(x)];
    511       path.set(b,true);
    512       b=_mate[b];
    513       while ( b!=INVALID ) {
    514         b=rep[blossom.find(ear[b])];
    515         path.set(b,true);
    516         b=_mate[b];
    517       } //we go until the root through bases of blossoms and odd vertices
    518 
    519       Node top=y;
    520       Node middle=rep[blossom.find(top)];
    521       Node bottom=x;
    522       while ( !path[middle] )
    523         shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
    524       //Until we arrive to a node on the path, we update blossom, tree
    525       //and the positions of the odd nodes.
    526 
    527       Node base=middle;
    528       top=x;
    529       middle=rep[blossom.find(top)];
    530       bottom=y;
    531       Node blossom_base=rep[blossom.find(base)];
    532       while ( middle!=blossom_base )
    533         shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
    534       //Until we arrive to a node on the path, we update blossom, tree
    535       //and the positions of the odd nodes.
    536 
    537       rep[blossom.find(base)] = base;
    538     }
    539 
    540     void shrinkStep(Node& top, Node& middle, Node& bottom,
    541                     typename Graph::template NodeMap<Node>& ear,
    542                     UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
    543       //We traverse a blossom and update everything.
    544 
    545       ear.set(top,bottom);
    546       Node t=top;
    547       while ( t!=middle ) {
    548         Node u=_mate[t];
    549         t=ear[u];
    550         ear.set(t,u);
    551       }
    552       bottom=_mate[middle];
    553       position.set(bottom,D);
    554       Q.push(bottom);
    555       top=ear[bottom];
    556       Node oldmiddle=middle;
    557       middle=rep[blossom.find(top)];
    558       tree.erase(bottom);
    559       tree.erase(oldmiddle);
    560       blossom.insert(bottom);
    561       blossom.join(bottom, oldmiddle);
    562       blossom.join(top, oldmiddle);
    563     }
    564 
    565 
    566 
    567     bool growOrAugment(Node& y, Node& x, typename Graph::template
    568                        NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
    569                        std::queue<Node>& Q) {
    570       //x is in a blossom in the tree, y is outside. If y is covered by
    571       //the matching we grow, otherwise we augment. In this case we
    572       //return 1.
    573 
    574       if ( _mate[y]!=INVALID ) {       //grow
    575         ear.set(y,x);
    576         Node w=_mate[y];
    577         rep[blossom.insert(w)] = w;
    578         position.set(y,A);
    579         position.set(w,D);
    580         int t = tree.find(rep[blossom.find(x)]);
    581         tree.insert(y,t);
    582         tree.insert(w,t);
    583         Q.push(w);
    584       } else {                      //augment
    585         augment(x, ear, blossom, rep, tree);
    586         _mate.set(x,y);
    587         _mate.set(y,x);
    588         return true;
    589       }
    590       return false;
    591     }
    592 
    593     void augment(Node x, typename Graph::template NodeMap<Node>& ear,
    594                  UFE& blossom, NV& rep, EFE& tree) {
    595       Node v=_mate[x];
    596       while ( v!=INVALID ) {
    597 
    598         Node u=ear[v];
    599         _mate.set(v,u);
    600         Node tmp=v;
    601         v=_mate[u];
    602         _mate.set(u,tmp);
    603       }
    604       int y = tree.find(rep[blossom.find(x)]);
    605       for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
    606         if ( position[tit] == D ) {
    607           int b = blossom.find(tit);
    608           for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
    609             position.set(bit, C);
    610           }
    611           blossom.eraseClass(b);
    612         } else position.set(tit, C);
    613       }
    614       tree.eraseClass(y);
    615 
    616     }
     602    bool barrier(const Node& n) const {
     603      return (*_status)[n] == ODD;
     604    }
     605
     606    /// @}
    617607
    618608  };
     
    628618  ///
    629619  /// The maximum weighted matching problem is to find undirected
    630   /// arcs in the digraph with maximum overall weight and no two of
    631   /// them shares their endpoints. The problem can be formulated with
    632   /// the next linear program:
     620  /// edges in the graph with maximum overall weight and no two of
     621  /// them shares their ends. The problem can be formulated with the
     622  /// following linear program.
    633623  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
    634   ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
     624  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
     625      \quad \forall B\in\mathcal{O}\f] */
    635626  /// \f[x_e \ge 0\quad \forall e\in E\f]
    636627  /// \f[\max \sum_{e\in E}x_ew_e\f]
    637   /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
    638   /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
    639   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
    640   /// the nodes.
     628  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     629  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
     630  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
     631  /// subsets of the nodes.
    641632  ///
    642633  /// The algorithm calculates an optimal matching and a proof of the
    643634  /// optimality. The solution of the dual problem can be used to check
    644   /// the result of the algorithm. The dual linear problem is the next:
    645   /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
     635  /// the result of the algorithm. The dual linear problem is the
     636  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
     637      z_B \ge w_{uv} \quad \forall uv\in E\f] */
    646638  /// \f[y_u \ge 0 \quad \forall u \in V\f]
    647639  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
    648   /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
     640  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
     641      \frac{\vert B \vert - 1}{2}z_B\f] */
    649642  ///
    650643  /// The algorithm can be executed with \c run() or the \c init() and
     
    706699    int _blossom_num;
    707700
    708     typedef typename Graph::template NodeMap<int> NodeIntMap;
    709     typedef typename Graph::template ArcMap<int> ArcIntMap;
    710     typedef typename Graph::template EdgeMap<int> EdgeIntMap;
    711701    typedef RangeMap<int> IntIntMap;
    712702
     
    715705    };
    716706
    717     typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
     707    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
    718708    struct BlossomData {
    719709      int tree;
     
    724714    };
    725715
    726     NodeIntMap *_blossom_index;
     716    IntNodeMap *_blossom_index;
    727717    BlossomSet *_blossom_set;
    728718    RangeMap<BlossomData>* _blossom_data;
    729719
    730     NodeIntMap *_node_index;
    731     ArcIntMap *_node_heap_index;
     720    IntNodeMap *_node_index;
     721    IntArcMap *_node_heap_index;
    732722
    733723    struct NodeData {
    734724
    735       NodeData(ArcIntMap& node_heap_index)
     725      NodeData(IntArcMap& node_heap_index)
    736726        : heap(node_heap_index) {}
    737727
    738728      int blossom;
    739729      Value pot;
    740       BinHeap<Value, ArcIntMap> heap;
     730      BinHeap<Value, IntArcMap> heap;
    741731      std::map<int, Arc> heap_index;
    742732
     
    751741    TreeSet *_tree_set;
    752742
    753     NodeIntMap *_delta1_index;
    754     BinHeap<Value, NodeIntMap> *_delta1;
     743    IntNodeMap *_delta1_index;
     744    BinHeap<Value, IntNodeMap> *_delta1;
    755745
    756746    IntIntMap *_delta2_index;
    757747    BinHeap<Value, IntIntMap> *_delta2;
    758748
    759     EdgeIntMap *_delta3_index;
    760     BinHeap<Value, EdgeIntMap> *_delta3;
     749    IntEdgeMap *_delta3_index;
     750    BinHeap<Value, IntEdgeMap> *_delta3;
    761751
    762752    IntIntMap *_delta4_index;
     
    776766      }
    777767      if (!_blossom_set) {
    778         _blossom_index = new NodeIntMap(_graph);
     768        _blossom_index = new IntNodeMap(_graph);
    779769        _blossom_set = new BlossomSet(*_blossom_index);
    780770        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
     
    782772
    783773      if (!_node_index) {
    784         _node_index = new NodeIntMap(_graph);
    785         _node_heap_index = new ArcIntMap(_graph);
     774        _node_index = new IntNodeMap(_graph);
     775        _node_heap_index = new IntArcMap(_graph);
    786776        _node_data = new RangeMap<NodeData>(_node_num,
    787777                                              NodeData(*_node_heap_index));
     
    793783      }
    794784      if (!_delta1) {
    795         _delta1_index = new NodeIntMap(_graph);
    796         _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
     785        _delta1_index = new IntNodeMap(_graph);
     786        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
    797787      }
    798788      if (!_delta2) {
     
    801791      }
    802792      if (!_delta3) {
    803         _delta3_index = new EdgeIntMap(_graph);
    804         _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
     793        _delta3_index = new IntEdgeMap(_graph);
     794        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
    805795      }
    806796      if (!_delta4) {
     
    12671257
    12681258
    1269     void augmentOnArc(const Edge& arc) {
    1270 
    1271       int left = _blossom_set->find(_graph.u(arc));
    1272       int right = _blossom_set->find(_graph.v(arc));
     1259    void augmentOnEdge(const Edge& edge) {
     1260
     1261      int left = _blossom_set->find(_graph.u(edge));
     1262      int right = _blossom_set->find(_graph.v(edge));
    12731263
    12741264      if ((*_blossom_data)[left].status == EVEN) {
     
    12901280      }
    12911281
    1292       (*_blossom_data)[left].next = _graph.direct(arc, true);
    1293       (*_blossom_data)[right].next = _graph.direct(arc, false);
     1282      (*_blossom_data)[left].next = _graph.direct(edge, true);
     1283      (*_blossom_data)[right].next = _graph.direct(edge, false);
    12941284    }
    12951285
     
    13111301    }
    13121302
    1313     void shrinkOnArc(const Edge& edge, int tree) {
     1303    void shrinkOnEdge(const Edge& edge, int tree) {
    13141304      int nca = -1;
    13151305      std::vector<int> left_path, right_path;
     
    16531643
    16541644      for (ArcIt e(_graph); e != INVALID; ++e) {
    1655         _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
     1645        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
    16561646      }
    16571647      for (NodeIt n(_graph); n != INVALID; ++n) {
     
    17701760
    17711761              if (left_tree == right_tree) {
    1772                 shrinkOnArc(e, left_tree);
     1762                shrinkOnEdge(e, left_tree);
    17731763              } else {
    1774                 augmentOnArc(e);
     1764                augmentOnEdge(e);
    17751765                unmatched -= 2;
    17761766              }
     
    18191809    }
    18201810
    1821     /// \brief Returns true when the arc is in the matching.
    1822     ///
    1823     /// Returns true when the arc is in the matching.
    1824     bool matching(const Edge& arc) const {
    1825       return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
     1811    /// \brief Returns the cardinality of the matching.
     1812    ///
     1813    /// Returns the cardinality of the matching.
     1814    int matchingSize() const {
     1815      int num = 0;
     1816      for (NodeIt n(_graph); n != INVALID; ++n) {
     1817        if ((*_matching)[n] != INVALID) {
     1818          ++num;
     1819        }
     1820      }
     1821      return num /= 2;
     1822    }
     1823
     1824    /// \brief Returns true when the edge is in the matching.
     1825    ///
     1826    /// Returns true when the edge is in the matching.
     1827    bool matching(const Edge& edge) const {
     1828      return edge == (*_matching)[_graph.u(edge)];
    18261829    }
    18271830
     
    19141917      }
    19151918
    1916       /// \brief Invalid constructor.
    1917       ///
    1918       /// Invalid constructor.
    1919       BlossomIt(Invalid) : _index(-1) {}
    1920 
    19211919      /// \brief Conversion to node.
    19221920      ///
    19231921      /// Conversion to node.
    19241922      operator Node() const {
    1925         return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
     1923        return _algorithm->_blossom_node_list[_index];
    19261924      }
    19271925
     
    19311929      BlossomIt& operator++() {
    19321930        ++_index;
    1933         if (_index == _last) {
    1934           _index = -1;
    1935         }
    19361931        return *this;
    19371932      }
    19381933
    1939       bool operator==(const BlossomIt& it) const {
    1940         return _index == it._index;
    1941       }
    1942       bool operator!=(const BlossomIt& it) const {
    1943         return _index != it._index;
    1944       }
     1934      /// \brief Validity checking
     1935      ///
     1936      /// Checks whether the iterator is invalid.
     1937      bool operator==(Invalid) const { return _index == _last; }
     1938
     1939      /// \brief Validity checking
     1940      ///
     1941      /// Checks whether the iterator is valid.
     1942      bool operator!=(Invalid) const { return _index != _last; }
    19451943
    19461944    private:
     
    19591957  ///
    19601958  /// This class provides an efficient implementation of Edmond's
    1961   /// maximum weighted perfecr matching algorithm. The implementation
     1959  /// maximum weighted perfect matching algorithm. The implementation
    19621960  /// is based on extensive use of priority queues and provides
    19631961  /// \f$O(nm\log(n))\f$ time complexity.
    19641962  ///
    19651963  /// The maximum weighted matching problem is to find undirected
    1966   /// arcs in the digraph with maximum overall weight and no two of
    1967   /// them shares their endpoints and covers all nodes. The problem
    1968   /// can be formulated with the next linear program:
     1964  /// edges in the graph with maximum overall weight and no two of
     1965  /// them shares their ends and covers all nodes. The problem can be
     1966  /// formulated with the following linear program.
    19691967  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
    1970   ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
     1968  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
     1969      \quad \forall B\in\mathcal{O}\f] */
    19711970  /// \f[x_e \ge 0\quad \forall e\in E\f]
    19721971  /// \f[\max \sum_{e\in E}x_ew_e\f]
    1973   /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
    1974   /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
    1975   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
    1976   /// the nodes.
     1972  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     1973  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
     1974  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
     1975  /// subsets of the nodes.
    19771976  ///
    19781977  /// The algorithm calculates an optimal matching and a proof of the
    19791978  /// optimality. The solution of the dual problem can be used to check
    1980   /// the result of the algorithm. The dual linear problem is the next:
    1981   /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
     1979  /// the result of the algorithm. The dual linear problem is the
     1980  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
     1981      w_{uv} \quad \forall uv\in E\f] */
    19821982  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
    1983   /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
     1983  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
     1984      \frac{\vert B \vert - 1}{2}z_B\f] */
    19841985  ///
    19851986  /// The algorithm can be executed with \c run() or the \c init() and
     
    20412042    int _blossom_num;
    20422043
    2043     typedef typename Graph::template NodeMap<int> NodeIntMap;
    2044     typedef typename Graph::template ArcMap<int> ArcIntMap;
    2045     typedef typename Graph::template EdgeMap<int> EdgeIntMap;
    20462044    typedef RangeMap<int> IntIntMap;
    20472045
     
    20502048    };
    20512049
    2052     typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
     2050    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
    20532051    struct BlossomData {
    20542052      int tree;
     
    20582056    };
    20592057
    2060     NodeIntMap *_blossom_index;
     2058    IntNodeMap *_blossom_index;
    20612059    BlossomSet *_blossom_set;
    20622060    RangeMap<BlossomData>* _blossom_data;
    20632061
    2064     NodeIntMap *_node_index;
    2065     ArcIntMap *_node_heap_index;
     2062    IntNodeMap *_node_index;
     2063    IntArcMap *_node_heap_index;
    20662064
    20672065    struct NodeData {
    20682066
    2069       NodeData(ArcIntMap& node_heap_index)
     2067      NodeData(IntArcMap& node_heap_index)
    20702068        : heap(node_heap_index) {}
    20712069
    20722070      int blossom;
    20732071      Value pot;
    2074       BinHeap<Value, ArcIntMap> heap;
     2072      BinHeap<Value, IntArcMap> heap;
    20752073      std::map<int, Arc> heap_index;
    20762074
     
    20882086    BinHeap<Value, IntIntMap> *_delta2;
    20892087
    2090     EdgeIntMap *_delta3_index;
    2091     BinHeap<Value, EdgeIntMap> *_delta3;
     2088    IntEdgeMap *_delta3_index;
     2089    BinHeap<Value, IntEdgeMap> *_delta3;
    20922090
    20932091    IntIntMap *_delta4_index;
     
    21072105      }
    21082106      if (!_blossom_set) {
    2109         _blossom_index = new NodeIntMap(_graph);
     2107        _blossom_index = new IntNodeMap(_graph);
    21102108        _blossom_set = new BlossomSet(*_blossom_index);
    21112109        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
     
    21132111
    21142112      if (!_node_index) {
    2115         _node_index = new NodeIntMap(_graph);
    2116         _node_heap_index = new ArcIntMap(_graph);
     2113        _node_index = new IntNodeMap(_graph);
     2114        _node_heap_index = new IntArcMap(_graph);
    21172115        _node_data = new RangeMap<NodeData>(_node_num,
    2118                                               NodeData(*_node_heap_index));
     2116                                            NodeData(*_node_heap_index));
    21192117      }
    21202118
     
    21282126      }
    21292127      if (!_delta3) {
    2130         _delta3_index = new EdgeIntMap(_graph);
    2131         _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
     2128        _delta3_index = new IntEdgeMap(_graph);
     2129        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
    21322130      }
    21332131      if (!_delta4) {
     
    24622460    }
    24632461
    2464     void augmentOnArc(const Edge& arc) {
    2465 
    2466       int left = _blossom_set->find(_graph.u(arc));
    2467       int right = _blossom_set->find(_graph.v(arc));
     2462    void augmentOnEdge(const Edge& edge) {
     2463
     2464      int left = _blossom_set->find(_graph.u(edge));
     2465      int right = _blossom_set->find(_graph.v(edge));
    24682466
    24692467      int left_tree = _tree_set->find(left);
     
    24752473      destroyTree(right_tree);
    24762474
    2477       (*_blossom_data)[left].next = _graph.direct(arc, true);
    2478       (*_blossom_data)[right].next = _graph.direct(arc, false);
     2475      (*_blossom_data)[left].next = _graph.direct(edge, true);
     2476      (*_blossom_data)[right].next = _graph.direct(edge, false);
    24792477    }
    24802478
     
    24962494    }
    24972495
    2498     void shrinkOnArc(const Edge& edge, int tree) {
     2496    void shrinkOnEdge(const Edge& edge, int tree) {
    24992497      int nca = -1;
    25002498      std::vector<int> left_path, right_path;
     
    28322830
    28332831      for (ArcIt e(_graph); e != INVALID; ++e) {
    2834         _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
     2832        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
    28352833      }
    28362834      for (EdgeIt e(_graph); e != INVALID; ++e) {
     
    29252923
    29262924              if (left_tree == right_tree) {
    2927                 shrinkOnArc(e, left_tree);
     2925                shrinkOnEdge(e, left_tree);
    29282926              } else {
    2929                 augmentOnArc(e);
     2927                augmentOnEdge(e);
    29302928                unmatched -= 2;
    29312929              }
     
    29752973    }
    29762974
    2977     /// \brief Returns true when the arc is in the matching.
    2978     ///
    2979     /// Returns true when the arc is in the matching.
    2980     bool matching(const Edge& arc) const {
    2981       return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
    2982     }
    2983 
    2984     /// \brief Returns the incident matching arc.
    2985     ///
    2986     /// Returns the incident matching arc from given node.
     2975    /// \brief Returns true when the edge is in the matching.
     2976    ///
     2977    /// Returns true when the edge is in the matching.
     2978    bool matching(const Edge& edge) const {
     2979      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
     2980    }
     2981
     2982    /// \brief Returns the incident matching edge.
     2983    ///
     2984    /// Returns the incident matching arc from given edge.
    29872985    Arc matching(const Node& node) const {
    29882986      return (*_matching)[node];
     
    30673065      }
    30683066
    3069       /// \brief Invalid constructor.
    3070       ///
    3071       /// Invalid constructor.
    3072       BlossomIt(Invalid) : _index(-1) {}
    3073 
    30743067      /// \brief Conversion to node.
    30753068      ///
    30763069      /// Conversion to node.
    30773070      operator Node() const {
    3078         return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
     3071        return _algorithm->_blossom_node_list[_index];
    30793072      }
    30803073
     
    30843077      BlossomIt& operator++() {
    30853078        ++_index;
    3086         if (_index == _last) {
    3087           _index = -1;
    3088         }
    30893079        return *this;
    30903080      }
    30913081
    3092       bool operator==(const BlossomIt& it) const {
    3093         return _index == it._index;
    3094       }
    3095       bool operator!=(const BlossomIt& it) const {
    3096         return _index != it._index;
    3097       }
     3082      /// \brief Validity checking
     3083      ///
     3084      /// Checks whether the iterator is invalid.
     3085      bool operator==(Invalid) const { return _index == _last; }
     3086
     3087      /// \brief Validity checking
     3088      ///
     3089      /// Checks whether the iterator is valid.
     3090      bool operator!=(Invalid) const { return _index != _last; }
    30983091
    30993092    private:
  • test/CMakeLists.txt

    r326 r327  
    1818  maps_test
    1919  max_matching_test
    20   max_weighted_matching_test
    2120  random_test
    2221  path_test
  • test/Makefile.am

    r326 r327  
    2121        test/maps_test \
    2222        test/max_matching_test \
    23         test/max_weighted_matching_test \
    2423        test/random_test \
    2524        test/path_test \
     
    4645test_maps_test_SOURCES = test/maps_test.cc
    4746test_max_matching_test_SOURCES = test/max_matching_test.cc
    48 test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
    4947test_path_test_SOURCES = test/path_test.cc
    5048test_random_test_SOURCES = test/random_test.cc
  • test/max_matching_test.cc

    r326 r327  
    1818
    1919#include <iostream>
     20#include <sstream>
    2021#include <vector>
    2122#include <queue>
     
    2324#include <cstdlib>
    2425
     26#include <lemon/max_matching.h>
     27#include <lemon/smart_graph.h>
     28#include <lemon/lgf_reader.h>
     29
    2530#include "test_tools.h"
    26 #include <lemon/list_graph.h>
    27 #include <lemon/max_matching.h>
    2831
    2932using namespace std;
    3033using namespace lemon;
    3134
     35GRAPH_TYPEDEFS(SmartGraph);
     36
     37
     38const int lgfn = 3;
     39const std::string lgf[lgfn] = {
     40  "@nodes\n"
     41  "label\n"
     42  "0\n"
     43  "1\n"
     44  "2\n"
     45  "3\n"
     46  "4\n"
     47  "5\n"
     48  "6\n"
     49  "7\n"
     50  "@edges\n"
     51  "     label  weight\n"
     52  "7 4  0      984\n"
     53  "0 7  1      73\n"
     54  "7 1  2      204\n"
     55  "2 3  3      583\n"
     56  "2 7  4      565\n"
     57  "2 1  5      582\n"
     58  "0 4  6      551\n"
     59  "2 5  7      385\n"
     60  "1 5  8      561\n"
     61  "5 3  9      484\n"
     62  "7 5  10     904\n"
     63  "3 6  11     47\n"
     64  "7 6  12     888\n"
     65  "3 0  13     747\n"
     66  "6 1  14     310\n",
     67
     68  "@nodes\n"
     69  "label\n"
     70  "0\n"
     71  "1\n"
     72  "2\n"
     73  "3\n"
     74  "4\n"
     75  "5\n"
     76  "6\n"
     77  "7\n"
     78  "@edges\n"
     79  "     label  weight\n"
     80  "2 5  0      710\n"
     81  "0 5  1      241\n"
     82  "2 4  2      856\n"
     83  "2 6  3      762\n"
     84  "4 1  4      747\n"
     85  "6 1  5      962\n"
     86  "4 7  6      723\n"
     87  "1 7  7      661\n"
     88  "2 3  8      376\n"
     89  "1 0  9      416\n"
     90  "6 7  10     391\n",
     91
     92  "@nodes\n"
     93  "label\n"
     94  "0\n"
     95  "1\n"
     96  "2\n"
     97  "3\n"
     98  "4\n"
     99  "5\n"
     100  "6\n"
     101  "7\n"
     102  "@edges\n"
     103  "     label  weight\n"
     104  "6 2  0      553\n"
     105  "0 7  1      653\n"
     106  "6 3  2      22\n"
     107  "4 7  3      846\n"
     108  "7 2  4      981\n"
     109  "7 6  5      250\n"
     110  "5 2  6      539\n",
     111};
     112
     113void checkMatching(const SmartGraph& graph,
     114                   const MaxMatching<SmartGraph>& mm) {
     115  int num = 0;
     116
     117  IntNodeMap comp_index(graph);
     118  UnionFind<IntNodeMap> comp(comp_index);
     119
     120  int barrier_num = 0;
     121
     122  for (NodeIt n(graph); n != INVALID; ++n) {
     123    check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
     124          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
     125    if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
     126      ++barrier_num;
     127    } else {
     128      comp.insert(n);
     129    }
     130  }
     131
     132  for (EdgeIt e(graph); e != INVALID; ++e) {
     133    if (mm.matching(e)) {
     134      check(e == mm.matching(graph.u(e)), "Wrong matching");
     135      check(e == mm.matching(graph.v(e)), "Wrong matching");
     136      ++num;
     137    }
     138    check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
     139          mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
     140          "Wrong Gallai-Edmonds decomposition");
     141
     142    check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
     143          mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
     144          "Wrong Gallai-Edmonds decomposition");
     145
     146    if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
     147        mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
     148      comp.join(graph.u(e), graph.v(e));
     149    }
     150  }
     151
     152  std::set<int> comp_root;
     153  int odd_comp_num = 0;
     154  for (NodeIt n(graph); n != INVALID; ++n) {
     155    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
     156      int root = comp.find(n);
     157      if (comp_root.find(root) == comp_root.end()) {
     158        comp_root.insert(root);
     159        if (comp.size(n) % 2 == 1) {
     160          ++odd_comp_num;
     161        }
     162      }
     163    }
     164  }
     165
     166  check(mm.matchingSize() == num, "Wrong matching");
     167  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
     168         "Wrong matching");
     169  return;
     170}
     171
     172void checkWeightedMatching(const SmartGraph& graph,
     173                   const SmartGraph::EdgeMap<int>& weight,
     174                   const MaxWeightedMatching<SmartGraph>& mwm) {
     175  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     176    if (graph.u(e) == graph.v(e)) continue;
     177    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
     178
     179    for (int i = 0; i < mwm.blossomNum(); ++i) {
     180      bool s = false, t = false;
     181      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
     182           n != INVALID; ++n) {
     183        if (graph.u(e) == n) s = true;
     184        if (graph.v(e) == n) t = true;
     185      }
     186      if (s == true && t == true) {
     187        rw += mwm.blossomValue(i);
     188      }
     189    }
     190    rw -= weight[e] * mwm.dualScale;
     191
     192    check(rw >= 0, "Negative reduced weight");
     193    check(rw == 0 || !mwm.matching(e),
     194          "Non-zero reduced weight on matching edge");
     195  }
     196
     197  int pv = 0;
     198  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     199    if (mwm.matching(n) != INVALID) {
     200      check(mwm.nodeValue(n) >= 0, "Invalid node value");
     201      pv += weight[mwm.matching(n)];
     202      SmartGraph::Node o = graph.target(mwm.matching(n));
     203      check(mwm.mate(n) == o, "Invalid matching");
     204      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
     205            "Invalid matching");
     206    } else {
     207      check(mwm.mate(n) == INVALID, "Invalid matching");
     208      check(mwm.nodeValue(n) == 0, "Invalid matching");
     209    }
     210  }
     211
     212  int dv = 0;
     213  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     214    dv += mwm.nodeValue(n);
     215  }
     216
     217  for (int i = 0; i < mwm.blossomNum(); ++i) {
     218    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
     219    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
     220    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
     221  }
     222
     223  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
     224
     225  return;
     226}
     227
     228void checkWeightedPerfectMatching(const SmartGraph& graph,
     229                          const SmartGraph::EdgeMap<int>& weight,
     230                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
     231  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     232    if (graph.u(e) == graph.v(e)) continue;
     233    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
     234
     235    for (int i = 0; i < mwpm.blossomNum(); ++i) {
     236      bool s = false, t = false;
     237      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
     238           n != INVALID; ++n) {
     239        if (graph.u(e) == n) s = true;
     240        if (graph.v(e) == n) t = true;
     241      }
     242      if (s == true && t == true) {
     243        rw += mwpm.blossomValue(i);
     244      }
     245    }
     246    rw -= weight[e] * mwpm.dualScale;
     247
     248    check(rw >= 0, "Negative reduced weight");
     249    check(rw == 0 || !mwpm.matching(e),
     250          "Non-zero reduced weight on matching edge");
     251  }
     252
     253  int pv = 0;
     254  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     255    check(mwpm.matching(n) != INVALID, "Non perfect");
     256    pv += weight[mwpm.matching(n)];
     257    SmartGraph::Node o = graph.target(mwpm.matching(n));
     258    check(mwpm.mate(n) == o, "Invalid matching");
     259    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
     260          "Invalid matching");
     261  }
     262
     263  int dv = 0;
     264  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     265    dv += mwpm.nodeValue(n);
     266  }
     267
     268  for (int i = 0; i < mwpm.blossomNum(); ++i) {
     269    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
     270    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
     271    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
     272  }
     273
     274  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
     275
     276  return;
     277}
     278
     279
    32280int main() {
    33281
    34   typedef ListGraph Graph;
    35 
    36   GRAPH_TYPEDEFS(Graph);
    37 
    38   Graph g;
    39   g.clear();
    40 
    41   std::vector<Graph::Node> nodes;
    42   for (int i=0; i<13; ++i)
    43       nodes.push_back(g.addNode());
    44 
    45   g.addEdge(nodes[0], nodes[0]);
    46   g.addEdge(nodes[6], nodes[10]);
    47   g.addEdge(nodes[5], nodes[10]);
    48   g.addEdge(nodes[4], nodes[10]);
    49   g.addEdge(nodes[3], nodes[11]);
    50   g.addEdge(nodes[1], nodes[6]);
    51   g.addEdge(nodes[4], nodes[7]);
    52   g.addEdge(nodes[1], nodes[8]);
    53   g.addEdge(nodes[0], nodes[8]);
    54   g.addEdge(nodes[3], nodes[12]);
    55   g.addEdge(nodes[6], nodes[9]);
    56   g.addEdge(nodes[9], nodes[11]);
    57   g.addEdge(nodes[2], nodes[10]);
    58   g.addEdge(nodes[10], nodes[8]);
    59   g.addEdge(nodes[5], nodes[8]);
    60   g.addEdge(nodes[6], nodes[3]);
    61   g.addEdge(nodes[0], nodes[5]);
    62   g.addEdge(nodes[6], nodes[12]);
    63 
    64   MaxMatching<Graph> max_matching(g);
    65   max_matching.init();
    66   max_matching.startDense();
    67 
    68   int s=0;
    69   Graph::NodeMap<Node> mate(g,INVALID);
    70   max_matching.mateMap(mate);
    71   for(NodeIt v(g); v!=INVALID; ++v) {
    72     if ( mate[v]!=INVALID ) ++s;
    73   }
    74   int size=int(s/2);  //size will be used as the size of a maxmatching
    75 
    76   for(NodeIt v(g); v!=INVALID; ++v) {
    77     max_matching.mate(v);
    78   }
    79 
    80   check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
    81 
    82   Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
    83   max_matching.decomposition(pos0);
    84 
    85   max_matching.init();
    86   max_matching.startSparse();
    87   s=0;
    88   max_matching.mateMap(mate);
    89   for(NodeIt v(g); v!=INVALID; ++v) {
    90     if ( mate[v]!=INVALID ) ++s;
    91   }
    92   check ( int(s/2) == size, "The size does not equal!" );
    93 
    94   Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
    95   max_matching.decomposition(pos1);
    96 
    97   max_matching.run();
    98   s=0;
    99   max_matching.mateMap(mate);
    100   for(NodeIt v(g); v!=INVALID; ++v) {
    101     if ( mate[v]!=INVALID ) ++s;
    102   }
    103   check ( int(s/2) == size, "The size does not equal!" );
    104 
    105   Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
    106   max_matching.decomposition(pos2);
    107 
    108   max_matching.run();
    109   s=0;
    110   max_matching.mateMap(mate);
    111   for(NodeIt v(g); v!=INVALID; ++v) {
    112     if ( mate[v]!=INVALID ) ++s;
    113   }
    114   check ( int(s/2) == size, "The size does not equal!" );
    115 
    116   Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
    117   max_matching.decomposition(pos);
    118 
    119   bool ismatching=true;
    120   for(NodeIt v(g); v!=INVALID; ++v) {
    121     if ( mate[v]!=INVALID ) {
    122       Node u=mate[v];
    123       if (mate[u]!=v) ismatching=false;
    124     }
    125   }
    126   check ( ismatching, "It is not a matching!" );
    127 
    128   bool coincide=true;
    129   for(NodeIt v(g); v!=INVALID; ++v) {
    130    if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
    131      coincide=false;
    132     }
    133   }
    134   check ( coincide, "The decompositions do not coincide! " );
    135 
    136   bool noarc=true;
    137   for(EdgeIt e(g); e!=INVALID; ++e) {
    138    if ( (pos[g.v(e)]==max_matching.C &&
    139          pos[g.u(e)]==max_matching.D) ||
    140          (pos[g.v(e)]==max_matching.D &&
    141           pos[g.u(e)]==max_matching.C) )
    142       noarc=false;
    143   }
    144   check ( noarc, "There are arcs between D and C!" );
    145 
    146   bool oddcomp=true;
    147   Graph::NodeMap<bool> todo(g,true);
    148   int num_comp=0;
    149   for(NodeIt v(g); v!=INVALID; ++v) {
    150    if ( pos[v]==max_matching.D && todo[v] ) {
    151       int comp_size=1;
    152       ++num_comp;
    153       std::queue<Node> Q;
    154       Q.push(v);
    155       todo.set(v,false);
    156       while (!Q.empty()) {
    157         Node w=Q.front();
    158         Q.pop();
    159         for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
    160           Node u=g.runningNode(e);
    161           if ( pos[u]==max_matching.D && todo[u] ) {
    162             ++comp_size;
    163             Q.push(u);
    164             todo.set(u,false);
    165           }
    166         }
    167       }
    168       if ( !(comp_size % 2) ) oddcomp=false;
    169     }
    170   }
    171   check ( oddcomp, "A component of g[D] is not odd." );
    172 
    173   int barrier=0;
    174   for(NodeIt v(g); v!=INVALID; ++v) {
    175     if ( pos[v]==max_matching.A ) ++barrier;
    176   }
    177   int expected_size=int( countNodes(g)-num_comp+barrier)/2;
    178   check ( size==expected_size, "The size of the matching is wrong." );
     282  for (int i = 0; i < lgfn; ++i) {
     283    SmartGraph graph;
     284    SmartGraph::EdgeMap<int> weight(graph);
     285
     286    istringstream lgfs(lgf[i]);
     287    graphReader(graph, lgfs).
     288      edgeMap("weight", weight).run();
     289
     290    MaxMatching<SmartGraph> mm(graph);
     291    mm.run();
     292    checkMatching(graph, mm);
     293
     294    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
     295    mwm.run();
     296    checkWeightedMatching(graph, weight, mwm);
     297
     298    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
     299    bool perfect = mwpm.run();
     300
     301    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
     302          "Perfect matching found");
     303
     304    if (perfect) {
     305      checkWeightedPerfectMatching(graph, weight, mwpm);
     306    }
     307  }
    179308
    180309  return 0;
Note: See TracChangeset for help on using the changeset viewer.