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