1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/fractional_matching.h Fri Sep 25 21:51:36 2009 +0200
1.3 @@ -0,0 +1,2135 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2009
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_FRACTIONAL_MATCHING_H
1.23 +#define LEMON_FRACTIONAL_MATCHING_H
1.24 +
1.25 +#include <vector>
1.26 +#include <queue>
1.27 +#include <set>
1.28 +#include <limits>
1.29 +
1.30 +#include <lemon/core.h>
1.31 +#include <lemon/unionfind.h>
1.32 +#include <lemon/bin_heap.h>
1.33 +#include <lemon/maps.h>
1.34 +#include <lemon/assert.h>
1.35 +#include <lemon/elevator.h>
1.36 +
1.37 +///\ingroup matching
1.38 +///\file
1.39 +///\brief Fractional matching algorithms in general graphs.
1.40 +
1.41 +namespace lemon {
1.42 +
1.43 + /// \brief Default traits class of MaxFractionalMatching class.
1.44 + ///
1.45 + /// Default traits class of MaxFractionalMatching class.
1.46 + /// \tparam GR Graph type.
1.47 + template <typename GR>
1.48 + struct MaxFractionalMatchingDefaultTraits {
1.49 +
1.50 + /// \brief The type of the graph the algorithm runs on.
1.51 + typedef GR Graph;
1.52 +
1.53 + /// \brief The type of the map that stores the matching.
1.54 + ///
1.55 + /// The type of the map that stores the matching arcs.
1.56 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.57 + typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
1.58 +
1.59 + /// \brief Instantiates a MatchingMap.
1.60 + ///
1.61 + /// This function instantiates a \ref MatchingMap.
1.62 + /// \param graph The graph for which we would like to define
1.63 + /// the matching map.
1.64 + static MatchingMap* createMatchingMap(const Graph& graph) {
1.65 + return new MatchingMap(graph);
1.66 + }
1.67 +
1.68 + /// \brief The elevator type used by MaxFractionalMatching algorithm.
1.69 + ///
1.70 + /// The elevator type used by MaxFractionalMatching algorithm.
1.71 + ///
1.72 + /// \sa Elevator
1.73 + /// \sa LinkedElevator
1.74 + typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
1.75 +
1.76 + /// \brief Instantiates an Elevator.
1.77 + ///
1.78 + /// This function instantiates an \ref Elevator.
1.79 + /// \param graph The graph for which we would like to define
1.80 + /// the elevator.
1.81 + /// \param max_level The maximum level of the elevator.
1.82 + static Elevator* createElevator(const Graph& graph, int max_level) {
1.83 + return new Elevator(graph, max_level);
1.84 + }
1.85 + };
1.86 +
1.87 + /// \ingroup matching
1.88 + ///
1.89 + /// \brief Max cardinality fractional matching
1.90 + ///
1.91 + /// This class provides an implementation of fractional matching
1.92 + /// algorithm based on push-relabel principle.
1.93 + ///
1.94 + /// The maximum cardinality fractional matching is a relaxation of the
1.95 + /// maximum cardinality matching problem where the odd set constraints
1.96 + /// are omitted.
1.97 + /// It can be formulated with the following linear program.
1.98 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.99 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.100 + /// \f[\max \sum_{e\in E}x_e\f]
1.101 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.102 + /// \f$X\f$. The result can be represented as the union of a
1.103 + /// matching with one value edges and a set of odd length cycles
1.104 + /// with half value edges.
1.105 + ///
1.106 + /// The algorithm calculates an optimal fractional matching and a
1.107 + /// barrier. The number of adjacents of any node set minus the size
1.108 + /// of node set is a lower bound on the uncovered nodes in the
1.109 + /// graph. For maximum matching a barrier is computed which
1.110 + /// maximizes this difference.
1.111 + ///
1.112 + /// The algorithm can be executed with the run() function. After it
1.113 + /// the matching (the primal solution) and the barrier (the dual
1.114 + /// solution) can be obtained using the query functions.
1.115 + ///
1.116 + /// The primal solution is multiplied by
1.117 + /// \ref MaxWeightedMatching::primalScale "2".
1.118 + ///
1.119 + /// \tparam GR The undirected graph type the algorithm runs on.
1.120 +#ifdef DOXYGEN
1.121 + template <typename GR, typename TR>
1.122 +#else
1.123 + template <typename GR,
1.124 + typename TR = MaxFractionalMatchingDefaultTraits<GR> >
1.125 +#endif
1.126 + class MaxFractionalMatching {
1.127 + public:
1.128 +
1.129 + /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
1.130 + /// class" of the algorithm.
1.131 + typedef TR Traits;
1.132 + /// The type of the graph the algorithm runs on.
1.133 + typedef typename TR::Graph Graph;
1.134 + /// The type of the matching map.
1.135 + typedef typename TR::MatchingMap MatchingMap;
1.136 + /// The type of the elevator.
1.137 + typedef typename TR::Elevator Elevator;
1.138 +
1.139 + /// \brief Scaling factor for primal solution
1.140 + ///
1.141 + /// Scaling factor for primal solution.
1.142 + static const int primalScale = 2;
1.143 +
1.144 + private:
1.145 +
1.146 + const Graph &_graph;
1.147 + int _node_num;
1.148 + bool _allow_loops;
1.149 + int _empty_level;
1.150 +
1.151 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.152 +
1.153 + bool _local_matching;
1.154 + MatchingMap *_matching;
1.155 +
1.156 + bool _local_level;
1.157 + Elevator *_level;
1.158 +
1.159 + typedef typename Graph::template NodeMap<int> InDegMap;
1.160 + InDegMap *_indeg;
1.161 +
1.162 + void createStructures() {
1.163 + _node_num = countNodes(_graph);
1.164 +
1.165 + if (!_matching) {
1.166 + _local_matching = true;
1.167 + _matching = Traits::createMatchingMap(_graph);
1.168 + }
1.169 + if (!_level) {
1.170 + _local_level = true;
1.171 + _level = Traits::createElevator(_graph, _node_num);
1.172 + }
1.173 + if (!_indeg) {
1.174 + _indeg = new InDegMap(_graph);
1.175 + }
1.176 + }
1.177 +
1.178 + void destroyStructures() {
1.179 + if (_local_matching) {
1.180 + delete _matching;
1.181 + }
1.182 + if (_local_level) {
1.183 + delete _level;
1.184 + }
1.185 + if (_indeg) {
1.186 + delete _indeg;
1.187 + }
1.188 + }
1.189 +
1.190 + void postprocessing() {
1.191 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.192 + if ((*_indeg)[n] != 0) continue;
1.193 + _indeg->set(n, -1);
1.194 + Node u = n;
1.195 + while ((*_matching)[u] != INVALID) {
1.196 + Node v = _graph.target((*_matching)[u]);
1.197 + _indeg->set(v, -1);
1.198 + Arc a = _graph.oppositeArc((*_matching)[u]);
1.199 + u = _graph.target((*_matching)[v]);
1.200 + _indeg->set(u, -1);
1.201 + _matching->set(v, a);
1.202 + }
1.203 + }
1.204 +
1.205 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.206 + if ((*_indeg)[n] != 1) continue;
1.207 + _indeg->set(n, -1);
1.208 +
1.209 + int num = 1;
1.210 + Node u = _graph.target((*_matching)[n]);
1.211 + while (u != n) {
1.212 + _indeg->set(u, -1);
1.213 + u = _graph.target((*_matching)[u]);
1.214 + ++num;
1.215 + }
1.216 + if (num % 2 == 0 && num > 2) {
1.217 + Arc prev = _graph.oppositeArc((*_matching)[n]);
1.218 + Node v = _graph.target((*_matching)[n]);
1.219 + u = _graph.target((*_matching)[v]);
1.220 + _matching->set(v, prev);
1.221 + while (u != n) {
1.222 + prev = _graph.oppositeArc((*_matching)[u]);
1.223 + v = _graph.target((*_matching)[u]);
1.224 + u = _graph.target((*_matching)[v]);
1.225 + _matching->set(v, prev);
1.226 + }
1.227 + }
1.228 + }
1.229 + }
1.230 +
1.231 + public:
1.232 +
1.233 + typedef MaxFractionalMatching Create;
1.234 +
1.235 + ///\name Named Template Parameters
1.236 +
1.237 + ///@{
1.238 +
1.239 + template <typename T>
1.240 + struct SetMatchingMapTraits : public Traits {
1.241 + typedef T MatchingMap;
1.242 + static MatchingMap *createMatchingMap(const Graph&) {
1.243 + LEMON_ASSERT(false, "MatchingMap is not initialized");
1.244 + return 0; // ignore warnings
1.245 + }
1.246 + };
1.247 +
1.248 + /// \brief \ref named-templ-param "Named parameter" for setting
1.249 + /// MatchingMap type
1.250 + ///
1.251 + /// \ref named-templ-param "Named parameter" for setting MatchingMap
1.252 + /// type.
1.253 + template <typename T>
1.254 + struct SetMatchingMap
1.255 + : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
1.256 + typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
1.257 + };
1.258 +
1.259 + template <typename T>
1.260 + struct SetElevatorTraits : public Traits {
1.261 + typedef T Elevator;
1.262 + static Elevator *createElevator(const Graph&, int) {
1.263 + LEMON_ASSERT(false, "Elevator is not initialized");
1.264 + return 0; // ignore warnings
1.265 + }
1.266 + };
1.267 +
1.268 + /// \brief \ref named-templ-param "Named parameter" for setting
1.269 + /// Elevator type
1.270 + ///
1.271 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.272 + /// type. If this named parameter is used, then an external
1.273 + /// elevator object must be passed to the algorithm using the
1.274 + /// \ref elevator(Elevator&) "elevator()" function before calling
1.275 + /// \ref run() or \ref init().
1.276 + /// \sa SetStandardElevator
1.277 + template <typename T>
1.278 + struct SetElevator
1.279 + : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
1.280 + typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
1.281 + };
1.282 +
1.283 + template <typename T>
1.284 + struct SetStandardElevatorTraits : public Traits {
1.285 + typedef T Elevator;
1.286 + static Elevator *createElevator(const Graph& graph, int max_level) {
1.287 + return new Elevator(graph, max_level);
1.288 + }
1.289 + };
1.290 +
1.291 + /// \brief \ref named-templ-param "Named parameter" for setting
1.292 + /// Elevator type with automatic allocation
1.293 + ///
1.294 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.295 + /// type with automatic allocation.
1.296 + /// The Elevator should have standard constructor interface to be
1.297 + /// able to automatically created by the algorithm (i.e. the
1.298 + /// graph and the maximum level should be passed to it).
1.299 + /// However an external elevator object could also be passed to the
1.300 + /// algorithm with the \ref elevator(Elevator&) "elevator()" function
1.301 + /// before calling \ref run() or \ref init().
1.302 + /// \sa SetElevator
1.303 + template <typename T>
1.304 + struct SetStandardElevator
1.305 + : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
1.306 + typedef MaxFractionalMatching<Graph,
1.307 + SetStandardElevatorTraits<T> > Create;
1.308 + };
1.309 +
1.310 + /// @}
1.311 +
1.312 + protected:
1.313 +
1.314 + MaxFractionalMatching() {}
1.315 +
1.316 + public:
1.317 +
1.318 + /// \brief Constructor
1.319 + ///
1.320 + /// Constructor.
1.321 + ///
1.322 + MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
1.323 + : _graph(graph), _allow_loops(allow_loops),
1.324 + _local_matching(false), _matching(0),
1.325 + _local_level(false), _level(0), _indeg(0)
1.326 + {}
1.327 +
1.328 + ~MaxFractionalMatching() {
1.329 + destroyStructures();
1.330 + }
1.331 +
1.332 + /// \brief Sets the matching map.
1.333 + ///
1.334 + /// Sets the matching map.
1.335 + /// If you don't use this function before calling \ref run() or
1.336 + /// \ref init(), an instance will be allocated automatically.
1.337 + /// The destructor deallocates this automatically allocated map,
1.338 + /// of course.
1.339 + /// \return <tt>(*this)</tt>
1.340 + MaxFractionalMatching& matchingMap(MatchingMap& map) {
1.341 + if (_local_matching) {
1.342 + delete _matching;
1.343 + _local_matching = false;
1.344 + }
1.345 + _matching = ↦
1.346 + return *this;
1.347 + }
1.348 +
1.349 + /// \brief Sets the elevator used by algorithm.
1.350 + ///
1.351 + /// Sets the elevator used by algorithm.
1.352 + /// If you don't use this function before calling \ref run() or
1.353 + /// \ref init(), an instance will be allocated automatically.
1.354 + /// The destructor deallocates this automatically allocated elevator,
1.355 + /// of course.
1.356 + /// \return <tt>(*this)</tt>
1.357 + MaxFractionalMatching& elevator(Elevator& elevator) {
1.358 + if (_local_level) {
1.359 + delete _level;
1.360 + _local_level = false;
1.361 + }
1.362 + _level = &elevator;
1.363 + return *this;
1.364 + }
1.365 +
1.366 + /// \brief Returns a const reference to the elevator.
1.367 + ///
1.368 + /// Returns a const reference to the elevator.
1.369 + ///
1.370 + /// \pre Either \ref run() or \ref init() must be called before
1.371 + /// using this function.
1.372 + const Elevator& elevator() const {
1.373 + return *_level;
1.374 + }
1.375 +
1.376 + /// \name Execution control
1.377 + /// The simplest way to execute the algorithm is to use one of the
1.378 + /// member functions called \c run(). \n
1.379 + /// If you need more control on the execution, first
1.380 + /// you must call \ref init() and then one variant of the start()
1.381 + /// member.
1.382 +
1.383 + /// @{
1.384 +
1.385 + /// \brief Initializes the internal data structures.
1.386 + ///
1.387 + /// Initializes the internal data structures and sets the initial
1.388 + /// matching.
1.389 + void init() {
1.390 + createStructures();
1.391 +
1.392 + _level->initStart();
1.393 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.394 + _indeg->set(n, 0);
1.395 + _matching->set(n, INVALID);
1.396 + _level->initAddItem(n);
1.397 + }
1.398 + _level->initFinish();
1.399 +
1.400 + _empty_level = _node_num;
1.401 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.402 + for (OutArcIt a(_graph, n); a != INVALID; ++a) {
1.403 + if (_graph.target(a) == n && !_allow_loops) continue;
1.404 + _matching->set(n, a);
1.405 + Node v = _graph.target((*_matching)[n]);
1.406 + _indeg->set(v, (*_indeg)[v] + 1);
1.407 + break;
1.408 + }
1.409 + }
1.410 +
1.411 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.412 + if ((*_indeg)[n] == 0) {
1.413 + _level->activate(n);
1.414 + }
1.415 + }
1.416 + }
1.417 +
1.418 + /// \brief Starts the algorithm and computes a fractional matching
1.419 + ///
1.420 + /// The algorithm computes a maximum fractional matching.
1.421 + ///
1.422 + /// \param postprocess The algorithm computes first a matching
1.423 + /// which is a union of a matching with one value edges, cycles
1.424 + /// with half value edges and even length paths with half value
1.425 + /// edges. If the parameter is true, then after the push-relabel
1.426 + /// algorithm it postprocesses the matching to contain only
1.427 + /// matching edges and half value odd cycles.
1.428 + void start(bool postprocess = true) {
1.429 + Node n;
1.430 + while ((n = _level->highestActive()) != INVALID) {
1.431 + int level = _level->highestActiveLevel();
1.432 + int new_level = _level->maxLevel();
1.433 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
1.434 + Node u = _graph.source(a);
1.435 + if (n == u && !_allow_loops) continue;
1.436 + Node v = _graph.target((*_matching)[u]);
1.437 + if ((*_level)[v] < level) {
1.438 + _indeg->set(v, (*_indeg)[v] - 1);
1.439 + if ((*_indeg)[v] == 0) {
1.440 + _level->activate(v);
1.441 + }
1.442 + _matching->set(u, a);
1.443 + _indeg->set(n, (*_indeg)[n] + 1);
1.444 + _level->deactivate(n);
1.445 + goto no_more_push;
1.446 + } else if (new_level > (*_level)[v]) {
1.447 + new_level = (*_level)[v];
1.448 + }
1.449 + }
1.450 +
1.451 + if (new_level + 1 < _level->maxLevel()) {
1.452 + _level->liftHighestActive(new_level + 1);
1.453 + } else {
1.454 + _level->liftHighestActiveToTop();
1.455 + }
1.456 + if (_level->emptyLevel(level)) {
1.457 + _level->liftToTop(level);
1.458 + }
1.459 + no_more_push:
1.460 + ;
1.461 + }
1.462 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.463 + if ((*_matching)[n] == INVALID) continue;
1.464 + Node u = _graph.target((*_matching)[n]);
1.465 + if ((*_indeg)[u] > 1) {
1.466 + _indeg->set(u, (*_indeg)[u] - 1);
1.467 + _matching->set(n, INVALID);
1.468 + }
1.469 + }
1.470 + if (postprocess) {
1.471 + postprocessing();
1.472 + }
1.473 + }
1.474 +
1.475 + /// \brief Starts the algorithm and computes a perfect fractional
1.476 + /// matching
1.477 + ///
1.478 + /// The algorithm computes a perfect fractional matching. If it
1.479 + /// does not exists, then the algorithm returns false and the
1.480 + /// matching is undefined and the barrier.
1.481 + ///
1.482 + /// \param postprocess The algorithm computes first a matching
1.483 + /// which is a union of a matching with one value edges, cycles
1.484 + /// with half value edges and even length paths with half value
1.485 + /// edges. If the parameter is true, then after the push-relabel
1.486 + /// algorithm it postprocesses the matching to contain only
1.487 + /// matching edges and half value odd cycles.
1.488 + bool startPerfect(bool postprocess = true) {
1.489 + Node n;
1.490 + while ((n = _level->highestActive()) != INVALID) {
1.491 + int level = _level->highestActiveLevel();
1.492 + int new_level = _level->maxLevel();
1.493 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
1.494 + Node u = _graph.source(a);
1.495 + if (n == u && !_allow_loops) continue;
1.496 + Node v = _graph.target((*_matching)[u]);
1.497 + if ((*_level)[v] < level) {
1.498 + _indeg->set(v, (*_indeg)[v] - 1);
1.499 + if ((*_indeg)[v] == 0) {
1.500 + _level->activate(v);
1.501 + }
1.502 + _matching->set(u, a);
1.503 + _indeg->set(n, (*_indeg)[n] + 1);
1.504 + _level->deactivate(n);
1.505 + goto no_more_push;
1.506 + } else if (new_level > (*_level)[v]) {
1.507 + new_level = (*_level)[v];
1.508 + }
1.509 + }
1.510 +
1.511 + if (new_level + 1 < _level->maxLevel()) {
1.512 + _level->liftHighestActive(new_level + 1);
1.513 + } else {
1.514 + _level->liftHighestActiveToTop();
1.515 + _empty_level = _level->maxLevel() - 1;
1.516 + return false;
1.517 + }
1.518 + if (_level->emptyLevel(level)) {
1.519 + _level->liftToTop(level);
1.520 + _empty_level = level;
1.521 + return false;
1.522 + }
1.523 + no_more_push:
1.524 + ;
1.525 + }
1.526 + if (postprocess) {
1.527 + postprocessing();
1.528 + }
1.529 + return true;
1.530 + }
1.531 +
1.532 + /// \brief Runs the algorithm
1.533 + ///
1.534 + /// Just a shortcut for the next code:
1.535 + ///\code
1.536 + /// init();
1.537 + /// start();
1.538 + ///\endcode
1.539 + void run(bool postprocess = true) {
1.540 + init();
1.541 + start(postprocess);
1.542 + }
1.543 +
1.544 + /// \brief Runs the algorithm to find a perfect fractional matching
1.545 + ///
1.546 + /// Just a shortcut for the next code:
1.547 + ///\code
1.548 + /// init();
1.549 + /// startPerfect();
1.550 + ///\endcode
1.551 + bool runPerfect(bool postprocess = true) {
1.552 + init();
1.553 + return startPerfect(postprocess);
1.554 + }
1.555 +
1.556 + ///@}
1.557 +
1.558 + /// \name Query Functions
1.559 + /// The result of the %Matching algorithm can be obtained using these
1.560 + /// functions.\n
1.561 + /// Before the use of these functions,
1.562 + /// either run() or start() must be called.
1.563 + ///@{
1.564 +
1.565 +
1.566 + /// \brief Return the number of covered nodes in the matching.
1.567 + ///
1.568 + /// This function returns the number of covered nodes in the matching.
1.569 + ///
1.570 + /// \pre Either run() or start() must be called before using this function.
1.571 + int matchingSize() const {
1.572 + int num = 0;
1.573 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.574 + if ((*_matching)[n] != INVALID) {
1.575 + ++num;
1.576 + }
1.577 + }
1.578 + return num;
1.579 + }
1.580 +
1.581 + /// \brief Returns a const reference to the matching map.
1.582 + ///
1.583 + /// Returns a const reference to the node map storing the found
1.584 + /// fractional matching. This method can be called after
1.585 + /// running the algorithm.
1.586 + ///
1.587 + /// \pre Either \ref run() or \ref init() must be called before
1.588 + /// using this function.
1.589 + const MatchingMap& matchingMap() const {
1.590 + return *_matching;
1.591 + }
1.592 +
1.593 + /// \brief Return \c true if the given edge is in the matching.
1.594 + ///
1.595 + /// This function returns \c true if the given edge is in the
1.596 + /// found matching. The result is scaled by \ref primalScale
1.597 + /// "primal scale".
1.598 + ///
1.599 + /// \pre Either run() or start() must be called before using this function.
1.600 + int matching(const Edge& edge) const {
1.601 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
1.602 + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
1.603 + }
1.604 +
1.605 + /// \brief Return the fractional matching arc (or edge) incident
1.606 + /// to the given node.
1.607 + ///
1.608 + /// This function returns one of the fractional matching arc (or
1.609 + /// edge) incident to the given node in the found matching or \c
1.610 + /// INVALID if the node is not covered by the matching or if the
1.611 + /// node is on an odd length cycle then it is the successor edge
1.612 + /// on the cycle.
1.613 + ///
1.614 + /// \pre Either run() or start() must be called before using this function.
1.615 + Arc matching(const Node& node) const {
1.616 + return (*_matching)[node];
1.617 + }
1.618 +
1.619 + /// \brief Returns true if the node is in the barrier
1.620 + ///
1.621 + /// The barrier is a subset of the nodes. If the nodes in the
1.622 + /// barrier have less adjacent nodes than the size of the barrier,
1.623 + /// then at least as much nodes cannot be covered as the
1.624 + /// difference of the two subsets.
1.625 + bool barrier(const Node& node) const {
1.626 + return (*_level)[node] >= _empty_level;
1.627 + }
1.628 +
1.629 + /// @}
1.630 +
1.631 + };
1.632 +
1.633 + /// \ingroup matching
1.634 + ///
1.635 + /// \brief Weighted fractional matching in general graphs
1.636 + ///
1.637 + /// This class provides an efficient implementation of fractional
1.638 + /// matching algorithm. The implementation is based on extensive use
1.639 + /// of priority queues and provides \f$O(nm\log n)\f$ time
1.640 + /// complexity.
1.641 + ///
1.642 + /// The maximum weighted fractional matching is a relaxation of the
1.643 + /// maximum weighted matching problem where the odd set constraints
1.644 + /// are omitted.
1.645 + /// It can be formulated with the following linear program.
1.646 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.647 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.648 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.649 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.650 + /// \f$X\f$. The result must be the union of a matching with one
1.651 + /// value edges and a set of odd length cycles with half value edges.
1.652 + ///
1.653 + /// The algorithm calculates an optimal fractional matching and a
1.654 + /// proof of the optimality. The solution of the dual problem can be
1.655 + /// used to check the result of the algorithm. The dual linear
1.656 + /// problem is the following.
1.657 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1.658 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.659 + /// \f[\min \sum_{u \in V}y_u \f] ///
1.660 + ///
1.661 + /// The algorithm can be executed with the run() function.
1.662 + /// After it the matching (the primal solution) and the dual solution
1.663 + /// can be obtained using the query functions.
1.664 + ///
1.665 + /// If the value type is integer, then the primal and the dual
1.666 + /// solutions are multiplied by
1.667 + /// \ref MaxWeightedMatching::primalScale "2" and
1.668 + /// \ref MaxWeightedMatching::dualScale "4" respectively.
1.669 + ///
1.670 + /// \tparam GR The undirected graph type the algorithm runs on.
1.671 + /// \tparam WM The type edge weight map. The default type is
1.672 + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.673 +#ifdef DOXYGEN
1.674 + template <typename GR, typename WM>
1.675 +#else
1.676 + template <typename GR,
1.677 + typename WM = typename GR::template EdgeMap<int> >
1.678 +#endif
1.679 + class MaxWeightedFractionalMatching {
1.680 + public:
1.681 +
1.682 + /// The graph type of the algorithm
1.683 + typedef GR Graph;
1.684 + /// The type of the edge weight map
1.685 + typedef WM WeightMap;
1.686 + /// The value type of the edge weights
1.687 + typedef typename WeightMap::Value Value;
1.688 +
1.689 + /// The type of the matching map
1.690 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.691 + MatchingMap;
1.692 +
1.693 + /// \brief Scaling factor for primal solution
1.694 + ///
1.695 + /// Scaling factor for primal solution. It is equal to 2 or 1
1.696 + /// according to the value type.
1.697 + static const int primalScale =
1.698 + std::numeric_limits<Value>::is_integer ? 2 : 1;
1.699 +
1.700 + /// \brief Scaling factor for dual solution
1.701 + ///
1.702 + /// Scaling factor for dual solution. It is equal to 4 or 1
1.703 + /// according to the value type.
1.704 + static const int dualScale =
1.705 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.706 +
1.707 + private:
1.708 +
1.709 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.710 +
1.711 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.712 +
1.713 + const Graph& _graph;
1.714 + const WeightMap& _weight;
1.715 +
1.716 + MatchingMap* _matching;
1.717 + NodePotential* _node_potential;
1.718 +
1.719 + int _node_num;
1.720 + bool _allow_loops;
1.721 +
1.722 + enum Status {
1.723 + EVEN = -1, MATCHED = 0, ODD = 1
1.724 + };
1.725 +
1.726 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.727 + StatusMap* _status;
1.728 +
1.729 + typedef typename Graph::template NodeMap<Arc> PredMap;
1.730 + PredMap* _pred;
1.731 +
1.732 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.733 +
1.734 + IntNodeMap *_tree_set_index;
1.735 + TreeSet *_tree_set;
1.736 +
1.737 + IntNodeMap *_delta1_index;
1.738 + BinHeap<Value, IntNodeMap> *_delta1;
1.739 +
1.740 + IntNodeMap *_delta2_index;
1.741 + BinHeap<Value, IntNodeMap> *_delta2;
1.742 +
1.743 + IntEdgeMap *_delta3_index;
1.744 + BinHeap<Value, IntEdgeMap> *_delta3;
1.745 +
1.746 + Value _delta_sum;
1.747 +
1.748 + void createStructures() {
1.749 + _node_num = countNodes(_graph);
1.750 +
1.751 + if (!_matching) {
1.752 + _matching = new MatchingMap(_graph);
1.753 + }
1.754 + if (!_node_potential) {
1.755 + _node_potential = new NodePotential(_graph);
1.756 + }
1.757 + if (!_status) {
1.758 + _status = new StatusMap(_graph);
1.759 + }
1.760 + if (!_pred) {
1.761 + _pred = new PredMap(_graph);
1.762 + }
1.763 + if (!_tree_set) {
1.764 + _tree_set_index = new IntNodeMap(_graph);
1.765 + _tree_set = new TreeSet(*_tree_set_index);
1.766 + }
1.767 + if (!_delta1) {
1.768 + _delta1_index = new IntNodeMap(_graph);
1.769 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
1.770 + }
1.771 + if (!_delta2) {
1.772 + _delta2_index = new IntNodeMap(_graph);
1.773 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1.774 + }
1.775 + if (!_delta3) {
1.776 + _delta3_index = new IntEdgeMap(_graph);
1.777 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.778 + }
1.779 + }
1.780 +
1.781 + void destroyStructures() {
1.782 + if (_matching) {
1.783 + delete _matching;
1.784 + }
1.785 + if (_node_potential) {
1.786 + delete _node_potential;
1.787 + }
1.788 + if (_status) {
1.789 + delete _status;
1.790 + }
1.791 + if (_pred) {
1.792 + delete _pred;
1.793 + }
1.794 + if (_tree_set) {
1.795 + delete _tree_set_index;
1.796 + delete _tree_set;
1.797 + }
1.798 + if (_delta1) {
1.799 + delete _delta1_index;
1.800 + delete _delta1;
1.801 + }
1.802 + if (_delta2) {
1.803 + delete _delta2_index;
1.804 + delete _delta2;
1.805 + }
1.806 + if (_delta3) {
1.807 + delete _delta3_index;
1.808 + delete _delta3;
1.809 + }
1.810 + }
1.811 +
1.812 + void matchedToEven(Node node, int tree) {
1.813 + _tree_set->insert(node, tree);
1.814 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.815 + _delta1->push(node, (*_node_potential)[node]);
1.816 +
1.817 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.818 + _delta2->erase(node);
1.819 + }
1.820 +
1.821 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.822 + Node v = _graph.source(a);
1.823 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.824 + dualScale * _weight[a];
1.825 + if (node == v) {
1.826 + if (_allow_loops && _graph.direction(a)) {
1.827 + _delta3->push(a, rw / 2);
1.828 + }
1.829 + } else if ((*_status)[v] == EVEN) {
1.830 + _delta3->push(a, rw / 2);
1.831 + } else if ((*_status)[v] == MATCHED) {
1.832 + if (_delta2->state(v) != _delta2->IN_HEAP) {
1.833 + _pred->set(v, a);
1.834 + _delta2->push(v, rw);
1.835 + } else if ((*_delta2)[v] > rw) {
1.836 + _pred->set(v, a);
1.837 + _delta2->decrease(v, rw);
1.838 + }
1.839 + }
1.840 + }
1.841 + }
1.842 +
1.843 + void matchedToOdd(Node node, int tree) {
1.844 + _tree_set->insert(node, tree);
1.845 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.846 +
1.847 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.848 + _delta2->erase(node);
1.849 + }
1.850 + }
1.851 +
1.852 + void evenToMatched(Node node, int tree) {
1.853 + _delta1->erase(node);
1.854 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.855 + Arc min = INVALID;
1.856 + Value minrw = std::numeric_limits<Value>::max();
1.857 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.858 + Node v = _graph.source(a);
1.859 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.860 + dualScale * _weight[a];
1.861 +
1.862 + if (node == v) {
1.863 + if (_allow_loops && _graph.direction(a)) {
1.864 + _delta3->erase(a);
1.865 + }
1.866 + } else if ((*_status)[v] == EVEN) {
1.867 + _delta3->erase(a);
1.868 + if (minrw > rw) {
1.869 + min = _graph.oppositeArc(a);
1.870 + minrw = rw;
1.871 + }
1.872 + } else if ((*_status)[v] == MATCHED) {
1.873 + if ((*_pred)[v] == a) {
1.874 + Arc mina = INVALID;
1.875 + Value minrwa = std::numeric_limits<Value>::max();
1.876 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1.877 + Node va = _graph.target(aa);
1.878 + if ((*_status)[va] != EVEN ||
1.879 + _tree_set->find(va) == tree) continue;
1.880 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1.881 + dualScale * _weight[aa];
1.882 + if (minrwa > rwa) {
1.883 + minrwa = rwa;
1.884 + mina = aa;
1.885 + }
1.886 + }
1.887 + if (mina != INVALID) {
1.888 + _pred->set(v, mina);
1.889 + _delta2->increase(v, minrwa);
1.890 + } else {
1.891 + _pred->set(v, INVALID);
1.892 + _delta2->erase(v);
1.893 + }
1.894 + }
1.895 + }
1.896 + }
1.897 + if (min != INVALID) {
1.898 + _pred->set(node, min);
1.899 + _delta2->push(node, minrw);
1.900 + } else {
1.901 + _pred->set(node, INVALID);
1.902 + }
1.903 + }
1.904 +
1.905 + void oddToMatched(Node node) {
1.906 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.907 + Arc min = INVALID;
1.908 + Value minrw = std::numeric_limits<Value>::max();
1.909 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.910 + Node v = _graph.source(a);
1.911 + if ((*_status)[v] != EVEN) continue;
1.912 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.913 + dualScale * _weight[a];
1.914 +
1.915 + if (minrw > rw) {
1.916 + min = _graph.oppositeArc(a);
1.917 + minrw = rw;
1.918 + }
1.919 + }
1.920 + if (min != INVALID) {
1.921 + _pred->set(node, min);
1.922 + _delta2->push(node, minrw);
1.923 + } else {
1.924 + _pred->set(node, INVALID);
1.925 + }
1.926 + }
1.927 +
1.928 + void alternatePath(Node even, int tree) {
1.929 + Node odd;
1.930 +
1.931 + _status->set(even, MATCHED);
1.932 + evenToMatched(even, tree);
1.933 +
1.934 + Arc prev = (*_matching)[even];
1.935 + while (prev != INVALID) {
1.936 + odd = _graph.target(prev);
1.937 + even = _graph.target((*_pred)[odd]);
1.938 + _matching->set(odd, (*_pred)[odd]);
1.939 + _status->set(odd, MATCHED);
1.940 + oddToMatched(odd);
1.941 +
1.942 + prev = (*_matching)[even];
1.943 + _status->set(even, MATCHED);
1.944 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.945 + evenToMatched(even, tree);
1.946 + }
1.947 + }
1.948 +
1.949 + void destroyTree(int tree) {
1.950 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1.951 + if ((*_status)[n] == EVEN) {
1.952 + _status->set(n, MATCHED);
1.953 + evenToMatched(n, tree);
1.954 + } else if ((*_status)[n] == ODD) {
1.955 + _status->set(n, MATCHED);
1.956 + oddToMatched(n);
1.957 + }
1.958 + }
1.959 + _tree_set->eraseClass(tree);
1.960 + }
1.961 +
1.962 +
1.963 + void unmatchNode(const Node& node) {
1.964 + int tree = _tree_set->find(node);
1.965 +
1.966 + alternatePath(node, tree);
1.967 + destroyTree(tree);
1.968 +
1.969 + _matching->set(node, INVALID);
1.970 + }
1.971 +
1.972 +
1.973 + void augmentOnEdge(const Edge& edge) {
1.974 + Node left = _graph.u(edge);
1.975 + int left_tree = _tree_set->find(left);
1.976 +
1.977 + alternatePath(left, left_tree);
1.978 + destroyTree(left_tree);
1.979 + _matching->set(left, _graph.direct(edge, true));
1.980 +
1.981 + Node right = _graph.v(edge);
1.982 + int right_tree = _tree_set->find(right);
1.983 +
1.984 + alternatePath(right, right_tree);
1.985 + destroyTree(right_tree);
1.986 + _matching->set(right, _graph.direct(edge, false));
1.987 + }
1.988 +
1.989 + void augmentOnArc(const Arc& arc) {
1.990 + Node left = _graph.source(arc);
1.991 + _status->set(left, MATCHED);
1.992 + _matching->set(left, arc);
1.993 + _pred->set(left, arc);
1.994 +
1.995 + Node right = _graph.target(arc);
1.996 + int right_tree = _tree_set->find(right);
1.997 +
1.998 + alternatePath(right, right_tree);
1.999 + destroyTree(right_tree);
1.1000 + _matching->set(right, _graph.oppositeArc(arc));
1.1001 + }
1.1002 +
1.1003 + void extendOnArc(const Arc& arc) {
1.1004 + Node base = _graph.target(arc);
1.1005 + int tree = _tree_set->find(base);
1.1006 +
1.1007 + Node odd = _graph.source(arc);
1.1008 + _tree_set->insert(odd, tree);
1.1009 + _status->set(odd, ODD);
1.1010 + matchedToOdd(odd, tree);
1.1011 + _pred->set(odd, arc);
1.1012 +
1.1013 + Node even = _graph.target((*_matching)[odd]);
1.1014 + _tree_set->insert(even, tree);
1.1015 + _status->set(even, EVEN);
1.1016 + matchedToEven(even, tree);
1.1017 + }
1.1018 +
1.1019 + void cycleOnEdge(const Edge& edge, int tree) {
1.1020 + Node nca = INVALID;
1.1021 + std::vector<Node> left_path, right_path;
1.1022 +
1.1023 + {
1.1024 + std::set<Node> left_set, right_set;
1.1025 + Node left = _graph.u(edge);
1.1026 + left_path.push_back(left);
1.1027 + left_set.insert(left);
1.1028 +
1.1029 + Node right = _graph.v(edge);
1.1030 + right_path.push_back(right);
1.1031 + right_set.insert(right);
1.1032 +
1.1033 + while (true) {
1.1034 +
1.1035 + if (left_set.find(right) != left_set.end()) {
1.1036 + nca = right;
1.1037 + break;
1.1038 + }
1.1039 +
1.1040 + if ((*_matching)[left] == INVALID) break;
1.1041 +
1.1042 + left = _graph.target((*_matching)[left]);
1.1043 + left_path.push_back(left);
1.1044 + left = _graph.target((*_pred)[left]);
1.1045 + left_path.push_back(left);
1.1046 +
1.1047 + left_set.insert(left);
1.1048 +
1.1049 + if (right_set.find(left) != right_set.end()) {
1.1050 + nca = left;
1.1051 + break;
1.1052 + }
1.1053 +
1.1054 + if ((*_matching)[right] == INVALID) break;
1.1055 +
1.1056 + right = _graph.target((*_matching)[right]);
1.1057 + right_path.push_back(right);
1.1058 + right = _graph.target((*_pred)[right]);
1.1059 + right_path.push_back(right);
1.1060 +
1.1061 + right_set.insert(right);
1.1062 +
1.1063 + }
1.1064 +
1.1065 + if (nca == INVALID) {
1.1066 + if ((*_matching)[left] == INVALID) {
1.1067 + nca = right;
1.1068 + while (left_set.find(nca) == left_set.end()) {
1.1069 + nca = _graph.target((*_matching)[nca]);
1.1070 + right_path.push_back(nca);
1.1071 + nca = _graph.target((*_pred)[nca]);
1.1072 + right_path.push_back(nca);
1.1073 + }
1.1074 + } else {
1.1075 + nca = left;
1.1076 + while (right_set.find(nca) == right_set.end()) {
1.1077 + nca = _graph.target((*_matching)[nca]);
1.1078 + left_path.push_back(nca);
1.1079 + nca = _graph.target((*_pred)[nca]);
1.1080 + left_path.push_back(nca);
1.1081 + }
1.1082 + }
1.1083 + }
1.1084 + }
1.1085 +
1.1086 + alternatePath(nca, tree);
1.1087 + Arc prev;
1.1088 +
1.1089 + prev = _graph.direct(edge, true);
1.1090 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1091 + _matching->set(left_path[i], prev);
1.1092 + _status->set(left_path[i], MATCHED);
1.1093 + evenToMatched(left_path[i], tree);
1.1094 +
1.1095 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1.1096 + _status->set(left_path[i + 1], MATCHED);
1.1097 + oddToMatched(left_path[i + 1]);
1.1098 + }
1.1099 + _matching->set(nca, prev);
1.1100 +
1.1101 + for (int i = 0; right_path[i] != nca; i += 2) {
1.1102 + _status->set(right_path[i], MATCHED);
1.1103 + evenToMatched(right_path[i], tree);
1.1104 +
1.1105 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1.1106 + _status->set(right_path[i + 1], MATCHED);
1.1107 + oddToMatched(right_path[i + 1]);
1.1108 + }
1.1109 +
1.1110 + destroyTree(tree);
1.1111 + }
1.1112 +
1.1113 + void extractCycle(const Arc &arc) {
1.1114 + Node left = _graph.source(arc);
1.1115 + Node odd = _graph.target((*_matching)[left]);
1.1116 + Arc prev;
1.1117 + while (odd != left) {
1.1118 + Node even = _graph.target((*_matching)[odd]);
1.1119 + prev = (*_matching)[odd];
1.1120 + odd = _graph.target((*_matching)[even]);
1.1121 + _matching->set(even, _graph.oppositeArc(prev));
1.1122 + }
1.1123 + _matching->set(left, arc);
1.1124 +
1.1125 + Node right = _graph.target(arc);
1.1126 + int right_tree = _tree_set->find(right);
1.1127 + alternatePath(right, right_tree);
1.1128 + destroyTree(right_tree);
1.1129 + _matching->set(right, _graph.oppositeArc(arc));
1.1130 + }
1.1131 +
1.1132 + public:
1.1133 +
1.1134 + /// \brief Constructor
1.1135 + ///
1.1136 + /// Constructor.
1.1137 + MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
1.1138 + bool allow_loops = true)
1.1139 + : _graph(graph), _weight(weight), _matching(0),
1.1140 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1.1141 + _status(0), _pred(0),
1.1142 + _tree_set_index(0), _tree_set(0),
1.1143 +
1.1144 + _delta1_index(0), _delta1(0),
1.1145 + _delta2_index(0), _delta2(0),
1.1146 + _delta3_index(0), _delta3(0),
1.1147 +
1.1148 + _delta_sum() {}
1.1149 +
1.1150 + ~MaxWeightedFractionalMatching() {
1.1151 + destroyStructures();
1.1152 + }
1.1153 +
1.1154 + /// \name Execution Control
1.1155 + /// The simplest way to execute the algorithm is to use the
1.1156 + /// \ref run() member function.
1.1157 +
1.1158 + ///@{
1.1159 +
1.1160 + /// \brief Initialize the algorithm
1.1161 + ///
1.1162 + /// This function initializes the algorithm.
1.1163 + void init() {
1.1164 + createStructures();
1.1165 +
1.1166 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1167 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
1.1168 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
1.1169 + }
1.1170 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1171 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.1172 + }
1.1173 +
1.1174 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1175 + Value max = 0;
1.1176 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1177 + if (_graph.target(e) == n && !_allow_loops) continue;
1.1178 + if ((dualScale * _weight[e]) / 2 > max) {
1.1179 + max = (dualScale * _weight[e]) / 2;
1.1180 + }
1.1181 + }
1.1182 + _node_potential->set(n, max);
1.1183 + _delta1->push(n, max);
1.1184 +
1.1185 + _tree_set->insert(n);
1.1186 +
1.1187 + _matching->set(n, INVALID);
1.1188 + _status->set(n, EVEN);
1.1189 + }
1.1190 +
1.1191 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1192 + Node left = _graph.u(e);
1.1193 + Node right = _graph.v(e);
1.1194 + if (left == right && !_allow_loops) continue;
1.1195 + _delta3->push(e, ((*_node_potential)[left] +
1.1196 + (*_node_potential)[right] -
1.1197 + dualScale * _weight[e]) / 2);
1.1198 + }
1.1199 + }
1.1200 +
1.1201 + /// \brief Start the algorithm
1.1202 + ///
1.1203 + /// This function starts the algorithm.
1.1204 + ///
1.1205 + /// \pre \ref init() must be called before using this function.
1.1206 + void start() {
1.1207 + enum OpType {
1.1208 + D1, D2, D3
1.1209 + };
1.1210 +
1.1211 + int unmatched = _node_num;
1.1212 + while (unmatched > 0) {
1.1213 + Value d1 = !_delta1->empty() ?
1.1214 + _delta1->prio() : std::numeric_limits<Value>::max();
1.1215 +
1.1216 + Value d2 = !_delta2->empty() ?
1.1217 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1218 +
1.1219 + Value d3 = !_delta3->empty() ?
1.1220 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1221 +
1.1222 + _delta_sum = d3; OpType ot = D3;
1.1223 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1.1224 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1225 +
1.1226 + switch (ot) {
1.1227 + case D1:
1.1228 + {
1.1229 + Node n = _delta1->top();
1.1230 + unmatchNode(n);
1.1231 + --unmatched;
1.1232 + }
1.1233 + break;
1.1234 + case D2:
1.1235 + {
1.1236 + Node n = _delta2->top();
1.1237 + Arc a = (*_pred)[n];
1.1238 + if ((*_matching)[n] == INVALID) {
1.1239 + augmentOnArc(a);
1.1240 + --unmatched;
1.1241 + } else {
1.1242 + Node v = _graph.target((*_matching)[n]);
1.1243 + if ((*_matching)[n] !=
1.1244 + _graph.oppositeArc((*_matching)[v])) {
1.1245 + extractCycle(a);
1.1246 + --unmatched;
1.1247 + } else {
1.1248 + extendOnArc(a);
1.1249 + }
1.1250 + }
1.1251 + } break;
1.1252 + case D3:
1.1253 + {
1.1254 + Edge e = _delta3->top();
1.1255 +
1.1256 + Node left = _graph.u(e);
1.1257 + Node right = _graph.v(e);
1.1258 +
1.1259 + int left_tree = _tree_set->find(left);
1.1260 + int right_tree = _tree_set->find(right);
1.1261 +
1.1262 + if (left_tree == right_tree) {
1.1263 + cycleOnEdge(e, left_tree);
1.1264 + --unmatched;
1.1265 + } else {
1.1266 + augmentOnEdge(e);
1.1267 + unmatched -= 2;
1.1268 + }
1.1269 + } break;
1.1270 + }
1.1271 + }
1.1272 + }
1.1273 +
1.1274 + /// \brief Run the algorithm.
1.1275 + ///
1.1276 + /// This method runs the \c %MaxWeightedMatching algorithm.
1.1277 + ///
1.1278 + /// \note mwfm.run() is just a shortcut of the following code.
1.1279 + /// \code
1.1280 + /// mwfm.init();
1.1281 + /// mwfm.start();
1.1282 + /// \endcode
1.1283 + void run() {
1.1284 + init();
1.1285 + start();
1.1286 + }
1.1287 +
1.1288 + /// @}
1.1289 +
1.1290 + /// \name Primal Solution
1.1291 + /// Functions to get the primal solution, i.e. the maximum weighted
1.1292 + /// matching.\n
1.1293 + /// Either \ref run() or \ref start() function should be called before
1.1294 + /// using them.
1.1295 +
1.1296 + /// @{
1.1297 +
1.1298 + /// \brief Return the weight of the matching.
1.1299 + ///
1.1300 + /// This function returns the weight of the found matching. This
1.1301 + /// value is scaled by \ref primalScale "primal scale".
1.1302 + ///
1.1303 + /// \pre Either run() or start() must be called before using this function.
1.1304 + Value matchingWeight() const {
1.1305 + Value sum = 0;
1.1306 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1307 + if ((*_matching)[n] != INVALID) {
1.1308 + sum += _weight[(*_matching)[n]];
1.1309 + }
1.1310 + }
1.1311 + return sum * primalScale / 2;
1.1312 + }
1.1313 +
1.1314 + /// \brief Return the number of covered nodes in the matching.
1.1315 + ///
1.1316 + /// This function returns the number of covered nodes in the matching.
1.1317 + ///
1.1318 + /// \pre Either run() or start() must be called before using this function.
1.1319 + int matchingSize() const {
1.1320 + int num = 0;
1.1321 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1322 + if ((*_matching)[n] != INVALID) {
1.1323 + ++num;
1.1324 + }
1.1325 + }
1.1326 + return num;
1.1327 + }
1.1328 +
1.1329 + /// \brief Return \c true if the given edge is in the matching.
1.1330 + ///
1.1331 + /// This function returns \c true if the given edge is in the
1.1332 + /// found matching. The result is scaled by \ref primalScale
1.1333 + /// "primal scale".
1.1334 + ///
1.1335 + /// \pre Either run() or start() must be called before using this function.
1.1336 + Value matching(const Edge& edge) const {
1.1337 + return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1.1338 + * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
1.1339 + * primalScale / 2;
1.1340 + }
1.1341 +
1.1342 + /// \brief Return the fractional matching arc (or edge) incident
1.1343 + /// to the given node.
1.1344 + ///
1.1345 + /// This function returns one of the fractional matching arc (or
1.1346 + /// edge) incident to the given node in the found matching or \c
1.1347 + /// INVALID if the node is not covered by the matching or if the
1.1348 + /// node is on an odd length cycle then it is the successor edge
1.1349 + /// on the cycle.
1.1350 + ///
1.1351 + /// \pre Either run() or start() must be called before using this function.
1.1352 + Arc matching(const Node& node) const {
1.1353 + return (*_matching)[node];
1.1354 + }
1.1355 +
1.1356 + /// \brief Return a const reference to the matching map.
1.1357 + ///
1.1358 + /// This function returns a const reference to a node map that stores
1.1359 + /// the matching arc (or edge) incident to each node.
1.1360 + const MatchingMap& matchingMap() const {
1.1361 + return *_matching;
1.1362 + }
1.1363 +
1.1364 + /// @}
1.1365 +
1.1366 + /// \name Dual Solution
1.1367 + /// Functions to get the dual solution.\n
1.1368 + /// Either \ref run() or \ref start() function should be called before
1.1369 + /// using them.
1.1370 +
1.1371 + /// @{
1.1372 +
1.1373 + /// \brief Return the value of the dual solution.
1.1374 + ///
1.1375 + /// This function returns the value of the dual solution.
1.1376 + /// It should be equal to the primal value scaled by \ref dualScale
1.1377 + /// "dual scale".
1.1378 + ///
1.1379 + /// \pre Either run() or start() must be called before using this function.
1.1380 + Value dualValue() const {
1.1381 + Value sum = 0;
1.1382 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1383 + sum += nodeValue(n);
1.1384 + }
1.1385 + return sum;
1.1386 + }
1.1387 +
1.1388 + /// \brief Return the dual value (potential) of the given node.
1.1389 + ///
1.1390 + /// This function returns the dual value (potential) of the given node.
1.1391 + ///
1.1392 + /// \pre Either run() or start() must be called before using this function.
1.1393 + Value nodeValue(const Node& n) const {
1.1394 + return (*_node_potential)[n];
1.1395 + }
1.1396 +
1.1397 + /// @}
1.1398 +
1.1399 + };
1.1400 +
1.1401 + /// \ingroup matching
1.1402 + ///
1.1403 + /// \brief Weighted fractional perfect matching in general graphs
1.1404 + ///
1.1405 + /// This class provides an efficient implementation of fractional
1.1406 + /// matching algorithm. The implementation is based on extensive use
1.1407 + /// of priority queues and provides \f$O(nm\log n)\f$ time
1.1408 + /// complexity.
1.1409 + ///
1.1410 + /// The maximum weighted fractional perfect matching is a relaxation
1.1411 + /// of the maximum weighted perfect matching problem where the odd
1.1412 + /// set constraints are omitted.
1.1413 + /// It can be formulated with the following linear program.
1.1414 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1415 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1416 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1417 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1418 + /// \f$X\f$. The result must be the union of a matching with one
1.1419 + /// value edges and a set of odd length cycles with half value edges.
1.1420 + ///
1.1421 + /// The algorithm calculates an optimal fractional matching and a
1.1422 + /// proof of the optimality. The solution of the dual problem can be
1.1423 + /// used to check the result of the algorithm. The dual linear
1.1424 + /// problem is the following.
1.1425 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1.1426 + /// \f[\min \sum_{u \in V}y_u \f] ///
1.1427 + ///
1.1428 + /// The algorithm can be executed with the run() function.
1.1429 + /// After it the matching (the primal solution) and the dual solution
1.1430 + /// can be obtained using the query functions.
1.1431 +
1.1432 + /// If the value type is integer, then the primal and the dual
1.1433 + /// solutions are multiplied by
1.1434 + /// \ref MaxWeightedMatching::primalScale "2" and
1.1435 + /// \ref MaxWeightedMatching::dualScale "4" respectively.
1.1436 + ///
1.1437 + /// \tparam GR The undirected graph type the algorithm runs on.
1.1438 + /// \tparam WM The type edge weight map. The default type is
1.1439 + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.1440 +#ifdef DOXYGEN
1.1441 + template <typename GR, typename WM>
1.1442 +#else
1.1443 + template <typename GR,
1.1444 + typename WM = typename GR::template EdgeMap<int> >
1.1445 +#endif
1.1446 + class MaxWeightedPerfectFractionalMatching {
1.1447 + public:
1.1448 +
1.1449 + /// The graph type of the algorithm
1.1450 + typedef GR Graph;
1.1451 + /// The type of the edge weight map
1.1452 + typedef WM WeightMap;
1.1453 + /// The value type of the edge weights
1.1454 + typedef typename WeightMap::Value Value;
1.1455 +
1.1456 + /// The type of the matching map
1.1457 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.1458 + MatchingMap;
1.1459 +
1.1460 + /// \brief Scaling factor for primal solution
1.1461 + ///
1.1462 + /// Scaling factor for primal solution. It is equal to 2 or 1
1.1463 + /// according to the value type.
1.1464 + static const int primalScale =
1.1465 + std::numeric_limits<Value>::is_integer ? 2 : 1;
1.1466 +
1.1467 + /// \brief Scaling factor for dual solution
1.1468 + ///
1.1469 + /// Scaling factor for dual solution. It is equal to 4 or 1
1.1470 + /// according to the value type.
1.1471 + static const int dualScale =
1.1472 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.1473 +
1.1474 + private:
1.1475 +
1.1476 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.1477 +
1.1478 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.1479 +
1.1480 + const Graph& _graph;
1.1481 + const WeightMap& _weight;
1.1482 +
1.1483 + MatchingMap* _matching;
1.1484 + NodePotential* _node_potential;
1.1485 +
1.1486 + int _node_num;
1.1487 + bool _allow_loops;
1.1488 +
1.1489 + enum Status {
1.1490 + EVEN = -1, MATCHED = 0, ODD = 1
1.1491 + };
1.1492 +
1.1493 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.1494 + StatusMap* _status;
1.1495 +
1.1496 + typedef typename Graph::template NodeMap<Arc> PredMap;
1.1497 + PredMap* _pred;
1.1498 +
1.1499 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.1500 +
1.1501 + IntNodeMap *_tree_set_index;
1.1502 + TreeSet *_tree_set;
1.1503 +
1.1504 + IntNodeMap *_delta2_index;
1.1505 + BinHeap<Value, IntNodeMap> *_delta2;
1.1506 +
1.1507 + IntEdgeMap *_delta3_index;
1.1508 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1509 +
1.1510 + Value _delta_sum;
1.1511 +
1.1512 + void createStructures() {
1.1513 + _node_num = countNodes(_graph);
1.1514 +
1.1515 + if (!_matching) {
1.1516 + _matching = new MatchingMap(_graph);
1.1517 + }
1.1518 + if (!_node_potential) {
1.1519 + _node_potential = new NodePotential(_graph);
1.1520 + }
1.1521 + if (!_status) {
1.1522 + _status = new StatusMap(_graph);
1.1523 + }
1.1524 + if (!_pred) {
1.1525 + _pred = new PredMap(_graph);
1.1526 + }
1.1527 + if (!_tree_set) {
1.1528 + _tree_set_index = new IntNodeMap(_graph);
1.1529 + _tree_set = new TreeSet(*_tree_set_index);
1.1530 + }
1.1531 + if (!_delta2) {
1.1532 + _delta2_index = new IntNodeMap(_graph);
1.1533 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1.1534 + }
1.1535 + if (!_delta3) {
1.1536 + _delta3_index = new IntEdgeMap(_graph);
1.1537 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1538 + }
1.1539 + }
1.1540 +
1.1541 + void destroyStructures() {
1.1542 + if (_matching) {
1.1543 + delete _matching;
1.1544 + }
1.1545 + if (_node_potential) {
1.1546 + delete _node_potential;
1.1547 + }
1.1548 + if (_status) {
1.1549 + delete _status;
1.1550 + }
1.1551 + if (_pred) {
1.1552 + delete _pred;
1.1553 + }
1.1554 + if (_tree_set) {
1.1555 + delete _tree_set_index;
1.1556 + delete _tree_set;
1.1557 + }
1.1558 + if (_delta2) {
1.1559 + delete _delta2_index;
1.1560 + delete _delta2;
1.1561 + }
1.1562 + if (_delta3) {
1.1563 + delete _delta3_index;
1.1564 + delete _delta3;
1.1565 + }
1.1566 + }
1.1567 +
1.1568 + void matchedToEven(Node node, int tree) {
1.1569 + _tree_set->insert(node, tree);
1.1570 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.1571 +
1.1572 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.1573 + _delta2->erase(node);
1.1574 + }
1.1575 +
1.1576 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1577 + Node v = _graph.source(a);
1.1578 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1579 + dualScale * _weight[a];
1.1580 + if (node == v) {
1.1581 + if (_allow_loops && _graph.direction(a)) {
1.1582 + _delta3->push(a, rw / 2);
1.1583 + }
1.1584 + } else if ((*_status)[v] == EVEN) {
1.1585 + _delta3->push(a, rw / 2);
1.1586 + } else if ((*_status)[v] == MATCHED) {
1.1587 + if (_delta2->state(v) != _delta2->IN_HEAP) {
1.1588 + _pred->set(v, a);
1.1589 + _delta2->push(v, rw);
1.1590 + } else if ((*_delta2)[v] > rw) {
1.1591 + _pred->set(v, a);
1.1592 + _delta2->decrease(v, rw);
1.1593 + }
1.1594 + }
1.1595 + }
1.1596 + }
1.1597 +
1.1598 + void matchedToOdd(Node node, int tree) {
1.1599 + _tree_set->insert(node, tree);
1.1600 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.1601 +
1.1602 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.1603 + _delta2->erase(node);
1.1604 + }
1.1605 + }
1.1606 +
1.1607 + void evenToMatched(Node node, int tree) {
1.1608 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.1609 + Arc min = INVALID;
1.1610 + Value minrw = std::numeric_limits<Value>::max();
1.1611 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1612 + Node v = _graph.source(a);
1.1613 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1614 + dualScale * _weight[a];
1.1615 +
1.1616 + if (node == v) {
1.1617 + if (_allow_loops && _graph.direction(a)) {
1.1618 + _delta3->erase(a);
1.1619 + }
1.1620 + } else if ((*_status)[v] == EVEN) {
1.1621 + _delta3->erase(a);
1.1622 + if (minrw > rw) {
1.1623 + min = _graph.oppositeArc(a);
1.1624 + minrw = rw;
1.1625 + }
1.1626 + } else if ((*_status)[v] == MATCHED) {
1.1627 + if ((*_pred)[v] == a) {
1.1628 + Arc mina = INVALID;
1.1629 + Value minrwa = std::numeric_limits<Value>::max();
1.1630 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1.1631 + Node va = _graph.target(aa);
1.1632 + if ((*_status)[va] != EVEN ||
1.1633 + _tree_set->find(va) == tree) continue;
1.1634 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1.1635 + dualScale * _weight[aa];
1.1636 + if (minrwa > rwa) {
1.1637 + minrwa = rwa;
1.1638 + mina = aa;
1.1639 + }
1.1640 + }
1.1641 + if (mina != INVALID) {
1.1642 + _pred->set(v, mina);
1.1643 + _delta2->increase(v, minrwa);
1.1644 + } else {
1.1645 + _pred->set(v, INVALID);
1.1646 + _delta2->erase(v);
1.1647 + }
1.1648 + }
1.1649 + }
1.1650 + }
1.1651 + if (min != INVALID) {
1.1652 + _pred->set(node, min);
1.1653 + _delta2->push(node, minrw);
1.1654 + } else {
1.1655 + _pred->set(node, INVALID);
1.1656 + }
1.1657 + }
1.1658 +
1.1659 + void oddToMatched(Node node) {
1.1660 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.1661 + Arc min = INVALID;
1.1662 + Value minrw = std::numeric_limits<Value>::max();
1.1663 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1664 + Node v = _graph.source(a);
1.1665 + if ((*_status)[v] != EVEN) continue;
1.1666 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1667 + dualScale * _weight[a];
1.1668 +
1.1669 + if (minrw > rw) {
1.1670 + min = _graph.oppositeArc(a);
1.1671 + minrw = rw;
1.1672 + }
1.1673 + }
1.1674 + if (min != INVALID) {
1.1675 + _pred->set(node, min);
1.1676 + _delta2->push(node, minrw);
1.1677 + } else {
1.1678 + _pred->set(node, INVALID);
1.1679 + }
1.1680 + }
1.1681 +
1.1682 + void alternatePath(Node even, int tree) {
1.1683 + Node odd;
1.1684 +
1.1685 + _status->set(even, MATCHED);
1.1686 + evenToMatched(even, tree);
1.1687 +
1.1688 + Arc prev = (*_matching)[even];
1.1689 + while (prev != INVALID) {
1.1690 + odd = _graph.target(prev);
1.1691 + even = _graph.target((*_pred)[odd]);
1.1692 + _matching->set(odd, (*_pred)[odd]);
1.1693 + _status->set(odd, MATCHED);
1.1694 + oddToMatched(odd);
1.1695 +
1.1696 + prev = (*_matching)[even];
1.1697 + _status->set(even, MATCHED);
1.1698 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.1699 + evenToMatched(even, tree);
1.1700 + }
1.1701 + }
1.1702 +
1.1703 + void destroyTree(int tree) {
1.1704 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1.1705 + if ((*_status)[n] == EVEN) {
1.1706 + _status->set(n, MATCHED);
1.1707 + evenToMatched(n, tree);
1.1708 + } else if ((*_status)[n] == ODD) {
1.1709 + _status->set(n, MATCHED);
1.1710 + oddToMatched(n);
1.1711 + }
1.1712 + }
1.1713 + _tree_set->eraseClass(tree);
1.1714 + }
1.1715 +
1.1716 + void augmentOnEdge(const Edge& edge) {
1.1717 + Node left = _graph.u(edge);
1.1718 + int left_tree = _tree_set->find(left);
1.1719 +
1.1720 + alternatePath(left, left_tree);
1.1721 + destroyTree(left_tree);
1.1722 + _matching->set(left, _graph.direct(edge, true));
1.1723 +
1.1724 + Node right = _graph.v(edge);
1.1725 + int right_tree = _tree_set->find(right);
1.1726 +
1.1727 + alternatePath(right, right_tree);
1.1728 + destroyTree(right_tree);
1.1729 + _matching->set(right, _graph.direct(edge, false));
1.1730 + }
1.1731 +
1.1732 + void augmentOnArc(const Arc& arc) {
1.1733 + Node left = _graph.source(arc);
1.1734 + _status->set(left, MATCHED);
1.1735 + _matching->set(left, arc);
1.1736 + _pred->set(left, arc);
1.1737 +
1.1738 + Node right = _graph.target(arc);
1.1739 + int right_tree = _tree_set->find(right);
1.1740 +
1.1741 + alternatePath(right, right_tree);
1.1742 + destroyTree(right_tree);
1.1743 + _matching->set(right, _graph.oppositeArc(arc));
1.1744 + }
1.1745 +
1.1746 + void extendOnArc(const Arc& arc) {
1.1747 + Node base = _graph.target(arc);
1.1748 + int tree = _tree_set->find(base);
1.1749 +
1.1750 + Node odd = _graph.source(arc);
1.1751 + _tree_set->insert(odd, tree);
1.1752 + _status->set(odd, ODD);
1.1753 + matchedToOdd(odd, tree);
1.1754 + _pred->set(odd, arc);
1.1755 +
1.1756 + Node even = _graph.target((*_matching)[odd]);
1.1757 + _tree_set->insert(even, tree);
1.1758 + _status->set(even, EVEN);
1.1759 + matchedToEven(even, tree);
1.1760 + }
1.1761 +
1.1762 + void cycleOnEdge(const Edge& edge, int tree) {
1.1763 + Node nca = INVALID;
1.1764 + std::vector<Node> left_path, right_path;
1.1765 +
1.1766 + {
1.1767 + std::set<Node> left_set, right_set;
1.1768 + Node left = _graph.u(edge);
1.1769 + left_path.push_back(left);
1.1770 + left_set.insert(left);
1.1771 +
1.1772 + Node right = _graph.v(edge);
1.1773 + right_path.push_back(right);
1.1774 + right_set.insert(right);
1.1775 +
1.1776 + while (true) {
1.1777 +
1.1778 + if (left_set.find(right) != left_set.end()) {
1.1779 + nca = right;
1.1780 + break;
1.1781 + }
1.1782 +
1.1783 + if ((*_matching)[left] == INVALID) break;
1.1784 +
1.1785 + left = _graph.target((*_matching)[left]);
1.1786 + left_path.push_back(left);
1.1787 + left = _graph.target((*_pred)[left]);
1.1788 + left_path.push_back(left);
1.1789 +
1.1790 + left_set.insert(left);
1.1791 +
1.1792 + if (right_set.find(left) != right_set.end()) {
1.1793 + nca = left;
1.1794 + break;
1.1795 + }
1.1796 +
1.1797 + if ((*_matching)[right] == INVALID) break;
1.1798 +
1.1799 + right = _graph.target((*_matching)[right]);
1.1800 + right_path.push_back(right);
1.1801 + right = _graph.target((*_pred)[right]);
1.1802 + right_path.push_back(right);
1.1803 +
1.1804 + right_set.insert(right);
1.1805 +
1.1806 + }
1.1807 +
1.1808 + if (nca == INVALID) {
1.1809 + if ((*_matching)[left] == INVALID) {
1.1810 + nca = right;
1.1811 + while (left_set.find(nca) == left_set.end()) {
1.1812 + nca = _graph.target((*_matching)[nca]);
1.1813 + right_path.push_back(nca);
1.1814 + nca = _graph.target((*_pred)[nca]);
1.1815 + right_path.push_back(nca);
1.1816 + }
1.1817 + } else {
1.1818 + nca = left;
1.1819 + while (right_set.find(nca) == right_set.end()) {
1.1820 + nca = _graph.target((*_matching)[nca]);
1.1821 + left_path.push_back(nca);
1.1822 + nca = _graph.target((*_pred)[nca]);
1.1823 + left_path.push_back(nca);
1.1824 + }
1.1825 + }
1.1826 + }
1.1827 + }
1.1828 +
1.1829 + alternatePath(nca, tree);
1.1830 + Arc prev;
1.1831 +
1.1832 + prev = _graph.direct(edge, true);
1.1833 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1834 + _matching->set(left_path[i], prev);
1.1835 + _status->set(left_path[i], MATCHED);
1.1836 + evenToMatched(left_path[i], tree);
1.1837 +
1.1838 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1.1839 + _status->set(left_path[i + 1], MATCHED);
1.1840 + oddToMatched(left_path[i + 1]);
1.1841 + }
1.1842 + _matching->set(nca, prev);
1.1843 +
1.1844 + for (int i = 0; right_path[i] != nca; i += 2) {
1.1845 + _status->set(right_path[i], MATCHED);
1.1846 + evenToMatched(right_path[i], tree);
1.1847 +
1.1848 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1.1849 + _status->set(right_path[i + 1], MATCHED);
1.1850 + oddToMatched(right_path[i + 1]);
1.1851 + }
1.1852 +
1.1853 + destroyTree(tree);
1.1854 + }
1.1855 +
1.1856 + void extractCycle(const Arc &arc) {
1.1857 + Node left = _graph.source(arc);
1.1858 + Node odd = _graph.target((*_matching)[left]);
1.1859 + Arc prev;
1.1860 + while (odd != left) {
1.1861 + Node even = _graph.target((*_matching)[odd]);
1.1862 + prev = (*_matching)[odd];
1.1863 + odd = _graph.target((*_matching)[even]);
1.1864 + _matching->set(even, _graph.oppositeArc(prev));
1.1865 + }
1.1866 + _matching->set(left, arc);
1.1867 +
1.1868 + Node right = _graph.target(arc);
1.1869 + int right_tree = _tree_set->find(right);
1.1870 + alternatePath(right, right_tree);
1.1871 + destroyTree(right_tree);
1.1872 + _matching->set(right, _graph.oppositeArc(arc));
1.1873 + }
1.1874 +
1.1875 + public:
1.1876 +
1.1877 + /// \brief Constructor
1.1878 + ///
1.1879 + /// Constructor.
1.1880 + MaxWeightedPerfectFractionalMatching(const Graph& graph,
1.1881 + const WeightMap& weight,
1.1882 + bool allow_loops = true)
1.1883 + : _graph(graph), _weight(weight), _matching(0),
1.1884 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1.1885 + _status(0), _pred(0),
1.1886 + _tree_set_index(0), _tree_set(0),
1.1887 +
1.1888 + _delta2_index(0), _delta2(0),
1.1889 + _delta3_index(0), _delta3(0),
1.1890 +
1.1891 + _delta_sum() {}
1.1892 +
1.1893 + ~MaxWeightedPerfectFractionalMatching() {
1.1894 + destroyStructures();
1.1895 + }
1.1896 +
1.1897 + /// \name Execution Control
1.1898 + /// The simplest way to execute the algorithm is to use the
1.1899 + /// \ref run() member function.
1.1900 +
1.1901 + ///@{
1.1902 +
1.1903 + /// \brief Initialize the algorithm
1.1904 + ///
1.1905 + /// This function initializes the algorithm.
1.1906 + void init() {
1.1907 + createStructures();
1.1908 +
1.1909 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1910 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
1.1911 + }
1.1912 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1913 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.1914 + }
1.1915 +
1.1916 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1917 + Value max = - std::numeric_limits<Value>::max();
1.1918 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1919 + if (_graph.target(e) == n && !_allow_loops) continue;
1.1920 + if ((dualScale * _weight[e]) / 2 > max) {
1.1921 + max = (dualScale * _weight[e]) / 2;
1.1922 + }
1.1923 + }
1.1924 + _node_potential->set(n, max);
1.1925 +
1.1926 + _tree_set->insert(n);
1.1927 +
1.1928 + _matching->set(n, INVALID);
1.1929 + _status->set(n, EVEN);
1.1930 + }
1.1931 +
1.1932 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1933 + Node left = _graph.u(e);
1.1934 + Node right = _graph.v(e);
1.1935 + if (left == right && !_allow_loops) continue;
1.1936 + _delta3->push(e, ((*_node_potential)[left] +
1.1937 + (*_node_potential)[right] -
1.1938 + dualScale * _weight[e]) / 2);
1.1939 + }
1.1940 + }
1.1941 +
1.1942 + /// \brief Start the algorithm
1.1943 + ///
1.1944 + /// This function starts the algorithm.
1.1945 + ///
1.1946 + /// \pre \ref init() must be called before using this function.
1.1947 + bool start() {
1.1948 + enum OpType {
1.1949 + D2, D3
1.1950 + };
1.1951 +
1.1952 + int unmatched = _node_num;
1.1953 + while (unmatched > 0) {
1.1954 + Value d2 = !_delta2->empty() ?
1.1955 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1956 +
1.1957 + Value d3 = !_delta3->empty() ?
1.1958 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1959 +
1.1960 + _delta_sum = d3; OpType ot = D3;
1.1961 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1962 +
1.1963 + if (_delta_sum == std::numeric_limits<Value>::max()) {
1.1964 + return false;
1.1965 + }
1.1966 +
1.1967 + switch (ot) {
1.1968 + case D2:
1.1969 + {
1.1970 + Node n = _delta2->top();
1.1971 + Arc a = (*_pred)[n];
1.1972 + if ((*_matching)[n] == INVALID) {
1.1973 + augmentOnArc(a);
1.1974 + --unmatched;
1.1975 + } else {
1.1976 + Node v = _graph.target((*_matching)[n]);
1.1977 + if ((*_matching)[n] !=
1.1978 + _graph.oppositeArc((*_matching)[v])) {
1.1979 + extractCycle(a);
1.1980 + --unmatched;
1.1981 + } else {
1.1982 + extendOnArc(a);
1.1983 + }
1.1984 + }
1.1985 + } break;
1.1986 + case D3:
1.1987 + {
1.1988 + Edge e = _delta3->top();
1.1989 +
1.1990 + Node left = _graph.u(e);
1.1991 + Node right = _graph.v(e);
1.1992 +
1.1993 + int left_tree = _tree_set->find(left);
1.1994 + int right_tree = _tree_set->find(right);
1.1995 +
1.1996 + if (left_tree == right_tree) {
1.1997 + cycleOnEdge(e, left_tree);
1.1998 + --unmatched;
1.1999 + } else {
1.2000 + augmentOnEdge(e);
1.2001 + unmatched -= 2;
1.2002 + }
1.2003 + } break;
1.2004 + }
1.2005 + }
1.2006 + return true;
1.2007 + }
1.2008 +
1.2009 + /// \brief Run the algorithm.
1.2010 + ///
1.2011 + /// This method runs the \c %MaxWeightedMatching algorithm.
1.2012 + ///
1.2013 + /// \note mwfm.run() is just a shortcut of the following code.
1.2014 + /// \code
1.2015 + /// mwpfm.init();
1.2016 + /// mwpfm.start();
1.2017 + /// \endcode
1.2018 + bool run() {
1.2019 + init();
1.2020 + return start();
1.2021 + }
1.2022 +
1.2023 + /// @}
1.2024 +
1.2025 + /// \name Primal Solution
1.2026 + /// Functions to get the primal solution, i.e. the maximum weighted
1.2027 + /// matching.\n
1.2028 + /// Either \ref run() or \ref start() function should be called before
1.2029 + /// using them.
1.2030 +
1.2031 + /// @{
1.2032 +
1.2033 + /// \brief Return the weight of the matching.
1.2034 + ///
1.2035 + /// This function returns the weight of the found matching. This
1.2036 + /// value is scaled by \ref primalScale "primal scale".
1.2037 + ///
1.2038 + /// \pre Either run() or start() must be called before using this function.
1.2039 + Value matchingWeight() const {
1.2040 + Value sum = 0;
1.2041 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2042 + if ((*_matching)[n] != INVALID) {
1.2043 + sum += _weight[(*_matching)[n]];
1.2044 + }
1.2045 + }
1.2046 + return sum * primalScale / 2;
1.2047 + }
1.2048 +
1.2049 + /// \brief Return the number of covered nodes in the matching.
1.2050 + ///
1.2051 + /// This function returns the number of covered nodes in the matching.
1.2052 + ///
1.2053 + /// \pre Either run() or start() must be called before using this function.
1.2054 + int matchingSize() const {
1.2055 + int num = 0;
1.2056 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2057 + if ((*_matching)[n] != INVALID) {
1.2058 + ++num;
1.2059 + }
1.2060 + }
1.2061 + return num;
1.2062 + }
1.2063 +
1.2064 + /// \brief Return \c true if the given edge is in the matching.
1.2065 + ///
1.2066 + /// This function returns \c true if the given edge is in the
1.2067 + /// found matching. The result is scaled by \ref primalScale
1.2068 + /// "primal scale".
1.2069 + ///
1.2070 + /// \pre Either run() or start() must be called before using this function.
1.2071 + Value matching(const Edge& edge) const {
1.2072 + return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1.2073 + * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
1.2074 + * primalScale / 2;
1.2075 + }
1.2076 +
1.2077 + /// \brief Return the fractional matching arc (or edge) incident
1.2078 + /// to the given node.
1.2079 + ///
1.2080 + /// This function returns one of the fractional matching arc (or
1.2081 + /// edge) incident to the given node in the found matching or \c
1.2082 + /// INVALID if the node is not covered by the matching or if the
1.2083 + /// node is on an odd length cycle then it is the successor edge
1.2084 + /// on the cycle.
1.2085 + ///
1.2086 + /// \pre Either run() or start() must be called before using this function.
1.2087 + Arc matching(const Node& node) const {
1.2088 + return (*_matching)[node];
1.2089 + }
1.2090 +
1.2091 + /// \brief Return a const reference to the matching map.
1.2092 + ///
1.2093 + /// This function returns a const reference to a node map that stores
1.2094 + /// the matching arc (or edge) incident to each node.
1.2095 + const MatchingMap& matchingMap() const {
1.2096 + return *_matching;
1.2097 + }
1.2098 +
1.2099 + /// @}
1.2100 +
1.2101 + /// \name Dual Solution
1.2102 + /// Functions to get the dual solution.\n
1.2103 + /// Either \ref run() or \ref start() function should be called before
1.2104 + /// using them.
1.2105 +
1.2106 + /// @{
1.2107 +
1.2108 + /// \brief Return the value of the dual solution.
1.2109 + ///
1.2110 + /// This function returns the value of the dual solution.
1.2111 + /// It should be equal to the primal value scaled by \ref dualScale
1.2112 + /// "dual scale".
1.2113 + ///
1.2114 + /// \pre Either run() or start() must be called before using this function.
1.2115 + Value dualValue() const {
1.2116 + Value sum = 0;
1.2117 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2118 + sum += nodeValue(n);
1.2119 + }
1.2120 + return sum;
1.2121 + }
1.2122 +
1.2123 + /// \brief Return the dual value (potential) of the given node.
1.2124 + ///
1.2125 + /// This function returns the dual value (potential) of the given node.
1.2126 + ///
1.2127 + /// \pre Either run() or start() must be called before using this function.
1.2128 + Value nodeValue(const Node& n) const {
1.2129 + return (*_node_potential)[n];
1.2130 + }
1.2131 +
1.2132 + /// @}
1.2133 +
1.2134 + };
1.2135 +
1.2136 +} //END OF NAMESPACE LEMON
1.2137 +
1.2138 +#endif //LEMON_FRACTIONAL_MATCHING_H