2140 lines
63 KiB
C++
Executable File
2140 lines
63 KiB
C++
Executable File
/* -*- mode: C++; indent-tabs-mode: nil; -*-
|
|
*
|
|
* This file is a part of LEMON, a generic C++ optimization library.
|
|
*
|
|
* Copyright (C) 2003-2013
|
|
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
|
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
|
*
|
|
* Permission to use, modify and distribute this software is granted
|
|
* provided that this copyright notice appears in all copies. For
|
|
* precise terms see the accompanying LICENSE file.
|
|
*
|
|
* This software is provided "AS IS" with no warranty of any kind,
|
|
* express or implied, and with no claim as to its suitability for any
|
|
* purpose.
|
|
*
|
|
*/
|
|
|
|
#ifndef LEMON_FRACTIONAL_MATCHING_H
|
|
#define LEMON_FRACTIONAL_MATCHING_H
|
|
|
|
#include <vector>
|
|
#include <queue>
|
|
#include <set>
|
|
#include <limits>
|
|
|
|
#include <lemon/core.h>
|
|
#include <lemon/unionfind.h>
|
|
#include <lemon/bin_heap.h>
|
|
#include <lemon/maps.h>
|
|
#include <lemon/assert.h>
|
|
#include <lemon/elevator.h>
|
|
|
|
///\ingroup matching
|
|
///\file
|
|
///\brief Fractional matching algorithms in general graphs.
|
|
|
|
namespace lemon {
|
|
|
|
/// \brief Default traits class of MaxFractionalMatching class.
|
|
///
|
|
/// Default traits class of MaxFractionalMatching class.
|
|
/// \tparam GR Graph type.
|
|
template <typename GR>
|
|
struct MaxFractionalMatchingDefaultTraits {
|
|
|
|
/// \brief The type of the graph the algorithm runs on.
|
|
typedef GR Graph;
|
|
|
|
/// \brief The type of the map that stores the matching.
|
|
///
|
|
/// The type of the map that stores the matching arcs.
|
|
/// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
|
|
typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
|
|
|
|
/// \brief Instantiates a MatchingMap.
|
|
///
|
|
/// This function instantiates a \ref MatchingMap.
|
|
/// \param graph The graph for which we would like to define
|
|
/// the matching map.
|
|
static MatchingMap* createMatchingMap(const Graph& graph) {
|
|
return new MatchingMap(graph);
|
|
}
|
|
|
|
/// \brief The elevator type used by MaxFractionalMatching algorithm.
|
|
///
|
|
/// The elevator type used by MaxFractionalMatching algorithm.
|
|
///
|
|
/// \sa Elevator
|
|
/// \sa LinkedElevator
|
|
typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
|
|
|
|
/// \brief Instantiates an Elevator.
|
|
///
|
|
/// This function instantiates an \ref Elevator.
|
|
/// \param graph The graph for which we would like to define
|
|
/// the elevator.
|
|
/// \param max_level The maximum level of the elevator.
|
|
static Elevator* createElevator(const Graph& graph, int max_level) {
|
|
return new Elevator(graph, max_level);
|
|
}
|
|
};
|
|
|
|
/// \ingroup matching
|
|
///
|
|
/// \brief Max cardinality fractional matching
|
|
///
|
|
/// This class provides an implementation of fractional matching
|
|
/// algorithm based on push-relabel principle.
|
|
///
|
|
/// The maximum cardinality fractional matching is a relaxation of the
|
|
/// maximum cardinality matching problem where the odd set constraints
|
|
/// are omitted.
|
|
/// It can be formulated with the following linear program.
|
|
/// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
|
|
/// \f[x_e \ge 0\quad \forall e\in E\f]
|
|
/// \f[\max \sum_{e\in E}x_e\f]
|
|
/// where \f$\delta(X)\f$ is the set of edges incident to a node in
|
|
/// \f$X\f$. The result can be represented as the union of a
|
|
/// matching with one value edges and a set of odd length cycles
|
|
/// with half value edges.
|
|
///
|
|
/// The algorithm calculates an optimal fractional matching and a
|
|
/// barrier. The number of adjacents of any node set minus the size
|
|
/// of node set is a lower bound on the uncovered nodes in the
|
|
/// graph. For maximum matching a barrier is computed which
|
|
/// maximizes this difference.
|
|
///
|
|
/// The algorithm can be executed with the run() function. After it
|
|
/// the matching (the primal solution) and the barrier (the dual
|
|
/// solution) can be obtained using the query functions.
|
|
///
|
|
/// The primal solution is multiplied by
|
|
/// \ref MaxFractionalMatching::primalScale "2".
|
|
///
|
|
/// \tparam GR The undirected graph type the algorithm runs on.
|
|
#ifdef DOXYGEN
|
|
template <typename GR, typename TR>
|
|
#else
|
|
template <typename GR,
|
|
typename TR = MaxFractionalMatchingDefaultTraits<GR> >
|
|
#endif
|
|
class MaxFractionalMatching {
|
|
public:
|
|
|
|
/// \brief The \ref lemon::MaxFractionalMatchingDefaultTraits
|
|
/// "traits class" of the algorithm.
|
|
typedef TR Traits;
|
|
/// The type of the graph the algorithm runs on.
|
|
typedef typename TR::Graph Graph;
|
|
/// The type of the matching map.
|
|
typedef typename TR::MatchingMap MatchingMap;
|
|
/// The type of the elevator.
|
|
typedef typename TR::Elevator Elevator;
|
|
|
|
/// \brief Scaling factor for primal solution
|
|
///
|
|
/// Scaling factor for primal solution.
|
|
static const int primalScale = 2;
|
|
|
|
private:
|
|
|
|
const Graph &_graph;
|
|
int _node_num;
|
|
bool _allow_loops;
|
|
int _empty_level;
|
|
|
|
TEMPLATE_GRAPH_TYPEDEFS(Graph);
|
|
|
|
bool _local_matching;
|
|
MatchingMap *_matching;
|
|
|
|
bool _local_level;
|
|
Elevator *_level;
|
|
|
|
typedef typename Graph::template NodeMap<int> InDegMap;
|
|
InDegMap *_indeg;
|
|
|
|
void createStructures() {
|
|
_node_num = countNodes(_graph);
|
|
|
|
if (!_matching) {
|
|
_local_matching = true;
|
|
_matching = Traits::createMatchingMap(_graph);
|
|
}
|
|
if (!_level) {
|
|
_local_level = true;
|
|
_level = Traits::createElevator(_graph, _node_num);
|
|
}
|
|
if (!_indeg) {
|
|
_indeg = new InDegMap(_graph);
|
|
}
|
|
}
|
|
|
|
void destroyStructures() {
|
|
if (_local_matching) {
|
|
delete _matching;
|
|
}
|
|
if (_local_level) {
|
|
delete _level;
|
|
}
|
|
if (_indeg) {
|
|
delete _indeg;
|
|
}
|
|
}
|
|
|
|
void postprocessing() {
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_indeg)[n] != 0) continue;
|
|
_indeg->set(n, -1);
|
|
Node u = n;
|
|
while ((*_matching)[u] != INVALID) {
|
|
Node v = _graph.target((*_matching)[u]);
|
|
_indeg->set(v, -1);
|
|
Arc a = _graph.oppositeArc((*_matching)[u]);
|
|
u = _graph.target((*_matching)[v]);
|
|
_indeg->set(u, -1);
|
|
_matching->set(v, a);
|
|
}
|
|
}
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_indeg)[n] != 1) continue;
|
|
_indeg->set(n, -1);
|
|
|
|
int num = 1;
|
|
Node u = _graph.target((*_matching)[n]);
|
|
while (u != n) {
|
|
_indeg->set(u, -1);
|
|
u = _graph.target((*_matching)[u]);
|
|
++num;
|
|
}
|
|
if (num % 2 == 0 && num > 2) {
|
|
Arc prev = _graph.oppositeArc((*_matching)[n]);
|
|
Node v = _graph.target((*_matching)[n]);
|
|
u = _graph.target((*_matching)[v]);
|
|
_matching->set(v, prev);
|
|
while (u != n) {
|
|
prev = _graph.oppositeArc((*_matching)[u]);
|
|
v = _graph.target((*_matching)[u]);
|
|
u = _graph.target((*_matching)[v]);
|
|
_matching->set(v, prev);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
public:
|
|
|
|
typedef MaxFractionalMatching Create;
|
|
|
|
///\name Named Template Parameters
|
|
|
|
///@{
|
|
|
|
template <typename T>
|
|
struct SetMatchingMapTraits : public Traits {
|
|
typedef T MatchingMap;
|
|
static MatchingMap *createMatchingMap(const Graph&) {
|
|
LEMON_ASSERT(false, "MatchingMap is not initialized");
|
|
return 0; // ignore warnings
|
|
}
|
|
};
|
|
|
|
/// \brief \ref named-templ-param "Named parameter" for setting
|
|
/// MatchingMap type
|
|
///
|
|
/// \ref named-templ-param "Named parameter" for setting MatchingMap
|
|
/// type.
|
|
template <typename T>
|
|
struct SetMatchingMap
|
|
: public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
|
|
typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
|
|
};
|
|
|
|
template <typename T>
|
|
struct SetElevatorTraits : public Traits {
|
|
typedef T Elevator;
|
|
static Elevator *createElevator(const Graph&, int) {
|
|
LEMON_ASSERT(false, "Elevator is not initialized");
|
|
return 0; // ignore warnings
|
|
}
|
|
};
|
|
|
|
/// \brief \ref named-templ-param "Named parameter" for setting
|
|
/// Elevator type
|
|
///
|
|
/// \ref named-templ-param "Named parameter" for setting Elevator
|
|
/// type. If this named parameter is used, then an external
|
|
/// elevator object must be passed to the algorithm using the
|
|
/// \ref elevator(Elevator&) "elevator()" function before calling
|
|
/// \ref run() or \ref init().
|
|
/// \sa SetStandardElevator
|
|
template <typename T>
|
|
struct SetElevator
|
|
: public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
|
|
typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
|
|
};
|
|
|
|
template <typename T>
|
|
struct SetStandardElevatorTraits : public Traits {
|
|
typedef T Elevator;
|
|
static Elevator *createElevator(const Graph& graph, int max_level) {
|
|
return new Elevator(graph, max_level);
|
|
}
|
|
};
|
|
|
|
/// \brief \ref named-templ-param "Named parameter" for setting
|
|
/// Elevator type with automatic allocation
|
|
///
|
|
/// \ref named-templ-param "Named parameter" for setting Elevator
|
|
/// type with automatic allocation.
|
|
/// The Elevator should have standard constructor interface to be
|
|
/// able to automatically created by the algorithm (i.e. the
|
|
/// graph and the maximum level should be passed to it).
|
|
/// However an external elevator object could also be passed to the
|
|
/// algorithm with the \ref elevator(Elevator&) "elevator()" function
|
|
/// before calling \ref run() or \ref init().
|
|
/// \sa SetElevator
|
|
template <typename T>
|
|
struct SetStandardElevator
|
|
: public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
|
|
typedef MaxFractionalMatching<Graph,
|
|
SetStandardElevatorTraits<T> > Create;
|
|
};
|
|
|
|
/// @}
|
|
|
|
protected:
|
|
|
|
MaxFractionalMatching() {}
|
|
|
|
public:
|
|
|
|
/// \brief Constructor
|
|
///
|
|
/// Constructor.
|
|
///
|
|
MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
|
|
: _graph(graph), _allow_loops(allow_loops),
|
|
_local_matching(false), _matching(0),
|
|
_local_level(false), _level(0), _indeg(0)
|
|
{}
|
|
|
|
~MaxFractionalMatching() {
|
|
destroyStructures();
|
|
}
|
|
|
|
/// \brief Sets the matching map.
|
|
///
|
|
/// Sets the matching map.
|
|
/// If you don't use this function before calling \ref run() or
|
|
/// \ref init(), an instance will be allocated automatically.
|
|
/// The destructor deallocates this automatically allocated map,
|
|
/// of course.
|
|
/// \return <tt>(*this)</tt>
|
|
MaxFractionalMatching& matchingMap(MatchingMap& map) {
|
|
if (_local_matching) {
|
|
delete _matching;
|
|
_local_matching = false;
|
|
}
|
|
_matching = ↦
|
|
return *this;
|
|
}
|
|
|
|
/// \brief Sets the elevator used by algorithm.
|
|
///
|
|
/// Sets the elevator used by algorithm.
|
|
/// If you don't use this function before calling \ref run() or
|
|
/// \ref init(), an instance will be allocated automatically.
|
|
/// The destructor deallocates this automatically allocated elevator,
|
|
/// of course.
|
|
/// \return <tt>(*this)</tt>
|
|
MaxFractionalMatching& elevator(Elevator& elevator) {
|
|
if (_local_level) {
|
|
delete _level;
|
|
_local_level = false;
|
|
}
|
|
_level = &elevator;
|
|
return *this;
|
|
}
|
|
|
|
/// \brief Returns a const reference to the elevator.
|
|
///
|
|
/// Returns a const reference to the elevator.
|
|
///
|
|
/// \pre Either \ref run() or \ref init() must be called before
|
|
/// using this function.
|
|
const Elevator& elevator() const {
|
|
return *_level;
|
|
}
|
|
|
|
/// \name Execution control
|
|
/// The simplest way to execute the algorithm is to use one of the
|
|
/// member functions called \c run(). \n
|
|
/// If you need more control on the execution, first
|
|
/// you must call \ref init() and then one variant of the start()
|
|
/// member.
|
|
|
|
/// @{
|
|
|
|
/// \brief Initializes the internal data structures.
|
|
///
|
|
/// Initializes the internal data structures and sets the initial
|
|
/// matching.
|
|
void init() {
|
|
createStructures();
|
|
|
|
_level->initStart();
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
_indeg->set(n, 0);
|
|
_matching->set(n, INVALID);
|
|
_level->initAddItem(n);
|
|
}
|
|
_level->initFinish();
|
|
|
|
_empty_level = _node_num;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
for (OutArcIt a(_graph, n); a != INVALID; ++a) {
|
|
if (_graph.target(a) == n && !_allow_loops) continue;
|
|
_matching->set(n, a);
|
|
Node v = _graph.target((*_matching)[n]);
|
|
_indeg->set(v, (*_indeg)[v] + 1);
|
|
break;
|
|
}
|
|
}
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_indeg)[n] == 0) {
|
|
_level->activate(n);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// \brief Starts the algorithm and computes a fractional matching
|
|
///
|
|
/// The algorithm computes a maximum fractional matching.
|
|
///
|
|
/// \param postprocess The algorithm computes first a matching
|
|
/// which is a union of a matching with one value edges, cycles
|
|
/// with half value edges and even length paths with half value
|
|
/// edges. If the parameter is true, then after the push-relabel
|
|
/// algorithm it postprocesses the matching to contain only
|
|
/// matching edges and half value odd cycles.
|
|
void start(bool postprocess = true) {
|
|
Node n;
|
|
while ((n = _level->highestActive()) != INVALID) {
|
|
int level = _level->highestActiveLevel();
|
|
int new_level = _level->maxLevel();
|
|
for (InArcIt a(_graph, n); a != INVALID; ++a) {
|
|
Node u = _graph.source(a);
|
|
if (n == u && !_allow_loops) continue;
|
|
Node v = _graph.target((*_matching)[u]);
|
|
if ((*_level)[v] < level) {
|
|
_indeg->set(v, (*_indeg)[v] - 1);
|
|
if ((*_indeg)[v] == 0) {
|
|
_level->activate(v);
|
|
}
|
|
_matching->set(u, a);
|
|
_indeg->set(n, (*_indeg)[n] + 1);
|
|
_level->deactivate(n);
|
|
goto no_more_push;
|
|
} else if (new_level > (*_level)[v]) {
|
|
new_level = (*_level)[v];
|
|
}
|
|
}
|
|
|
|
if (new_level + 1 < _level->maxLevel()) {
|
|
_level->liftHighestActive(new_level + 1);
|
|
} else {
|
|
_level->liftHighestActiveToTop();
|
|
}
|
|
if (_level->emptyLevel(level)) {
|
|
_level->liftToTop(level);
|
|
}
|
|
no_more_push:
|
|
;
|
|
}
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] == INVALID) continue;
|
|
Node u = _graph.target((*_matching)[n]);
|
|
if ((*_indeg)[u] > 1) {
|
|
_indeg->set(u, (*_indeg)[u] - 1);
|
|
_matching->set(n, INVALID);
|
|
}
|
|
}
|
|
if (postprocess) {
|
|
postprocessing();
|
|
}
|
|
}
|
|
|
|
/// \brief Starts the algorithm and computes a perfect fractional
|
|
/// matching
|
|
///
|
|
/// The algorithm computes a perfect fractional matching. If it
|
|
/// does not exists, then the algorithm returns false and the
|
|
/// matching is undefined and the barrier.
|
|
///
|
|
/// \param postprocess The algorithm computes first a matching
|
|
/// which is a union of a matching with one value edges, cycles
|
|
/// with half value edges and even length paths with half value
|
|
/// edges. If the parameter is true, then after the push-relabel
|
|
/// algorithm it postprocesses the matching to contain only
|
|
/// matching edges and half value odd cycles.
|
|
bool startPerfect(bool postprocess = true) {
|
|
Node n;
|
|
while ((n = _level->highestActive()) != INVALID) {
|
|
int level = _level->highestActiveLevel();
|
|
int new_level = _level->maxLevel();
|
|
for (InArcIt a(_graph, n); a != INVALID; ++a) {
|
|
Node u = _graph.source(a);
|
|
if (n == u && !_allow_loops) continue;
|
|
Node v = _graph.target((*_matching)[u]);
|
|
if ((*_level)[v] < level) {
|
|
_indeg->set(v, (*_indeg)[v] - 1);
|
|
if ((*_indeg)[v] == 0) {
|
|
_level->activate(v);
|
|
}
|
|
_matching->set(u, a);
|
|
_indeg->set(n, (*_indeg)[n] + 1);
|
|
_level->deactivate(n);
|
|
goto no_more_push;
|
|
} else if (new_level > (*_level)[v]) {
|
|
new_level = (*_level)[v];
|
|
}
|
|
}
|
|
|
|
if (new_level + 1 < _level->maxLevel()) {
|
|
_level->liftHighestActive(new_level + 1);
|
|
} else {
|
|
_level->liftHighestActiveToTop();
|
|
_empty_level = _level->maxLevel() - 1;
|
|
return false;
|
|
}
|
|
if (_level->emptyLevel(level)) {
|
|
_level->liftToTop(level);
|
|
_empty_level = level;
|
|
return false;
|
|
}
|
|
no_more_push:
|
|
;
|
|
}
|
|
if (postprocess) {
|
|
postprocessing();
|
|
}
|
|
return true;
|
|
}
|
|
|
|
/// \brief Runs the algorithm
|
|
///
|
|
/// Just a shortcut for the next code:
|
|
///\code
|
|
/// init();
|
|
/// start();
|
|
///\endcode
|
|
void run(bool postprocess = true) {
|
|
init();
|
|
start(postprocess);
|
|
}
|
|
|
|
/// \brief Runs the algorithm to find a perfect fractional matching
|
|
///
|
|
/// Just a shortcut for the next code:
|
|
///\code
|
|
/// init();
|
|
/// startPerfect();
|
|
///\endcode
|
|
bool runPerfect(bool postprocess = true) {
|
|
init();
|
|
return startPerfect(postprocess);
|
|
}
|
|
|
|
///@}
|
|
|
|
/// \name Query Functions
|
|
/// The result of the %Matching algorithm can be obtained using these
|
|
/// functions.\n
|
|
/// Before the use of these functions,
|
|
/// either run() or start() must be called.
|
|
///@{
|
|
|
|
|
|
/// \brief Return the number of covered nodes in the matching.
|
|
///
|
|
/// This function returns the number of covered nodes in the matching.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matchingSize() const {
|
|
int num = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] != INVALID) {
|
|
++num;
|
|
}
|
|
}
|
|
return num;
|
|
}
|
|
|
|
/// \brief Returns a const reference to the matching map.
|
|
///
|
|
/// Returns a const reference to the node map storing the found
|
|
/// fractional matching. This method can be called after
|
|
/// running the algorithm.
|
|
///
|
|
/// \pre Either \ref run() or \ref init() must be called before
|
|
/// using this function.
|
|
const MatchingMap& matchingMap() const {
|
|
return *_matching;
|
|
}
|
|
|
|
/// \brief Return \c true if the given edge is in the matching.
|
|
///
|
|
/// This function returns \c true if the given edge is in the
|
|
/// found matching. The result is scaled by \ref primalScale
|
|
/// "primal scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matching(const Edge& edge) const {
|
|
return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
|
|
(edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
|
|
}
|
|
|
|
/// \brief Return the fractional matching arc (or edge) incident
|
|
/// to the given node.
|
|
///
|
|
/// This function returns one of the fractional matching arc (or
|
|
/// edge) incident to the given node in the found matching or \c
|
|
/// INVALID if the node is not covered by the matching or if the
|
|
/// node is on an odd length cycle then it is the successor edge
|
|
/// on the cycle.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Arc matching(const Node& node) const {
|
|
return (*_matching)[node];
|
|
}
|
|
|
|
/// \brief Returns true if the node is in the barrier
|
|
///
|
|
/// The barrier is a subset of the nodes. If the nodes in the
|
|
/// barrier have less adjacent nodes than the size of the barrier,
|
|
/// then at least as much nodes cannot be covered as the
|
|
/// difference of the two subsets.
|
|
bool barrier(const Node& node) const {
|
|
return (*_level)[node] >= _empty_level;
|
|
}
|
|
|
|
/// @}
|
|
|
|
};
|
|
|
|
/// \ingroup matching
|
|
///
|
|
/// \brief Weighted fractional matching in general graphs
|
|
///
|
|
/// This class provides an efficient implementation of fractional
|
|
/// matching algorithm. The implementation uses priority queues and
|
|
/// provides \f$O(nm\log n)\f$ time complexity.
|
|
///
|
|
/// The maximum weighted fractional matching is a relaxation of the
|
|
/// maximum weighted matching problem where the odd set constraints
|
|
/// are omitted.
|
|
/// It can be formulated with the following linear program.
|
|
/// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
|
|
/// \f[x_e \ge 0\quad \forall e\in E\f]
|
|
/// \f[\max \sum_{e\in E}x_ew_e\f]
|
|
/// where \f$\delta(X)\f$ is the set of edges incident to a node in
|
|
/// \f$X\f$. The result must be the union of a matching with one
|
|
/// value edges and a set of odd length cycles with half value edges.
|
|
///
|
|
/// The algorithm calculates an optimal fractional matching and a
|
|
/// proof of the optimality. The solution of the dual problem can be
|
|
/// used to check the result of the algorithm. The dual linear
|
|
/// problem is the following.
|
|
/// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
|
|
/// \f[y_u \ge 0 \quad \forall u \in V\f]
|
|
/// \f[\min \sum_{u \in V}y_u \f]
|
|
///
|
|
/// The algorithm can be executed with the run() function.
|
|
/// After it the matching (the primal solution) and the dual solution
|
|
/// can be obtained using the query functions.
|
|
///
|
|
/// The primal solution is multiplied by
|
|
/// \ref MaxWeightedFractionalMatching::primalScale "2".
|
|
/// If the value type is integer, then the dual
|
|
/// solution is scaled by
|
|
/// \ref MaxWeightedFractionalMatching::dualScale "4".
|
|
///
|
|
/// \tparam GR The undirected graph type the algorithm runs on.
|
|
/// \tparam WM The type edge weight map. The default type is
|
|
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
|
|
#ifdef DOXYGEN
|
|
template <typename GR, typename WM>
|
|
#else
|
|
template <typename GR,
|
|
typename WM = typename GR::template EdgeMap<int> >
|
|
#endif
|
|
class MaxWeightedFractionalMatching {
|
|
public:
|
|
|
|
/// The graph type of the algorithm
|
|
typedef GR Graph;
|
|
/// The type of the edge weight map
|
|
typedef WM WeightMap;
|
|
/// The value type of the edge weights
|
|
typedef typename WeightMap::Value Value;
|
|
|
|
/// The type of the matching map
|
|
typedef typename Graph::template NodeMap<typename Graph::Arc>
|
|
MatchingMap;
|
|
|
|
/// \brief Scaling factor for primal solution
|
|
///
|
|
/// Scaling factor for primal solution.
|
|
static const int primalScale = 2;
|
|
|
|
/// \brief Scaling factor for dual solution
|
|
///
|
|
/// Scaling factor for dual solution. It is equal to 4 or 1
|
|
/// according to the value type.
|
|
static const int dualScale =
|
|
std::numeric_limits<Value>::is_integer ? 4 : 1;
|
|
|
|
private:
|
|
|
|
TEMPLATE_GRAPH_TYPEDEFS(Graph);
|
|
|
|
typedef typename Graph::template NodeMap<Value> NodePotential;
|
|
|
|
const Graph& _graph;
|
|
const WeightMap& _weight;
|
|
|
|
MatchingMap* _matching;
|
|
NodePotential* _node_potential;
|
|
|
|
int _node_num;
|
|
bool _allow_loops;
|
|
|
|
enum Status {
|
|
EVEN = -1, MATCHED = 0, ODD = 1
|
|
};
|
|
|
|
typedef typename Graph::template NodeMap<Status> StatusMap;
|
|
StatusMap* _status;
|
|
|
|
typedef typename Graph::template NodeMap<Arc> PredMap;
|
|
PredMap* _pred;
|
|
|
|
typedef ExtendFindEnum<IntNodeMap> TreeSet;
|
|
|
|
IntNodeMap *_tree_set_index;
|
|
TreeSet *_tree_set;
|
|
|
|
IntNodeMap *_delta1_index;
|
|
BinHeap<Value, IntNodeMap> *_delta1;
|
|
|
|
IntNodeMap *_delta2_index;
|
|
BinHeap<Value, IntNodeMap> *_delta2;
|
|
|
|
IntEdgeMap *_delta3_index;
|
|
BinHeap<Value, IntEdgeMap> *_delta3;
|
|
|
|
Value _delta_sum;
|
|
|
|
void createStructures() {
|
|
_node_num = countNodes(_graph);
|
|
|
|
if (!_matching) {
|
|
_matching = new MatchingMap(_graph);
|
|
}
|
|
if (!_node_potential) {
|
|
_node_potential = new NodePotential(_graph);
|
|
}
|
|
if (!_status) {
|
|
_status = new StatusMap(_graph);
|
|
}
|
|
if (!_pred) {
|
|
_pred = new PredMap(_graph);
|
|
}
|
|
if (!_tree_set) {
|
|
_tree_set_index = new IntNodeMap(_graph);
|
|
_tree_set = new TreeSet(*_tree_set_index);
|
|
}
|
|
if (!_delta1) {
|
|
_delta1_index = new IntNodeMap(_graph);
|
|
_delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
|
|
}
|
|
if (!_delta2) {
|
|
_delta2_index = new IntNodeMap(_graph);
|
|
_delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
|
|
}
|
|
if (!_delta3) {
|
|
_delta3_index = new IntEdgeMap(_graph);
|
|
_delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
|
|
}
|
|
}
|
|
|
|
void destroyStructures() {
|
|
if (_matching) {
|
|
delete _matching;
|
|
}
|
|
if (_node_potential) {
|
|
delete _node_potential;
|
|
}
|
|
if (_status) {
|
|
delete _status;
|
|
}
|
|
if (_pred) {
|
|
delete _pred;
|
|
}
|
|
if (_tree_set) {
|
|
delete _tree_set_index;
|
|
delete _tree_set;
|
|
}
|
|
if (_delta1) {
|
|
delete _delta1_index;
|
|
delete _delta1;
|
|
}
|
|
if (_delta2) {
|
|
delete _delta2_index;
|
|
delete _delta2;
|
|
}
|
|
if (_delta3) {
|
|
delete _delta3_index;
|
|
delete _delta3;
|
|
}
|
|
}
|
|
|
|
void matchedToEven(Node node, int tree) {
|
|
_tree_set->insert(node, tree);
|
|
_node_potential->set(node, (*_node_potential)[node] + _delta_sum);
|
|
_delta1->push(node, (*_node_potential)[node]);
|
|
|
|
if (_delta2->state(node) == _delta2->IN_HEAP) {
|
|
_delta2->erase(node);
|
|
}
|
|
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
if (node == v) {
|
|
if (_allow_loops && _graph.direction(a)) {
|
|
_delta3->push(a, rw / 2);
|
|
}
|
|
} else if ((*_status)[v] == EVEN) {
|
|
_delta3->push(a, rw / 2);
|
|
} else if ((*_status)[v] == MATCHED) {
|
|
if (_delta2->state(v) != _delta2->IN_HEAP) {
|
|
_pred->set(v, a);
|
|
_delta2->push(v, rw);
|
|
} else if ((*_delta2)[v] > rw) {
|
|
_pred->set(v, a);
|
|
_delta2->decrease(v, rw);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void matchedToOdd(Node node, int tree) {
|
|
_tree_set->insert(node, tree);
|
|
_node_potential->set(node, (*_node_potential)[node] - _delta_sum);
|
|
|
|
if (_delta2->state(node) == _delta2->IN_HEAP) {
|
|
_delta2->erase(node);
|
|
}
|
|
}
|
|
|
|
void evenToMatched(Node node, int tree) {
|
|
_delta1->erase(node);
|
|
_node_potential->set(node, (*_node_potential)[node] - _delta_sum);
|
|
Arc min = INVALID;
|
|
Value minrw = std::numeric_limits<Value>::max();
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
|
|
if (node == v) {
|
|
if (_allow_loops && _graph.direction(a)) {
|
|
_delta3->erase(a);
|
|
}
|
|
} else if ((*_status)[v] == EVEN) {
|
|
_delta3->erase(a);
|
|
if (minrw > rw) {
|
|
min = _graph.oppositeArc(a);
|
|
minrw = rw;
|
|
}
|
|
} else if ((*_status)[v] == MATCHED) {
|
|
if ((*_pred)[v] == a) {
|
|
Arc mina = INVALID;
|
|
Value minrwa = std::numeric_limits<Value>::max();
|
|
for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
|
|
Node va = _graph.target(aa);
|
|
if ((*_status)[va] != EVEN ||
|
|
_tree_set->find(va) == tree) continue;
|
|
Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
|
|
dualScale * _weight[aa];
|
|
if (minrwa > rwa) {
|
|
minrwa = rwa;
|
|
mina = aa;
|
|
}
|
|
}
|
|
if (mina != INVALID) {
|
|
_pred->set(v, mina);
|
|
_delta2->increase(v, minrwa);
|
|
} else {
|
|
_pred->set(v, INVALID);
|
|
_delta2->erase(v);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (min != INVALID) {
|
|
_pred->set(node, min);
|
|
_delta2->push(node, minrw);
|
|
} else {
|
|
_pred->set(node, INVALID);
|
|
}
|
|
}
|
|
|
|
void oddToMatched(Node node) {
|
|
_node_potential->set(node, (*_node_potential)[node] + _delta_sum);
|
|
Arc min = INVALID;
|
|
Value minrw = std::numeric_limits<Value>::max();
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
if ((*_status)[v] != EVEN) continue;
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
|
|
if (minrw > rw) {
|
|
min = _graph.oppositeArc(a);
|
|
minrw = rw;
|
|
}
|
|
}
|
|
if (min != INVALID) {
|
|
_pred->set(node, min);
|
|
_delta2->push(node, minrw);
|
|
} else {
|
|
_pred->set(node, INVALID);
|
|
}
|
|
}
|
|
|
|
void alternatePath(Node even, int tree) {
|
|
Node odd;
|
|
|
|
_status->set(even, MATCHED);
|
|
evenToMatched(even, tree);
|
|
|
|
Arc prev = (*_matching)[even];
|
|
while (prev != INVALID) {
|
|
odd = _graph.target(prev);
|
|
even = _graph.target((*_pred)[odd]);
|
|
_matching->set(odd, (*_pred)[odd]);
|
|
_status->set(odd, MATCHED);
|
|
oddToMatched(odd);
|
|
|
|
prev = (*_matching)[even];
|
|
_status->set(even, MATCHED);
|
|
_matching->set(even, _graph.oppositeArc((*_matching)[odd]));
|
|
evenToMatched(even, tree);
|
|
}
|
|
}
|
|
|
|
void destroyTree(int tree) {
|
|
for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
|
|
if ((*_status)[n] == EVEN) {
|
|
_status->set(n, MATCHED);
|
|
evenToMatched(n, tree);
|
|
} else if ((*_status)[n] == ODD) {
|
|
_status->set(n, MATCHED);
|
|
oddToMatched(n);
|
|
}
|
|
}
|
|
_tree_set->eraseClass(tree);
|
|
}
|
|
|
|
|
|
void unmatchNode(const Node& node) {
|
|
int tree = _tree_set->find(node);
|
|
|
|
alternatePath(node, tree);
|
|
destroyTree(tree);
|
|
|
|
_matching->set(node, INVALID);
|
|
}
|
|
|
|
|
|
void augmentOnEdge(const Edge& edge) {
|
|
Node left = _graph.u(edge);
|
|
int left_tree = _tree_set->find(left);
|
|
|
|
alternatePath(left, left_tree);
|
|
destroyTree(left_tree);
|
|
_matching->set(left, _graph.direct(edge, true));
|
|
|
|
Node right = _graph.v(edge);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.direct(edge, false));
|
|
}
|
|
|
|
void augmentOnArc(const Arc& arc) {
|
|
Node left = _graph.source(arc);
|
|
_status->set(left, MATCHED);
|
|
_matching->set(left, arc);
|
|
_pred->set(left, arc);
|
|
|
|
Node right = _graph.target(arc);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.oppositeArc(arc));
|
|
}
|
|
|
|
void extendOnArc(const Arc& arc) {
|
|
Node base = _graph.target(arc);
|
|
int tree = _tree_set->find(base);
|
|
|
|
Node odd = _graph.source(arc);
|
|
_tree_set->insert(odd, tree);
|
|
_status->set(odd, ODD);
|
|
matchedToOdd(odd, tree);
|
|
_pred->set(odd, arc);
|
|
|
|
Node even = _graph.target((*_matching)[odd]);
|
|
_tree_set->insert(even, tree);
|
|
_status->set(even, EVEN);
|
|
matchedToEven(even, tree);
|
|
}
|
|
|
|
void cycleOnEdge(const Edge& edge, int tree) {
|
|
Node nca = INVALID;
|
|
std::vector<Node> left_path, right_path;
|
|
|
|
{
|
|
std::set<Node> left_set, right_set;
|
|
Node left = _graph.u(edge);
|
|
left_path.push_back(left);
|
|
left_set.insert(left);
|
|
|
|
Node right = _graph.v(edge);
|
|
right_path.push_back(right);
|
|
right_set.insert(right);
|
|
|
|
while (true) {
|
|
|
|
if (left_set.find(right) != left_set.end()) {
|
|
nca = right;
|
|
break;
|
|
}
|
|
|
|
if ((*_matching)[left] == INVALID) break;
|
|
|
|
left = _graph.target((*_matching)[left]);
|
|
left_path.push_back(left);
|
|
left = _graph.target((*_pred)[left]);
|
|
left_path.push_back(left);
|
|
|
|
left_set.insert(left);
|
|
|
|
if (right_set.find(left) != right_set.end()) {
|
|
nca = left;
|
|
break;
|
|
}
|
|
|
|
if ((*_matching)[right] == INVALID) break;
|
|
|
|
right = _graph.target((*_matching)[right]);
|
|
right_path.push_back(right);
|
|
right = _graph.target((*_pred)[right]);
|
|
right_path.push_back(right);
|
|
|
|
right_set.insert(right);
|
|
|
|
}
|
|
|
|
if (nca == INVALID) {
|
|
if ((*_matching)[left] == INVALID) {
|
|
nca = right;
|
|
while (left_set.find(nca) == left_set.end()) {
|
|
nca = _graph.target((*_matching)[nca]);
|
|
right_path.push_back(nca);
|
|
nca = _graph.target((*_pred)[nca]);
|
|
right_path.push_back(nca);
|
|
}
|
|
} else {
|
|
nca = left;
|
|
while (right_set.find(nca) == right_set.end()) {
|
|
nca = _graph.target((*_matching)[nca]);
|
|
left_path.push_back(nca);
|
|
nca = _graph.target((*_pred)[nca]);
|
|
left_path.push_back(nca);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
alternatePath(nca, tree);
|
|
Arc prev;
|
|
|
|
prev = _graph.direct(edge, true);
|
|
for (int i = 0; left_path[i] != nca; i += 2) {
|
|
_matching->set(left_path[i], prev);
|
|
_status->set(left_path[i], MATCHED);
|
|
evenToMatched(left_path[i], tree);
|
|
|
|
prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
|
|
_status->set(left_path[i + 1], MATCHED);
|
|
oddToMatched(left_path[i + 1]);
|
|
}
|
|
_matching->set(nca, prev);
|
|
|
|
for (int i = 0; right_path[i] != nca; i += 2) {
|
|
_status->set(right_path[i], MATCHED);
|
|
evenToMatched(right_path[i], tree);
|
|
|
|
_matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
|
|
_status->set(right_path[i + 1], MATCHED);
|
|
oddToMatched(right_path[i + 1]);
|
|
}
|
|
|
|
destroyTree(tree);
|
|
}
|
|
|
|
void extractCycle(const Arc &arc) {
|
|
Node left = _graph.source(arc);
|
|
Node odd = _graph.target((*_matching)[left]);
|
|
Arc prev;
|
|
while (odd != left) {
|
|
Node even = _graph.target((*_matching)[odd]);
|
|
prev = (*_matching)[odd];
|
|
odd = _graph.target((*_matching)[even]);
|
|
_matching->set(even, _graph.oppositeArc(prev));
|
|
}
|
|
_matching->set(left, arc);
|
|
|
|
Node right = _graph.target(arc);
|
|
int right_tree = _tree_set->find(right);
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.oppositeArc(arc));
|
|
}
|
|
|
|
public:
|
|
|
|
/// \brief Constructor
|
|
///
|
|
/// Constructor.
|
|
MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
|
|
bool allow_loops = true)
|
|
: _graph(graph), _weight(weight), _matching(0),
|
|
_node_potential(0), _node_num(0), _allow_loops(allow_loops),
|
|
_status(0), _pred(0),
|
|
_tree_set_index(0), _tree_set(0),
|
|
|
|
_delta1_index(0), _delta1(0),
|
|
_delta2_index(0), _delta2(0),
|
|
_delta3_index(0), _delta3(0),
|
|
|
|
_delta_sum() {}
|
|
|
|
~MaxWeightedFractionalMatching() {
|
|
destroyStructures();
|
|
}
|
|
|
|
/// \name Execution Control
|
|
/// The simplest way to execute the algorithm is to use the
|
|
/// \ref run() member function.
|
|
|
|
///@{
|
|
|
|
/// \brief Initialize the algorithm
|
|
///
|
|
/// This function initializes the algorithm.
|
|
void init() {
|
|
createStructures();
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
(*_delta1_index)[n] = _delta1->PRE_HEAP;
|
|
(*_delta2_index)[n] = _delta2->PRE_HEAP;
|
|
}
|
|
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
|
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
|
}
|
|
|
|
_delta1->clear();
|
|
_delta2->clear();
|
|
_delta3->clear();
|
|
_tree_set->clear();
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
Value max = 0;
|
|
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
|
if (_graph.target(e) == n && !_allow_loops) continue;
|
|
if ((dualScale * _weight[e]) / 2 > max) {
|
|
max = (dualScale * _weight[e]) / 2;
|
|
}
|
|
}
|
|
_node_potential->set(n, max);
|
|
_delta1->push(n, max);
|
|
|
|
_tree_set->insert(n);
|
|
|
|
_matching->set(n, INVALID);
|
|
_status->set(n, EVEN);
|
|
}
|
|
|
|
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
|
Node left = _graph.u(e);
|
|
Node right = _graph.v(e);
|
|
if (left == right && !_allow_loops) continue;
|
|
_delta3->push(e, ((*_node_potential)[left] +
|
|
(*_node_potential)[right] -
|
|
dualScale * _weight[e]) / 2);
|
|
}
|
|
}
|
|
|
|
/// \brief Start the algorithm
|
|
///
|
|
/// This function starts the algorithm.
|
|
///
|
|
/// \pre \ref init() must be called before using this function.
|
|
void start() {
|
|
enum OpType {
|
|
D1, D2, D3
|
|
};
|
|
|
|
int unmatched = _node_num;
|
|
while (unmatched > 0) {
|
|
Value d1 = !_delta1->empty() ?
|
|
_delta1->prio() : std::numeric_limits<Value>::max();
|
|
|
|
Value d2 = !_delta2->empty() ?
|
|
_delta2->prio() : std::numeric_limits<Value>::max();
|
|
|
|
Value d3 = !_delta3->empty() ?
|
|
_delta3->prio() : std::numeric_limits<Value>::max();
|
|
|
|
_delta_sum = d3; OpType ot = D3;
|
|
if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
|
|
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
|
|
|
|
switch (ot) {
|
|
case D1:
|
|
{
|
|
Node n = _delta1->top();
|
|
unmatchNode(n);
|
|
--unmatched;
|
|
}
|
|
break;
|
|
case D2:
|
|
{
|
|
Node n = _delta2->top();
|
|
Arc a = (*_pred)[n];
|
|
if ((*_matching)[n] == INVALID) {
|
|
augmentOnArc(a);
|
|
--unmatched;
|
|
} else {
|
|
Node v = _graph.target((*_matching)[n]);
|
|
if ((*_matching)[n] !=
|
|
_graph.oppositeArc((*_matching)[v])) {
|
|
extractCycle(a);
|
|
--unmatched;
|
|
} else {
|
|
extendOnArc(a);
|
|
}
|
|
}
|
|
} break;
|
|
case D3:
|
|
{
|
|
Edge e = _delta3->top();
|
|
|
|
Node left = _graph.u(e);
|
|
Node right = _graph.v(e);
|
|
|
|
int left_tree = _tree_set->find(left);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
if (left_tree == right_tree) {
|
|
cycleOnEdge(e, left_tree);
|
|
--unmatched;
|
|
} else {
|
|
augmentOnEdge(e);
|
|
unmatched -= 2;
|
|
}
|
|
} break;
|
|
}
|
|
}
|
|
}
|
|
|
|
/// \brief Run the algorithm.
|
|
///
|
|
/// This method runs the \c %MaxWeightedFractionalMatching algorithm.
|
|
///
|
|
/// \note mwfm.run() is just a shortcut of the following code.
|
|
/// \code
|
|
/// mwfm.init();
|
|
/// mwfm.start();
|
|
/// \endcode
|
|
void run() {
|
|
init();
|
|
start();
|
|
}
|
|
|
|
/// @}
|
|
|
|
/// \name Primal Solution
|
|
/// Functions to get the primal solution, i.e. the maximum weighted
|
|
/// matching.\n
|
|
/// Either \ref run() or \ref start() function should be called before
|
|
/// using them.
|
|
|
|
/// @{
|
|
|
|
/// \brief Return the weight of the matching.
|
|
///
|
|
/// This function returns the weight of the found matching. This
|
|
/// value is scaled by \ref primalScale "primal scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value matchingWeight() const {
|
|
Value sum = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] != INVALID) {
|
|
sum += _weight[(*_matching)[n]];
|
|
}
|
|
}
|
|
return sum * primalScale / 2;
|
|
}
|
|
|
|
/// \brief Return the number of covered nodes in the matching.
|
|
///
|
|
/// This function returns the number of covered nodes in the matching.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matchingSize() const {
|
|
int num = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] != INVALID) {
|
|
++num;
|
|
}
|
|
}
|
|
return num;
|
|
}
|
|
|
|
/// \brief Return \c true if the given edge is in the matching.
|
|
///
|
|
/// This function returns \c true if the given edge is in the
|
|
/// found matching. The result is scaled by \ref primalScale
|
|
/// "primal scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matching(const Edge& edge) const {
|
|
return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
|
|
+ (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
|
|
}
|
|
|
|
/// \brief Return the fractional matching arc (or edge) incident
|
|
/// to the given node.
|
|
///
|
|
/// This function returns one of the fractional matching arc (or
|
|
/// edge) incident to the given node in the found matching or \c
|
|
/// INVALID if the node is not covered by the matching or if the
|
|
/// node is on an odd length cycle then it is the successor edge
|
|
/// on the cycle.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Arc matching(const Node& node) const {
|
|
return (*_matching)[node];
|
|
}
|
|
|
|
/// \brief Return a const reference to the matching map.
|
|
///
|
|
/// This function returns a const reference to a node map that stores
|
|
/// the matching arc (or edge) incident to each node.
|
|
const MatchingMap& matchingMap() const {
|
|
return *_matching;
|
|
}
|
|
|
|
/// @}
|
|
|
|
/// \name Dual Solution
|
|
/// Functions to get the dual solution.\n
|
|
/// Either \ref run() or \ref start() function should be called before
|
|
/// using them.
|
|
|
|
/// @{
|
|
|
|
/// \brief Return the value of the dual solution.
|
|
///
|
|
/// This function returns the value of the dual solution.
|
|
/// It should be equal to the primal value scaled by \ref dualScale
|
|
/// "dual scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value dualValue() const {
|
|
Value sum = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
sum += nodeValue(n);
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
/// \brief Return the dual value (potential) of the given node.
|
|
///
|
|
/// This function returns the dual value (potential) of the given node.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value nodeValue(const Node& n) const {
|
|
return (*_node_potential)[n];
|
|
}
|
|
|
|
/// @}
|
|
|
|
};
|
|
|
|
/// \ingroup matching
|
|
///
|
|
/// \brief Weighted fractional perfect matching in general graphs
|
|
///
|
|
/// This class provides an efficient implementation of fractional
|
|
/// matching algorithm. The implementation uses priority queues and
|
|
/// provides \f$O(nm\log n)\f$ time complexity.
|
|
///
|
|
/// The maximum weighted fractional perfect matching is a relaxation
|
|
/// of the maximum weighted perfect matching problem where the odd
|
|
/// set constraints are omitted.
|
|
/// It can be formulated with the following linear program.
|
|
/// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
|
|
/// \f[x_e \ge 0\quad \forall e\in E\f]
|
|
/// \f[\max \sum_{e\in E}x_ew_e\f]
|
|
/// where \f$\delta(X)\f$ is the set of edges incident to a node in
|
|
/// \f$X\f$. The result must be the union of a matching with one
|
|
/// value edges and a set of odd length cycles with half value edges.
|
|
///
|
|
/// The algorithm calculates an optimal fractional matching and a
|
|
/// proof of the optimality. The solution of the dual problem can be
|
|
/// used to check the result of the algorithm. The dual linear
|
|
/// problem is the following.
|
|
/// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
|
|
/// \f[\min \sum_{u \in V}y_u \f]
|
|
///
|
|
/// The algorithm can be executed with the run() function.
|
|
/// After it the matching (the primal solution) and the dual solution
|
|
/// can be obtained using the query functions.
|
|
///
|
|
/// The primal solution is multiplied by
|
|
/// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
|
|
/// If the value type is integer, then the dual
|
|
/// solution is scaled by
|
|
/// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
|
|
///
|
|
/// \tparam GR The undirected graph type the algorithm runs on.
|
|
/// \tparam WM The type edge weight map. The default type is
|
|
/// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
|
|
#ifdef DOXYGEN
|
|
template <typename GR, typename WM>
|
|
#else
|
|
template <typename GR,
|
|
typename WM = typename GR::template EdgeMap<int> >
|
|
#endif
|
|
class MaxWeightedPerfectFractionalMatching {
|
|
public:
|
|
|
|
/// The graph type of the algorithm
|
|
typedef GR Graph;
|
|
/// The type of the edge weight map
|
|
typedef WM WeightMap;
|
|
/// The value type of the edge weights
|
|
typedef typename WeightMap::Value Value;
|
|
|
|
/// The type of the matching map
|
|
typedef typename Graph::template NodeMap<typename Graph::Arc>
|
|
MatchingMap;
|
|
|
|
/// \brief Scaling factor for primal solution
|
|
///
|
|
/// Scaling factor for primal solution.
|
|
static const int primalScale = 2;
|
|
|
|
/// \brief Scaling factor for dual solution
|
|
///
|
|
/// Scaling factor for dual solution. It is equal to 4 or 1
|
|
/// according to the value type.
|
|
static const int dualScale =
|
|
std::numeric_limits<Value>::is_integer ? 4 : 1;
|
|
|
|
private:
|
|
|
|
TEMPLATE_GRAPH_TYPEDEFS(Graph);
|
|
|
|
typedef typename Graph::template NodeMap<Value> NodePotential;
|
|
|
|
const Graph& _graph;
|
|
const WeightMap& _weight;
|
|
|
|
MatchingMap* _matching;
|
|
NodePotential* _node_potential;
|
|
|
|
int _node_num;
|
|
bool _allow_loops;
|
|
|
|
enum Status {
|
|
EVEN = -1, MATCHED = 0, ODD = 1
|
|
};
|
|
|
|
typedef typename Graph::template NodeMap<Status> StatusMap;
|
|
StatusMap* _status;
|
|
|
|
typedef typename Graph::template NodeMap<Arc> PredMap;
|
|
PredMap* _pred;
|
|
|
|
typedef ExtendFindEnum<IntNodeMap> TreeSet;
|
|
|
|
IntNodeMap *_tree_set_index;
|
|
TreeSet *_tree_set;
|
|
|
|
IntNodeMap *_delta2_index;
|
|
BinHeap<Value, IntNodeMap> *_delta2;
|
|
|
|
IntEdgeMap *_delta3_index;
|
|
BinHeap<Value, IntEdgeMap> *_delta3;
|
|
|
|
Value _delta_sum;
|
|
|
|
void createStructures() {
|
|
_node_num = countNodes(_graph);
|
|
|
|
if (!_matching) {
|
|
_matching = new MatchingMap(_graph);
|
|
}
|
|
if (!_node_potential) {
|
|
_node_potential = new NodePotential(_graph);
|
|
}
|
|
if (!_status) {
|
|
_status = new StatusMap(_graph);
|
|
}
|
|
if (!_pred) {
|
|
_pred = new PredMap(_graph);
|
|
}
|
|
if (!_tree_set) {
|
|
_tree_set_index = new IntNodeMap(_graph);
|
|
_tree_set = new TreeSet(*_tree_set_index);
|
|
}
|
|
if (!_delta2) {
|
|
_delta2_index = new IntNodeMap(_graph);
|
|
_delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
|
|
}
|
|
if (!_delta3) {
|
|
_delta3_index = new IntEdgeMap(_graph);
|
|
_delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
|
|
}
|
|
}
|
|
|
|
void destroyStructures() {
|
|
if (_matching) {
|
|
delete _matching;
|
|
}
|
|
if (_node_potential) {
|
|
delete _node_potential;
|
|
}
|
|
if (_status) {
|
|
delete _status;
|
|
}
|
|
if (_pred) {
|
|
delete _pred;
|
|
}
|
|
if (_tree_set) {
|
|
delete _tree_set_index;
|
|
delete _tree_set;
|
|
}
|
|
if (_delta2) {
|
|
delete _delta2_index;
|
|
delete _delta2;
|
|
}
|
|
if (_delta3) {
|
|
delete _delta3_index;
|
|
delete _delta3;
|
|
}
|
|
}
|
|
|
|
void matchedToEven(Node node, int tree) {
|
|
_tree_set->insert(node, tree);
|
|
_node_potential->set(node, (*_node_potential)[node] + _delta_sum);
|
|
|
|
if (_delta2->state(node) == _delta2->IN_HEAP) {
|
|
_delta2->erase(node);
|
|
}
|
|
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
if (node == v) {
|
|
if (_allow_loops && _graph.direction(a)) {
|
|
_delta3->push(a, rw / 2);
|
|
}
|
|
} else if ((*_status)[v] == EVEN) {
|
|
_delta3->push(a, rw / 2);
|
|
} else if ((*_status)[v] == MATCHED) {
|
|
if (_delta2->state(v) != _delta2->IN_HEAP) {
|
|
_pred->set(v, a);
|
|
_delta2->push(v, rw);
|
|
} else if ((*_delta2)[v] > rw) {
|
|
_pred->set(v, a);
|
|
_delta2->decrease(v, rw);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void matchedToOdd(Node node, int tree) {
|
|
_tree_set->insert(node, tree);
|
|
_node_potential->set(node, (*_node_potential)[node] - _delta_sum);
|
|
|
|
if (_delta2->state(node) == _delta2->IN_HEAP) {
|
|
_delta2->erase(node);
|
|
}
|
|
}
|
|
|
|
void evenToMatched(Node node, int tree) {
|
|
_node_potential->set(node, (*_node_potential)[node] - _delta_sum);
|
|
Arc min = INVALID;
|
|
Value minrw = std::numeric_limits<Value>::max();
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
|
|
if (node == v) {
|
|
if (_allow_loops && _graph.direction(a)) {
|
|
_delta3->erase(a);
|
|
}
|
|
} else if ((*_status)[v] == EVEN) {
|
|
_delta3->erase(a);
|
|
if (minrw > rw) {
|
|
min = _graph.oppositeArc(a);
|
|
minrw = rw;
|
|
}
|
|
} else if ((*_status)[v] == MATCHED) {
|
|
if ((*_pred)[v] == a) {
|
|
Arc mina = INVALID;
|
|
Value minrwa = std::numeric_limits<Value>::max();
|
|
for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
|
|
Node va = _graph.target(aa);
|
|
if ((*_status)[va] != EVEN ||
|
|
_tree_set->find(va) == tree) continue;
|
|
Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
|
|
dualScale * _weight[aa];
|
|
if (minrwa > rwa) {
|
|
minrwa = rwa;
|
|
mina = aa;
|
|
}
|
|
}
|
|
if (mina != INVALID) {
|
|
_pred->set(v, mina);
|
|
_delta2->increase(v, minrwa);
|
|
} else {
|
|
_pred->set(v, INVALID);
|
|
_delta2->erase(v);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (min != INVALID) {
|
|
_pred->set(node, min);
|
|
_delta2->push(node, minrw);
|
|
} else {
|
|
_pred->set(node, INVALID);
|
|
}
|
|
}
|
|
|
|
void oddToMatched(Node node) {
|
|
_node_potential->set(node, (*_node_potential)[node] + _delta_sum);
|
|
Arc min = INVALID;
|
|
Value minrw = std::numeric_limits<Value>::max();
|
|
for (InArcIt a(_graph, node); a != INVALID; ++a) {
|
|
Node v = _graph.source(a);
|
|
if ((*_status)[v] != EVEN) continue;
|
|
Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
|
|
dualScale * _weight[a];
|
|
|
|
if (minrw > rw) {
|
|
min = _graph.oppositeArc(a);
|
|
minrw = rw;
|
|
}
|
|
}
|
|
if (min != INVALID) {
|
|
_pred->set(node, min);
|
|
_delta2->push(node, minrw);
|
|
} else {
|
|
_pred->set(node, INVALID);
|
|
}
|
|
}
|
|
|
|
void alternatePath(Node even, int tree) {
|
|
Node odd;
|
|
|
|
_status->set(even, MATCHED);
|
|
evenToMatched(even, tree);
|
|
|
|
Arc prev = (*_matching)[even];
|
|
while (prev != INVALID) {
|
|
odd = _graph.target(prev);
|
|
even = _graph.target((*_pred)[odd]);
|
|
_matching->set(odd, (*_pred)[odd]);
|
|
_status->set(odd, MATCHED);
|
|
oddToMatched(odd);
|
|
|
|
prev = (*_matching)[even];
|
|
_status->set(even, MATCHED);
|
|
_matching->set(even, _graph.oppositeArc((*_matching)[odd]));
|
|
evenToMatched(even, tree);
|
|
}
|
|
}
|
|
|
|
void destroyTree(int tree) {
|
|
for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
|
|
if ((*_status)[n] == EVEN) {
|
|
_status->set(n, MATCHED);
|
|
evenToMatched(n, tree);
|
|
} else if ((*_status)[n] == ODD) {
|
|
_status->set(n, MATCHED);
|
|
oddToMatched(n);
|
|
}
|
|
}
|
|
_tree_set->eraseClass(tree);
|
|
}
|
|
|
|
void augmentOnEdge(const Edge& edge) {
|
|
Node left = _graph.u(edge);
|
|
int left_tree = _tree_set->find(left);
|
|
|
|
alternatePath(left, left_tree);
|
|
destroyTree(left_tree);
|
|
_matching->set(left, _graph.direct(edge, true));
|
|
|
|
Node right = _graph.v(edge);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.direct(edge, false));
|
|
}
|
|
|
|
void augmentOnArc(const Arc& arc) {
|
|
Node left = _graph.source(arc);
|
|
_status->set(left, MATCHED);
|
|
_matching->set(left, arc);
|
|
_pred->set(left, arc);
|
|
|
|
Node right = _graph.target(arc);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.oppositeArc(arc));
|
|
}
|
|
|
|
void extendOnArc(const Arc& arc) {
|
|
Node base = _graph.target(arc);
|
|
int tree = _tree_set->find(base);
|
|
|
|
Node odd = _graph.source(arc);
|
|
_tree_set->insert(odd, tree);
|
|
_status->set(odd, ODD);
|
|
matchedToOdd(odd, tree);
|
|
_pred->set(odd, arc);
|
|
|
|
Node even = _graph.target((*_matching)[odd]);
|
|
_tree_set->insert(even, tree);
|
|
_status->set(even, EVEN);
|
|
matchedToEven(even, tree);
|
|
}
|
|
|
|
void cycleOnEdge(const Edge& edge, int tree) {
|
|
Node nca = INVALID;
|
|
std::vector<Node> left_path, right_path;
|
|
|
|
{
|
|
std::set<Node> left_set, right_set;
|
|
Node left = _graph.u(edge);
|
|
left_path.push_back(left);
|
|
left_set.insert(left);
|
|
|
|
Node right = _graph.v(edge);
|
|
right_path.push_back(right);
|
|
right_set.insert(right);
|
|
|
|
while (true) {
|
|
|
|
if (left_set.find(right) != left_set.end()) {
|
|
nca = right;
|
|
break;
|
|
}
|
|
|
|
if ((*_matching)[left] == INVALID) break;
|
|
|
|
left = _graph.target((*_matching)[left]);
|
|
left_path.push_back(left);
|
|
left = _graph.target((*_pred)[left]);
|
|
left_path.push_back(left);
|
|
|
|
left_set.insert(left);
|
|
|
|
if (right_set.find(left) != right_set.end()) {
|
|
nca = left;
|
|
break;
|
|
}
|
|
|
|
if ((*_matching)[right] == INVALID) break;
|
|
|
|
right = _graph.target((*_matching)[right]);
|
|
right_path.push_back(right);
|
|
right = _graph.target((*_pred)[right]);
|
|
right_path.push_back(right);
|
|
|
|
right_set.insert(right);
|
|
|
|
}
|
|
|
|
if (nca == INVALID) {
|
|
if ((*_matching)[left] == INVALID) {
|
|
nca = right;
|
|
while (left_set.find(nca) == left_set.end()) {
|
|
nca = _graph.target((*_matching)[nca]);
|
|
right_path.push_back(nca);
|
|
nca = _graph.target((*_pred)[nca]);
|
|
right_path.push_back(nca);
|
|
}
|
|
} else {
|
|
nca = left;
|
|
while (right_set.find(nca) == right_set.end()) {
|
|
nca = _graph.target((*_matching)[nca]);
|
|
left_path.push_back(nca);
|
|
nca = _graph.target((*_pred)[nca]);
|
|
left_path.push_back(nca);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
alternatePath(nca, tree);
|
|
Arc prev;
|
|
|
|
prev = _graph.direct(edge, true);
|
|
for (int i = 0; left_path[i] != nca; i += 2) {
|
|
_matching->set(left_path[i], prev);
|
|
_status->set(left_path[i], MATCHED);
|
|
evenToMatched(left_path[i], tree);
|
|
|
|
prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
|
|
_status->set(left_path[i + 1], MATCHED);
|
|
oddToMatched(left_path[i + 1]);
|
|
}
|
|
_matching->set(nca, prev);
|
|
|
|
for (int i = 0; right_path[i] != nca; i += 2) {
|
|
_status->set(right_path[i], MATCHED);
|
|
evenToMatched(right_path[i], tree);
|
|
|
|
_matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
|
|
_status->set(right_path[i + 1], MATCHED);
|
|
oddToMatched(right_path[i + 1]);
|
|
}
|
|
|
|
destroyTree(tree);
|
|
}
|
|
|
|
void extractCycle(const Arc &arc) {
|
|
Node left = _graph.source(arc);
|
|
Node odd = _graph.target((*_matching)[left]);
|
|
Arc prev;
|
|
while (odd != left) {
|
|
Node even = _graph.target((*_matching)[odd]);
|
|
prev = (*_matching)[odd];
|
|
odd = _graph.target((*_matching)[even]);
|
|
_matching->set(even, _graph.oppositeArc(prev));
|
|
}
|
|
_matching->set(left, arc);
|
|
|
|
Node right = _graph.target(arc);
|
|
int right_tree = _tree_set->find(right);
|
|
alternatePath(right, right_tree);
|
|
destroyTree(right_tree);
|
|
_matching->set(right, _graph.oppositeArc(arc));
|
|
}
|
|
|
|
public:
|
|
|
|
/// \brief Constructor
|
|
///
|
|
/// Constructor.
|
|
MaxWeightedPerfectFractionalMatching(const Graph& graph,
|
|
const WeightMap& weight,
|
|
bool allow_loops = true)
|
|
: _graph(graph), _weight(weight), _matching(0),
|
|
_node_potential(0), _node_num(0), _allow_loops(allow_loops),
|
|
_status(0), _pred(0),
|
|
_tree_set_index(0), _tree_set(0),
|
|
|
|
_delta2_index(0), _delta2(0),
|
|
_delta3_index(0), _delta3(0),
|
|
|
|
_delta_sum() {}
|
|
|
|
~MaxWeightedPerfectFractionalMatching() {
|
|
destroyStructures();
|
|
}
|
|
|
|
/// \name Execution Control
|
|
/// The simplest way to execute the algorithm is to use the
|
|
/// \ref run() member function.
|
|
|
|
///@{
|
|
|
|
/// \brief Initialize the algorithm
|
|
///
|
|
/// This function initializes the algorithm.
|
|
void init() {
|
|
createStructures();
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
(*_delta2_index)[n] = _delta2->PRE_HEAP;
|
|
}
|
|
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
|
(*_delta3_index)[e] = _delta3->PRE_HEAP;
|
|
}
|
|
|
|
_delta2->clear();
|
|
_delta3->clear();
|
|
_tree_set->clear();
|
|
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
Value max = - std::numeric_limits<Value>::max();
|
|
for (OutArcIt e(_graph, n); e != INVALID; ++e) {
|
|
if (_graph.target(e) == n && !_allow_loops) continue;
|
|
if ((dualScale * _weight[e]) / 2 > max) {
|
|
max = (dualScale * _weight[e]) / 2;
|
|
}
|
|
}
|
|
_node_potential->set(n, max);
|
|
|
|
_tree_set->insert(n);
|
|
|
|
_matching->set(n, INVALID);
|
|
_status->set(n, EVEN);
|
|
}
|
|
|
|
for (EdgeIt e(_graph); e != INVALID; ++e) {
|
|
Node left = _graph.u(e);
|
|
Node right = _graph.v(e);
|
|
if (left == right && !_allow_loops) continue;
|
|
_delta3->push(e, ((*_node_potential)[left] +
|
|
(*_node_potential)[right] -
|
|
dualScale * _weight[e]) / 2);
|
|
}
|
|
}
|
|
|
|
/// \brief Start the algorithm
|
|
///
|
|
/// This function starts the algorithm.
|
|
///
|
|
/// \pre \ref init() must be called before using this function.
|
|
bool start() {
|
|
enum OpType {
|
|
D2, D3
|
|
};
|
|
|
|
int unmatched = _node_num;
|
|
while (unmatched > 0) {
|
|
Value d2 = !_delta2->empty() ?
|
|
_delta2->prio() : std::numeric_limits<Value>::max();
|
|
|
|
Value d3 = !_delta3->empty() ?
|
|
_delta3->prio() : std::numeric_limits<Value>::max();
|
|
|
|
_delta_sum = d3; OpType ot = D3;
|
|
if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
|
|
|
|
if (_delta_sum == std::numeric_limits<Value>::max()) {
|
|
return false;
|
|
}
|
|
|
|
switch (ot) {
|
|
case D2:
|
|
{
|
|
Node n = _delta2->top();
|
|
Arc a = (*_pred)[n];
|
|
if ((*_matching)[n] == INVALID) {
|
|
augmentOnArc(a);
|
|
--unmatched;
|
|
} else {
|
|
Node v = _graph.target((*_matching)[n]);
|
|
if ((*_matching)[n] !=
|
|
_graph.oppositeArc((*_matching)[v])) {
|
|
extractCycle(a);
|
|
--unmatched;
|
|
} else {
|
|
extendOnArc(a);
|
|
}
|
|
}
|
|
} break;
|
|
case D3:
|
|
{
|
|
Edge e = _delta3->top();
|
|
|
|
Node left = _graph.u(e);
|
|
Node right = _graph.v(e);
|
|
|
|
int left_tree = _tree_set->find(left);
|
|
int right_tree = _tree_set->find(right);
|
|
|
|
if (left_tree == right_tree) {
|
|
cycleOnEdge(e, left_tree);
|
|
--unmatched;
|
|
} else {
|
|
augmentOnEdge(e);
|
|
unmatched -= 2;
|
|
}
|
|
} break;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
/// \brief Run the algorithm.
|
|
///
|
|
/// This method runs the \c %MaxWeightedPerfectFractionalMatching
|
|
/// algorithm.
|
|
///
|
|
/// \note mwfm.run() is just a shortcut of the following code.
|
|
/// \code
|
|
/// mwpfm.init();
|
|
/// mwpfm.start();
|
|
/// \endcode
|
|
bool run() {
|
|
init();
|
|
return start();
|
|
}
|
|
|
|
/// @}
|
|
|
|
/// \name Primal Solution
|
|
/// Functions to get the primal solution, i.e. the maximum weighted
|
|
/// matching.\n
|
|
/// Either \ref run() or \ref start() function should be called before
|
|
/// using them.
|
|
|
|
/// @{
|
|
|
|
/// \brief Return the weight of the matching.
|
|
///
|
|
/// This function returns the weight of the found matching. This
|
|
/// value is scaled by \ref primalScale "primal scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value matchingWeight() const {
|
|
Value sum = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] != INVALID) {
|
|
sum += _weight[(*_matching)[n]];
|
|
}
|
|
}
|
|
return sum * primalScale / 2;
|
|
}
|
|
|
|
/// \brief Return the number of covered nodes in the matching.
|
|
///
|
|
/// This function returns the number of covered nodes in the matching.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matchingSize() const {
|
|
int num = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
if ((*_matching)[n] != INVALID) {
|
|
++num;
|
|
}
|
|
}
|
|
return num;
|
|
}
|
|
|
|
/// \brief Return \c true if the given edge is in the matching.
|
|
///
|
|
/// This function returns \c true if the given edge is in the
|
|
/// found matching. The result is scaled by \ref primalScale
|
|
/// "primal scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
int matching(const Edge& edge) const {
|
|
return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
|
|
+ (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
|
|
}
|
|
|
|
/// \brief Return the fractional matching arc (or edge) incident
|
|
/// to the given node.
|
|
///
|
|
/// This function returns one of the fractional matching arc (or
|
|
/// edge) incident to the given node in the found matching or \c
|
|
/// INVALID if the node is not covered by the matching or if the
|
|
/// node is on an odd length cycle then it is the successor edge
|
|
/// on the cycle.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Arc matching(const Node& node) const {
|
|
return (*_matching)[node];
|
|
}
|
|
|
|
/// \brief Return a const reference to the matching map.
|
|
///
|
|
/// This function returns a const reference to a node map that stores
|
|
/// the matching arc (or edge) incident to each node.
|
|
const MatchingMap& matchingMap() const {
|
|
return *_matching;
|
|
}
|
|
|
|
/// @}
|
|
|
|
/// \name Dual Solution
|
|
/// Functions to get the dual solution.\n
|
|
/// Either \ref run() or \ref start() function should be called before
|
|
/// using them.
|
|
|
|
/// @{
|
|
|
|
/// \brief Return the value of the dual solution.
|
|
///
|
|
/// This function returns the value of the dual solution.
|
|
/// It should be equal to the primal value scaled by \ref dualScale
|
|
/// "dual scale".
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value dualValue() const {
|
|
Value sum = 0;
|
|
for (NodeIt n(_graph); n != INVALID; ++n) {
|
|
sum += nodeValue(n);
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
/// \brief Return the dual value (potential) of the given node.
|
|
///
|
|
/// This function returns the dual value (potential) of the given node.
|
|
///
|
|
/// \pre Either run() or start() must be called before using this function.
|
|
Value nodeValue(const Node& n) const {
|
|
return (*_node_potential)[n];
|
|
}
|
|
|
|
/// @}
|
|
|
|
};
|
|
|
|
} //END OF NAMESPACE LEMON
|
|
|
|
#endif //LEMON_FRACTIONAL_MATCHING_H
|