COIN-OR::LEMON - Graph Library

source: lemon/lemon/fractional_matching.h @ 949:61120524af27

Last change on this file since 949:61120524af27 was 948:636dadefe1e6, checked in by Balazs Dezso <deba@…>, 15 years ago

Add fractional matching algorithms (#314)

File size: 63.2 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2009
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_FRACTIONAL_MATCHING_H
20#define LEMON_FRACTIONAL_MATCHING_H
21
22#include <vector>
23#include <queue>
24#include <set>
25#include <limits>
26
27#include <lemon/core.h>
28#include <lemon/unionfind.h>
29#include <lemon/bin_heap.h>
30#include <lemon/maps.h>
31#include <lemon/assert.h>
32#include <lemon/elevator.h>
33
34///\ingroup matching
35///\file
36///\brief Fractional matching algorithms in general graphs.
37
38namespace lemon {
39
40  /// \brief Default traits class of MaxFractionalMatching class.
41  ///
42  /// Default traits class of MaxFractionalMatching class.
43  /// \tparam GR Graph type.
44  template <typename GR>
45  struct MaxFractionalMatchingDefaultTraits {
46
47    /// \brief The type of the graph the algorithm runs on.
48    typedef GR Graph;
49
50    /// \brief The type of the map that stores the matching.
51    ///
52    /// The type of the map that stores the matching arcs.
53    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
54    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
55
56    /// \brief Instantiates a MatchingMap.
57    ///
58    /// This function instantiates a \ref MatchingMap.
59    /// \param graph The graph for which we would like to define
60    /// the matching map.
61    static MatchingMap* createMatchingMap(const Graph& graph) {
62      return new MatchingMap(graph);
63    }
64
65    /// \brief The elevator type used by MaxFractionalMatching algorithm.
66    ///
67    /// The elevator type used by MaxFractionalMatching algorithm.
68    ///
69    /// \sa Elevator
70    /// \sa LinkedElevator
71    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
72
73    /// \brief Instantiates an Elevator.
74    ///
75    /// This function instantiates an \ref Elevator.
76    /// \param graph The graph for which we would like to define
77    /// the elevator.
78    /// \param max_level The maximum level of the elevator.
79    static Elevator* createElevator(const Graph& graph, int max_level) {
80      return new Elevator(graph, max_level);
81    }
82  };
83
84  /// \ingroup matching
85  ///
86  /// \brief Max cardinality fractional matching
87  ///
88  /// This class provides an implementation of fractional matching
89  /// algorithm based on push-relabel principle.
90  ///
91  /// The maximum cardinality fractional matching is a relaxation of the
92  /// maximum cardinality matching problem where the odd set constraints
93  /// are omitted.
94  /// It can be formulated with the following linear program.
95  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
96  /// \f[x_e \ge 0\quad \forall e\in E\f]
97  /// \f[\max \sum_{e\in E}x_e\f]
98  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
99  /// \f$X\f$. The result can be represented as the union of a
100  /// matching with one value edges and a set of odd length cycles
101  /// with half value edges.
102  ///
103  /// The algorithm calculates an optimal fractional matching and a
104  /// barrier. The number of adjacents of any node set minus the size
105  /// of node set is a lower bound on the uncovered nodes in the
106  /// graph. For maximum matching a barrier is computed which
107  /// maximizes this difference.
108  ///
109  /// The algorithm can be executed with the run() function.  After it
110  /// the matching (the primal solution) and the barrier (the dual
111  /// solution) can be obtained using the query functions.
112  ///
113  /// The primal solution is multiplied by
114  /// \ref MaxWeightedMatching::primalScale "2".
115  ///
116  /// \tparam GR The undirected graph type the algorithm runs on.
117#ifdef DOXYGEN
118  template <typename GR, typename TR>
119#else
120  template <typename GR,
121            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
122#endif
123  class MaxFractionalMatching {
124  public:
125
126    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
127    /// class" of the algorithm.
128    typedef TR Traits;
129    /// The type of the graph the algorithm runs on.
130    typedef typename TR::Graph Graph;
131    /// The type of the matching map.
132    typedef typename TR::MatchingMap MatchingMap;
133    /// The type of the elevator.
134    typedef typename TR::Elevator Elevator;
135
136    /// \brief Scaling factor for primal solution
137    ///
138    /// Scaling factor for primal solution.
139    static const int primalScale = 2;
140
141  private:
142
143    const Graph &_graph;
144    int _node_num;
145    bool _allow_loops;
146    int _empty_level;
147
148    TEMPLATE_GRAPH_TYPEDEFS(Graph);
149
150    bool _local_matching;
151    MatchingMap *_matching;
152
153    bool _local_level;
154    Elevator *_level;
155
156    typedef typename Graph::template NodeMap<int> InDegMap;
157    InDegMap *_indeg;
158
159    void createStructures() {
160      _node_num = countNodes(_graph);
161
162      if (!_matching) {
163        _local_matching = true;
164        _matching = Traits::createMatchingMap(_graph);
165      }
166      if (!_level) {
167        _local_level = true;
168        _level = Traits::createElevator(_graph, _node_num);
169      }
170      if (!_indeg) {
171        _indeg = new InDegMap(_graph);
172      }
173    }
174
175    void destroyStructures() {
176      if (_local_matching) {
177        delete _matching;
178      }
179      if (_local_level) {
180        delete _level;
181      }
182      if (_indeg) {
183        delete _indeg;
184      }
185    }
186
187    void postprocessing() {
188      for (NodeIt n(_graph); n != INVALID; ++n) {
189        if ((*_indeg)[n] != 0) continue;
190        _indeg->set(n, -1);
191        Node u = n;
192        while ((*_matching)[u] != INVALID) {
193          Node v = _graph.target((*_matching)[u]);
194          _indeg->set(v, -1);
195          Arc a = _graph.oppositeArc((*_matching)[u]);
196          u = _graph.target((*_matching)[v]);
197          _indeg->set(u, -1);
198          _matching->set(v, a);
199        }
200      }
201
202      for (NodeIt n(_graph); n != INVALID; ++n) {
203        if ((*_indeg)[n] != 1) continue;
204        _indeg->set(n, -1);
205
206        int num = 1;
207        Node u = _graph.target((*_matching)[n]);
208        while (u != n) {
209          _indeg->set(u, -1);
210          u = _graph.target((*_matching)[u]);
211          ++num;
212        }
213        if (num % 2 == 0 && num > 2) {
214          Arc prev = _graph.oppositeArc((*_matching)[n]);
215          Node v = _graph.target((*_matching)[n]);
216          u = _graph.target((*_matching)[v]);
217          _matching->set(v, prev);
218          while (u != n) {
219            prev = _graph.oppositeArc((*_matching)[u]);
220            v = _graph.target((*_matching)[u]);
221            u = _graph.target((*_matching)[v]);
222            _matching->set(v, prev);
223          }
224        }
225      }
226    }
227
228  public:
229
230    typedef MaxFractionalMatching Create;
231
232    ///\name Named Template Parameters
233
234    ///@{
235
236    template <typename T>
237    struct SetMatchingMapTraits : public Traits {
238      typedef T MatchingMap;
239      static MatchingMap *createMatchingMap(const Graph&) {
240        LEMON_ASSERT(false, "MatchingMap is not initialized");
241        return 0; // ignore warnings
242      }
243    };
244
245    /// \brief \ref named-templ-param "Named parameter" for setting
246    /// MatchingMap type
247    ///
248    /// \ref named-templ-param "Named parameter" for setting MatchingMap
249    /// type.
250    template <typename T>
251    struct SetMatchingMap
252      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
253      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
254    };
255
256    template <typename T>
257    struct SetElevatorTraits : public Traits {
258      typedef T Elevator;
259      static Elevator *createElevator(const Graph&, int) {
260        LEMON_ASSERT(false, "Elevator is not initialized");
261        return 0; // ignore warnings
262      }
263    };
264
265    /// \brief \ref named-templ-param "Named parameter" for setting
266    /// Elevator type
267    ///
268    /// \ref named-templ-param "Named parameter" for setting Elevator
269    /// type. If this named parameter is used, then an external
270    /// elevator object must be passed to the algorithm using the
271    /// \ref elevator(Elevator&) "elevator()" function before calling
272    /// \ref run() or \ref init().
273    /// \sa SetStandardElevator
274    template <typename T>
275    struct SetElevator
276      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
277      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
278    };
279
280    template <typename T>
281    struct SetStandardElevatorTraits : public Traits {
282      typedef T Elevator;
283      static Elevator *createElevator(const Graph& graph, int max_level) {
284        return new Elevator(graph, max_level);
285      }
286    };
287
288    /// \brief \ref named-templ-param "Named parameter" for setting
289    /// Elevator type with automatic allocation
290    ///
291    /// \ref named-templ-param "Named parameter" for setting Elevator
292    /// type with automatic allocation.
293    /// The Elevator should have standard constructor interface to be
294    /// able to automatically created by the algorithm (i.e. the
295    /// graph and the maximum level should be passed to it).
296    /// However an external elevator object could also be passed to the
297    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
298    /// before calling \ref run() or \ref init().
299    /// \sa SetElevator
300    template <typename T>
301    struct SetStandardElevator
302      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
303      typedef MaxFractionalMatching<Graph,
304                                    SetStandardElevatorTraits<T> > Create;
305    };
306
307    /// @}
308
309  protected:
310
311    MaxFractionalMatching() {}
312
313  public:
314
315    /// \brief Constructor
316    ///
317    /// Constructor.
318    ///
319    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
320      : _graph(graph), _allow_loops(allow_loops),
321        _local_matching(false), _matching(0),
322        _local_level(false), _level(0),  _indeg(0)
323    {}
324
325    ~MaxFractionalMatching() {
326      destroyStructures();
327    }
328
329    /// \brief Sets the matching map.
330    ///
331    /// Sets the matching map.
332    /// If you don't use this function before calling \ref run() or
333    /// \ref init(), an instance will be allocated automatically.
334    /// The destructor deallocates this automatically allocated map,
335    /// of course.
336    /// \return <tt>(*this)</tt>
337    MaxFractionalMatching& matchingMap(MatchingMap& map) {
338      if (_local_matching) {
339        delete _matching;
340        _local_matching = false;
341      }
342      _matching = &map;
343      return *this;
344    }
345
346    /// \brief Sets the elevator used by algorithm.
347    ///
348    /// Sets the elevator used by algorithm.
349    /// If you don't use this function before calling \ref run() or
350    /// \ref init(), an instance will be allocated automatically.
351    /// The destructor deallocates this automatically allocated elevator,
352    /// of course.
353    /// \return <tt>(*this)</tt>
354    MaxFractionalMatching& elevator(Elevator& elevator) {
355      if (_local_level) {
356        delete _level;
357        _local_level = false;
358      }
359      _level = &elevator;
360      return *this;
361    }
362
363    /// \brief Returns a const reference to the elevator.
364    ///
365    /// Returns a const reference to the elevator.
366    ///
367    /// \pre Either \ref run() or \ref init() must be called before
368    /// using this function.
369    const Elevator& elevator() const {
370      return *_level;
371    }
372
373    /// \name Execution control
374    /// The simplest way to execute the algorithm is to use one of the
375    /// member functions called \c run(). \n
376    /// If you need more control on the execution, first
377    /// you must call \ref init() and then one variant of the start()
378    /// member.
379
380    /// @{
381
382    /// \brief Initializes the internal data structures.
383    ///
384    /// Initializes the internal data structures and sets the initial
385    /// matching.
386    void init() {
387      createStructures();
388
389      _level->initStart();
390      for (NodeIt n(_graph); n != INVALID; ++n) {
391        _indeg->set(n, 0);
392        _matching->set(n, INVALID);
393        _level->initAddItem(n);
394      }
395      _level->initFinish();
396
397      _empty_level = _node_num;
398      for (NodeIt n(_graph); n != INVALID; ++n) {
399        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
400          if (_graph.target(a) == n && !_allow_loops) continue;
401          _matching->set(n, a);
402          Node v = _graph.target((*_matching)[n]);
403          _indeg->set(v, (*_indeg)[v] + 1);
404          break;
405        }
406      }
407
408      for (NodeIt n(_graph); n != INVALID; ++n) {
409        if ((*_indeg)[n] == 0) {
410          _level->activate(n);
411        }
412      }
413    }
414
415    /// \brief Starts the algorithm and computes a fractional matching
416    ///
417    /// The algorithm computes a maximum fractional matching.
418    ///
419    /// \param postprocess The algorithm computes first a matching
420    /// which is a union of a matching with one value edges, cycles
421    /// with half value edges and even length paths with half value
422    /// edges. If the parameter is true, then after the push-relabel
423    /// algorithm it postprocesses the matching to contain only
424    /// matching edges and half value odd cycles.
425    void start(bool postprocess = true) {
426      Node n;
427      while ((n = _level->highestActive()) != INVALID) {
428        int level = _level->highestActiveLevel();
429        int new_level = _level->maxLevel();
430        for (InArcIt a(_graph, n); a != INVALID; ++a) {
431          Node u = _graph.source(a);
432          if (n == u && !_allow_loops) continue;
433          Node v = _graph.target((*_matching)[u]);
434          if ((*_level)[v] < level) {
435            _indeg->set(v, (*_indeg)[v] - 1);
436            if ((*_indeg)[v] == 0) {
437              _level->activate(v);
438            }
439            _matching->set(u, a);
440            _indeg->set(n, (*_indeg)[n] + 1);
441            _level->deactivate(n);
442            goto no_more_push;
443          } else if (new_level > (*_level)[v]) {
444            new_level = (*_level)[v];
445          }
446        }
447
448        if (new_level + 1 < _level->maxLevel()) {
449          _level->liftHighestActive(new_level + 1);
450        } else {
451          _level->liftHighestActiveToTop();
452        }
453        if (_level->emptyLevel(level)) {
454          _level->liftToTop(level);
455        }
456      no_more_push:
457        ;
458      }
459      for (NodeIt n(_graph); n != INVALID; ++n) {
460        if ((*_matching)[n] == INVALID) continue;
461        Node u = _graph.target((*_matching)[n]);
462        if ((*_indeg)[u] > 1) {
463          _indeg->set(u, (*_indeg)[u] - 1);
464          _matching->set(n, INVALID);
465        }
466      }
467      if (postprocess) {
468        postprocessing();
469      }
470    }
471
472    /// \brief Starts the algorithm and computes a perfect fractional
473    /// matching
474    ///
475    /// The algorithm computes a perfect fractional matching. If it
476    /// does not exists, then the algorithm returns false and the
477    /// matching is undefined and the barrier.
478    ///
479    /// \param postprocess The algorithm computes first a matching
480    /// which is a union of a matching with one value edges, cycles
481    /// with half value edges and even length paths with half value
482    /// edges. If the parameter is true, then after the push-relabel
483    /// algorithm it postprocesses the matching to contain only
484    /// matching edges and half value odd cycles.
485    bool startPerfect(bool postprocess = true) {
486      Node n;
487      while ((n = _level->highestActive()) != INVALID) {
488        int level = _level->highestActiveLevel();
489        int new_level = _level->maxLevel();
490        for (InArcIt a(_graph, n); a != INVALID; ++a) {
491          Node u = _graph.source(a);
492          if (n == u && !_allow_loops) continue;
493          Node v = _graph.target((*_matching)[u]);
494          if ((*_level)[v] < level) {
495            _indeg->set(v, (*_indeg)[v] - 1);
496            if ((*_indeg)[v] == 0) {
497              _level->activate(v);
498            }
499            _matching->set(u, a);
500            _indeg->set(n, (*_indeg)[n] + 1);
501            _level->deactivate(n);
502            goto no_more_push;
503          } else if (new_level > (*_level)[v]) {
504            new_level = (*_level)[v];
505          }
506        }
507
508        if (new_level + 1 < _level->maxLevel()) {
509          _level->liftHighestActive(new_level + 1);
510        } else {
511          _level->liftHighestActiveToTop();
512          _empty_level = _level->maxLevel() - 1;
513          return false;
514        }
515        if (_level->emptyLevel(level)) {
516          _level->liftToTop(level);
517          _empty_level = level;
518          return false;
519        }
520      no_more_push:
521        ;
522      }
523      if (postprocess) {
524        postprocessing();
525      }
526      return true;
527    }
528
529    /// \brief Runs the algorithm
530    ///
531    /// Just a shortcut for the next code:
532    ///\code
533    /// init();
534    /// start();
535    ///\endcode
536    void run(bool postprocess = true) {
537      init();
538      start(postprocess);
539    }
540
541    /// \brief Runs the algorithm to find a perfect fractional matching
542    ///
543    /// Just a shortcut for the next code:
544    ///\code
545    /// init();
546    /// startPerfect();
547    ///\endcode
548    bool runPerfect(bool postprocess = true) {
549      init();
550      return startPerfect(postprocess);
551    }
552
553    ///@}
554
555    /// \name Query Functions
556    /// The result of the %Matching algorithm can be obtained using these
557    /// functions.\n
558    /// Before the use of these functions,
559    /// either run() or start() must be called.
560    ///@{
561
562
563    /// \brief Return the number of covered nodes in the matching.
564    ///
565    /// This function returns the number of covered nodes in the matching.
566    ///
567    /// \pre Either run() or start() must be called before using this function.
568    int matchingSize() const {
569      int num = 0;
570      for (NodeIt n(_graph); n != INVALID; ++n) {
571        if ((*_matching)[n] != INVALID) {
572          ++num;
573        }
574      }
575      return num;
576    }
577
578    /// \brief Returns a const reference to the matching map.
579    ///
580    /// Returns a const reference to the node map storing the found
581    /// fractional matching. This method can be called after
582    /// running the algorithm.
583    ///
584    /// \pre Either \ref run() or \ref init() must be called before
585    /// using this function.
586    const MatchingMap& matchingMap() const {
587      return *_matching;
588    }
589
590    /// \brief Return \c true if the given edge is in the matching.
591    ///
592    /// This function returns \c true if the given edge is in the
593    /// found matching. The result is scaled by \ref primalScale
594    /// "primal scale".
595    ///
596    /// \pre Either run() or start() must be called before using this function.
597    int matching(const Edge& edge) const {
598      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
599        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
600    }
601
602    /// \brief Return the fractional matching arc (or edge) incident
603    /// to the given node.
604    ///
605    /// This function returns one of the fractional matching arc (or
606    /// edge) incident to the given node in the found matching or \c
607    /// INVALID if the node is not covered by the matching or if the
608    /// node is on an odd length cycle then it is the successor edge
609    /// on the cycle.
610    ///
611    /// \pre Either run() or start() must be called before using this function.
612    Arc matching(const Node& node) const {
613      return (*_matching)[node];
614    }
615
616    /// \brief Returns true if the node is in the barrier
617    ///
618    /// The barrier is a subset of the nodes. If the nodes in the
619    /// barrier have less adjacent nodes than the size of the barrier,
620    /// then at least as much nodes cannot be covered as the
621    /// difference of the two subsets.
622    bool barrier(const Node& node) const {
623      return (*_level)[node] >= _empty_level;
624    }
625
626    /// @}
627
628  };
629
630  /// \ingroup matching
631  ///
632  /// \brief Weighted fractional matching in general graphs
633  ///
634  /// This class provides an efficient implementation of fractional
635  /// matching algorithm. The implementation is based on extensive use
636  /// of priority queues and provides \f$O(nm\log n)\f$ time
637  /// complexity.
638  ///
639  /// The maximum weighted fractional matching is a relaxation of the
640  /// maximum weighted matching problem where the odd set constraints
641  /// are omitted.
642  /// It can be formulated with the following linear program.
643  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
644  /// \f[x_e \ge 0\quad \forall e\in E\f]
645  /// \f[\max \sum_{e\in E}x_ew_e\f]
646  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
647  /// \f$X\f$. The result must be the union of a matching with one
648  /// value edges and a set of odd length cycles with half value edges.
649  ///
650  /// The algorithm calculates an optimal fractional matching and a
651  /// proof of the optimality. The solution of the dual problem can be
652  /// used to check the result of the algorithm. The dual linear
653  /// problem is the following.
654  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
655  /// \f[y_u \ge 0 \quad \forall u \in V\f]
656  /// \f[\min \sum_{u \in V}y_u \f] ///
657  ///
658  /// The algorithm can be executed with the run() function.
659  /// After it the matching (the primal solution) and the dual solution
660  /// can be obtained using the query functions.
661  ///
662  /// If the value type is integer, then the primal and the dual
663  /// solutions are multiplied by
664  /// \ref MaxWeightedMatching::primalScale "2" and
665  /// \ref MaxWeightedMatching::dualScale "4" respectively.
666  ///
667  /// \tparam GR The undirected graph type the algorithm runs on.
668  /// \tparam WM The type edge weight map. The default type is
669  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
670#ifdef DOXYGEN
671  template <typename GR, typename WM>
672#else
673  template <typename GR,
674            typename WM = typename GR::template EdgeMap<int> >
675#endif
676  class MaxWeightedFractionalMatching {
677  public:
678
679    /// The graph type of the algorithm
680    typedef GR Graph;
681    /// The type of the edge weight map
682    typedef WM WeightMap;
683    /// The value type of the edge weights
684    typedef typename WeightMap::Value Value;
685
686    /// The type of the matching map
687    typedef typename Graph::template NodeMap<typename Graph::Arc>
688    MatchingMap;
689
690    /// \brief Scaling factor for primal solution
691    ///
692    /// Scaling factor for primal solution. It is equal to 2 or 1
693    /// according to the value type.
694    static const int primalScale =
695      std::numeric_limits<Value>::is_integer ? 2 : 1;
696
697    /// \brief Scaling factor for dual solution
698    ///
699    /// Scaling factor for dual solution. It is equal to 4 or 1
700    /// according to the value type.
701    static const int dualScale =
702      std::numeric_limits<Value>::is_integer ? 4 : 1;
703
704  private:
705
706    TEMPLATE_GRAPH_TYPEDEFS(Graph);
707
708    typedef typename Graph::template NodeMap<Value> NodePotential;
709
710    const Graph& _graph;
711    const WeightMap& _weight;
712
713    MatchingMap* _matching;
714    NodePotential* _node_potential;
715
716    int _node_num;
717    bool _allow_loops;
718
719    enum Status {
720      EVEN = -1, MATCHED = 0, ODD = 1
721    };
722
723    typedef typename Graph::template NodeMap<Status> StatusMap;
724    StatusMap* _status;
725
726    typedef typename Graph::template NodeMap<Arc> PredMap;
727    PredMap* _pred;
728
729    typedef ExtendFindEnum<IntNodeMap> TreeSet;
730
731    IntNodeMap *_tree_set_index;
732    TreeSet *_tree_set;
733
734    IntNodeMap *_delta1_index;
735    BinHeap<Value, IntNodeMap> *_delta1;
736
737    IntNodeMap *_delta2_index;
738    BinHeap<Value, IntNodeMap> *_delta2;
739
740    IntEdgeMap *_delta3_index;
741    BinHeap<Value, IntEdgeMap> *_delta3;
742
743    Value _delta_sum;
744
745    void createStructures() {
746      _node_num = countNodes(_graph);
747
748      if (!_matching) {
749        _matching = new MatchingMap(_graph);
750      }
751      if (!_node_potential) {
752        _node_potential = new NodePotential(_graph);
753      }
754      if (!_status) {
755        _status = new StatusMap(_graph);
756      }
757      if (!_pred) {
758        _pred = new PredMap(_graph);
759      }
760      if (!_tree_set) {
761        _tree_set_index = new IntNodeMap(_graph);
762        _tree_set = new TreeSet(*_tree_set_index);
763      }
764      if (!_delta1) {
765        _delta1_index = new IntNodeMap(_graph);
766        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
767      }
768      if (!_delta2) {
769        _delta2_index = new IntNodeMap(_graph);
770        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
771      }
772      if (!_delta3) {
773        _delta3_index = new IntEdgeMap(_graph);
774        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
775      }
776    }
777
778    void destroyStructures() {
779      if (_matching) {
780        delete _matching;
781      }
782      if (_node_potential) {
783        delete _node_potential;
784      }
785      if (_status) {
786        delete _status;
787      }
788      if (_pred) {
789        delete _pred;
790      }
791      if (_tree_set) {
792        delete _tree_set_index;
793        delete _tree_set;
794      }
795      if (_delta1) {
796        delete _delta1_index;
797        delete _delta1;
798      }
799      if (_delta2) {
800        delete _delta2_index;
801        delete _delta2;
802      }
803      if (_delta3) {
804        delete _delta3_index;
805        delete _delta3;
806      }
807    }
808
809    void matchedToEven(Node node, int tree) {
810      _tree_set->insert(node, tree);
811      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
812      _delta1->push(node, (*_node_potential)[node]);
813
814      if (_delta2->state(node) == _delta2->IN_HEAP) {
815        _delta2->erase(node);
816      }
817
818      for (InArcIt a(_graph, node); a != INVALID; ++a) {
819        Node v = _graph.source(a);
820        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
821          dualScale * _weight[a];
822        if (node == v) {
823          if (_allow_loops && _graph.direction(a)) {
824            _delta3->push(a, rw / 2);
825          }
826        } else if ((*_status)[v] == EVEN) {
827          _delta3->push(a, rw / 2);
828        } else if ((*_status)[v] == MATCHED) {
829          if (_delta2->state(v) != _delta2->IN_HEAP) {
830            _pred->set(v, a);
831            _delta2->push(v, rw);
832          } else if ((*_delta2)[v] > rw) {
833            _pred->set(v, a);
834            _delta2->decrease(v, rw);
835          }
836        }
837      }
838    }
839
840    void matchedToOdd(Node node, int tree) {
841      _tree_set->insert(node, tree);
842      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
843
844      if (_delta2->state(node) == _delta2->IN_HEAP) {
845        _delta2->erase(node);
846      }
847    }
848
849    void evenToMatched(Node node, int tree) {
850      _delta1->erase(node);
851      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
852      Arc min = INVALID;
853      Value minrw = std::numeric_limits<Value>::max();
854      for (InArcIt a(_graph, node); a != INVALID; ++a) {
855        Node v = _graph.source(a);
856        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
857          dualScale * _weight[a];
858
859        if (node == v) {
860          if (_allow_loops && _graph.direction(a)) {
861            _delta3->erase(a);
862          }
863        } else if ((*_status)[v] == EVEN) {
864          _delta3->erase(a);
865          if (minrw > rw) {
866            min = _graph.oppositeArc(a);
867            minrw = rw;
868          }
869        } else if ((*_status)[v]  == MATCHED) {
870          if ((*_pred)[v] == a) {
871            Arc mina = INVALID;
872            Value minrwa = std::numeric_limits<Value>::max();
873            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
874              Node va = _graph.target(aa);
875              if ((*_status)[va] != EVEN ||
876                  _tree_set->find(va) == tree) continue;
877              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
878                dualScale * _weight[aa];
879              if (minrwa > rwa) {
880                minrwa = rwa;
881                mina = aa;
882              }
883            }
884            if (mina != INVALID) {
885              _pred->set(v, mina);
886              _delta2->increase(v, minrwa);
887            } else {
888              _pred->set(v, INVALID);
889              _delta2->erase(v);
890            }
891          }
892        }
893      }
894      if (min != INVALID) {
895        _pred->set(node, min);
896        _delta2->push(node, minrw);
897      } else {
898        _pred->set(node, INVALID);
899      }
900    }
901
902    void oddToMatched(Node node) {
903      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
904      Arc min = INVALID;
905      Value minrw = std::numeric_limits<Value>::max();
906      for (InArcIt a(_graph, node); a != INVALID; ++a) {
907        Node v = _graph.source(a);
908        if ((*_status)[v] != EVEN) continue;
909        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
910          dualScale * _weight[a];
911
912        if (minrw > rw) {
913          min = _graph.oppositeArc(a);
914          minrw = rw;
915        }
916      }
917      if (min != INVALID) {
918        _pred->set(node, min);
919        _delta2->push(node, minrw);
920      } else {
921        _pred->set(node, INVALID);
922      }
923    }
924
925    void alternatePath(Node even, int tree) {
926      Node odd;
927
928      _status->set(even, MATCHED);
929      evenToMatched(even, tree);
930
931      Arc prev = (*_matching)[even];
932      while (prev != INVALID) {
933        odd = _graph.target(prev);
934        even = _graph.target((*_pred)[odd]);
935        _matching->set(odd, (*_pred)[odd]);
936        _status->set(odd, MATCHED);
937        oddToMatched(odd);
938
939        prev = (*_matching)[even];
940        _status->set(even, MATCHED);
941        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
942        evenToMatched(even, tree);
943      }
944    }
945
946    void destroyTree(int tree) {
947      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
948        if ((*_status)[n] == EVEN) {
949          _status->set(n, MATCHED);
950          evenToMatched(n, tree);
951        } else if ((*_status)[n] == ODD) {
952          _status->set(n, MATCHED);
953          oddToMatched(n);
954        }
955      }
956      _tree_set->eraseClass(tree);
957    }
958
959
960    void unmatchNode(const Node& node) {
961      int tree = _tree_set->find(node);
962
963      alternatePath(node, tree);
964      destroyTree(tree);
965
966      _matching->set(node, INVALID);
967    }
968
969
970    void augmentOnEdge(const Edge& edge) {
971      Node left = _graph.u(edge);
972      int left_tree = _tree_set->find(left);
973
974      alternatePath(left, left_tree);
975      destroyTree(left_tree);
976      _matching->set(left, _graph.direct(edge, true));
977
978      Node right = _graph.v(edge);
979      int right_tree = _tree_set->find(right);
980
981      alternatePath(right, right_tree);
982      destroyTree(right_tree);
983      _matching->set(right, _graph.direct(edge, false));
984    }
985
986    void augmentOnArc(const Arc& arc) {
987      Node left = _graph.source(arc);
988      _status->set(left, MATCHED);
989      _matching->set(left, arc);
990      _pred->set(left, arc);
991
992      Node right = _graph.target(arc);
993      int right_tree = _tree_set->find(right);
994
995      alternatePath(right, right_tree);
996      destroyTree(right_tree);
997      _matching->set(right, _graph.oppositeArc(arc));
998    }
999
1000    void extendOnArc(const Arc& arc) {
1001      Node base = _graph.target(arc);
1002      int tree = _tree_set->find(base);
1003
1004      Node odd = _graph.source(arc);
1005      _tree_set->insert(odd, tree);
1006      _status->set(odd, ODD);
1007      matchedToOdd(odd, tree);
1008      _pred->set(odd, arc);
1009
1010      Node even = _graph.target((*_matching)[odd]);
1011      _tree_set->insert(even, tree);
1012      _status->set(even, EVEN);
1013      matchedToEven(even, tree);
1014    }
1015
1016    void cycleOnEdge(const Edge& edge, int tree) {
1017      Node nca = INVALID;
1018      std::vector<Node> left_path, right_path;
1019
1020      {
1021        std::set<Node> left_set, right_set;
1022        Node left = _graph.u(edge);
1023        left_path.push_back(left);
1024        left_set.insert(left);
1025
1026        Node right = _graph.v(edge);
1027        right_path.push_back(right);
1028        right_set.insert(right);
1029
1030        while (true) {
1031
1032          if (left_set.find(right) != left_set.end()) {
1033            nca = right;
1034            break;
1035          }
1036
1037          if ((*_matching)[left] == INVALID) break;
1038
1039          left = _graph.target((*_matching)[left]);
1040          left_path.push_back(left);
1041          left = _graph.target((*_pred)[left]);
1042          left_path.push_back(left);
1043
1044          left_set.insert(left);
1045
1046          if (right_set.find(left) != right_set.end()) {
1047            nca = left;
1048            break;
1049          }
1050
1051          if ((*_matching)[right] == INVALID) break;
1052
1053          right = _graph.target((*_matching)[right]);
1054          right_path.push_back(right);
1055          right = _graph.target((*_pred)[right]);
1056          right_path.push_back(right);
1057
1058          right_set.insert(right);
1059
1060        }
1061
1062        if (nca == INVALID) {
1063          if ((*_matching)[left] == INVALID) {
1064            nca = right;
1065            while (left_set.find(nca) == left_set.end()) {
1066              nca = _graph.target((*_matching)[nca]);
1067              right_path.push_back(nca);
1068              nca = _graph.target((*_pred)[nca]);
1069              right_path.push_back(nca);
1070            }
1071          } else {
1072            nca = left;
1073            while (right_set.find(nca) == right_set.end()) {
1074              nca = _graph.target((*_matching)[nca]);
1075              left_path.push_back(nca);
1076              nca = _graph.target((*_pred)[nca]);
1077              left_path.push_back(nca);
1078            }
1079          }
1080        }
1081      }
1082
1083      alternatePath(nca, tree);
1084      Arc prev;
1085
1086      prev = _graph.direct(edge, true);
1087      for (int i = 0; left_path[i] != nca; i += 2) {
1088        _matching->set(left_path[i], prev);
1089        _status->set(left_path[i], MATCHED);
1090        evenToMatched(left_path[i], tree);
1091
1092        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1093        _status->set(left_path[i + 1], MATCHED);
1094        oddToMatched(left_path[i + 1]);
1095      }
1096      _matching->set(nca, prev);
1097
1098      for (int i = 0; right_path[i] != nca; i += 2) {
1099        _status->set(right_path[i], MATCHED);
1100        evenToMatched(right_path[i], tree);
1101
1102        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1103        _status->set(right_path[i + 1], MATCHED);
1104        oddToMatched(right_path[i + 1]);
1105      }
1106
1107      destroyTree(tree);
1108    }
1109
1110    void extractCycle(const Arc &arc) {
1111      Node left = _graph.source(arc);
1112      Node odd = _graph.target((*_matching)[left]);
1113      Arc prev;
1114      while (odd != left) {
1115        Node even = _graph.target((*_matching)[odd]);
1116        prev = (*_matching)[odd];
1117        odd = _graph.target((*_matching)[even]);
1118        _matching->set(even, _graph.oppositeArc(prev));
1119      }
1120      _matching->set(left, arc);
1121
1122      Node right = _graph.target(arc);
1123      int right_tree = _tree_set->find(right);
1124      alternatePath(right, right_tree);
1125      destroyTree(right_tree);
1126      _matching->set(right, _graph.oppositeArc(arc));
1127    }
1128
1129  public:
1130
1131    /// \brief Constructor
1132    ///
1133    /// Constructor.
1134    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
1135                                  bool allow_loops = true)
1136      : _graph(graph), _weight(weight), _matching(0),
1137      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1138      _status(0),  _pred(0),
1139      _tree_set_index(0), _tree_set(0),
1140
1141      _delta1_index(0), _delta1(0),
1142      _delta2_index(0), _delta2(0),
1143      _delta3_index(0), _delta3(0),
1144
1145      _delta_sum() {}
1146
1147    ~MaxWeightedFractionalMatching() {
1148      destroyStructures();
1149    }
1150
1151    /// \name Execution Control
1152    /// The simplest way to execute the algorithm is to use the
1153    /// \ref run() member function.
1154
1155    ///@{
1156
1157    /// \brief Initialize the algorithm
1158    ///
1159    /// This function initializes the algorithm.
1160    void init() {
1161      createStructures();
1162
1163      for (NodeIt n(_graph); n != INVALID; ++n) {
1164        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1165        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1166      }
1167      for (EdgeIt e(_graph); e != INVALID; ++e) {
1168        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1169      }
1170
1171      for (NodeIt n(_graph); n != INVALID; ++n) {
1172        Value max = 0;
1173        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1174          if (_graph.target(e) == n && !_allow_loops) continue;
1175          if ((dualScale * _weight[e]) / 2 > max) {
1176            max = (dualScale * _weight[e]) / 2;
1177          }
1178        }
1179        _node_potential->set(n, max);
1180        _delta1->push(n, max);
1181
1182        _tree_set->insert(n);
1183
1184        _matching->set(n, INVALID);
1185        _status->set(n, EVEN);
1186      }
1187
1188      for (EdgeIt e(_graph); e != INVALID; ++e) {
1189        Node left = _graph.u(e);
1190        Node right = _graph.v(e);
1191        if (left == right && !_allow_loops) continue;
1192        _delta3->push(e, ((*_node_potential)[left] +
1193                          (*_node_potential)[right] -
1194                          dualScale * _weight[e]) / 2);
1195      }
1196    }
1197
1198    /// \brief Start the algorithm
1199    ///
1200    /// This function starts the algorithm.
1201    ///
1202    /// \pre \ref init() must be called before using this function.
1203    void start() {
1204      enum OpType {
1205        D1, D2, D3
1206      };
1207
1208      int unmatched = _node_num;
1209      while (unmatched > 0) {
1210        Value d1 = !_delta1->empty() ?
1211          _delta1->prio() : std::numeric_limits<Value>::max();
1212
1213        Value d2 = !_delta2->empty() ?
1214          _delta2->prio() : std::numeric_limits<Value>::max();
1215
1216        Value d3 = !_delta3->empty() ?
1217          _delta3->prio() : std::numeric_limits<Value>::max();
1218
1219        _delta_sum = d3; OpType ot = D3;
1220        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1221        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1222
1223        switch (ot) {
1224        case D1:
1225          {
1226            Node n = _delta1->top();
1227            unmatchNode(n);
1228            --unmatched;
1229          }
1230          break;
1231        case D2:
1232          {
1233            Node n = _delta2->top();
1234            Arc a = (*_pred)[n];
1235            if ((*_matching)[n] == INVALID) {
1236              augmentOnArc(a);
1237              --unmatched;
1238            } else {
1239              Node v = _graph.target((*_matching)[n]);
1240              if ((*_matching)[n] !=
1241                  _graph.oppositeArc((*_matching)[v])) {
1242                extractCycle(a);
1243                --unmatched;
1244              } else {
1245                extendOnArc(a);
1246              }
1247            }
1248          } break;
1249        case D3:
1250          {
1251            Edge e = _delta3->top();
1252
1253            Node left = _graph.u(e);
1254            Node right = _graph.v(e);
1255
1256            int left_tree = _tree_set->find(left);
1257            int right_tree = _tree_set->find(right);
1258
1259            if (left_tree == right_tree) {
1260              cycleOnEdge(e, left_tree);
1261              --unmatched;
1262            } else {
1263              augmentOnEdge(e);
1264              unmatched -= 2;
1265            }
1266          } break;
1267        }
1268      }
1269    }
1270
1271    /// \brief Run the algorithm.
1272    ///
1273    /// This method runs the \c %MaxWeightedMatching algorithm.
1274    ///
1275    /// \note mwfm.run() is just a shortcut of the following code.
1276    /// \code
1277    ///   mwfm.init();
1278    ///   mwfm.start();
1279    /// \endcode
1280    void run() {
1281      init();
1282      start();
1283    }
1284
1285    /// @}
1286
1287    /// \name Primal Solution
1288    /// Functions to get the primal solution, i.e. the maximum weighted
1289    /// matching.\n
1290    /// Either \ref run() or \ref start() function should be called before
1291    /// using them.
1292
1293    /// @{
1294
1295    /// \brief Return the weight of the matching.
1296    ///
1297    /// This function returns the weight of the found matching. This
1298    /// value is scaled by \ref primalScale "primal scale".
1299    ///
1300    /// \pre Either run() or start() must be called before using this function.
1301    Value matchingWeight() const {
1302      Value sum = 0;
1303      for (NodeIt n(_graph); n != INVALID; ++n) {
1304        if ((*_matching)[n] != INVALID) {
1305          sum += _weight[(*_matching)[n]];
1306        }
1307      }
1308      return sum * primalScale / 2;
1309    }
1310
1311    /// \brief Return the number of covered nodes in the matching.
1312    ///
1313    /// This function returns the number of covered nodes in the matching.
1314    ///
1315    /// \pre Either run() or start() must be called before using this function.
1316    int matchingSize() const {
1317      int num = 0;
1318      for (NodeIt n(_graph); n != INVALID; ++n) {
1319        if ((*_matching)[n] != INVALID) {
1320          ++num;
1321        }
1322      }
1323      return num;
1324    }
1325
1326    /// \brief Return \c true if the given edge is in the matching.
1327    ///
1328    /// This function returns \c true if the given edge is in the
1329    /// found matching. The result is scaled by \ref primalScale
1330    /// "primal scale".
1331    ///
1332    /// \pre Either run() or start() must be called before using this function.
1333    Value matching(const Edge& edge) const {
1334      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1335        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
1336        * primalScale / 2;
1337    }
1338
1339    /// \brief Return the fractional matching arc (or edge) incident
1340    /// to the given node.
1341    ///
1342    /// This function returns one of the fractional matching arc (or
1343    /// edge) incident to the given node in the found matching or \c
1344    /// INVALID if the node is not covered by the matching or if the
1345    /// node is on an odd length cycle then it is the successor edge
1346    /// on the cycle.
1347    ///
1348    /// \pre Either run() or start() must be called before using this function.
1349    Arc matching(const Node& node) const {
1350      return (*_matching)[node];
1351    }
1352
1353    /// \brief Return a const reference to the matching map.
1354    ///
1355    /// This function returns a const reference to a node map that stores
1356    /// the matching arc (or edge) incident to each node.
1357    const MatchingMap& matchingMap() const {
1358      return *_matching;
1359    }
1360
1361    /// @}
1362
1363    /// \name Dual Solution
1364    /// Functions to get the dual solution.\n
1365    /// Either \ref run() or \ref start() function should be called before
1366    /// using them.
1367
1368    /// @{
1369
1370    /// \brief Return the value of the dual solution.
1371    ///
1372    /// This function returns the value of the dual solution.
1373    /// It should be equal to the primal value scaled by \ref dualScale
1374    /// "dual scale".
1375    ///
1376    /// \pre Either run() or start() must be called before using this function.
1377    Value dualValue() const {
1378      Value sum = 0;
1379      for (NodeIt n(_graph); n != INVALID; ++n) {
1380        sum += nodeValue(n);
1381      }
1382      return sum;
1383    }
1384
1385    /// \brief Return the dual value (potential) of the given node.
1386    ///
1387    /// This function returns the dual value (potential) of the given node.
1388    ///
1389    /// \pre Either run() or start() must be called before using this function.
1390    Value nodeValue(const Node& n) const {
1391      return (*_node_potential)[n];
1392    }
1393
1394    /// @}
1395
1396  };
1397
1398  /// \ingroup matching
1399  ///
1400  /// \brief Weighted fractional perfect matching in general graphs
1401  ///
1402  /// This class provides an efficient implementation of fractional
1403  /// matching algorithm. The implementation is based on extensive use
1404  /// of priority queues and provides \f$O(nm\log n)\f$ time
1405  /// complexity.
1406  ///
1407  /// The maximum weighted fractional perfect matching is a relaxation
1408  /// of the maximum weighted perfect matching problem where the odd
1409  /// set constraints are omitted.
1410  /// It can be formulated with the following linear program.
1411  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1412  /// \f[x_e \ge 0\quad \forall e\in E\f]
1413  /// \f[\max \sum_{e\in E}x_ew_e\f]
1414  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1415  /// \f$X\f$. The result must be the union of a matching with one
1416  /// value edges and a set of odd length cycles with half value edges.
1417  ///
1418  /// The algorithm calculates an optimal fractional matching and a
1419  /// proof of the optimality. The solution of the dual problem can be
1420  /// used to check the result of the algorithm. The dual linear
1421  /// problem is the following.
1422  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1423  /// \f[\min \sum_{u \in V}y_u \f] ///
1424  ///
1425  /// The algorithm can be executed with the run() function.
1426  /// After it the matching (the primal solution) and the dual solution
1427  /// can be obtained using the query functions.
1428
1429  /// If the value type is integer, then the primal and the dual
1430  /// solutions are multiplied by
1431  /// \ref MaxWeightedMatching::primalScale "2" and
1432  /// \ref MaxWeightedMatching::dualScale "4" respectively.
1433  ///
1434  /// \tparam GR The undirected graph type the algorithm runs on.
1435  /// \tparam WM The type edge weight map. The default type is
1436  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1437#ifdef DOXYGEN
1438  template <typename GR, typename WM>
1439#else
1440  template <typename GR,
1441            typename WM = typename GR::template EdgeMap<int> >
1442#endif
1443  class MaxWeightedPerfectFractionalMatching {
1444  public:
1445
1446    /// The graph type of the algorithm
1447    typedef GR Graph;
1448    /// The type of the edge weight map
1449    typedef WM WeightMap;
1450    /// The value type of the edge weights
1451    typedef typename WeightMap::Value Value;
1452
1453    /// The type of the matching map
1454    typedef typename Graph::template NodeMap<typename Graph::Arc>
1455    MatchingMap;
1456
1457    /// \brief Scaling factor for primal solution
1458    ///
1459    /// Scaling factor for primal solution. It is equal to 2 or 1
1460    /// according to the value type.
1461    static const int primalScale =
1462      std::numeric_limits<Value>::is_integer ? 2 : 1;
1463
1464    /// \brief Scaling factor for dual solution
1465    ///
1466    /// Scaling factor for dual solution. It is equal to 4 or 1
1467    /// according to the value type.
1468    static const int dualScale =
1469      std::numeric_limits<Value>::is_integer ? 4 : 1;
1470
1471  private:
1472
1473    TEMPLATE_GRAPH_TYPEDEFS(Graph);
1474
1475    typedef typename Graph::template NodeMap<Value> NodePotential;
1476
1477    const Graph& _graph;
1478    const WeightMap& _weight;
1479
1480    MatchingMap* _matching;
1481    NodePotential* _node_potential;
1482
1483    int _node_num;
1484    bool _allow_loops;
1485
1486    enum Status {
1487      EVEN = -1, MATCHED = 0, ODD = 1
1488    };
1489
1490    typedef typename Graph::template NodeMap<Status> StatusMap;
1491    StatusMap* _status;
1492
1493    typedef typename Graph::template NodeMap<Arc> PredMap;
1494    PredMap* _pred;
1495
1496    typedef ExtendFindEnum<IntNodeMap> TreeSet;
1497
1498    IntNodeMap *_tree_set_index;
1499    TreeSet *_tree_set;
1500
1501    IntNodeMap *_delta2_index;
1502    BinHeap<Value, IntNodeMap> *_delta2;
1503
1504    IntEdgeMap *_delta3_index;
1505    BinHeap<Value, IntEdgeMap> *_delta3;
1506
1507    Value _delta_sum;
1508
1509    void createStructures() {
1510      _node_num = countNodes(_graph);
1511
1512      if (!_matching) {
1513        _matching = new MatchingMap(_graph);
1514      }
1515      if (!_node_potential) {
1516        _node_potential = new NodePotential(_graph);
1517      }
1518      if (!_status) {
1519        _status = new StatusMap(_graph);
1520      }
1521      if (!_pred) {
1522        _pred = new PredMap(_graph);
1523      }
1524      if (!_tree_set) {
1525        _tree_set_index = new IntNodeMap(_graph);
1526        _tree_set = new TreeSet(*_tree_set_index);
1527      }
1528      if (!_delta2) {
1529        _delta2_index = new IntNodeMap(_graph);
1530        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1531      }
1532      if (!_delta3) {
1533        _delta3_index = new IntEdgeMap(_graph);
1534        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1535      }
1536    }
1537
1538    void destroyStructures() {
1539      if (_matching) {
1540        delete _matching;
1541      }
1542      if (_node_potential) {
1543        delete _node_potential;
1544      }
1545      if (_status) {
1546        delete _status;
1547      }
1548      if (_pred) {
1549        delete _pred;
1550      }
1551      if (_tree_set) {
1552        delete _tree_set_index;
1553        delete _tree_set;
1554      }
1555      if (_delta2) {
1556        delete _delta2_index;
1557        delete _delta2;
1558      }
1559      if (_delta3) {
1560        delete _delta3_index;
1561        delete _delta3;
1562      }
1563    }
1564
1565    void matchedToEven(Node node, int tree) {
1566      _tree_set->insert(node, tree);
1567      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1568
1569      if (_delta2->state(node) == _delta2->IN_HEAP) {
1570        _delta2->erase(node);
1571      }
1572
1573      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1574        Node v = _graph.source(a);
1575        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1576          dualScale * _weight[a];
1577        if (node == v) {
1578          if (_allow_loops && _graph.direction(a)) {
1579            _delta3->push(a, rw / 2);
1580          }
1581        } else if ((*_status)[v] == EVEN) {
1582          _delta3->push(a, rw / 2);
1583        } else if ((*_status)[v] == MATCHED) {
1584          if (_delta2->state(v) != _delta2->IN_HEAP) {
1585            _pred->set(v, a);
1586            _delta2->push(v, rw);
1587          } else if ((*_delta2)[v] > rw) {
1588            _pred->set(v, a);
1589            _delta2->decrease(v, rw);
1590          }
1591        }
1592      }
1593    }
1594
1595    void matchedToOdd(Node node, int tree) {
1596      _tree_set->insert(node, tree);
1597      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1598
1599      if (_delta2->state(node) == _delta2->IN_HEAP) {
1600        _delta2->erase(node);
1601      }
1602    }
1603
1604    void evenToMatched(Node node, int tree) {
1605      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1606      Arc min = INVALID;
1607      Value minrw = std::numeric_limits<Value>::max();
1608      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1609        Node v = _graph.source(a);
1610        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1611          dualScale * _weight[a];
1612
1613        if (node == v) {
1614          if (_allow_loops && _graph.direction(a)) {
1615            _delta3->erase(a);
1616          }
1617        } else if ((*_status)[v] == EVEN) {
1618          _delta3->erase(a);
1619          if (minrw > rw) {
1620            min = _graph.oppositeArc(a);
1621            minrw = rw;
1622          }
1623        } else if ((*_status)[v]  == MATCHED) {
1624          if ((*_pred)[v] == a) {
1625            Arc mina = INVALID;
1626            Value minrwa = std::numeric_limits<Value>::max();
1627            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1628              Node va = _graph.target(aa);
1629              if ((*_status)[va] != EVEN ||
1630                  _tree_set->find(va) == tree) continue;
1631              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1632                dualScale * _weight[aa];
1633              if (minrwa > rwa) {
1634                minrwa = rwa;
1635                mina = aa;
1636              }
1637            }
1638            if (mina != INVALID) {
1639              _pred->set(v, mina);
1640              _delta2->increase(v, minrwa);
1641            } else {
1642              _pred->set(v, INVALID);
1643              _delta2->erase(v);
1644            }
1645          }
1646        }
1647      }
1648      if (min != INVALID) {
1649        _pred->set(node, min);
1650        _delta2->push(node, minrw);
1651      } else {
1652        _pred->set(node, INVALID);
1653      }
1654    }
1655
1656    void oddToMatched(Node node) {
1657      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1658      Arc min = INVALID;
1659      Value minrw = std::numeric_limits<Value>::max();
1660      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1661        Node v = _graph.source(a);
1662        if ((*_status)[v] != EVEN) continue;
1663        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1664          dualScale * _weight[a];
1665
1666        if (minrw > rw) {
1667          min = _graph.oppositeArc(a);
1668          minrw = rw;
1669        }
1670      }
1671      if (min != INVALID) {
1672        _pred->set(node, min);
1673        _delta2->push(node, minrw);
1674      } else {
1675        _pred->set(node, INVALID);
1676      }
1677    }
1678
1679    void alternatePath(Node even, int tree) {
1680      Node odd;
1681
1682      _status->set(even, MATCHED);
1683      evenToMatched(even, tree);
1684
1685      Arc prev = (*_matching)[even];
1686      while (prev != INVALID) {
1687        odd = _graph.target(prev);
1688        even = _graph.target((*_pred)[odd]);
1689        _matching->set(odd, (*_pred)[odd]);
1690        _status->set(odd, MATCHED);
1691        oddToMatched(odd);
1692
1693        prev = (*_matching)[even];
1694        _status->set(even, MATCHED);
1695        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1696        evenToMatched(even, tree);
1697      }
1698    }
1699
1700    void destroyTree(int tree) {
1701      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1702        if ((*_status)[n] == EVEN) {
1703          _status->set(n, MATCHED);
1704          evenToMatched(n, tree);
1705        } else if ((*_status)[n] == ODD) {
1706          _status->set(n, MATCHED);
1707          oddToMatched(n);
1708        }
1709      }
1710      _tree_set->eraseClass(tree);
1711    }
1712
1713    void augmentOnEdge(const Edge& edge) {
1714      Node left = _graph.u(edge);
1715      int left_tree = _tree_set->find(left);
1716
1717      alternatePath(left, left_tree);
1718      destroyTree(left_tree);
1719      _matching->set(left, _graph.direct(edge, true));
1720
1721      Node right = _graph.v(edge);
1722      int right_tree = _tree_set->find(right);
1723
1724      alternatePath(right, right_tree);
1725      destroyTree(right_tree);
1726      _matching->set(right, _graph.direct(edge, false));
1727    }
1728
1729    void augmentOnArc(const Arc& arc) {
1730      Node left = _graph.source(arc);
1731      _status->set(left, MATCHED);
1732      _matching->set(left, arc);
1733      _pred->set(left, arc);
1734
1735      Node right = _graph.target(arc);
1736      int right_tree = _tree_set->find(right);
1737
1738      alternatePath(right, right_tree);
1739      destroyTree(right_tree);
1740      _matching->set(right, _graph.oppositeArc(arc));
1741    }
1742
1743    void extendOnArc(const Arc& arc) {
1744      Node base = _graph.target(arc);
1745      int tree = _tree_set->find(base);
1746
1747      Node odd = _graph.source(arc);
1748      _tree_set->insert(odd, tree);
1749      _status->set(odd, ODD);
1750      matchedToOdd(odd, tree);
1751      _pred->set(odd, arc);
1752
1753      Node even = _graph.target((*_matching)[odd]);
1754      _tree_set->insert(even, tree);
1755      _status->set(even, EVEN);
1756      matchedToEven(even, tree);
1757    }
1758
1759    void cycleOnEdge(const Edge& edge, int tree) {
1760      Node nca = INVALID;
1761      std::vector<Node> left_path, right_path;
1762
1763      {
1764        std::set<Node> left_set, right_set;
1765        Node left = _graph.u(edge);
1766        left_path.push_back(left);
1767        left_set.insert(left);
1768
1769        Node right = _graph.v(edge);
1770        right_path.push_back(right);
1771        right_set.insert(right);
1772
1773        while (true) {
1774
1775          if (left_set.find(right) != left_set.end()) {
1776            nca = right;
1777            break;
1778          }
1779
1780          if ((*_matching)[left] == INVALID) break;
1781
1782          left = _graph.target((*_matching)[left]);
1783          left_path.push_back(left);
1784          left = _graph.target((*_pred)[left]);
1785          left_path.push_back(left);
1786
1787          left_set.insert(left);
1788
1789          if (right_set.find(left) != right_set.end()) {
1790            nca = left;
1791            break;
1792          }
1793
1794          if ((*_matching)[right] == INVALID) break;
1795
1796          right = _graph.target((*_matching)[right]);
1797          right_path.push_back(right);
1798          right = _graph.target((*_pred)[right]);
1799          right_path.push_back(right);
1800
1801          right_set.insert(right);
1802
1803        }
1804
1805        if (nca == INVALID) {
1806          if ((*_matching)[left] == INVALID) {
1807            nca = right;
1808            while (left_set.find(nca) == left_set.end()) {
1809              nca = _graph.target((*_matching)[nca]);
1810              right_path.push_back(nca);
1811              nca = _graph.target((*_pred)[nca]);
1812              right_path.push_back(nca);
1813            }
1814          } else {
1815            nca = left;
1816            while (right_set.find(nca) == right_set.end()) {
1817              nca = _graph.target((*_matching)[nca]);
1818              left_path.push_back(nca);
1819              nca = _graph.target((*_pred)[nca]);
1820              left_path.push_back(nca);
1821            }
1822          }
1823        }
1824      }
1825
1826      alternatePath(nca, tree);
1827      Arc prev;
1828
1829      prev = _graph.direct(edge, true);
1830      for (int i = 0; left_path[i] != nca; i += 2) {
1831        _matching->set(left_path[i], prev);
1832        _status->set(left_path[i], MATCHED);
1833        evenToMatched(left_path[i], tree);
1834
1835        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1836        _status->set(left_path[i + 1], MATCHED);
1837        oddToMatched(left_path[i + 1]);
1838      }
1839      _matching->set(nca, prev);
1840
1841      for (int i = 0; right_path[i] != nca; i += 2) {
1842        _status->set(right_path[i], MATCHED);
1843        evenToMatched(right_path[i], tree);
1844
1845        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1846        _status->set(right_path[i + 1], MATCHED);
1847        oddToMatched(right_path[i + 1]);
1848      }
1849
1850      destroyTree(tree);
1851    }
1852
1853    void extractCycle(const Arc &arc) {
1854      Node left = _graph.source(arc);
1855      Node odd = _graph.target((*_matching)[left]);
1856      Arc prev;
1857      while (odd != left) {
1858        Node even = _graph.target((*_matching)[odd]);
1859        prev = (*_matching)[odd];
1860        odd = _graph.target((*_matching)[even]);
1861        _matching->set(even, _graph.oppositeArc(prev));
1862      }
1863      _matching->set(left, arc);
1864
1865      Node right = _graph.target(arc);
1866      int right_tree = _tree_set->find(right);
1867      alternatePath(right, right_tree);
1868      destroyTree(right_tree);
1869      _matching->set(right, _graph.oppositeArc(arc));
1870    }
1871
1872  public:
1873
1874    /// \brief Constructor
1875    ///
1876    /// Constructor.
1877    MaxWeightedPerfectFractionalMatching(const Graph& graph,
1878                                         const WeightMap& weight,
1879                                         bool allow_loops = true)
1880      : _graph(graph), _weight(weight), _matching(0),
1881      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1882      _status(0),  _pred(0),
1883      _tree_set_index(0), _tree_set(0),
1884
1885      _delta2_index(0), _delta2(0),
1886      _delta3_index(0), _delta3(0),
1887
1888      _delta_sum() {}
1889
1890    ~MaxWeightedPerfectFractionalMatching() {
1891      destroyStructures();
1892    }
1893
1894    /// \name Execution Control
1895    /// The simplest way to execute the algorithm is to use the
1896    /// \ref run() member function.
1897
1898    ///@{
1899
1900    /// \brief Initialize the algorithm
1901    ///
1902    /// This function initializes the algorithm.
1903    void init() {
1904      createStructures();
1905
1906      for (NodeIt n(_graph); n != INVALID; ++n) {
1907        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1908      }
1909      for (EdgeIt e(_graph); e != INVALID; ++e) {
1910        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1911      }
1912
1913      for (NodeIt n(_graph); n != INVALID; ++n) {
1914        Value max = - std::numeric_limits<Value>::max();
1915        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1916          if (_graph.target(e) == n && !_allow_loops) continue;
1917          if ((dualScale * _weight[e]) / 2 > max) {
1918            max = (dualScale * _weight[e]) / 2;
1919          }
1920        }
1921        _node_potential->set(n, max);
1922
1923        _tree_set->insert(n);
1924
1925        _matching->set(n, INVALID);
1926        _status->set(n, EVEN);
1927      }
1928
1929      for (EdgeIt e(_graph); e != INVALID; ++e) {
1930        Node left = _graph.u(e);
1931        Node right = _graph.v(e);
1932        if (left == right && !_allow_loops) continue;
1933        _delta3->push(e, ((*_node_potential)[left] +
1934                          (*_node_potential)[right] -
1935                          dualScale * _weight[e]) / 2);
1936      }
1937    }
1938
1939    /// \brief Start the algorithm
1940    ///
1941    /// This function starts the algorithm.
1942    ///
1943    /// \pre \ref init() must be called before using this function.
1944    bool start() {
1945      enum OpType {
1946        D2, D3
1947      };
1948
1949      int unmatched = _node_num;
1950      while (unmatched > 0) {
1951        Value d2 = !_delta2->empty() ?
1952          _delta2->prio() : std::numeric_limits<Value>::max();
1953
1954        Value d3 = !_delta3->empty() ?
1955          _delta3->prio() : std::numeric_limits<Value>::max();
1956
1957        _delta_sum = d3; OpType ot = D3;
1958        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1959
1960        if (_delta_sum == std::numeric_limits<Value>::max()) {
1961          return false;
1962        }
1963
1964        switch (ot) {
1965        case D2:
1966          {
1967            Node n = _delta2->top();
1968            Arc a = (*_pred)[n];
1969            if ((*_matching)[n] == INVALID) {
1970              augmentOnArc(a);
1971              --unmatched;
1972            } else {
1973              Node v = _graph.target((*_matching)[n]);
1974              if ((*_matching)[n] !=
1975                  _graph.oppositeArc((*_matching)[v])) {
1976                extractCycle(a);
1977                --unmatched;
1978              } else {
1979                extendOnArc(a);
1980              }
1981            }
1982          } break;
1983        case D3:
1984          {
1985            Edge e = _delta3->top();
1986
1987            Node left = _graph.u(e);
1988            Node right = _graph.v(e);
1989
1990            int left_tree = _tree_set->find(left);
1991            int right_tree = _tree_set->find(right);
1992
1993            if (left_tree == right_tree) {
1994              cycleOnEdge(e, left_tree);
1995              --unmatched;
1996            } else {
1997              augmentOnEdge(e);
1998              unmatched -= 2;
1999            }
2000          } break;
2001        }
2002      }
2003      return true;
2004    }
2005
2006    /// \brief Run the algorithm.
2007    ///
2008    /// This method runs the \c %MaxWeightedMatching algorithm.
2009    ///
2010    /// \note mwfm.run() is just a shortcut of the following code.
2011    /// \code
2012    ///   mwpfm.init();
2013    ///   mwpfm.start();
2014    /// \endcode
2015    bool run() {
2016      init();
2017      return start();
2018    }
2019
2020    /// @}
2021
2022    /// \name Primal Solution
2023    /// Functions to get the primal solution, i.e. the maximum weighted
2024    /// matching.\n
2025    /// Either \ref run() or \ref start() function should be called before
2026    /// using them.
2027
2028    /// @{
2029
2030    /// \brief Return the weight of the matching.
2031    ///
2032    /// This function returns the weight of the found matching. This
2033    /// value is scaled by \ref primalScale "primal scale".
2034    ///
2035    /// \pre Either run() or start() must be called before using this function.
2036    Value matchingWeight() const {
2037      Value sum = 0;
2038      for (NodeIt n(_graph); n != INVALID; ++n) {
2039        if ((*_matching)[n] != INVALID) {
2040          sum += _weight[(*_matching)[n]];
2041        }
2042      }
2043      return sum * primalScale / 2;
2044    }
2045
2046    /// \brief Return the number of covered nodes in the matching.
2047    ///
2048    /// This function returns the number of covered nodes in the matching.
2049    ///
2050    /// \pre Either run() or start() must be called before using this function.
2051    int matchingSize() const {
2052      int num = 0;
2053      for (NodeIt n(_graph); n != INVALID; ++n) {
2054        if ((*_matching)[n] != INVALID) {
2055          ++num;
2056        }
2057      }
2058      return num;
2059    }
2060
2061    /// \brief Return \c true if the given edge is in the matching.
2062    ///
2063    /// This function returns \c true if the given edge is in the
2064    /// found matching. The result is scaled by \ref primalScale
2065    /// "primal scale".
2066    ///
2067    /// \pre Either run() or start() must be called before using this function.
2068    Value matching(const Edge& edge) const {
2069      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
2070        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
2071        * primalScale / 2;
2072    }
2073
2074    /// \brief Return the fractional matching arc (or edge) incident
2075    /// to the given node.
2076    ///
2077    /// This function returns one of the fractional matching arc (or
2078    /// edge) incident to the given node in the found matching or \c
2079    /// INVALID if the node is not covered by the matching or if the
2080    /// node is on an odd length cycle then it is the successor edge
2081    /// on the cycle.
2082    ///
2083    /// \pre Either run() or start() must be called before using this function.
2084    Arc matching(const Node& node) const {
2085      return (*_matching)[node];
2086    }
2087
2088    /// \brief Return a const reference to the matching map.
2089    ///
2090    /// This function returns a const reference to a node map that stores
2091    /// the matching arc (or edge) incident to each node.
2092    const MatchingMap& matchingMap() const {
2093      return *_matching;
2094    }
2095
2096    /// @}
2097
2098    /// \name Dual Solution
2099    /// Functions to get the dual solution.\n
2100    /// Either \ref run() or \ref start() function should be called before
2101    /// using them.
2102
2103    /// @{
2104
2105    /// \brief Return the value of the dual solution.
2106    ///
2107    /// This function returns the value of the dual solution.
2108    /// It should be equal to the primal value scaled by \ref dualScale
2109    /// "dual scale".
2110    ///
2111    /// \pre Either run() or start() must be called before using this function.
2112    Value dualValue() const {
2113      Value sum = 0;
2114      for (NodeIt n(_graph); n != INVALID; ++n) {
2115        sum += nodeValue(n);
2116      }
2117      return sum;
2118    }
2119
2120    /// \brief Return the dual value (potential) of the given node.
2121    ///
2122    /// This function returns the dual value (potential) of the given node.
2123    ///
2124    /// \pre Either run() or start() must be called before using this function.
2125    Value nodeValue(const Node& n) const {
2126      return (*_node_potential)[n];
2127    }
2128
2129    /// @}
2130
2131  };
2132
2133} //END OF NAMESPACE LEMON
2134
2135#endif //LEMON_FRACTIONAL_MATCHING_H
Note: See TracBrowser for help on using the repository browser.