1.1 --- a/doc/groups.dox Tue Mar 16 21:18:39 2010 +0100
1.2 +++ b/doc/groups.dox Tue Mar 16 21:27:35 2010 +0100
1.3 @@ -386,7 +386,7 @@
1.4 also provide functions to query the minimum cut, which is the dual
1.5 problem of maximum flow.
1.6
1.7 -\ref Circulation is a preflow push-relabel algorithm implemented directly
1.8 +\ref Circulation is a preflow push-relabel algorithm implemented directly
1.9 for finding feasible circulations, which is a somewhat different problem,
1.10 but it is strongly related to maximum flow.
1.11 For more information, see \ref Circulation.
1.12 @@ -522,6 +522,13 @@
1.13 - \ref MaxWeightedPerfectMatching
1.14 Edmond's blossom shrinking algorithm for calculating maximum weighted
1.15 perfect matching in general graphs.
1.16 +- \ref MaxFractionalMatching Push-relabel algorithm for calculating
1.17 + maximum cardinality fractional matching in general graphs.
1.18 +- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
1.19 + maximum weighted fractional matching in general graphs.
1.20 +- \ref MaxWeightedPerfectFractionalMatching
1.21 + Augmenting path algorithm for calculating maximum weighted
1.22 + perfect fractional matching in general graphs.
1.23
1.24 \image html matching.png
1.25 \image latex matching.eps "Min Cost Perfect Matching" width=\textwidth
2.1 --- a/lemon/Makefile.am Tue Mar 16 21:18:39 2010 +0100
2.2 +++ b/lemon/Makefile.am Tue Mar 16 21:27:35 2010 +0100
2.3 @@ -84,6 +84,7 @@
2.4 lemon/error.h \
2.5 lemon/euler.h \
2.6 lemon/fib_heap.h \
2.7 + lemon/fractional_matching.h \
2.8 lemon/full_graph.h \
2.9 lemon/glpk.h \
2.10 lemon/gomory_hu.h \
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/lemon/fractional_matching.h Tue Mar 16 21:27:35 2010 +0100
3.3 @@ -0,0 +1,2130 @@
3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
3.5 + *
3.6 + * This file is a part of LEMON, a generic C++ optimization library.
3.7 + *
3.8 + * Copyright (C) 2003-2009
3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
3.11 + *
3.12 + * Permission to use, modify and distribute this software is granted
3.13 + * provided that this copyright notice appears in all copies. For
3.14 + * precise terms see the accompanying LICENSE file.
3.15 + *
3.16 + * This software is provided "AS IS" with no warranty of any kind,
3.17 + * express or implied, and with no claim as to its suitability for any
3.18 + * purpose.
3.19 + *
3.20 + */
3.21 +
3.22 +#ifndef LEMON_FRACTIONAL_MATCHING_H
3.23 +#define LEMON_FRACTIONAL_MATCHING_H
3.24 +
3.25 +#include <vector>
3.26 +#include <queue>
3.27 +#include <set>
3.28 +#include <limits>
3.29 +
3.30 +#include <lemon/core.h>
3.31 +#include <lemon/unionfind.h>
3.32 +#include <lemon/bin_heap.h>
3.33 +#include <lemon/maps.h>
3.34 +#include <lemon/assert.h>
3.35 +#include <lemon/elevator.h>
3.36 +
3.37 +///\ingroup matching
3.38 +///\file
3.39 +///\brief Fractional matching algorithms in general graphs.
3.40 +
3.41 +namespace lemon {
3.42 +
3.43 + /// \brief Default traits class of MaxFractionalMatching class.
3.44 + ///
3.45 + /// Default traits class of MaxFractionalMatching class.
3.46 + /// \tparam GR Graph type.
3.47 + template <typename GR>
3.48 + struct MaxFractionalMatchingDefaultTraits {
3.49 +
3.50 + /// \brief The type of the graph the algorithm runs on.
3.51 + typedef GR Graph;
3.52 +
3.53 + /// \brief The type of the map that stores the matching.
3.54 + ///
3.55 + /// The type of the map that stores the matching arcs.
3.56 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
3.57 + typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
3.58 +
3.59 + /// \brief Instantiates a MatchingMap.
3.60 + ///
3.61 + /// This function instantiates a \ref MatchingMap.
3.62 + /// \param graph The graph for which we would like to define
3.63 + /// the matching map.
3.64 + static MatchingMap* createMatchingMap(const Graph& graph) {
3.65 + return new MatchingMap(graph);
3.66 + }
3.67 +
3.68 + /// \brief The elevator type used by MaxFractionalMatching algorithm.
3.69 + ///
3.70 + /// The elevator type used by MaxFractionalMatching algorithm.
3.71 + ///
3.72 + /// \sa Elevator
3.73 + /// \sa LinkedElevator
3.74 + typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
3.75 +
3.76 + /// \brief Instantiates an Elevator.
3.77 + ///
3.78 + /// This function instantiates an \ref Elevator.
3.79 + /// \param graph The graph for which we would like to define
3.80 + /// the elevator.
3.81 + /// \param max_level The maximum level of the elevator.
3.82 + static Elevator* createElevator(const Graph& graph, int max_level) {
3.83 + return new Elevator(graph, max_level);
3.84 + }
3.85 + };
3.86 +
3.87 + /// \ingroup matching
3.88 + ///
3.89 + /// \brief Max cardinality fractional matching
3.90 + ///
3.91 + /// This class provides an implementation of fractional matching
3.92 + /// algorithm based on push-relabel principle.
3.93 + ///
3.94 + /// The maximum cardinality fractional matching is a relaxation of the
3.95 + /// maximum cardinality matching problem where the odd set constraints
3.96 + /// are omitted.
3.97 + /// It can be formulated with the following linear program.
3.98 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
3.99 + /// \f[x_e \ge 0\quad \forall e\in E\f]
3.100 + /// \f[\max \sum_{e\in E}x_e\f]
3.101 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.102 + /// \f$X\f$. The result can be represented as the union of a
3.103 + /// matching with one value edges and a set of odd length cycles
3.104 + /// with half value edges.
3.105 + ///
3.106 + /// The algorithm calculates an optimal fractional matching and a
3.107 + /// barrier. The number of adjacents of any node set minus the size
3.108 + /// of node set is a lower bound on the uncovered nodes in the
3.109 + /// graph. For maximum matching a barrier is computed which
3.110 + /// maximizes this difference.
3.111 + ///
3.112 + /// The algorithm can be executed with the run() function. After it
3.113 + /// the matching (the primal solution) and the barrier (the dual
3.114 + /// solution) can be obtained using the query functions.
3.115 + ///
3.116 + /// The primal solution is multiplied by
3.117 + /// \ref MaxFractionalMatching::primalScale "2".
3.118 + ///
3.119 + /// \tparam GR The undirected graph type the algorithm runs on.
3.120 +#ifdef DOXYGEN
3.121 + template <typename GR, typename TR>
3.122 +#else
3.123 + template <typename GR,
3.124 + typename TR = MaxFractionalMatchingDefaultTraits<GR> >
3.125 +#endif
3.126 + class MaxFractionalMatching {
3.127 + public:
3.128 +
3.129 + /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
3.130 + /// class" of the algorithm.
3.131 + typedef TR Traits;
3.132 + /// The type of the graph the algorithm runs on.
3.133 + typedef typename TR::Graph Graph;
3.134 + /// The type of the matching map.
3.135 + typedef typename TR::MatchingMap MatchingMap;
3.136 + /// The type of the elevator.
3.137 + typedef typename TR::Elevator Elevator;
3.138 +
3.139 + /// \brief Scaling factor for primal solution
3.140 + ///
3.141 + /// Scaling factor for primal solution.
3.142 + static const int primalScale = 2;
3.143 +
3.144 + private:
3.145 +
3.146 + const Graph &_graph;
3.147 + int _node_num;
3.148 + bool _allow_loops;
3.149 + int _empty_level;
3.150 +
3.151 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
3.152 +
3.153 + bool _local_matching;
3.154 + MatchingMap *_matching;
3.155 +
3.156 + bool _local_level;
3.157 + Elevator *_level;
3.158 +
3.159 + typedef typename Graph::template NodeMap<int> InDegMap;
3.160 + InDegMap *_indeg;
3.161 +
3.162 + void createStructures() {
3.163 + _node_num = countNodes(_graph);
3.164 +
3.165 + if (!_matching) {
3.166 + _local_matching = true;
3.167 + _matching = Traits::createMatchingMap(_graph);
3.168 + }
3.169 + if (!_level) {
3.170 + _local_level = true;
3.171 + _level = Traits::createElevator(_graph, _node_num);
3.172 + }
3.173 + if (!_indeg) {
3.174 + _indeg = new InDegMap(_graph);
3.175 + }
3.176 + }
3.177 +
3.178 + void destroyStructures() {
3.179 + if (_local_matching) {
3.180 + delete _matching;
3.181 + }
3.182 + if (_local_level) {
3.183 + delete _level;
3.184 + }
3.185 + if (_indeg) {
3.186 + delete _indeg;
3.187 + }
3.188 + }
3.189 +
3.190 + void postprocessing() {
3.191 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.192 + if ((*_indeg)[n] != 0) continue;
3.193 + _indeg->set(n, -1);
3.194 + Node u = n;
3.195 + while ((*_matching)[u] != INVALID) {
3.196 + Node v = _graph.target((*_matching)[u]);
3.197 + _indeg->set(v, -1);
3.198 + Arc a = _graph.oppositeArc((*_matching)[u]);
3.199 + u = _graph.target((*_matching)[v]);
3.200 + _indeg->set(u, -1);
3.201 + _matching->set(v, a);
3.202 + }
3.203 + }
3.204 +
3.205 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.206 + if ((*_indeg)[n] != 1) continue;
3.207 + _indeg->set(n, -1);
3.208 +
3.209 + int num = 1;
3.210 + Node u = _graph.target((*_matching)[n]);
3.211 + while (u != n) {
3.212 + _indeg->set(u, -1);
3.213 + u = _graph.target((*_matching)[u]);
3.214 + ++num;
3.215 + }
3.216 + if (num % 2 == 0 && num > 2) {
3.217 + Arc prev = _graph.oppositeArc((*_matching)[n]);
3.218 + Node v = _graph.target((*_matching)[n]);
3.219 + u = _graph.target((*_matching)[v]);
3.220 + _matching->set(v, prev);
3.221 + while (u != n) {
3.222 + prev = _graph.oppositeArc((*_matching)[u]);
3.223 + v = _graph.target((*_matching)[u]);
3.224 + u = _graph.target((*_matching)[v]);
3.225 + _matching->set(v, prev);
3.226 + }
3.227 + }
3.228 + }
3.229 + }
3.230 +
3.231 + public:
3.232 +
3.233 + typedef MaxFractionalMatching Create;
3.234 +
3.235 + ///\name Named Template Parameters
3.236 +
3.237 + ///@{
3.238 +
3.239 + template <typename T>
3.240 + struct SetMatchingMapTraits : public Traits {
3.241 + typedef T MatchingMap;
3.242 + static MatchingMap *createMatchingMap(const Graph&) {
3.243 + LEMON_ASSERT(false, "MatchingMap is not initialized");
3.244 + return 0; // ignore warnings
3.245 + }
3.246 + };
3.247 +
3.248 + /// \brief \ref named-templ-param "Named parameter" for setting
3.249 + /// MatchingMap type
3.250 + ///
3.251 + /// \ref named-templ-param "Named parameter" for setting MatchingMap
3.252 + /// type.
3.253 + template <typename T>
3.254 + struct SetMatchingMap
3.255 + : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
3.256 + typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
3.257 + };
3.258 +
3.259 + template <typename T>
3.260 + struct SetElevatorTraits : public Traits {
3.261 + typedef T Elevator;
3.262 + static Elevator *createElevator(const Graph&, int) {
3.263 + LEMON_ASSERT(false, "Elevator is not initialized");
3.264 + return 0; // ignore warnings
3.265 + }
3.266 + };
3.267 +
3.268 + /// \brief \ref named-templ-param "Named parameter" for setting
3.269 + /// Elevator type
3.270 + ///
3.271 + /// \ref named-templ-param "Named parameter" for setting Elevator
3.272 + /// type. If this named parameter is used, then an external
3.273 + /// elevator object must be passed to the algorithm using the
3.274 + /// \ref elevator(Elevator&) "elevator()" function before calling
3.275 + /// \ref run() or \ref init().
3.276 + /// \sa SetStandardElevator
3.277 + template <typename T>
3.278 + struct SetElevator
3.279 + : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
3.280 + typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
3.281 + };
3.282 +
3.283 + template <typename T>
3.284 + struct SetStandardElevatorTraits : public Traits {
3.285 + typedef T Elevator;
3.286 + static Elevator *createElevator(const Graph& graph, int max_level) {
3.287 + return new Elevator(graph, max_level);
3.288 + }
3.289 + };
3.290 +
3.291 + /// \brief \ref named-templ-param "Named parameter" for setting
3.292 + /// Elevator type with automatic allocation
3.293 + ///
3.294 + /// \ref named-templ-param "Named parameter" for setting Elevator
3.295 + /// type with automatic allocation.
3.296 + /// The Elevator should have standard constructor interface to be
3.297 + /// able to automatically created by the algorithm (i.e. the
3.298 + /// graph and the maximum level should be passed to it).
3.299 + /// However an external elevator object could also be passed to the
3.300 + /// algorithm with the \ref elevator(Elevator&) "elevator()" function
3.301 + /// before calling \ref run() or \ref init().
3.302 + /// \sa SetElevator
3.303 + template <typename T>
3.304 + struct SetStandardElevator
3.305 + : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
3.306 + typedef MaxFractionalMatching<Graph,
3.307 + SetStandardElevatorTraits<T> > Create;
3.308 + };
3.309 +
3.310 + /// @}
3.311 +
3.312 + protected:
3.313 +
3.314 + MaxFractionalMatching() {}
3.315 +
3.316 + public:
3.317 +
3.318 + /// \brief Constructor
3.319 + ///
3.320 + /// Constructor.
3.321 + ///
3.322 + MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
3.323 + : _graph(graph), _allow_loops(allow_loops),
3.324 + _local_matching(false), _matching(0),
3.325 + _local_level(false), _level(0), _indeg(0)
3.326 + {}
3.327 +
3.328 + ~MaxFractionalMatching() {
3.329 + destroyStructures();
3.330 + }
3.331 +
3.332 + /// \brief Sets the matching map.
3.333 + ///
3.334 + /// Sets the matching map.
3.335 + /// If you don't use this function before calling \ref run() or
3.336 + /// \ref init(), an instance will be allocated automatically.
3.337 + /// The destructor deallocates this automatically allocated map,
3.338 + /// of course.
3.339 + /// \return <tt>(*this)</tt>
3.340 + MaxFractionalMatching& matchingMap(MatchingMap& map) {
3.341 + if (_local_matching) {
3.342 + delete _matching;
3.343 + _local_matching = false;
3.344 + }
3.345 + _matching = ↦
3.346 + return *this;
3.347 + }
3.348 +
3.349 + /// \brief Sets the elevator used by algorithm.
3.350 + ///
3.351 + /// Sets the elevator used by algorithm.
3.352 + /// If you don't use this function before calling \ref run() or
3.353 + /// \ref init(), an instance will be allocated automatically.
3.354 + /// The destructor deallocates this automatically allocated elevator,
3.355 + /// of course.
3.356 + /// \return <tt>(*this)</tt>
3.357 + MaxFractionalMatching& elevator(Elevator& elevator) {
3.358 + if (_local_level) {
3.359 + delete _level;
3.360 + _local_level = false;
3.361 + }
3.362 + _level = &elevator;
3.363 + return *this;
3.364 + }
3.365 +
3.366 + /// \brief Returns a const reference to the elevator.
3.367 + ///
3.368 + /// Returns a const reference to the elevator.
3.369 + ///
3.370 + /// \pre Either \ref run() or \ref init() must be called before
3.371 + /// using this function.
3.372 + const Elevator& elevator() const {
3.373 + return *_level;
3.374 + }
3.375 +
3.376 + /// \name Execution control
3.377 + /// The simplest way to execute the algorithm is to use one of the
3.378 + /// member functions called \c run(). \n
3.379 + /// If you need more control on the execution, first
3.380 + /// you must call \ref init() and then one variant of the start()
3.381 + /// member.
3.382 +
3.383 + /// @{
3.384 +
3.385 + /// \brief Initializes the internal data structures.
3.386 + ///
3.387 + /// Initializes the internal data structures and sets the initial
3.388 + /// matching.
3.389 + void init() {
3.390 + createStructures();
3.391 +
3.392 + _level->initStart();
3.393 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.394 + _indeg->set(n, 0);
3.395 + _matching->set(n, INVALID);
3.396 + _level->initAddItem(n);
3.397 + }
3.398 + _level->initFinish();
3.399 +
3.400 + _empty_level = _node_num;
3.401 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.402 + for (OutArcIt a(_graph, n); a != INVALID; ++a) {
3.403 + if (_graph.target(a) == n && !_allow_loops) continue;
3.404 + _matching->set(n, a);
3.405 + Node v = _graph.target((*_matching)[n]);
3.406 + _indeg->set(v, (*_indeg)[v] + 1);
3.407 + break;
3.408 + }
3.409 + }
3.410 +
3.411 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.412 + if ((*_indeg)[n] == 0) {
3.413 + _level->activate(n);
3.414 + }
3.415 + }
3.416 + }
3.417 +
3.418 + /// \brief Starts the algorithm and computes a fractional matching
3.419 + ///
3.420 + /// The algorithm computes a maximum fractional matching.
3.421 + ///
3.422 + /// \param postprocess The algorithm computes first a matching
3.423 + /// which is a union of a matching with one value edges, cycles
3.424 + /// with half value edges and even length paths with half value
3.425 + /// edges. If the parameter is true, then after the push-relabel
3.426 + /// algorithm it postprocesses the matching to contain only
3.427 + /// matching edges and half value odd cycles.
3.428 + void start(bool postprocess = true) {
3.429 + Node n;
3.430 + while ((n = _level->highestActive()) != INVALID) {
3.431 + int level = _level->highestActiveLevel();
3.432 + int new_level = _level->maxLevel();
3.433 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
3.434 + Node u = _graph.source(a);
3.435 + if (n == u && !_allow_loops) continue;
3.436 + Node v = _graph.target((*_matching)[u]);
3.437 + if ((*_level)[v] < level) {
3.438 + _indeg->set(v, (*_indeg)[v] - 1);
3.439 + if ((*_indeg)[v] == 0) {
3.440 + _level->activate(v);
3.441 + }
3.442 + _matching->set(u, a);
3.443 + _indeg->set(n, (*_indeg)[n] + 1);
3.444 + _level->deactivate(n);
3.445 + goto no_more_push;
3.446 + } else if (new_level > (*_level)[v]) {
3.447 + new_level = (*_level)[v];
3.448 + }
3.449 + }
3.450 +
3.451 + if (new_level + 1 < _level->maxLevel()) {
3.452 + _level->liftHighestActive(new_level + 1);
3.453 + } else {
3.454 + _level->liftHighestActiveToTop();
3.455 + }
3.456 + if (_level->emptyLevel(level)) {
3.457 + _level->liftToTop(level);
3.458 + }
3.459 + no_more_push:
3.460 + ;
3.461 + }
3.462 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.463 + if ((*_matching)[n] == INVALID) continue;
3.464 + Node u = _graph.target((*_matching)[n]);
3.465 + if ((*_indeg)[u] > 1) {
3.466 + _indeg->set(u, (*_indeg)[u] - 1);
3.467 + _matching->set(n, INVALID);
3.468 + }
3.469 + }
3.470 + if (postprocess) {
3.471 + postprocessing();
3.472 + }
3.473 + }
3.474 +
3.475 + /// \brief Starts the algorithm and computes a perfect fractional
3.476 + /// matching
3.477 + ///
3.478 + /// The algorithm computes a perfect fractional matching. If it
3.479 + /// does not exists, then the algorithm returns false and the
3.480 + /// matching is undefined and the barrier.
3.481 + ///
3.482 + /// \param postprocess The algorithm computes first a matching
3.483 + /// which is a union of a matching with one value edges, cycles
3.484 + /// with half value edges and even length paths with half value
3.485 + /// edges. If the parameter is true, then after the push-relabel
3.486 + /// algorithm it postprocesses the matching to contain only
3.487 + /// matching edges and half value odd cycles.
3.488 + bool startPerfect(bool postprocess = true) {
3.489 + Node n;
3.490 + while ((n = _level->highestActive()) != INVALID) {
3.491 + int level = _level->highestActiveLevel();
3.492 + int new_level = _level->maxLevel();
3.493 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
3.494 + Node u = _graph.source(a);
3.495 + if (n == u && !_allow_loops) continue;
3.496 + Node v = _graph.target((*_matching)[u]);
3.497 + if ((*_level)[v] < level) {
3.498 + _indeg->set(v, (*_indeg)[v] - 1);
3.499 + if ((*_indeg)[v] == 0) {
3.500 + _level->activate(v);
3.501 + }
3.502 + _matching->set(u, a);
3.503 + _indeg->set(n, (*_indeg)[n] + 1);
3.504 + _level->deactivate(n);
3.505 + goto no_more_push;
3.506 + } else if (new_level > (*_level)[v]) {
3.507 + new_level = (*_level)[v];
3.508 + }
3.509 + }
3.510 +
3.511 + if (new_level + 1 < _level->maxLevel()) {
3.512 + _level->liftHighestActive(new_level + 1);
3.513 + } else {
3.514 + _level->liftHighestActiveToTop();
3.515 + _empty_level = _level->maxLevel() - 1;
3.516 + return false;
3.517 + }
3.518 + if (_level->emptyLevel(level)) {
3.519 + _level->liftToTop(level);
3.520 + _empty_level = level;
3.521 + return false;
3.522 + }
3.523 + no_more_push:
3.524 + ;
3.525 + }
3.526 + if (postprocess) {
3.527 + postprocessing();
3.528 + }
3.529 + return true;
3.530 + }
3.531 +
3.532 + /// \brief Runs the algorithm
3.533 + ///
3.534 + /// Just a shortcut for the next code:
3.535 + ///\code
3.536 + /// init();
3.537 + /// start();
3.538 + ///\endcode
3.539 + void run(bool postprocess = true) {
3.540 + init();
3.541 + start(postprocess);
3.542 + }
3.543 +
3.544 + /// \brief Runs the algorithm to find a perfect fractional matching
3.545 + ///
3.546 + /// Just a shortcut for the next code:
3.547 + ///\code
3.548 + /// init();
3.549 + /// startPerfect();
3.550 + ///\endcode
3.551 + bool runPerfect(bool postprocess = true) {
3.552 + init();
3.553 + return startPerfect(postprocess);
3.554 + }
3.555 +
3.556 + ///@}
3.557 +
3.558 + /// \name Query Functions
3.559 + /// The result of the %Matching algorithm can be obtained using these
3.560 + /// functions.\n
3.561 + /// Before the use of these functions,
3.562 + /// either run() or start() must be called.
3.563 + ///@{
3.564 +
3.565 +
3.566 + /// \brief Return the number of covered nodes in the matching.
3.567 + ///
3.568 + /// This function returns the number of covered nodes in the matching.
3.569 + ///
3.570 + /// \pre Either run() or start() must be called before using this function.
3.571 + int matchingSize() const {
3.572 + int num = 0;
3.573 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.574 + if ((*_matching)[n] != INVALID) {
3.575 + ++num;
3.576 + }
3.577 + }
3.578 + return num;
3.579 + }
3.580 +
3.581 + /// \brief Returns a const reference to the matching map.
3.582 + ///
3.583 + /// Returns a const reference to the node map storing the found
3.584 + /// fractional matching. This method can be called after
3.585 + /// running the algorithm.
3.586 + ///
3.587 + /// \pre Either \ref run() or \ref init() must be called before
3.588 + /// using this function.
3.589 + const MatchingMap& matchingMap() const {
3.590 + return *_matching;
3.591 + }
3.592 +
3.593 + /// \brief Return \c true if the given edge is in the matching.
3.594 + ///
3.595 + /// This function returns \c true if the given edge is in the
3.596 + /// found matching. The result is scaled by \ref primalScale
3.597 + /// "primal scale".
3.598 + ///
3.599 + /// \pre Either run() or start() must be called before using this function.
3.600 + int matching(const Edge& edge) const {
3.601 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
3.602 + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
3.603 + }
3.604 +
3.605 + /// \brief Return the fractional matching arc (or edge) incident
3.606 + /// to the given node.
3.607 + ///
3.608 + /// This function returns one of the fractional matching arc (or
3.609 + /// edge) incident to the given node in the found matching or \c
3.610 + /// INVALID if the node is not covered by the matching or if the
3.611 + /// node is on an odd length cycle then it is the successor edge
3.612 + /// on the cycle.
3.613 + ///
3.614 + /// \pre Either run() or start() must be called before using this function.
3.615 + Arc matching(const Node& node) const {
3.616 + return (*_matching)[node];
3.617 + }
3.618 +
3.619 + /// \brief Returns true if the node is in the barrier
3.620 + ///
3.621 + /// The barrier is a subset of the nodes. If the nodes in the
3.622 + /// barrier have less adjacent nodes than the size of the barrier,
3.623 + /// then at least as much nodes cannot be covered as the
3.624 + /// difference of the two subsets.
3.625 + bool barrier(const Node& node) const {
3.626 + return (*_level)[node] >= _empty_level;
3.627 + }
3.628 +
3.629 + /// @}
3.630 +
3.631 + };
3.632 +
3.633 + /// \ingroup matching
3.634 + ///
3.635 + /// \brief Weighted fractional matching in general graphs
3.636 + ///
3.637 + /// This class provides an efficient implementation of fractional
3.638 + /// matching algorithm. The implementation uses priority queues and
3.639 + /// provides \f$O(nm\log n)\f$ time complexity.
3.640 + ///
3.641 + /// The maximum weighted fractional matching is a relaxation of the
3.642 + /// maximum weighted matching problem where the odd set constraints
3.643 + /// are omitted.
3.644 + /// It can be formulated with the following linear program.
3.645 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
3.646 + /// \f[x_e \ge 0\quad \forall e\in E\f]
3.647 + /// \f[\max \sum_{e\in E}x_ew_e\f]
3.648 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.649 + /// \f$X\f$. The result must be the union of a matching with one
3.650 + /// value edges and a set of odd length cycles with half value edges.
3.651 + ///
3.652 + /// The algorithm calculates an optimal fractional matching and a
3.653 + /// proof of the optimality. The solution of the dual problem can be
3.654 + /// used to check the result of the algorithm. The dual linear
3.655 + /// problem is the following.
3.656 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
3.657 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
3.658 + /// \f[\min \sum_{u \in V}y_u \f]
3.659 + ///
3.660 + /// The algorithm can be executed with the run() function.
3.661 + /// After it the matching (the primal solution) and the dual solution
3.662 + /// can be obtained using the query functions.
3.663 + ///
3.664 + /// The primal solution is multiplied by
3.665 + /// \ref MaxWeightedFractionalMatching::primalScale "2".
3.666 + /// If the value type is integer, then the dual
3.667 + /// solution is scaled by
3.668 + /// \ref MaxWeightedFractionalMatching::dualScale "4".
3.669 + ///
3.670 + /// \tparam GR The undirected graph type the algorithm runs on.
3.671 + /// \tparam WM The type edge weight map. The default type is
3.672 + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
3.673 +#ifdef DOXYGEN
3.674 + template <typename GR, typename WM>
3.675 +#else
3.676 + template <typename GR,
3.677 + typename WM = typename GR::template EdgeMap<int> >
3.678 +#endif
3.679 + class MaxWeightedFractionalMatching {
3.680 + public:
3.681 +
3.682 + /// The graph type of the algorithm
3.683 + typedef GR Graph;
3.684 + /// The type of the edge weight map
3.685 + typedef WM WeightMap;
3.686 + /// The value type of the edge weights
3.687 + typedef typename WeightMap::Value Value;
3.688 +
3.689 + /// The type of the matching map
3.690 + typedef typename Graph::template NodeMap<typename Graph::Arc>
3.691 + MatchingMap;
3.692 +
3.693 + /// \brief Scaling factor for primal solution
3.694 + ///
3.695 + /// Scaling factor for primal solution.
3.696 + static const int primalScale = 2;
3.697 +
3.698 + /// \brief Scaling factor for dual solution
3.699 + ///
3.700 + /// Scaling factor for dual solution. It is equal to 4 or 1
3.701 + /// according to the value type.
3.702 + static const int dualScale =
3.703 + std::numeric_limits<Value>::is_integer ? 4 : 1;
3.704 +
3.705 + private:
3.706 +
3.707 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
3.708 +
3.709 + typedef typename Graph::template NodeMap<Value> NodePotential;
3.710 +
3.711 + const Graph& _graph;
3.712 + const WeightMap& _weight;
3.713 +
3.714 + MatchingMap* _matching;
3.715 + NodePotential* _node_potential;
3.716 +
3.717 + int _node_num;
3.718 + bool _allow_loops;
3.719 +
3.720 + enum Status {
3.721 + EVEN = -1, MATCHED = 0, ODD = 1
3.722 + };
3.723 +
3.724 + typedef typename Graph::template NodeMap<Status> StatusMap;
3.725 + StatusMap* _status;
3.726 +
3.727 + typedef typename Graph::template NodeMap<Arc> PredMap;
3.728 + PredMap* _pred;
3.729 +
3.730 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
3.731 +
3.732 + IntNodeMap *_tree_set_index;
3.733 + TreeSet *_tree_set;
3.734 +
3.735 + IntNodeMap *_delta1_index;
3.736 + BinHeap<Value, IntNodeMap> *_delta1;
3.737 +
3.738 + IntNodeMap *_delta2_index;
3.739 + BinHeap<Value, IntNodeMap> *_delta2;
3.740 +
3.741 + IntEdgeMap *_delta3_index;
3.742 + BinHeap<Value, IntEdgeMap> *_delta3;
3.743 +
3.744 + Value _delta_sum;
3.745 +
3.746 + void createStructures() {
3.747 + _node_num = countNodes(_graph);
3.748 +
3.749 + if (!_matching) {
3.750 + _matching = new MatchingMap(_graph);
3.751 + }
3.752 + if (!_node_potential) {
3.753 + _node_potential = new NodePotential(_graph);
3.754 + }
3.755 + if (!_status) {
3.756 + _status = new StatusMap(_graph);
3.757 + }
3.758 + if (!_pred) {
3.759 + _pred = new PredMap(_graph);
3.760 + }
3.761 + if (!_tree_set) {
3.762 + _tree_set_index = new IntNodeMap(_graph);
3.763 + _tree_set = new TreeSet(*_tree_set_index);
3.764 + }
3.765 + if (!_delta1) {
3.766 + _delta1_index = new IntNodeMap(_graph);
3.767 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
3.768 + }
3.769 + if (!_delta2) {
3.770 + _delta2_index = new IntNodeMap(_graph);
3.771 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
3.772 + }
3.773 + if (!_delta3) {
3.774 + _delta3_index = new IntEdgeMap(_graph);
3.775 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
3.776 + }
3.777 + }
3.778 +
3.779 + void destroyStructures() {
3.780 + if (_matching) {
3.781 + delete _matching;
3.782 + }
3.783 + if (_node_potential) {
3.784 + delete _node_potential;
3.785 + }
3.786 + if (_status) {
3.787 + delete _status;
3.788 + }
3.789 + if (_pred) {
3.790 + delete _pred;
3.791 + }
3.792 + if (_tree_set) {
3.793 + delete _tree_set_index;
3.794 + delete _tree_set;
3.795 + }
3.796 + if (_delta1) {
3.797 + delete _delta1_index;
3.798 + delete _delta1;
3.799 + }
3.800 + if (_delta2) {
3.801 + delete _delta2_index;
3.802 + delete _delta2;
3.803 + }
3.804 + if (_delta3) {
3.805 + delete _delta3_index;
3.806 + delete _delta3;
3.807 + }
3.808 + }
3.809 +
3.810 + void matchedToEven(Node node, int tree) {
3.811 + _tree_set->insert(node, tree);
3.812 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.813 + _delta1->push(node, (*_node_potential)[node]);
3.814 +
3.815 + if (_delta2->state(node) == _delta2->IN_HEAP) {
3.816 + _delta2->erase(node);
3.817 + }
3.818 +
3.819 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.820 + Node v = _graph.source(a);
3.821 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.822 + dualScale * _weight[a];
3.823 + if (node == v) {
3.824 + if (_allow_loops && _graph.direction(a)) {
3.825 + _delta3->push(a, rw / 2);
3.826 + }
3.827 + } else if ((*_status)[v] == EVEN) {
3.828 + _delta3->push(a, rw / 2);
3.829 + } else if ((*_status)[v] == MATCHED) {
3.830 + if (_delta2->state(v) != _delta2->IN_HEAP) {
3.831 + _pred->set(v, a);
3.832 + _delta2->push(v, rw);
3.833 + } else if ((*_delta2)[v] > rw) {
3.834 + _pred->set(v, a);
3.835 + _delta2->decrease(v, rw);
3.836 + }
3.837 + }
3.838 + }
3.839 + }
3.840 +
3.841 + void matchedToOdd(Node node, int tree) {
3.842 + _tree_set->insert(node, tree);
3.843 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.844 +
3.845 + if (_delta2->state(node) == _delta2->IN_HEAP) {
3.846 + _delta2->erase(node);
3.847 + }
3.848 + }
3.849 +
3.850 + void evenToMatched(Node node, int tree) {
3.851 + _delta1->erase(node);
3.852 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.853 + Arc min = INVALID;
3.854 + Value minrw = std::numeric_limits<Value>::max();
3.855 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.856 + Node v = _graph.source(a);
3.857 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.858 + dualScale * _weight[a];
3.859 +
3.860 + if (node == v) {
3.861 + if (_allow_loops && _graph.direction(a)) {
3.862 + _delta3->erase(a);
3.863 + }
3.864 + } else if ((*_status)[v] == EVEN) {
3.865 + _delta3->erase(a);
3.866 + if (minrw > rw) {
3.867 + min = _graph.oppositeArc(a);
3.868 + minrw = rw;
3.869 + }
3.870 + } else if ((*_status)[v] == MATCHED) {
3.871 + if ((*_pred)[v] == a) {
3.872 + Arc mina = INVALID;
3.873 + Value minrwa = std::numeric_limits<Value>::max();
3.874 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
3.875 + Node va = _graph.target(aa);
3.876 + if ((*_status)[va] != EVEN ||
3.877 + _tree_set->find(va) == tree) continue;
3.878 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
3.879 + dualScale * _weight[aa];
3.880 + if (minrwa > rwa) {
3.881 + minrwa = rwa;
3.882 + mina = aa;
3.883 + }
3.884 + }
3.885 + if (mina != INVALID) {
3.886 + _pred->set(v, mina);
3.887 + _delta2->increase(v, minrwa);
3.888 + } else {
3.889 + _pred->set(v, INVALID);
3.890 + _delta2->erase(v);
3.891 + }
3.892 + }
3.893 + }
3.894 + }
3.895 + if (min != INVALID) {
3.896 + _pred->set(node, min);
3.897 + _delta2->push(node, minrw);
3.898 + } else {
3.899 + _pred->set(node, INVALID);
3.900 + }
3.901 + }
3.902 +
3.903 + void oddToMatched(Node node) {
3.904 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.905 + Arc min = INVALID;
3.906 + Value minrw = std::numeric_limits<Value>::max();
3.907 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.908 + Node v = _graph.source(a);
3.909 + if ((*_status)[v] != EVEN) continue;
3.910 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.911 + dualScale * _weight[a];
3.912 +
3.913 + if (minrw > rw) {
3.914 + min = _graph.oppositeArc(a);
3.915 + minrw = rw;
3.916 + }
3.917 + }
3.918 + if (min != INVALID) {
3.919 + _pred->set(node, min);
3.920 + _delta2->push(node, minrw);
3.921 + } else {
3.922 + _pred->set(node, INVALID);
3.923 + }
3.924 + }
3.925 +
3.926 + void alternatePath(Node even, int tree) {
3.927 + Node odd;
3.928 +
3.929 + _status->set(even, MATCHED);
3.930 + evenToMatched(even, tree);
3.931 +
3.932 + Arc prev = (*_matching)[even];
3.933 + while (prev != INVALID) {
3.934 + odd = _graph.target(prev);
3.935 + even = _graph.target((*_pred)[odd]);
3.936 + _matching->set(odd, (*_pred)[odd]);
3.937 + _status->set(odd, MATCHED);
3.938 + oddToMatched(odd);
3.939 +
3.940 + prev = (*_matching)[even];
3.941 + _status->set(even, MATCHED);
3.942 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
3.943 + evenToMatched(even, tree);
3.944 + }
3.945 + }
3.946 +
3.947 + void destroyTree(int tree) {
3.948 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
3.949 + if ((*_status)[n] == EVEN) {
3.950 + _status->set(n, MATCHED);
3.951 + evenToMatched(n, tree);
3.952 + } else if ((*_status)[n] == ODD) {
3.953 + _status->set(n, MATCHED);
3.954 + oddToMatched(n);
3.955 + }
3.956 + }
3.957 + _tree_set->eraseClass(tree);
3.958 + }
3.959 +
3.960 +
3.961 + void unmatchNode(const Node& node) {
3.962 + int tree = _tree_set->find(node);
3.963 +
3.964 + alternatePath(node, tree);
3.965 + destroyTree(tree);
3.966 +
3.967 + _matching->set(node, INVALID);
3.968 + }
3.969 +
3.970 +
3.971 + void augmentOnEdge(const Edge& edge) {
3.972 + Node left = _graph.u(edge);
3.973 + int left_tree = _tree_set->find(left);
3.974 +
3.975 + alternatePath(left, left_tree);
3.976 + destroyTree(left_tree);
3.977 + _matching->set(left, _graph.direct(edge, true));
3.978 +
3.979 + Node right = _graph.v(edge);
3.980 + int right_tree = _tree_set->find(right);
3.981 +
3.982 + alternatePath(right, right_tree);
3.983 + destroyTree(right_tree);
3.984 + _matching->set(right, _graph.direct(edge, false));
3.985 + }
3.986 +
3.987 + void augmentOnArc(const Arc& arc) {
3.988 + Node left = _graph.source(arc);
3.989 + _status->set(left, MATCHED);
3.990 + _matching->set(left, arc);
3.991 + _pred->set(left, arc);
3.992 +
3.993 + Node right = _graph.target(arc);
3.994 + int right_tree = _tree_set->find(right);
3.995 +
3.996 + alternatePath(right, right_tree);
3.997 + destroyTree(right_tree);
3.998 + _matching->set(right, _graph.oppositeArc(arc));
3.999 + }
3.1000 +
3.1001 + void extendOnArc(const Arc& arc) {
3.1002 + Node base = _graph.target(arc);
3.1003 + int tree = _tree_set->find(base);
3.1004 +
3.1005 + Node odd = _graph.source(arc);
3.1006 + _tree_set->insert(odd, tree);
3.1007 + _status->set(odd, ODD);
3.1008 + matchedToOdd(odd, tree);
3.1009 + _pred->set(odd, arc);
3.1010 +
3.1011 + Node even = _graph.target((*_matching)[odd]);
3.1012 + _tree_set->insert(even, tree);
3.1013 + _status->set(even, EVEN);
3.1014 + matchedToEven(even, tree);
3.1015 + }
3.1016 +
3.1017 + void cycleOnEdge(const Edge& edge, int tree) {
3.1018 + Node nca = INVALID;
3.1019 + std::vector<Node> left_path, right_path;
3.1020 +
3.1021 + {
3.1022 + std::set<Node> left_set, right_set;
3.1023 + Node left = _graph.u(edge);
3.1024 + left_path.push_back(left);
3.1025 + left_set.insert(left);
3.1026 +
3.1027 + Node right = _graph.v(edge);
3.1028 + right_path.push_back(right);
3.1029 + right_set.insert(right);
3.1030 +
3.1031 + while (true) {
3.1032 +
3.1033 + if (left_set.find(right) != left_set.end()) {
3.1034 + nca = right;
3.1035 + break;
3.1036 + }
3.1037 +
3.1038 + if ((*_matching)[left] == INVALID) break;
3.1039 +
3.1040 + left = _graph.target((*_matching)[left]);
3.1041 + left_path.push_back(left);
3.1042 + left = _graph.target((*_pred)[left]);
3.1043 + left_path.push_back(left);
3.1044 +
3.1045 + left_set.insert(left);
3.1046 +
3.1047 + if (right_set.find(left) != right_set.end()) {
3.1048 + nca = left;
3.1049 + break;
3.1050 + }
3.1051 +
3.1052 + if ((*_matching)[right] == INVALID) break;
3.1053 +
3.1054 + right = _graph.target((*_matching)[right]);
3.1055 + right_path.push_back(right);
3.1056 + right = _graph.target((*_pred)[right]);
3.1057 + right_path.push_back(right);
3.1058 +
3.1059 + right_set.insert(right);
3.1060 +
3.1061 + }
3.1062 +
3.1063 + if (nca == INVALID) {
3.1064 + if ((*_matching)[left] == INVALID) {
3.1065 + nca = right;
3.1066 + while (left_set.find(nca) == left_set.end()) {
3.1067 + nca = _graph.target((*_matching)[nca]);
3.1068 + right_path.push_back(nca);
3.1069 + nca = _graph.target((*_pred)[nca]);
3.1070 + right_path.push_back(nca);
3.1071 + }
3.1072 + } else {
3.1073 + nca = left;
3.1074 + while (right_set.find(nca) == right_set.end()) {
3.1075 + nca = _graph.target((*_matching)[nca]);
3.1076 + left_path.push_back(nca);
3.1077 + nca = _graph.target((*_pred)[nca]);
3.1078 + left_path.push_back(nca);
3.1079 + }
3.1080 + }
3.1081 + }
3.1082 + }
3.1083 +
3.1084 + alternatePath(nca, tree);
3.1085 + Arc prev;
3.1086 +
3.1087 + prev = _graph.direct(edge, true);
3.1088 + for (int i = 0; left_path[i] != nca; i += 2) {
3.1089 + _matching->set(left_path[i], prev);
3.1090 + _status->set(left_path[i], MATCHED);
3.1091 + evenToMatched(left_path[i], tree);
3.1092 +
3.1093 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
3.1094 + _status->set(left_path[i + 1], MATCHED);
3.1095 + oddToMatched(left_path[i + 1]);
3.1096 + }
3.1097 + _matching->set(nca, prev);
3.1098 +
3.1099 + for (int i = 0; right_path[i] != nca; i += 2) {
3.1100 + _status->set(right_path[i], MATCHED);
3.1101 + evenToMatched(right_path[i], tree);
3.1102 +
3.1103 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
3.1104 + _status->set(right_path[i + 1], MATCHED);
3.1105 + oddToMatched(right_path[i + 1]);
3.1106 + }
3.1107 +
3.1108 + destroyTree(tree);
3.1109 + }
3.1110 +
3.1111 + void extractCycle(const Arc &arc) {
3.1112 + Node left = _graph.source(arc);
3.1113 + Node odd = _graph.target((*_matching)[left]);
3.1114 + Arc prev;
3.1115 + while (odd != left) {
3.1116 + Node even = _graph.target((*_matching)[odd]);
3.1117 + prev = (*_matching)[odd];
3.1118 + odd = _graph.target((*_matching)[even]);
3.1119 + _matching->set(even, _graph.oppositeArc(prev));
3.1120 + }
3.1121 + _matching->set(left, arc);
3.1122 +
3.1123 + Node right = _graph.target(arc);
3.1124 + int right_tree = _tree_set->find(right);
3.1125 + alternatePath(right, right_tree);
3.1126 + destroyTree(right_tree);
3.1127 + _matching->set(right, _graph.oppositeArc(arc));
3.1128 + }
3.1129 +
3.1130 + public:
3.1131 +
3.1132 + /// \brief Constructor
3.1133 + ///
3.1134 + /// Constructor.
3.1135 + MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
3.1136 + bool allow_loops = true)
3.1137 + : _graph(graph), _weight(weight), _matching(0),
3.1138 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
3.1139 + _status(0), _pred(0),
3.1140 + _tree_set_index(0), _tree_set(0),
3.1141 +
3.1142 + _delta1_index(0), _delta1(0),
3.1143 + _delta2_index(0), _delta2(0),
3.1144 + _delta3_index(0), _delta3(0),
3.1145 +
3.1146 + _delta_sum() {}
3.1147 +
3.1148 + ~MaxWeightedFractionalMatching() {
3.1149 + destroyStructures();
3.1150 + }
3.1151 +
3.1152 + /// \name Execution Control
3.1153 + /// The simplest way to execute the algorithm is to use the
3.1154 + /// \ref run() member function.
3.1155 +
3.1156 + ///@{
3.1157 +
3.1158 + /// \brief Initialize the algorithm
3.1159 + ///
3.1160 + /// This function initializes the algorithm.
3.1161 + void init() {
3.1162 + createStructures();
3.1163 +
3.1164 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1165 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
3.1166 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
3.1167 + }
3.1168 + for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1169 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
3.1170 + }
3.1171 +
3.1172 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1173 + Value max = 0;
3.1174 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3.1175 + if (_graph.target(e) == n && !_allow_loops) continue;
3.1176 + if ((dualScale * _weight[e]) / 2 > max) {
3.1177 + max = (dualScale * _weight[e]) / 2;
3.1178 + }
3.1179 + }
3.1180 + _node_potential->set(n, max);
3.1181 + _delta1->push(n, max);
3.1182 +
3.1183 + _tree_set->insert(n);
3.1184 +
3.1185 + _matching->set(n, INVALID);
3.1186 + _status->set(n, EVEN);
3.1187 + }
3.1188 +
3.1189 + for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1190 + Node left = _graph.u(e);
3.1191 + Node right = _graph.v(e);
3.1192 + if (left == right && !_allow_loops) continue;
3.1193 + _delta3->push(e, ((*_node_potential)[left] +
3.1194 + (*_node_potential)[right] -
3.1195 + dualScale * _weight[e]) / 2);
3.1196 + }
3.1197 + }
3.1198 +
3.1199 + /// \brief Start the algorithm
3.1200 + ///
3.1201 + /// This function starts the algorithm.
3.1202 + ///
3.1203 + /// \pre \ref init() must be called before using this function.
3.1204 + void start() {
3.1205 + enum OpType {
3.1206 + D1, D2, D3
3.1207 + };
3.1208 +
3.1209 + int unmatched = _node_num;
3.1210 + while (unmatched > 0) {
3.1211 + Value d1 = !_delta1->empty() ?
3.1212 + _delta1->prio() : std::numeric_limits<Value>::max();
3.1213 +
3.1214 + Value d2 = !_delta2->empty() ?
3.1215 + _delta2->prio() : std::numeric_limits<Value>::max();
3.1216 +
3.1217 + Value d3 = !_delta3->empty() ?
3.1218 + _delta3->prio() : std::numeric_limits<Value>::max();
3.1219 +
3.1220 + _delta_sum = d3; OpType ot = D3;
3.1221 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
3.1222 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3.1223 +
3.1224 + switch (ot) {
3.1225 + case D1:
3.1226 + {
3.1227 + Node n = _delta1->top();
3.1228 + unmatchNode(n);
3.1229 + --unmatched;
3.1230 + }
3.1231 + break;
3.1232 + case D2:
3.1233 + {
3.1234 + Node n = _delta2->top();
3.1235 + Arc a = (*_pred)[n];
3.1236 + if ((*_matching)[n] == INVALID) {
3.1237 + augmentOnArc(a);
3.1238 + --unmatched;
3.1239 + } else {
3.1240 + Node v = _graph.target((*_matching)[n]);
3.1241 + if ((*_matching)[n] !=
3.1242 + _graph.oppositeArc((*_matching)[v])) {
3.1243 + extractCycle(a);
3.1244 + --unmatched;
3.1245 + } else {
3.1246 + extendOnArc(a);
3.1247 + }
3.1248 + }
3.1249 + } break;
3.1250 + case D3:
3.1251 + {
3.1252 + Edge e = _delta3->top();
3.1253 +
3.1254 + Node left = _graph.u(e);
3.1255 + Node right = _graph.v(e);
3.1256 +
3.1257 + int left_tree = _tree_set->find(left);
3.1258 + int right_tree = _tree_set->find(right);
3.1259 +
3.1260 + if (left_tree == right_tree) {
3.1261 + cycleOnEdge(e, left_tree);
3.1262 + --unmatched;
3.1263 + } else {
3.1264 + augmentOnEdge(e);
3.1265 + unmatched -= 2;
3.1266 + }
3.1267 + } break;
3.1268 + }
3.1269 + }
3.1270 + }
3.1271 +
3.1272 + /// \brief Run the algorithm.
3.1273 + ///
3.1274 + /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
3.1275 + ///
3.1276 + /// \note mwfm.run() is just a shortcut of the following code.
3.1277 + /// \code
3.1278 + /// mwfm.init();
3.1279 + /// mwfm.start();
3.1280 + /// \endcode
3.1281 + void run() {
3.1282 + init();
3.1283 + start();
3.1284 + }
3.1285 +
3.1286 + /// @}
3.1287 +
3.1288 + /// \name Primal Solution
3.1289 + /// Functions to get the primal solution, i.e. the maximum weighted
3.1290 + /// matching.\n
3.1291 + /// Either \ref run() or \ref start() function should be called before
3.1292 + /// using them.
3.1293 +
3.1294 + /// @{
3.1295 +
3.1296 + /// \brief Return the weight of the matching.
3.1297 + ///
3.1298 + /// This function returns the weight of the found matching. This
3.1299 + /// value is scaled by \ref primalScale "primal scale".
3.1300 + ///
3.1301 + /// \pre Either run() or start() must be called before using this function.
3.1302 + Value matchingWeight() const {
3.1303 + Value sum = 0;
3.1304 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1305 + if ((*_matching)[n] != INVALID) {
3.1306 + sum += _weight[(*_matching)[n]];
3.1307 + }
3.1308 + }
3.1309 + return sum * primalScale / 2;
3.1310 + }
3.1311 +
3.1312 + /// \brief Return the number of covered nodes in the matching.
3.1313 + ///
3.1314 + /// This function returns the number of covered nodes in the matching.
3.1315 + ///
3.1316 + /// \pre Either run() or start() must be called before using this function.
3.1317 + int matchingSize() const {
3.1318 + int num = 0;
3.1319 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1320 + if ((*_matching)[n] != INVALID) {
3.1321 + ++num;
3.1322 + }
3.1323 + }
3.1324 + return num;
3.1325 + }
3.1326 +
3.1327 + /// \brief Return \c true if the given edge is in the matching.
3.1328 + ///
3.1329 + /// This function returns \c true if the given edge is in the
3.1330 + /// found matching. The result is scaled by \ref primalScale
3.1331 + /// "primal scale".
3.1332 + ///
3.1333 + /// \pre Either run() or start() must be called before using this function.
3.1334 + int matching(const Edge& edge) const {
3.1335 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
3.1336 + + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
3.1337 + }
3.1338 +
3.1339 + /// \brief Return the fractional matching arc (or edge) incident
3.1340 + /// to the given node.
3.1341 + ///
3.1342 + /// This function returns one of the fractional matching arc (or
3.1343 + /// edge) incident to the given node in the found matching or \c
3.1344 + /// INVALID if the node is not covered by the matching or if the
3.1345 + /// node is on an odd length cycle then it is the successor edge
3.1346 + /// on the cycle.
3.1347 + ///
3.1348 + /// \pre Either run() or start() must be called before using this function.
3.1349 + Arc matching(const Node& node) const {
3.1350 + return (*_matching)[node];
3.1351 + }
3.1352 +
3.1353 + /// \brief Return a const reference to the matching map.
3.1354 + ///
3.1355 + /// This function returns a const reference to a node map that stores
3.1356 + /// the matching arc (or edge) incident to each node.
3.1357 + const MatchingMap& matchingMap() const {
3.1358 + return *_matching;
3.1359 + }
3.1360 +
3.1361 + /// @}
3.1362 +
3.1363 + /// \name Dual Solution
3.1364 + /// Functions to get the dual solution.\n
3.1365 + /// Either \ref run() or \ref start() function should be called before
3.1366 + /// using them.
3.1367 +
3.1368 + /// @{
3.1369 +
3.1370 + /// \brief Return the value of the dual solution.
3.1371 + ///
3.1372 + /// This function returns the value of the dual solution.
3.1373 + /// It should be equal to the primal value scaled by \ref dualScale
3.1374 + /// "dual scale".
3.1375 + ///
3.1376 + /// \pre Either run() or start() must be called before using this function.
3.1377 + Value dualValue() const {
3.1378 + Value sum = 0;
3.1379 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1380 + sum += nodeValue(n);
3.1381 + }
3.1382 + return sum;
3.1383 + }
3.1384 +
3.1385 + /// \brief Return the dual value (potential) of the given node.
3.1386 + ///
3.1387 + /// This function returns the dual value (potential) of the given node.
3.1388 + ///
3.1389 + /// \pre Either run() or start() must be called before using this function.
3.1390 + Value nodeValue(const Node& n) const {
3.1391 + return (*_node_potential)[n];
3.1392 + }
3.1393 +
3.1394 + /// @}
3.1395 +
3.1396 + };
3.1397 +
3.1398 + /// \ingroup matching
3.1399 + ///
3.1400 + /// \brief Weighted fractional perfect matching in general graphs
3.1401 + ///
3.1402 + /// This class provides an efficient implementation of fractional
3.1403 + /// matching algorithm. The implementation uses priority queues and
3.1404 + /// provides \f$O(nm\log n)\f$ time complexity.
3.1405 + ///
3.1406 + /// The maximum weighted fractional perfect matching is a relaxation
3.1407 + /// of the maximum weighted perfect matching problem where the odd
3.1408 + /// set constraints are omitted.
3.1409 + /// It can be formulated with the following linear program.
3.1410 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
3.1411 + /// \f[x_e \ge 0\quad \forall e\in E\f]
3.1412 + /// \f[\max \sum_{e\in E}x_ew_e\f]
3.1413 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
3.1414 + /// \f$X\f$. The result must be the union of a matching with one
3.1415 + /// value edges and a set of odd length cycles with half value edges.
3.1416 + ///
3.1417 + /// The algorithm calculates an optimal fractional matching and a
3.1418 + /// proof of the optimality. The solution of the dual problem can be
3.1419 + /// used to check the result of the algorithm. The dual linear
3.1420 + /// problem is the following.
3.1421 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
3.1422 + /// \f[\min \sum_{u \in V}y_u \f]
3.1423 + ///
3.1424 + /// The algorithm can be executed with the run() function.
3.1425 + /// After it the matching (the primal solution) and the dual solution
3.1426 + /// can be obtained using the query functions.
3.1427 + ///
3.1428 + /// The primal solution is multiplied by
3.1429 + /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
3.1430 + /// If the value type is integer, then the dual
3.1431 + /// solution is scaled by
3.1432 + /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
3.1433 + ///
3.1434 + /// \tparam GR The undirected graph type the algorithm runs on.
3.1435 + /// \tparam WM The type edge weight map. The default type is
3.1436 + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
3.1437 +#ifdef DOXYGEN
3.1438 + template <typename GR, typename WM>
3.1439 +#else
3.1440 + template <typename GR,
3.1441 + typename WM = typename GR::template EdgeMap<int> >
3.1442 +#endif
3.1443 + class MaxWeightedPerfectFractionalMatching {
3.1444 + public:
3.1445 +
3.1446 + /// The graph type of the algorithm
3.1447 + typedef GR Graph;
3.1448 + /// The type of the edge weight map
3.1449 + typedef WM WeightMap;
3.1450 + /// The value type of the edge weights
3.1451 + typedef typename WeightMap::Value Value;
3.1452 +
3.1453 + /// The type of the matching map
3.1454 + typedef typename Graph::template NodeMap<typename Graph::Arc>
3.1455 + MatchingMap;
3.1456 +
3.1457 + /// \brief Scaling factor for primal solution
3.1458 + ///
3.1459 + /// Scaling factor for primal solution.
3.1460 + static const int primalScale = 2;
3.1461 +
3.1462 + /// \brief Scaling factor for dual solution
3.1463 + ///
3.1464 + /// Scaling factor for dual solution. It is equal to 4 or 1
3.1465 + /// according to the value type.
3.1466 + static const int dualScale =
3.1467 + std::numeric_limits<Value>::is_integer ? 4 : 1;
3.1468 +
3.1469 + private:
3.1470 +
3.1471 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
3.1472 +
3.1473 + typedef typename Graph::template NodeMap<Value> NodePotential;
3.1474 +
3.1475 + const Graph& _graph;
3.1476 + const WeightMap& _weight;
3.1477 +
3.1478 + MatchingMap* _matching;
3.1479 + NodePotential* _node_potential;
3.1480 +
3.1481 + int _node_num;
3.1482 + bool _allow_loops;
3.1483 +
3.1484 + enum Status {
3.1485 + EVEN = -1, MATCHED = 0, ODD = 1
3.1486 + };
3.1487 +
3.1488 + typedef typename Graph::template NodeMap<Status> StatusMap;
3.1489 + StatusMap* _status;
3.1490 +
3.1491 + typedef typename Graph::template NodeMap<Arc> PredMap;
3.1492 + PredMap* _pred;
3.1493 +
3.1494 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
3.1495 +
3.1496 + IntNodeMap *_tree_set_index;
3.1497 + TreeSet *_tree_set;
3.1498 +
3.1499 + IntNodeMap *_delta2_index;
3.1500 + BinHeap<Value, IntNodeMap> *_delta2;
3.1501 +
3.1502 + IntEdgeMap *_delta3_index;
3.1503 + BinHeap<Value, IntEdgeMap> *_delta3;
3.1504 +
3.1505 + Value _delta_sum;
3.1506 +
3.1507 + void createStructures() {
3.1508 + _node_num = countNodes(_graph);
3.1509 +
3.1510 + if (!_matching) {
3.1511 + _matching = new MatchingMap(_graph);
3.1512 + }
3.1513 + if (!_node_potential) {
3.1514 + _node_potential = new NodePotential(_graph);
3.1515 + }
3.1516 + if (!_status) {
3.1517 + _status = new StatusMap(_graph);
3.1518 + }
3.1519 + if (!_pred) {
3.1520 + _pred = new PredMap(_graph);
3.1521 + }
3.1522 + if (!_tree_set) {
3.1523 + _tree_set_index = new IntNodeMap(_graph);
3.1524 + _tree_set = new TreeSet(*_tree_set_index);
3.1525 + }
3.1526 + if (!_delta2) {
3.1527 + _delta2_index = new IntNodeMap(_graph);
3.1528 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
3.1529 + }
3.1530 + if (!_delta3) {
3.1531 + _delta3_index = new IntEdgeMap(_graph);
3.1532 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
3.1533 + }
3.1534 + }
3.1535 +
3.1536 + void destroyStructures() {
3.1537 + if (_matching) {
3.1538 + delete _matching;
3.1539 + }
3.1540 + if (_node_potential) {
3.1541 + delete _node_potential;
3.1542 + }
3.1543 + if (_status) {
3.1544 + delete _status;
3.1545 + }
3.1546 + if (_pred) {
3.1547 + delete _pred;
3.1548 + }
3.1549 + if (_tree_set) {
3.1550 + delete _tree_set_index;
3.1551 + delete _tree_set;
3.1552 + }
3.1553 + if (_delta2) {
3.1554 + delete _delta2_index;
3.1555 + delete _delta2;
3.1556 + }
3.1557 + if (_delta3) {
3.1558 + delete _delta3_index;
3.1559 + delete _delta3;
3.1560 + }
3.1561 + }
3.1562 +
3.1563 + void matchedToEven(Node node, int tree) {
3.1564 + _tree_set->insert(node, tree);
3.1565 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.1566 +
3.1567 + if (_delta2->state(node) == _delta2->IN_HEAP) {
3.1568 + _delta2->erase(node);
3.1569 + }
3.1570 +
3.1571 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1572 + Node v = _graph.source(a);
3.1573 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1574 + dualScale * _weight[a];
3.1575 + if (node == v) {
3.1576 + if (_allow_loops && _graph.direction(a)) {
3.1577 + _delta3->push(a, rw / 2);
3.1578 + }
3.1579 + } else if ((*_status)[v] == EVEN) {
3.1580 + _delta3->push(a, rw / 2);
3.1581 + } else if ((*_status)[v] == MATCHED) {
3.1582 + if (_delta2->state(v) != _delta2->IN_HEAP) {
3.1583 + _pred->set(v, a);
3.1584 + _delta2->push(v, rw);
3.1585 + } else if ((*_delta2)[v] > rw) {
3.1586 + _pred->set(v, a);
3.1587 + _delta2->decrease(v, rw);
3.1588 + }
3.1589 + }
3.1590 + }
3.1591 + }
3.1592 +
3.1593 + void matchedToOdd(Node node, int tree) {
3.1594 + _tree_set->insert(node, tree);
3.1595 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.1596 +
3.1597 + if (_delta2->state(node) == _delta2->IN_HEAP) {
3.1598 + _delta2->erase(node);
3.1599 + }
3.1600 + }
3.1601 +
3.1602 + void evenToMatched(Node node, int tree) {
3.1603 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
3.1604 + Arc min = INVALID;
3.1605 + Value minrw = std::numeric_limits<Value>::max();
3.1606 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1607 + Node v = _graph.source(a);
3.1608 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1609 + dualScale * _weight[a];
3.1610 +
3.1611 + if (node == v) {
3.1612 + if (_allow_loops && _graph.direction(a)) {
3.1613 + _delta3->erase(a);
3.1614 + }
3.1615 + } else if ((*_status)[v] == EVEN) {
3.1616 + _delta3->erase(a);
3.1617 + if (minrw > rw) {
3.1618 + min = _graph.oppositeArc(a);
3.1619 + minrw = rw;
3.1620 + }
3.1621 + } else if ((*_status)[v] == MATCHED) {
3.1622 + if ((*_pred)[v] == a) {
3.1623 + Arc mina = INVALID;
3.1624 + Value minrwa = std::numeric_limits<Value>::max();
3.1625 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
3.1626 + Node va = _graph.target(aa);
3.1627 + if ((*_status)[va] != EVEN ||
3.1628 + _tree_set->find(va) == tree) continue;
3.1629 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
3.1630 + dualScale * _weight[aa];
3.1631 + if (minrwa > rwa) {
3.1632 + minrwa = rwa;
3.1633 + mina = aa;
3.1634 + }
3.1635 + }
3.1636 + if (mina != INVALID) {
3.1637 + _pred->set(v, mina);
3.1638 + _delta2->increase(v, minrwa);
3.1639 + } else {
3.1640 + _pred->set(v, INVALID);
3.1641 + _delta2->erase(v);
3.1642 + }
3.1643 + }
3.1644 + }
3.1645 + }
3.1646 + if (min != INVALID) {
3.1647 + _pred->set(node, min);
3.1648 + _delta2->push(node, minrw);
3.1649 + } else {
3.1650 + _pred->set(node, INVALID);
3.1651 + }
3.1652 + }
3.1653 +
3.1654 + void oddToMatched(Node node) {
3.1655 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
3.1656 + Arc min = INVALID;
3.1657 + Value minrw = std::numeric_limits<Value>::max();
3.1658 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
3.1659 + Node v = _graph.source(a);
3.1660 + if ((*_status)[v] != EVEN) continue;
3.1661 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
3.1662 + dualScale * _weight[a];
3.1663 +
3.1664 + if (minrw > rw) {
3.1665 + min = _graph.oppositeArc(a);
3.1666 + minrw = rw;
3.1667 + }
3.1668 + }
3.1669 + if (min != INVALID) {
3.1670 + _pred->set(node, min);
3.1671 + _delta2->push(node, minrw);
3.1672 + } else {
3.1673 + _pred->set(node, INVALID);
3.1674 + }
3.1675 + }
3.1676 +
3.1677 + void alternatePath(Node even, int tree) {
3.1678 + Node odd;
3.1679 +
3.1680 + _status->set(even, MATCHED);
3.1681 + evenToMatched(even, tree);
3.1682 +
3.1683 + Arc prev = (*_matching)[even];
3.1684 + while (prev != INVALID) {
3.1685 + odd = _graph.target(prev);
3.1686 + even = _graph.target((*_pred)[odd]);
3.1687 + _matching->set(odd, (*_pred)[odd]);
3.1688 + _status->set(odd, MATCHED);
3.1689 + oddToMatched(odd);
3.1690 +
3.1691 + prev = (*_matching)[even];
3.1692 + _status->set(even, MATCHED);
3.1693 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
3.1694 + evenToMatched(even, tree);
3.1695 + }
3.1696 + }
3.1697 +
3.1698 + void destroyTree(int tree) {
3.1699 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
3.1700 + if ((*_status)[n] == EVEN) {
3.1701 + _status->set(n, MATCHED);
3.1702 + evenToMatched(n, tree);
3.1703 + } else if ((*_status)[n] == ODD) {
3.1704 + _status->set(n, MATCHED);
3.1705 + oddToMatched(n);
3.1706 + }
3.1707 + }
3.1708 + _tree_set->eraseClass(tree);
3.1709 + }
3.1710 +
3.1711 + void augmentOnEdge(const Edge& edge) {
3.1712 + Node left = _graph.u(edge);
3.1713 + int left_tree = _tree_set->find(left);
3.1714 +
3.1715 + alternatePath(left, left_tree);
3.1716 + destroyTree(left_tree);
3.1717 + _matching->set(left, _graph.direct(edge, true));
3.1718 +
3.1719 + Node right = _graph.v(edge);
3.1720 + int right_tree = _tree_set->find(right);
3.1721 +
3.1722 + alternatePath(right, right_tree);
3.1723 + destroyTree(right_tree);
3.1724 + _matching->set(right, _graph.direct(edge, false));
3.1725 + }
3.1726 +
3.1727 + void augmentOnArc(const Arc& arc) {
3.1728 + Node left = _graph.source(arc);
3.1729 + _status->set(left, MATCHED);
3.1730 + _matching->set(left, arc);
3.1731 + _pred->set(left, arc);
3.1732 +
3.1733 + Node right = _graph.target(arc);
3.1734 + int right_tree = _tree_set->find(right);
3.1735 +
3.1736 + alternatePath(right, right_tree);
3.1737 + destroyTree(right_tree);
3.1738 + _matching->set(right, _graph.oppositeArc(arc));
3.1739 + }
3.1740 +
3.1741 + void extendOnArc(const Arc& arc) {
3.1742 + Node base = _graph.target(arc);
3.1743 + int tree = _tree_set->find(base);
3.1744 +
3.1745 + Node odd = _graph.source(arc);
3.1746 + _tree_set->insert(odd, tree);
3.1747 + _status->set(odd, ODD);
3.1748 + matchedToOdd(odd, tree);
3.1749 + _pred->set(odd, arc);
3.1750 +
3.1751 + Node even = _graph.target((*_matching)[odd]);
3.1752 + _tree_set->insert(even, tree);
3.1753 + _status->set(even, EVEN);
3.1754 + matchedToEven(even, tree);
3.1755 + }
3.1756 +
3.1757 + void cycleOnEdge(const Edge& edge, int tree) {
3.1758 + Node nca = INVALID;
3.1759 + std::vector<Node> left_path, right_path;
3.1760 +
3.1761 + {
3.1762 + std::set<Node> left_set, right_set;
3.1763 + Node left = _graph.u(edge);
3.1764 + left_path.push_back(left);
3.1765 + left_set.insert(left);
3.1766 +
3.1767 + Node right = _graph.v(edge);
3.1768 + right_path.push_back(right);
3.1769 + right_set.insert(right);
3.1770 +
3.1771 + while (true) {
3.1772 +
3.1773 + if (left_set.find(right) != left_set.end()) {
3.1774 + nca = right;
3.1775 + break;
3.1776 + }
3.1777 +
3.1778 + if ((*_matching)[left] == INVALID) break;
3.1779 +
3.1780 + left = _graph.target((*_matching)[left]);
3.1781 + left_path.push_back(left);
3.1782 + left = _graph.target((*_pred)[left]);
3.1783 + left_path.push_back(left);
3.1784 +
3.1785 + left_set.insert(left);
3.1786 +
3.1787 + if (right_set.find(left) != right_set.end()) {
3.1788 + nca = left;
3.1789 + break;
3.1790 + }
3.1791 +
3.1792 + if ((*_matching)[right] == INVALID) break;
3.1793 +
3.1794 + right = _graph.target((*_matching)[right]);
3.1795 + right_path.push_back(right);
3.1796 + right = _graph.target((*_pred)[right]);
3.1797 + right_path.push_back(right);
3.1798 +
3.1799 + right_set.insert(right);
3.1800 +
3.1801 + }
3.1802 +
3.1803 + if (nca == INVALID) {
3.1804 + if ((*_matching)[left] == INVALID) {
3.1805 + nca = right;
3.1806 + while (left_set.find(nca) == left_set.end()) {
3.1807 + nca = _graph.target((*_matching)[nca]);
3.1808 + right_path.push_back(nca);
3.1809 + nca = _graph.target((*_pred)[nca]);
3.1810 + right_path.push_back(nca);
3.1811 + }
3.1812 + } else {
3.1813 + nca = left;
3.1814 + while (right_set.find(nca) == right_set.end()) {
3.1815 + nca = _graph.target((*_matching)[nca]);
3.1816 + left_path.push_back(nca);
3.1817 + nca = _graph.target((*_pred)[nca]);
3.1818 + left_path.push_back(nca);
3.1819 + }
3.1820 + }
3.1821 + }
3.1822 + }
3.1823 +
3.1824 + alternatePath(nca, tree);
3.1825 + Arc prev;
3.1826 +
3.1827 + prev = _graph.direct(edge, true);
3.1828 + for (int i = 0; left_path[i] != nca; i += 2) {
3.1829 + _matching->set(left_path[i], prev);
3.1830 + _status->set(left_path[i], MATCHED);
3.1831 + evenToMatched(left_path[i], tree);
3.1832 +
3.1833 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
3.1834 + _status->set(left_path[i + 1], MATCHED);
3.1835 + oddToMatched(left_path[i + 1]);
3.1836 + }
3.1837 + _matching->set(nca, prev);
3.1838 +
3.1839 + for (int i = 0; right_path[i] != nca; i += 2) {
3.1840 + _status->set(right_path[i], MATCHED);
3.1841 + evenToMatched(right_path[i], tree);
3.1842 +
3.1843 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
3.1844 + _status->set(right_path[i + 1], MATCHED);
3.1845 + oddToMatched(right_path[i + 1]);
3.1846 + }
3.1847 +
3.1848 + destroyTree(tree);
3.1849 + }
3.1850 +
3.1851 + void extractCycle(const Arc &arc) {
3.1852 + Node left = _graph.source(arc);
3.1853 + Node odd = _graph.target((*_matching)[left]);
3.1854 + Arc prev;
3.1855 + while (odd != left) {
3.1856 + Node even = _graph.target((*_matching)[odd]);
3.1857 + prev = (*_matching)[odd];
3.1858 + odd = _graph.target((*_matching)[even]);
3.1859 + _matching->set(even, _graph.oppositeArc(prev));
3.1860 + }
3.1861 + _matching->set(left, arc);
3.1862 +
3.1863 + Node right = _graph.target(arc);
3.1864 + int right_tree = _tree_set->find(right);
3.1865 + alternatePath(right, right_tree);
3.1866 + destroyTree(right_tree);
3.1867 + _matching->set(right, _graph.oppositeArc(arc));
3.1868 + }
3.1869 +
3.1870 + public:
3.1871 +
3.1872 + /// \brief Constructor
3.1873 + ///
3.1874 + /// Constructor.
3.1875 + MaxWeightedPerfectFractionalMatching(const Graph& graph,
3.1876 + const WeightMap& weight,
3.1877 + bool allow_loops = true)
3.1878 + : _graph(graph), _weight(weight), _matching(0),
3.1879 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
3.1880 + _status(0), _pred(0),
3.1881 + _tree_set_index(0), _tree_set(0),
3.1882 +
3.1883 + _delta2_index(0), _delta2(0),
3.1884 + _delta3_index(0), _delta3(0),
3.1885 +
3.1886 + _delta_sum() {}
3.1887 +
3.1888 + ~MaxWeightedPerfectFractionalMatching() {
3.1889 + destroyStructures();
3.1890 + }
3.1891 +
3.1892 + /// \name Execution Control
3.1893 + /// The simplest way to execute the algorithm is to use the
3.1894 + /// \ref run() member function.
3.1895 +
3.1896 + ///@{
3.1897 +
3.1898 + /// \brief Initialize the algorithm
3.1899 + ///
3.1900 + /// This function initializes the algorithm.
3.1901 + void init() {
3.1902 + createStructures();
3.1903 +
3.1904 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1905 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
3.1906 + }
3.1907 + for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1908 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
3.1909 + }
3.1910 +
3.1911 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.1912 + Value max = - std::numeric_limits<Value>::max();
3.1913 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
3.1914 + if (_graph.target(e) == n && !_allow_loops) continue;
3.1915 + if ((dualScale * _weight[e]) / 2 > max) {
3.1916 + max = (dualScale * _weight[e]) / 2;
3.1917 + }
3.1918 + }
3.1919 + _node_potential->set(n, max);
3.1920 +
3.1921 + _tree_set->insert(n);
3.1922 +
3.1923 + _matching->set(n, INVALID);
3.1924 + _status->set(n, EVEN);
3.1925 + }
3.1926 +
3.1927 + for (EdgeIt e(_graph); e != INVALID; ++e) {
3.1928 + Node left = _graph.u(e);
3.1929 + Node right = _graph.v(e);
3.1930 + if (left == right && !_allow_loops) continue;
3.1931 + _delta3->push(e, ((*_node_potential)[left] +
3.1932 + (*_node_potential)[right] -
3.1933 + dualScale * _weight[e]) / 2);
3.1934 + }
3.1935 + }
3.1936 +
3.1937 + /// \brief Start the algorithm
3.1938 + ///
3.1939 + /// This function starts the algorithm.
3.1940 + ///
3.1941 + /// \pre \ref init() must be called before using this function.
3.1942 + bool start() {
3.1943 + enum OpType {
3.1944 + D2, D3
3.1945 + };
3.1946 +
3.1947 + int unmatched = _node_num;
3.1948 + while (unmatched > 0) {
3.1949 + Value d2 = !_delta2->empty() ?
3.1950 + _delta2->prio() : std::numeric_limits<Value>::max();
3.1951 +
3.1952 + Value d3 = !_delta3->empty() ?
3.1953 + _delta3->prio() : std::numeric_limits<Value>::max();
3.1954 +
3.1955 + _delta_sum = d3; OpType ot = D3;
3.1956 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
3.1957 +
3.1958 + if (_delta_sum == std::numeric_limits<Value>::max()) {
3.1959 + return false;
3.1960 + }
3.1961 +
3.1962 + switch (ot) {
3.1963 + case D2:
3.1964 + {
3.1965 + Node n = _delta2->top();
3.1966 + Arc a = (*_pred)[n];
3.1967 + if ((*_matching)[n] == INVALID) {
3.1968 + augmentOnArc(a);
3.1969 + --unmatched;
3.1970 + } else {
3.1971 + Node v = _graph.target((*_matching)[n]);
3.1972 + if ((*_matching)[n] !=
3.1973 + _graph.oppositeArc((*_matching)[v])) {
3.1974 + extractCycle(a);
3.1975 + --unmatched;
3.1976 + } else {
3.1977 + extendOnArc(a);
3.1978 + }
3.1979 + }
3.1980 + } break;
3.1981 + case D3:
3.1982 + {
3.1983 + Edge e = _delta3->top();
3.1984 +
3.1985 + Node left = _graph.u(e);
3.1986 + Node right = _graph.v(e);
3.1987 +
3.1988 + int left_tree = _tree_set->find(left);
3.1989 + int right_tree = _tree_set->find(right);
3.1990 +
3.1991 + if (left_tree == right_tree) {
3.1992 + cycleOnEdge(e, left_tree);
3.1993 + --unmatched;
3.1994 + } else {
3.1995 + augmentOnEdge(e);
3.1996 + unmatched -= 2;
3.1997 + }
3.1998 + } break;
3.1999 + }
3.2000 + }
3.2001 + return true;
3.2002 + }
3.2003 +
3.2004 + /// \brief Run the algorithm.
3.2005 + ///
3.2006 + /// This method runs the \c %MaxWeightedPerfectFractionalMatching
3.2007 + /// algorithm.
3.2008 + ///
3.2009 + /// \note mwfm.run() is just a shortcut of the following code.
3.2010 + /// \code
3.2011 + /// mwpfm.init();
3.2012 + /// mwpfm.start();
3.2013 + /// \endcode
3.2014 + bool run() {
3.2015 + init();
3.2016 + return start();
3.2017 + }
3.2018 +
3.2019 + /// @}
3.2020 +
3.2021 + /// \name Primal Solution
3.2022 + /// Functions to get the primal solution, i.e. the maximum weighted
3.2023 + /// matching.\n
3.2024 + /// Either \ref run() or \ref start() function should be called before
3.2025 + /// using them.
3.2026 +
3.2027 + /// @{
3.2028 +
3.2029 + /// \brief Return the weight of the matching.
3.2030 + ///
3.2031 + /// This function returns the weight of the found matching. This
3.2032 + /// value is scaled by \ref primalScale "primal scale".
3.2033 + ///
3.2034 + /// \pre Either run() or start() must be called before using this function.
3.2035 + Value matchingWeight() const {
3.2036 + Value sum = 0;
3.2037 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.2038 + if ((*_matching)[n] != INVALID) {
3.2039 + sum += _weight[(*_matching)[n]];
3.2040 + }
3.2041 + }
3.2042 + return sum * primalScale / 2;
3.2043 + }
3.2044 +
3.2045 + /// \brief Return the number of covered nodes in the matching.
3.2046 + ///
3.2047 + /// This function returns the number of covered nodes in the matching.
3.2048 + ///
3.2049 + /// \pre Either run() or start() must be called before using this function.
3.2050 + int matchingSize() const {
3.2051 + int num = 0;
3.2052 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.2053 + if ((*_matching)[n] != INVALID) {
3.2054 + ++num;
3.2055 + }
3.2056 + }
3.2057 + return num;
3.2058 + }
3.2059 +
3.2060 + /// \brief Return \c true if the given edge is in the matching.
3.2061 + ///
3.2062 + /// This function returns \c true if the given edge is in the
3.2063 + /// found matching. The result is scaled by \ref primalScale
3.2064 + /// "primal scale".
3.2065 + ///
3.2066 + /// \pre Either run() or start() must be called before using this function.
3.2067 + int matching(const Edge& edge) const {
3.2068 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
3.2069 + + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
3.2070 + }
3.2071 +
3.2072 + /// \brief Return the fractional matching arc (or edge) incident
3.2073 + /// to the given node.
3.2074 + ///
3.2075 + /// This function returns one of the fractional matching arc (or
3.2076 + /// edge) incident to the given node in the found matching or \c
3.2077 + /// INVALID if the node is not covered by the matching or if the
3.2078 + /// node is on an odd length cycle then it is the successor edge
3.2079 + /// on the cycle.
3.2080 + ///
3.2081 + /// \pre Either run() or start() must be called before using this function.
3.2082 + Arc matching(const Node& node) const {
3.2083 + return (*_matching)[node];
3.2084 + }
3.2085 +
3.2086 + /// \brief Return a const reference to the matching map.
3.2087 + ///
3.2088 + /// This function returns a const reference to a node map that stores
3.2089 + /// the matching arc (or edge) incident to each node.
3.2090 + const MatchingMap& matchingMap() const {
3.2091 + return *_matching;
3.2092 + }
3.2093 +
3.2094 + /// @}
3.2095 +
3.2096 + /// \name Dual Solution
3.2097 + /// Functions to get the dual solution.\n
3.2098 + /// Either \ref run() or \ref start() function should be called before
3.2099 + /// using them.
3.2100 +
3.2101 + /// @{
3.2102 +
3.2103 + /// \brief Return the value of the dual solution.
3.2104 + ///
3.2105 + /// This function returns the value of the dual solution.
3.2106 + /// It should be equal to the primal value scaled by \ref dualScale
3.2107 + /// "dual scale".
3.2108 + ///
3.2109 + /// \pre Either run() or start() must be called before using this function.
3.2110 + Value dualValue() const {
3.2111 + Value sum = 0;
3.2112 + for (NodeIt n(_graph); n != INVALID; ++n) {
3.2113 + sum += nodeValue(n);
3.2114 + }
3.2115 + return sum;
3.2116 + }
3.2117 +
3.2118 + /// \brief Return the dual value (potential) of the given node.
3.2119 + ///
3.2120 + /// This function returns the dual value (potential) of the given node.
3.2121 + ///
3.2122 + /// \pre Either run() or start() must be called before using this function.
3.2123 + Value nodeValue(const Node& n) const {
3.2124 + return (*_node_potential)[n];
3.2125 + }
3.2126 +
3.2127 + /// @}
3.2128 +
3.2129 + };
3.2130 +
3.2131 +} //END OF NAMESPACE LEMON
3.2132 +
3.2133 +#endif //LEMON_FRACTIONAL_MATCHING_H
4.1 --- a/lemon/matching.h Tue Mar 16 21:18:39 2010 +0100
4.2 +++ b/lemon/matching.h Tue Mar 16 21:27:35 2010 +0100
4.3 @@ -16,8 +16,8 @@
4.4 *
4.5 */
4.6
4.7 -#ifndef LEMON_MAX_MATCHING_H
4.8 -#define LEMON_MAX_MATCHING_H
4.9 +#ifndef LEMON_MATCHING_H
4.10 +#define LEMON_MATCHING_H
4.11
4.12 #include <vector>
4.13 #include <queue>
4.14 @@ -28,6 +28,7 @@
4.15 #include <lemon/unionfind.h>
4.16 #include <lemon/bin_heap.h>
4.17 #include <lemon/maps.h>
4.18 +#include <lemon/fractional_matching.h>
4.19
4.20 ///\ingroup matching
4.21 ///\file
4.22 @@ -41,7 +42,7 @@
4.23 ///
4.24 /// This class implements Edmonds' alternating forest matching algorithm
4.25 /// for finding a maximum cardinality matching in a general undirected graph.
4.26 - /// It can be started from an arbitrary initial matching
4.27 + /// It can be started from an arbitrary initial matching
4.28 /// (the default is the empty one).
4.29 ///
4.30 /// The dual solution of the problem is a map of the nodes to
4.31 @@ -69,11 +70,11 @@
4.32
4.33 ///\brief Status constants for Gallai-Edmonds decomposition.
4.34 ///
4.35 - ///These constants are used for indicating the Gallai-Edmonds
4.36 + ///These constants are used for indicating the Gallai-Edmonds
4.37 ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
4.38 ///induce a subgraph with factor-critical components, the nodes with
4.39 ///status \c ODD (or \c A) form the canonical barrier, and the nodes
4.40 - ///with status \c MATCHED (or \c C) induce a subgraph having a
4.41 + ///with status \c MATCHED (or \c C) induce a subgraph having a
4.42 ///perfect matching.
4.43 enum Status {
4.44 EVEN = 1, ///< = 1. (\c D is an alias for \c EVEN.)
4.45 @@ -512,7 +513,7 @@
4.46 }
4.47 }
4.48
4.49 - /// \brief Start Edmonds' algorithm with a heuristic improvement
4.50 + /// \brief Start Edmonds' algorithm with a heuristic improvement
4.51 /// for dense graphs
4.52 ///
4.53 /// This function runs Edmonds' algorithm with a heuristic of postponing
4.54 @@ -534,8 +535,8 @@
4.55
4.56 /// \brief Run Edmonds' algorithm
4.57 ///
4.58 - /// This function runs Edmonds' algorithm. An additional heuristic of
4.59 - /// postponing shrinks is used for relatively dense graphs
4.60 + /// This function runs Edmonds' algorithm. An additional heuristic of
4.61 + /// postponing shrinks is used for relatively dense graphs
4.62 /// (for which <tt>m>=2*n</tt> holds).
4.63 void run() {
4.64 if (countEdges(_graph) < 2 * countNodes(_graph)) {
4.65 @@ -556,7 +557,7 @@
4.66
4.67 /// \brief Return the size (cardinality) of the matching.
4.68 ///
4.69 - /// This function returns the size (cardinality) of the current matching.
4.70 + /// This function returns the size (cardinality) of the current matching.
4.71 /// After run() it returns the size of the maximum matching in the graph.
4.72 int matchingSize() const {
4.73 int size = 0;
4.74 @@ -570,7 +571,7 @@
4.75
4.76 /// \brief Return \c true if the given edge is in the matching.
4.77 ///
4.78 - /// This function returns \c true if the given edge is in the current
4.79 + /// This function returns \c true if the given edge is in the current
4.80 /// matching.
4.81 bool matching(const Edge& edge) const {
4.82 return edge == (*_matching)[_graph.u(edge)];
4.83 @@ -579,7 +580,7 @@
4.84 /// \brief Return the matching arc (or edge) incident to the given node.
4.85 ///
4.86 /// This function returns the matching arc (or edge) incident to the
4.87 - /// given node in the current matching or \c INVALID if the node is
4.88 + /// given node in the current matching or \c INVALID if the node is
4.89 /// not covered by the matching.
4.90 Arc matching(const Node& n) const {
4.91 return (*_matching)[n];
4.92 @@ -595,7 +596,7 @@
4.93
4.94 /// \brief Return the mate of the given node.
4.95 ///
4.96 - /// This function returns the mate of the given node in the current
4.97 + /// This function returns the mate of the given node in the current
4.98 /// matching or \c INVALID if the node is not covered by the matching.
4.99 Node mate(const Node& n) const {
4.100 return (*_matching)[n] != INVALID ?
4.101 @@ -605,7 +606,7 @@
4.102 /// @}
4.103
4.104 /// \name Dual Solution
4.105 - /// Functions to get the dual solution, i.e. the Gallai-Edmonds
4.106 + /// Functions to get the dual solution, i.e. the Gallai-Edmonds
4.107 /// decomposition.
4.108
4.109 /// @{
4.110 @@ -648,8 +649,8 @@
4.111 /// on extensive use of priority queues and provides
4.112 /// \f$O(nm\log n)\f$ time complexity.
4.113 ///
4.114 - /// The maximum weighted matching problem is to find a subset of the
4.115 - /// edges in an undirected graph with maximum overall weight for which
4.116 + /// The maximum weighted matching problem is to find a subset of the
4.117 + /// edges in an undirected graph with maximum overall weight for which
4.118 /// each node has at most one incident edge.
4.119 /// It can be formulated with the following linear program.
4.120 /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
4.121 @@ -673,16 +674,16 @@
4.122 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
4.123 \frac{\vert B \vert - 1}{2}z_B\f] */
4.124 ///
4.125 - /// The algorithm can be executed with the run() function.
4.126 + /// The algorithm can be executed with the run() function.
4.127 /// After it the matching (the primal solution) and the dual solution
4.128 - /// can be obtained using the query functions and the
4.129 - /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
4.130 - /// which is able to iterate on the nodes of a blossom.
4.131 + /// can be obtained using the query functions and the
4.132 + /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
4.133 + /// which is able to iterate on the nodes of a blossom.
4.134 /// If the value type is integer, then the dual solution is multiplied
4.135 /// by \ref MaxWeightedMatching::dualScale "4".
4.136 ///
4.137 /// \tparam GR The undirected graph type the algorithm runs on.
4.138 - /// \tparam WM The type edge weight map. The default type is
4.139 + /// \tparam WM The type edge weight map. The default type is
4.140 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
4.141 #ifdef DOXYGEN
4.142 template <typename GR, typename WM>
4.143 @@ -745,7 +746,7 @@
4.144 typedef RangeMap<int> IntIntMap;
4.145
4.146 enum Status {
4.147 - EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
4.148 + EVEN = -1, MATCHED = 0, ODD = 1
4.149 };
4.150
4.151 typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
4.152 @@ -797,6 +798,10 @@
4.153 BinHeap<Value, IntIntMap> *_delta4;
4.154
4.155 Value _delta_sum;
4.156 + int _unmatched;
4.157 +
4.158 + typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
4.159 + FractionalMatching *_fractional;
4.160
4.161 void createStructures() {
4.162 _node_num = countNodes(_graph);
4.163 @@ -844,9 +849,6 @@
4.164 }
4.165
4.166 void destroyStructures() {
4.167 - _node_num = countNodes(_graph);
4.168 - _blossom_num = _node_num * 3 / 2;
4.169 -
4.170 if (_matching) {
4.171 delete _matching;
4.172 }
4.173 @@ -922,10 +924,6 @@
4.174 if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
4.175 _delta3->push(e, rw / 2);
4.176 }
4.177 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
4.178 - if (_delta3->state(e) != _delta3->IN_HEAP) {
4.179 - _delta3->push(e, rw);
4.180 - }
4.181 } else {
4.182 typename std::map<int, Arc>::iterator it =
4.183 (*_node_data)[vi].heap_index.find(tree);
4.184 @@ -949,202 +947,6 @@
4.185 _delta2->push(vb, _blossom_set->classPrio(vb) -
4.186 (*_blossom_data)[vb].offset);
4.187 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
4.188 - (*_blossom_data)[vb].offset){
4.189 - _delta2->decrease(vb, _blossom_set->classPrio(vb) -
4.190 - (*_blossom_data)[vb].offset);
4.191 - }
4.192 - }
4.193 - }
4.194 - }
4.195 - }
4.196 - }
4.197 - (*_blossom_data)[blossom].offset = 0;
4.198 - }
4.199 -
4.200 - void matchedToOdd(int blossom) {
4.201 - if (_delta2->state(blossom) == _delta2->IN_HEAP) {
4.202 - _delta2->erase(blossom);
4.203 - }
4.204 - (*_blossom_data)[blossom].offset += _delta_sum;
4.205 - if (!_blossom_set->trivial(blossom)) {
4.206 - _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
4.207 - (*_blossom_data)[blossom].offset);
4.208 - }
4.209 - }
4.210 -
4.211 - void evenToMatched(int blossom, int tree) {
4.212 - if (!_blossom_set->trivial(blossom)) {
4.213 - (*_blossom_data)[blossom].pot += 2 * _delta_sum;
4.214 - }
4.215 -
4.216 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
4.217 - n != INVALID; ++n) {
4.218 - int ni = (*_node_index)[n];
4.219 - (*_node_data)[ni].pot -= _delta_sum;
4.220 -
4.221 - _delta1->erase(n);
4.222 -
4.223 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
4.224 - Node v = _graph.source(e);
4.225 - int vb = _blossom_set->find(v);
4.226 - int vi = (*_node_index)[v];
4.227 -
4.228 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.229 - dualScale * _weight[e];
4.230 -
4.231 - if (vb == blossom) {
4.232 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.233 - _delta3->erase(e);
4.234 - }
4.235 - } else if ((*_blossom_data)[vb].status == EVEN) {
4.236 -
4.237 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.238 - _delta3->erase(e);
4.239 - }
4.240 -
4.241 - int vt = _tree_set->find(vb);
4.242 -
4.243 - if (vt != tree) {
4.244 -
4.245 - Arc r = _graph.oppositeArc(e);
4.246 -
4.247 - typename std::map<int, Arc>::iterator it =
4.248 - (*_node_data)[ni].heap_index.find(vt);
4.249 -
4.250 - if (it != (*_node_data)[ni].heap_index.end()) {
4.251 - if ((*_node_data)[ni].heap[it->second] > rw) {
4.252 - (*_node_data)[ni].heap.replace(it->second, r);
4.253 - (*_node_data)[ni].heap.decrease(r, rw);
4.254 - it->second = r;
4.255 - }
4.256 - } else {
4.257 - (*_node_data)[ni].heap.push(r, rw);
4.258 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
4.259 - }
4.260 -
4.261 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
4.262 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
4.263 -
4.264 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
4.265 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
4.266 - (*_blossom_data)[blossom].offset);
4.267 - } else if ((*_delta2)[blossom] >
4.268 - _blossom_set->classPrio(blossom) -
4.269 - (*_blossom_data)[blossom].offset){
4.270 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
4.271 - (*_blossom_data)[blossom].offset);
4.272 - }
4.273 - }
4.274 - }
4.275 -
4.276 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
4.277 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.278 - _delta3->erase(e);
4.279 - }
4.280 - } else {
4.281 -
4.282 - typename std::map<int, Arc>::iterator it =
4.283 - (*_node_data)[vi].heap_index.find(tree);
4.284 -
4.285 - if (it != (*_node_data)[vi].heap_index.end()) {
4.286 - (*_node_data)[vi].heap.erase(it->second);
4.287 - (*_node_data)[vi].heap_index.erase(it);
4.288 - if ((*_node_data)[vi].heap.empty()) {
4.289 - _blossom_set->increase(v, std::numeric_limits<Value>::max());
4.290 - } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
4.291 - _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
4.292 - }
4.293 -
4.294 - if ((*_blossom_data)[vb].status == MATCHED) {
4.295 - if (_blossom_set->classPrio(vb) ==
4.296 - std::numeric_limits<Value>::max()) {
4.297 - _delta2->erase(vb);
4.298 - } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
4.299 - (*_blossom_data)[vb].offset) {
4.300 - _delta2->increase(vb, _blossom_set->classPrio(vb) -
4.301 - (*_blossom_data)[vb].offset);
4.302 - }
4.303 - }
4.304 - }
4.305 - }
4.306 - }
4.307 - }
4.308 - }
4.309 -
4.310 - void oddToMatched(int blossom) {
4.311 - (*_blossom_data)[blossom].offset -= _delta_sum;
4.312 -
4.313 - if (_blossom_set->classPrio(blossom) !=
4.314 - std::numeric_limits<Value>::max()) {
4.315 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
4.316 - (*_blossom_data)[blossom].offset);
4.317 - }
4.318 -
4.319 - if (!_blossom_set->trivial(blossom)) {
4.320 - _delta4->erase(blossom);
4.321 - }
4.322 - }
4.323 -
4.324 - void oddToEven(int blossom, int tree) {
4.325 - if (!_blossom_set->trivial(blossom)) {
4.326 - _delta4->erase(blossom);
4.327 - (*_blossom_data)[blossom].pot -=
4.328 - 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
4.329 - }
4.330 -
4.331 - for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
4.332 - n != INVALID; ++n) {
4.333 - int ni = (*_node_index)[n];
4.334 -
4.335 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
4.336 -
4.337 - (*_node_data)[ni].heap.clear();
4.338 - (*_node_data)[ni].heap_index.clear();
4.339 - (*_node_data)[ni].pot +=
4.340 - 2 * _delta_sum - (*_blossom_data)[blossom].offset;
4.341 -
4.342 - _delta1->push(n, (*_node_data)[ni].pot);
4.343 -
4.344 - for (InArcIt e(_graph, n); e != INVALID; ++e) {
4.345 - Node v = _graph.source(e);
4.346 - int vb = _blossom_set->find(v);
4.347 - int vi = (*_node_index)[v];
4.348 -
4.349 - Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.350 - dualScale * _weight[e];
4.351 -
4.352 - if ((*_blossom_data)[vb].status == EVEN) {
4.353 - if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
4.354 - _delta3->push(e, rw / 2);
4.355 - }
4.356 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
4.357 - if (_delta3->state(e) != _delta3->IN_HEAP) {
4.358 - _delta3->push(e, rw);
4.359 - }
4.360 - } else {
4.361 -
4.362 - typename std::map<int, Arc>::iterator it =
4.363 - (*_node_data)[vi].heap_index.find(tree);
4.364 -
4.365 - if (it != (*_node_data)[vi].heap_index.end()) {
4.366 - if ((*_node_data)[vi].heap[it->second] > rw) {
4.367 - (*_node_data)[vi].heap.replace(it->second, e);
4.368 - (*_node_data)[vi].heap.decrease(e, rw);
4.369 - it->second = e;
4.370 - }
4.371 - } else {
4.372 - (*_node_data)[vi].heap.push(e, rw);
4.373 - (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
4.374 - }
4.375 -
4.376 - if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
4.377 - _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
4.378 -
4.379 - if ((*_blossom_data)[vb].status == MATCHED) {
4.380 - if (_delta2->state(vb) != _delta2->IN_HEAP) {
4.381 - _delta2->push(vb, _blossom_set->classPrio(vb) -
4.382 - (*_blossom_data)[vb].offset);
4.383 - } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
4.384 (*_blossom_data)[vb].offset) {
4.385 _delta2->decrease(vb, _blossom_set->classPrio(vb) -
4.386 (*_blossom_data)[vb].offset);
4.387 @@ -1157,43 +959,145 @@
4.388 (*_blossom_data)[blossom].offset = 0;
4.389 }
4.390
4.391 -
4.392 - void matchedToUnmatched(int blossom) {
4.393 + void matchedToOdd(int blossom) {
4.394 if (_delta2->state(blossom) == _delta2->IN_HEAP) {
4.395 _delta2->erase(blossom);
4.396 }
4.397 + (*_blossom_data)[blossom].offset += _delta_sum;
4.398 + if (!_blossom_set->trivial(blossom)) {
4.399 + _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
4.400 + (*_blossom_data)[blossom].offset);
4.401 + }
4.402 + }
4.403 +
4.404 + void evenToMatched(int blossom, int tree) {
4.405 + if (!_blossom_set->trivial(blossom)) {
4.406 + (*_blossom_data)[blossom].pot += 2 * _delta_sum;
4.407 + }
4.408
4.409 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
4.410 n != INVALID; ++n) {
4.411 int ni = (*_node_index)[n];
4.412 -
4.413 - _blossom_set->increase(n, std::numeric_limits<Value>::max());
4.414 -
4.415 - (*_node_data)[ni].heap.clear();
4.416 - (*_node_data)[ni].heap_index.clear();
4.417 -
4.418 - for (OutArcIt e(_graph, n); e != INVALID; ++e) {
4.419 - Node v = _graph.target(e);
4.420 + (*_node_data)[ni].pot -= _delta_sum;
4.421 +
4.422 + _delta1->erase(n);
4.423 +
4.424 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
4.425 + Node v = _graph.source(e);
4.426 int vb = _blossom_set->find(v);
4.427 int vi = (*_node_index)[v];
4.428
4.429 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.430 dualScale * _weight[e];
4.431
4.432 - if ((*_blossom_data)[vb].status == EVEN) {
4.433 - if (_delta3->state(e) != _delta3->IN_HEAP) {
4.434 - _delta3->push(e, rw);
4.435 + if (vb == blossom) {
4.436 + if (_delta3->state(e) == _delta3->IN_HEAP) {
4.437 + _delta3->erase(e);
4.438 + }
4.439 + } else if ((*_blossom_data)[vb].status == EVEN) {
4.440 +
4.441 + if (_delta3->state(e) == _delta3->IN_HEAP) {
4.442 + _delta3->erase(e);
4.443 + }
4.444 +
4.445 + int vt = _tree_set->find(vb);
4.446 +
4.447 + if (vt != tree) {
4.448 +
4.449 + Arc r = _graph.oppositeArc(e);
4.450 +
4.451 + typename std::map<int, Arc>::iterator it =
4.452 + (*_node_data)[ni].heap_index.find(vt);
4.453 +
4.454 + if (it != (*_node_data)[ni].heap_index.end()) {
4.455 + if ((*_node_data)[ni].heap[it->second] > rw) {
4.456 + (*_node_data)[ni].heap.replace(it->second, r);
4.457 + (*_node_data)[ni].heap.decrease(r, rw);
4.458 + it->second = r;
4.459 + }
4.460 + } else {
4.461 + (*_node_data)[ni].heap.push(r, rw);
4.462 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
4.463 + }
4.464 +
4.465 + if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
4.466 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
4.467 +
4.468 + if (_delta2->state(blossom) != _delta2->IN_HEAP) {
4.469 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
4.470 + (*_blossom_data)[blossom].offset);
4.471 + } else if ((*_delta2)[blossom] >
4.472 + _blossom_set->classPrio(blossom) -
4.473 + (*_blossom_data)[blossom].offset){
4.474 + _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
4.475 + (*_blossom_data)[blossom].offset);
4.476 + }
4.477 + }
4.478 + }
4.479 + } else {
4.480 +
4.481 + typename std::map<int, Arc>::iterator it =
4.482 + (*_node_data)[vi].heap_index.find(tree);
4.483 +
4.484 + if (it != (*_node_data)[vi].heap_index.end()) {
4.485 + (*_node_data)[vi].heap.erase(it->second);
4.486 + (*_node_data)[vi].heap_index.erase(it);
4.487 + if ((*_node_data)[vi].heap.empty()) {
4.488 + _blossom_set->increase(v, std::numeric_limits<Value>::max());
4.489 + } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
4.490 + _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
4.491 + }
4.492 +
4.493 + if ((*_blossom_data)[vb].status == MATCHED) {
4.494 + if (_blossom_set->classPrio(vb) ==
4.495 + std::numeric_limits<Value>::max()) {
4.496 + _delta2->erase(vb);
4.497 + } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
4.498 + (*_blossom_data)[vb].offset) {
4.499 + _delta2->increase(vb, _blossom_set->classPrio(vb) -
4.500 + (*_blossom_data)[vb].offset);
4.501 + }
4.502 + }
4.503 }
4.504 }
4.505 }
4.506 }
4.507 }
4.508
4.509 - void unmatchedToMatched(int blossom) {
4.510 + void oddToMatched(int blossom) {
4.511 + (*_blossom_data)[blossom].offset -= _delta_sum;
4.512 +
4.513 + if (_blossom_set->classPrio(blossom) !=
4.514 + std::numeric_limits<Value>::max()) {
4.515 + _delta2->push(blossom, _blossom_set->classPrio(blossom) -
4.516 + (*_blossom_data)[blossom].offset);
4.517 + }
4.518 +
4.519 + if (!_blossom_set->trivial(blossom)) {
4.520 + _delta4->erase(blossom);
4.521 + }
4.522 + }
4.523 +
4.524 + void oddToEven(int blossom, int tree) {
4.525 + if (!_blossom_set->trivial(blossom)) {
4.526 + _delta4->erase(blossom);
4.527 + (*_blossom_data)[blossom].pot -=
4.528 + 2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
4.529 + }
4.530 +
4.531 for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
4.532 n != INVALID; ++n) {
4.533 int ni = (*_node_index)[n];
4.534
4.535 + _blossom_set->increase(n, std::numeric_limits<Value>::max());
4.536 +
4.537 + (*_node_data)[ni].heap.clear();
4.538 + (*_node_data)[ni].heap_index.clear();
4.539 + (*_node_data)[ni].pot +=
4.540 + 2 * _delta_sum - (*_blossom_data)[blossom].offset;
4.541 +
4.542 + _delta1->push(n, (*_node_data)[ni].pot);
4.543 +
4.544 for (InArcIt e(_graph, n); e != INVALID; ++e) {
4.545 Node v = _graph.source(e);
4.546 int vb = _blossom_set->find(v);
4.547 @@ -1202,54 +1106,44 @@
4.548 Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.549 dualScale * _weight[e];
4.550
4.551 - if (vb == blossom) {
4.552 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.553 - _delta3->erase(e);
4.554 + if ((*_blossom_data)[vb].status == EVEN) {
4.555 + if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
4.556 + _delta3->push(e, rw / 2);
4.557 }
4.558 - } else if ((*_blossom_data)[vb].status == EVEN) {
4.559 -
4.560 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.561 - _delta3->erase(e);
4.562 - }
4.563 -
4.564 - int vt = _tree_set->find(vb);
4.565 -
4.566 - Arc r = _graph.oppositeArc(e);
4.567 + } else {
4.568
4.569 typename std::map<int, Arc>::iterator it =
4.570 - (*_node_data)[ni].heap_index.find(vt);
4.571 -
4.572 - if (it != (*_node_data)[ni].heap_index.end()) {
4.573 - if ((*_node_data)[ni].heap[it->second] > rw) {
4.574 - (*_node_data)[ni].heap.replace(it->second, r);
4.575 - (*_node_data)[ni].heap.decrease(r, rw);
4.576 - it->second = r;
4.577 + (*_node_data)[vi].heap_index.find(tree);
4.578 +
4.579 + if (it != (*_node_data)[vi].heap_index.end()) {
4.580 + if ((*_node_data)[vi].heap[it->second] > rw) {
4.581 + (*_node_data)[vi].heap.replace(it->second, e);
4.582 + (*_node_data)[vi].heap.decrease(e, rw);
4.583 + it->second = e;
4.584 }
4.585 } else {
4.586 - (*_node_data)[ni].heap.push(r, rw);
4.587 - (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
4.588 + (*_node_data)[vi].heap.push(e, rw);
4.589 + (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
4.590 }
4.591
4.592 - if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
4.593 - _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
4.594 -
4.595 - if (_delta2->state(blossom) != _delta2->IN_HEAP) {
4.596 - _delta2->push(blossom, _blossom_set->classPrio(blossom) -
4.597 - (*_blossom_data)[blossom].offset);
4.598 - } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
4.599 - (*_blossom_data)[blossom].offset){
4.600 - _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
4.601 - (*_blossom_data)[blossom].offset);
4.602 + if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
4.603 + _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
4.604 +
4.605 + if ((*_blossom_data)[vb].status == MATCHED) {
4.606 + if (_delta2->state(vb) != _delta2->IN_HEAP) {
4.607 + _delta2->push(vb, _blossom_set->classPrio(vb) -
4.608 + (*_blossom_data)[vb].offset);
4.609 + } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
4.610 + (*_blossom_data)[vb].offset) {
4.611 + _delta2->decrease(vb, _blossom_set->classPrio(vb) -
4.612 + (*_blossom_data)[vb].offset);
4.613 + }
4.614 }
4.615 }
4.616 -
4.617 - } else if ((*_blossom_data)[vb].status == UNMATCHED) {
4.618 - if (_delta3->state(e) == _delta3->IN_HEAP) {
4.619 - _delta3->erase(e);
4.620 - }
4.621 }
4.622 }
4.623 }
4.624 + (*_blossom_data)[blossom].offset = 0;
4.625 }
4.626
4.627 void alternatePath(int even, int tree) {
4.628 @@ -1294,39 +1188,42 @@
4.629 alternatePath(blossom, tree);
4.630 destroyTree(tree);
4.631
4.632 - (*_blossom_data)[blossom].status = UNMATCHED;
4.633 (*_blossom_data)[blossom].base = node;
4.634 - matchedToUnmatched(blossom);
4.635 + (*_blossom_data)[blossom].next = INVALID;
4.636 }
4.637
4.638 -
4.639 void augmentOnEdge(const Edge& edge) {
4.640
4.641 int left = _blossom_set->find(_graph.u(edge));
4.642 int right = _blossom_set->find(_graph.v(edge));
4.643
4.644 - if ((*_blossom_data)[left].status == EVEN) {
4.645 - int left_tree = _tree_set->find(left);
4.646 - alternatePath(left, left_tree);
4.647 - destroyTree(left_tree);
4.648 - } else {
4.649 - (*_blossom_data)[left].status = MATCHED;
4.650 - unmatchedToMatched(left);
4.651 - }
4.652 -
4.653 - if ((*_blossom_data)[right].status == EVEN) {
4.654 - int right_tree = _tree_set->find(right);
4.655 - alternatePath(right, right_tree);
4.656 - destroyTree(right_tree);
4.657 - } else {
4.658 - (*_blossom_data)[right].status = MATCHED;
4.659 - unmatchedToMatched(right);
4.660 - }
4.661 + int left_tree = _tree_set->find(left);
4.662 + alternatePath(left, left_tree);
4.663 + destroyTree(left_tree);
4.664 +
4.665 + int right_tree = _tree_set->find(right);
4.666 + alternatePath(right, right_tree);
4.667 + destroyTree(right_tree);
4.668
4.669 (*_blossom_data)[left].next = _graph.direct(edge, true);
4.670 (*_blossom_data)[right].next = _graph.direct(edge, false);
4.671 }
4.672
4.673 + void augmentOnArc(const Arc& arc) {
4.674 +
4.675 + int left = _blossom_set->find(_graph.source(arc));
4.676 + int right = _blossom_set->find(_graph.target(arc));
4.677 +
4.678 + (*_blossom_data)[left].status = MATCHED;
4.679 +
4.680 + int right_tree = _tree_set->find(right);
4.681 + alternatePath(right, right_tree);
4.682 + destroyTree(right_tree);
4.683 +
4.684 + (*_blossom_data)[left].next = arc;
4.685 + (*_blossom_data)[right].next = _graph.oppositeArc(arc);
4.686 + }
4.687 +
4.688 void extendOnArc(const Arc& arc) {
4.689 int base = _blossom_set->find(_graph.target(arc));
4.690 int tree = _tree_set->find(base);
4.691 @@ -1529,7 +1426,7 @@
4.692 _tree_set->insert(sb, tree);
4.693 (*_blossom_data)[sb].pred = pred;
4.694 (*_blossom_data)[sb].next =
4.695 - _graph.oppositeArc((*_blossom_data)[tb].next);
4.696 + _graph.oppositeArc((*_blossom_data)[tb].next);
4.697
4.698 pred = (*_blossom_data)[ub].next;
4.699
4.700 @@ -1629,7 +1526,7 @@
4.701 }
4.702
4.703 for (int i = 0; i < int(blossoms.size()); ++i) {
4.704 - if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
4.705 + if ((*_blossom_data)[blossoms[i]].next != INVALID) {
4.706
4.707 Value offset = (*_blossom_data)[blossoms[i]].offset;
4.708 (*_blossom_data)[blossoms[i]].pot += 2 * offset;
4.709 @@ -1667,10 +1564,16 @@
4.710 _delta3_index(0), _delta3(0),
4.711 _delta4_index(0), _delta4(0),
4.712
4.713 - _delta_sum() {}
4.714 + _delta_sum(), _unmatched(0),
4.715 +
4.716 + _fractional(0)
4.717 + {}
4.718
4.719 ~MaxWeightedMatching() {
4.720 destroyStructures();
4.721 + if (_fractional) {
4.722 + delete _fractional;
4.723 + }
4.724 }
4.725
4.726 /// \name Execution Control
4.727 @@ -1699,6 +1602,8 @@
4.728 (*_delta4_index)[i] = _delta4->PRE_HEAP;
4.729 }
4.730
4.731 + _unmatched = _node_num;
4.732 +
4.733 int index = 0;
4.734 for (NodeIt n(_graph); n != INVALID; ++n) {
4.735 Value max = 0;
4.736 @@ -1733,18 +1638,155 @@
4.737 }
4.738 }
4.739
4.740 + /// \brief Initialize the algorithm with fractional matching
4.741 + ///
4.742 + /// This function initializes the algorithm with a fractional
4.743 + /// matching. This initialization is also called jumpstart heuristic.
4.744 + void fractionalInit() {
4.745 + createStructures();
4.746 +
4.747 + if (_fractional == 0) {
4.748 + _fractional = new FractionalMatching(_graph, _weight, false);
4.749 + }
4.750 + _fractional->run();
4.751 +
4.752 + for (ArcIt e(_graph); e != INVALID; ++e) {
4.753 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
4.754 + }
4.755 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.756 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
4.757 + }
4.758 + for (EdgeIt e(_graph); e != INVALID; ++e) {
4.759 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
4.760 + }
4.761 + for (int i = 0; i < _blossom_num; ++i) {
4.762 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
4.763 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
4.764 + }
4.765 +
4.766 + _unmatched = 0;
4.767 +
4.768 + int index = 0;
4.769 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.770 + Value pot = _fractional->nodeValue(n);
4.771 + (*_node_index)[n] = index;
4.772 + (*_node_data)[index].pot = pot;
4.773 + int blossom =
4.774 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
4.775 +
4.776 + (*_blossom_data)[blossom].status = MATCHED;
4.777 + (*_blossom_data)[blossom].pred = INVALID;
4.778 + (*_blossom_data)[blossom].next = _fractional->matching(n);
4.779 + if (_fractional->matching(n) == INVALID) {
4.780 + (*_blossom_data)[blossom].base = n;
4.781 + }
4.782 + (*_blossom_data)[blossom].pot = 0;
4.783 + (*_blossom_data)[blossom].offset = 0;
4.784 + ++index;
4.785 + }
4.786 +
4.787 + typename Graph::template NodeMap<bool> processed(_graph, false);
4.788 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.789 + if (processed[n]) continue;
4.790 + processed[n] = true;
4.791 + if (_fractional->matching(n) == INVALID) continue;
4.792 + int num = 1;
4.793 + Node v = _graph.target(_fractional->matching(n));
4.794 + while (n != v) {
4.795 + processed[v] = true;
4.796 + v = _graph.target(_fractional->matching(v));
4.797 + ++num;
4.798 + }
4.799 +
4.800 + if (num % 2 == 1) {
4.801 + std::vector<int> subblossoms(num);
4.802 +
4.803 + subblossoms[--num] = _blossom_set->find(n);
4.804 + _delta1->push(n, _fractional->nodeValue(n));
4.805 + v = _graph.target(_fractional->matching(n));
4.806 + while (n != v) {
4.807 + subblossoms[--num] = _blossom_set->find(v);
4.808 + _delta1->push(v, _fractional->nodeValue(v));
4.809 + v = _graph.target(_fractional->matching(v));
4.810 + }
4.811 +
4.812 + int surface =
4.813 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
4.814 + (*_blossom_data)[surface].status = EVEN;
4.815 + (*_blossom_data)[surface].pred = INVALID;
4.816 + (*_blossom_data)[surface].next = INVALID;
4.817 + (*_blossom_data)[surface].pot = 0;
4.818 + (*_blossom_data)[surface].offset = 0;
4.819 +
4.820 + _tree_set->insert(surface);
4.821 + ++_unmatched;
4.822 + }
4.823 + }
4.824 +
4.825 + for (EdgeIt e(_graph); e != INVALID; ++e) {
4.826 + int si = (*_node_index)[_graph.u(e)];
4.827 + int sb = _blossom_set->find(_graph.u(e));
4.828 + int ti = (*_node_index)[_graph.v(e)];
4.829 + int tb = _blossom_set->find(_graph.v(e));
4.830 + if ((*_blossom_data)[sb].status == EVEN &&
4.831 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
4.832 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
4.833 + dualScale * _weight[e]) / 2);
4.834 + }
4.835 + }
4.836 +
4.837 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.838 + int nb = _blossom_set->find(n);
4.839 + if ((*_blossom_data)[nb].status != MATCHED) continue;
4.840 + int ni = (*_node_index)[n];
4.841 +
4.842 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
4.843 + Node v = _graph.target(e);
4.844 + int vb = _blossom_set->find(v);
4.845 + int vi = (*_node_index)[v];
4.846 +
4.847 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.848 + dualScale * _weight[e];
4.849 +
4.850 + if ((*_blossom_data)[vb].status == EVEN) {
4.851 +
4.852 + int vt = _tree_set->find(vb);
4.853 +
4.854 + typename std::map<int, Arc>::iterator it =
4.855 + (*_node_data)[ni].heap_index.find(vt);
4.856 +
4.857 + if (it != (*_node_data)[ni].heap_index.end()) {
4.858 + if ((*_node_data)[ni].heap[it->second] > rw) {
4.859 + (*_node_data)[ni].heap.replace(it->second, e);
4.860 + (*_node_data)[ni].heap.decrease(e, rw);
4.861 + it->second = e;
4.862 + }
4.863 + } else {
4.864 + (*_node_data)[ni].heap.push(e, rw);
4.865 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
4.866 + }
4.867 + }
4.868 + }
4.869 +
4.870 + if (!(*_node_data)[ni].heap.empty()) {
4.871 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
4.872 + _delta2->push(nb, _blossom_set->classPrio(nb));
4.873 + }
4.874 + }
4.875 + }
4.876 +
4.877 /// \brief Start the algorithm
4.878 ///
4.879 /// This function starts the algorithm.
4.880 ///
4.881 - /// \pre \ref init() must be called before using this function.
4.882 + /// \pre \ref init() or \ref fractionalInit() must be called
4.883 + /// before using this function.
4.884 void start() {
4.885 enum OpType {
4.886 D1, D2, D3, D4
4.887 };
4.888
4.889 - int unmatched = _node_num;
4.890 - while (unmatched > 0) {
4.891 + while (_unmatched > 0) {
4.892 Value d1 = !_delta1->empty() ?
4.893 _delta1->prio() : std::numeric_limits<Value>::max();
4.894
4.895 @@ -1757,26 +1799,30 @@
4.896 Value d4 = !_delta4->empty() ?
4.897 _delta4->prio() : std::numeric_limits<Value>::max();
4.898
4.899 - _delta_sum = d1; OpType ot = D1;
4.900 + _delta_sum = d3; OpType ot = D3;
4.901 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
4.902 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
4.903 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
4.904 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
4.905
4.906 -
4.907 switch (ot) {
4.908 case D1:
4.909 {
4.910 Node n = _delta1->top();
4.911 unmatchNode(n);
4.912 - --unmatched;
4.913 + --_unmatched;
4.914 }
4.915 break;
4.916 case D2:
4.917 {
4.918 int blossom = _delta2->top();
4.919 Node n = _blossom_set->classTop(blossom);
4.920 - Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
4.921 - extendOnArc(e);
4.922 + Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
4.923 + if ((*_blossom_data)[blossom].next == INVALID) {
4.924 + augmentOnArc(a);
4.925 + --_unmatched;
4.926 + } else {
4.927 + extendOnArc(a);
4.928 + }
4.929 }
4.930 break;
4.931 case D3:
4.932 @@ -1789,26 +1835,14 @@
4.933 if (left_blossom == right_blossom) {
4.934 _delta3->pop();
4.935 } else {
4.936 - int left_tree;
4.937 - if ((*_blossom_data)[left_blossom].status == EVEN) {
4.938 - left_tree = _tree_set->find(left_blossom);
4.939 - } else {
4.940 - left_tree = -1;
4.941 - ++unmatched;
4.942 - }
4.943 - int right_tree;
4.944 - if ((*_blossom_data)[right_blossom].status == EVEN) {
4.945 - right_tree = _tree_set->find(right_blossom);
4.946 - } else {
4.947 - right_tree = -1;
4.948 - ++unmatched;
4.949 - }
4.950 + int left_tree = _tree_set->find(left_blossom);
4.951 + int right_tree = _tree_set->find(right_blossom);
4.952
4.953 if (left_tree == right_tree) {
4.954 shrinkOnEdge(e, left_tree);
4.955 } else {
4.956 augmentOnEdge(e);
4.957 - unmatched -= 2;
4.958 + _unmatched -= 2;
4.959 }
4.960 }
4.961 } break;
4.962 @@ -1826,18 +1860,18 @@
4.963 ///
4.964 /// \note mwm.run() is just a shortcut of the following code.
4.965 /// \code
4.966 - /// mwm.init();
4.967 + /// mwm.fractionalInit();
4.968 /// mwm.start();
4.969 /// \endcode
4.970 void run() {
4.971 - init();
4.972 + fractionalInit();
4.973 start();
4.974 }
4.975
4.976 /// @}
4.977
4.978 /// \name Primal Solution
4.979 - /// Functions to get the primal solution, i.e. the maximum weighted
4.980 + /// Functions to get the primal solution, i.e. the maximum weighted
4.981 /// matching.\n
4.982 /// Either \ref run() or \ref start() function should be called before
4.983 /// using them.
4.984 @@ -1856,7 +1890,7 @@
4.985 sum += _weight[(*_matching)[n]];
4.986 }
4.987 }
4.988 - return sum /= 2;
4.989 + return sum / 2;
4.990 }
4.991
4.992 /// \brief Return the size (cardinality) of the matching.
4.993 @@ -1876,7 +1910,7 @@
4.994
4.995 /// \brief Return \c true if the given edge is in the matching.
4.996 ///
4.997 - /// This function returns \c true if the given edge is in the found
4.998 + /// This function returns \c true if the given edge is in the found
4.999 /// matching.
4.1000 ///
4.1001 /// \pre Either run() or start() must be called before using this function.
4.1002 @@ -1887,7 +1921,7 @@
4.1003 /// \brief Return the matching arc (or edge) incident to the given node.
4.1004 ///
4.1005 /// This function returns the matching arc (or edge) incident to the
4.1006 - /// given node in the found matching or \c INVALID if the node is
4.1007 + /// given node in the found matching or \c INVALID if the node is
4.1008 /// not covered by the matching.
4.1009 ///
4.1010 /// \pre Either run() or start() must be called before using this function.
4.1011 @@ -1905,7 +1939,7 @@
4.1012
4.1013 /// \brief Return the mate of the given node.
4.1014 ///
4.1015 - /// This function returns the mate of the given node in the found
4.1016 + /// This function returns the mate of the given node in the found
4.1017 /// matching or \c INVALID if the node is not covered by the matching.
4.1018 ///
4.1019 /// \pre Either run() or start() must be called before using this function.
4.1020 @@ -1925,8 +1959,8 @@
4.1021
4.1022 /// \brief Return the value of the dual solution.
4.1023 ///
4.1024 - /// This function returns the value of the dual solution.
4.1025 - /// It should be equal to the primal value scaled by \ref dualScale
4.1026 + /// This function returns the value of the dual solution.
4.1027 + /// It should be equal to the primal value scaled by \ref dualScale
4.1028 /// "dual scale".
4.1029 ///
4.1030 /// \pre Either run() or start() must be called before using this function.
4.1031 @@ -1981,9 +2015,9 @@
4.1032
4.1033 /// \brief Iterator for obtaining the nodes of a blossom.
4.1034 ///
4.1035 - /// This class provides an iterator for obtaining the nodes of the
4.1036 + /// This class provides an iterator for obtaining the nodes of the
4.1037 /// given blossom. It lists a subset of the nodes.
4.1038 - /// Before using this iterator, you must allocate a
4.1039 + /// Before using this iterator, you must allocate a
4.1040 /// MaxWeightedMatching class and execute it.
4.1041 class BlossomIt {
4.1042 public:
4.1043 @@ -1992,8 +2026,8 @@
4.1044 ///
4.1045 /// Constructor to get the nodes of the given variable.
4.1046 ///
4.1047 - /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
4.1048 - /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
4.1049 + /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
4.1050 + /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
4.1051 /// called before initializing this iterator.
4.1052 BlossomIt(const MaxWeightedMatching& algorithm, int variable)
4.1053 : _algorithm(&algorithm)
4.1054 @@ -2046,8 +2080,8 @@
4.1055 /// is based on extensive use of priority queues and provides
4.1056 /// \f$O(nm\log n)\f$ time complexity.
4.1057 ///
4.1058 - /// The maximum weighted perfect matching problem is to find a subset of
4.1059 - /// the edges in an undirected graph with maximum overall weight for which
4.1060 + /// The maximum weighted perfect matching problem is to find a subset of
4.1061 + /// the edges in an undirected graph with maximum overall weight for which
4.1062 /// each node has exactly one incident edge.
4.1063 /// It can be formulated with the following linear program.
4.1064 /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
4.1065 @@ -2070,16 +2104,16 @@
4.1066 /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
4.1067 \frac{\vert B \vert - 1}{2}z_B\f] */
4.1068 ///
4.1069 - /// The algorithm can be executed with the run() function.
4.1070 + /// The algorithm can be executed with the run() function.
4.1071 /// After it the matching (the primal solution) and the dual solution
4.1072 - /// can be obtained using the query functions and the
4.1073 - /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
4.1074 - /// which is able to iterate on the nodes of a blossom.
4.1075 + /// can be obtained using the query functions and the
4.1076 + /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
4.1077 + /// which is able to iterate on the nodes of a blossom.
4.1078 /// If the value type is integer, then the dual solution is multiplied
4.1079 /// by \ref MaxWeightedMatching::dualScale "4".
4.1080 ///
4.1081 /// \tparam GR The undirected graph type the algorithm runs on.
4.1082 - /// \tparam WM The type edge weight map. The default type is
4.1083 + /// \tparam WM The type edge weight map. The default type is
4.1084 /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
4.1085 #ifdef DOXYGEN
4.1086 template <typename GR, typename WM>
4.1087 @@ -2190,6 +2224,11 @@
4.1088 BinHeap<Value, IntIntMap> *_delta4;
4.1089
4.1090 Value _delta_sum;
4.1091 + int _unmatched;
4.1092 +
4.1093 + typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
4.1094 + FractionalMatching;
4.1095 + FractionalMatching *_fractional;
4.1096
4.1097 void createStructures() {
4.1098 _node_num = countNodes(_graph);
4.1099 @@ -2233,9 +2272,6 @@
4.1100 }
4.1101
4.1102 void destroyStructures() {
4.1103 - _node_num = countNodes(_graph);
4.1104 - _blossom_num = _node_num * 3 / 2;
4.1105 -
4.1106 if (_matching) {
4.1107 delete _matching;
4.1108 }
4.1109 @@ -2908,10 +2944,16 @@
4.1110 _delta3_index(0), _delta3(0),
4.1111 _delta4_index(0), _delta4(0),
4.1112
4.1113 - _delta_sum() {}
4.1114 + _delta_sum(), _unmatched(0),
4.1115 +
4.1116 + _fractional(0)
4.1117 + {}
4.1118
4.1119 ~MaxWeightedPerfectMatching() {
4.1120 destroyStructures();
4.1121 + if (_fractional) {
4.1122 + delete _fractional;
4.1123 + }
4.1124 }
4.1125
4.1126 /// \name Execution Control
4.1127 @@ -2937,6 +2979,8 @@
4.1128 (*_delta4_index)[i] = _delta4->PRE_HEAP;
4.1129 }
4.1130
4.1131 + _unmatched = _node_num;
4.1132 +
4.1133 int index = 0;
4.1134 for (NodeIt n(_graph); n != INVALID; ++n) {
4.1135 Value max = - std::numeric_limits<Value>::max();
4.1136 @@ -2970,18 +3014,152 @@
4.1137 }
4.1138 }
4.1139
4.1140 + /// \brief Initialize the algorithm with fractional matching
4.1141 + ///
4.1142 + /// This function initializes the algorithm with a fractional
4.1143 + /// matching. This initialization is also called jumpstart heuristic.
4.1144 + void fractionalInit() {
4.1145 + createStructures();
4.1146 +
4.1147 + if (_fractional == 0) {
4.1148 + _fractional = new FractionalMatching(_graph, _weight, false);
4.1149 + }
4.1150 + if (!_fractional->run()) {
4.1151 + _unmatched = -1;
4.1152 + return;
4.1153 + }
4.1154 +
4.1155 + for (ArcIt e(_graph); e != INVALID; ++e) {
4.1156 + (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
4.1157 + }
4.1158 + for (EdgeIt e(_graph); e != INVALID; ++e) {
4.1159 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
4.1160 + }
4.1161 + for (int i = 0; i < _blossom_num; ++i) {
4.1162 + (*_delta2_index)[i] = _delta2->PRE_HEAP;
4.1163 + (*_delta4_index)[i] = _delta4->PRE_HEAP;
4.1164 + }
4.1165 +
4.1166 + _unmatched = 0;
4.1167 +
4.1168 + int index = 0;
4.1169 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.1170 + Value pot = _fractional->nodeValue(n);
4.1171 + (*_node_index)[n] = index;
4.1172 + (*_node_data)[index].pot = pot;
4.1173 + int blossom =
4.1174 + _blossom_set->insert(n, std::numeric_limits<Value>::max());
4.1175 +
4.1176 + (*_blossom_data)[blossom].status = MATCHED;
4.1177 + (*_blossom_data)[blossom].pred = INVALID;
4.1178 + (*_blossom_data)[blossom].next = _fractional->matching(n);
4.1179 + (*_blossom_data)[blossom].pot = 0;
4.1180 + (*_blossom_data)[blossom].offset = 0;
4.1181 + ++index;
4.1182 + }
4.1183 +
4.1184 + typename Graph::template NodeMap<bool> processed(_graph, false);
4.1185 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.1186 + if (processed[n]) continue;
4.1187 + processed[n] = true;
4.1188 + if (_fractional->matching(n) == INVALID) continue;
4.1189 + int num = 1;
4.1190 + Node v = _graph.target(_fractional->matching(n));
4.1191 + while (n != v) {
4.1192 + processed[v] = true;
4.1193 + v = _graph.target(_fractional->matching(v));
4.1194 + ++num;
4.1195 + }
4.1196 +
4.1197 + if (num % 2 == 1) {
4.1198 + std::vector<int> subblossoms(num);
4.1199 +
4.1200 + subblossoms[--num] = _blossom_set->find(n);
4.1201 + v = _graph.target(_fractional->matching(n));
4.1202 + while (n != v) {
4.1203 + subblossoms[--num] = _blossom_set->find(v);
4.1204 + v = _graph.target(_fractional->matching(v));
4.1205 + }
4.1206 +
4.1207 + int surface =
4.1208 + _blossom_set->join(subblossoms.begin(), subblossoms.end());
4.1209 + (*_blossom_data)[surface].status = EVEN;
4.1210 + (*_blossom_data)[surface].pred = INVALID;
4.1211 + (*_blossom_data)[surface].next = INVALID;
4.1212 + (*_blossom_data)[surface].pot = 0;
4.1213 + (*_blossom_data)[surface].offset = 0;
4.1214 +
4.1215 + _tree_set->insert(surface);
4.1216 + ++_unmatched;
4.1217 + }
4.1218 + }
4.1219 +
4.1220 + for (EdgeIt e(_graph); e != INVALID; ++e) {
4.1221 + int si = (*_node_index)[_graph.u(e)];
4.1222 + int sb = _blossom_set->find(_graph.u(e));
4.1223 + int ti = (*_node_index)[_graph.v(e)];
4.1224 + int tb = _blossom_set->find(_graph.v(e));
4.1225 + if ((*_blossom_data)[sb].status == EVEN &&
4.1226 + (*_blossom_data)[tb].status == EVEN && sb != tb) {
4.1227 + _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
4.1228 + dualScale * _weight[e]) / 2);
4.1229 + }
4.1230 + }
4.1231 +
4.1232 + for (NodeIt n(_graph); n != INVALID; ++n) {
4.1233 + int nb = _blossom_set->find(n);
4.1234 + if ((*_blossom_data)[nb].status != MATCHED) continue;
4.1235 + int ni = (*_node_index)[n];
4.1236 +
4.1237 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
4.1238 + Node v = _graph.target(e);
4.1239 + int vb = _blossom_set->find(v);
4.1240 + int vi = (*_node_index)[v];
4.1241 +
4.1242 + Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
4.1243 + dualScale * _weight[e];
4.1244 +
4.1245 + if ((*_blossom_data)[vb].status == EVEN) {
4.1246 +
4.1247 + int vt = _tree_set->find(vb);
4.1248 +
4.1249 + typename std::map<int, Arc>::iterator it =
4.1250 + (*_node_data)[ni].heap_index.find(vt);
4.1251 +
4.1252 + if (it != (*_node_data)[ni].heap_index.end()) {
4.1253 + if ((*_node_data)[ni].heap[it->second] > rw) {
4.1254 + (*_node_data)[ni].heap.replace(it->second, e);
4.1255 + (*_node_data)[ni].heap.decrease(e, rw);
4.1256 + it->second = e;
4.1257 + }
4.1258 + } else {
4.1259 + (*_node_data)[ni].heap.push(e, rw);
4.1260 + (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
4.1261 + }
4.1262 + }
4.1263 + }
4.1264 +
4.1265 + if (!(*_node_data)[ni].heap.empty()) {
4.1266 + _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
4.1267 + _delta2->push(nb, _blossom_set->classPrio(nb));
4.1268 + }
4.1269 + }
4.1270 + }
4.1271 +
4.1272 /// \brief Start the algorithm
4.1273 ///
4.1274 /// This function starts the algorithm.
4.1275 ///
4.1276 - /// \pre \ref init() must be called before using this function.
4.1277 + /// \pre \ref init() or \ref fractionalInit() must be called before
4.1278 + /// using this function.
4.1279 bool start() {
4.1280 enum OpType {
4.1281 D2, D3, D4
4.1282 };
4.1283
4.1284 - int unmatched = _node_num;
4.1285 - while (unmatched > 0) {
4.1286 + if (_unmatched == -1) return false;
4.1287 +
4.1288 + while (_unmatched > 0) {
4.1289 Value d2 = !_delta2->empty() ?
4.1290 _delta2->prio() : std::numeric_limits<Value>::max();
4.1291
4.1292 @@ -2991,8 +3169,8 @@
4.1293 Value d4 = !_delta4->empty() ?
4.1294 _delta4->prio() : std::numeric_limits<Value>::max();
4.1295
4.1296 - _delta_sum = d2; OpType ot = D2;
4.1297 - if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
4.1298 + _delta_sum = d3; OpType ot = D3;
4.1299 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
4.1300 if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
4.1301
4.1302 if (_delta_sum == std::numeric_limits<Value>::max()) {
4.1303 @@ -3025,7 +3203,7 @@
4.1304 shrinkOnEdge(e, left_tree);
4.1305 } else {
4.1306 augmentOnEdge(e);
4.1307 - unmatched -= 2;
4.1308 + _unmatched -= 2;
4.1309 }
4.1310 }
4.1311 } break;
4.1312 @@ -3044,18 +3222,18 @@
4.1313 ///
4.1314 /// \note mwpm.run() is just a shortcut of the following code.
4.1315 /// \code
4.1316 - /// mwpm.init();
4.1317 + /// mwpm.fractionalInit();
4.1318 /// mwpm.start();
4.1319 /// \endcode
4.1320 bool run() {
4.1321 - init();
4.1322 + fractionalInit();
4.1323 return start();
4.1324 }
4.1325
4.1326 /// @}
4.1327
4.1328 /// \name Primal Solution
4.1329 - /// Functions to get the primal solution, i.e. the maximum weighted
4.1330 + /// Functions to get the primal solution, i.e. the maximum weighted
4.1331 /// perfect matching.\n
4.1332 /// Either \ref run() or \ref start() function should be called before
4.1333 /// using them.
4.1334 @@ -3074,12 +3252,12 @@
4.1335 sum += _weight[(*_matching)[n]];
4.1336 }
4.1337 }
4.1338 - return sum /= 2;
4.1339 + return sum / 2;
4.1340 }
4.1341
4.1342 /// \brief Return \c true if the given edge is in the matching.
4.1343 ///
4.1344 - /// This function returns \c true if the given edge is in the found
4.1345 + /// This function returns \c true if the given edge is in the found
4.1346 /// matching.
4.1347 ///
4.1348 /// \pre Either run() or start() must be called before using this function.
4.1349 @@ -3090,7 +3268,7 @@
4.1350 /// \brief Return the matching arc (or edge) incident to the given node.
4.1351 ///
4.1352 /// This function returns the matching arc (or edge) incident to the
4.1353 - /// given node in the found matching or \c INVALID if the node is
4.1354 + /// given node in the found matching or \c INVALID if the node is
4.1355 /// not covered by the matching.
4.1356 ///
4.1357 /// \pre Either run() or start() must be called before using this function.
4.1358 @@ -3108,7 +3286,7 @@
4.1359
4.1360 /// \brief Return the mate of the given node.
4.1361 ///
4.1362 - /// This function returns the mate of the given node in the found
4.1363 + /// This function returns the mate of the given node in the found
4.1364 /// matching or \c INVALID if the node is not covered by the matching.
4.1365 ///
4.1366 /// \pre Either run() or start() must be called before using this function.
4.1367 @@ -3127,8 +3305,8 @@
4.1368
4.1369 /// \brief Return the value of the dual solution.
4.1370 ///
4.1371 - /// This function returns the value of the dual solution.
4.1372 - /// It should be equal to the primal value scaled by \ref dualScale
4.1373 + /// This function returns the value of the dual solution.
4.1374 + /// It should be equal to the primal value scaled by \ref dualScale
4.1375 /// "dual scale".
4.1376 ///
4.1377 /// \pre Either run() or start() must be called before using this function.
4.1378 @@ -3183,9 +3361,9 @@
4.1379
4.1380 /// \brief Iterator for obtaining the nodes of a blossom.
4.1381 ///
4.1382 - /// This class provides an iterator for obtaining the nodes of the
4.1383 + /// This class provides an iterator for obtaining the nodes of the
4.1384 /// given blossom. It lists a subset of the nodes.
4.1385 - /// Before using this iterator, you must allocate a
4.1386 + /// Before using this iterator, you must allocate a
4.1387 /// MaxWeightedPerfectMatching class and execute it.
4.1388 class BlossomIt {
4.1389 public:
4.1390 @@ -3194,8 +3372,8 @@
4.1391 ///
4.1392 /// Constructor to get the nodes of the given variable.
4.1393 ///
4.1394 - /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
4.1395 - /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
4.1396 + /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
4.1397 + /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
4.1398 /// must be called before initializing this iterator.
4.1399 BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
4.1400 : _algorithm(&algorithm)
4.1401 @@ -3241,4 +3419,4 @@
4.1402
4.1403 } //END OF NAMESPACE LEMON
4.1404
4.1405 -#endif //LEMON_MAX_MATCHING_H
4.1406 +#endif //LEMON_MATCHING_H
5.1 --- a/test/CMakeLists.txt Tue Mar 16 21:18:39 2010 +0100
5.2 +++ b/test/CMakeLists.txt Tue Mar 16 21:27:35 2010 +0100
5.3 @@ -21,6 +21,7 @@
5.4 edge_set_test
5.5 error_test
5.6 euler_test
5.7 + fractional_matching_test
5.8 gomory_hu_test
5.9 graph_copy_test
5.10 graph_test
6.1 --- a/test/Makefile.am Tue Mar 16 21:18:39 2010 +0100
6.2 +++ b/test/Makefile.am Tue Mar 16 21:27:35 2010 +0100
6.3 @@ -23,6 +23,7 @@
6.4 test/edge_set_test \
6.5 test/error_test \
6.6 test/euler_test \
6.7 + test/fractional_matching_test \
6.8 test/gomory_hu_test \
6.9 test/graph_copy_test \
6.10 test/graph_test \
6.11 @@ -71,6 +72,7 @@
6.12 test_edge_set_test_SOURCES = test/edge_set_test.cc
6.13 test_error_test_SOURCES = test/error_test.cc
6.14 test_euler_test_SOURCES = test/euler_test.cc
6.15 +test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
6.16 test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
6.17 test_graph_copy_test_SOURCES = test/graph_copy_test.cc
6.18 test_graph_test_SOURCES = test/graph_test.cc
7.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
7.2 +++ b/test/fractional_matching_test.cc Tue Mar 16 21:27:35 2010 +0100
7.3 @@ -0,0 +1,525 @@
7.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
7.5 + *
7.6 + * This file is a part of LEMON, a generic C++ optimization library.
7.7 + *
7.8 + * Copyright (C) 2003-2009
7.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
7.11 + *
7.12 + * Permission to use, modify and distribute this software is granted
7.13 + * provided that this copyright notice appears in all copies. For
7.14 + * precise terms see the accompanying LICENSE file.
7.15 + *
7.16 + * This software is provided "AS IS" with no warranty of any kind,
7.17 + * express or implied, and with no claim as to its suitability for any
7.18 + * purpose.
7.19 + *
7.20 + */
7.21 +
7.22 +#include <iostream>
7.23 +#include <sstream>
7.24 +#include <vector>
7.25 +#include <queue>
7.26 +#include <cstdlib>
7.27 +
7.28 +#include <lemon/fractional_matching.h>
7.29 +#include <lemon/smart_graph.h>
7.30 +#include <lemon/concepts/graph.h>
7.31 +#include <lemon/concepts/maps.h>
7.32 +#include <lemon/lgf_reader.h>
7.33 +#include <lemon/math.h>
7.34 +
7.35 +#include "test_tools.h"
7.36 +
7.37 +using namespace std;
7.38 +using namespace lemon;
7.39 +
7.40 +GRAPH_TYPEDEFS(SmartGraph);
7.41 +
7.42 +
7.43 +const int lgfn = 4;
7.44 +const std::string lgf[lgfn] = {
7.45 + "@nodes\n"
7.46 + "label\n"
7.47 + "0\n"
7.48 + "1\n"
7.49 + "2\n"
7.50 + "3\n"
7.51 + "4\n"
7.52 + "5\n"
7.53 + "6\n"
7.54 + "7\n"
7.55 + "@edges\n"
7.56 + " label weight\n"
7.57 + "7 4 0 984\n"
7.58 + "0 7 1 73\n"
7.59 + "7 1 2 204\n"
7.60 + "2 3 3 583\n"
7.61 + "2 7 4 565\n"
7.62 + "2 1 5 582\n"
7.63 + "0 4 6 551\n"
7.64 + "2 5 7 385\n"
7.65 + "1 5 8 561\n"
7.66 + "5 3 9 484\n"
7.67 + "7 5 10 904\n"
7.68 + "3 6 11 47\n"
7.69 + "7 6 12 888\n"
7.70 + "3 0 13 747\n"
7.71 + "6 1 14 310\n",
7.72 +
7.73 + "@nodes\n"
7.74 + "label\n"
7.75 + "0\n"
7.76 + "1\n"
7.77 + "2\n"
7.78 + "3\n"
7.79 + "4\n"
7.80 + "5\n"
7.81 + "6\n"
7.82 + "7\n"
7.83 + "@edges\n"
7.84 + " label weight\n"
7.85 + "2 5 0 710\n"
7.86 + "0 5 1 241\n"
7.87 + "2 4 2 856\n"
7.88 + "2 6 3 762\n"
7.89 + "4 1 4 747\n"
7.90 + "6 1 5 962\n"
7.91 + "4 7 6 723\n"
7.92 + "1 7 7 661\n"
7.93 + "2 3 8 376\n"
7.94 + "1 0 9 416\n"
7.95 + "6 7 10 391\n",
7.96 +
7.97 + "@nodes\n"
7.98 + "label\n"
7.99 + "0\n"
7.100 + "1\n"
7.101 + "2\n"
7.102 + "3\n"
7.103 + "4\n"
7.104 + "5\n"
7.105 + "6\n"
7.106 + "7\n"
7.107 + "@edges\n"
7.108 + " label weight\n"
7.109 + "6 2 0 553\n"
7.110 + "0 7 1 653\n"
7.111 + "6 3 2 22\n"
7.112 + "4 7 3 846\n"
7.113 + "7 2 4 981\n"
7.114 + "7 6 5 250\n"
7.115 + "5 2 6 539\n",
7.116 +
7.117 + "@nodes\n"
7.118 + "label\n"
7.119 + "0\n"
7.120 + "@edges\n"
7.121 + " label weight\n"
7.122 + "0 0 0 100\n"
7.123 +};
7.124 +
7.125 +void checkMaxFractionalMatchingCompile()
7.126 +{
7.127 + typedef concepts::Graph Graph;
7.128 + typedef Graph::Node Node;
7.129 + typedef Graph::Edge Edge;
7.130 +
7.131 + Graph g;
7.132 + Node n;
7.133 + Edge e;
7.134 +
7.135 + MaxFractionalMatching<Graph> mat_test(g);
7.136 + const MaxFractionalMatching<Graph>&
7.137 + const_mat_test = mat_test;
7.138 +
7.139 + mat_test.init();
7.140 + mat_test.start();
7.141 + mat_test.start(true);
7.142 + mat_test.startPerfect();
7.143 + mat_test.startPerfect(true);
7.144 + mat_test.run();
7.145 + mat_test.run(true);
7.146 + mat_test.runPerfect();
7.147 + mat_test.runPerfect(true);
7.148 +
7.149 + const_mat_test.matchingSize();
7.150 + const_mat_test.matching(e);
7.151 + const_mat_test.matching(n);
7.152 + const MaxFractionalMatching<Graph>::MatchingMap& mmap =
7.153 + const_mat_test.matchingMap();
7.154 + e = mmap[n];
7.155 +
7.156 + const_mat_test.barrier(n);
7.157 +}
7.158 +
7.159 +void checkMaxWeightedFractionalMatchingCompile()
7.160 +{
7.161 + typedef concepts::Graph Graph;
7.162 + typedef Graph::Node Node;
7.163 + typedef Graph::Edge Edge;
7.164 + typedef Graph::EdgeMap<int> WeightMap;
7.165 +
7.166 + Graph g;
7.167 + Node n;
7.168 + Edge e;
7.169 + WeightMap w(g);
7.170 +
7.171 + MaxWeightedFractionalMatching<Graph> mat_test(g, w);
7.172 + const MaxWeightedFractionalMatching<Graph>&
7.173 + const_mat_test = mat_test;
7.174 +
7.175 + mat_test.init();
7.176 + mat_test.start();
7.177 + mat_test.run();
7.178 +
7.179 + const_mat_test.matchingWeight();
7.180 + const_mat_test.matchingSize();
7.181 + const_mat_test.matching(e);
7.182 + const_mat_test.matching(n);
7.183 + const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
7.184 + const_mat_test.matchingMap();
7.185 + e = mmap[n];
7.186 +
7.187 + const_mat_test.dualValue();
7.188 + const_mat_test.nodeValue(n);
7.189 +}
7.190 +
7.191 +void checkMaxWeightedPerfectFractionalMatchingCompile()
7.192 +{
7.193 + typedef concepts::Graph Graph;
7.194 + typedef Graph::Node Node;
7.195 + typedef Graph::Edge Edge;
7.196 + typedef Graph::EdgeMap<int> WeightMap;
7.197 +
7.198 + Graph g;
7.199 + Node n;
7.200 + Edge e;
7.201 + WeightMap w(g);
7.202 +
7.203 + MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
7.204 + const MaxWeightedPerfectFractionalMatching<Graph>&
7.205 + const_mat_test = mat_test;
7.206 +
7.207 + mat_test.init();
7.208 + mat_test.start();
7.209 + mat_test.run();
7.210 +
7.211 + const_mat_test.matchingWeight();
7.212 + const_mat_test.matching(e);
7.213 + const_mat_test.matching(n);
7.214 + const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
7.215 + const_mat_test.matchingMap();
7.216 + e = mmap[n];
7.217 +
7.218 + const_mat_test.dualValue();
7.219 + const_mat_test.nodeValue(n);
7.220 +}
7.221 +
7.222 +void checkFractionalMatching(const SmartGraph& graph,
7.223 + const MaxFractionalMatching<SmartGraph>& mfm,
7.224 + bool allow_loops = true) {
7.225 + int pv = 0;
7.226 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.227 + int indeg = 0;
7.228 + for (InArcIt a(graph, n); a != INVALID; ++a) {
7.229 + if (mfm.matching(graph.source(a)) == a) {
7.230 + ++indeg;
7.231 + }
7.232 + }
7.233 + if (mfm.matching(n) != INVALID) {
7.234 + check(indeg == 1, "Invalid matching");
7.235 + ++pv;
7.236 + } else {
7.237 + check(indeg == 0, "Invalid matching");
7.238 + }
7.239 + }
7.240 + check(pv == mfm.matchingSize(), "Wrong matching size");
7.241 +
7.242 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.243 + check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
7.244 + (e == mfm.matching(graph.v(e)) ? 1 : 0) ==
7.245 + mfm.matching(e), "Invalid matching");
7.246 + }
7.247 +
7.248 + SmartGraph::NodeMap<bool> processed(graph, false);
7.249 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.250 + if (processed[n]) continue;
7.251 + processed[n] = true;
7.252 + if (mfm.matching(n) == INVALID) continue;
7.253 + int num = 1;
7.254 + Node v = graph.target(mfm.matching(n));
7.255 + while (v != n) {
7.256 + processed[v] = true;
7.257 + ++num;
7.258 + v = graph.target(mfm.matching(v));
7.259 + }
7.260 + check(num == 2 || num % 2 == 1, "Wrong cycle size");
7.261 + check(allow_loops || num != 1, "Wrong cycle size");
7.262 + }
7.263 +
7.264 + int anum = 0, bnum = 0;
7.265 + SmartGraph::NodeMap<bool> neighbours(graph, false);
7.266 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.267 + if (!mfm.barrier(n)) continue;
7.268 + ++anum;
7.269 + for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
7.270 + Node u = graph.source(a);
7.271 + if (!allow_loops && u == n) continue;
7.272 + if (!neighbours[u]) {
7.273 + neighbours[u] = true;
7.274 + ++bnum;
7.275 + }
7.276 + }
7.277 + }
7.278 + check(anum - bnum + mfm.matchingSize() == countNodes(graph),
7.279 + "Wrong barrier");
7.280 +}
7.281 +
7.282 +void checkPerfectFractionalMatching(const SmartGraph& graph,
7.283 + const MaxFractionalMatching<SmartGraph>& mfm,
7.284 + bool perfect, bool allow_loops = true) {
7.285 + if (perfect) {
7.286 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.287 + int indeg = 0;
7.288 + for (InArcIt a(graph, n); a != INVALID; ++a) {
7.289 + if (mfm.matching(graph.source(a)) == a) {
7.290 + ++indeg;
7.291 + }
7.292 + }
7.293 + check(mfm.matching(n) != INVALID, "Invalid matching");
7.294 + check(indeg == 1, "Invalid matching");
7.295 + }
7.296 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.297 + check((e == mfm.matching(graph.u(e)) ? 1 : 0) +
7.298 + (e == mfm.matching(graph.v(e)) ? 1 : 0) ==
7.299 + mfm.matching(e), "Invalid matching");
7.300 + }
7.301 + } else {
7.302 + int anum = 0, bnum = 0;
7.303 + SmartGraph::NodeMap<bool> neighbours(graph, false);
7.304 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.305 + if (!mfm.barrier(n)) continue;
7.306 + ++anum;
7.307 + for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
7.308 + Node u = graph.source(a);
7.309 + if (!allow_loops && u == n) continue;
7.310 + if (!neighbours[u]) {
7.311 + neighbours[u] = true;
7.312 + ++bnum;
7.313 + }
7.314 + }
7.315 + }
7.316 + check(anum - bnum > 0, "Wrong barrier");
7.317 + }
7.318 +}
7.319 +
7.320 +void checkWeightedFractionalMatching(const SmartGraph& graph,
7.321 + const SmartGraph::EdgeMap<int>& weight,
7.322 + const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
7.323 + bool allow_loops = true) {
7.324 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.325 + if (graph.u(e) == graph.v(e) && !allow_loops) continue;
7.326 + int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
7.327 + - weight[e] * mwfm.dualScale;
7.328 +
7.329 + check(rw >= 0, "Negative reduced weight");
7.330 + check(rw == 0 || !mwfm.matching(e),
7.331 + "Non-zero reduced weight on matching edge");
7.332 + }
7.333 +
7.334 + int pv = 0;
7.335 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.336 + int indeg = 0;
7.337 + for (InArcIt a(graph, n); a != INVALID; ++a) {
7.338 + if (mwfm.matching(graph.source(a)) == a) {
7.339 + ++indeg;
7.340 + }
7.341 + }
7.342 + check(indeg <= 1, "Invalid matching");
7.343 + if (mwfm.matching(n) != INVALID) {
7.344 + check(mwfm.nodeValue(n) >= 0, "Invalid node value");
7.345 + check(indeg == 1, "Invalid matching");
7.346 + pv += weight[mwfm.matching(n)];
7.347 + SmartGraph::Node o = graph.target(mwfm.matching(n));
7.348 + } else {
7.349 + check(mwfm.nodeValue(n) == 0, "Invalid matching");
7.350 + check(indeg == 0, "Invalid matching");
7.351 + }
7.352 + }
7.353 +
7.354 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.355 + check((e == mwfm.matching(graph.u(e)) ? 1 : 0) +
7.356 + (e == mwfm.matching(graph.v(e)) ? 1 : 0) ==
7.357 + mwfm.matching(e), "Invalid matching");
7.358 + }
7.359 +
7.360 + int dv = 0;
7.361 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.362 + dv += mwfm.nodeValue(n);
7.363 + }
7.364 +
7.365 + check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
7.366 +
7.367 + SmartGraph::NodeMap<bool> processed(graph, false);
7.368 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.369 + if (processed[n]) continue;
7.370 + processed[n] = true;
7.371 + if (mwfm.matching(n) == INVALID) continue;
7.372 + int num = 1;
7.373 + Node v = graph.target(mwfm.matching(n));
7.374 + while (v != n) {
7.375 + processed[v] = true;
7.376 + ++num;
7.377 + v = graph.target(mwfm.matching(v));
7.378 + }
7.379 + check(num == 2 || num % 2 == 1, "Wrong cycle size");
7.380 + check(allow_loops || num != 1, "Wrong cycle size");
7.381 + }
7.382 +
7.383 + return;
7.384 +}
7.385 +
7.386 +void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
7.387 + const SmartGraph::EdgeMap<int>& weight,
7.388 + const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
7.389 + bool allow_loops = true) {
7.390 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.391 + if (graph.u(e) == graph.v(e) && !allow_loops) continue;
7.392 + int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
7.393 + - weight[e] * mwpfm.dualScale;
7.394 +
7.395 + check(rw >= 0, "Negative reduced weight");
7.396 + check(rw == 0 || !mwpfm.matching(e),
7.397 + "Non-zero reduced weight on matching edge");
7.398 + }
7.399 +
7.400 + int pv = 0;
7.401 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.402 + int indeg = 0;
7.403 + for (InArcIt a(graph, n); a != INVALID; ++a) {
7.404 + if (mwpfm.matching(graph.source(a)) == a) {
7.405 + ++indeg;
7.406 + }
7.407 + }
7.408 + check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
7.409 + check(indeg == 1, "Invalid perfect matching");
7.410 + pv += weight[mwpfm.matching(n)];
7.411 + SmartGraph::Node o = graph.target(mwpfm.matching(n));
7.412 + }
7.413 +
7.414 + for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
7.415 + check((e == mwpfm.matching(graph.u(e)) ? 1 : 0) +
7.416 + (e == mwpfm.matching(graph.v(e)) ? 1 : 0) ==
7.417 + mwpfm.matching(e), "Invalid matching");
7.418 + }
7.419 +
7.420 + int dv = 0;
7.421 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.422 + dv += mwpfm.nodeValue(n);
7.423 + }
7.424 +
7.425 + check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
7.426 +
7.427 + SmartGraph::NodeMap<bool> processed(graph, false);
7.428 + for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
7.429 + if (processed[n]) continue;
7.430 + processed[n] = true;
7.431 + if (mwpfm.matching(n) == INVALID) continue;
7.432 + int num = 1;
7.433 + Node v = graph.target(mwpfm.matching(n));
7.434 + while (v != n) {
7.435 + processed[v] = true;
7.436 + ++num;
7.437 + v = graph.target(mwpfm.matching(v));
7.438 + }
7.439 + check(num == 2 || num % 2 == 1, "Wrong cycle size");
7.440 + check(allow_loops || num != 1, "Wrong cycle size");
7.441 + }
7.442 +
7.443 + return;
7.444 +}
7.445 +
7.446 +
7.447 +int main() {
7.448 +
7.449 + for (int i = 0; i < lgfn; ++i) {
7.450 + SmartGraph graph;
7.451 + SmartGraph::EdgeMap<int> weight(graph);
7.452 +
7.453 + istringstream lgfs(lgf[i]);
7.454 + graphReader(graph, lgfs).
7.455 + edgeMap("weight", weight).run();
7.456 +
7.457 + bool perfect_with_loops;
7.458 + {
7.459 + MaxFractionalMatching<SmartGraph> mfm(graph, true);
7.460 + mfm.run();
7.461 + checkFractionalMatching(graph, mfm, true);
7.462 + perfect_with_loops = mfm.matchingSize() == countNodes(graph);
7.463 + }
7.464 +
7.465 + bool perfect_without_loops;
7.466 + {
7.467 + MaxFractionalMatching<SmartGraph> mfm(graph, false);
7.468 + mfm.run();
7.469 + checkFractionalMatching(graph, mfm, false);
7.470 + perfect_without_loops = mfm.matchingSize() == countNodes(graph);
7.471 + }
7.472 +
7.473 + {
7.474 + MaxFractionalMatching<SmartGraph> mfm(graph, true);
7.475 + bool result = mfm.runPerfect();
7.476 + checkPerfectFractionalMatching(graph, mfm, result, true);
7.477 + check(result == perfect_with_loops, "Wrong perfect matching");
7.478 + }
7.479 +
7.480 + {
7.481 + MaxFractionalMatching<SmartGraph> mfm(graph, false);
7.482 + bool result = mfm.runPerfect();
7.483 + checkPerfectFractionalMatching(graph, mfm, result, false);
7.484 + check(result == perfect_without_loops, "Wrong perfect matching");
7.485 + }
7.486 +
7.487 + {
7.488 + MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
7.489 + mwfm.run();
7.490 + checkWeightedFractionalMatching(graph, weight, mwfm, true);
7.491 + }
7.492 +
7.493 + {
7.494 + MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
7.495 + mwfm.run();
7.496 + checkWeightedFractionalMatching(graph, weight, mwfm, false);
7.497 + }
7.498 +
7.499 + {
7.500 + MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
7.501 + true);
7.502 + bool perfect = mwpfm.run();
7.503 + check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
7.504 + "Perfect matching found");
7.505 + check(perfect == perfect_with_loops, "Wrong perfect matching");
7.506 +
7.507 + if (perfect) {
7.508 + checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
7.509 + }
7.510 + }
7.511 +
7.512 + {
7.513 + MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
7.514 + false);
7.515 + bool perfect = mwpfm.run();
7.516 + check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
7.517 + "Perfect matching found");
7.518 + check(perfect == perfect_without_loops, "Wrong perfect matching");
7.519 +
7.520 + if (perfect) {
7.521 + checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
7.522 + }
7.523 + }
7.524 +
7.525 + }
7.526 +
7.527 + return 0;
7.528 +}
8.1 --- a/test/matching_test.cc Tue Mar 16 21:18:39 2010 +0100
8.2 +++ b/test/matching_test.cc Tue Mar 16 21:27:35 2010 +0100
8.3 @@ -401,22 +401,46 @@
8.4 graphReader(graph, lgfs).
8.5 edgeMap("weight", weight).run();
8.6
8.7 - MaxMatching<SmartGraph> mm(graph);
8.8 - mm.run();
8.9 - checkMatching(graph, mm);
8.10 + bool perfect;
8.11 + {
8.12 + MaxMatching<SmartGraph> mm(graph);
8.13 + mm.run();
8.14 + checkMatching(graph, mm);
8.15 + perfect = 2 * mm.matchingSize() == countNodes(graph);
8.16 + }
8.17
8.18 - MaxWeightedMatching<SmartGraph> mwm(graph, weight);
8.19 - mwm.run();
8.20 - checkWeightedMatching(graph, weight, mwm);
8.21 + {
8.22 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
8.23 + mwm.run();
8.24 + checkWeightedMatching(graph, weight, mwm);
8.25 + }
8.26
8.27 - MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
8.28 - bool perfect = mwpm.run();
8.29 + {
8.30 + MaxWeightedMatching<SmartGraph> mwm(graph, weight);
8.31 + mwm.init();
8.32 + mwm.start();
8.33 + checkWeightedMatching(graph, weight, mwm);
8.34 + }
8.35
8.36 - check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
8.37 - "Perfect matching found");
8.38 + {
8.39 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
8.40 + bool result = mwpm.run();
8.41 +
8.42 + check(result == perfect, "Perfect matching found");
8.43 + if (perfect) {
8.44 + checkWeightedPerfectMatching(graph, weight, mwpm);
8.45 + }
8.46 + }
8.47
8.48 - if (perfect) {
8.49 - checkWeightedPerfectMatching(graph, weight, mwpm);
8.50 + {
8.51 + MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
8.52 + mwpm.init();
8.53 + bool result = mwpm.start();
8.54 +
8.55 + check(result == perfect, "Perfect matching found");
8.56 + if (perfect) {
8.57 + checkWeightedPerfectMatching(graph, weight, mwpm);
8.58 + }
8.59 }
8.60 }
8.61