1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/fractional_matching.h Tue Mar 16 21:27:35 2010 +0100
1.3 @@ -0,0 +1,2130 @@
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 MaxFractionalMatching::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 uses priority queues and
1.639 + /// provides \f$O(nm\log n)\f$ time complexity.
1.640 + ///
1.641 + /// The maximum weighted fractional matching is a relaxation of the
1.642 + /// maximum weighted matching problem where the odd set constraints
1.643 + /// are omitted.
1.644 + /// It can be formulated with the following linear program.
1.645 + /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
1.646 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.647 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.648 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.649 + /// \f$X\f$. The result must be the union of a matching with one
1.650 + /// value edges and a set of odd length cycles with half value edges.
1.651 + ///
1.652 + /// The algorithm calculates an optimal fractional matching and a
1.653 + /// proof of the optimality. The solution of the dual problem can be
1.654 + /// used to check the result of the algorithm. The dual linear
1.655 + /// problem is the following.
1.656 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1.657 + /// \f[y_u \ge 0 \quad \forall u \in V\f]
1.658 + /// \f[\min \sum_{u \in V}y_u \f]
1.659 + ///
1.660 + /// The algorithm can be executed with the run() function.
1.661 + /// After it the matching (the primal solution) and the dual solution
1.662 + /// can be obtained using the query functions.
1.663 + ///
1.664 + /// The primal solution is multiplied by
1.665 + /// \ref MaxWeightedFractionalMatching::primalScale "2".
1.666 + /// If the value type is integer, then the dual
1.667 + /// solution is scaled by
1.668 + /// \ref MaxWeightedFractionalMatching::dualScale "4".
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.
1.696 + static const int primalScale = 2;
1.697 +
1.698 + /// \brief Scaling factor for dual solution
1.699 + ///
1.700 + /// Scaling factor for dual solution. It is equal to 4 or 1
1.701 + /// according to the value type.
1.702 + static const int dualScale =
1.703 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.704 +
1.705 + private:
1.706 +
1.707 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.708 +
1.709 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.710 +
1.711 + const Graph& _graph;
1.712 + const WeightMap& _weight;
1.713 +
1.714 + MatchingMap* _matching;
1.715 + NodePotential* _node_potential;
1.716 +
1.717 + int _node_num;
1.718 + bool _allow_loops;
1.719 +
1.720 + enum Status {
1.721 + EVEN = -1, MATCHED = 0, ODD = 1
1.722 + };
1.723 +
1.724 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.725 + StatusMap* _status;
1.726 +
1.727 + typedef typename Graph::template NodeMap<Arc> PredMap;
1.728 + PredMap* _pred;
1.729 +
1.730 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.731 +
1.732 + IntNodeMap *_tree_set_index;
1.733 + TreeSet *_tree_set;
1.734 +
1.735 + IntNodeMap *_delta1_index;
1.736 + BinHeap<Value, IntNodeMap> *_delta1;
1.737 +
1.738 + IntNodeMap *_delta2_index;
1.739 + BinHeap<Value, IntNodeMap> *_delta2;
1.740 +
1.741 + IntEdgeMap *_delta3_index;
1.742 + BinHeap<Value, IntEdgeMap> *_delta3;
1.743 +
1.744 + Value _delta_sum;
1.745 +
1.746 + void createStructures() {
1.747 + _node_num = countNodes(_graph);
1.748 +
1.749 + if (!_matching) {
1.750 + _matching = new MatchingMap(_graph);
1.751 + }
1.752 + if (!_node_potential) {
1.753 + _node_potential = new NodePotential(_graph);
1.754 + }
1.755 + if (!_status) {
1.756 + _status = new StatusMap(_graph);
1.757 + }
1.758 + if (!_pred) {
1.759 + _pred = new PredMap(_graph);
1.760 + }
1.761 + if (!_tree_set) {
1.762 + _tree_set_index = new IntNodeMap(_graph);
1.763 + _tree_set = new TreeSet(*_tree_set_index);
1.764 + }
1.765 + if (!_delta1) {
1.766 + _delta1_index = new IntNodeMap(_graph);
1.767 + _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
1.768 + }
1.769 + if (!_delta2) {
1.770 + _delta2_index = new IntNodeMap(_graph);
1.771 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1.772 + }
1.773 + if (!_delta3) {
1.774 + _delta3_index = new IntEdgeMap(_graph);
1.775 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.776 + }
1.777 + }
1.778 +
1.779 + void destroyStructures() {
1.780 + if (_matching) {
1.781 + delete _matching;
1.782 + }
1.783 + if (_node_potential) {
1.784 + delete _node_potential;
1.785 + }
1.786 + if (_status) {
1.787 + delete _status;
1.788 + }
1.789 + if (_pred) {
1.790 + delete _pred;
1.791 + }
1.792 + if (_tree_set) {
1.793 + delete _tree_set_index;
1.794 + delete _tree_set;
1.795 + }
1.796 + if (_delta1) {
1.797 + delete _delta1_index;
1.798 + delete _delta1;
1.799 + }
1.800 + if (_delta2) {
1.801 + delete _delta2_index;
1.802 + delete _delta2;
1.803 + }
1.804 + if (_delta3) {
1.805 + delete _delta3_index;
1.806 + delete _delta3;
1.807 + }
1.808 + }
1.809 +
1.810 + void matchedToEven(Node node, int tree) {
1.811 + _tree_set->insert(node, tree);
1.812 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.813 + _delta1->push(node, (*_node_potential)[node]);
1.814 +
1.815 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.816 + _delta2->erase(node);
1.817 + }
1.818 +
1.819 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.820 + Node v = _graph.source(a);
1.821 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.822 + dualScale * _weight[a];
1.823 + if (node == v) {
1.824 + if (_allow_loops && _graph.direction(a)) {
1.825 + _delta3->push(a, rw / 2);
1.826 + }
1.827 + } else if ((*_status)[v] == EVEN) {
1.828 + _delta3->push(a, rw / 2);
1.829 + } else if ((*_status)[v] == MATCHED) {
1.830 + if (_delta2->state(v) != _delta2->IN_HEAP) {
1.831 + _pred->set(v, a);
1.832 + _delta2->push(v, rw);
1.833 + } else if ((*_delta2)[v] > rw) {
1.834 + _pred->set(v, a);
1.835 + _delta2->decrease(v, rw);
1.836 + }
1.837 + }
1.838 + }
1.839 + }
1.840 +
1.841 + void matchedToOdd(Node node, int tree) {
1.842 + _tree_set->insert(node, tree);
1.843 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.844 +
1.845 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.846 + _delta2->erase(node);
1.847 + }
1.848 + }
1.849 +
1.850 + void evenToMatched(Node node, int tree) {
1.851 + _delta1->erase(node);
1.852 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.853 + Arc min = INVALID;
1.854 + Value minrw = std::numeric_limits<Value>::max();
1.855 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.856 + Node v = _graph.source(a);
1.857 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.858 + dualScale * _weight[a];
1.859 +
1.860 + if (node == v) {
1.861 + if (_allow_loops && _graph.direction(a)) {
1.862 + _delta3->erase(a);
1.863 + }
1.864 + } else if ((*_status)[v] == EVEN) {
1.865 + _delta3->erase(a);
1.866 + if (minrw > rw) {
1.867 + min = _graph.oppositeArc(a);
1.868 + minrw = rw;
1.869 + }
1.870 + } else if ((*_status)[v] == MATCHED) {
1.871 + if ((*_pred)[v] == a) {
1.872 + Arc mina = INVALID;
1.873 + Value minrwa = std::numeric_limits<Value>::max();
1.874 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1.875 + Node va = _graph.target(aa);
1.876 + if ((*_status)[va] != EVEN ||
1.877 + _tree_set->find(va) == tree) continue;
1.878 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1.879 + dualScale * _weight[aa];
1.880 + if (minrwa > rwa) {
1.881 + minrwa = rwa;
1.882 + mina = aa;
1.883 + }
1.884 + }
1.885 + if (mina != INVALID) {
1.886 + _pred->set(v, mina);
1.887 + _delta2->increase(v, minrwa);
1.888 + } else {
1.889 + _pred->set(v, INVALID);
1.890 + _delta2->erase(v);
1.891 + }
1.892 + }
1.893 + }
1.894 + }
1.895 + if (min != INVALID) {
1.896 + _pred->set(node, min);
1.897 + _delta2->push(node, minrw);
1.898 + } else {
1.899 + _pred->set(node, INVALID);
1.900 + }
1.901 + }
1.902 +
1.903 + void oddToMatched(Node node) {
1.904 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.905 + Arc min = INVALID;
1.906 + Value minrw = std::numeric_limits<Value>::max();
1.907 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.908 + Node v = _graph.source(a);
1.909 + if ((*_status)[v] != EVEN) continue;
1.910 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.911 + dualScale * _weight[a];
1.912 +
1.913 + if (minrw > rw) {
1.914 + min = _graph.oppositeArc(a);
1.915 + minrw = rw;
1.916 + }
1.917 + }
1.918 + if (min != INVALID) {
1.919 + _pred->set(node, min);
1.920 + _delta2->push(node, minrw);
1.921 + } else {
1.922 + _pred->set(node, INVALID);
1.923 + }
1.924 + }
1.925 +
1.926 + void alternatePath(Node even, int tree) {
1.927 + Node odd;
1.928 +
1.929 + _status->set(even, MATCHED);
1.930 + evenToMatched(even, tree);
1.931 +
1.932 + Arc prev = (*_matching)[even];
1.933 + while (prev != INVALID) {
1.934 + odd = _graph.target(prev);
1.935 + even = _graph.target((*_pred)[odd]);
1.936 + _matching->set(odd, (*_pred)[odd]);
1.937 + _status->set(odd, MATCHED);
1.938 + oddToMatched(odd);
1.939 +
1.940 + prev = (*_matching)[even];
1.941 + _status->set(even, MATCHED);
1.942 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.943 + evenToMatched(even, tree);
1.944 + }
1.945 + }
1.946 +
1.947 + void destroyTree(int tree) {
1.948 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1.949 + if ((*_status)[n] == EVEN) {
1.950 + _status->set(n, MATCHED);
1.951 + evenToMatched(n, tree);
1.952 + } else if ((*_status)[n] == ODD) {
1.953 + _status->set(n, MATCHED);
1.954 + oddToMatched(n);
1.955 + }
1.956 + }
1.957 + _tree_set->eraseClass(tree);
1.958 + }
1.959 +
1.960 +
1.961 + void unmatchNode(const Node& node) {
1.962 + int tree = _tree_set->find(node);
1.963 +
1.964 + alternatePath(node, tree);
1.965 + destroyTree(tree);
1.966 +
1.967 + _matching->set(node, INVALID);
1.968 + }
1.969 +
1.970 +
1.971 + void augmentOnEdge(const Edge& edge) {
1.972 + Node left = _graph.u(edge);
1.973 + int left_tree = _tree_set->find(left);
1.974 +
1.975 + alternatePath(left, left_tree);
1.976 + destroyTree(left_tree);
1.977 + _matching->set(left, _graph.direct(edge, true));
1.978 +
1.979 + Node right = _graph.v(edge);
1.980 + int right_tree = _tree_set->find(right);
1.981 +
1.982 + alternatePath(right, right_tree);
1.983 + destroyTree(right_tree);
1.984 + _matching->set(right, _graph.direct(edge, false));
1.985 + }
1.986 +
1.987 + void augmentOnArc(const Arc& arc) {
1.988 + Node left = _graph.source(arc);
1.989 + _status->set(left, MATCHED);
1.990 + _matching->set(left, arc);
1.991 + _pred->set(left, arc);
1.992 +
1.993 + Node right = _graph.target(arc);
1.994 + int right_tree = _tree_set->find(right);
1.995 +
1.996 + alternatePath(right, right_tree);
1.997 + destroyTree(right_tree);
1.998 + _matching->set(right, _graph.oppositeArc(arc));
1.999 + }
1.1000 +
1.1001 + void extendOnArc(const Arc& arc) {
1.1002 + Node base = _graph.target(arc);
1.1003 + int tree = _tree_set->find(base);
1.1004 +
1.1005 + Node odd = _graph.source(arc);
1.1006 + _tree_set->insert(odd, tree);
1.1007 + _status->set(odd, ODD);
1.1008 + matchedToOdd(odd, tree);
1.1009 + _pred->set(odd, arc);
1.1010 +
1.1011 + Node even = _graph.target((*_matching)[odd]);
1.1012 + _tree_set->insert(even, tree);
1.1013 + _status->set(even, EVEN);
1.1014 + matchedToEven(even, tree);
1.1015 + }
1.1016 +
1.1017 + void cycleOnEdge(const Edge& edge, int tree) {
1.1018 + Node nca = INVALID;
1.1019 + std::vector<Node> left_path, right_path;
1.1020 +
1.1021 + {
1.1022 + std::set<Node> left_set, right_set;
1.1023 + Node left = _graph.u(edge);
1.1024 + left_path.push_back(left);
1.1025 + left_set.insert(left);
1.1026 +
1.1027 + Node right = _graph.v(edge);
1.1028 + right_path.push_back(right);
1.1029 + right_set.insert(right);
1.1030 +
1.1031 + while (true) {
1.1032 +
1.1033 + if (left_set.find(right) != left_set.end()) {
1.1034 + nca = right;
1.1035 + break;
1.1036 + }
1.1037 +
1.1038 + if ((*_matching)[left] == INVALID) break;
1.1039 +
1.1040 + left = _graph.target((*_matching)[left]);
1.1041 + left_path.push_back(left);
1.1042 + left = _graph.target((*_pred)[left]);
1.1043 + left_path.push_back(left);
1.1044 +
1.1045 + left_set.insert(left);
1.1046 +
1.1047 + if (right_set.find(left) != right_set.end()) {
1.1048 + nca = left;
1.1049 + break;
1.1050 + }
1.1051 +
1.1052 + if ((*_matching)[right] == INVALID) break;
1.1053 +
1.1054 + right = _graph.target((*_matching)[right]);
1.1055 + right_path.push_back(right);
1.1056 + right = _graph.target((*_pred)[right]);
1.1057 + right_path.push_back(right);
1.1058 +
1.1059 + right_set.insert(right);
1.1060 +
1.1061 + }
1.1062 +
1.1063 + if (nca == INVALID) {
1.1064 + if ((*_matching)[left] == INVALID) {
1.1065 + nca = right;
1.1066 + while (left_set.find(nca) == left_set.end()) {
1.1067 + nca = _graph.target((*_matching)[nca]);
1.1068 + right_path.push_back(nca);
1.1069 + nca = _graph.target((*_pred)[nca]);
1.1070 + right_path.push_back(nca);
1.1071 + }
1.1072 + } else {
1.1073 + nca = left;
1.1074 + while (right_set.find(nca) == right_set.end()) {
1.1075 + nca = _graph.target((*_matching)[nca]);
1.1076 + left_path.push_back(nca);
1.1077 + nca = _graph.target((*_pred)[nca]);
1.1078 + left_path.push_back(nca);
1.1079 + }
1.1080 + }
1.1081 + }
1.1082 + }
1.1083 +
1.1084 + alternatePath(nca, tree);
1.1085 + Arc prev;
1.1086 +
1.1087 + prev = _graph.direct(edge, true);
1.1088 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1089 + _matching->set(left_path[i], prev);
1.1090 + _status->set(left_path[i], MATCHED);
1.1091 + evenToMatched(left_path[i], tree);
1.1092 +
1.1093 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1.1094 + _status->set(left_path[i + 1], MATCHED);
1.1095 + oddToMatched(left_path[i + 1]);
1.1096 + }
1.1097 + _matching->set(nca, prev);
1.1098 +
1.1099 + for (int i = 0; right_path[i] != nca; i += 2) {
1.1100 + _status->set(right_path[i], MATCHED);
1.1101 + evenToMatched(right_path[i], tree);
1.1102 +
1.1103 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1.1104 + _status->set(right_path[i + 1], MATCHED);
1.1105 + oddToMatched(right_path[i + 1]);
1.1106 + }
1.1107 +
1.1108 + destroyTree(tree);
1.1109 + }
1.1110 +
1.1111 + void extractCycle(const Arc &arc) {
1.1112 + Node left = _graph.source(arc);
1.1113 + Node odd = _graph.target((*_matching)[left]);
1.1114 + Arc prev;
1.1115 + while (odd != left) {
1.1116 + Node even = _graph.target((*_matching)[odd]);
1.1117 + prev = (*_matching)[odd];
1.1118 + odd = _graph.target((*_matching)[even]);
1.1119 + _matching->set(even, _graph.oppositeArc(prev));
1.1120 + }
1.1121 + _matching->set(left, arc);
1.1122 +
1.1123 + Node right = _graph.target(arc);
1.1124 + int right_tree = _tree_set->find(right);
1.1125 + alternatePath(right, right_tree);
1.1126 + destroyTree(right_tree);
1.1127 + _matching->set(right, _graph.oppositeArc(arc));
1.1128 + }
1.1129 +
1.1130 + public:
1.1131 +
1.1132 + /// \brief Constructor
1.1133 + ///
1.1134 + /// Constructor.
1.1135 + MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
1.1136 + bool allow_loops = true)
1.1137 + : _graph(graph), _weight(weight), _matching(0),
1.1138 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1.1139 + _status(0), _pred(0),
1.1140 + _tree_set_index(0), _tree_set(0),
1.1141 +
1.1142 + _delta1_index(0), _delta1(0),
1.1143 + _delta2_index(0), _delta2(0),
1.1144 + _delta3_index(0), _delta3(0),
1.1145 +
1.1146 + _delta_sum() {}
1.1147 +
1.1148 + ~MaxWeightedFractionalMatching() {
1.1149 + destroyStructures();
1.1150 + }
1.1151 +
1.1152 + /// \name Execution Control
1.1153 + /// The simplest way to execute the algorithm is to use the
1.1154 + /// \ref run() member function.
1.1155 +
1.1156 + ///@{
1.1157 +
1.1158 + /// \brief Initialize the algorithm
1.1159 + ///
1.1160 + /// This function initializes the algorithm.
1.1161 + void init() {
1.1162 + createStructures();
1.1163 +
1.1164 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1165 + (*_delta1_index)[n] = _delta1->PRE_HEAP;
1.1166 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
1.1167 + }
1.1168 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1169 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.1170 + }
1.1171 +
1.1172 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1173 + Value max = 0;
1.1174 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1175 + if (_graph.target(e) == n && !_allow_loops) continue;
1.1176 + if ((dualScale * _weight[e]) / 2 > max) {
1.1177 + max = (dualScale * _weight[e]) / 2;
1.1178 + }
1.1179 + }
1.1180 + _node_potential->set(n, max);
1.1181 + _delta1->push(n, max);
1.1182 +
1.1183 + _tree_set->insert(n);
1.1184 +
1.1185 + _matching->set(n, INVALID);
1.1186 + _status->set(n, EVEN);
1.1187 + }
1.1188 +
1.1189 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1190 + Node left = _graph.u(e);
1.1191 + Node right = _graph.v(e);
1.1192 + if (left == right && !_allow_loops) continue;
1.1193 + _delta3->push(e, ((*_node_potential)[left] +
1.1194 + (*_node_potential)[right] -
1.1195 + dualScale * _weight[e]) / 2);
1.1196 + }
1.1197 + }
1.1198 +
1.1199 + /// \brief Start the algorithm
1.1200 + ///
1.1201 + /// This function starts the algorithm.
1.1202 + ///
1.1203 + /// \pre \ref init() must be called before using this function.
1.1204 + void start() {
1.1205 + enum OpType {
1.1206 + D1, D2, D3
1.1207 + };
1.1208 +
1.1209 + int unmatched = _node_num;
1.1210 + while (unmatched > 0) {
1.1211 + Value d1 = !_delta1->empty() ?
1.1212 + _delta1->prio() : std::numeric_limits<Value>::max();
1.1213 +
1.1214 + Value d2 = !_delta2->empty() ?
1.1215 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1216 +
1.1217 + Value d3 = !_delta3->empty() ?
1.1218 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1219 +
1.1220 + _delta_sum = d3; OpType ot = D3;
1.1221 + if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1.1222 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1223 +
1.1224 + switch (ot) {
1.1225 + case D1:
1.1226 + {
1.1227 + Node n = _delta1->top();
1.1228 + unmatchNode(n);
1.1229 + --unmatched;
1.1230 + }
1.1231 + break;
1.1232 + case D2:
1.1233 + {
1.1234 + Node n = _delta2->top();
1.1235 + Arc a = (*_pred)[n];
1.1236 + if ((*_matching)[n] == INVALID) {
1.1237 + augmentOnArc(a);
1.1238 + --unmatched;
1.1239 + } else {
1.1240 + Node v = _graph.target((*_matching)[n]);
1.1241 + if ((*_matching)[n] !=
1.1242 + _graph.oppositeArc((*_matching)[v])) {
1.1243 + extractCycle(a);
1.1244 + --unmatched;
1.1245 + } else {
1.1246 + extendOnArc(a);
1.1247 + }
1.1248 + }
1.1249 + } break;
1.1250 + case D3:
1.1251 + {
1.1252 + Edge e = _delta3->top();
1.1253 +
1.1254 + Node left = _graph.u(e);
1.1255 + Node right = _graph.v(e);
1.1256 +
1.1257 + int left_tree = _tree_set->find(left);
1.1258 + int right_tree = _tree_set->find(right);
1.1259 +
1.1260 + if (left_tree == right_tree) {
1.1261 + cycleOnEdge(e, left_tree);
1.1262 + --unmatched;
1.1263 + } else {
1.1264 + augmentOnEdge(e);
1.1265 + unmatched -= 2;
1.1266 + }
1.1267 + } break;
1.1268 + }
1.1269 + }
1.1270 + }
1.1271 +
1.1272 + /// \brief Run the algorithm.
1.1273 + ///
1.1274 + /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
1.1275 + ///
1.1276 + /// \note mwfm.run() is just a shortcut of the following code.
1.1277 + /// \code
1.1278 + /// mwfm.init();
1.1279 + /// mwfm.start();
1.1280 + /// \endcode
1.1281 + void run() {
1.1282 + init();
1.1283 + start();
1.1284 + }
1.1285 +
1.1286 + /// @}
1.1287 +
1.1288 + /// \name Primal Solution
1.1289 + /// Functions to get the primal solution, i.e. the maximum weighted
1.1290 + /// matching.\n
1.1291 + /// Either \ref run() or \ref start() function should be called before
1.1292 + /// using them.
1.1293 +
1.1294 + /// @{
1.1295 +
1.1296 + /// \brief Return the weight of the matching.
1.1297 + ///
1.1298 + /// This function returns the weight of the found matching. This
1.1299 + /// value is scaled by \ref primalScale "primal scale".
1.1300 + ///
1.1301 + /// \pre Either run() or start() must be called before using this function.
1.1302 + Value matchingWeight() const {
1.1303 + Value sum = 0;
1.1304 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1305 + if ((*_matching)[n] != INVALID) {
1.1306 + sum += _weight[(*_matching)[n]];
1.1307 + }
1.1308 + }
1.1309 + return sum * primalScale / 2;
1.1310 + }
1.1311 +
1.1312 + /// \brief Return the number of covered nodes in the matching.
1.1313 + ///
1.1314 + /// This function returns the number of covered nodes in the matching.
1.1315 + ///
1.1316 + /// \pre Either run() or start() must be called before using this function.
1.1317 + int matchingSize() const {
1.1318 + int num = 0;
1.1319 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1320 + if ((*_matching)[n] != INVALID) {
1.1321 + ++num;
1.1322 + }
1.1323 + }
1.1324 + return num;
1.1325 + }
1.1326 +
1.1327 + /// \brief Return \c true if the given edge is in the matching.
1.1328 + ///
1.1329 + /// This function returns \c true if the given edge is in the
1.1330 + /// found matching. The result is scaled by \ref primalScale
1.1331 + /// "primal scale".
1.1332 + ///
1.1333 + /// \pre Either run() or start() must be called before using this function.
1.1334 + int matching(const Edge& edge) const {
1.1335 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1.1336 + + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
1.1337 + }
1.1338 +
1.1339 + /// \brief Return the fractional matching arc (or edge) incident
1.1340 + /// to the given node.
1.1341 + ///
1.1342 + /// This function returns one of the fractional matching arc (or
1.1343 + /// edge) incident to the given node in the found matching or \c
1.1344 + /// INVALID if the node is not covered by the matching or if the
1.1345 + /// node is on an odd length cycle then it is the successor edge
1.1346 + /// on the cycle.
1.1347 + ///
1.1348 + /// \pre Either run() or start() must be called before using this function.
1.1349 + Arc matching(const Node& node) const {
1.1350 + return (*_matching)[node];
1.1351 + }
1.1352 +
1.1353 + /// \brief Return a const reference to the matching map.
1.1354 + ///
1.1355 + /// This function returns a const reference to a node map that stores
1.1356 + /// the matching arc (or edge) incident to each node.
1.1357 + const MatchingMap& matchingMap() const {
1.1358 + return *_matching;
1.1359 + }
1.1360 +
1.1361 + /// @}
1.1362 +
1.1363 + /// \name Dual Solution
1.1364 + /// Functions to get the dual solution.\n
1.1365 + /// Either \ref run() or \ref start() function should be called before
1.1366 + /// using them.
1.1367 +
1.1368 + /// @{
1.1369 +
1.1370 + /// \brief Return the value of the dual solution.
1.1371 + ///
1.1372 + /// This function returns the value of the dual solution.
1.1373 + /// It should be equal to the primal value scaled by \ref dualScale
1.1374 + /// "dual scale".
1.1375 + ///
1.1376 + /// \pre Either run() or start() must be called before using this function.
1.1377 + Value dualValue() const {
1.1378 + Value sum = 0;
1.1379 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1380 + sum += nodeValue(n);
1.1381 + }
1.1382 + return sum;
1.1383 + }
1.1384 +
1.1385 + /// \brief Return the dual value (potential) of the given node.
1.1386 + ///
1.1387 + /// This function returns the dual value (potential) of the given node.
1.1388 + ///
1.1389 + /// \pre Either run() or start() must be called before using this function.
1.1390 + Value nodeValue(const Node& n) const {
1.1391 + return (*_node_potential)[n];
1.1392 + }
1.1393 +
1.1394 + /// @}
1.1395 +
1.1396 + };
1.1397 +
1.1398 + /// \ingroup matching
1.1399 + ///
1.1400 + /// \brief Weighted fractional perfect matching in general graphs
1.1401 + ///
1.1402 + /// This class provides an efficient implementation of fractional
1.1403 + /// matching algorithm. The implementation uses priority queues and
1.1404 + /// provides \f$O(nm\log n)\f$ time complexity.
1.1405 + ///
1.1406 + /// The maximum weighted fractional perfect matching is a relaxation
1.1407 + /// of the maximum weighted perfect matching problem where the odd
1.1408 + /// set constraints are omitted.
1.1409 + /// It can be formulated with the following linear program.
1.1410 + /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1.1411 + /// \f[x_e \ge 0\quad \forall e\in E\f]
1.1412 + /// \f[\max \sum_{e\in E}x_ew_e\f]
1.1413 + /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1.1414 + /// \f$X\f$. The result must be the union of a matching with one
1.1415 + /// value edges and a set of odd length cycles with half value edges.
1.1416 + ///
1.1417 + /// The algorithm calculates an optimal fractional matching and a
1.1418 + /// proof of the optimality. The solution of the dual problem can be
1.1419 + /// used to check the result of the algorithm. The dual linear
1.1420 + /// problem is the following.
1.1421 + /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1.1422 + /// \f[\min \sum_{u \in V}y_u \f]
1.1423 + ///
1.1424 + /// The algorithm can be executed with the run() function.
1.1425 + /// After it the matching (the primal solution) and the dual solution
1.1426 + /// can be obtained using the query functions.
1.1427 + ///
1.1428 + /// The primal solution is multiplied by
1.1429 + /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
1.1430 + /// If the value type is integer, then the dual
1.1431 + /// solution is scaled by
1.1432 + /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
1.1433 + ///
1.1434 + /// \tparam GR The undirected graph type the algorithm runs on.
1.1435 + /// \tparam WM The type edge weight map. The default type is
1.1436 + /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.1437 +#ifdef DOXYGEN
1.1438 + template <typename GR, typename WM>
1.1439 +#else
1.1440 + template <typename GR,
1.1441 + typename WM = typename GR::template EdgeMap<int> >
1.1442 +#endif
1.1443 + class MaxWeightedPerfectFractionalMatching {
1.1444 + public:
1.1445 +
1.1446 + /// The graph type of the algorithm
1.1447 + typedef GR Graph;
1.1448 + /// The type of the edge weight map
1.1449 + typedef WM WeightMap;
1.1450 + /// The value type of the edge weights
1.1451 + typedef typename WeightMap::Value Value;
1.1452 +
1.1453 + /// The type of the matching map
1.1454 + typedef typename Graph::template NodeMap<typename Graph::Arc>
1.1455 + MatchingMap;
1.1456 +
1.1457 + /// \brief Scaling factor for primal solution
1.1458 + ///
1.1459 + /// Scaling factor for primal solution.
1.1460 + static const int primalScale = 2;
1.1461 +
1.1462 + /// \brief Scaling factor for dual solution
1.1463 + ///
1.1464 + /// Scaling factor for dual solution. It is equal to 4 or 1
1.1465 + /// according to the value type.
1.1466 + static const int dualScale =
1.1467 + std::numeric_limits<Value>::is_integer ? 4 : 1;
1.1468 +
1.1469 + private:
1.1470 +
1.1471 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.1472 +
1.1473 + typedef typename Graph::template NodeMap<Value> NodePotential;
1.1474 +
1.1475 + const Graph& _graph;
1.1476 + const WeightMap& _weight;
1.1477 +
1.1478 + MatchingMap* _matching;
1.1479 + NodePotential* _node_potential;
1.1480 +
1.1481 + int _node_num;
1.1482 + bool _allow_loops;
1.1483 +
1.1484 + enum Status {
1.1485 + EVEN = -1, MATCHED = 0, ODD = 1
1.1486 + };
1.1487 +
1.1488 + typedef typename Graph::template NodeMap<Status> StatusMap;
1.1489 + StatusMap* _status;
1.1490 +
1.1491 + typedef typename Graph::template NodeMap<Arc> PredMap;
1.1492 + PredMap* _pred;
1.1493 +
1.1494 + typedef ExtendFindEnum<IntNodeMap> TreeSet;
1.1495 +
1.1496 + IntNodeMap *_tree_set_index;
1.1497 + TreeSet *_tree_set;
1.1498 +
1.1499 + IntNodeMap *_delta2_index;
1.1500 + BinHeap<Value, IntNodeMap> *_delta2;
1.1501 +
1.1502 + IntEdgeMap *_delta3_index;
1.1503 + BinHeap<Value, IntEdgeMap> *_delta3;
1.1504 +
1.1505 + Value _delta_sum;
1.1506 +
1.1507 + void createStructures() {
1.1508 + _node_num = countNodes(_graph);
1.1509 +
1.1510 + if (!_matching) {
1.1511 + _matching = new MatchingMap(_graph);
1.1512 + }
1.1513 + if (!_node_potential) {
1.1514 + _node_potential = new NodePotential(_graph);
1.1515 + }
1.1516 + if (!_status) {
1.1517 + _status = new StatusMap(_graph);
1.1518 + }
1.1519 + if (!_pred) {
1.1520 + _pred = new PredMap(_graph);
1.1521 + }
1.1522 + if (!_tree_set) {
1.1523 + _tree_set_index = new IntNodeMap(_graph);
1.1524 + _tree_set = new TreeSet(*_tree_set_index);
1.1525 + }
1.1526 + if (!_delta2) {
1.1527 + _delta2_index = new IntNodeMap(_graph);
1.1528 + _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1.1529 + }
1.1530 + if (!_delta3) {
1.1531 + _delta3_index = new IntEdgeMap(_graph);
1.1532 + _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1.1533 + }
1.1534 + }
1.1535 +
1.1536 + void destroyStructures() {
1.1537 + if (_matching) {
1.1538 + delete _matching;
1.1539 + }
1.1540 + if (_node_potential) {
1.1541 + delete _node_potential;
1.1542 + }
1.1543 + if (_status) {
1.1544 + delete _status;
1.1545 + }
1.1546 + if (_pred) {
1.1547 + delete _pred;
1.1548 + }
1.1549 + if (_tree_set) {
1.1550 + delete _tree_set_index;
1.1551 + delete _tree_set;
1.1552 + }
1.1553 + if (_delta2) {
1.1554 + delete _delta2_index;
1.1555 + delete _delta2;
1.1556 + }
1.1557 + if (_delta3) {
1.1558 + delete _delta3_index;
1.1559 + delete _delta3;
1.1560 + }
1.1561 + }
1.1562 +
1.1563 + void matchedToEven(Node node, int tree) {
1.1564 + _tree_set->insert(node, tree);
1.1565 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.1566 +
1.1567 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.1568 + _delta2->erase(node);
1.1569 + }
1.1570 +
1.1571 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1572 + Node v = _graph.source(a);
1.1573 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1574 + dualScale * _weight[a];
1.1575 + if (node == v) {
1.1576 + if (_allow_loops && _graph.direction(a)) {
1.1577 + _delta3->push(a, rw / 2);
1.1578 + }
1.1579 + } else if ((*_status)[v] == EVEN) {
1.1580 + _delta3->push(a, rw / 2);
1.1581 + } else if ((*_status)[v] == MATCHED) {
1.1582 + if (_delta2->state(v) != _delta2->IN_HEAP) {
1.1583 + _pred->set(v, a);
1.1584 + _delta2->push(v, rw);
1.1585 + } else if ((*_delta2)[v] > rw) {
1.1586 + _pred->set(v, a);
1.1587 + _delta2->decrease(v, rw);
1.1588 + }
1.1589 + }
1.1590 + }
1.1591 + }
1.1592 +
1.1593 + void matchedToOdd(Node node, int tree) {
1.1594 + _tree_set->insert(node, tree);
1.1595 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.1596 +
1.1597 + if (_delta2->state(node) == _delta2->IN_HEAP) {
1.1598 + _delta2->erase(node);
1.1599 + }
1.1600 + }
1.1601 +
1.1602 + void evenToMatched(Node node, int tree) {
1.1603 + _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1.1604 + Arc min = INVALID;
1.1605 + Value minrw = std::numeric_limits<Value>::max();
1.1606 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1607 + Node v = _graph.source(a);
1.1608 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1609 + dualScale * _weight[a];
1.1610 +
1.1611 + if (node == v) {
1.1612 + if (_allow_loops && _graph.direction(a)) {
1.1613 + _delta3->erase(a);
1.1614 + }
1.1615 + } else if ((*_status)[v] == EVEN) {
1.1616 + _delta3->erase(a);
1.1617 + if (minrw > rw) {
1.1618 + min = _graph.oppositeArc(a);
1.1619 + minrw = rw;
1.1620 + }
1.1621 + } else if ((*_status)[v] == MATCHED) {
1.1622 + if ((*_pred)[v] == a) {
1.1623 + Arc mina = INVALID;
1.1624 + Value minrwa = std::numeric_limits<Value>::max();
1.1625 + for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1.1626 + Node va = _graph.target(aa);
1.1627 + if ((*_status)[va] != EVEN ||
1.1628 + _tree_set->find(va) == tree) continue;
1.1629 + Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1.1630 + dualScale * _weight[aa];
1.1631 + if (minrwa > rwa) {
1.1632 + minrwa = rwa;
1.1633 + mina = aa;
1.1634 + }
1.1635 + }
1.1636 + if (mina != INVALID) {
1.1637 + _pred->set(v, mina);
1.1638 + _delta2->increase(v, minrwa);
1.1639 + } else {
1.1640 + _pred->set(v, INVALID);
1.1641 + _delta2->erase(v);
1.1642 + }
1.1643 + }
1.1644 + }
1.1645 + }
1.1646 + if (min != INVALID) {
1.1647 + _pred->set(node, min);
1.1648 + _delta2->push(node, minrw);
1.1649 + } else {
1.1650 + _pred->set(node, INVALID);
1.1651 + }
1.1652 + }
1.1653 +
1.1654 + void oddToMatched(Node node) {
1.1655 + _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1.1656 + Arc min = INVALID;
1.1657 + Value minrw = std::numeric_limits<Value>::max();
1.1658 + for (InArcIt a(_graph, node); a != INVALID; ++a) {
1.1659 + Node v = _graph.source(a);
1.1660 + if ((*_status)[v] != EVEN) continue;
1.1661 + Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1.1662 + dualScale * _weight[a];
1.1663 +
1.1664 + if (minrw > rw) {
1.1665 + min = _graph.oppositeArc(a);
1.1666 + minrw = rw;
1.1667 + }
1.1668 + }
1.1669 + if (min != INVALID) {
1.1670 + _pred->set(node, min);
1.1671 + _delta2->push(node, minrw);
1.1672 + } else {
1.1673 + _pred->set(node, INVALID);
1.1674 + }
1.1675 + }
1.1676 +
1.1677 + void alternatePath(Node even, int tree) {
1.1678 + Node odd;
1.1679 +
1.1680 + _status->set(even, MATCHED);
1.1681 + evenToMatched(even, tree);
1.1682 +
1.1683 + Arc prev = (*_matching)[even];
1.1684 + while (prev != INVALID) {
1.1685 + odd = _graph.target(prev);
1.1686 + even = _graph.target((*_pred)[odd]);
1.1687 + _matching->set(odd, (*_pred)[odd]);
1.1688 + _status->set(odd, MATCHED);
1.1689 + oddToMatched(odd);
1.1690 +
1.1691 + prev = (*_matching)[even];
1.1692 + _status->set(even, MATCHED);
1.1693 + _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1.1694 + evenToMatched(even, tree);
1.1695 + }
1.1696 + }
1.1697 +
1.1698 + void destroyTree(int tree) {
1.1699 + for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1.1700 + if ((*_status)[n] == EVEN) {
1.1701 + _status->set(n, MATCHED);
1.1702 + evenToMatched(n, tree);
1.1703 + } else if ((*_status)[n] == ODD) {
1.1704 + _status->set(n, MATCHED);
1.1705 + oddToMatched(n);
1.1706 + }
1.1707 + }
1.1708 + _tree_set->eraseClass(tree);
1.1709 + }
1.1710 +
1.1711 + void augmentOnEdge(const Edge& edge) {
1.1712 + Node left = _graph.u(edge);
1.1713 + int left_tree = _tree_set->find(left);
1.1714 +
1.1715 + alternatePath(left, left_tree);
1.1716 + destroyTree(left_tree);
1.1717 + _matching->set(left, _graph.direct(edge, true));
1.1718 +
1.1719 + Node right = _graph.v(edge);
1.1720 + int right_tree = _tree_set->find(right);
1.1721 +
1.1722 + alternatePath(right, right_tree);
1.1723 + destroyTree(right_tree);
1.1724 + _matching->set(right, _graph.direct(edge, false));
1.1725 + }
1.1726 +
1.1727 + void augmentOnArc(const Arc& arc) {
1.1728 + Node left = _graph.source(arc);
1.1729 + _status->set(left, MATCHED);
1.1730 + _matching->set(left, arc);
1.1731 + _pred->set(left, arc);
1.1732 +
1.1733 + Node right = _graph.target(arc);
1.1734 + int right_tree = _tree_set->find(right);
1.1735 +
1.1736 + alternatePath(right, right_tree);
1.1737 + destroyTree(right_tree);
1.1738 + _matching->set(right, _graph.oppositeArc(arc));
1.1739 + }
1.1740 +
1.1741 + void extendOnArc(const Arc& arc) {
1.1742 + Node base = _graph.target(arc);
1.1743 + int tree = _tree_set->find(base);
1.1744 +
1.1745 + Node odd = _graph.source(arc);
1.1746 + _tree_set->insert(odd, tree);
1.1747 + _status->set(odd, ODD);
1.1748 + matchedToOdd(odd, tree);
1.1749 + _pred->set(odd, arc);
1.1750 +
1.1751 + Node even = _graph.target((*_matching)[odd]);
1.1752 + _tree_set->insert(even, tree);
1.1753 + _status->set(even, EVEN);
1.1754 + matchedToEven(even, tree);
1.1755 + }
1.1756 +
1.1757 + void cycleOnEdge(const Edge& edge, int tree) {
1.1758 + Node nca = INVALID;
1.1759 + std::vector<Node> left_path, right_path;
1.1760 +
1.1761 + {
1.1762 + std::set<Node> left_set, right_set;
1.1763 + Node left = _graph.u(edge);
1.1764 + left_path.push_back(left);
1.1765 + left_set.insert(left);
1.1766 +
1.1767 + Node right = _graph.v(edge);
1.1768 + right_path.push_back(right);
1.1769 + right_set.insert(right);
1.1770 +
1.1771 + while (true) {
1.1772 +
1.1773 + if (left_set.find(right) != left_set.end()) {
1.1774 + nca = right;
1.1775 + break;
1.1776 + }
1.1777 +
1.1778 + if ((*_matching)[left] == INVALID) break;
1.1779 +
1.1780 + left = _graph.target((*_matching)[left]);
1.1781 + left_path.push_back(left);
1.1782 + left = _graph.target((*_pred)[left]);
1.1783 + left_path.push_back(left);
1.1784 +
1.1785 + left_set.insert(left);
1.1786 +
1.1787 + if (right_set.find(left) != right_set.end()) {
1.1788 + nca = left;
1.1789 + break;
1.1790 + }
1.1791 +
1.1792 + if ((*_matching)[right] == INVALID) break;
1.1793 +
1.1794 + right = _graph.target((*_matching)[right]);
1.1795 + right_path.push_back(right);
1.1796 + right = _graph.target((*_pred)[right]);
1.1797 + right_path.push_back(right);
1.1798 +
1.1799 + right_set.insert(right);
1.1800 +
1.1801 + }
1.1802 +
1.1803 + if (nca == INVALID) {
1.1804 + if ((*_matching)[left] == INVALID) {
1.1805 + nca = right;
1.1806 + while (left_set.find(nca) == left_set.end()) {
1.1807 + nca = _graph.target((*_matching)[nca]);
1.1808 + right_path.push_back(nca);
1.1809 + nca = _graph.target((*_pred)[nca]);
1.1810 + right_path.push_back(nca);
1.1811 + }
1.1812 + } else {
1.1813 + nca = left;
1.1814 + while (right_set.find(nca) == right_set.end()) {
1.1815 + nca = _graph.target((*_matching)[nca]);
1.1816 + left_path.push_back(nca);
1.1817 + nca = _graph.target((*_pred)[nca]);
1.1818 + left_path.push_back(nca);
1.1819 + }
1.1820 + }
1.1821 + }
1.1822 + }
1.1823 +
1.1824 + alternatePath(nca, tree);
1.1825 + Arc prev;
1.1826 +
1.1827 + prev = _graph.direct(edge, true);
1.1828 + for (int i = 0; left_path[i] != nca; i += 2) {
1.1829 + _matching->set(left_path[i], prev);
1.1830 + _status->set(left_path[i], MATCHED);
1.1831 + evenToMatched(left_path[i], tree);
1.1832 +
1.1833 + prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1.1834 + _status->set(left_path[i + 1], MATCHED);
1.1835 + oddToMatched(left_path[i + 1]);
1.1836 + }
1.1837 + _matching->set(nca, prev);
1.1838 +
1.1839 + for (int i = 0; right_path[i] != nca; i += 2) {
1.1840 + _status->set(right_path[i], MATCHED);
1.1841 + evenToMatched(right_path[i], tree);
1.1842 +
1.1843 + _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1.1844 + _status->set(right_path[i + 1], MATCHED);
1.1845 + oddToMatched(right_path[i + 1]);
1.1846 + }
1.1847 +
1.1848 + destroyTree(tree);
1.1849 + }
1.1850 +
1.1851 + void extractCycle(const Arc &arc) {
1.1852 + Node left = _graph.source(arc);
1.1853 + Node odd = _graph.target((*_matching)[left]);
1.1854 + Arc prev;
1.1855 + while (odd != left) {
1.1856 + Node even = _graph.target((*_matching)[odd]);
1.1857 + prev = (*_matching)[odd];
1.1858 + odd = _graph.target((*_matching)[even]);
1.1859 + _matching->set(even, _graph.oppositeArc(prev));
1.1860 + }
1.1861 + _matching->set(left, arc);
1.1862 +
1.1863 + Node right = _graph.target(arc);
1.1864 + int right_tree = _tree_set->find(right);
1.1865 + alternatePath(right, right_tree);
1.1866 + destroyTree(right_tree);
1.1867 + _matching->set(right, _graph.oppositeArc(arc));
1.1868 + }
1.1869 +
1.1870 + public:
1.1871 +
1.1872 + /// \brief Constructor
1.1873 + ///
1.1874 + /// Constructor.
1.1875 + MaxWeightedPerfectFractionalMatching(const Graph& graph,
1.1876 + const WeightMap& weight,
1.1877 + bool allow_loops = true)
1.1878 + : _graph(graph), _weight(weight), _matching(0),
1.1879 + _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1.1880 + _status(0), _pred(0),
1.1881 + _tree_set_index(0), _tree_set(0),
1.1882 +
1.1883 + _delta2_index(0), _delta2(0),
1.1884 + _delta3_index(0), _delta3(0),
1.1885 +
1.1886 + _delta_sum() {}
1.1887 +
1.1888 + ~MaxWeightedPerfectFractionalMatching() {
1.1889 + destroyStructures();
1.1890 + }
1.1891 +
1.1892 + /// \name Execution Control
1.1893 + /// The simplest way to execute the algorithm is to use the
1.1894 + /// \ref run() member function.
1.1895 +
1.1896 + ///@{
1.1897 +
1.1898 + /// \brief Initialize the algorithm
1.1899 + ///
1.1900 + /// This function initializes the algorithm.
1.1901 + void init() {
1.1902 + createStructures();
1.1903 +
1.1904 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1905 + (*_delta2_index)[n] = _delta2->PRE_HEAP;
1.1906 + }
1.1907 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1908 + (*_delta3_index)[e] = _delta3->PRE_HEAP;
1.1909 + }
1.1910 +
1.1911 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.1912 + Value max = - std::numeric_limits<Value>::max();
1.1913 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.1914 + if (_graph.target(e) == n && !_allow_loops) continue;
1.1915 + if ((dualScale * _weight[e]) / 2 > max) {
1.1916 + max = (dualScale * _weight[e]) / 2;
1.1917 + }
1.1918 + }
1.1919 + _node_potential->set(n, max);
1.1920 +
1.1921 + _tree_set->insert(n);
1.1922 +
1.1923 + _matching->set(n, INVALID);
1.1924 + _status->set(n, EVEN);
1.1925 + }
1.1926 +
1.1927 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.1928 + Node left = _graph.u(e);
1.1929 + Node right = _graph.v(e);
1.1930 + if (left == right && !_allow_loops) continue;
1.1931 + _delta3->push(e, ((*_node_potential)[left] +
1.1932 + (*_node_potential)[right] -
1.1933 + dualScale * _weight[e]) / 2);
1.1934 + }
1.1935 + }
1.1936 +
1.1937 + /// \brief Start the algorithm
1.1938 + ///
1.1939 + /// This function starts the algorithm.
1.1940 + ///
1.1941 + /// \pre \ref init() must be called before using this function.
1.1942 + bool start() {
1.1943 + enum OpType {
1.1944 + D2, D3
1.1945 + };
1.1946 +
1.1947 + int unmatched = _node_num;
1.1948 + while (unmatched > 0) {
1.1949 + Value d2 = !_delta2->empty() ?
1.1950 + _delta2->prio() : std::numeric_limits<Value>::max();
1.1951 +
1.1952 + Value d3 = !_delta3->empty() ?
1.1953 + _delta3->prio() : std::numeric_limits<Value>::max();
1.1954 +
1.1955 + _delta_sum = d3; OpType ot = D3;
1.1956 + if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1.1957 +
1.1958 + if (_delta_sum == std::numeric_limits<Value>::max()) {
1.1959 + return false;
1.1960 + }
1.1961 +
1.1962 + switch (ot) {
1.1963 + case D2:
1.1964 + {
1.1965 + Node n = _delta2->top();
1.1966 + Arc a = (*_pred)[n];
1.1967 + if ((*_matching)[n] == INVALID) {
1.1968 + augmentOnArc(a);
1.1969 + --unmatched;
1.1970 + } else {
1.1971 + Node v = _graph.target((*_matching)[n]);
1.1972 + if ((*_matching)[n] !=
1.1973 + _graph.oppositeArc((*_matching)[v])) {
1.1974 + extractCycle(a);
1.1975 + --unmatched;
1.1976 + } else {
1.1977 + extendOnArc(a);
1.1978 + }
1.1979 + }
1.1980 + } break;
1.1981 + case D3:
1.1982 + {
1.1983 + Edge e = _delta3->top();
1.1984 +
1.1985 + Node left = _graph.u(e);
1.1986 + Node right = _graph.v(e);
1.1987 +
1.1988 + int left_tree = _tree_set->find(left);
1.1989 + int right_tree = _tree_set->find(right);
1.1990 +
1.1991 + if (left_tree == right_tree) {
1.1992 + cycleOnEdge(e, left_tree);
1.1993 + --unmatched;
1.1994 + } else {
1.1995 + augmentOnEdge(e);
1.1996 + unmatched -= 2;
1.1997 + }
1.1998 + } break;
1.1999 + }
1.2000 + }
1.2001 + return true;
1.2002 + }
1.2003 +
1.2004 + /// \brief Run the algorithm.
1.2005 + ///
1.2006 + /// This method runs the \c %MaxWeightedPerfectFractionalMatching
1.2007 + /// algorithm.
1.2008 + ///
1.2009 + /// \note mwfm.run() is just a shortcut of the following code.
1.2010 + /// \code
1.2011 + /// mwpfm.init();
1.2012 + /// mwpfm.start();
1.2013 + /// \endcode
1.2014 + bool run() {
1.2015 + init();
1.2016 + return start();
1.2017 + }
1.2018 +
1.2019 + /// @}
1.2020 +
1.2021 + /// \name Primal Solution
1.2022 + /// Functions to get the primal solution, i.e. the maximum weighted
1.2023 + /// matching.\n
1.2024 + /// Either \ref run() or \ref start() function should be called before
1.2025 + /// using them.
1.2026 +
1.2027 + /// @{
1.2028 +
1.2029 + /// \brief Return the weight of the matching.
1.2030 + ///
1.2031 + /// This function returns the weight of the found matching. This
1.2032 + /// value is scaled by \ref primalScale "primal scale".
1.2033 + ///
1.2034 + /// \pre Either run() or start() must be called before using this function.
1.2035 + Value matchingWeight() const {
1.2036 + Value sum = 0;
1.2037 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2038 + if ((*_matching)[n] != INVALID) {
1.2039 + sum += _weight[(*_matching)[n]];
1.2040 + }
1.2041 + }
1.2042 + return sum * primalScale / 2;
1.2043 + }
1.2044 +
1.2045 + /// \brief Return the number of covered nodes in the matching.
1.2046 + ///
1.2047 + /// This function returns the number of covered nodes in the matching.
1.2048 + ///
1.2049 + /// \pre Either run() or start() must be called before using this function.
1.2050 + int matchingSize() const {
1.2051 + int num = 0;
1.2052 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2053 + if ((*_matching)[n] != INVALID) {
1.2054 + ++num;
1.2055 + }
1.2056 + }
1.2057 + return num;
1.2058 + }
1.2059 +
1.2060 + /// \brief Return \c true if the given edge is in the matching.
1.2061 + ///
1.2062 + /// This function returns \c true if the given edge is in the
1.2063 + /// found matching. The result is scaled by \ref primalScale
1.2064 + /// "primal scale".
1.2065 + ///
1.2066 + /// \pre Either run() or start() must be called before using this function.
1.2067 + int matching(const Edge& edge) const {
1.2068 + return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1.2069 + + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
1.2070 + }
1.2071 +
1.2072 + /// \brief Return the fractional matching arc (or edge) incident
1.2073 + /// to the given node.
1.2074 + ///
1.2075 + /// This function returns one of the fractional matching arc (or
1.2076 + /// edge) incident to the given node in the found matching or \c
1.2077 + /// INVALID if the node is not covered by the matching or if the
1.2078 + /// node is on an odd length cycle then it is the successor edge
1.2079 + /// on the cycle.
1.2080 + ///
1.2081 + /// \pre Either run() or start() must be called before using this function.
1.2082 + Arc matching(const Node& node) const {
1.2083 + return (*_matching)[node];
1.2084 + }
1.2085 +
1.2086 + /// \brief Return a const reference to the matching map.
1.2087 + ///
1.2088 + /// This function returns a const reference to a node map that stores
1.2089 + /// the matching arc (or edge) incident to each node.
1.2090 + const MatchingMap& matchingMap() const {
1.2091 + return *_matching;
1.2092 + }
1.2093 +
1.2094 + /// @}
1.2095 +
1.2096 + /// \name Dual Solution
1.2097 + /// Functions to get the dual solution.\n
1.2098 + /// Either \ref run() or \ref start() function should be called before
1.2099 + /// using them.
1.2100 +
1.2101 + /// @{
1.2102 +
1.2103 + /// \brief Return the value of the dual solution.
1.2104 + ///
1.2105 + /// This function returns the value of the dual solution.
1.2106 + /// It should be equal to the primal value scaled by \ref dualScale
1.2107 + /// "dual scale".
1.2108 + ///
1.2109 + /// \pre Either run() or start() must be called before using this function.
1.2110 + Value dualValue() const {
1.2111 + Value sum = 0;
1.2112 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.2113 + sum += nodeValue(n);
1.2114 + }
1.2115 + return sum;
1.2116 + }
1.2117 +
1.2118 + /// \brief Return the dual value (potential) of the given node.
1.2119 + ///
1.2120 + /// This function returns the dual value (potential) of the given node.
1.2121 + ///
1.2122 + /// \pre Either run() or start() must be called before using this function.
1.2123 + Value nodeValue(const Node& n) const {
1.2124 + return (*_node_potential)[n];
1.2125 + }
1.2126 +
1.2127 + /// @}
1.2128 +
1.2129 + };
1.2130 +
1.2131 +} //END OF NAMESPACE LEMON
1.2132 +
1.2133 +#endif //LEMON_FRACTIONAL_MATCHING_H