COIN-OR::LEMON - Graph Library

Ticket #168: 3180fd18651b.patch

File 3180fd18651b.patch, 105.7 KB (added by Balazs Dezso, 2 years ago)
  • doc/groups.dox

    # HG changeset patch
    # User Balazs Dezso <deba@google.com>
    # Date 1543258354 -3600
    #      Mon Nov 26 19:52:34 2018 +0100
    # Node ID 3180fd18651baf7c21facac7b7952aa1704d9afd
    # Parent  4fd76139b69ec938db4d67ce0927ad5923647ed2
    MaxWeightBpMatching algorithm variants
    
    diff -r 4fd76139b69e -r 3180fd18651b doc/groups.dox
    a b  
    523523  for calculating maximum cardinality matching in bipartite graphs.
    524524- \ref PrBipartiteMatching Push-relabel algorithm
    525525  for calculating maximum cardinality matching in bipartite graphs.
    526 - \ref MaxWeightedBipartiteMatching
    527   Successive shortest path algorithm for calculating maximum weighted
    528   matching and maximum weighted bipartite matching in bipartite graphs.
     526- \ref MaxWeightedBpMatching Multiple search tree augmenting path algorithm for
     527  calculating maximum weighted matching in bipartite graphs.
    529528- \ref MinCostMaxBipartiteMatching
    530529  Successive shortest path algorithm for calculating minimum cost maximum
    531530  matching in bipartite graphs.
  • new file lemon/bpmatching.h

    diff -r 4fd76139b69e -r 3180fd18651b lemon/bpmatching.h
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2018
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_BPMATCHING_H
     20#define LEMON_BPMATCHING_H
     21
     22#include <limits>
     23
     24#include <lemon/core.h>
     25#include <lemon/unionfind.h>
     26#include <lemon/bin_heap.h>
     27#include <lemon/maps.h>
     28
     29///\ingroup matching
     30///\file
     31///\brief Maximum matching algorithms in bipartite graphs.
     32
     33namespace lemon {
     34
     35  /// \ingroup matching
     36  ///
     37  /// \brief Weighted matching in bipartite graphs
     38  ///
     39  /// This class provides an efficient implementation of multiple search tree
     40  /// augmenting path matching algorithm. The implementation is based on
     41  /// extensive use of priority queues and provides \f$O(nm\log n)\f$ time
     42  /// complexity.
     43  ///
     44  /// The maximum weighted matching problem is to find a subset of the
     45  /// edges in a bipartite graph with maximum overall weight for which
     46  /// each node has at most one incident edge.
     47  /// It can be formulated with the following linear program.
     48  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
     49  /// \f[x_e \ge 0\quad \forall e\in E\f]
     50  /// \f[\max \sum_{e\in E}x_ew_e\f]
     51  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
     52  /// \f$X\f$.
     53  ///
     54  /// The algorithm calculates an optimal matching and a proof of the
     55  /// optimality. The solution of the dual problem can be used to check
     56  /// the result of the algorithm. The dual linear problem is the
     57  /// following.
     58  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
     59  /// \f[y_u \ge 0 \quad \forall u \in V\f]
     60  /// \f[\min \sum_{u \in V}y_u \f]
     61  /// \tparam BPGR The bipartite graph type the algorithm runs on.
     62  /// \tparam WM The type edge weight map. The default type is
     63  /// \ref concepts::BpGraph::EdgeMap "BPGR::EdgeMap<int>".
     64#ifdef DOXYGEN
     65  template <typename BPGR, typename WM>
     66#else
     67  template <typename BPGR,
     68            typename WM = typename BPGR::template EdgeMap<int> >
     69#endif
     70  class MaxWeightedBpMatching1 {
     71  public:
     72
     73    /// The graph type of the algorithm
     74    typedef BPGR BpGraph;
     75    /// The type of the edge weight map
     76    typedef WM WeightMap;
     77    /// The value type of the edge weights
     78    typedef typename WeightMap::Value Value;
     79
     80    /// The type of the matching map
     81    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     82    MatchingMap;
     83
     84    /// \brief Scaling factor for dual solution
     85    ///
     86    /// Scaling factor for dual solution. It is equal to 4 or 1
     87    /// according to the value type.
     88    static const int dualScale =
     89      std::numeric_limits<Value>::is_integer ? 4 : 1;
     90
     91  private:
     92
     93    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     94
     95    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     96
     97    const BpGraph& _bpgraph;
     98    const WeightMap& _weight;
     99
     100    MatchingMap* _matching;
     101
     102    NodePotential* _node_potential;
     103
     104    int _node_num;
     105
     106    enum Status {
     107      EVEN = -1, MATCHED = 0, ODD = 1
     108    };
     109
     110    typedef typename BpGraph::template NodeMap<Status> StatusMap;
     111    StatusMap* _status;
     112
     113    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     114    PredMap* _pred;
     115
     116    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     117    IntNodeMap *_tree_set_index;
     118    TreeSet *_tree_set;
     119
     120    IntNodeMap *_delta1_index;
     121    BinHeap<Value, IntNodeMap> *_delta1;
     122
     123    IntNodeMap *_delta2_index;
     124    BinHeap<Value, IntNodeMap> *_delta2;
     125
     126    IntEdgeMap *_delta3_index;
     127    BinHeap<Value, IntEdgeMap> *_delta3;
     128
     129    Value _delta_sum;
     130    int _unmatched;
     131
     132    void createStructures() {
     133      _node_num = countNodes(_bpgraph);
     134
     135      if (!_matching) {
     136        _matching = new MatchingMap(_bpgraph);
     137      }
     138
     139      if (!_node_potential) {
     140        _node_potential = new NodePotential(_bpgraph);
     141      }
     142
     143      if (!_status) {
     144        _status = new StatusMap(_bpgraph);
     145      }
     146
     147      if (!_pred) {
     148        _pred = new PredMap(_bpgraph);
     149      }
     150
     151      if (!_tree_set) {
     152        _tree_set_index = new IntNodeMap(_bpgraph);
     153        _tree_set = new TreeSet(*_tree_set_index);
     154      }
     155
     156      if (!_delta1) {
     157        _delta1_index = new IntNodeMap(_bpgraph);
     158        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
     159      }
     160
     161      if (!_delta2) {
     162        _delta2_index = new IntNodeMap(_bpgraph);
     163        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
     164      }
     165
     166      if (!_delta3) {
     167        _delta3_index = new IntEdgeMap(_bpgraph);
     168        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
     169      }
     170    }
     171
     172    void destroyStructures() {
     173      if (_matching) {
     174        delete _matching;
     175      }
     176      if (_node_potential) {
     177        delete _node_potential;
     178      }
     179      if (_status) {
     180        delete _status;
     181      }
     182      if (_pred) {
     183        delete _pred;
     184      }
     185      if (_tree_set) {
     186        delete _tree_set_index;
     187        delete _tree_set;
     188      }
     189      if (_delta2) {
     190        delete _delta2_index;
     191        delete _delta2;
     192      }
     193      if (_delta3) {
     194        delete _delta3_index;
     195        delete _delta3;
     196      }
     197    }
     198
     199    void matchedToEven(Node node, int tree) {
     200      _tree_set->insert(node, tree);
     201      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     202      _delta1->push(node, (*_node_potential)[node]);
     203
     204      if (_delta2->state(node) == _delta2->IN_HEAP) {
     205        _delta2->erase(node);
     206      }
     207
     208      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     209        Node v = _bpgraph.source(a);
     210        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     211          dualScale * _weight[a];
     212        if ((*_status)[v] == EVEN) {
     213          _delta3->push(a, rw / 2);
     214        } else if ((*_status)[v] == MATCHED) {
     215          if (_delta2->state(v) != _delta2->IN_HEAP) {
     216            _pred->set(v, a);
     217            _delta2->push(v, rw);
     218          } else if ((*_delta2)[v] > rw) {
     219            _pred->set(v, a);
     220            _delta2->decrease(v, rw);
     221          }
     222        }
     223      }
     224    }
     225
     226    void matchedToOdd(Node node, int tree) {
     227      _tree_set->insert(node, tree);
     228      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     229
     230      if (_delta2->state(node) == _delta2->IN_HEAP) {
     231        _delta2->erase(node);
     232      }
     233    }
     234
     235    void evenToMatched(Node node, int tree) {
     236      _delta1->erase(node);
     237      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     238      Arc min = INVALID;
     239      Value minrw = std::numeric_limits<Value>::max();
     240      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     241        Node v = _bpgraph.source(a);
     242        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     243          dualScale * _weight[a];
     244
     245        if ((*_status)[v] == EVEN) {
     246          _delta3->erase(a);
     247          if (minrw > rw) {
     248            min = _bpgraph.oppositeArc(a);
     249            minrw = rw;
     250          }
     251        } else if ((*_status)[v]  == MATCHED) {
     252          if ((*_pred)[v] == a) {
     253            Arc mina = INVALID;
     254            Value minrwa = std::numeric_limits<Value>::max();
     255            for (OutArcIt aa(_bpgraph, v); aa != INVALID; ++aa) {
     256              Node va = _bpgraph.target(aa);
     257              if ((*_status)[va] != EVEN ||
     258                  _tree_set->find(va) == tree) continue;
     259              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
     260                dualScale * _weight[aa];
     261              if (minrwa > rwa) {
     262                minrwa = rwa;
     263                mina = aa;
     264              }
     265            }
     266            if (mina != INVALID) {
     267              _pred->set(v, mina);
     268              _delta2->increase(v, minrwa);
     269            } else {
     270              _pred->set(v, INVALID);
     271              _delta2->erase(v);
     272            }
     273          }
     274        }
     275      }
     276      if (min != INVALID) {
     277        _pred->set(node, min);
     278        _delta2->push(node, minrw);
     279      } else {
     280        _pred->set(node, INVALID);
     281      }
     282    }
     283
     284    void oddToMatched(Node node) {
     285      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     286      Arc min = INVALID;
     287      Value minrw = std::numeric_limits<Value>::max();
     288      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     289        Node v = _bpgraph.source(a);
     290        if ((*_status)[v] != EVEN) continue;
     291        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     292          dualScale * _weight[a];
     293
     294        if (minrw > rw) {
     295          min = _bpgraph.oppositeArc(a);
     296          minrw = rw;
     297        }
     298      }
     299      if (min != INVALID) {
     300        _pred->set(node, min);
     301        _delta2->push(node, minrw);
     302      } else {
     303        _pred->set(node, INVALID);
     304      }
     305    }
     306
     307    void alternatePath(Node even, int tree) {
     308      Node odd;
     309
     310      _status->set(even, MATCHED);
     311      evenToMatched(even, tree);
     312
     313      Arc prev = (*_matching)[even];
     314      while (prev != INVALID) {
     315        odd = _bpgraph.target(prev);
     316        even = _bpgraph.target((*_pred)[odd]);
     317        _matching->set(odd, (*_pred)[odd]);
     318        _status->set(odd, MATCHED);
     319        oddToMatched(odd);
     320
     321        prev = (*_matching)[even];
     322        _status->set(even, MATCHED);
     323        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     324        evenToMatched(even, tree);
     325      }
     326    }
     327
     328    void destroyTree(int tree) {
     329      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     330        if ((*_status)[n] == EVEN) {
     331          _status->set(n, MATCHED);
     332          evenToMatched(n, tree);
     333        } else if ((*_status)[n] == ODD) {
     334          _status->set(n, MATCHED);
     335          oddToMatched(n);
     336        }
     337      }
     338      _tree_set->eraseClass(tree);
     339    }
     340
     341    void unmatchNode(const Node& node) {
     342      int tree = _tree_set->find(node);
     343
     344      alternatePath(node, tree);
     345      destroyTree(tree);
     346
     347      _matching->set(node, INVALID);
     348    }
     349
     350    void augmentOnEdge(const Edge& edge) {
     351      Node left = _bpgraph.u(edge);
     352      int left_tree = _tree_set->find(left);
     353
     354      alternatePath(left, left_tree);
     355      destroyTree(left_tree);
     356      _matching->set(left, _bpgraph.direct(edge, true));
     357
     358      Node right = _bpgraph.v(edge);
     359      int right_tree = _tree_set->find(right);
     360
     361      alternatePath(right, right_tree);
     362      destroyTree(right_tree);
     363      _matching->set(right, _bpgraph.direct(edge, false));
     364    }
     365
     366    void augmentOnArc(const Arc& arc) {
     367      Node left = _bpgraph.source(arc);
     368      _status->set(left, MATCHED);
     369      _matching->set(left, arc);
     370      _pred->set(left, arc);
     371
     372      Node right = _bpgraph.target(arc);
     373      int right_tree = _tree_set->find(right);
     374
     375      alternatePath(right, right_tree);
     376      destroyTree(right_tree);
     377      _matching->set(right, _bpgraph.oppositeArc(arc));
     378    }
     379
     380    void extendOnArc(const Arc& arc) {
     381      Node base = _bpgraph.target(arc);
     382      int tree = _tree_set->find(base);
     383
     384      Node odd = _bpgraph.source(arc);
     385      _tree_set->insert(odd, tree);
     386      _status->set(odd, ODD);
     387      matchedToOdd(odd, tree);
     388      _pred->set(odd, arc);
     389
     390      Node even = _bpgraph.target((*_matching)[odd]);
     391      _tree_set->insert(even, tree);
     392      _status->set(even, EVEN);
     393      matchedToEven(even, tree);
     394    }
     395
     396  public:
     397
     398    /// \brief Constructor
     399    ///
     400    /// Constructor.
     401    MaxWeightedBpMatching1(const BpGraph& bpgraph, const WeightMap& weight)
     402      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     403        _node_potential(0), _node_num(0),
     404        _status(0), _pred(0),
     405        _tree_set_index(0), _tree_set(0),
     406
     407        _delta1_index(0), _delta1(0),
     408        _delta2_index(0), _delta2(0),
     409        _delta3_index(0), _delta3(0),
     410
     411        _delta_sum(), _unmatched(0)
     412    {}
     413
     414    ~MaxWeightedBpMatching1() {
     415      destroyStructures();
     416    }
     417
     418    /// \name Execution Control
     419    /// The simplest way to execute the algorithm is to use the
     420    /// \ref run() member function.
     421
     422    ///@{
     423
     424    /// \brief Initialize the algorithm
     425    ///
     426    /// This function initializes the algorithm.
     427    void init() {
     428      createStructures();
     429
     430      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     431        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     432        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     433      }
     434      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     435        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     436      }
     437
     438      _delta1->clear();
     439      _delta2->clear();
     440      _delta3->clear();
     441      _tree_set->clear();
     442
     443      _unmatched = _node_num;
     444
     445      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     446        Value max = 0;
     447        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     448          if ((dualScale * _weight[a]) / 2 > max) {
     449            max = (dualScale * _weight[a]) / 2;
     450          }
     451        }
     452        _node_potential->set(n, max);
     453        _delta1->push(n, max);
     454
     455        _tree_set->insert(n);
     456
     457        _matching->set(n, INVALID);
     458        _status->set(n, EVEN);
     459      }
     460      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     461        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
     462                          (*_node_potential)[_bpgraph.v(e)] -
     463                          dualScale * _weight[e]) / 2);
     464      }
     465    }
     466
     467    /// \brief Initialize the algorithm
     468    ///
     469    /// This function initializes the algorithm.
     470    void redRootInit() {
     471      createStructures();
     472
     473      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     474        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     475        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     476      }
     477      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     478        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     479      }
     480
     481      _delta1->clear();
     482      _delta2->clear();
     483      _delta3->clear();
     484      _tree_set->clear();
     485
     486      _unmatched = 0;
     487      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     488        Value max = 0;
     489        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     490          if ((dualScale * _weight[a]) > max) {
     491            max = dualScale * _weight[a];
     492          }
     493        }
     494        _node_potential->set(n, max);
     495        _delta1->push(n, max);
     496
     497        _tree_set->insert(n);
     498
     499        _matching->set(n, INVALID);
     500        _status->set(n, EVEN);
     501
     502        ++_unmatched;
     503      }
     504      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     505        _matching->set(n, INVALID);
     506        _status->set(n, MATCHED);
     507
     508        Arc min = INVALID;
     509        Value minrw = std::numeric_limits<Value>::max();
     510        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
     511          Node v = _bpgraph.source(a);
     512          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     513                     dualScale * _weight[a];
     514
     515          if (minrw > rw) {
     516            min = _bpgraph.oppositeArc(a);
     517            minrw = rw;
     518          }
     519        }
     520        if (min != INVALID) {
     521          _pred->set(n, min);
     522          _delta2->push(n, minrw);
     523        } else {
     524          _pred->set(n, INVALID);
     525        }
     526      }
     527    }
     528
     529    /// \brief Initialize the algorithm
     530    ///
     531    /// This function initializes the algorithm.
     532    void blueRootInit() {
     533      createStructures();
     534
     535      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     536        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     537        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     538      }
     539      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     540        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     541      }
     542
     543      _delta1->clear();
     544      _delta2->clear();
     545      _delta3->clear();
     546      _tree_set->clear();
     547
     548      _unmatched = 0;
     549      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     550        Value max = 0;
     551        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     552          if ((dualScale * _weight[a]) > max) {
     553            max = dualScale * _weight[a];
     554          }
     555        }
     556        _node_potential->set(n, max);
     557        _delta1->push(n, max);
     558
     559        _tree_set->insert(n);
     560
     561        _matching->set(n, INVALID);
     562        _status->set(n, EVEN);
     563
     564        ++_unmatched;
     565      }
     566      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     567        _matching->set(n, INVALID);
     568        _status->set(n, MATCHED);
     569
     570        Arc min = INVALID;
     571        Value minrw = std::numeric_limits<Value>::max();
     572        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
     573          Node v = _bpgraph.source(a);
     574          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     575                     dualScale * _weight[a];
     576
     577          if (minrw > rw) {
     578            min = _bpgraph.oppositeArc(a);
     579            minrw = rw;
     580          }
     581        }
     582        if (min != INVALID) {
     583          _pred->set(n, min);
     584          _delta2->push(n, minrw);
     585        } else {
     586          _pred->set(n, INVALID);
     587        }
     588      }
     589    }
     590
     591    /// \brief Start the algorithm
     592    ///
     593    /// This function starts the algorithm.
     594    ///
     595    /// \pre \ref init() must be called before using this function.
     596    void start() {
     597      enum OpType {
     598        D1, D2, D3
     599      };
     600
     601      while (_unmatched > 0) {
     602        Value d1 = !_delta1->empty() ?
     603          _delta1->prio() : std::numeric_limits<Value>::max();
     604
     605        Value d2 = !_delta2->empty() ?
     606          _delta2->prio() : std::numeric_limits<Value>::max();
     607
     608        Value d3 = !_delta3->empty() ?
     609          _delta3->prio() : std::numeric_limits<Value>::max();
     610
     611        _delta_sum = d3; OpType ot = D3;
     612        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
     613        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     614
     615        switch (ot) {
     616        case D1:
     617          {
     618            Node n = _delta1->top();
     619            unmatchNode(n);
     620            --_unmatched;
     621          }
     622          break;
     623        case D2:
     624          {
     625            Node n = _delta2->top();
     626            Arc a = (*_pred)[n];
     627            if ((*_matching)[n] == INVALID) {
     628              augmentOnArc(a);
     629              --_unmatched;
     630            } else {
     631              extendOnArc(a);
     632            }
     633          }
     634          break;
     635        case D3:
     636          {
     637            Edge e = _delta3->top();
     638            augmentOnEdge(e);
     639            _unmatched -= 2;
     640          }
     641          break;
     642        }
     643      }
     644    }
     645
     646    /// \brief Run the algorithm.
     647    ///
     648    /// This method runs the \c %MaxWeightedBpMatching1 algorithm.
     649    ///
     650    /// \note mwbpm.run() is just a shortcut of the following code.
     651    /// \code
     652    ///   mwbpm.init();
     653    ///   mwbpm.start();
     654    /// \endcode
     655    void run() {
     656      init();
     657      start();
     658    }
     659
     660    /// @}
     661
     662    /// \name Primal Solution
     663    /// Functions to get the primal solution, i.e. the maximum weighted
     664    /// bipartite matching.\n
     665    /// Either \ref run() or \ref start() function should be called before
     666    /// using them.
     667
     668    /// @{
     669
     670    /// \brief Return the weight of the matching.
     671    ///
     672    /// This function returns the weight of the found matching.
     673    ///
     674    /// \pre Either run() or start() must be called before using this function.
     675    Value matchingWeight() const {
     676      Value sum = 0;
     677      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     678        if ((*_matching)[n] != INVALID) {
     679          sum += _weight[(*_matching)[n]];
     680        }
     681      }
     682      return sum / 2;
     683    }
     684
     685    /// \brief Return the size (cardinality) of the matching.
     686    ///
     687    /// This function returns the size (cardinality) of the found matching.
     688    ///
     689    /// \pre Either run() or start() must be called before using this function.
     690    int matchingSize() const {
     691      int num = 0;
     692      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     693        if ((*_matching)[n] != INVALID) {
     694          ++num;
     695        }
     696      }
     697      return num / 2;
     698    }
     699
     700    /// \brief Return \c true if the given edge is in the matching.
     701    ///
     702    /// This function returns \c true if the given edge is in the found
     703    /// matching.
     704    ///
     705    /// \pre Either run() or start() must be called before using this function.
     706    bool matching(const Edge& edge) const {
     707      return edge == (*_matching)[_bpgraph.u(edge)];
     708    }
     709
     710    /// \brief Return the matching arc (or edge) incident to the given node.
     711    ///
     712    /// This function returns the matching arc (or edge) incident to the
     713    /// given node in the found matching or \c INVALID if the node is
     714    /// not covered by the matching.
     715    ///
     716    /// \pre Either run() or start() must be called before using this function.
     717    Arc matching(const Node& node) const {
     718      return (*_matching)[node];
     719    }
     720
     721    /// \brief Return the mate of the given node.
     722    ///
     723    /// This function returns the mate of the given node in the found
     724    /// matching or \c INVALID if the node is not covered by the matching.
     725    ///
     726    /// \pre Either run() or start() must be called before using this function.
     727    Node mate(const Node& node) const {
     728      return (*_matching)[node] != INVALID ?
     729        _bpgraph.target((*_matching)[node]) : INVALID;
     730    }
     731
     732    /// @}
     733
     734    /// \name Dual Solution
     735    /// Functions to get the dual solution.\n
     736    /// Either \ref run() or \ref start() function should be called before
     737    /// using them.
     738
     739    /// @{
     740
     741    /// \brief Return the value of the dual solution.
     742    ///
     743    /// This function returns the value of the dual solution.
     744    /// It should be equal to the primal value scaled by \ref dualScale
     745    /// "dual scale".
     746    ///
     747    /// \pre Either run() or start() must be called before using this function.
     748    Value dualValue() const {
     749      Value sum = 0;
     750      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     751        sum += nodeValue(n);
     752      }
     753      return sum;
     754    }
     755
     756    /// \brief Return the dual value (potential) of the given node.
     757    ///
     758    /// This function returns the dual value (potential) of the given node.
     759    ///
     760    /// \pre Either run() or start() must be called before using this function.
     761    Value nodeValue(const Node& n) const {
     762      return (*_node_potential)[n];
     763    }
     764
     765    /// @}
     766
     767  };
     768
     769#ifdef DOXYGEN
     770  template <typename BPGR, typename WM>
     771#else
     772  template <typename BPGR,
     773            typename WM = typename BPGR::template EdgeMap<int> >
     774#endif
     775  class MaxWeightedBpMatching2 {
     776  public:
     777
     778    /// The graph type of the algorithm
     779    typedef BPGR BpGraph;
     780    /// The type of the edge weight map
     781    typedef WM WeightMap;
     782    /// The value type of the edge weights
     783    typedef typename WeightMap::Value Value;
     784
     785    /// The type of the matching map
     786    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     787    MatchingMap;
     788
     789    /// \brief Scaling factor for dual solution
     790    ///
     791    /// Scaling factor for dual solution. It is equal to 4 or 1
     792    /// according to the value type.
     793    static const int dualScale = 1;
     794
     795  private:
     796
     797    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     798
     799    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     800
     801    const BpGraph& _bpgraph;
     802    const WeightMap& _weight;
     803
     804    MatchingMap* _matching;
     805
     806    NodePotential* _node_potential;
     807
     808    int _node_num;
     809
     810    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     811    PredMap* _pred;
     812
     813    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     814    IntNodeMap *_tree_set_index;
     815    TreeSet *_tree_set;
     816
     817    IntNodeMap *_heap_index;
     818    BinHeap<Value, IntNodeMap> *_heap;
     819
     820    Value _delta_sum;
     821    int _unmatched;
     822
     823    void createStructures() {
     824      _node_num = countNodes(_bpgraph);
     825
     826      if (!_matching) {
     827        _matching = new MatchingMap(_bpgraph);
     828      }
     829
     830      if (!_node_potential) {
     831        _node_potential = new NodePotential(_bpgraph);
     832      }
     833
     834      if (!_pred) {
     835        _pred = new PredMap(_bpgraph);
     836      }
     837
     838      if (!_tree_set) {
     839        _tree_set_index = new IntNodeMap(_bpgraph);
     840        _tree_set = new TreeSet(*_tree_set_index);
     841      }
     842
     843      if (!_heap) {
     844        _heap_index = new IntNodeMap(_bpgraph);
     845        _heap = new BinHeap<Value, IntNodeMap>(*_heap_index);
     846      }
     847
     848    }
     849
     850    void destroyStructures() {
     851      if (_matching) {
     852        delete _matching;
     853      }
     854      if (_node_potential) {
     855        delete _node_potential;
     856      }
     857      if (_pred) {
     858        delete _pred;
     859      }
     860      if (_tree_set) {
     861        delete _tree_set_index;
     862        delete _tree_set;
     863      }
     864      if (_heap) {
     865        delete _heap_index;
     866        delete _heap;
     867      }
     868    }
     869
     870    void matchedToEven(Node node, int tree) {
     871      _tree_set->insert(node, tree);
     872      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     873      _heap->push(node, (*_node_potential)[node]);
     874
     875      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     876        Node v = _bpgraph.source(a);
     877        Value rw = (*_node_potential)[node] + (*_node_potential)[v] - _weight[a];
     878        if (_heap->state(v) == _heap->PRE_HEAP) {
     879          _pred->set(v, a);
     880          _heap->push(v, rw);
     881        } else if (_heap->state(v) == _heap->IN_HEAP && (*_heap)[v] > rw) {
     882          _pred->set(v, a);
     883          _heap->decrease(v, rw);
     884        }
     885      }
     886    }
     887
     888    void matchedToOdd(Node node, int tree) {
     889      _tree_set->insert(node, tree);
     890      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     891
     892      if (_heap->state(node) == _heap->IN_HEAP) {
     893        _heap->erase(node);
     894      } else {
     895        _heap_index->set(node, _heap->POST_HEAP);
     896      }
     897    }
     898
     899    void evenToMatched(Node node, int tree) {
     900      _heap->erase(node);
     901      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     902      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     903        Node v = _bpgraph.source(a);
     904        if (_heap->state(v) == _heap->IN_HEAP && (*_pred)[v] == a) {
     905          Arc mina = INVALID;
     906          Value minrwa = std::numeric_limits<Value>::max();
     907          for (OutArcIt aa(_bpgraph, v); aa != INVALID; ++aa) {
     908            Node va = _bpgraph.target(aa);
     909            if (_heap->state(va) == _heap->POST_HEAP ||
     910                _tree_set->find(va) == tree) continue;
     911            Value rwa = (*_node_potential)[v] + (*_node_potential)[va] - _weight[aa];
     912            if (minrwa > rwa) {
     913              minrwa = rwa;
     914              mina = aa;
     915            }
     916          }
     917          if (mina != INVALID) {
     918            _pred->set(v, mina);
     919            _heap->increase(v, minrwa);
     920          } else {
     921            _pred->set(v, INVALID);
     922            _heap->erase(v);
     923            _heap_index->set(v, _heap->PRE_HEAP);
     924          }
     925        }
     926      }
     927    }
     928
     929    void oddToMatched(Node node) {
     930      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     931      Arc min = INVALID;
     932      Value minrw = std::numeric_limits<Value>::max();
     933      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     934        Node v = _bpgraph.source(a);
     935        if (_heap->state(v) == _heap->POST_HEAP) continue;
     936        Value rw = (*_node_potential)[node] + (*_node_potential)[v] - _weight[a];
     937
     938        if (minrw > rw) {
     939          min = _bpgraph.oppositeArc(a);
     940          minrw = rw;
     941        }
     942      }
     943      if (min != INVALID) {
     944        _pred->set(node, min);
     945        _heap->push(node, minrw);
     946      } else {
     947        _pred->set(node, INVALID);
     948        _heap_index->set(node, _heap->PRE_HEAP);
     949      }
     950    }
     951
     952    void alternatePath(Node even, int tree) {
     953      Node odd;
     954
     955      evenToMatched(even, tree);
     956
     957      Arc prev = (*_matching)[even];
     958      while (prev != INVALID) {
     959        odd = _bpgraph.target(prev);
     960        even = _bpgraph.target((*_pred)[odd]);
     961        _matching->set(odd, (*_pred)[odd]);
     962        oddToMatched(odd);
     963
     964        prev = (*_matching)[even];
     965        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     966        evenToMatched(even, tree);
     967      }
     968    }
     969
     970    void destroyTree(int tree) {
     971      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     972        if (_bpgraph.red(n)) {
     973          if (_heap->state(n) != _heap->POST_HEAP) {
     974            evenToMatched(n, tree);
     975          }
     976        } else {
     977          if (_heap->state(n) == _heap->POST_HEAP) {
     978            oddToMatched(n);
     979          }
     980        }
     981      }
     982      _tree_set->eraseClass(tree);
     983    }
     984
     985    void unmatchNode(Node node) {
     986      int tree = _tree_set->find(node);
     987
     988      alternatePath(node, tree);
     989      destroyTree(tree);
     990
     991      _matching->set(node, INVALID);
     992    }
     993
     994    void augmentOnArc(Arc arc) {
     995      Node left = _bpgraph.source(arc);
     996      _matching->set(left, arc);
     997
     998      Node right = _bpgraph.target(arc);
     999      int right_tree = _tree_set->find(right);
     1000
     1001      alternatePath(right, right_tree);
     1002      destroyTree(right_tree);
     1003      _matching->set(right, _bpgraph.oppositeArc(arc));
     1004    }
     1005
     1006    void extendOnArc(Arc arc) {
     1007      Node base = _bpgraph.target(arc);
     1008      int tree = _tree_set->find(base);
     1009
     1010      Node odd = _bpgraph.source(arc);
     1011      _tree_set->insert(odd, tree);
     1012      matchedToOdd(odd, tree);
     1013      _pred->set(odd, arc);
     1014
     1015      Node even = _bpgraph.target((*_matching)[odd]);
     1016      _tree_set->insert(even, tree);
     1017      matchedToEven(even, tree);
     1018    }
     1019
     1020  public:
     1021
     1022    /// \brief Constructor
     1023    ///
     1024    /// Constructor.
     1025    MaxWeightedBpMatching2(const BpGraph& bpgraph, const WeightMap& weight)
     1026      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     1027        _node_potential(0), _node_num(0), _pred(0),
     1028        _tree_set_index(0), _tree_set(0),
     1029
     1030        _heap_index(0), _heap(0),
     1031
     1032        _delta_sum(), _unmatched(0)
     1033    {}
     1034
     1035    ~MaxWeightedBpMatching2() {
     1036      destroyStructures();
     1037    }
     1038
     1039    /// \brief Initialize the algorithm
     1040    ///
     1041    /// This function initializes the algorithm.
     1042    void redRootInit() {
     1043      createStructures();
     1044
     1045      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1046        _heap_index->set(n, _heap->PRE_HEAP);
     1047      }
     1048
     1049      _heap->clear();
     1050      _tree_set->clear();
     1051
     1052      _unmatched = 0;
     1053      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     1054        _matching->set(n, INVALID);
     1055
     1056        Value max = 0;
     1057        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1058          if (_weight[a] > max) {
     1059            max = _weight[a];
     1060          }
     1061        }
     1062        _node_potential->set(n, max);
     1063        _heap->push(n, max);
     1064
     1065        _tree_set->insert(n);
     1066
     1067        ++_unmatched;
     1068      }
     1069      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     1070        _matching->set(n, INVALID);
     1071        _node_potential->set(n, 0);
     1072
     1073        Arc min = INVALID;
     1074        Value minrw = std::numeric_limits<Value>::max();
     1075        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1076          Node v = _bpgraph.target(a);
     1077          Value rw = (*_node_potential)[n] + (*_node_potential)[v] - _weight[a];
     1078
     1079          if (minrw > rw) {
     1080            min = a;
     1081            minrw = rw;
     1082          }
     1083        }
     1084        if (min != INVALID) {
     1085          _pred->set(n, min);
     1086          _heap->push(n, minrw);
     1087        } else {
     1088          _pred->set(n, INVALID);
     1089        }
     1090      }
     1091    }
     1092
     1093    /// \brief Start the algorithm
     1094    ///
     1095    /// This function starts the algorithm.
     1096    ///
     1097    /// \pre \ref init() must be called before using this function.
     1098    void start() {
     1099      while (_unmatched > 0) {
     1100        _delta_sum = _heap->prio();
     1101        Node n = _heap->top();
     1102        if (_bpgraph.red(n)) {
     1103          unmatchNode(n);
     1104          --_unmatched;
     1105        } else if ((*_matching)[n] == INVALID) {
     1106          augmentOnArc((*_pred)[n]);
     1107          --_unmatched;
     1108        } else {
     1109          extendOnArc((*_pred)[n]);
     1110        }
     1111      }
     1112    }
     1113
     1114    /// \brief Run the algorithm.
     1115    ///
     1116    /// This method runs the \c %MaxWeightedBpMatching algorithm.
     1117    ///
     1118    /// \note mwbpm.run() is just a shortcut of the following code.
     1119    /// \code
     1120    ///   mwbpm.init();
     1121    ///   mwbpm.start();
     1122    /// \endcode
     1123    void run() {
     1124      redRootInit();
     1125      start();
     1126    }
     1127
     1128    /// @}
     1129
     1130    /// \name Primal Solution
     1131    /// Functions to get the primal solution, i.e. the maximum weighted
     1132    /// bipartite matching.\n
     1133    /// Either \ref run() or \ref start() function should be called before
     1134    /// using them.
     1135
     1136    /// @{
     1137
     1138    /// \brief Return the weight of the matching.
     1139    ///
     1140    /// This function returns the weight of the found matching.
     1141    ///
     1142    /// \pre Either run() or start() must be called before using this function.
     1143    Value matchingWeight() const {
     1144      Value sum = 0;
     1145      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1146        if ((*_matching)[n] != INVALID) {
     1147          sum += _weight[(*_matching)[n]];
     1148        }
     1149      }
     1150      return sum / 2;
     1151    }
     1152
     1153    /// \brief Return the size (cardinality) of the matching.
     1154    ///
     1155    /// This function returns the size (cardinality) of the found matching.
     1156    ///
     1157    /// \pre Either run() or start() must be called before using this function.
     1158    int matchingSize() const {
     1159      int num = 0;
     1160      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1161        if ((*_matching)[n] != INVALID) {
     1162          ++num;
     1163        }
     1164      }
     1165      return num / 2;
     1166    }
     1167
     1168    /// \brief Return \c true if the given edge is in the matching.
     1169    ///
     1170    /// This function returns \c true if the given edge is in the found
     1171    /// matching.
     1172    ///
     1173    /// \pre Either run() or start() must be called before using this function.
     1174    bool matching(const Edge& edge) const {
     1175      return edge == (*_matching)[_bpgraph.u(edge)];
     1176    }
     1177
     1178    /// \brief Return the matching arc (or edge) incident to the given node.
     1179    ///
     1180    /// This function returns the matching arc (or edge) incident to the
     1181    /// given node in the found matching or \c INVALID if the node is
     1182    /// not covered by the matching.
     1183    ///
     1184    /// \pre Either run() or start() must be called before using this function.
     1185    Arc matching(const Node& node) const {
     1186      return (*_matching)[node];
     1187    }
     1188
     1189    /// \brief Return the mate of the given node.
     1190    ///
     1191    /// This function returns the mate of the given node in the found
     1192    /// matching or \c INVALID if the node is not covered by the matching.
     1193    ///
     1194    /// \pre Either run() or start() must be called before using this function.
     1195    Node mate(const Node& node) const {
     1196      return (*_matching)[node] != INVALID ?
     1197        _bpgraph.target((*_matching)[node]) : INVALID;
     1198    }
     1199
     1200    /// @}
     1201
     1202    /// \name Dual Solution
     1203    /// Functions to get the dual solution.\n
     1204    /// Either \ref run() or \ref start() function should be called before
     1205    /// using them.
     1206
     1207    /// @{
     1208
     1209    /// \brief Return the value of the dual solution.
     1210    ///
     1211    /// This function returns the value of the dual solution.
     1212    /// It should be equal to the primal value scaled by \ref dualScale
     1213    /// "dual scale".
     1214    ///
     1215    /// \pre Either run() or start() must be called before using this function.
     1216    Value dualValue() const {
     1217      Value sum = 0;
     1218      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1219        sum += nodeValue(n);
     1220      }
     1221      return sum;
     1222    }
     1223
     1224    /// \brief Return the dual value (potential) of the given node.
     1225    ///
     1226    /// This function returns the dual value (potential) of the given node.
     1227    ///
     1228    /// \pre Either run() or start() must be called before using this function.
     1229    Value nodeValue(const Node& n) const {
     1230      return (*_node_potential)[n];
     1231    }
     1232
     1233    /// @}
     1234
     1235  };
     1236
     1237#ifdef DOXYGEN
     1238  template <typename BPGR, typename WM>
     1239#else
     1240  template <typename BPGR,
     1241            typename WM = typename BPGR::template EdgeMap<int> >
     1242#endif
     1243  class MaxWeightedBpMatching3 {
     1244  public:
     1245
     1246    /// The graph type of the algorithm
     1247    typedef BPGR BpGraph;
     1248    /// The type of the edge weight map
     1249    typedef WM WeightMap;
     1250    /// The value type of the edge weights
     1251    typedef typename WeightMap::Value Value;
     1252
     1253    /// The type of the matching map
     1254    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     1255    MatchingMap;
     1256
     1257    /// \brief Scaling factor for dual solution
     1258    ///
     1259    /// Scaling factor for dual solution. It is equal to 4 or 1
     1260    /// according to the value type.
     1261    static const int dualScale =
     1262      std::numeric_limits<Value>::is_integer ? 4 : 1;
     1263
     1264  private:
     1265
     1266    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     1267
     1268    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     1269
     1270    const BpGraph& _bpgraph;
     1271    const WeightMap& _weight;
     1272
     1273    MatchingMap* _matching;
     1274
     1275    NodePotential* _node_potential;
     1276
     1277    int _node_num;
     1278
     1279    enum Status {
     1280      EVEN = -1, MATCHED = 0, ODD = 1
     1281    };
     1282
     1283    typedef typename BpGraph::template NodeMap<Status> StatusMap;
     1284    StatusMap* _status;
     1285
     1286    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     1287    PredMap* _pred;
     1288
     1289    IntArcMap* _node_heap_index;
     1290    typedef typename BpGraph::template NodeMap<BinHeap<Value, IntArcMap>*> NodeHeapMap;
     1291    NodeHeapMap* _node_heap;
     1292
     1293    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     1294    IntNodeMap *_tree_set_index;
     1295    TreeSet *_tree_set;
     1296
     1297    IntNodeMap *_delta1_index;
     1298    BinHeap<Value, IntNodeMap> *_delta1;
     1299
     1300    IntNodeMap *_delta2_index;
     1301    BinHeap<Value, IntNodeMap> *_delta2;
     1302
     1303    IntEdgeMap *_delta3_index;
     1304    BinHeap<Value, IntEdgeMap> *_delta3;
     1305
     1306    Value _delta_sum;
     1307    int _unmatched;
     1308
     1309    void createStructures() {
     1310      _node_num = countNodes(_bpgraph);
     1311
     1312      if (!_matching) {
     1313        _matching = new MatchingMap(_bpgraph);
     1314      }
     1315
     1316      if (!_node_potential) {
     1317        _node_potential = new NodePotential(_bpgraph);
     1318      }
     1319
     1320      if (!_status) {
     1321        _status = new StatusMap(_bpgraph);
     1322      }
     1323
     1324      if (!_pred) {
     1325        _pred = new PredMap(_bpgraph);
     1326      }
     1327
     1328      if (!_node_heap) {
     1329        _node_heap_index = new IntArcMap(_bpgraph);
     1330        _node_heap = new NodeHeapMap(_bpgraph, 0);
     1331      } else {
     1332        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1333          if ((*_status)[n] == MATCHED) {
     1334            delete (*_node_heap)[n];
     1335          }
     1336        }
     1337      }
     1338
     1339      if (!_tree_set) {
     1340        _tree_set_index = new IntNodeMap(_bpgraph);
     1341        _tree_set = new TreeSet(*_tree_set_index);
     1342      }
     1343
     1344      if (!_delta1) {
     1345        _delta1_index = new IntNodeMap(_bpgraph);
     1346        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
     1347      }
     1348
     1349      if (!_delta2) {
     1350        _delta2_index = new IntNodeMap(_bpgraph);
     1351        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
     1352      }
     1353
     1354      if (!_delta3) {
     1355        _delta3_index = new IntEdgeMap(_bpgraph);
     1356        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
     1357      }
     1358    }
     1359
     1360    void destroyStructures() {
     1361      if (_matching) {
     1362        delete _matching;
     1363      }
     1364      if (_node_potential) {
     1365        delete _node_potential;
     1366      }
     1367      if (_node_heap) {
     1368        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1369          if ((*_status)[n] == MATCHED) {
     1370            delete (*_node_heap)[n];
     1371          }
     1372        }
     1373        delete _node_heap_index;
     1374        delete _node_heap;
     1375      }
     1376      if (_status) {
     1377        delete _status;
     1378      }
     1379      if (_pred) {
     1380        delete _pred;
     1381      }
     1382      if (_tree_set) {
     1383        delete _tree_set_index;
     1384        delete _tree_set;
     1385      }
     1386      if (_delta2) {
     1387        delete _delta2_index;
     1388        delete _delta2;
     1389      }
     1390      if (_delta3) {
     1391        delete _delta3_index;
     1392        delete _delta3;
     1393      }
     1394    }
     1395
     1396    void matchedToEven(Node node, int tree) {
     1397      _tree_set->insert(node, tree);
     1398      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     1399      _delta1->push(node, (*_node_potential)[node]);
     1400
     1401      if (_delta2->state(node) == _delta2->IN_HEAP) {
     1402        _delta2->erase(node);
     1403      }
     1404      delete (*_node_heap)[node];
     1405
     1406      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     1407        Node v = _bpgraph.source(a);
     1408        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     1409          dualScale * _weight[a];
     1410        if ((*_status)[v] == EVEN) {
     1411          _delta3->push(a, rw / 2);
     1412        } else if ((*_status)[v] == MATCHED) {
     1413          (*_node_heap)[v]->push(a, rw);
     1414          if (_delta2->state(v) != _delta2->IN_HEAP) {
     1415            _delta2->push(v, rw);
     1416          } else if ((*_delta2)[v] > rw) {
     1417            _delta2->decrease(v, rw);
     1418          }
     1419        }
     1420      }
     1421    }
     1422
     1423    void matchedToOdd(Node node, int tree) {
     1424      _tree_set->insert(node, tree);
     1425      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     1426
     1427      if (_delta2->state(node) == _delta2->IN_HEAP) {
     1428        _delta2->erase(node);
     1429        delete (*_node_heap)[node];
     1430      }
     1431    }
     1432
     1433    void evenToMatched(Node node) {
     1434      _delta1->erase(node);
     1435      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     1436      (*_node_heap)[node] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     1437      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     1438        Node v = _bpgraph.source(a);
     1439        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     1440          dualScale * _weight[a];
     1441
     1442        if ((*_status)[v] == EVEN) {
     1443          _delta3->erase(a);
     1444
     1445          (*_node_heap)[node]->push(_bpgraph.oppositeArc(a), rw);
     1446        } else if ((*_status)[v]  == MATCHED) {
     1447          (*_node_heap)[v]->erase(a);
     1448          if ((*_node_heap)[v]->empty()) {
     1449            _delta2->erase(v);
     1450          } else if ((*_node_heap)[v]->prio() > (*_delta2)[v]) {
     1451            _delta2->increase(v, (*_node_heap)[v]->prio());
     1452          }
     1453        }
     1454      }
     1455      if (!(*_node_heap)[node]->empty()) {
     1456        _delta2->push(node, (*_node_heap)[node]->prio());
     1457      }
     1458    }
     1459
     1460    void oddToMatched(Node node) {
     1461      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     1462      (*_node_heap)[node] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     1463      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     1464        Node v = _bpgraph.source(a);
     1465        if ((*_status)[v] != EVEN) continue;
     1466        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
     1467          dualScale * _weight[a];
     1468        (*_node_heap)[node]->push(_bpgraph.oppositeArc(a), rw);
     1469      }
     1470      if (!(*_node_heap)[node]->empty()) {
     1471        _delta2->push(node, (*_node_heap)[node]->prio());
     1472      }
     1473    }
     1474
     1475    void alternatePath(Node even) {
     1476      Node odd;
     1477
     1478      _status->set(even, MATCHED);
     1479      evenToMatched(even);
     1480
     1481      Arc prev = (*_matching)[even];
     1482      while (prev != INVALID) {
     1483        odd = _bpgraph.target(prev);
     1484        even = _bpgraph.target((*_pred)[odd]);
     1485        _matching->set(odd, (*_pred)[odd]);
     1486        _status->set(odd, MATCHED);
     1487        oddToMatched(odd);
     1488
     1489        prev = (*_matching)[even];
     1490        _status->set(even, MATCHED);
     1491        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     1492        evenToMatched(even);
     1493      }
     1494    }
     1495
     1496    void destroyTree(int tree) {
     1497      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     1498        if ((*_status)[n] == EVEN) {
     1499          _status->set(n, MATCHED);
     1500          evenToMatched(n);
     1501        } else if ((*_status)[n] == ODD) {
     1502          _status->set(n, MATCHED);
     1503          oddToMatched(n);
     1504        }
     1505      }
     1506      _tree_set->eraseClass(tree);
     1507    }
     1508
     1509    void unmatchNode(const Node& node) {
     1510      int tree = _tree_set->find(node);
     1511
     1512      alternatePath(node);
     1513      destroyTree(tree);
     1514
     1515      _matching->set(node, INVALID);
     1516    }
     1517
     1518    void augmentOnEdge(const Edge& edge) {
     1519      Node left = _bpgraph.u(edge);
     1520      int left_tree = _tree_set->find(left);
     1521
     1522      alternatePath(left);
     1523      destroyTree(left_tree);
     1524      _matching->set(left, _bpgraph.direct(edge, true));
     1525
     1526      Node right = _bpgraph.v(edge);
     1527      int right_tree = _tree_set->find(right);
     1528
     1529      alternatePath(right);
     1530      destroyTree(right_tree);
     1531      _matching->set(right, _bpgraph.direct(edge, false));
     1532    }
     1533
     1534    void augmentOnArc(const Arc& arc) {
     1535      Node left = _bpgraph.source(arc);
     1536      _matching->set(left, arc);
     1537      _pred->set(left, arc);
     1538
     1539      Node right = _bpgraph.target(arc);
     1540      int right_tree = _tree_set->find(right);
     1541
     1542      alternatePath(right);
     1543      destroyTree(right_tree);
     1544      _matching->set(right, _bpgraph.oppositeArc(arc));
     1545    }
     1546
     1547    void extendOnArc(const Arc& arc) {
     1548      Node base = _bpgraph.target(arc);
     1549      int tree = _tree_set->find(base);
     1550
     1551      Node odd = _bpgraph.source(arc);
     1552      _tree_set->insert(odd, tree);
     1553      _status->set(odd, ODD);
     1554      matchedToOdd(odd, tree);
     1555      _pred->set(odd, arc);
     1556
     1557      Node even = _bpgraph.target((*_matching)[odd]);
     1558      _tree_set->insert(even, tree);
     1559      _status->set(even, EVEN);
     1560      matchedToEven(even, tree);
     1561    }
     1562
     1563  public:
     1564
     1565    /// \brief Constructor
     1566    ///
     1567    /// Constructor.
     1568    MaxWeightedBpMatching3(const BpGraph& bpgraph, const WeightMap& weight)
     1569      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     1570        _node_potential(0), _node_num(0),
     1571        _status(0), _pred(0), _node_heap_index(0), _node_heap(0),
     1572        _tree_set_index(0), _tree_set(0),
     1573
     1574        _delta1_index(0), _delta1(0),
     1575        _delta2_index(0), _delta2(0),
     1576        _delta3_index(0), _delta3(0),
     1577
     1578        _delta_sum(), _unmatched(0)
     1579    {}
     1580
     1581    ~MaxWeightedBpMatching3() {
     1582      destroyStructures();
     1583    }
     1584
     1585    /// \name Execution Control
     1586    /// The simplest way to execute the algorithm is to use the
     1587    /// \ref run() member function.
     1588
     1589    ///@{
     1590
     1591    /// \brief Initialize the algorithm
     1592    ///
     1593    /// This function initializes the algorithm.
     1594    void init() {
     1595      createStructures();
     1596
     1597      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1598        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     1599        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     1600      }
     1601      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1602        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     1603      }
     1604
     1605      _delta1->clear();
     1606      _delta2->clear();
     1607      _delta3->clear();
     1608      _tree_set->clear();
     1609
     1610      _unmatched = _node_num;
     1611
     1612      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1613        Value max = 0;
     1614        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1615          if ((dualScale * _weight[a]) / 2 > max) {
     1616            max = (dualScale * _weight[a]) / 2;
     1617          }
     1618        }
     1619        _node_potential->set(n, max);
     1620        _delta1->push(n, max);
     1621
     1622        _tree_set->insert(n);
     1623
     1624        _matching->set(n, INVALID);
     1625        _status->set(n, EVEN);
     1626      }
     1627      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1628        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
     1629                          (*_node_potential)[_bpgraph.v(e)] -
     1630                          dualScale * _weight[e]) / 2);
     1631      }
     1632    }
     1633
     1634    /// \brief Initialize the algorithm
     1635    ///
     1636    /// This function initializes the algorithm.
     1637    void redRootInit() {
     1638      createStructures();
     1639
     1640      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1641        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     1642        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     1643      }
     1644      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1645        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     1646      }
     1647
     1648      _delta1->clear();
     1649      _delta2->clear();
     1650      _delta3->clear();
     1651      _tree_set->clear();
     1652
     1653      _unmatched = 0;
     1654      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     1655        Value max = 0;
     1656        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1657          if ((dualScale * _weight[a]) > max) {
     1658            max = dualScale * _weight[a];
     1659          }
     1660        }
     1661        _node_potential->set(n, max);
     1662        _delta1->push(n, max);
     1663
     1664        _tree_set->insert(n);
     1665
     1666        _matching->set(n, INVALID);
     1667        _status->set(n, EVEN);
     1668
     1669        ++_unmatched;
     1670      }
     1671      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     1672        _matching->set(n, INVALID);
     1673        _status->set(n, MATCHED);
     1674        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     1675        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1676          Node v = _bpgraph.target(a);
     1677          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     1678                     dualScale * _weight[a];
     1679          (*_node_heap)[n]->push(a, rw);
     1680        }
     1681        if (!(*_node_heap)[n]->empty()) {
     1682          _delta2->push(n, (*_node_heap)[n]->prio());
     1683        }
     1684      }
     1685    }
     1686
     1687    /// \brief Initialize the algorithm
     1688    ///
     1689    /// This function initializes the algorithm.
     1690    void blueRootInit() {
     1691      createStructures();
     1692
     1693      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1694        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     1695        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     1696      }
     1697      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     1698        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     1699      }
     1700
     1701      _delta1->clear();
     1702      _delta2->clear();
     1703      _delta3->clear();
     1704      _tree_set->clear();
     1705
     1706      _unmatched = 0;
     1707      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     1708        Value max = 0;
     1709        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1710          if ((dualScale * _weight[a]) > max) {
     1711            max = dualScale * _weight[a];
     1712          }
     1713        }
     1714        _node_potential->set(n, max);
     1715        _delta1->push(n, max);
     1716
     1717        _tree_set->insert(n);
     1718
     1719        _matching->set(n, INVALID);
     1720        _status->set(n, EVEN);
     1721
     1722        ++_unmatched;
     1723      }
     1724      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     1725        _matching->set(n, INVALID);
     1726        _status->set(n, MATCHED);
     1727        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     1728        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     1729          Node v = _bpgraph.target(a);
     1730          Value rw = (*_node_potential)[n] + (*_node_potential)[v] -
     1731                     dualScale * _weight[a];
     1732          (*_node_heap)[n]->push(a, rw);
     1733        }
     1734        if (!(*_node_heap)[n]->empty()) {
     1735          _delta2->push(n, (*_node_heap)[n]->prio());
     1736        }
     1737      }
     1738    }
     1739
     1740    /// \brief Start the algorithm
     1741    ///
     1742    /// This function starts the algorithm.
     1743    ///
     1744    /// \pre \ref init() must be called before using this function.
     1745    void start() {
     1746      enum OpType {
     1747        D1, D2, D3
     1748      };
     1749
     1750      while (_unmatched > 0) {
     1751        Value d1 = !_delta1->empty() ?
     1752          _delta1->prio() : std::numeric_limits<Value>::max();
     1753
     1754        Value d2 = !_delta2->empty() ?
     1755          _delta2->prio() : std::numeric_limits<Value>::max();
     1756
     1757        Value d3 = !_delta3->empty() ?
     1758          _delta3->prio() : std::numeric_limits<Value>::max();
     1759
     1760        _delta_sum = d3; OpType ot = D3;
     1761        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
     1762        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     1763
     1764        switch (ot) {
     1765        case D1:
     1766          {
     1767            Node n = _delta1->top();
     1768            unmatchNode(n);
     1769            --_unmatched;
     1770          }
     1771          break;
     1772        case D2:
     1773          {
     1774            Node n = _delta2->top();
     1775            Arc a = (*_node_heap)[n]->top();
     1776            if ((*_matching)[n] == INVALID) {
     1777              augmentOnArc(a);
     1778              --_unmatched;
     1779            } else {
     1780              extendOnArc(a);
     1781            }
     1782          }
     1783          break;
     1784        case D3:
     1785          {
     1786            Edge e = _delta3->top();
     1787            augmentOnEdge(e);
     1788            _unmatched -= 2;
     1789          }
     1790          break;
     1791        }
     1792      }
     1793    }
     1794
     1795    /// \brief Run the algorithm.
     1796    ///
     1797    /// This method runs the \c %MaxWeightedBpMatching3 algorithm.
     1798    ///
     1799    /// \note mwbpm.run() is just a shortcut of the following code.
     1800    /// \code
     1801    ///   mwbpm.init();
     1802    ///   mwbpm.start();
     1803    /// \endcode
     1804    void run() {
     1805      init();
     1806      start();
     1807    }
     1808
     1809    /// @}
     1810
     1811    /// \name Primal Solution
     1812    /// Functions to get the primal solution, i.e. the maximum weighted
     1813    /// bipartite matching.\n
     1814    /// Either \ref run() or \ref start() function should be called before
     1815    /// using them.
     1816
     1817    /// @{
     1818
     1819    /// \brief Return the weight of the matching.
     1820    ///
     1821    /// This function returns the weight of the found matching.
     1822    ///
     1823    /// \pre Either run() or start() must be called before using this function.
     1824    Value matchingWeight() const {
     1825      Value sum = 0;
     1826      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1827        if ((*_matching)[n] != INVALID) {
     1828          sum += _weight[(*_matching)[n]];
     1829        }
     1830      }
     1831      return sum / 2;
     1832    }
     1833
     1834    /// \brief Return the size (cardinality) of the matching.
     1835    ///
     1836    /// This function returns the size (cardinality) of the found matching.
     1837    ///
     1838    /// \pre Either run() or start() must be called before using this function.
     1839    int matchingSize() const {
     1840      int num = 0;
     1841      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1842        if ((*_matching)[n] != INVALID) {
     1843          ++num;
     1844        }
     1845      }
     1846      return num / 2;
     1847    }
     1848
     1849    /// \brief Return \c true if the given edge is in the matching.
     1850    ///
     1851    /// This function returns \c true if the given edge is in the found
     1852    /// matching.
     1853    ///
     1854    /// \pre Either run() or start() must be called before using this function.
     1855    bool matching(const Edge& edge) const {
     1856      return edge == (*_matching)[_bpgraph.u(edge)];
     1857    }
     1858
     1859    /// \brief Return the matching arc (or edge) incident to the given node.
     1860    ///
     1861    /// This function returns the matching arc (or edge) incident to the
     1862    /// given node in the found matching or \c INVALID if the node is
     1863    /// not covered by the matching.
     1864    ///
     1865    /// \pre Either run() or start() must be called before using this function.
     1866    Arc matching(const Node& node) const {
     1867      return (*_matching)[node];
     1868    }
     1869
     1870    /// \brief Return the mate of the given node.
     1871    ///
     1872    /// This function returns the mate of the given node in the found
     1873    /// matching or \c INVALID if the node is not covered by the matching.
     1874    ///
     1875    /// \pre Either run() or start() must be called before using this function.
     1876    Node mate(const Node& node) const {
     1877      return (*_matching)[node] != INVALID ?
     1878        _bpgraph.target((*_matching)[node]) : INVALID;
     1879    }
     1880
     1881    /// @}
     1882
     1883    /// \name Dual Solution
     1884    /// Functions to get the dual solution.\n
     1885    /// Either \ref run() or \ref start() function should be called before
     1886    /// using them.
     1887
     1888    /// @{
     1889
     1890    /// \brief Return the value of the dual solution.
     1891    ///
     1892    /// This function returns the value of the dual solution.
     1893    /// It should be equal to the primal value scaled by \ref dualScale
     1894    /// "dual scale".
     1895    ///
     1896    /// \pre Either run() or start() must be called before using this function.
     1897    Value dualValue() const {
     1898      Value sum = 0;
     1899      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     1900        sum += nodeValue(n);
     1901      }
     1902      return sum;
     1903    }
     1904
     1905    /// \brief Return the dual value (potential) of the given node.
     1906    ///
     1907    /// This function returns the dual value (potential) of the given node.
     1908    ///
     1909    /// \pre Either run() or start() must be called before using this function.
     1910    Value nodeValue(const Node& n) const {
     1911      return (*_node_potential)[n];
     1912    }
     1913
     1914    /// @}
     1915
     1916  };
     1917
     1918#ifdef DOXYGEN
     1919  template <typename BPGR, typename WM>
     1920#else
     1921  template <typename BPGR,
     1922            typename WM = typename BPGR::template EdgeMap<int> >
     1923#endif
     1924  class MaxWeightedBpMatching4 {
     1925  public:
     1926
     1927    /// The graph type of the algorithm
     1928    typedef BPGR BpGraph;
     1929    /// The type of the edge weight map
     1930    typedef WM WeightMap;
     1931    /// The value type of the edge weights
     1932    typedef typename WeightMap::Value Value;
     1933
     1934    /// The type of the matching map
     1935    typedef typename BpGraph::template NodeMap<typename BpGraph::Arc>
     1936    MatchingMap;
     1937
     1938    /// \brief Scaling factor for dual solution
     1939    ///
     1940    /// Scaling factor for dual solution. It is equal to 4 or 1
     1941    /// according to the value type.
     1942    static const int dualScale =
     1943      std::numeric_limits<Value>::is_integer ? 4 : 1;
     1944
     1945  private:
     1946
     1947    TEMPLATE_BPGRAPH_TYPEDEFS(BpGraph);
     1948
     1949    typedef typename BpGraph::template NodeMap<Value> NodePotential;
     1950
     1951    const BpGraph& _bpgraph;
     1952    const WeightMap& _weight;
     1953
     1954    MatchingMap* _matching;
     1955
     1956    NodePotential* _node_potential;
     1957
     1958    int _node_num;
     1959
     1960    enum Status {
     1961      EVEN = -1, MATCHED = 0, ODD = 1
     1962    };
     1963
     1964    typedef typename BpGraph::template NodeMap<Status> StatusMap;
     1965    StatusMap* _status;
     1966
     1967    typedef typename BpGraph::template NodeMap<Arc> PredMap;
     1968    PredMap* _pred;
     1969
     1970    IntArcMap* _node_heap_index;
     1971    typedef typename BpGraph::template NodeMap<BinHeap<Value, IntArcMap>*> NodeHeapMap;
     1972    NodeHeapMap* _node_heap;
     1973
     1974    typedef ExtendFindEnum<IntNodeMap> TreeSet;
     1975    IntNodeMap *_tree_set_index;
     1976    TreeSet *_tree_set;
     1977
     1978    IntNodeMap *_delta1_index;
     1979    BinHeap<Value, IntNodeMap> *_delta1;
     1980
     1981    IntNodeMap *_delta2_index;
     1982    BinHeap<Value, IntNodeMap> *_delta2;
     1983
     1984    IntEdgeMap *_delta3_index;
     1985    BinHeap<Value, IntEdgeMap> *_delta3;
     1986
     1987    Value _delta_sum;
     1988    int _unmatched;
     1989
     1990    void createStructures() {
     1991      _node_num = countNodes(_bpgraph);
     1992
     1993      if (!_matching) {
     1994        _matching = new MatchingMap(_bpgraph);
     1995      }
     1996
     1997      if (!_node_potential) {
     1998        _node_potential = new NodePotential(_bpgraph);
     1999      }
     2000
     2001      if (!_status) {
     2002        _status = new StatusMap(_bpgraph);
     2003      }
     2004
     2005      if (!_pred) {
     2006        _pred = new PredMap(_bpgraph);
     2007      }
     2008
     2009      if (!_node_heap) {
     2010        _node_heap_index = new IntArcMap(_bpgraph);
     2011        _node_heap = new NodeHeapMap(_bpgraph, 0);
     2012      }
     2013
     2014      if (!_tree_set) {
     2015        _tree_set_index = new IntNodeMap(_bpgraph);
     2016        _tree_set = new TreeSet(*_tree_set_index);
     2017      }
     2018
     2019      if (!_delta1) {
     2020        _delta1_index = new IntNodeMap(_bpgraph);
     2021        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
     2022      }
     2023
     2024      if (!_delta2) {
     2025        _delta2_index = new IntNodeMap(_bpgraph);
     2026        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
     2027      }
     2028
     2029      if (!_delta3) {
     2030        _delta3_index = new IntEdgeMap(_bpgraph);
     2031        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
     2032      }
     2033    }
     2034
     2035    void destroyStructures() {
     2036      if (_matching) {
     2037        delete _matching;
     2038      }
     2039      if (_node_potential) {
     2040        delete _node_potential;
     2041      }
     2042      if (_node_heap) {
     2043        for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2044          delete (*_node_heap)[n];
     2045        }
     2046        delete _node_heap_index;
     2047        delete _node_heap;
     2048      }
     2049      if (_status) {
     2050        delete _status;
     2051      }
     2052      if (_pred) {
     2053        delete _pred;
     2054      }
     2055      if (_tree_set) {
     2056        delete _tree_set_index;
     2057        delete _tree_set;
     2058      }
     2059      if (_delta2) {
     2060        delete _delta2_index;
     2061        delete _delta2;
     2062      }
     2063      if (_delta3) {
     2064        delete _delta3_index;
     2065        delete _delta3;
     2066      }
     2067    }
     2068
     2069    void matchedToEven(Node node, int tree) {
     2070      _tree_set->insert(node, tree);
     2071      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     2072      _delta1->push(node, (*_node_potential)[node]);
     2073
     2074      if (_delta2->state(node) == _delta2->IN_HEAP) {
     2075        _delta2->erase(node);
     2076      }
     2077
     2078      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     2079        Node v = _bpgraph.source(a);
     2080        Value hrw = (*_node_potential)[node] - dualScale * _weight[a];
     2081        (*_node_heap)[v]->push(a, hrw);
     2082        Value rw = hrw + (*_node_potential)[v];
     2083        if ((*_status)[v] == EVEN) {
     2084          _delta3->push(a, rw / 2);
     2085        } else if ((*_status)[v] == MATCHED) {
     2086          if (_delta2->state(v) != _delta2->IN_HEAP) {
     2087            _delta2->push(v, rw);
     2088          } else if ((*_delta2)[v] > rw) {
     2089            _delta2->decrease(v, rw);
     2090          }
     2091        }
     2092      }
     2093    }
     2094
     2095    void matchedToOdd(Node node, int tree) {
     2096      _tree_set->insert(node, tree);
     2097      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     2098
     2099      if (_delta2->state(node) == _delta2->IN_HEAP) {
     2100        _delta2->erase(node);
     2101      }
     2102    }
     2103
     2104    void evenToMatched(Node node) {
     2105      _delta1->erase(node);
     2106      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
     2107      for (InArcIt a(_bpgraph, node); a != INVALID; ++a) {
     2108        Node v = _bpgraph.source(a);
     2109
     2110        (*_node_heap)[v]->erase(a);
     2111
     2112        if ((*_status)[v] == EVEN) {
     2113          _delta3->erase(a);
     2114
     2115        } else if ((*_status)[v]  == MATCHED) {
     2116          if ((*_node_heap)[v]->empty()) {
     2117            _delta2->erase(v);
     2118          } else if ((*_node_heap)[v]->prio() + (*_node_potential)[v] > (*_delta2)[v]) {
     2119            _delta2->increase(v, (*_node_heap)[v]->prio() + (*_node_potential)[v]);
     2120          }
     2121        }
     2122      }
     2123      if (!(*_node_heap)[node]->empty()) {
     2124        _delta2->push(
     2125            node, (*_node_heap)[node]->prio() + (*_node_potential)[node]);
     2126      }
     2127    }
     2128
     2129    void oddToMatched(Node node) {
     2130      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
     2131      if (!(*_node_heap)[node]->empty()) {
     2132        _delta2->push(node, (*_node_heap)[node]->prio() + (*_node_potential)[node]);
     2133      }
     2134    }
     2135
     2136    void alternatePath(Node even) {
     2137      Node odd;
     2138
     2139      _status->set(even, MATCHED);
     2140      evenToMatched(even);
     2141
     2142      Arc prev = (*_matching)[even];
     2143      while (prev != INVALID) {
     2144        odd = _bpgraph.target(prev);
     2145        even = _bpgraph.target((*_pred)[odd]);
     2146        _matching->set(odd, (*_pred)[odd]);
     2147        _status->set(odd, MATCHED);
     2148        oddToMatched(odd);
     2149
     2150        prev = (*_matching)[even];
     2151        _status->set(even, MATCHED);
     2152        _matching->set(even, _bpgraph.oppositeArc((*_matching)[odd]));
     2153        evenToMatched(even);
     2154      }
     2155    }
     2156
     2157    void destroyTree(int tree) {
     2158      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
     2159        if ((*_status)[n] == EVEN) {
     2160          _status->set(n, MATCHED);
     2161          evenToMatched(n);
     2162        } else if ((*_status)[n] == ODD) {
     2163          _status->set(n, MATCHED);
     2164          oddToMatched(n);
     2165        }
     2166      }
     2167      _tree_set->eraseClass(tree);
     2168    }
     2169
     2170    void unmatchNode(const Node& node) {
     2171      int tree = _tree_set->find(node);
     2172
     2173      alternatePath(node);
     2174      destroyTree(tree);
     2175
     2176      _matching->set(node, INVALID);
     2177    }
     2178
     2179    void augmentOnEdge(const Edge& edge) {
     2180      Node left = _bpgraph.u(edge);
     2181      int left_tree = _tree_set->find(left);
     2182
     2183      alternatePath(left);
     2184      destroyTree(left_tree);
     2185      _matching->set(left, _bpgraph.direct(edge, true));
     2186
     2187      Node right = _bpgraph.v(edge);
     2188      int right_tree = _tree_set->find(right);
     2189
     2190      alternatePath(right);
     2191      destroyTree(right_tree);
     2192      _matching->set(right, _bpgraph.direct(edge, false));
     2193    }
     2194
     2195    void augmentOnArc(const Arc& arc) {
     2196      Node left = _bpgraph.source(arc);
     2197      _matching->set(left, arc);
     2198      _pred->set(left, arc);
     2199
     2200      Node right = _bpgraph.target(arc);
     2201      int right_tree = _tree_set->find(right);
     2202
     2203      alternatePath(right);
     2204      destroyTree(right_tree);
     2205      _matching->set(right, _bpgraph.oppositeArc(arc));
     2206    }
     2207
     2208    void extendOnArc(const Arc& arc) {
     2209      Node base = _bpgraph.target(arc);
     2210      int tree = _tree_set->find(base);
     2211
     2212      Node odd = _bpgraph.source(arc);
     2213      _tree_set->insert(odd, tree);
     2214      _status->set(odd, ODD);
     2215      matchedToOdd(odd, tree);
     2216      _pred->set(odd, arc);
     2217
     2218      Node even = _bpgraph.target((*_matching)[odd]);
     2219      _tree_set->insert(even, tree);
     2220      _status->set(even, EVEN);
     2221      matchedToEven(even, tree);
     2222    }
     2223
     2224  public:
     2225
     2226    /// \brief Constructor
     2227    ///
     2228    /// Constructor.
     2229    MaxWeightedBpMatching4(const BpGraph& bpgraph, const WeightMap& weight)
     2230      : _bpgraph(bpgraph), _weight(weight), _matching(0),
     2231        _node_potential(0), _node_num(0),
     2232        _status(0), _pred(0), _node_heap_index(0), _node_heap(0),
     2233        _tree_set_index(0), _tree_set(0),
     2234
     2235        _delta1_index(0), _delta1(0),
     2236        _delta2_index(0), _delta2(0),
     2237        _delta3_index(0), _delta3(0),
     2238
     2239        _delta_sum(), _unmatched(0)
     2240    {}
     2241
     2242    ~MaxWeightedBpMatching4() {
     2243      destroyStructures();
     2244    }
     2245
     2246    /// \name Execution Control
     2247    /// The simplest way to execute the algorithm is to use the
     2248    /// \ref run() member function.
     2249
     2250    ///@{
     2251
     2252    /// \brief Initialize the algorithm
     2253    ///
     2254    /// This function initializes the algorithm.
     2255    void init() {
     2256      createStructures();
     2257
     2258      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2259        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     2260        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     2261        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     2262      }
     2263      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     2264        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     2265      }
     2266
     2267      _delta1->clear();
     2268      _delta2->clear();
     2269      _delta3->clear();
     2270      _tree_set->clear();
     2271
     2272      _unmatched = _node_num;
     2273
     2274      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2275        Value max = 0;
     2276        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2277          if ((dualScale * _weight[a]) / 2 > max) {
     2278            max = (dualScale * _weight[a]) / 2;
     2279          }
     2280        }
     2281        _node_potential->set(n, max);
     2282        _delta1->push(n, max);
     2283        _tree_set->insert(n);
     2284
     2285        for (InArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2286          (*_node_heap)[_bpgraph.source(a)]->push(
     2287              a, (*_node_potential)[n] - dualScale * _weight[a]);
     2288        }
     2289
     2290        _matching->set(n, INVALID);
     2291        _status->set(n, EVEN);
     2292      }
     2293      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     2294        _delta3->push(e, ((*_node_potential)[_bpgraph.u(e)] +
     2295                          (*_node_potential)[_bpgraph.v(e)] -
     2296                          dualScale * _weight[e]) / 2);
     2297      }
     2298    }
     2299
     2300    /// \brief Initialize the algorithm
     2301    ///
     2302    /// This function initializes the algorithm.
     2303    void redRootInit() {
     2304      createStructures();
     2305
     2306      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2307        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     2308        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     2309        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     2310      }
     2311      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     2312        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     2313      }
     2314
     2315      _delta1->clear();
     2316      _delta2->clear();
     2317      _delta3->clear();
     2318      _tree_set->clear();
     2319
     2320      _unmatched = 0;
     2321      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     2322        Value max = 0;
     2323        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2324          if ((dualScale * _weight[a]) > max) {
     2325            max = dualScale * _weight[a];
     2326          }
     2327        }
     2328        _node_potential->set(n, max);
     2329        _delta1->push(n, max);
     2330
     2331        _tree_set->insert(n);
     2332
     2333        _matching->set(n, INVALID);
     2334        _status->set(n, EVEN);
     2335
     2336        ++_unmatched;
     2337      }
     2338      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     2339        _matching->set(n, INVALID);
     2340        _status->set(n, MATCHED);
     2341        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2342          Node v = _bpgraph.target(a);
     2343          Value rw = (*_node_potential)[v] - dualScale * _weight[a];
     2344          (*_node_heap)[n]->push(a, rw);
     2345        }
     2346        if (!(*_node_heap)[n]->empty()) {
     2347          _delta2->push(n, (*_node_heap)[n]->prio() + (*_node_potential)[n]);
     2348        }
     2349      }
     2350    }
     2351
     2352    /// \brief Initialize the algorithm
     2353    ///
     2354    /// This function initializes the algorithm.
     2355    void blueRootInit() {
     2356      createStructures();
     2357
     2358      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2359        (*_delta1_index)[n] = _delta1->PRE_HEAP;
     2360        (*_delta2_index)[n] = _delta1->PRE_HEAP;
     2361        (*_node_heap)[n] = new BinHeap<Value, IntArcMap>(*_node_heap_index);
     2362      }
     2363      for (EdgeIt e(_bpgraph); e != INVALID; ++e) {
     2364        (*_delta3_index)[e] = _delta3->PRE_HEAP;
     2365      }
     2366
     2367      _delta1->clear();
     2368      _delta2->clear();
     2369      _delta3->clear();
     2370      _tree_set->clear();
     2371
     2372      _unmatched = 0;
     2373      for (BlueNodeIt n(_bpgraph); n != INVALID; ++n) {
     2374        Value max = 0;
     2375        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2376          if ((dualScale * _weight[a]) > max) {
     2377            max = dualScale * _weight[a];
     2378          }
     2379        }
     2380        _node_potential->set(n, max);
     2381        _delta1->push(n, max);
     2382
     2383        _tree_set->insert(n);
     2384
     2385        _matching->set(n, INVALID);
     2386        _status->set(n, EVEN);
     2387
     2388        ++_unmatched;
     2389      }
     2390      for (RedNodeIt n(_bpgraph); n != INVALID; ++n) {
     2391        _matching->set(n, INVALID);
     2392        _status->set(n, MATCHED);
     2393        for (OutArcIt a(_bpgraph, n); a != INVALID; ++a) {
     2394          Node v = _bpgraph.target(a);
     2395          Value rw = (*_node_potential)[v] - dualScale * _weight[a];
     2396          (*_node_heap)[n]->push(a, rw);
     2397        }
     2398        if (!(*_node_heap)[n]->empty()) {
     2399          _delta2->push(n, (*_node_heap)[n]->prio() + (*_node_potential)[n]);
     2400        }
     2401      }
     2402    }
     2403
     2404    /// \brief Start the algorithm
     2405    ///
     2406    /// This function starts the algorithm.
     2407    ///
     2408    /// \pre \ref init() must be called before using this function.
     2409    void start() {
     2410      enum OpType {
     2411        D1, D2, D3
     2412      };
     2413
     2414      while (_unmatched > 0) {
     2415        Value d1 = !_delta1->empty() ?
     2416          _delta1->prio() : std::numeric_limits<Value>::max();
     2417
     2418        Value d2 = !_delta2->empty() ?
     2419          _delta2->prio() : std::numeric_limits<Value>::max();
     2420
     2421        Value d3 = !_delta3->empty() ?
     2422          _delta3->prio() : std::numeric_limits<Value>::max();
     2423
     2424        _delta_sum = d3; OpType ot = D3;
     2425        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
     2426        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
     2427
     2428        switch (ot) {
     2429        case D1:
     2430          {
     2431            Node n = _delta1->top();
     2432            unmatchNode(n);
     2433            --_unmatched;
     2434          }
     2435          break;
     2436        case D2:
     2437          {
     2438            Node n = _delta2->top();
     2439            Arc a = (*_node_heap)[n]->top();
     2440            if ((*_matching)[n] == INVALID) {
     2441              augmentOnArc(a);
     2442              --_unmatched;
     2443            } else {
     2444              extendOnArc(a);
     2445            }
     2446          }
     2447          break;
     2448        case D3:
     2449          {
     2450            Edge e = _delta3->top();
     2451            augmentOnEdge(e);
     2452            _unmatched -= 2;
     2453          }
     2454          break;
     2455        }
     2456      }
     2457    }
     2458
     2459    /// \brief Run the algorithm.
     2460    ///
     2461    /// This method runs the \c %MaxWeightedBpMatching4 algorithm.
     2462    ///
     2463    /// \note mwbpm.run() is just a shortcut of the following code.
     2464    /// \code
     2465    ///   mwbpm.init();
     2466    ///   mwbpm.start();
     2467    /// \endcode
     2468    void run() {
     2469      init();
     2470      start();
     2471    }
     2472
     2473    /// @}
     2474
     2475    /// \name Primal Solution
     2476    /// Functions to get the primal solution, i.e. the maximum weighted
     2477    /// bipartite matching.\n
     2478    /// Either \ref run() or \ref start() function should be called before
     2479    /// using them.
     2480
     2481    /// @{
     2482
     2483    /// \brief Return the weight of the matching.
     2484    ///
     2485    /// This function returns the weight of the found matching.
     2486    ///
     2487    /// \pre Either run() or start() must be called before using this function.
     2488    Value matchingWeight() const {
     2489      Value sum = 0;
     2490      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2491        if ((*_matching)[n] != INVALID) {
     2492          sum += _weight[(*_matching)[n]];
     2493        }
     2494      }
     2495      return sum / 2;
     2496    }
     2497
     2498    /// \brief Return the size (cardinality) of the matching.
     2499    ///
     2500    /// This function returns the size (cardinality) of the found matching.
     2501    ///
     2502    /// \pre Either run() or start() must be called before using this function.
     2503    int matchingSize() const {
     2504      int num = 0;
     2505      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2506        if ((*_matching)[n] != INVALID) {
     2507          ++num;
     2508        }
     2509      }
     2510      return num / 2;
     2511    }
     2512
     2513    /// \brief Return \c true if the given edge is in the matching.
     2514    ///
     2515    /// This function returns \c true if the given edge is in the found
     2516    /// matching.
     2517    ///
     2518    /// \pre Either run() or start() must be called before using this function.
     2519    bool matching(const Edge& edge) const {
     2520      return edge == (*_matching)[_bpgraph.u(edge)];
     2521    }
     2522
     2523    /// \brief Return the matching arc (or edge) incident to the given node.
     2524    ///
     2525    /// This function returns the matching arc (or edge) incident to the
     2526    /// given node in the found matching or \c INVALID if the node is
     2527    /// not covered by the matching.
     2528    ///
     2529    /// \pre Either run() or start() must be called before using this function.
     2530    Arc matching(const Node& node) const {
     2531      return (*_matching)[node];
     2532    }
     2533
     2534    /// \brief Return the mate of the given node.
     2535    ///
     2536    /// This function returns the mate of the given node in the found
     2537    /// matching or \c INVALID if the node is not covered by the matching.
     2538    ///
     2539    /// \pre Either run() or start() must be called before using this function.
     2540    Node mate(const Node& node) const {
     2541      return (*_matching)[node] != INVALID ?
     2542        _bpgraph.target((*_matching)[node]) : INVALID;
     2543    }
     2544
     2545    /// @}
     2546
     2547    /// \name Dual Solution
     2548    /// Functions to get the dual solution.\n
     2549    /// Either \ref run() or \ref start() function should be called before
     2550    /// using them.
     2551
     2552    /// @{
     2553
     2554    /// \brief Return the value of the dual solution.
     2555    ///
     2556    /// This function returns the value of the dual solution.
     2557    /// It should be equal to the primal value scaled by \ref dualScale
     2558    /// "dual scale".
     2559    ///
     2560    /// \pre Either run() or start() must be called before using this function.
     2561    Value dualValue() const {
     2562      Value sum = 0;
     2563      for (NodeIt n(_bpgraph); n != INVALID; ++n) {
     2564        sum += nodeValue(n);
     2565      }
     2566      return sum;
     2567    }
     2568
     2569    /// \brief Return the dual value (potential) of the given node.
     2570    ///
     2571    /// This function returns the dual value (potential) of the given node.
     2572    ///
     2573    /// \pre Either run() or start() must be called before using this function.
     2574    Value nodeValue(const Node& n) const {
     2575      return (*_node_potential)[n];
     2576    }
     2577
     2578    /// @}
     2579
     2580  };
     2581
     2582} //END OF NAMESPACE LEMON
     2583
     2584#endif //LEMON_BPMATCHING_H
  • lemon/matching.h

    diff -r 4fd76139b69e -r 3180fd18651b lemon/matching.h
    a b  
    11831183        (*_blossom_data)[even].next =
    11841184          _graph.oppositeArc((*_blossom_data)[odd].pred);
    11851185      }
    1186 
    11871186    }
    11881187
    11891188    void destroyTree(int tree) {
     
    15861585        _delta_sum(), _unmatched(0),
    15871586
    15881587        _fractional(0)
     1588
    15891589    {}
    15901590
    15911591    ~MaxWeightedMatching() {
  • test/CMakeLists.txt

    diff -r 4fd76139b69e -r 3180fd18651b test/CMakeLists.txt
    a b  
    1717  bellman_ford_test
    1818  bfs_test
    1919  bpgraph_test
     20  bpmatching_test
    2021  circulation_test
    2122  connectivity_test
    2223  counter_test
  • new file test/bpmatching_test.cc

    diff -r 4fd76139b69e -r 3180fd18651b test/bpmatching_test.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2012
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <list>
     22
     23#include <lemon/random.h>
     24#include <lemon/bpmatching.h>
     25#include <lemon/matching.h>
     26#include <lemon/smart_graph.h>
     27#include <lemon/concepts/bpgraph.h>
     28#include <lemon/concepts/maps.h>
     29#include <lemon/lgf_reader.h>
     30#include <lemon/time_measure.h>
     31#include <lemon/network_simplex.h>
     32#include <lemon/cost_scaling.h>
     33#include <lemon/capacity_scaling.h>
     34
     35#include "test_tools.h"
     36
     37using namespace std;
     38using namespace lemon;
     39
     40const int lgfn = 6;
     41const string lgf[lgfn] = {
     42  // Empty graph
     43  "@red_nodes\n"
     44  "label\n"
     45  "\n"
     46  "@blue_nodes\n"
     47  "label\n"
     48  "\n"
     49  "@edges\n"
     50  "  label weight\n"
     51  "\n",
     52
     53  // No red nodes
     54  "@red_nodes\n"
     55  "label\n"
     56  "\n"
     57  "@blue_nodes\n"
     58  "label\n"
     59  "1\n"
     60  "@edges\n"
     61  "  label weight\n"
     62  "\n",
     63
     64  // No blue nodes
     65  "@red_nodes\n"
     66  "label\n"
     67  "1\n"
     68  "\n"
     69  "@blue_nodes\n"
     70  "label\n"
     71  "\n"
     72  "@edges\n"
     73  "  label weight\n"
     74  "\n",
     75
     76  "@red_nodes\n"
     77  "label\n"
     78  "1\n"
     79  "@blue_nodes\n"
     80  "label\n"
     81  "2\n"
     82  "\n"
     83  "@edges\n"
     84  "  label weight\n"
     85  "1 2 1 1\n"
     86  "\n",
     87
     88  "@red_nodes\n"
     89  "label\n"
     90  "1\n"
     91  "2\n"
     92  "\n"
     93  "@blue_nodes\n"
     94  "label\n"
     95  "3\n"
     96  "4\n"
     97  "\n"
     98  "@edges\n"
     99  "  label weight\n"
     100  "1 3 1 5\n"
     101  "1 4 2 3\n"
     102  "2 3 3 1\n"
     103  "2 4 4 2\n"
     104  "\n",
     105
     106  "@red_nodes\n"
     107  "label\n"
     108  "1\n"
     109  "2\n"
     110  "\n"
     111  "@blue_nodes\n"
     112  "label\n"
     113  "3\n"
     114  "4\n"
     115  "\n"
     116  "@edges\n"
     117  "  label weight\n"
     118  "1 3 1 3\n"
     119  "1 4 2 2\n"
     120  "2 3 3 2\n"
     121  "\n",
     122};
     123
     124void checkMaxWeightedBpMatchingCompiles() {
     125  typedef concepts::BpGraph BpGraph;
     126  typedef BpGraph::Node Node;
     127  typedef BpGraph::Edge Edge;
     128  typedef BpGraph::EdgeMap<int> IntEdgeMap;
     129  typedef MaxWeightedBpMatching1<BpGraph, IntEdgeMap> MWBPM;
     130
     131  int v;
     132  Node n;
     133  Edge e;
     134  bool b;
     135
     136  BpGraph graph;
     137  IntEdgeMap weight(graph);
     138  MWBPM mwbpm(graph, weight);
     139  const MWBPM& cmwbpm = mwbpm;
     140
     141  mwbpm.init();
     142  mwbpm.start();
     143  mwbpm.run();
     144
     145  v = cmwbpm.matchingWeight();
     146  v = cmwbpm.matchingSize();
     147  e = cmwbpm.matching(n);
     148  b = cmwbpm.matching(e);
     149  n = cmwbpm.mate(n);
     150  v = cmwbpm.dualValue();
     151  v = cmwbpm.nodeValue(n);
     152
     153  ::lemon::ignore_unused_variable_warning(v, b);
     154}
     155
     156
     157template <typename BpMatching>
     158void checkMaxWeightedBpMatching(const SmartBpGraph& bpgraph,
     159                                const SmartBpGraph::EdgeMap<int>& weight,
     160                                const BpMatching& mwbpm) {
     161  int pv = 0;
     162  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
     163    int rw = mwbpm.nodeValue(bpgraph.redNode(e)) +
     164             mwbpm.nodeValue(bpgraph.blueNode(e));
     165    rw -= weight[e] * mwbpm.dualScale;
     166    check(rw >= 0, "Negative reduced weight");
     167    check(rw == 0 || !mwbpm.matching(e),
     168          "Non-zero reduced weight on matching edge");
     169    if (mwbpm.matching(e)) {
     170      pv += weight[e];
     171    }
     172  }
     173  int dv = 0;
     174  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     175    if (mwbpm.matching(n) != INVALID) {
     176      check(mwbpm.nodeValue(n) >= 0, "Invalid node value");
     177      SmartBpGraph::Node o = bpgraph.target(mwbpm.matching(n));
     178      check(mwbpm.mate(n) == o, "Invalid matching");
     179      check(mwbpm.matching(n) == bpgraph.oppositeArc(mwbpm.matching(o)),
     180            "Invalid matching");
     181    } else {
     182      check(mwbpm.mate(n) == INVALID, "Invalid matching");
     183      check(mwbpm.nodeValue(n) == 0, "Invalid matching");
     184    }
     185  }
     186
     187  for (SmartBpGraph::NodeIt n(bpgraph); n != INVALID; ++n) {
     188    dv += mwbpm.nodeValue(n);
     189  }
     190
     191  check(pv * mwbpm.dualScale == dv, "Wrong duality");
     192}
     193
     194void generateRandomBpGraph(int redNum, int blueNum, int edgeNum, int min, int max,
     195                           SmartBpGraph& bpgraph, SmartBpGraph::EdgeMap<int>& weight) {
     196  std::vector<SmartBpGraph::RedNode> redNodes;
     197  for (int i = 0; i < redNum; ++i) {
     198    redNodes.push_back(bpgraph.addRedNode());
     199  }
     200  std::vector<SmartBpGraph::BlueNode> blueNodes;
     201  for (int i = 0; i < blueNum; ++i) {
     202    blueNodes.push_back(bpgraph.addBlueNode());
     203  }
     204  for (int i = 0; i < edgeNum; ++i) {
     205    SmartBpGraph::Edge e = bpgraph.addEdge(
     206        redNodes[lemon::rnd[redNum]], blueNodes[lemon::rnd[blueNum]]);
     207    weight[e] = lemon::rnd[max - min] + min;
     208  }
     209}
     210
     211void generateFullBpGraph(int redNum, int blueNum, int min, int max,
     212                         SmartBpGraph& bpgraph, SmartBpGraph::EdgeMap<int>& weight) {
     213  std::vector<SmartBpGraph::RedNode> redNodes;
     214  for (int i = 0; i < redNum; ++i) {
     215    redNodes.push_back(bpgraph.addRedNode());
     216  }
     217  std::vector<SmartBpGraph::BlueNode> blueNodes;
     218  for (int i = 0; i < blueNum; ++i) {
     219    blueNodes.push_back(bpgraph.addBlueNode());
     220  }
     221  for (int i = 0; i < redNum; ++i) {
     222    for (int j = 0; j < blueNum; ++j) {
     223      SmartBpGraph::Edge e = bpgraph.addEdge(
     224          redNodes[i], blueNodes[j]);
     225      weight[e] = lemon::rnd[max - min] + min;
     226    }
     227  }
     228}
     229
     230void generateSkewBpGraph(int partNodeNum, int deg,
     231                         SmartBpGraph& bpgraph,
     232                         SmartBpGraph::EdgeMap<int>& weight) {
     233  std::vector<SmartBpGraph::RedNode> redNodes;
     234  for (int i = 0; i < partNodeNum; ++i) {
     235    redNodes.push_back(bpgraph.addRedNode());
     236  }
     237  for (int i = 0; i < partNodeNum; ++i) {
     238    std::swap(redNodes[i], redNodes[lemon::rnd[partNodeNum - i] + i]);
     239  }
     240  std::vector<SmartBpGraph::BlueNode> blueNodes;
     241  for (int i = 0; i < partNodeNum; ++i) {
     242    blueNodes.push_back(bpgraph.addBlueNode());
     243  }
     244  for (int i = 0; i < partNodeNum; ++i) {
     245    std::swap(blueNodes[i], blueNodes[lemon::rnd[partNodeNum - i] + i]);
     246  }
     247  for (int d = 0; d < deg; ++d) {
     248    for (int i = d; i < partNodeNum; ++i) {
     249      SmartBpGraph::Edge e = bpgraph.addEdge(
     250          redNodes[i], blueNodes[i - d]);
     251      weight[e] = partNodeNum + d;
     252    }
     253  }
     254}
     255
     256void convertToMinCostFlow(const SmartBpGraph& bpgraph,
     257                          const SmartBpGraph::EdgeMap<int>& weight,
     258                          SmartDigraph& digraph,
     259                          SmartDigraph::NodeMap<int>& supply,
     260                          SmartDigraph::ArcMap<int>& cost) {
     261  SmartDigraph::Node outer = digraph.addNode();
     262  supply.set(outer, countBlueNodes(bpgraph) - countRedNodes(bpgraph));
     263
     264  SmartBpGraph::NodeMap<SmartDigraph::Node> nodeRef(bpgraph);
     265  for (SmartBpGraph::RedNodeIt n(bpgraph); n != INVALID; ++n) {
     266    nodeRef[n] = digraph.addNode();
     267    supply.set(nodeRef[n], 1);
     268    cost.set(digraph.addArc(nodeRef[n], outer), 0);
     269  }
     270  for (SmartBpGraph::BlueNodeIt n(bpgraph); n != INVALID; ++n) {
     271    nodeRef[n] = digraph.addNode();
     272    supply.set(nodeRef[n], -1);
     273    cost.set(digraph.addArc(outer, nodeRef[n]), 0);
     274  }
     275
     276  for (SmartBpGraph::EdgeIt e(bpgraph); e != INVALID; ++e) {
     277    cost.set(
     278        digraph.addArc(nodeRef[bpgraph.redNode(e)], nodeRef[bpgraph.blueNode(e)]),
     279        - weight[e]);
     280  }
     281}
     282
     283
     284void runAll(const SmartBpGraph& bpgraph,
     285            const SmartBpGraph::EdgeMap<int>& weight,
     286            bool includeMcf=true) {
     287  SmartDigraph mcfDigraph;
     288  SmartDigraph::NodeMap<int> mcfSupply(mcfDigraph);
     289  SmartDigraph::ArcMap<int> mcfCost(mcfDigraph);
     290  convertToMinCostFlow(bpgraph, weight, mcfDigraph, mcfSupply, mcfCost);
     291
     292  {
     293    Timer t;
     294    MaxWeightedMatching<SmartBpGraph> mwm(bpgraph, weight);
     295    mwm.init();
     296    mwm.start();
     297    std::cerr << "MaxWeightedMatching/init: " << t << std::endl;
     298  }
     299  {
     300    Timer t;
     301    MaxWeightedMatching<SmartBpGraph> mwm(bpgraph, weight);
     302    mwm.fractionalInit();
     303    mwm.start();
     304    std::cerr << "MaxWeightedMatching/fractionalInit: " << t << std::endl;
     305  }
     306  {
     307    Timer t;
     308    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
     309    mwbpm.init();
     310    mwbpm.start();
     311    std::cerr << "MaxWeightedBpMatching1/init: " << t << std::endl;
     312    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     313  }
     314  {
     315    Timer t;
     316    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
     317    mwbpm.redRootInit();
     318    mwbpm.start();
     319    std::cerr << "MaxWeightedBpMatching1/redRootInit: " << t << std::endl;
     320    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     321  }
     322  {
     323    Timer t;
     324    MaxWeightedBpMatching1<SmartBpGraph> mwbpm(bpgraph, weight);
     325    mwbpm.blueRootInit();
     326    mwbpm.start();
     327    std::cerr << "MaxWeightedBpMatching1/blueRootInit: " << t << std::endl;
     328    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     329  }
     330  {
     331    Timer t;
     332    MaxWeightedBpMatching2<SmartBpGraph> mwbpm(bpgraph, weight);
     333    mwbpm.redRootInit();
     334    mwbpm.start();
     335    std::cerr << "MaxWeightedBpMatching2/redRootInit: " << t << std::endl;
     336    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     337  }
     338  {
     339    Timer t;
     340    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
     341    mwbpm.init();
     342    mwbpm.start();
     343    std::cerr << "MaxWeightedBpMatching3/init: " << t << std::endl;
     344    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     345  }
     346  {
     347    Timer t;
     348    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
     349    mwbpm.redRootInit();
     350    mwbpm.start();
     351    std::cerr << "MaxWeightedBpMatching3/redRootInit: " << t << std::endl;
     352    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     353  }
     354  {
     355    Timer t;
     356    MaxWeightedBpMatching3<SmartBpGraph> mwbpm(bpgraph, weight);
     357    mwbpm.blueRootInit();
     358    mwbpm.start();
     359    std::cerr << "MaxWeightedBpMatching3/blueRootInit: " << t << std::endl;
     360    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     361  }
     362  {
     363    Timer t;
     364    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
     365    mwbpm.init();
     366    mwbpm.start();
     367    std::cerr << "MaxWeightedBpMatching4/init: " << t << std::endl;
     368    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     369  }
     370  {
     371    Timer t;
     372    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
     373    mwbpm.redRootInit();
     374    mwbpm.start();
     375    std::cerr << "MaxWeightedBpMatching4/redRootInit: " << t << std::endl;
     376    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     377  }
     378  {
     379    Timer t;
     380    MaxWeightedBpMatching4<SmartBpGraph> mwbpm(bpgraph, weight);
     381    mwbpm.blueRootInit();
     382    mwbpm.start();
     383    std::cerr << "MaxWeightedBpMatching4/blueRootInit: " << t << std::endl;
     384    checkMaxWeightedBpMatching(bpgraph, weight, mwbpm);
     385  }
     386  if (includeMcf) {
     387    Timer t;
     388    NetworkSimplex<SmartDigraph> ns(mcfDigraph);
     389    ns.costMap(mcfCost);
     390    ns.supplyMap(mcfSupply);
     391    ns.upperMap(constMap<SmartDigraph::Arc, int>(1));
     392    ns.run();
     393    std::cerr << "NetworkSimplex: " << t << std::endl;
     394  }
     395  if (includeMcf) {
     396    Timer t;
     397    NetworkSimplex<SmartDigraph> ns(mcfDigraph);
     398    ns.costMap(mcfCost);
     399    ns.supplyMap(mcfSupply);
     400    ns.upperMap(constMap<SmartDigraph::Arc, int>(1));
     401    ns.run(NetworkSimplex<SmartDigraph>::ALTERING_LIST);
     402    std::cerr << "NetworkSimplex/ALTERING_LIST: " << t << std::endl;
     403  }
     404  if (includeMcf) {
     405    Timer t;
     406    CapacityScaling<SmartDigraph> cs(mcfDigraph);
     407    cs.costMap(mcfCost);
     408    cs.supplyMap(mcfSupply);
     409    cs.upperMap(constMap<SmartDigraph::Arc, int>(1));
     410    cs.run();
     411    std::cerr << "CapacityScaling: " << t << std::endl;
     412  }
     413  if (includeMcf) {
     414    Timer t;
     415    CostScaling<SmartDigraph> cs(mcfDigraph);
     416    cs.costMap(mcfCost);
     417    cs.supplyMap(mcfSupply);
     418    cs.upperMap(constMap<SmartDigraph::Arc, int>(1));
     419    cs.run();
     420    std::cerr << "CostScaling: " << t << std::endl;
     421  }
     422}
     423
     424
     425int main() {
     426  // // Checking hard wired extremal graphs
     427  // for (int i=0; i<lgfn; ++i) {
     428  //   SmartBpGraph bpgraph;
     429  //   SmartBpGraph::EdgeMap<int> weight(bpgraph);
     430  //   istringstream lgfs(lgf[i]);
     431  //   bpGraphReader(bpgraph, lgfs)
     432  //       .edgeMap("weight", weight).run();
     433  //   runAll(bpgraph, weight);
     434  // }
     435
     436  std::cerr << "Random red=10000 blue=10000 edge=10000 range=0-1000" << std::endl;
     437  std::cerr << "=================================================================" << std::endl;
     438  for (int i = 0; i < 3; ++i) {
     439    std::cerr << "===== case " << i << ":" << std::endl;
     440    SmartBpGraph bpgraph;
     441    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     442    generateRandomBpGraph(10000, 10000, 10000, 0, 1000, bpgraph, weight);
     443    runAll(bpgraph, weight);
     444  }
     445  std::cerr << std::endl;
     446
     447  std::cerr << "Random red=20000 blue=20000 edge=20000 range=0-1000" << std::endl;
     448  std::cerr << "=================================================================" << std::endl;
     449  for (int i = 0; i < 3; ++i) {
     450    std::cerr << "===== case " << i << ":" << std::endl;
     451    SmartBpGraph bpgraph;
     452    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     453    generateRandomBpGraph(20000, 20000, 20000, 0, 1000, bpgraph, weight);
     454    runAll(bpgraph, weight);
     455  }
     456  std::cerr << std::endl;
     457
     458  std::cerr << "Random red=40000 blue=40000 edge=40000 range=0-1000" << std::endl;
     459  std::cerr << "=================================================================" << std::endl;
     460  for (int i = 0; i < 3; ++i) {
     461    std::cerr << "===== case " << i << ":" << std::endl;
     462    SmartBpGraph bpgraph;
     463    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     464    generateRandomBpGraph(40000, 40000, 40000, 0, 1000, bpgraph, weight);
     465    runAll(bpgraph, weight);
     466  }
     467  std::cerr << std::endl;
     468
     469  std::cerr << "Random red=100000 blue=100000 edge=100000 range=0-1000" << std::endl;
     470  std::cerr << "=================================================================" << std::endl;
     471  for (int i = 0; i < 3; ++i) {
     472    std::cerr << "===== case " << i << ":" << std::endl;
     473    SmartBpGraph bpgraph;
     474    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     475    generateRandomBpGraph(100000, 100000, 100000, 0, 1000, bpgraph, weight);
     476    runAll(bpgraph, weight, false);
     477  }
     478  std::cerr << std::endl;
     479
     480  std::cerr << "Random red=1000000 blue=1000000 edge=1000000 range=0-1000" << std::endl;
     481  for (int i = 0; i < 3; ++i) {
     482    std::cerr << "===== case " << i << ":" << std::endl;
     483    SmartBpGraph bpgraph;
     484    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     485    generateRandomBpGraph(1000000, 1000000, 1000000, 0, 1000, bpgraph, weight);
     486    runAll(bpgraph, weight, false);
     487  }
     488  std::cerr << std::endl;
     489
     490  std::cerr << "Random red=10000 blue=10000 edge=10000 range=1000-1100" << std::endl;
     491  std::cerr << "=================================================================" << std::endl;
     492  for (int i = 0; i < 3; ++i) {
     493    std::cerr << "===== case " << i << ":" << std::endl;
     494    SmartBpGraph bpgraph;
     495    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     496    generateRandomBpGraph(10000, 10000, 10000, 1000, 1100, bpgraph, weight);
     497    runAll(bpgraph, weight);
     498  }
     499  std::cerr << std::endl;
     500
     501  std::cerr << "Random red=20000 blue=20000 edge=20000 range=1000-1100" << std::endl;
     502  std::cerr << "=================================================================" << std::endl;
     503  for (int i = 0; i < 3; ++i) {
     504    std::cerr << "===== case " << i << ":" << std::endl;
     505    SmartBpGraph bpgraph;
     506    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     507    generateRandomBpGraph(20000, 20000, 20000, 1000, 1100, bpgraph, weight);
     508    runAll(bpgraph, weight);
     509  }
     510  std::cerr << std::endl;
     511
     512  std::cerr << "Random red=40000 blue=40000 edge=40000 range=1000-1100" << std::endl;
     513  std::cerr << "=================================================================" << std::endl;
     514  for (int i = 0; i < 3; ++i) {
     515    std::cerr << "===== case " << i << ":" << std::endl;
     516    SmartBpGraph bpgraph;
     517    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     518    generateRandomBpGraph(40000, 40000, 40000, 1000, 1100, bpgraph, weight);
     519    runAll(bpgraph, weight);
     520  }
     521  std::cerr << std::endl;
     522
     523  std::cerr << "Random red=100000 blue=100000 edge=100000 range=1000-1100" << std::endl;
     524  std::cerr << "=================================================================" << std::endl;
     525  for (int i = 0; i < 3; ++i) {
     526    std::cerr << "===== case " << i << ":" << std::endl;
     527    SmartBpGraph bpgraph;
     528    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     529    generateRandomBpGraph(100000, 100000, 100000, 1000, 1100, bpgraph, weight);
     530    runAll(bpgraph, weight, false);
     531  }
     532  std::cerr << std::endl;
     533
     534  std::cerr << "Random red=1000000 blue=1000000 edge=1000000 range=1000-1100" << std::endl;
     535  std::cerr << "=================================================================" << std::endl;
     536  for (int i = 0; i < 3; ++i) {
     537    std::cerr << "===== case " << i << ":" << std::endl;
     538    SmartBpGraph bpgraph;
     539    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     540    generateRandomBpGraph(1000000, 1000000, 1000000, 1000, 1100, bpgraph, weight);
     541    runAll(bpgraph, weight, false);
     542  }
     543  std::cerr << std::endl;
     544
     545  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
     546  std::cerr << "=================================================================" << std::endl;
     547  for (int i = 0; i < 3; ++i) {
     548    std::cerr << "===== case " << i << ":" << std::endl;
     549    SmartBpGraph bpgraph;
     550    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     551    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
     552    runAll(bpgraph, weight, false);
     553  }
     554  std::cerr << std::endl;
     555
     556  std::cerr << "Random red=10000 blue=10000 edge=20000 range=0-1000" << std::endl;
     557  std::cerr << "=================================================================" << std::endl;
     558  for (int i = 0; i < 3; ++i) {
     559    std::cerr << "===== case " << i << ":" << std::endl;
     560    SmartBpGraph bpgraph;
     561    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     562    generateRandomBpGraph(10000, 10000, 20000, 0, 1000, bpgraph, weight);
     563    runAll(bpgraph, weight);
     564  }
     565  std::cerr << std::endl;
     566
     567  std::cerr << "Random red=20000 blue=20000 edge=40000 range=0-1000" << std::endl;
     568  std::cerr << "=================================================================" << std::endl;
     569  for (int i = 0; i < 3; ++i) {
     570    std::cerr << "===== case " << i << ":" << std::endl;
     571    SmartBpGraph bpgraph;
     572    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     573    generateRandomBpGraph(20000, 20000, 40000, 0, 1000, bpgraph, weight);
     574    runAll(bpgraph, weight);
     575  }
     576  std::cerr << std::endl;
     577
     578  std::cerr << "Random red=100000 blue=100000 edge=200000 range=0-1000" << std::endl;
     579  std::cerr << "=================================================================" << std::endl;
     580  for (int i = 0; i < 3; ++i) {
     581    std::cerr << "===== case " << i << ":" << std::endl;
     582    SmartBpGraph bpgraph;
     583    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     584    generateRandomBpGraph(100000, 100000, 200000, 0, 1000, bpgraph, weight);
     585    runAll(bpgraph, weight, false);
     586  }
     587  std::cerr << std::endl;
     588
     589  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=0-1000" << std::endl;
     590  std::cerr << "=================================================================" << std::endl;
     591  for (int i = 0; i < 3; ++i) {
     592    std::cerr << "===== case " << i << ":" << std::endl;
     593    SmartBpGraph bpgraph;
     594    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     595    generateRandomBpGraph(1000000, 1000000, 2000000, 0, 1000, bpgraph, weight);
     596    runAll(bpgraph, weight, false);
     597  }
     598  std::cerr << std::endl;
     599
     600  std::cerr << "Random red=10000 blue=10000 edge=20000 range=1000-1100" << std::endl;
     601  std::cerr << "=================================================================" << std::endl;
     602  for (int i = 0; i < 3; ++i) {
     603    std::cerr << "===== case " << i << ":" << std::endl;
     604    SmartBpGraph bpgraph;
     605    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     606    generateRandomBpGraph(10000, 10000, 20000, 1000, 1100, bpgraph, weight);
     607    runAll(bpgraph, weight);
     608  }
     609  std::cerr << std::endl;
     610
     611  std::cerr << "Random red=20000 blue=20000 edge=40000 range=1000-1100" << std::endl;
     612  std::cerr << "=================================================================" << std::endl;
     613  for (int i = 0; i < 3; ++i) {
     614    std::cerr << "===== case " << i << ":" << std::endl;
     615    SmartBpGraph bpgraph;
     616    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     617    generateRandomBpGraph(20000, 20000, 40000, 1000, 1100, bpgraph, weight);
     618    runAll(bpgraph, weight);
     619  }
     620  std::cerr << std::endl;
     621
     622  std::cerr << "Random red=100000 blue=100000 edge=200000 range=1000-1100" << std::endl;
     623  std::cerr << "=================================================================" << std::endl;
     624  for (int i = 0; i < 3; ++i) {
     625    std::cerr << "===== case " << i << ":" << std::endl;
     626    SmartBpGraph bpgraph;
     627    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     628    generateRandomBpGraph(100000, 100000, 200000, 1000, 1100, bpgraph, weight);
     629    runAll(bpgraph, weight, false);
     630  }
     631  std::cerr << std::endl;
     632
     633  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
     634  std::cerr << "=================================================================" << std::endl;
     635  for (int i = 0; i < 3; ++i) {
     636    std::cerr << "===== case " << i << ":" << std::endl;
     637    SmartBpGraph bpgraph;
     638    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     639    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
     640    runAll(bpgraph, weight, false);
     641  }
     642  std::cerr << std::endl;
     643
     644  std::cerr << "Random red=1000000 blue=1000000 edge=2000000 range=1000-1100" << std::endl;
     645  std::cerr << "=================================================================" << std::endl;
     646  for (int i = 0; i < 3; ++i) {
     647    std::cerr << "===== case " << i << ":" << std::endl;
     648    SmartBpGraph bpgraph;
     649    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     650    generateRandomBpGraph(1000000, 1000000, 2000000, 1000, 1100, bpgraph, weight);
     651    runAll(bpgraph, weight, false);
     652  }
     653  std::cerr << std::endl;
     654
     655  std::cerr << "Random red=10000 blue=10000 edge=200000 range=0-10000" << std::endl;
     656  std::cerr << "=================================================================" << std::endl;
     657  for (int i = 0; i < 3; ++i) {
     658    std::cerr << "===== case " << i << ":" << std::endl;
     659    SmartBpGraph bpgraph;
     660    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     661    generateRandomBpGraph(10000, 10000, 200000, 0, 10000, bpgraph, weight);
     662    runAll(bpgraph, weight);
     663  }
     664  std::cerr << std::endl;
     665
     666  std::cerr << "Random red=20000 blue=20000 edge=400000 range=0-1000" << std::endl;
     667  std::cerr << "=================================================================" << std::endl;
     668  for (int i = 0; i < 3; ++i) {
     669    std::cerr << "===== case " << i << ":" << std::endl;
     670    SmartBpGraph bpgraph;
     671    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     672    generateRandomBpGraph(20000, 20000, 400000, 0, 1000, bpgraph, weight);
     673    runAll(bpgraph, weight);
     674  }
     675  std::cerr << std::endl;
     676
     677  std::cerr << "Random red=40000 blue=40000 edge=800000 range=0-1000" << std::endl;
     678  std::cerr << "=================================================================" << std::endl;
     679  for (int i = 0; i < 3; ++i) {
     680    std::cerr << "===== case " << i << ":" << std::endl;
     681    SmartBpGraph bpgraph;
     682    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     683    generateRandomBpGraph(40000, 40000, 800000, 0, 1000, bpgraph, weight);
     684    runAll(bpgraph, weight);
     685  }
     686  std::cerr << std::endl;
     687
     688  std::cerr << "Random red=100000 blue=100000 edge=2000000 range=0-10000" << std::endl;
     689  std::cerr << "=================================================================" << std::endl;
     690  for (int i = 0; i < 3; ++i) {
     691    std::cerr << "===== case " << i << ":" << std::endl;
     692    SmartBpGraph bpgraph;
     693    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     694    generateRandomBpGraph(100000, 100000, 2000000, 0, 10000, bpgraph, weight);
     695    runAll(bpgraph, weight, false);
     696  }
     697  std::cerr << std::endl;
     698
     699  std::cerr << "Skew red=10000 blue=10000 deg=2" << std::endl;
     700  std::cerr << "=================================================================" << std::endl;
     701  for (int i = 0; i < 3; ++i) {
     702    std::cerr << "===== case " << i << ":" << std::endl;
     703    SmartBpGraph bpgraph;
     704    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     705    generateSkewBpGraph(10000, 2, bpgraph, weight);
     706    runAll(bpgraph, weight);
     707  }
     708  std::cerr << std::endl;
     709
     710  std::cerr << "Skew red=20000 blue=20000 deg=2" << std::endl;
     711  std::cerr << "=================================================================" << std::endl;
     712  for (int i = 0; i < 3; ++i) {
     713    std::cerr << "===== case " << i << ":" << std::endl;
     714    SmartBpGraph bpgraph;
     715    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     716    generateSkewBpGraph(20000, 2, bpgraph, weight);
     717    runAll(bpgraph, weight);
     718  }
     719  std::cerr << std::endl;
     720
     721  std::cerr << "Skew red=100000 blue=100000 deg=2" << std::endl;
     722  std::cerr << "=================================================================" << std::endl;
     723  for (int i = 0; i < 3; ++i) {
     724    std::cerr << "===== case " << i << ":" << std::endl;
     725    SmartBpGraph bpgraph;
     726    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     727    generateSkewBpGraph(100000, 2, bpgraph, weight);
     728    runAll(bpgraph, weight, false);
     729  }
     730  std::cerr << std::endl;
     731
     732  std::cerr << "Skew red=1000000 blue=1000000 deg=2" << std::endl;
     733  std::cerr << "=================================================================" << std::endl;
     734  for (int i = 0; i < 3; ++i) {
     735    std::cerr << "===== case " << i << ":" << std::endl;
     736    SmartBpGraph bpgraph;
     737    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     738    generateSkewBpGraph(1000000, 2, bpgraph, weight);
     739    runAll(bpgraph, weight, false);
     740  }
     741  std::cerr << std::endl;
     742
     743  std::cerr << "Skew red=10000 blue=10000 deg=3" << std::endl;
     744  std::cerr << "=================================================================" << std::endl;
     745  for (int i = 0; i < 3; ++i) {
     746    std::cerr << "===== case " << i << ":" << std::endl;
     747    SmartBpGraph bpgraph;
     748    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     749    generateSkewBpGraph(10000, 3, bpgraph, weight);
     750    runAll(bpgraph, weight);
     751  }
     752  std::cerr << std::endl;
     753
     754  std::cerr << "Skew red=20000 blue=20000 deg=3" << std::endl;
     755  std::cerr << "=================================================================" << std::endl;
     756  for (int i = 0; i < 3; ++i) {
     757    std::cerr << "===== case " << i << ":" << std::endl;
     758    SmartBpGraph bpgraph;
     759    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     760    generateSkewBpGraph(20000, 3, bpgraph, weight);
     761    runAll(bpgraph, weight);
     762  }
     763  std::cerr << std::endl;
     764
     765  std::cerr << "Skew red=100000 blue=100000 deg=3" << std::endl;
     766  std::cerr << "=================================================================" << std::endl;
     767  for (int i = 0; i < 3; ++i) {
     768    std::cerr << "===== case " << i << ":" << std::endl;
     769    SmartBpGraph bpgraph;
     770    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     771    generateSkewBpGraph(100000, 3, bpgraph, weight);
     772    runAll(bpgraph, weight, false);
     773  }
     774  std::cerr << std::endl;
     775
     776  std::cerr << "Skew red=1000000 blue=1000000 deg=3" << std::endl;
     777  std::cerr << "=================================================================" << std::endl;
     778  for (int i = 0; i < 3; ++i) {
     779    std::cerr << "===== case " << i << ":" << std::endl;
     780    SmartBpGraph bpgraph;
     781    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     782    generateSkewBpGraph(1000000, 3, bpgraph, weight);
     783    runAll(bpgraph, weight, false);
     784  }
     785  std::cerr << std::endl;
     786
     787  std::cerr << "Skew red=10000 blue=10000 deg=20" << std::endl;
     788  std::cerr << "=================================================================" << std::endl;
     789  for (int i = 0; i < 3; ++i) {
     790    std::cerr << "===== case " << i << ":" << std::endl;
     791    SmartBpGraph bpgraph;
     792    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     793    generateSkewBpGraph(10000, 20, bpgraph, weight);
     794    runAll(bpgraph, weight);
     795  }
     796  std::cerr << std::endl;
     797
     798  std::cerr << "Skew red=20000 blue=20000 deg=20" << std::endl;
     799  std::cerr << "=================================================================" << std::endl;
     800  for (int i = 0; i < 3; ++i) {
     801    std::cerr << "===== case " << i << ":" << std::endl;
     802    SmartBpGraph bpgraph;
     803    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     804    generateSkewBpGraph(20000, 20, bpgraph, weight);
     805    runAll(bpgraph, weight);
     806  }
     807  std::cerr << std::endl;
     808
     809  std::cerr << "Skew red=100000 blue=100000 deg=20" << std::endl;
     810  std::cerr << "=================================================================" << std::endl;
     811  for (int i = 0; i < 3; ++i) {
     812    std::cerr << "===== case " << i << ":" << std::endl;
     813    SmartBpGraph bpgraph;
     814    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     815    generateSkewBpGraph(100000, 20, bpgraph, weight);
     816    runAll(bpgraph, weight, false);
     817  }
     818  std::cerr << std::endl;
     819
     820  std::cerr << "Skew red=1000 blue=1000 deg=1000" << std::endl;
     821  std::cerr << "=================================================================" << std::endl;
     822  for (int i = 0; i < 3; ++i) {
     823    std::cerr << "===== case " << i << ":" << std::endl;
     824    SmartBpGraph bpgraph;
     825    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     826    generateSkewBpGraph(1000, 1000, bpgraph, weight);
     827    runAll(bpgraph, weight);
     828  }
     829  std::cerr << std::endl;
     830
     831  std::cerr << "Full red=4000 blue=1000 range=10000-11000" << std::endl;
     832  std::cerr << "=================================================================" << std::endl;
     833  for (int i = 0; i < 3; ++i) {
     834    std::cerr << "===== case " << i << ":" << std::endl;
     835    SmartBpGraph bpgraph;
     836    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     837    generateFullBpGraph(4000, 1000, 10000, 11000, bpgraph, weight);
     838    runAll(bpgraph, weight);
     839  }
     840  std::cerr << std::endl;
     841
     842  std::cerr << "Full red=1000 blue=4000 range=10000-11000" << std::endl;
     843  std::cerr << "=================================================================" << std::endl;
     844  for (int i = 0; i < 3; ++i) {
     845    std::cerr << "===== case " << i << ":" << std::endl;
     846    SmartBpGraph bpgraph;
     847    SmartBpGraph::EdgeMap<int> weight(bpgraph);
     848    generateFullBpGraph(1000, 4000, 10000, 11000, bpgraph, weight);
     849    runAll(bpgraph, weight);
     850  }
     851  std::cerr << std::endl;
     852
     853  return 0;
     854}