gravatar
deba@inf.elte.hu
deba@inf.elte.hu
General improvements in weighted matching algorithms (#314) - Fix include guard - Uniform handling of MATCHED and UNMATCHED blossoms - Prefer operations which decrease the number of trees - Fix improper use of '/=' The solved problems did not cause wrong solution.
0 1 0
default
1 file changed with 38 insertions and 157 deletions:
↑ Collapse diff ↑
Show white space 12 line context
... ...
@@ -13,14 +13,14 @@
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19
#ifndef LEMON_MAX_MATCHING_H
20
#define LEMON_MAX_MATCHING_H
19
#ifndef LEMON_MATCHING_H
20
#define LEMON_MATCHING_H
21 21

	
22 22
#include <vector>
23 23
#include <queue>
24 24
#include <set>
25 25
#include <limits>
26 26

	
... ...
@@ -742,13 +742,13 @@
742 742
    int _node_num;
743 743
    int _blossom_num;
744 744

	
745 745
    typedef RangeMap<int> IntIntMap;
746 746

	
747 747
    enum Status {
748
      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
748
      EVEN = -1, MATCHED = 0, ODD = 1
749 749
    };
750 750

	
751 751
    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
752 752
    struct BlossomData {
753 753
      int tree;
754 754
      Status status;
... ...
@@ -841,15 +841,12 @@
841 841
        _delta4_index = new IntIntMap(_blossom_num);
842 842
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
843 843
      }
844 844
    }
845 845

	
846 846
    void destroyStructures() {
847
      _node_num = countNodes(_graph);
848
      _blossom_num = _node_num * 3 / 2;
849

	
850 847
      if (_matching) {
851 848
        delete _matching;
852 849
      }
853 850
      if (_node_potential) {
854 851
        delete _node_potential;
855 852
      }
... ...
@@ -919,16 +916,12 @@
919 916
            dualScale * _weight[e];
920 917

	
921 918
          if ((*_blossom_data)[vb].status == EVEN) {
922 919
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
923 920
              _delta3->push(e, rw / 2);
924 921
            }
925
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
926
            if (_delta3->state(e) != _delta3->IN_HEAP) {
927
              _delta3->push(e, rw);
928
            }
929 922
          } else {
930 923
            typename std::map<int, Arc>::iterator it =
931 924
              (*_node_data)[vi].heap_index.find(tree);
932 925

	
933 926
            if (it != (*_node_data)[vi].heap_index.end()) {
934 927
              if ((*_node_data)[vi].heap[it->second] > rw) {
... ...
@@ -1033,17 +1026,12 @@
1033 1026
                           (*_blossom_data)[blossom].offset){
1034 1027
                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1035 1028
                                   (*_blossom_data)[blossom].offset);
1036 1029
                }
1037 1030
              }
1038 1031
            }
1039

	
1040
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1041
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1042
              _delta3->erase(e);
1043
            }
1044 1032
          } else {
1045 1033

	
1046 1034
            typename std::map<int, Arc>::iterator it =
1047 1035
              (*_node_data)[vi].heap_index.find(tree);
1048 1036

	
1049 1037
            if (it != (*_node_data)[vi].heap_index.end()) {
... ...
@@ -1114,16 +1102,12 @@
1114 1102
            dualScale * _weight[e];
1115 1103

	
1116 1104
          if ((*_blossom_data)[vb].status == EVEN) {
1117 1105
            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
1118 1106
              _delta3->push(e, rw / 2);
1119 1107
            }
1120
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1121
            if (_delta3->state(e) != _delta3->IN_HEAP) {
1122
              _delta3->push(e, rw);
1123
            }
1124 1108
          } else {
1125 1109

	
1126 1110
            typename std::map<int, Arc>::iterator it =
1127 1111
              (*_node_data)[vi].heap_index.find(tree);
1128 1112

	
1129 1113
            if (it != (*_node_data)[vi].heap_index.end()) {
... ...
@@ -1154,107 +1138,12 @@
1154 1138
          }
1155 1139
        }
1156 1140
      }
1157 1141
      (*_blossom_data)[blossom].offset = 0;
1158 1142
    }
1159 1143

	
1160

	
1161
    void matchedToUnmatched(int blossom) {
1162
      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
1163
        _delta2->erase(blossom);
1164
      }
1165

	
1166
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1167
           n != INVALID; ++n) {
1168
        int ni = (*_node_index)[n];
1169

	
1170
        _blossom_set->increase(n, std::numeric_limits<Value>::max());
1171

	
1172
        (*_node_data)[ni].heap.clear();
1173
        (*_node_data)[ni].heap_index.clear();
1174

	
1175
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1176
          Node v = _graph.target(e);
1177
          int vb = _blossom_set->find(v);
1178
          int vi = (*_node_index)[v];
1179

	
1180
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1181
            dualScale * _weight[e];
1182

	
1183
          if ((*_blossom_data)[vb].status == EVEN) {
1184
            if (_delta3->state(e) != _delta3->IN_HEAP) {
1185
              _delta3->push(e, rw);
1186
            }
1187
          }
1188
        }
1189
      }
1190
    }
1191

	
1192
    void unmatchedToMatched(int blossom) {
1193
      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
1194
           n != INVALID; ++n) {
1195
        int ni = (*_node_index)[n];
1196

	
1197
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
1198
          Node v = _graph.source(e);
1199
          int vb = _blossom_set->find(v);
1200
          int vi = (*_node_index)[v];
1201

	
1202
          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
1203
            dualScale * _weight[e];
1204

	
1205
          if (vb == blossom) {
1206
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1207
              _delta3->erase(e);
1208
            }
1209
          } else if ((*_blossom_data)[vb].status == EVEN) {
1210

	
1211
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1212
              _delta3->erase(e);
1213
            }
1214

	
1215
            int vt = _tree_set->find(vb);
1216

	
1217
            Arc r = _graph.oppositeArc(e);
1218

	
1219
            typename std::map<int, Arc>::iterator it =
1220
              (*_node_data)[ni].heap_index.find(vt);
1221

	
1222
            if (it != (*_node_data)[ni].heap_index.end()) {
1223
              if ((*_node_data)[ni].heap[it->second] > rw) {
1224
                (*_node_data)[ni].heap.replace(it->second, r);
1225
                (*_node_data)[ni].heap.decrease(r, rw);
1226
                it->second = r;
1227
              }
1228
            } else {
1229
              (*_node_data)[ni].heap.push(r, rw);
1230
              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
1231
            }
1232

	
1233
            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
1234
              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
1235

	
1236
              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
1237
                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
1238
                             (*_blossom_data)[blossom].offset);
1239
              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
1240
                         (*_blossom_data)[blossom].offset){
1241
                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
1242
                                 (*_blossom_data)[blossom].offset);
1243
              }
1244
            }
1245

	
1246
          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
1247
            if (_delta3->state(e) == _delta3->IN_HEAP) {
1248
              _delta3->erase(e);
1249
            }
1250
          }
1251
        }
1252
      }
1253
    }
1254

	
1255 1144
    void alternatePath(int even, int tree) {
1256 1145
      int odd;
1257 1146

	
1258 1147
      evenToMatched(even, tree);
1259 1148
      (*_blossom_data)[even].status = MATCHED;
1260 1149

	
... ...
@@ -1291,45 +1180,48 @@
1291 1180
      int blossom = _blossom_set->find(node);
1292 1181
      int tree = _tree_set->find(blossom);
1293 1182

	
1294 1183
      alternatePath(blossom, tree);
1295 1184
      destroyTree(tree);
1296 1185

	
1297
      (*_blossom_data)[blossom].status = UNMATCHED;
1298 1186
      (*_blossom_data)[blossom].base = node;
1299
      matchedToUnmatched(blossom);
1300
    }
1301

	
1187
      (*_blossom_data)[blossom].next = INVALID;
1188
    }
1302 1189

	
1303 1190
    void augmentOnEdge(const Edge& edge) {
1304 1191

	
1305 1192
      int left = _blossom_set->find(_graph.u(edge));
1306 1193
      int right = _blossom_set->find(_graph.v(edge));
1307 1194

	
1308
      if ((*_blossom_data)[left].status == EVEN) {
1309 1195
        int left_tree = _tree_set->find(left);
1310 1196
        alternatePath(left, left_tree);
1311 1197
        destroyTree(left_tree);
1312
      } else {
1313
        (*_blossom_data)[left].status = MATCHED;
1314
        unmatchedToMatched(left);
1315
      }
1316

	
1317
      if ((*_blossom_data)[right].status == EVEN) {
1198

	
1318 1199
        int right_tree = _tree_set->find(right);
1319 1200
        alternatePath(right, right_tree);
1320 1201
        destroyTree(right_tree);
1321
      } else {
1322
        (*_blossom_data)[right].status = MATCHED;
1323
        unmatchedToMatched(right);
1324
      }
1325 1202

	
1326 1203
      (*_blossom_data)[left].next = _graph.direct(edge, true);
1327 1204
      (*_blossom_data)[right].next = _graph.direct(edge, false);
1328 1205
    }
1329 1206

	
1207
    void augmentOnArc(const Arc& arc) {
1208

	
1209
      int left = _blossom_set->find(_graph.source(arc));
1210
      int right = _blossom_set->find(_graph.target(arc));
1211

	
1212
      (*_blossom_data)[left].status = MATCHED;
1213

	
1214
      int right_tree = _tree_set->find(right);
1215
      alternatePath(right, right_tree);
1216
      destroyTree(right_tree);
1217

	
1218
      (*_blossom_data)[left].next = arc;
1219
      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
1220
    }
1221

	
1330 1222
    void extendOnArc(const Arc& arc) {
1331 1223
      int base = _blossom_set->find(_graph.target(arc));
1332 1224
      int tree = _tree_set->find(base);
1333 1225

	
1334 1226
      int odd = _blossom_set->find(_graph.source(arc));
1335 1227
      _tree_set->insert(odd, tree);
... ...
@@ -1626,13 +1518,13 @@
1626 1518
      std::vector<int> blossoms;
1627 1519
      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
1628 1520
        blossoms.push_back(c);
1629 1521
      }
1630 1522

	
1631 1523
      for (int i = 0; i < int(blossoms.size()); ++i) {
1632
        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
1524
        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
1633 1525

	
1634 1526
          Value offset = (*_blossom_data)[blossoms[i]].offset;
1635 1527
          (*_blossom_data)[blossoms[i]].pot += 2 * offset;
1636 1528
          for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
1637 1529
               n != INVALID; ++n) {
1638 1530
            (*_node_data)[(*_node_index)[n]].pot -= offset;
... ...
@@ -1754,58 +1646,50 @@
1754 1646
        Value d3 = !_delta3->empty() ?
1755 1647
          _delta3->prio() : std::numeric_limits<Value>::max();
1756 1648

	
1757 1649
        Value d4 = !_delta4->empty() ?
1758 1650
          _delta4->prio() : std::numeric_limits<Value>::max();
1759 1651

	
1760
        _delta_sum = d1; OpType ot = D1;
1652
        _delta_sum = d3; OpType ot = D3;
1653
        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1761 1654
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1762
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
1763 1655
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
1764 1656

	
1765

	
1766 1657
        switch (ot) {
1767 1658
        case D1:
1768 1659
          {
1769 1660
            Node n = _delta1->top();
1770 1661
            unmatchNode(n);
1771 1662
            --unmatched;
1772 1663
          }
1773 1664
          break;
1774 1665
        case D2:
1775 1666
          {
1776 1667
            int blossom = _delta2->top();
1777 1668
            Node n = _blossom_set->classTop(blossom);
1778
            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
1779
            extendOnArc(e);
1669
            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
1670
            if ((*_blossom_data)[blossom].next == INVALID) {
1671
              augmentOnArc(a);
1672
              --unmatched;
1673
            } else {
1674
              extendOnArc(a);
1675
            }
1780 1676
          }
1781 1677
          break;
1782 1678
        case D3:
1783 1679
          {
1784 1680
            Edge e = _delta3->top();
1785 1681

	
1786 1682
            int left_blossom = _blossom_set->find(_graph.u(e));
1787 1683
            int right_blossom = _blossom_set->find(_graph.v(e));
1788 1684

	
1789 1685
            if (left_blossom == right_blossom) {
1790 1686
              _delta3->pop();
1791 1687
            } else {
1792
              int left_tree;
1793
              if ((*_blossom_data)[left_blossom].status == EVEN) {
1794
                left_tree = _tree_set->find(left_blossom);
1795
              } else {
1796
                left_tree = -1;
1797
                ++unmatched;
1798
              }
1799
              int right_tree;
1800
              if ((*_blossom_data)[right_blossom].status == EVEN) {
1801
                right_tree = _tree_set->find(right_blossom);
1802
              } else {
1803
                right_tree = -1;
1804
                ++unmatched;
1805
              }
1688
              int left_tree = _tree_set->find(left_blossom);
1689
              int right_tree = _tree_set->find(right_blossom);
1806 1690

	
1807 1691
              if (left_tree == right_tree) {
1808 1692
                shrinkOnEdge(e, left_tree);
1809 1693
              } else {
1810 1694
                augmentOnEdge(e);
1811 1695
                unmatched -= 2;
... ...
@@ -1853,13 +1737,13 @@
1853 1737
      Value sum = 0;
1854 1738
      for (NodeIt n(_graph); n != INVALID; ++n) {
1855 1739
        if ((*_matching)[n] != INVALID) {
1856 1740
          sum += _weight[(*_matching)[n]];
1857 1741
        }
1858 1742
      }
1859
      return sum /= 2;
1743
      return sum / 2;
1860 1744
    }
1861 1745

	
1862 1746
    /// \brief Return the size (cardinality) of the matching.
1863 1747
    ///
1864 1748
    /// This function returns the size (cardinality) of the found matching.
1865 1749
    ///
... ...
@@ -2230,15 +2114,12 @@
2230 2114
        _delta4_index = new IntIntMap(_blossom_num);
2231 2115
        _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
2232 2116
      }
2233 2117
    }
2234 2118

	
2235 2119
    void destroyStructures() {
2236
      _node_num = countNodes(_graph);
2237
      _blossom_num = _node_num * 3 / 2;
2238

	
2239 2120
      if (_matching) {
2240 2121
        delete _matching;
2241 2122
      }
2242 2123
      if (_node_potential) {
2243 2124
        delete _node_potential;
2244 2125
      }
... ...
@@ -2988,14 +2869,14 @@
2988 2869
        Value d3 = !_delta3->empty() ?
2989 2870
          _delta3->prio() : std::numeric_limits<Value>::max();
2990 2871

	
2991 2872
        Value d4 = !_delta4->empty() ?
2992 2873
          _delta4->prio() : std::numeric_limits<Value>::max();
2993 2874

	
2994
        _delta_sum = d2; OpType ot = D2;
2995
        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
2875
        _delta_sum = d3; OpType ot = D3;
2876
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
2996 2877
        if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
2997 2878

	
2998 2879
        if (_delta_sum == std::numeric_limits<Value>::max()) {
2999 2880
          return false;
3000 2881
        }
3001 2882

	
... ...
@@ -3071,13 +2952,13 @@
3071 2952
      Value sum = 0;
3072 2953
      for (NodeIt n(_graph); n != INVALID; ++n) {
3073 2954
        if ((*_matching)[n] != INVALID) {
3074 2955
          sum += _weight[(*_matching)[n]];
3075 2956
        }
3076 2957
      }
3077
      return sum /= 2;
2958
      return sum / 2;
3078 2959
    }
3079 2960

	
3080 2961
    /// \brief Return \c true if the given edge is in the matching.
3081 2962
    ///
3082 2963
    /// This function returns \c true if the given edge is in the found 
3083 2964
    /// matching.
... ...
@@ -3238,7 +3119,7 @@
3238 3119
    /// @}
3239 3120

	
3240 3121
  };
3241 3122

	
3242 3123
} //END OF NAMESPACE LEMON
3243 3124

	
3244
#endif //LEMON_MAX_MATCHING_H
3125
#endif //LEMON_MATCHING_H
0 comments (0 inline)