LEMON/LEMON-official Changeset - r917:61120524af27

Repositories » LEMON » LEMON-official

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r917:61120524af27

« 916:636dad

918:86613a »

r917:61120524af27

Sat, 26 Sep 2009 10:17:31

[Not Reviewed]

0 comments (0 inline)

deba@inf.elte.hu
deba@inf.elte.hu

Fractional matching initialization of weighted matchings (#314)

0 2 0

default

2 files changed with 341 insertions and 20 deletions:

lemon/matching.h

313

test/matching_test.cc

↑ Collapse diff ↑

lemon/matching.h

Show inline comments

@@ -25,12 +25,13 @@
 		#include <limits>
 		#include <lemon/core.h>
 		#include <lemon/unionfind.h>
 		#include <lemon/bin_heap.h>
 		#include <lemon/maps.h>
 		#include <lemon/fractional_matching.h>
 		///\ingroup matching
 		///\file
 		///\brief Maximum matching algorithms in general graphs.
 		namespace lemon {
@@ -794,12 +795,16 @@
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    IntIntMap *_delta4_index;
 		    BinHeap<Value, IntIntMap> *_delta4;
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (!_matching) {
@@ -1556,16 +1561,22 @@
 		        _delta1_index(0), _delta1(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
@@ -1588,12 +1599,14 @@
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = 0;
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
@@ -1622,24 +1635,161 @@
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Initialize the algorithm with fractional matching
 		    ///
 		    /// This function initializes the algorithm with a fractional
 		    /// matching. This initialization is also called jumpstart heuristic.
 		    void fractionalInit() {
 		      createStructures();
 		      if (_fractional == 0) {
 		        _fractional = new FractionalMatching(_graph, _weight, false);
+		      }
 		      _fractional->run();
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = 0;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value pot = _fractional->nodeValue(n);
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = pot;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        (*_blossom_data)[blossom].status = MATCHED;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = _fractional->matching(n);
 		        if (_fractional->matching(n) == INVALID) {
 		          (*_blossom_data)[blossom].base = n;
+		        }
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      typename Graph::template NodeMap<bool> processed(_graph, false);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if (processed[n]) continue;
 		        processed[n] = true;
 		        if (_fractional->matching(n) == INVALID) continue;
 		        int num = 1;
 		        Node v = _graph.target(_fractional->matching(n));
 		        while (n != v) {
 		          processed[v] = true;
 		          v = _graph.target(_fractional->matching(v));
 		          ++num;
+		        }
 		        if (num % 2 == 1) {
 		          std::vector<int> subblossoms(num);
 		          subblossoms[--num] = _blossom_set->find(n);
 		          _delta1->push(n, _fractional->nodeValue(n));
 		          v = _graph.target(_fractional->matching(n));
 		          while (n != v) {
 		            subblossoms[--num] = _blossom_set->find(v);
 		            _delta1->push(v, _fractional->nodeValue(v));
 		            v = _graph.target(_fractional->matching(v));
+		          }
 		          int surface =
 		            _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		          (*_blossom_data)[surface].status = EVEN;
 		          (*_blossom_data)[surface].pred = INVALID;
 		          (*_blossom_data)[surface].next = INVALID;
 		          (*_blossom_data)[surface].pot = 0;
 		          (*_blossom_data)[surface].offset = 0;
 		          _tree_set->insert(surface);
 		          ++_unmatched;
+		        }
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int sb = _blossom_set->find(_graph.u(e));
 		        int ti = (*_node_index)[_graph.v(e)];
 		        int tb = _blossom_set->find(_graph.v(e));
 		        if ((*_blossom_data)[sb].status == EVEN &&
 		            (*_blossom_data)[tb].status == EVEN && sb != tb) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        int nb = _blossom_set->find(n);
 		        if ((*_blossom_data)[nb].status != MATCHED) continue;
 		        int ni = (*_node_index)[n];
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.target(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            int vt = _tree_set->find(vb);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, e);
 		                (*_node_data)[ni].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(e, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+		            }
+		          }
+		        }
 		        if (!(*_node_data)[ni].heap.empty()) {
 		          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		          _delta2->push(nb, _blossom_set->classPrio(nb));
+		        }
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
-		    /// \pre \ref init() must be called before using this function.
+		    /// \pre \ref init() or \ref fractionalInit() must be called
 		    /// before using this function.
 		    void start() {
 		      enum OpType {
 		        D1, D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		      while (_unmatched > 0) {
 		        Value d1 = !_delta1->empty() ?
 		          _delta1->prio() : std::numeric_limits<Value>::max();
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
@@ -1656,23 +1806,23 @@
 		        switch (ot) {
 		        case D1:
+		          {
 		            Node n = _delta1->top();
 		            unmatchNode(n);
-		            --unmatched;
+		            --_unmatched;
+		          }
 		          break;
 		        case D2:
+		          {
 		            int blossom = _delta2->top();
 		            Node n = _blossom_set->classTop(blossom);
 		            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
 		            if ((*_blossom_data)[blossom].next == INVALID) {
 		              augmentOnArc(a);
-		              --unmatched;
+		              --_unmatched;
 		            } else {
 		              extendOnArc(a);
+		            }
+		          }
 		          break;
 		        case D3:
@@ -1689,13 +1839,13 @@
 		              int right_tree = _tree_set->find(right_blossom);
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
@@ -1707,17 +1857,17 @@
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedMatching algorithm.
 		    ///
 		    /// \note mwm.run() is just a shortcut of the following code.
 		    /// \code
-		    ///   mwm.init();
+		    ///   mwm.fractionalInit();
 		    ///   mwm.start();
 		    /// \endcode
 		    void run() {
-		      init();
+		      fractionalInit();
 		      start();
+		    }
 		    /// @}
 		    /// \name Primal Solution
@@ -2071,12 +2221,17 @@
 		    BinHeap<Value, IntEdgeMap> *_delta3;
 		    IntIntMap *_delta4_index;
 		    BinHeap<Value, IntIntMap> *_delta4;
 		    Value _delta_sum;
 		    int _unmatched;
 		    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
 		    FractionalMatching;
 		    FractionalMatching *_fractional;
 		    void createStructures() {
 		      _node_num = countNodes(_graph);
 		      _blossom_num = _node_num * 3 / 2;
 		      if (!_matching) {
@@ -2786,16 +2941,22 @@
 		        _tree_set_index(0), _tree_set(0),
 		        _delta2_index(0), _delta2(0),
 		        _delta3_index(0), _delta3(0),
 		        _delta4_index(0), _delta4(0),
-		        _delta_sum() {}
+		        _delta_sum(), _unmatched(0),
 		        _fractional(0)
 		    {}
 		    ~MaxWeightedPerfectMatching() {
 		      destroyStructures();
 		      if (_fractional) {
 		        delete _fractional;
+		      }
+		    }
 		    /// \name Execution Control
 		    /// The simplest way to execute the algorithm is to use the
 		    /// \ref run() member function.
@@ -2815,12 +2976,14 @@
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = _node_num;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value max = - std::numeric_limits<Value>::max();
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          if (_graph.target(e) == n) continue;
 		          if ((dualScale * _weight[e]) / 2 > max) {
@@ -2848,24 +3011,158 @@
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
+		    }
 		    /// \brief Initialize the algorithm with fractional matching
 		    ///
 		    /// This function initializes the algorithm with a fractional
 		    /// matching. This initialization is also called jumpstart heuristic.
 		    void fractionalInit() {
 		      createStructures();
 		      if (_fractional == 0) {
 		        _fractional = new FractionalMatching(_graph, _weight, false);
+		      }
 		      if (!_fractional->run()) {
 		        _unmatched = -1;
 		        return;
+		      }
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+		      }
 		      for (int i = 0; i < _blossom_num; ++i) {
 		        (*_delta2_index)[i] = _delta2->PRE_HEAP;
 		        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+		      }
 		      _unmatched = 0;
 		      int index = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        Value pot = _fractional->nodeValue(n);
 		        (*_node_index)[n] = index;
 		        (*_node_data)[index].pot = pot;
 		        int blossom =
 		          _blossom_set->insert(n, std::numeric_limits<Value>::max());
 		        (*_blossom_data)[blossom].status = MATCHED;
 		        (*_blossom_data)[blossom].pred = INVALID;
 		        (*_blossom_data)[blossom].next = _fractional->matching(n);
 		        (*_blossom_data)[blossom].pot = 0;
 		        (*_blossom_data)[blossom].offset = 0;
 		        ++index;
+		      }
 		      typename Graph::template NodeMap<bool> processed(_graph, false);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        if (processed[n]) continue;
 		        processed[n] = true;
 		        if (_fractional->matching(n) == INVALID) continue;
 		        int num = 1;
 		        Node v = _graph.target(_fractional->matching(n));
 		        while (n != v) {
 		          processed[v] = true;
 		          v = _graph.target(_fractional->matching(v));
 		          ++num;
+		        }
 		        if (num % 2 == 1) {
 		          std::vector<int> subblossoms(num);
 		          subblossoms[--num] = _blossom_set->find(n);
 		          v = _graph.target(_fractional->matching(n));
 		          while (n != v) {
 		            subblossoms[--num] = _blossom_set->find(v);
 		            v = _graph.target(_fractional->matching(v));
+		          }
 		          int surface =
 		            _blossom_set->join(subblossoms.begin(), subblossoms.end());
 		          (*_blossom_data)[surface].status = EVEN;
 		          (*_blossom_data)[surface].pred = INVALID;
 		          (*_blossom_data)[surface].next = INVALID;
 		          (*_blossom_data)[surface].pot = 0;
 		          (*_blossom_data)[surface].offset = 0;
 		          _tree_set->insert(surface);
 		          ++_unmatched;
+		        }
+		      }
 		      for (EdgeIt e(_graph); e != INVALID; ++e) {
 		        int si = (*_node_index)[_graph.u(e)];
 		        int sb = _blossom_set->find(_graph.u(e));
 		        int ti = (*_node_index)[_graph.v(e)];
 		        int tb = _blossom_set->find(_graph.v(e));
 		        if ((*_blossom_data)[sb].status == EVEN &&
 		            (*_blossom_data)[tb].status == EVEN && sb != tb) {
 		          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
 		                            dualScale * _weight[e]) / 2);
+		        }
+		      }
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        int nb = _blossom_set->find(n);
 		        if ((*_blossom_data)[nb].status != MATCHED) continue;
 		        int ni = (*_node_index)[n];
 		        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
 		          Node v = _graph.target(e);
 		          int vb = _blossom_set->find(v);
 		          int vi = (*_node_index)[v];
 		          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
 		            dualScale * _weight[e];
 		          if ((*_blossom_data)[vb].status == EVEN) {
 		            int vt = _tree_set->find(vb);
 		            typename std::map<int, Arc>::iterator it =
 		              (*_node_data)[ni].heap_index.find(vt);
 		            if (it != (*_node_data)[ni].heap_index.end()) {
 		              if ((*_node_data)[ni].heap[it->second] > rw) {
 		                (*_node_data)[ni].heap.replace(it->second, e);
 		                (*_node_data)[ni].heap.decrease(e, rw);
 		                it->second = e;
+		              }
 		            } else {
 		              (*_node_data)[ni].heap.push(e, rw);
 		              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+		            }
+		          }
+		        }
 		        if (!(*_node_data)[ni].heap.empty()) {
 		          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
 		          _delta2->push(nb, _blossom_set->classPrio(nb));
+		        }
+		      }
+		    }
 		    /// \brief Start the algorithm
 		    ///
 		    /// This function starts the algorithm.
 		    ///
-		    /// \pre \ref init() must be called before using this function.
+		    /// \pre \ref init() or \ref fractionalInit() must be called before
 		    /// using this function.
 		    bool start() {
 		      enum OpType {
 		        D2, D3, D4
 		      };
 		      int unmatched = _node_num;
 		      while (unmatched > 0) {
 		      if (_unmatched == -1) return false;
 		      while (_unmatched > 0) {
 		        Value d2 = !_delta2->empty() ?
 		          _delta2->prio() : std::numeric_limits<Value>::max();
 		        Value d3 = !_delta3->empty() ?
 		          _delta3->prio() : std::numeric_limits<Value>::max();
@@ -2903,13 +3200,13 @@
 		              int right_tree = _tree_set->find(right_blossom);
 		              if (left_tree == right_tree) {
 		                shrinkOnEdge(e, left_tree);
 		              } else {
 		                augmentOnEdge(e);
-		                unmatched -= 2;
+		                _unmatched -= 2;
+		              }
+		            }
 		          } break;
 		        case D4:
 		          splitBlossom(_delta4->top());
 		          break;
@@ -2922,17 +3219,17 @@
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
 		    ///
 		    /// \note mwpm.run() is just a shortcut of the following code.
 		    /// \code
-		    ///   mwpm.init();
+		    ///   mwpm.fractionalInit();
 		    ///   mwpm.start();
 		    /// \endcode
 		    bool run() {
-		      init();
+		      fractionalInit();
 		      return start();
+		    }
 		    /// @}
 		    /// \name Primal Solution

test/matching_test.cc

Show inline comments

@@ -398,27 +398,51 @@
 		    SmartGraph::EdgeMap<int> weight(graph);
 		    istringstream lgfs(lgf[i]);
 		    graphReader(graph, lgfs).
 		      edgeMap("weight", weight).run();
 		    bool perfect;
+		    {
 		    MaxMatching<SmartGraph> mm(graph);
 		    mm.run();
 		    checkMatching(graph, mm);
 		      perfect = 2 * mm.matchingSize() == countNodes(graph);
+		    }
+		    {
 		    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
 		    mwm.run();
 		    checkWeightedMatching(graph, weight, mwm);
+		    }
+		    {
 		      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
 		      mwm.init();
 		      mwm.start();
 		      checkWeightedMatching(graph, weight, mwm);
+		    }
+		    {
 		    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
-		    bool perfect = mwpm.run();
+		      bool result = mwpm.run();
 		    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
 		          "Perfect matching found");
+		      check(result == perfect, "Perfect matching found");
 		    if (perfect) {
 		      checkWeightedPerfectMatching(graph, weight, mwpm);
+		    }
+		  }
+		    {
 		      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
 		      mwpm.init();
 		      bool result = mwpm.start();
 		      check(result == perfect, "Perfect matching found");
 		      if (perfect) {
 		        checkWeightedPerfectMatching(graph, weight, mwpm);
+		      }
+		    }
+		  }
 		  return 0;
+		}

0 comments (0 inline)

RhodeCode

Login to your account