1.1 --- a/lemon/bipartite_matching.h Fri Apr 14 18:05:02 2006 +0000
1.2 +++ b/lemon/bipartite_matching.h Fri Apr 14 18:07:33 2006 +0000
1.3 @@ -19,8 +19,10 @@
1.4 #ifndef LEMON_BIPARTITE_MATCHING
1.5 #define LEMON_BIPARTITE_MATCHING
1.6
1.7 -#include <lemon/bpugraph_adaptor.h>
1.8 -#include <lemon/bfs.h>
1.9 +#include <functional>
1.10 +
1.11 +#include <lemon/bin_heap.h>
1.12 +#include <lemon/maps.h>
1.13
1.14 #include <iostream>
1.15
1.16 @@ -35,7 +37,7 @@
1.17 /// \brief Bipartite Max Cardinality Matching algorithm
1.18 ///
1.19 /// Bipartite Max Cardinality Matching algorithm. This class implements
1.20 - /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
1.21 + /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
1.22 /// complexity.
1.23 template <typename BpUGraph>
1.24 class MaxBipartiteMatching {
1.25 @@ -83,7 +85,7 @@
1.26 for (BNodeIt it(*graph); it != INVALID; ++it) {
1.27 bnode_matching[it] = INVALID;
1.28 }
1.29 - matching_value = 0;
1.30 + matching_size = 0;
1.31 }
1.32
1.33 /// \brief Initalize the data structures.
1.34 @@ -92,7 +94,7 @@
1.35 /// matching. From this matching sometimes it is faster to get
1.36 /// the matching than from the initial empty matching.
1.37 void greedyInit() {
1.38 - matching_value = 0;
1.39 + matching_size = 0;
1.40 for (BNodeIt it(*graph); it != INVALID; ++it) {
1.41 bnode_matching[it] = INVALID;
1.42 }
1.43 @@ -102,7 +104,7 @@
1.44 if (bnode_matching[graph->bNode(jt)] == INVALID) {
1.45 anode_matching[it] = jt;
1.46 bnode_matching[graph->bNode(jt)] = jt;
1.47 - ++matching_value;
1.48 + ++matching_size;
1.49 break;
1.50 }
1.51 }
1.52 @@ -120,10 +122,10 @@
1.53 for (BNodeIt it(*graph); it != INVALID; ++it) {
1.54 bnode_matching[it] = INVALID;
1.55 }
1.56 - matching_value = 0;
1.57 + matching_size = 0;
1.58 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.59 if (matching[it]) {
1.60 - ++matching_value;
1.61 + ++matching_size;
1.62 anode_matching[graph->aNode(it)] = it;
1.63 bnode_matching[graph->bNode(it)] = it;
1.64 }
1.65 @@ -142,10 +144,10 @@
1.66 for (BNodeIt it(*graph); it != INVALID; ++it) {
1.67 bnode_matching[it] = INVALID;
1.68 }
1.69 - matching_value = 0;
1.70 + matching_size = 0;
1.71 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.72 if (matching[it]) {
1.73 - ++matching_value;
1.74 + ++matching_size;
1.75 if (anode_matching[graph->aNode(it)] != INVALID) {
1.76 return false;
1.77 }
1.78 @@ -246,7 +248,7 @@
1.79
1.80 aused[anode] = true;
1.81 }
1.82 - ++matching_value;
1.83 + ++matching_size;
1.84
1.85 }
1.86 }
1.87 @@ -297,7 +299,7 @@
1.88 anode_matching[anode] = uedge;
1.89
1.90 }
1.91 - ++matching_value;
1.92 + ++matching_size;
1.93 return true;
1.94 } else {
1.95 Node newanode = graph->aNode(bnode_matching[bnode]);
1.96 @@ -407,7 +409,7 @@
1.97 matching[anode_matching[it]] = true;
1.98 }
1.99 }
1.100 - return matching_value;
1.101 + return matching_size;
1.102 }
1.103
1.104 /// \brief Set true all matching uedge in the map and the others to false.
1.105 @@ -419,7 +421,7 @@
1.106 for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.107 matching[it] = it == anode_matching[graph->aNode(it)];
1.108 }
1.109 - return matching_value;
1.110 + return matching_size;
1.111 }
1.112
1.113
1.114 @@ -445,8 +447,8 @@
1.115 /// \brief Gives back the number of the matching edges.
1.116 ///
1.117 /// Gives back the number of the matching edges.
1.118 - int matchingValue() const {
1.119 - return matching_value;
1.120 + int matchingSize() const {
1.121 + return matching_size;
1.122 }
1.123
1.124 /// @}
1.125 @@ -457,7 +459,1022 @@
1.126 BNodeMatchingMap bnode_matching;
1.127 const Graph *graph;
1.128
1.129 - int matching_value;
1.130 + int matching_size;
1.131 +
1.132 + };
1.133 +
1.134 + /// \brief Default traits class for weighted bipartite matching algoritms.
1.135 + ///
1.136 + /// Default traits class for weighted bipartite matching algoritms.
1.137 + /// \param _BpUGraph The bipartite undirected graph type.
1.138 + /// \param _WeightMap Type of weight map.
1.139 + template <typename _BpUGraph, typename _WeightMap>
1.140 + struct WeightedBipartiteMatchingDefaultTraits {
1.141 + /// \brief The type of the weight of the undirected edges.
1.142 + typedef typename _WeightMap::Value Value;
1.143 +
1.144 + /// The undirected bipartite graph type the algorithm runs on.
1.145 + typedef _BpUGraph BpUGraph;
1.146 +
1.147 + /// The map of the edges weights
1.148 + typedef _WeightMap WeightMap;
1.149 +
1.150 + /// \brief The cross reference type used by heap.
1.151 + ///
1.152 + /// The cross reference type used by heap.
1.153 + /// Usually it is \c Graph::NodeMap<int>.
1.154 + typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1.155 +
1.156 + /// \brief Instantiates a HeapCrossRef.
1.157 + ///
1.158 + /// This function instantiates a \ref HeapCrossRef.
1.159 + /// \param graph is the graph, to which we would like to define the
1.160 + /// HeapCrossRef.
1.161 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1.162 + return new HeapCrossRef(graph);
1.163 + }
1.164 +
1.165 + /// \brief The heap type used by weighted matching algorithms.
1.166 + ///
1.167 + /// The heap type used by weighted matching algorithms. It should
1.168 + /// minimize the priorities and the heap's key type is the graph's
1.169 + /// anode graph's node.
1.170 + ///
1.171 + /// \sa BinHeap
1.172 + typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
1.173 +
1.174 + /// \brief Instantiates a Heap.
1.175 + ///
1.176 + /// This function instantiates a \ref Heap.
1.177 + /// \param crossref The cross reference of the heap.
1.178 + static Heap *createHeap(HeapCrossRef& crossref) {
1.179 + return new Heap(crossref);
1.180 + }
1.181 +
1.182 + };
1.183 +
1.184 +
1.185 + /// \ingroup matching
1.186 + ///
1.187 + /// \brief Bipartite Max Weighted Matching algorithm
1.188 + ///
1.189 + /// This class implements the bipartite Max Weighted Matching
1.190 + /// algorithm. It uses the successive shortest path algorithm to
1.191 + /// calculate the maximum weighted matching in the bipartite
1.192 + /// graph. The algorithm can be used also to calculate the maximum
1.193 + /// cardinality maximum weighted matching. The time complexity
1.194 + /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
1.195 + /// heap implementation but this can be improved to
1.196 + /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1.197 + ///
1.198 + /// The algorithm also provides a potential function on the nodes
1.199 + /// which a dual solution of the matching algorithm and it can be
1.200 + /// used to proof the optimality of the given pimal solution.
1.201 +#ifdef DOXYGEN
1.202 + template <typename _BpUGraph, typename _WeightMap, typename _Traits>
1.203 +#else
1.204 + template <typename _BpUGraph,
1.205 + typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
1.206 + typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
1.207 +#endif
1.208 + class MaxWeightedBipartiteMatching {
1.209 + public:
1.210 +
1.211 + typedef _Traits Traits;
1.212 + typedef typename Traits::BpUGraph BpUGraph;
1.213 + typedef typename Traits::WeightMap WeightMap;
1.214 + typedef typename Traits::Value Value;
1.215 +
1.216 + protected:
1.217 +
1.218 + typedef typename Traits::HeapCrossRef HeapCrossRef;
1.219 + typedef typename Traits::Heap Heap;
1.220 +
1.221 +
1.222 + typedef typename BpUGraph::Node Node;
1.223 + typedef typename BpUGraph::ANodeIt ANodeIt;
1.224 + typedef typename BpUGraph::BNodeIt BNodeIt;
1.225 + typedef typename BpUGraph::UEdge UEdge;
1.226 + typedef typename BpUGraph::UEdgeIt UEdgeIt;
1.227 + typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1.228 +
1.229 + typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1.230 + typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1.231 +
1.232 + typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1.233 + typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1.234 +
1.235 +
1.236 + public:
1.237 +
1.238 + /// \brief \ref Exception for uninitialized parameters.
1.239 + ///
1.240 + /// This error represents problems in the initialization
1.241 + /// of the parameters of the algorithms.
1.242 + class UninitializedParameter : public lemon::UninitializedParameter {
1.243 + public:
1.244 + virtual const char* exceptionName() const {
1.245 + return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
1.246 + }
1.247 + };
1.248 +
1.249 + ///\name Named template parameters
1.250 +
1.251 + ///@{
1.252 +
1.253 + template <class H, class CR>
1.254 + struct DefHeapTraits : public Traits {
1.255 + typedef CR HeapCrossRef;
1.256 + typedef H Heap;
1.257 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1.258 + throw UninitializedParameter();
1.259 + }
1.260 + static Heap *createHeap(HeapCrossRef &) {
1.261 + throw UninitializedParameter();
1.262 + }
1.263 + };
1.264 +
1.265 + /// \brief \ref named-templ-param "Named parameter" for setting heap
1.266 + /// and cross reference type
1.267 + ///
1.268 + /// \ref named-templ-param "Named parameter" for setting heap and cross
1.269 + /// reference type
1.270 + template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1.271 + struct DefHeap
1.272 + : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
1.273 + DefHeapTraits<H, CR> > {
1.274 + typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
1.275 + DefHeapTraits<H, CR> > Create;
1.276 + };
1.277 +
1.278 + template <class H, class CR>
1.279 + struct DefStandardHeapTraits : public Traits {
1.280 + typedef CR HeapCrossRef;
1.281 + typedef H Heap;
1.282 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1.283 + return new HeapCrossRef(graph);
1.284 + }
1.285 + static Heap *createHeap(HeapCrossRef &crossref) {
1.286 + return new Heap(crossref);
1.287 + }
1.288 + };
1.289 +
1.290 + /// \brief \ref named-templ-param "Named parameter" for setting heap and
1.291 + /// cross reference type with automatic allocation
1.292 + ///
1.293 + /// \ref named-templ-param "Named parameter" for setting heap and cross
1.294 + /// reference type. It can allocate the heap and the cross reference
1.295 + /// object if the cross reference's constructor waits for the graph as
1.296 + /// parameter and the heap's constructor waits for the cross reference.
1.297 + template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1.298 + struct DefStandardHeap
1.299 + : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
1.300 + DefStandardHeapTraits<H, CR> > {
1.301 + typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
1.302 + DefStandardHeapTraits<H, CR> >
1.303 + Create;
1.304 + };
1.305 +
1.306 + ///@}
1.307 +
1.308 +
1.309 + /// \brief Constructor.
1.310 + ///
1.311 + /// Constructor of the algorithm.
1.312 + MaxWeightedBipartiteMatching(const BpUGraph& _graph,
1.313 + const WeightMap& _weight)
1.314 + : graph(&_graph), weight(&_weight),
1.315 + anode_matching(_graph), bnode_matching(_graph),
1.316 + anode_potential(_graph), bnode_potential(_graph),
1.317 + _heap_cross_ref(0), local_heap_cross_ref(false),
1.318 + _heap(0), local_heap(0) {}
1.319 +
1.320 + /// \brief Destructor.
1.321 + ///
1.322 + /// Destructor of the algorithm.
1.323 + ~MaxWeightedBipartiteMatching() {
1.324 + destroyStructures();
1.325 + }
1.326 +
1.327 + /// \brief Sets the heap and the cross reference used by algorithm.
1.328 + ///
1.329 + /// Sets the heap and the cross reference used by algorithm.
1.330 + /// If you don't use this function before calling \ref run(),
1.331 + /// it will allocate one. The destuctor deallocates this
1.332 + /// automatically allocated map, of course.
1.333 + /// \return \c (*this)
1.334 + MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
1.335 + if(local_heap_cross_ref) {
1.336 + delete _heap_cross_ref;
1.337 + local_heap_cross_ref = false;
1.338 + }
1.339 + _heap_cross_ref = &crossRef;
1.340 + if(local_heap) {
1.341 + delete _heap;
1.342 + local_heap = false;
1.343 + }
1.344 + _heap = &heap;
1.345 + return *this;
1.346 + }
1.347 +
1.348 + /// \name Execution control
1.349 + /// The simplest way to execute the algorithm is to use
1.350 + /// one of the member functions called \c run().
1.351 + /// \n
1.352 + /// If you need more control on the execution,
1.353 + /// first you must call \ref init() or one alternative for it.
1.354 + /// Finally \ref start() will perform the matching computation or
1.355 + /// with step-by-step execution you can augment the solution.
1.356 +
1.357 + /// @{
1.358 +
1.359 + /// \brief Initalize the data structures.
1.360 + ///
1.361 + /// It initalizes the data structures and creates an empty matching.
1.362 + void init() {
1.363 + initStructures();
1.364 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.365 + anode_matching[it] = INVALID;
1.366 + anode_potential[it] = 0;
1.367 + }
1.368 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.369 + bnode_matching[it] = INVALID;
1.370 + bnode_potential[it] = 0;
1.371 + for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
1.372 + if (-(*weight)[jt] <= bnode_potential[it]) {
1.373 + bnode_potential[it] = -(*weight)[jt];
1.374 + }
1.375 + }
1.376 + }
1.377 + matching_value = 0;
1.378 + matching_size = 0;
1.379 + }
1.380 +
1.381 +
1.382 + /// \brief An augmenting phase of the weighted matching algorithm
1.383 + ///
1.384 + /// It runs an augmenting phase of the weighted matching
1.385 + /// algorithm. The phase finds the best augmenting path and
1.386 + /// augments only on this paths.
1.387 + ///
1.388 + /// The algorithm consists at most
1.389 + /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1.390 + /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1.391 + /// with binary heap.
1.392 + /// \param decrease If the given parameter true the matching value
1.393 + /// can be decreased in the augmenting phase. If we would like
1.394 + /// to calculate the maximum cardinality maximum weighted matching
1.395 + /// then we should let the algorithm to decrease the matching
1.396 + /// value in order to increase the number of the matching edges.
1.397 + bool augment(bool decrease = false) {
1.398 +
1.399 + typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1.400 + typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1.401 +
1.402 + Node bestNode = INVALID;
1.403 + Value bestValue = 0;
1.404 +
1.405 + _heap->clear();
1.406 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.407 + (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1.408 + }
1.409 +
1.410 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.411 + if (anode_matching[it] == INVALID) {
1.412 + _heap->push(it, 0);
1.413 + }
1.414 + }
1.415 +
1.416 + Value bdistMax = 0;
1.417 + while (!_heap->empty()) {
1.418 + Node anode = _heap->top();
1.419 + Value avalue = _heap->prio();
1.420 + _heap->pop();
1.421 + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1.422 + if (jt == anode_matching[anode]) continue;
1.423 + Node bnode = graph->bNode(jt);
1.424 + Value bvalue = avalue + anode_potential[anode] -
1.425 + bnode_potential[bnode] - (*weight)[jt];
1.426 + if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1.427 + bdist[bnode] = bvalue;
1.428 + bpred[bnode] = jt;
1.429 + }
1.430 + if (bvalue > bdistMax) {
1.431 + bdistMax = bvalue;
1.432 + }
1.433 + if (bnode_matching[bnode] != INVALID) {
1.434 + Node newanode = graph->aNode(bnode_matching[bnode]);
1.435 + switch (_heap->state(newanode)) {
1.436 + case Heap::PRE_HEAP:
1.437 + _heap->push(newanode, bvalue);
1.438 + break;
1.439 + case Heap::IN_HEAP:
1.440 + if (bvalue < (*_heap)[newanode]) {
1.441 + _heap->decrease(newanode, bvalue);
1.442 + }
1.443 + break;
1.444 + case Heap::POST_HEAP:
1.445 + break;
1.446 + }
1.447 + } else {
1.448 + if (bestNode == INVALID ||
1.449 + - bvalue - bnode_potential[bnode] > bestValue) {
1.450 + bestValue = - bvalue - bnode_potential[bnode];
1.451 + bestNode = bnode;
1.452 + }
1.453 + }
1.454 + }
1.455 + }
1.456 +
1.457 + if (bestNode == INVALID || (!decrease && bestValue < 0)) {
1.458 + return false;
1.459 + }
1.460 +
1.461 + matching_value += bestValue;
1.462 + ++matching_size;
1.463 +
1.464 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.465 + if (bpred[it] != INVALID) {
1.466 + bnode_potential[it] += bdist[it];
1.467 + } else {
1.468 + bnode_potential[it] += bdistMax;
1.469 + }
1.470 + }
1.471 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.472 + if (anode_matching[it] != INVALID) {
1.473 + Node bnode = graph->bNode(anode_matching[it]);
1.474 + if (bpred[bnode] != INVALID) {
1.475 + anode_potential[it] += bdist[bnode];
1.476 + } else {
1.477 + anode_potential[it] += bdistMax;
1.478 + }
1.479 + }
1.480 + }
1.481 +
1.482 + while (bestNode != INVALID) {
1.483 + UEdge uedge = bpred[bestNode];
1.484 + Node anode = graph->aNode(uedge);
1.485 +
1.486 + bnode_matching[bestNode] = uedge;
1.487 + if (anode_matching[anode] != INVALID) {
1.488 + bestNode = graph->bNode(anode_matching[anode]);
1.489 + } else {
1.490 + bestNode = INVALID;
1.491 + }
1.492 + anode_matching[anode] = uedge;
1.493 + }
1.494 +
1.495 +
1.496 + return true;
1.497 + }
1.498 +
1.499 + /// \brief Starts the algorithm.
1.500 + ///
1.501 + /// Starts the algorithm. It runs augmenting phases until the
1.502 + /// optimal solution reached.
1.503 + ///
1.504 + /// \param maxCardinality If the given value is true it will
1.505 + /// calculate the maximum cardinality maximum matching instead of
1.506 + /// the maximum matching.
1.507 + void start(bool maxCardinality = false) {
1.508 + while (augment(maxCardinality)) {}
1.509 + }
1.510 +
1.511 + /// \brief Runs the algorithm.
1.512 + ///
1.513 + /// It just initalize the algorithm and then start it.
1.514 + ///
1.515 + /// \param maxCardinality If the given value is true it will
1.516 + /// calculate the maximum cardinality maximum matching instead of
1.517 + /// the maximum matching.
1.518 + void run(bool maxCardinality = false) {
1.519 + init();
1.520 + start(maxCardinality);
1.521 + }
1.522 +
1.523 + /// @}
1.524 +
1.525 + /// \name Query Functions
1.526 + /// The result of the %Matching algorithm can be obtained using these
1.527 + /// functions.\n
1.528 + /// Before the use of these functions,
1.529 + /// either run() or start() must be called.
1.530 +
1.531 + ///@{
1.532 +
1.533 + /// \brief Gives back the potential in the NodeMap
1.534 + ///
1.535 + /// Gives back the potential in the NodeMap. The potential
1.536 + /// is feasible if \f$ \pi(a) - \pi(b) - w(ab) = 0 \f$ for
1.537 + /// each matching edges and \f$ \pi(a) - \pi(b) - w(ab) \ge 0 \f$
1.538 + /// for each edges.
1.539 + template <typename PotentialMap>
1.540 + void potential(PotentialMap& potential) {
1.541 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.542 + potential[it] = anode_potential[it];
1.543 + }
1.544 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.545 + potential[it] = bnode_potential[it];
1.546 + }
1.547 + }
1.548 +
1.549 + /// \brief Set true all matching uedge in the map.
1.550 + ///
1.551 + /// Set true all matching uedge in the map. It does not change the
1.552 + /// value mapped to the other uedges.
1.553 + /// \return The number of the matching edges.
1.554 + template <typename MatchingMap>
1.555 + int quickMatching(MatchingMap& matching) {
1.556 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.557 + if (anode_matching[it] != INVALID) {
1.558 + matching[anode_matching[it]] = true;
1.559 + }
1.560 + }
1.561 + return matching_size;
1.562 + }
1.563 +
1.564 + /// \brief Set true all matching uedge in the map and the others to false.
1.565 + ///
1.566 + /// Set true all matching uedge in the map and the others to false.
1.567 + /// \return The number of the matching edges.
1.568 + template <typename MatchingMap>
1.569 + int matching(MatchingMap& matching) {
1.570 + for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.571 + matching[it] = it == anode_matching[graph->aNode(it)];
1.572 + }
1.573 + return matching_size;
1.574 + }
1.575 +
1.576 +
1.577 + /// \brief Return true if the given uedge is in the matching.
1.578 + ///
1.579 + /// It returns true if the given uedge is in the matching.
1.580 + bool matchingEdge(const UEdge& edge) {
1.581 + return anode_matching[graph->aNode(edge)] == edge;
1.582 + }
1.583 +
1.584 + /// \brief Returns the matching edge from the node.
1.585 + ///
1.586 + /// Returns the matching edge from the node. If there is not such
1.587 + /// edge it gives back \c INVALID.
1.588 + UEdge matchingEdge(const Node& node) {
1.589 + if (graph->aNode(node)) {
1.590 + return anode_matching[node];
1.591 + } else {
1.592 + return bnode_matching[node];
1.593 + }
1.594 + }
1.595 +
1.596 + /// \brief Gives back the sum of weights of the matching edges.
1.597 + ///
1.598 + /// Gives back the sum of weights of the matching edges.
1.599 + Value matchingValue() const {
1.600 + return matching_value;
1.601 + }
1.602 +
1.603 + /// \brief Gives back the number of the matching edges.
1.604 + ///
1.605 + /// Gives back the number of the matching edges.
1.606 + int matchingSize() const {
1.607 + return matching_size;
1.608 + }
1.609 +
1.610 + /// @}
1.611 +
1.612 + private:
1.613 +
1.614 + void initStructures() {
1.615 + if (!_heap_cross_ref) {
1.616 + local_heap_cross_ref = true;
1.617 + _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1.618 + }
1.619 + if (!_heap) {
1.620 + local_heap = true;
1.621 + _heap = Traits::createHeap(*_heap_cross_ref);
1.622 + }
1.623 + }
1.624 +
1.625 + void destroyStructures() {
1.626 + if (local_heap_cross_ref) delete _heap_cross_ref;
1.627 + if (local_heap) delete _heap;
1.628 + }
1.629 +
1.630 +
1.631 + private:
1.632 +
1.633 + const BpUGraph *graph;
1.634 + const WeightMap* weight;
1.635 +
1.636 + ANodeMatchingMap anode_matching;
1.637 + BNodeMatchingMap bnode_matching;
1.638 +
1.639 + ANodePotentialMap anode_potential;
1.640 + BNodePotentialMap bnode_potential;
1.641 +
1.642 + Value matching_value;
1.643 + int matching_size;
1.644 +
1.645 + HeapCrossRef *_heap_cross_ref;
1.646 + bool local_heap_cross_ref;
1.647 +
1.648 + Heap *_heap;
1.649 + bool local_heap;
1.650 +
1.651 + };
1.652 +
1.653 + /// \brief Default traits class for minimum cost bipartite matching
1.654 + /// algoritms.
1.655 + ///
1.656 + /// Default traits class for minimum cost bipartite matching
1.657 + /// algoritms.
1.658 + ///
1.659 + /// \param _BpUGraph The bipartite undirected graph
1.660 + /// type.
1.661 + ///
1.662 + /// \param _CostMap Type of cost map.
1.663 + template <typename _BpUGraph, typename _CostMap>
1.664 + struct MinCostMaxBipartiteMatchingDefaultTraits {
1.665 + /// \brief The type of the cost of the undirected edges.
1.666 + typedef typename _CostMap::Value Value;
1.667 +
1.668 + /// The undirected bipartite graph type the algorithm runs on.
1.669 + typedef _BpUGraph BpUGraph;
1.670 +
1.671 + /// The map of the edges costs
1.672 + typedef _CostMap CostMap;
1.673 +
1.674 + /// \brief The cross reference type used by heap.
1.675 + ///
1.676 + /// The cross reference type used by heap.
1.677 + /// Usually it is \c Graph::NodeMap<int>.
1.678 + typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
1.679 +
1.680 + /// \brief Instantiates a HeapCrossRef.
1.681 + ///
1.682 + /// This function instantiates a \ref HeapCrossRef.
1.683 + /// \param graph is the graph, to which we would like to define the
1.684 + /// HeapCrossRef.
1.685 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1.686 + return new HeapCrossRef(graph);
1.687 + }
1.688 +
1.689 + /// \brief The heap type used by costed matching algorithms.
1.690 + ///
1.691 + /// The heap type used by costed matching algorithms. It should
1.692 + /// minimize the priorities and the heap's key type is the graph's
1.693 + /// anode graph's node.
1.694 + ///
1.695 + /// \sa BinHeap
1.696 + typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
1.697 +
1.698 + /// \brief Instantiates a Heap.
1.699 + ///
1.700 + /// This function instantiates a \ref Heap.
1.701 + /// \param crossref The cross reference of the heap.
1.702 + static Heap *createHeap(HeapCrossRef& crossref) {
1.703 + return new Heap(crossref);
1.704 + }
1.705 +
1.706 + };
1.707 +
1.708 +
1.709 + /// \ingroup matching
1.710 + ///
1.711 + /// \brief Bipartite Min Cost Matching algorithm
1.712 + ///
1.713 + /// This class implements the bipartite Min Cost Matching algorithm.
1.714 + /// It uses the successive shortest path algorithm to calculate the
1.715 + /// minimum cost maximum matching in the bipartite graph. The time
1.716 + /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
1.717 + /// default binary heap implementation but this can be improved to
1.718 + /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
1.719 + ///
1.720 + /// The algorithm also provides a potential function on the nodes
1.721 + /// which a dual solution of the matching algorithm and it can be
1.722 + /// used to proof the optimality of the given pimal solution.
1.723 +#ifdef DOXYGEN
1.724 + template <typename _BpUGraph, typename _CostMap, typename _Traits>
1.725 +#else
1.726 + template <typename _BpUGraph,
1.727 + typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
1.728 + typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
1.729 +#endif
1.730 + class MinCostMaxBipartiteMatching {
1.731 + public:
1.732 +
1.733 + typedef _Traits Traits;
1.734 + typedef typename Traits::BpUGraph BpUGraph;
1.735 + typedef typename Traits::CostMap CostMap;
1.736 + typedef typename Traits::Value Value;
1.737 +
1.738 + protected:
1.739 +
1.740 + typedef typename Traits::HeapCrossRef HeapCrossRef;
1.741 + typedef typename Traits::Heap Heap;
1.742 +
1.743 +
1.744 + typedef typename BpUGraph::Node Node;
1.745 + typedef typename BpUGraph::ANodeIt ANodeIt;
1.746 + typedef typename BpUGraph::BNodeIt BNodeIt;
1.747 + typedef typename BpUGraph::UEdge UEdge;
1.748 + typedef typename BpUGraph::UEdgeIt UEdgeIt;
1.749 + typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
1.750 +
1.751 + typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
1.752 + typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
1.753 +
1.754 + typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
1.755 + typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
1.756 +
1.757 +
1.758 + public:
1.759 +
1.760 + /// \brief \ref Exception for uninitialized parameters.
1.761 + ///
1.762 + /// This error represents problems in the initialization
1.763 + /// of the parameters of the algorithms.
1.764 + class UninitializedParameter : public lemon::UninitializedParameter {
1.765 + public:
1.766 + virtual const char* exceptionName() const {
1.767 + return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
1.768 + }
1.769 + };
1.770 +
1.771 + ///\name Named template parameters
1.772 +
1.773 + ///@{
1.774 +
1.775 + template <class H, class CR>
1.776 + struct DefHeapTraits : public Traits {
1.777 + typedef CR HeapCrossRef;
1.778 + typedef H Heap;
1.779 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
1.780 + throw UninitializedParameter();
1.781 + }
1.782 + static Heap *createHeap(HeapCrossRef &) {
1.783 + throw UninitializedParameter();
1.784 + }
1.785 + };
1.786 +
1.787 + /// \brief \ref named-templ-param "Named parameter" for setting heap
1.788 + /// and cross reference type
1.789 + ///
1.790 + /// \ref named-templ-param "Named parameter" for setting heap and cross
1.791 + /// reference type
1.792 + template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1.793 + struct DefHeap
1.794 + : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1.795 + DefHeapTraits<H, CR> > {
1.796 + typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1.797 + DefHeapTraits<H, CR> > Create;
1.798 + };
1.799 +
1.800 + template <class H, class CR>
1.801 + struct DefStandardHeapTraits : public Traits {
1.802 + typedef CR HeapCrossRef;
1.803 + typedef H Heap;
1.804 + static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
1.805 + return new HeapCrossRef(graph);
1.806 + }
1.807 + static Heap *createHeap(HeapCrossRef &crossref) {
1.808 + return new Heap(crossref);
1.809 + }
1.810 + };
1.811 +
1.812 + /// \brief \ref named-templ-param "Named parameter" for setting heap and
1.813 + /// cross reference type with automatic allocation
1.814 + ///
1.815 + /// \ref named-templ-param "Named parameter" for setting heap and cross
1.816 + /// reference type. It can allocate the heap and the cross reference
1.817 + /// object if the cross reference's constructor waits for the graph as
1.818 + /// parameter and the heap's constructor waits for the cross reference.
1.819 + template <class H, class CR = typename BpUGraph::template NodeMap<int> >
1.820 + struct DefStandardHeap
1.821 + : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1.822 + DefStandardHeapTraits<H, CR> > {
1.823 + typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
1.824 + DefStandardHeapTraits<H, CR> >
1.825 + Create;
1.826 + };
1.827 +
1.828 + ///@}
1.829 +
1.830 +
1.831 + /// \brief Constructor.
1.832 + ///
1.833 + /// Constructor of the algorithm.
1.834 + MinCostMaxBipartiteMatching(const BpUGraph& _graph,
1.835 + const CostMap& _cost)
1.836 + : graph(&_graph), cost(&_cost),
1.837 + anode_matching(_graph), bnode_matching(_graph),
1.838 + anode_potential(_graph), bnode_potential(_graph),
1.839 + _heap_cross_ref(0), local_heap_cross_ref(false),
1.840 + _heap(0), local_heap(0) {}
1.841 +
1.842 + /// \brief Destructor.
1.843 + ///
1.844 + /// Destructor of the algorithm.
1.845 + ~MinCostMaxBipartiteMatching() {
1.846 + destroyStructures();
1.847 + }
1.848 +
1.849 + /// \brief Sets the heap and the cross reference used by algorithm.
1.850 + ///
1.851 + /// Sets the heap and the cross reference used by algorithm.
1.852 + /// If you don't use this function before calling \ref run(),
1.853 + /// it will allocate one. The destuctor deallocates this
1.854 + /// automatically allocated map, of course.
1.855 + /// \return \c (*this)
1.856 + MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
1.857 + if(local_heap_cross_ref) {
1.858 + delete _heap_cross_ref;
1.859 + local_heap_cross_ref = false;
1.860 + }
1.861 + _heap_cross_ref = &crossRef;
1.862 + if(local_heap) {
1.863 + delete _heap;
1.864 + local_heap = false;
1.865 + }
1.866 + _heap = &heap;
1.867 + return *this;
1.868 + }
1.869 +
1.870 + /// \name Execution control
1.871 + /// The simplest way to execute the algorithm is to use
1.872 + /// one of the member functions called \c run().
1.873 + /// \n
1.874 + /// If you need more control on the execution,
1.875 + /// first you must call \ref init() or one alternative for it.
1.876 + /// Finally \ref start() will perform the matching computation or
1.877 + /// with step-by-step execution you can augment the solution.
1.878 +
1.879 + /// @{
1.880 +
1.881 + /// \brief Initalize the data structures.
1.882 + ///
1.883 + /// It initalizes the data structures and creates an empty matching.
1.884 + void init() {
1.885 + initStructures();
1.886 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.887 + anode_matching[it] = INVALID;
1.888 + anode_potential[it] = 0;
1.889 + }
1.890 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.891 + bnode_matching[it] = INVALID;
1.892 + bnode_potential[it] = 0;
1.893 + }
1.894 + matching_cost = 0;
1.895 + matching_size = 0;
1.896 + }
1.897 +
1.898 +
1.899 + /// \brief An augmenting phase of the costed matching algorithm
1.900 + ///
1.901 + /// It runs an augmenting phase of the matching algorithm. The
1.902 + /// phase finds the best augmenting path and augments only on this
1.903 + /// paths.
1.904 + ///
1.905 + /// The algorithm consists at most
1.906 + /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
1.907 + /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
1.908 + /// with binary heap.
1.909 + bool augment() {
1.910 +
1.911 + typename BpUGraph::template BNodeMap<Value> bdist(*graph);
1.912 + typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
1.913 +
1.914 + Node bestNode = INVALID;
1.915 + Value bestValue = 0;
1.916 +
1.917 + _heap->clear();
1.918 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.919 + (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
1.920 + }
1.921 +
1.922 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.923 + if (anode_matching[it] == INVALID) {
1.924 + _heap->push(it, 0);
1.925 + }
1.926 + }
1.927 +
1.928 + while (!_heap->empty()) {
1.929 + Node anode = _heap->top();
1.930 + Value avalue = _heap->prio();
1.931 + _heap->pop();
1.932 + for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
1.933 + if (jt == anode_matching[anode]) continue;
1.934 + Node bnode = graph->bNode(jt);
1.935 + Value bvalue = avalue + (*cost)[jt] +
1.936 + anode_potential[anode] - bnode_potential[bnode];
1.937 + if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
1.938 + bdist[bnode] = bvalue;
1.939 + bpred[bnode] = jt;
1.940 + }
1.941 + if (bnode_matching[bnode] != INVALID) {
1.942 + Node newanode = graph->aNode(bnode_matching[bnode]);
1.943 + switch (_heap->state(newanode)) {
1.944 + case Heap::PRE_HEAP:
1.945 + _heap->push(newanode, bvalue);
1.946 + break;
1.947 + case Heap::IN_HEAP:
1.948 + if (bvalue < (*_heap)[newanode]) {
1.949 + _heap->decrease(newanode, bvalue);
1.950 + }
1.951 + break;
1.952 + case Heap::POST_HEAP:
1.953 + break;
1.954 + }
1.955 + } else {
1.956 + if (bestNode == INVALID ||
1.957 + bvalue + bnode_potential[bnode] < bestValue) {
1.958 + bestValue = bvalue + bnode_potential[bnode];
1.959 + bestNode = bnode;
1.960 + }
1.961 + }
1.962 + }
1.963 + }
1.964 +
1.965 + if (bestNode == INVALID) {
1.966 + return false;
1.967 + }
1.968 +
1.969 + matching_cost += bestValue;
1.970 + ++matching_size;
1.971 +
1.972 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.973 + if (bpred[it] != INVALID) {
1.974 + bnode_potential[it] += bdist[it];
1.975 + }
1.976 + }
1.977 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.978 + if (anode_matching[it] != INVALID) {
1.979 + Node bnode = graph->bNode(anode_matching[it]);
1.980 + if (bpred[bnode] != INVALID) {
1.981 + anode_potential[it] += bdist[bnode];
1.982 + }
1.983 + }
1.984 + }
1.985 +
1.986 + while (bestNode != INVALID) {
1.987 + UEdge uedge = bpred[bestNode];
1.988 + Node anode = graph->aNode(uedge);
1.989 +
1.990 + bnode_matching[bestNode] = uedge;
1.991 + if (anode_matching[anode] != INVALID) {
1.992 + bestNode = graph->bNode(anode_matching[anode]);
1.993 + } else {
1.994 + bestNode = INVALID;
1.995 + }
1.996 + anode_matching[anode] = uedge;
1.997 + }
1.998 +
1.999 +
1.1000 + return true;
1.1001 + }
1.1002 +
1.1003 + /// \brief Starts the algorithm.
1.1004 + ///
1.1005 + /// Starts the algorithm. It runs augmenting phases until the
1.1006 + /// optimal solution reached.
1.1007 + void start() {
1.1008 + while (augment()) {}
1.1009 + }
1.1010 +
1.1011 + /// \brief Runs the algorithm.
1.1012 + ///
1.1013 + /// It just initalize the algorithm and then start it.
1.1014 + void run() {
1.1015 + init();
1.1016 + start();
1.1017 + }
1.1018 +
1.1019 + /// @}
1.1020 +
1.1021 + /// \name Query Functions
1.1022 + /// The result of the %Matching algorithm can be obtained using these
1.1023 + /// functions.\n
1.1024 + /// Before the use of these functions,
1.1025 + /// either run() or start() must be called.
1.1026 +
1.1027 + ///@{
1.1028 +
1.1029 + /// \brief Gives back the potential in the NodeMap
1.1030 + ///
1.1031 + /// Gives back the potential in the NodeMap. The potential
1.1032 + /// is feasible if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
1.1033 + /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
1.1034 + /// for each edges.
1.1035 + template <typename PotentialMap>
1.1036 + void potential(PotentialMap& potential) {
1.1037 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.1038 + potential[it] = anode_potential[it];
1.1039 + }
1.1040 + for (BNodeIt it(*graph); it != INVALID; ++it) {
1.1041 + potential[it] = bnode_potential[it];
1.1042 + }
1.1043 + }
1.1044 +
1.1045 + /// \brief Set true all matching uedge in the map.
1.1046 + ///
1.1047 + /// Set true all matching uedge in the map. It does not change the
1.1048 + /// value mapped to the other uedges.
1.1049 + /// \return The number of the matching edges.
1.1050 + template <typename MatchingMap>
1.1051 + int quickMatching(MatchingMap& matching) {
1.1052 + for (ANodeIt it(*graph); it != INVALID; ++it) {
1.1053 + if (anode_matching[it] != INVALID) {
1.1054 + matching[anode_matching[it]] = true;
1.1055 + }
1.1056 + }
1.1057 + return matching_size;
1.1058 + }
1.1059 +
1.1060 + /// \brief Set true all matching uedge in the map and the others to false.
1.1061 + ///
1.1062 + /// Set true all matching uedge in the map and the others to false.
1.1063 + /// \return The number of the matching edges.
1.1064 + template <typename MatchingMap>
1.1065 + int matching(MatchingMap& matching) {
1.1066 + for (UEdgeIt it(*graph); it != INVALID; ++it) {
1.1067 + matching[it] = it == anode_matching[graph->aNode(it)];
1.1068 + }
1.1069 + return matching_size;
1.1070 + }
1.1071 +
1.1072 +
1.1073 + /// \brief Return true if the given uedge is in the matching.
1.1074 + ///
1.1075 + /// It returns true if the given uedge is in the matching.
1.1076 + bool matchingEdge(const UEdge& edge) {
1.1077 + return anode_matching[graph->aNode(edge)] == edge;
1.1078 + }
1.1079 +
1.1080 + /// \brief Returns the matching edge from the node.
1.1081 + ///
1.1082 + /// Returns the matching edge from the node. If there is not such
1.1083 + /// edge it gives back \c INVALID.
1.1084 + UEdge matchingEdge(const Node& node) {
1.1085 + if (graph->aNode(node)) {
1.1086 + return anode_matching[node];
1.1087 + } else {
1.1088 + return bnode_matching[node];
1.1089 + }
1.1090 + }
1.1091 +
1.1092 + /// \brief Gives back the sum of costs of the matching edges.
1.1093 + ///
1.1094 + /// Gives back the sum of costs of the matching edges.
1.1095 + Value matchingCost() const {
1.1096 + return matching_cost;
1.1097 + }
1.1098 +
1.1099 + /// \brief Gives back the number of the matching edges.
1.1100 + ///
1.1101 + /// Gives back the number of the matching edges.
1.1102 + int matchingSize() const {
1.1103 + return matching_size;
1.1104 + }
1.1105 +
1.1106 + /// @}
1.1107 +
1.1108 + private:
1.1109 +
1.1110 + void initStructures() {
1.1111 + if (!_heap_cross_ref) {
1.1112 + local_heap_cross_ref = true;
1.1113 + _heap_cross_ref = Traits::createHeapCrossRef(*graph);
1.1114 + }
1.1115 + if (!_heap) {
1.1116 + local_heap = true;
1.1117 + _heap = Traits::createHeap(*_heap_cross_ref);
1.1118 + }
1.1119 + }
1.1120 +
1.1121 + void destroyStructures() {
1.1122 + if (local_heap_cross_ref) delete _heap_cross_ref;
1.1123 + if (local_heap) delete _heap;
1.1124 + }
1.1125 +
1.1126 +
1.1127 + private:
1.1128 +
1.1129 + const BpUGraph *graph;
1.1130 + const CostMap* cost;
1.1131 +
1.1132 + ANodeMatchingMap anode_matching;
1.1133 + BNodeMatchingMap bnode_matching;
1.1134 +
1.1135 + ANodePotentialMap anode_potential;
1.1136 + BNodePotentialMap bnode_potential;
1.1137 +
1.1138 + Value matching_cost;
1.1139 + int matching_size;
1.1140 +
1.1141 + HeapCrossRef *_heap_cross_ref;
1.1142 + bool local_heap_cross_ref;
1.1143 +
1.1144 + Heap *_heap;
1.1145 + bool local_heap;
1.1146
1.1147 };
1.1148
2.1 --- a/test/bipartite_matching_test.cc Fri Apr 14 18:05:02 2006 +0000
2.2 +++ b/test/bipartite_matching_test.cc Fri Apr 14 18:07:33 2006 +0000
2.3 @@ -9,6 +9,8 @@
2.4
2.5 #include <lemon/graph_utils.h>
2.6
2.7 +#include <lemon/maps.h>
2.8 +
2.9 #include "test_tools.h"
2.10
2.11 using namespace std;
2.12 @@ -24,6 +26,13 @@
2.13 int n = argc > 1 ? atoi(argv[1]) : 100;
2.14 int m = argc > 2 ? atoi(argv[2]) : 100;
2.15 int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m));
2.16 + int c = argc > 4 ? atoi(argv[4]) : 100;
2.17 +
2.18 + Graph::UEdgeMap<int> weight(graph);
2.19 +
2.20 + int max_cardinality;
2.21 + int max_weight;
2.22 + int max_cardinality_max_weight;
2.23
2.24 for (int i = 0; i < n; ++i) {
2.25 Node node = graph.addANode();
2.26 @@ -36,7 +45,8 @@
2.27 for (int i = 0; i < e; ++i) {
2.28 Node aNode = aNodes[urandom(n)];
2.29 Node bNode = bNodes[urandom(m)];
2.30 - graph.addEdge(aNode, bNode);
2.31 + UEdge uedge = graph.addEdge(aNode, bNode);
2.32 + weight[uedge] = urandom(c);
2.33 }
2.34
2.35 {
2.36 @@ -57,9 +67,10 @@
2.37 int num = 0;
2.38 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
2.39 if (mm[jt]) ++num;
2.40 - }
2.41 + }
2.42 check(num <= 1, "INVALID PRIMAL");
2.43 }
2.44 + max_cardinality = bpmatch.matchingSize();
2.45 }
2.46
2.47 {
2.48 @@ -111,20 +122,63 @@
2.49 }
2.50
2.51 {
2.52 - SwapBpUGraphAdaptor<Graph> swapgraph(graph);
2.53 - MaxBipartiteMatching<SwapBpUGraphAdaptor<Graph> > bpmatch(swapgraph);
2.54 + MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
2.55 +
2.56 + bpmatch.init();
2.57 +
2.58 + max_weight = 0;
2.59 + while (bpmatch.augment(true)) {
2.60 +
2.61 + Graph::UEdgeMap<bool> mm(graph);
2.62 + Graph::NodeMap<int> pm(graph);
2.63 +
2.64 + bpmatch.matching(mm);
2.65 + bpmatch.potential(pm);
2.66 +
2.67 + for (UEdgeIt it(graph); it != INVALID; ++it) {
2.68 + if (mm[it]) {
2.69 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
2.70 + "INVALID DUAL");
2.71 + } else {
2.72 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
2.73 + "INVALID DUAL");
2.74 + }
2.75 + }
2.76 +
2.77 + for (ANodeIt it(graph); it != INVALID; ++it) {
2.78 + int num = 0;
2.79 + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
2.80 + if (mm[jt]) ++num;
2.81 + }
2.82 + check(num <= 1, "INVALID PRIMAL");
2.83 + }
2.84 + if (bpmatch.matchingValue() > max_weight) {
2.85 + max_weight = bpmatch.matchingValue();
2.86 + }
2.87 + }
2.88 + }
2.89 +
2.90 + {
2.91 + MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
2.92
2.93 bpmatch.run();
2.94 +
2.95 + Graph::UEdgeMap<bool> mm(graph);
2.96 + Graph::NodeMap<int> pm(graph);
2.97
2.98 - Graph::UEdgeMap<bool> mm(graph);
2.99 - Graph::NodeMap<bool> cs(graph);
2.100 -
2.101 - check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL");
2.102 + bpmatch.matching(mm);
2.103 + bpmatch.potential(pm);
2.104
2.105 for (UEdgeIt it(graph); it != INVALID; ++it) {
2.106 - check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL");
2.107 + if (mm[it]) {
2.108 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
2.109 + "INVALID DUAL");
2.110 + } else {
2.111 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
2.112 + "INVALID DUAL");
2.113 + }
2.114 }
2.115 -
2.116 +
2.117 for (ANodeIt it(graph); it != INVALID; ++it) {
2.118 int num = 0;
2.119 for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
2.120 @@ -132,6 +186,79 @@
2.121 }
2.122 check(num <= 1, "INVALID PRIMAL");
2.123 }
2.124 + check(max_weight == bpmatch.matchingValue(), "WRONG WEIGHT");
2.125 + }
2.126 +
2.127 + {
2.128 + MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
2.129 +
2.130 + bpmatch.run(true);
2.131 +
2.132 + Graph::UEdgeMap<bool> mm(graph);
2.133 + Graph::NodeMap<int> pm(graph);
2.134 +
2.135 + bpmatch.matching(mm);
2.136 + bpmatch.potential(pm);
2.137 +
2.138 + for (UEdgeIt it(graph); it != INVALID; ++it) {
2.139 + if (mm[it]) {
2.140 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
2.141 + "INVALID DUAL");
2.142 + } else {
2.143 + check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
2.144 + "INVALID DUAL");
2.145 + }
2.146 + }
2.147 +
2.148 + for (ANodeIt it(graph); it != INVALID; ++it) {
2.149 + int num = 0;
2.150 + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
2.151 + if (mm[jt]) ++num;
2.152 + }
2.153 + check(num <= 1, "INVALID PRIMAL");
2.154 + }
2.155 + check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
2.156 + max_cardinality_max_weight = bpmatch.matchingValue();
2.157 +
2.158 + }
2.159 +
2.160 + {
2.161 + Graph::UEdgeMap<int> cost(graph);
2.162 +
2.163 + cost = subMap(constMap<UEdge>(c), weight);
2.164 +
2.165 + MinCostMaxBipartiteMatching<Graph> bpmatch(graph, cost);
2.166 +
2.167 + bpmatch.run();
2.168 +
2.169 + Graph::UEdgeMap<bool> mm(graph);
2.170 + Graph::NodeMap<int> pm(graph);
2.171 +
2.172 + bpmatch.matching(mm);
2.173 + bpmatch.potential(pm);
2.174 +
2.175 + for (UEdgeIt it(graph); it != INVALID; ++it) {
2.176 + if (mm[it]) {
2.177 + check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] == 0,
2.178 + "INVALID DUAL");
2.179 + } else {
2.180 + check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] >= 0,
2.181 + "INVALID DUAL");
2.182 + }
2.183 + }
2.184 +
2.185 + for (ANodeIt it(graph); it != INVALID; ++it) {
2.186 + int num = 0;
2.187 + for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
2.188 + if (mm[jt]) ++num;
2.189 + }
2.190 + check(num <= 1, "INVALID PRIMAL");
2.191 + }
2.192 +
2.193 + check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
2.194 + check(max_cardinality * c - max_cardinality_max_weight
2.195 + == bpmatch.matchingCost(), "WRONG SIZE");
2.196 +
2.197 }
2.198
2.199 return 0;