Fractional matching initialization of weighted matchings (#314)
authorBalazs Dezso <deba@inf.elte.hu>
Sat, 26 Sep 2009 10:17:31 +0200 (2009-09-26)
changeset 87061120524af27
parent 869 636dadefe1e6
child 871 86613aa28a0c
Fractional matching initialization of weighted matchings (#314)
lemon/matching.h
test/matching_test.cc
     1.1 --- a/lemon/matching.h	Fri Sep 25 21:51:36 2009 +0200
     1.2 +++ b/lemon/matching.h	Sat Sep 26 10:17:31 2009 +0200
     1.3 @@ -28,6 +28,7 @@
     1.4  #include <lemon/unionfind.h>
     1.5  #include <lemon/bin_heap.h>
     1.6  #include <lemon/maps.h>
     1.7 +#include <lemon/fractional_matching.h>
     1.8  
     1.9  ///\ingroup matching
    1.10  ///\file
    1.11 @@ -797,6 +798,10 @@
    1.12      BinHeap<Value, IntIntMap> *_delta4;
    1.13  
    1.14      Value _delta_sum;
    1.15 +    int _unmatched;
    1.16 +
    1.17 +    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
    1.18 +    FractionalMatching *_fractional;
    1.19  
    1.20      void createStructures() {
    1.21        _node_num = countNodes(_graph);
    1.22 @@ -1559,10 +1564,16 @@
    1.23          _delta3_index(0), _delta3(0),
    1.24          _delta4_index(0), _delta4(0),
    1.25  
    1.26 -        _delta_sum() {}
    1.27 +        _delta_sum(), _unmatched(0),
    1.28 +
    1.29 +        _fractional(0)
    1.30 +    {}
    1.31  
    1.32      ~MaxWeightedMatching() {
    1.33        destroyStructures();
    1.34 +      if (_fractional) {
    1.35 +        delete _fractional;
    1.36 +      }
    1.37      }
    1.38  
    1.39      /// \name Execution Control
    1.40 @@ -1591,6 +1602,8 @@
    1.41          (*_delta4_index)[i] = _delta4->PRE_HEAP;
    1.42        }
    1.43  
    1.44 +      _unmatched = _node_num;
    1.45 +
    1.46        int index = 0;
    1.47        for (NodeIt n(_graph); n != INVALID; ++n) {
    1.48          Value max = 0;
    1.49 @@ -1625,18 +1638,155 @@
    1.50        }
    1.51      }
    1.52  
    1.53 +    /// \brief Initialize the algorithm with fractional matching
    1.54 +    ///
    1.55 +    /// This function initializes the algorithm with a fractional
    1.56 +    /// matching. This initialization is also called jumpstart heuristic.
    1.57 +    void fractionalInit() {
    1.58 +      createStructures();
    1.59 +      
    1.60 +      if (_fractional == 0) {
    1.61 +        _fractional = new FractionalMatching(_graph, _weight, false);
    1.62 +      }
    1.63 +      _fractional->run();
    1.64 +
    1.65 +      for (ArcIt e(_graph); e != INVALID; ++e) {
    1.66 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
    1.67 +      }
    1.68 +      for (NodeIt n(_graph); n != INVALID; ++n) {
    1.69 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
    1.70 +      }
    1.71 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
    1.72 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
    1.73 +      }
    1.74 +      for (int i = 0; i < _blossom_num; ++i) {
    1.75 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
    1.76 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
    1.77 +      }
    1.78 +
    1.79 +      _unmatched = 0;
    1.80 +
    1.81 +      int index = 0;
    1.82 +      for (NodeIt n(_graph); n != INVALID; ++n) {
    1.83 +        Value pot = _fractional->nodeValue(n);
    1.84 +        (*_node_index)[n] = index;
    1.85 +        (*_node_data)[index].pot = pot;
    1.86 +        int blossom =
    1.87 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
    1.88 +
    1.89 +        (*_blossom_data)[blossom].status = MATCHED;
    1.90 +        (*_blossom_data)[blossom].pred = INVALID;
    1.91 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
    1.92 +        if (_fractional->matching(n) == INVALID) {
    1.93 +          (*_blossom_data)[blossom].base = n;
    1.94 +        }
    1.95 +        (*_blossom_data)[blossom].pot = 0;
    1.96 +        (*_blossom_data)[blossom].offset = 0;
    1.97 +        ++index;
    1.98 +      }
    1.99 +
   1.100 +      typename Graph::template NodeMap<bool> processed(_graph, false);
   1.101 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.102 +        if (processed[n]) continue;
   1.103 +        processed[n] = true;
   1.104 +        if (_fractional->matching(n) == INVALID) continue;
   1.105 +        int num = 1;
   1.106 +        Node v = _graph.target(_fractional->matching(n));
   1.107 +        while (n != v) {
   1.108 +          processed[v] = true;
   1.109 +          v = _graph.target(_fractional->matching(v));
   1.110 +          ++num;
   1.111 +        }
   1.112 +
   1.113 +        if (num % 2 == 1) {
   1.114 +          std::vector<int> subblossoms(num);
   1.115 +
   1.116 +          subblossoms[--num] = _blossom_set->find(n);
   1.117 +          _delta1->push(n, _fractional->nodeValue(n));
   1.118 +          v = _graph.target(_fractional->matching(n));
   1.119 +          while (n != v) {
   1.120 +            subblossoms[--num] = _blossom_set->find(v);
   1.121 +            _delta1->push(v, _fractional->nodeValue(v));
   1.122 +            v = _graph.target(_fractional->matching(v));            
   1.123 +          }
   1.124 +          
   1.125 +          int surface = 
   1.126 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
   1.127 +          (*_blossom_data)[surface].status = EVEN;
   1.128 +          (*_blossom_data)[surface].pred = INVALID;
   1.129 +          (*_blossom_data)[surface].next = INVALID;
   1.130 +          (*_blossom_data)[surface].pot = 0;
   1.131 +          (*_blossom_data)[surface].offset = 0;
   1.132 +          
   1.133 +          _tree_set->insert(surface);
   1.134 +          ++_unmatched;
   1.135 +        }
   1.136 +      }
   1.137 +
   1.138 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.139 +        int si = (*_node_index)[_graph.u(e)];
   1.140 +        int sb = _blossom_set->find(_graph.u(e));
   1.141 +        int ti = (*_node_index)[_graph.v(e)];
   1.142 +        int tb = _blossom_set->find(_graph.v(e));
   1.143 +        if ((*_blossom_data)[sb].status == EVEN &&
   1.144 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
   1.145 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
   1.146 +                            dualScale * _weight[e]) / 2);
   1.147 +        }
   1.148 +      }
   1.149 +
   1.150 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.151 +        int nb = _blossom_set->find(n);
   1.152 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
   1.153 +        int ni = (*_node_index)[n];
   1.154 +
   1.155 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.156 +          Node v = _graph.target(e);
   1.157 +          int vb = _blossom_set->find(v);
   1.158 +          int vi = (*_node_index)[v];
   1.159 +
   1.160 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.161 +            dualScale * _weight[e];
   1.162 +
   1.163 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.164 +
   1.165 +            int vt = _tree_set->find(vb);
   1.166 +
   1.167 +            typename std::map<int, Arc>::iterator it =
   1.168 +              (*_node_data)[ni].heap_index.find(vt);
   1.169 +
   1.170 +            if (it != (*_node_data)[ni].heap_index.end()) {
   1.171 +              if ((*_node_data)[ni].heap[it->second] > rw) {
   1.172 +                (*_node_data)[ni].heap.replace(it->second, e);
   1.173 +                (*_node_data)[ni].heap.decrease(e, rw);
   1.174 +                it->second = e;
   1.175 +              }
   1.176 +            } else {
   1.177 +              (*_node_data)[ni].heap.push(e, rw);
   1.178 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
   1.179 +            }
   1.180 +          }
   1.181 +        }
   1.182 +            
   1.183 +        if (!(*_node_data)[ni].heap.empty()) {
   1.184 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.185 +          _delta2->push(nb, _blossom_set->classPrio(nb));
   1.186 +        }
   1.187 +      }
   1.188 +    }
   1.189 +
   1.190      /// \brief Start the algorithm
   1.191      ///
   1.192      /// This function starts the algorithm.
   1.193      ///
   1.194 -    /// \pre \ref init() must be called before using this function.
   1.195 +    /// \pre \ref init() or \ref fractionalInit() must be called
   1.196 +    /// before using this function.
   1.197      void start() {
   1.198        enum OpType {
   1.199          D1, D2, D3, D4
   1.200        };
   1.201  
   1.202 -      int unmatched = _node_num;
   1.203 -      while (unmatched > 0) {
   1.204 +      while (_unmatched > 0) {
   1.205          Value d1 = !_delta1->empty() ?
   1.206            _delta1->prio() : std::numeric_limits<Value>::max();
   1.207  
   1.208 @@ -1659,7 +1809,7 @@
   1.209            {
   1.210              Node n = _delta1->top();
   1.211              unmatchNode(n);
   1.212 -            --unmatched;
   1.213 +            --_unmatched;
   1.214            }
   1.215            break;
   1.216          case D2:
   1.217 @@ -1669,7 +1819,7 @@
   1.218              Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
   1.219              if ((*_blossom_data)[blossom].next == INVALID) {
   1.220                augmentOnArc(a);
   1.221 -              --unmatched;
   1.222 +              --_unmatched;
   1.223              } else {
   1.224                extendOnArc(a);
   1.225              }
   1.226 @@ -1692,7 +1842,7 @@
   1.227                  shrinkOnEdge(e, left_tree);
   1.228                } else {
   1.229                  augmentOnEdge(e);
   1.230 -                unmatched -= 2;
   1.231 +                _unmatched -= 2;
   1.232                }
   1.233              }
   1.234            } break;
   1.235 @@ -1710,11 +1860,11 @@
   1.236      ///
   1.237      /// \note mwm.run() is just a shortcut of the following code.
   1.238      /// \code
   1.239 -    ///   mwm.init();
   1.240 +    ///   mwm.fractionalInit();
   1.241      ///   mwm.start();
   1.242      /// \endcode
   1.243      void run() {
   1.244 -      init();
   1.245 +      fractionalInit();
   1.246        start();
   1.247      }
   1.248  
   1.249 @@ -2074,6 +2224,11 @@
   1.250      BinHeap<Value, IntIntMap> *_delta4;
   1.251  
   1.252      Value _delta_sum;
   1.253 +    int _unmatched;
   1.254 +
   1.255 +    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 
   1.256 +    FractionalMatching;
   1.257 +    FractionalMatching *_fractional;
   1.258  
   1.259      void createStructures() {
   1.260        _node_num = countNodes(_graph);
   1.261 @@ -2789,10 +2944,16 @@
   1.262          _delta3_index(0), _delta3(0),
   1.263          _delta4_index(0), _delta4(0),
   1.264  
   1.265 -        _delta_sum() {}
   1.266 +        _delta_sum(), _unmatched(0),
   1.267 +
   1.268 +        _fractional(0)
   1.269 +    {}
   1.270  
   1.271      ~MaxWeightedPerfectMatching() {
   1.272        destroyStructures();
   1.273 +      if (_fractional) {
   1.274 +        delete _fractional;
   1.275 +      }
   1.276      }
   1.277  
   1.278      /// \name Execution Control
   1.279 @@ -2818,6 +2979,8 @@
   1.280          (*_delta4_index)[i] = _delta4->PRE_HEAP;
   1.281        }
   1.282  
   1.283 +      _unmatched = _node_num;
   1.284 +
   1.285        int index = 0;
   1.286        for (NodeIt n(_graph); n != INVALID; ++n) {
   1.287          Value max = - std::numeric_limits<Value>::max();
   1.288 @@ -2851,18 +3014,152 @@
   1.289        }
   1.290      }
   1.291  
   1.292 +    /// \brief Initialize the algorithm with fractional matching
   1.293 +    ///
   1.294 +    /// This function initializes the algorithm with a fractional
   1.295 +    /// matching. This initialization is also called jumpstart heuristic.
   1.296 +    void fractionalInit() {
   1.297 +      createStructures();
   1.298 +      
   1.299 +      if (_fractional == 0) {
   1.300 +        _fractional = new FractionalMatching(_graph, _weight, false);
   1.301 +      }
   1.302 +      if (!_fractional->run()) {
   1.303 +        _unmatched = -1;
   1.304 +        return;
   1.305 +      }
   1.306 +
   1.307 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   1.308 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
   1.309 +      }
   1.310 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.311 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
   1.312 +      }
   1.313 +      for (int i = 0; i < _blossom_num; ++i) {
   1.314 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
   1.315 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
   1.316 +      }
   1.317 +
   1.318 +      _unmatched = 0;
   1.319 +
   1.320 +      int index = 0;
   1.321 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.322 +        Value pot = _fractional->nodeValue(n);
   1.323 +        (*_node_index)[n] = index;
   1.324 +        (*_node_data)[index].pot = pot;
   1.325 +        int blossom =
   1.326 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
   1.327 +
   1.328 +        (*_blossom_data)[blossom].status = MATCHED;
   1.329 +        (*_blossom_data)[blossom].pred = INVALID;
   1.330 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
   1.331 +        (*_blossom_data)[blossom].pot = 0;
   1.332 +        (*_blossom_data)[blossom].offset = 0;
   1.333 +        ++index;
   1.334 +      }
   1.335 +
   1.336 +      typename Graph::template NodeMap<bool> processed(_graph, false);
   1.337 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.338 +        if (processed[n]) continue;
   1.339 +        processed[n] = true;
   1.340 +        if (_fractional->matching(n) == INVALID) continue;
   1.341 +        int num = 1;
   1.342 +        Node v = _graph.target(_fractional->matching(n));
   1.343 +        while (n != v) {
   1.344 +          processed[v] = true;
   1.345 +          v = _graph.target(_fractional->matching(v));
   1.346 +          ++num;
   1.347 +        }
   1.348 +
   1.349 +        if (num % 2 == 1) {
   1.350 +          std::vector<int> subblossoms(num);
   1.351 +
   1.352 +          subblossoms[--num] = _blossom_set->find(n);
   1.353 +          v = _graph.target(_fractional->matching(n));
   1.354 +          while (n != v) {
   1.355 +            subblossoms[--num] = _blossom_set->find(v);
   1.356 +            v = _graph.target(_fractional->matching(v));            
   1.357 +          }
   1.358 +          
   1.359 +          int surface = 
   1.360 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
   1.361 +          (*_blossom_data)[surface].status = EVEN;
   1.362 +          (*_blossom_data)[surface].pred = INVALID;
   1.363 +          (*_blossom_data)[surface].next = INVALID;
   1.364 +          (*_blossom_data)[surface].pot = 0;
   1.365 +          (*_blossom_data)[surface].offset = 0;
   1.366 +          
   1.367 +          _tree_set->insert(surface);
   1.368 +          ++_unmatched;
   1.369 +        }
   1.370 +      }
   1.371 +
   1.372 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.373 +        int si = (*_node_index)[_graph.u(e)];
   1.374 +        int sb = _blossom_set->find(_graph.u(e));
   1.375 +        int ti = (*_node_index)[_graph.v(e)];
   1.376 +        int tb = _blossom_set->find(_graph.v(e));
   1.377 +        if ((*_blossom_data)[sb].status == EVEN &&
   1.378 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
   1.379 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
   1.380 +                            dualScale * _weight[e]) / 2);
   1.381 +        }
   1.382 +      }
   1.383 +
   1.384 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.385 +        int nb = _blossom_set->find(n);
   1.386 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
   1.387 +        int ni = (*_node_index)[n];
   1.388 +
   1.389 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.390 +          Node v = _graph.target(e);
   1.391 +          int vb = _blossom_set->find(v);
   1.392 +          int vi = (*_node_index)[v];
   1.393 +
   1.394 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.395 +            dualScale * _weight[e];
   1.396 +
   1.397 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.398 +
   1.399 +            int vt = _tree_set->find(vb);
   1.400 +
   1.401 +            typename std::map<int, Arc>::iterator it =
   1.402 +              (*_node_data)[ni].heap_index.find(vt);
   1.403 +
   1.404 +            if (it != (*_node_data)[ni].heap_index.end()) {
   1.405 +              if ((*_node_data)[ni].heap[it->second] > rw) {
   1.406 +                (*_node_data)[ni].heap.replace(it->second, e);
   1.407 +                (*_node_data)[ni].heap.decrease(e, rw);
   1.408 +                it->second = e;
   1.409 +              }
   1.410 +            } else {
   1.411 +              (*_node_data)[ni].heap.push(e, rw);
   1.412 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
   1.413 +            }
   1.414 +          }
   1.415 +        }
   1.416 +            
   1.417 +        if (!(*_node_data)[ni].heap.empty()) {
   1.418 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.419 +          _delta2->push(nb, _blossom_set->classPrio(nb));
   1.420 +        }
   1.421 +      }
   1.422 +    }
   1.423 +
   1.424      /// \brief Start the algorithm
   1.425      ///
   1.426      /// This function starts the algorithm.
   1.427      ///
   1.428 -    /// \pre \ref init() must be called before using this function.
   1.429 +    /// \pre \ref init() or \ref fractionalInit() must be called before
   1.430 +    /// using this function.
   1.431      bool start() {
   1.432        enum OpType {
   1.433          D2, D3, D4
   1.434        };
   1.435  
   1.436 -      int unmatched = _node_num;
   1.437 -      while (unmatched > 0) {
   1.438 +      if (_unmatched == -1) return false;
   1.439 +
   1.440 +      while (_unmatched > 0) {
   1.441          Value d2 = !_delta2->empty() ?
   1.442            _delta2->prio() : std::numeric_limits<Value>::max();
   1.443  
   1.444 @@ -2906,7 +3203,7 @@
   1.445                  shrinkOnEdge(e, left_tree);
   1.446                } else {
   1.447                  augmentOnEdge(e);
   1.448 -                unmatched -= 2;
   1.449 +                _unmatched -= 2;
   1.450                }
   1.451              }
   1.452            } break;
   1.453 @@ -2925,11 +3222,11 @@
   1.454      ///
   1.455      /// \note mwpm.run() is just a shortcut of the following code.
   1.456      /// \code
   1.457 -    ///   mwpm.init();
   1.458 +    ///   mwpm.fractionalInit();
   1.459      ///   mwpm.start();
   1.460      /// \endcode
   1.461      bool run() {
   1.462 -      init();
   1.463 +      fractionalInit();
   1.464        return start();
   1.465      }
   1.466  
     2.1 --- a/test/matching_test.cc	Fri Sep 25 21:51:36 2009 +0200
     2.2 +++ b/test/matching_test.cc	Sat Sep 26 10:17:31 2009 +0200
     2.3 @@ -401,22 +401,46 @@
     2.4      graphReader(graph, lgfs).
     2.5        edgeMap("weight", weight).run();
     2.6  
     2.7 -    MaxMatching<SmartGraph> mm(graph);
     2.8 -    mm.run();
     2.9 -    checkMatching(graph, mm);
    2.10 +    bool perfect;
    2.11 +    {
    2.12 +      MaxMatching<SmartGraph> mm(graph);
    2.13 +      mm.run();
    2.14 +      checkMatching(graph, mm);
    2.15 +      perfect = 2 * mm.matchingSize() == countNodes(graph);
    2.16 +    }
    2.17  
    2.18 -    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    2.19 -    mwm.run();
    2.20 -    checkWeightedMatching(graph, weight, mwm);
    2.21 +    {
    2.22 +      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    2.23 +      mwm.run();
    2.24 +      checkWeightedMatching(graph, weight, mwm);
    2.25 +    }
    2.26  
    2.27 -    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    2.28 -    bool perfect = mwpm.run();
    2.29 +    {
    2.30 +      MaxWeightedMatching<SmartGraph> mwm(graph, weight);
    2.31 +      mwm.init();
    2.32 +      mwm.start();
    2.33 +      checkWeightedMatching(graph, weight, mwm);
    2.34 +    }
    2.35  
    2.36 -    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
    2.37 -          "Perfect matching found");
    2.38 +    {
    2.39 +      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    2.40 +      bool result = mwpm.run();
    2.41 +      
    2.42 +      check(result == perfect, "Perfect matching found");
    2.43 +      if (perfect) {
    2.44 +        checkWeightedPerfectMatching(graph, weight, mwpm);
    2.45 +      }
    2.46 +    }
    2.47  
    2.48 -    if (perfect) {
    2.49 -      checkWeightedPerfectMatching(graph, weight, mwpm);
    2.50 +    {
    2.51 +      MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
    2.52 +      mwpm.init();
    2.53 +      bool result = mwpm.start();
    2.54 +      
    2.55 +      check(result == perfect, "Perfect matching found");
    2.56 +      if (perfect) {
    2.57 +        checkWeightedPerfectMatching(graph, weight, mwpm);
    2.58 +      }
    2.59      }
    2.60    }
    2.61