COIN-OR::LEMON - Graph Library

source: lemon/lemon/fractional_matching.h @ 954:07ec2b52e53d

Last change on this file since 954:07ec2b52e53d was 951:41d7ac528c3a, checked in by Balazs Dezso <deba@…>, 10 years ago

Uniforming primal scale to 2 (#314)

File size: 62.9 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 MaxFractionalMatching::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 uses priority queues and
636  /// provides \f$O(nm\log n)\f$ time complexity.
637  ///
638  /// The maximum weighted fractional matching is a relaxation of the
639  /// maximum weighted matching problem where the odd set constraints
640  /// are omitted.
641  /// It can be formulated with the following linear program.
642  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
643  /// \f[x_e \ge 0\quad \forall e\in E\f]
644  /// \f[\max \sum_{e\in E}x_ew_e\f]
645  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
646  /// \f$X\f$. The result must be the union of a matching with one
647  /// value edges and a set of odd length cycles with half value edges.
648  ///
649  /// The algorithm calculates an optimal fractional matching and a
650  /// proof of the optimality. The solution of the dual problem can be
651  /// used to check the result of the algorithm. The dual linear
652  /// problem is the following.
653  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
654  /// \f[y_u \ge 0 \quad \forall u \in V\f]
655  /// \f[\min \sum_{u \in V}y_u \f]
656  ///
657  /// The algorithm can be executed with the run() function.
658  /// After it the matching (the primal solution) and the dual solution
659  /// can be obtained using the query functions.
660  ///
661  /// The primal solution is multiplied by
662  /// \ref MaxWeightedFractionalMatching::primalScale "2".
663  /// If the value type is integer, then the dual
664  /// solution is scaled by
665  /// \ref MaxWeightedFractionalMatching::dualScale "4".
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.
693    static const int primalScale = 2;
694
695    /// \brief Scaling factor for dual solution
696    ///
697    /// Scaling factor for dual solution. It is equal to 4 or 1
698    /// according to the value type.
699    static const int dualScale =
700      std::numeric_limits<Value>::is_integer ? 4 : 1;
701
702  private:
703
704    TEMPLATE_GRAPH_TYPEDEFS(Graph);
705
706    typedef typename Graph::template NodeMap<Value> NodePotential;
707
708    const Graph& _graph;
709    const WeightMap& _weight;
710
711    MatchingMap* _matching;
712    NodePotential* _node_potential;
713
714    int _node_num;
715    bool _allow_loops;
716
717    enum Status {
718      EVEN = -1, MATCHED = 0, ODD = 1
719    };
720
721    typedef typename Graph::template NodeMap<Status> StatusMap;
722    StatusMap* _status;
723
724    typedef typename Graph::template NodeMap<Arc> PredMap;
725    PredMap* _pred;
726
727    typedef ExtendFindEnum<IntNodeMap> TreeSet;
728
729    IntNodeMap *_tree_set_index;
730    TreeSet *_tree_set;
731
732    IntNodeMap *_delta1_index;
733    BinHeap<Value, IntNodeMap> *_delta1;
734
735    IntNodeMap *_delta2_index;
736    BinHeap<Value, IntNodeMap> *_delta2;
737
738    IntEdgeMap *_delta3_index;
739    BinHeap<Value, IntEdgeMap> *_delta3;
740
741    Value _delta_sum;
742
743    void createStructures() {
744      _node_num = countNodes(_graph);
745
746      if (!_matching) {
747        _matching = new MatchingMap(_graph);
748      }
749      if (!_node_potential) {
750        _node_potential = new NodePotential(_graph);
751      }
752      if (!_status) {
753        _status = new StatusMap(_graph);
754      }
755      if (!_pred) {
756        _pred = new PredMap(_graph);
757      }
758      if (!_tree_set) {
759        _tree_set_index = new IntNodeMap(_graph);
760        _tree_set = new TreeSet(*_tree_set_index);
761      }
762      if (!_delta1) {
763        _delta1_index = new IntNodeMap(_graph);
764        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
765      }
766      if (!_delta2) {
767        _delta2_index = new IntNodeMap(_graph);
768        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
769      }
770      if (!_delta3) {
771        _delta3_index = new IntEdgeMap(_graph);
772        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
773      }
774    }
775
776    void destroyStructures() {
777      if (_matching) {
778        delete _matching;
779      }
780      if (_node_potential) {
781        delete _node_potential;
782      }
783      if (_status) {
784        delete _status;
785      }
786      if (_pred) {
787        delete _pred;
788      }
789      if (_tree_set) {
790        delete _tree_set_index;
791        delete _tree_set;
792      }
793      if (_delta1) {
794        delete _delta1_index;
795        delete _delta1;
796      }
797      if (_delta2) {
798        delete _delta2_index;
799        delete _delta2;
800      }
801      if (_delta3) {
802        delete _delta3_index;
803        delete _delta3;
804      }
805    }
806
807    void matchedToEven(Node node, int tree) {
808      _tree_set->insert(node, tree);
809      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
810      _delta1->push(node, (*_node_potential)[node]);
811
812      if (_delta2->state(node) == _delta2->IN_HEAP) {
813        _delta2->erase(node);
814      }
815
816      for (InArcIt a(_graph, node); a != INVALID; ++a) {
817        Node v = _graph.source(a);
818        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
819          dualScale * _weight[a];
820        if (node == v) {
821          if (_allow_loops && _graph.direction(a)) {
822            _delta3->push(a, rw / 2);
823          }
824        } else if ((*_status)[v] == EVEN) {
825          _delta3->push(a, rw / 2);
826        } else if ((*_status)[v] == MATCHED) {
827          if (_delta2->state(v) != _delta2->IN_HEAP) {
828            _pred->set(v, a);
829            _delta2->push(v, rw);
830          } else if ((*_delta2)[v] > rw) {
831            _pred->set(v, a);
832            _delta2->decrease(v, rw);
833          }
834        }
835      }
836    }
837
838    void matchedToOdd(Node node, int tree) {
839      _tree_set->insert(node, tree);
840      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
841
842      if (_delta2->state(node) == _delta2->IN_HEAP) {
843        _delta2->erase(node);
844      }
845    }
846
847    void evenToMatched(Node node, int tree) {
848      _delta1->erase(node);
849      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
850      Arc min = INVALID;
851      Value minrw = std::numeric_limits<Value>::max();
852      for (InArcIt a(_graph, node); a != INVALID; ++a) {
853        Node v = _graph.source(a);
854        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
855          dualScale * _weight[a];
856
857        if (node == v) {
858          if (_allow_loops && _graph.direction(a)) {
859            _delta3->erase(a);
860          }
861        } else if ((*_status)[v] == EVEN) {
862          _delta3->erase(a);
863          if (minrw > rw) {
864            min = _graph.oppositeArc(a);
865            minrw = rw;
866          }
867        } else if ((*_status)[v]  == MATCHED) {
868          if ((*_pred)[v] == a) {
869            Arc mina = INVALID;
870            Value minrwa = std::numeric_limits<Value>::max();
871            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
872              Node va = _graph.target(aa);
873              if ((*_status)[va] != EVEN ||
874                  _tree_set->find(va) == tree) continue;
875              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
876                dualScale * _weight[aa];
877              if (minrwa > rwa) {
878                minrwa = rwa;
879                mina = aa;
880              }
881            }
882            if (mina != INVALID) {
883              _pred->set(v, mina);
884              _delta2->increase(v, minrwa);
885            } else {
886              _pred->set(v, INVALID);
887              _delta2->erase(v);
888            }
889          }
890        }
891      }
892      if (min != INVALID) {
893        _pred->set(node, min);
894        _delta2->push(node, minrw);
895      } else {
896        _pred->set(node, INVALID);
897      }
898    }
899
900    void oddToMatched(Node node) {
901      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
902      Arc min = INVALID;
903      Value minrw = std::numeric_limits<Value>::max();
904      for (InArcIt a(_graph, node); a != INVALID; ++a) {
905        Node v = _graph.source(a);
906        if ((*_status)[v] != EVEN) continue;
907        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
908          dualScale * _weight[a];
909
910        if (minrw > rw) {
911          min = _graph.oppositeArc(a);
912          minrw = rw;
913        }
914      }
915      if (min != INVALID) {
916        _pred->set(node, min);
917        _delta2->push(node, minrw);
918      } else {
919        _pred->set(node, INVALID);
920      }
921    }
922
923    void alternatePath(Node even, int tree) {
924      Node odd;
925
926      _status->set(even, MATCHED);
927      evenToMatched(even, tree);
928
929      Arc prev = (*_matching)[even];
930      while (prev != INVALID) {
931        odd = _graph.target(prev);
932        even = _graph.target((*_pred)[odd]);
933        _matching->set(odd, (*_pred)[odd]);
934        _status->set(odd, MATCHED);
935        oddToMatched(odd);
936
937        prev = (*_matching)[even];
938        _status->set(even, MATCHED);
939        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
940        evenToMatched(even, tree);
941      }
942    }
943
944    void destroyTree(int tree) {
945      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
946        if ((*_status)[n] == EVEN) {
947          _status->set(n, MATCHED);
948          evenToMatched(n, tree);
949        } else if ((*_status)[n] == ODD) {
950          _status->set(n, MATCHED);
951          oddToMatched(n);
952        }
953      }
954      _tree_set->eraseClass(tree);
955    }
956
957
958    void unmatchNode(const Node& node) {
959      int tree = _tree_set->find(node);
960
961      alternatePath(node, tree);
962      destroyTree(tree);
963
964      _matching->set(node, INVALID);
965    }
966
967
968    void augmentOnEdge(const Edge& edge) {
969      Node left = _graph.u(edge);
970      int left_tree = _tree_set->find(left);
971
972      alternatePath(left, left_tree);
973      destroyTree(left_tree);
974      _matching->set(left, _graph.direct(edge, true));
975
976      Node right = _graph.v(edge);
977      int right_tree = _tree_set->find(right);
978
979      alternatePath(right, right_tree);
980      destroyTree(right_tree);
981      _matching->set(right, _graph.direct(edge, false));
982    }
983
984    void augmentOnArc(const Arc& arc) {
985      Node left = _graph.source(arc);
986      _status->set(left, MATCHED);
987      _matching->set(left, arc);
988      _pred->set(left, arc);
989
990      Node right = _graph.target(arc);
991      int right_tree = _tree_set->find(right);
992
993      alternatePath(right, right_tree);
994      destroyTree(right_tree);
995      _matching->set(right, _graph.oppositeArc(arc));
996    }
997
998    void extendOnArc(const Arc& arc) {
999      Node base = _graph.target(arc);
1000      int tree = _tree_set->find(base);
1001
1002      Node odd = _graph.source(arc);
1003      _tree_set->insert(odd, tree);
1004      _status->set(odd, ODD);
1005      matchedToOdd(odd, tree);
1006      _pred->set(odd, arc);
1007
1008      Node even = _graph.target((*_matching)[odd]);
1009      _tree_set->insert(even, tree);
1010      _status->set(even, EVEN);
1011      matchedToEven(even, tree);
1012    }
1013
1014    void cycleOnEdge(const Edge& edge, int tree) {
1015      Node nca = INVALID;
1016      std::vector<Node> left_path, right_path;
1017
1018      {
1019        std::set<Node> left_set, right_set;
1020        Node left = _graph.u(edge);
1021        left_path.push_back(left);
1022        left_set.insert(left);
1023
1024        Node right = _graph.v(edge);
1025        right_path.push_back(right);
1026        right_set.insert(right);
1027
1028        while (true) {
1029
1030          if (left_set.find(right) != left_set.end()) {
1031            nca = right;
1032            break;
1033          }
1034
1035          if ((*_matching)[left] == INVALID) break;
1036
1037          left = _graph.target((*_matching)[left]);
1038          left_path.push_back(left);
1039          left = _graph.target((*_pred)[left]);
1040          left_path.push_back(left);
1041
1042          left_set.insert(left);
1043
1044          if (right_set.find(left) != right_set.end()) {
1045            nca = left;
1046            break;
1047          }
1048
1049          if ((*_matching)[right] == INVALID) break;
1050
1051          right = _graph.target((*_matching)[right]);
1052          right_path.push_back(right);
1053          right = _graph.target((*_pred)[right]);
1054          right_path.push_back(right);
1055
1056          right_set.insert(right);
1057
1058        }
1059
1060        if (nca == INVALID) {
1061          if ((*_matching)[left] == INVALID) {
1062            nca = right;
1063            while (left_set.find(nca) == left_set.end()) {
1064              nca = _graph.target((*_matching)[nca]);
1065              right_path.push_back(nca);
1066              nca = _graph.target((*_pred)[nca]);
1067              right_path.push_back(nca);
1068            }
1069          } else {
1070            nca = left;
1071            while (right_set.find(nca) == right_set.end()) {
1072              nca = _graph.target((*_matching)[nca]);
1073              left_path.push_back(nca);
1074              nca = _graph.target((*_pred)[nca]);
1075              left_path.push_back(nca);
1076            }
1077          }
1078        }
1079      }
1080
1081      alternatePath(nca, tree);
1082      Arc prev;
1083
1084      prev = _graph.direct(edge, true);
1085      for (int i = 0; left_path[i] != nca; i += 2) {
1086        _matching->set(left_path[i], prev);
1087        _status->set(left_path[i], MATCHED);
1088        evenToMatched(left_path[i], tree);
1089
1090        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1091        _status->set(left_path[i + 1], MATCHED);
1092        oddToMatched(left_path[i + 1]);
1093      }
1094      _matching->set(nca, prev);
1095
1096      for (int i = 0; right_path[i] != nca; i += 2) {
1097        _status->set(right_path[i], MATCHED);
1098        evenToMatched(right_path[i], tree);
1099
1100        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1101        _status->set(right_path[i + 1], MATCHED);
1102        oddToMatched(right_path[i + 1]);
1103      }
1104
1105      destroyTree(tree);
1106    }
1107
1108    void extractCycle(const Arc &arc) {
1109      Node left = _graph.source(arc);
1110      Node odd = _graph.target((*_matching)[left]);
1111      Arc prev;
1112      while (odd != left) {
1113        Node even = _graph.target((*_matching)[odd]);
1114        prev = (*_matching)[odd];
1115        odd = _graph.target((*_matching)[even]);
1116        _matching->set(even, _graph.oppositeArc(prev));
1117      }
1118      _matching->set(left, arc);
1119
1120      Node right = _graph.target(arc);
1121      int right_tree = _tree_set->find(right);
1122      alternatePath(right, right_tree);
1123      destroyTree(right_tree);
1124      _matching->set(right, _graph.oppositeArc(arc));
1125    }
1126
1127  public:
1128
1129    /// \brief Constructor
1130    ///
1131    /// Constructor.
1132    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
1133                                  bool allow_loops = true)
1134      : _graph(graph), _weight(weight), _matching(0),
1135      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1136      _status(0),  _pred(0),
1137      _tree_set_index(0), _tree_set(0),
1138
1139      _delta1_index(0), _delta1(0),
1140      _delta2_index(0), _delta2(0),
1141      _delta3_index(0), _delta3(0),
1142
1143      _delta_sum() {}
1144
1145    ~MaxWeightedFractionalMatching() {
1146      destroyStructures();
1147    }
1148
1149    /// \name Execution Control
1150    /// The simplest way to execute the algorithm is to use the
1151    /// \ref run() member function.
1152
1153    ///@{
1154
1155    /// \brief Initialize the algorithm
1156    ///
1157    /// This function initializes the algorithm.
1158    void init() {
1159      createStructures();
1160
1161      for (NodeIt n(_graph); n != INVALID; ++n) {
1162        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1163        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1164      }
1165      for (EdgeIt e(_graph); e != INVALID; ++e) {
1166        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1167      }
1168
1169      for (NodeIt n(_graph); n != INVALID; ++n) {
1170        Value max = 0;
1171        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1172          if (_graph.target(e) == n && !_allow_loops) continue;
1173          if ((dualScale * _weight[e]) / 2 > max) {
1174            max = (dualScale * _weight[e]) / 2;
1175          }
1176        }
1177        _node_potential->set(n, max);
1178        _delta1->push(n, max);
1179
1180        _tree_set->insert(n);
1181
1182        _matching->set(n, INVALID);
1183        _status->set(n, EVEN);
1184      }
1185
1186      for (EdgeIt e(_graph); e != INVALID; ++e) {
1187        Node left = _graph.u(e);
1188        Node right = _graph.v(e);
1189        if (left == right && !_allow_loops) continue;
1190        _delta3->push(e, ((*_node_potential)[left] +
1191                          (*_node_potential)[right] -
1192                          dualScale * _weight[e]) / 2);
1193      }
1194    }
1195
1196    /// \brief Start the algorithm
1197    ///
1198    /// This function starts the algorithm.
1199    ///
1200    /// \pre \ref init() must be called before using this function.
1201    void start() {
1202      enum OpType {
1203        D1, D2, D3
1204      };
1205
1206      int unmatched = _node_num;
1207      while (unmatched > 0) {
1208        Value d1 = !_delta1->empty() ?
1209          _delta1->prio() : std::numeric_limits<Value>::max();
1210
1211        Value d2 = !_delta2->empty() ?
1212          _delta2->prio() : std::numeric_limits<Value>::max();
1213
1214        Value d3 = !_delta3->empty() ?
1215          _delta3->prio() : std::numeric_limits<Value>::max();
1216
1217        _delta_sum = d3; OpType ot = D3;
1218        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1219        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1220
1221        switch (ot) {
1222        case D1:
1223          {
1224            Node n = _delta1->top();
1225            unmatchNode(n);
1226            --unmatched;
1227          }
1228          break;
1229        case D2:
1230          {
1231            Node n = _delta2->top();
1232            Arc a = (*_pred)[n];
1233            if ((*_matching)[n] == INVALID) {
1234              augmentOnArc(a);
1235              --unmatched;
1236            } else {
1237              Node v = _graph.target((*_matching)[n]);
1238              if ((*_matching)[n] !=
1239                  _graph.oppositeArc((*_matching)[v])) {
1240                extractCycle(a);
1241                --unmatched;
1242              } else {
1243                extendOnArc(a);
1244              }
1245            }
1246          } break;
1247        case D3:
1248          {
1249            Edge e = _delta3->top();
1250
1251            Node left = _graph.u(e);
1252            Node right = _graph.v(e);
1253
1254            int left_tree = _tree_set->find(left);
1255            int right_tree = _tree_set->find(right);
1256
1257            if (left_tree == right_tree) {
1258              cycleOnEdge(e, left_tree);
1259              --unmatched;
1260            } else {
1261              augmentOnEdge(e);
1262              unmatched -= 2;
1263            }
1264          } break;
1265        }
1266      }
1267    }
1268
1269    /// \brief Run the algorithm.
1270    ///
1271    /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
1272    ///
1273    /// \note mwfm.run() is just a shortcut of the following code.
1274    /// \code
1275    ///   mwfm.init();
1276    ///   mwfm.start();
1277    /// \endcode
1278    void run() {
1279      init();
1280      start();
1281    }
1282
1283    /// @}
1284
1285    /// \name Primal Solution
1286    /// Functions to get the primal solution, i.e. the maximum weighted
1287    /// matching.\n
1288    /// Either \ref run() or \ref start() function should be called before
1289    /// using them.
1290
1291    /// @{
1292
1293    /// \brief Return the weight of the matching.
1294    ///
1295    /// This function returns the weight of the found matching. This
1296    /// value is scaled by \ref primalScale "primal scale".
1297    ///
1298    /// \pre Either run() or start() must be called before using this function.
1299    Value matchingWeight() const {
1300      Value sum = 0;
1301      for (NodeIt n(_graph); n != INVALID; ++n) {
1302        if ((*_matching)[n] != INVALID) {
1303          sum += _weight[(*_matching)[n]];
1304        }
1305      }
1306      return sum * primalScale / 2;
1307    }
1308
1309    /// \brief Return the number of covered nodes in the matching.
1310    ///
1311    /// This function returns the number of covered nodes in the matching.
1312    ///
1313    /// \pre Either run() or start() must be called before using this function.
1314    int matchingSize() const {
1315      int num = 0;
1316      for (NodeIt n(_graph); n != INVALID; ++n) {
1317        if ((*_matching)[n] != INVALID) {
1318          ++num;
1319        }
1320      }
1321      return num;
1322    }
1323
1324    /// \brief Return \c true if the given edge is in the matching.
1325    ///
1326    /// This function returns \c true if the given edge is in the
1327    /// found matching. The result is scaled by \ref primalScale
1328    /// "primal scale".
1329    ///
1330    /// \pre Either run() or start() must be called before using this function.
1331    int matching(const Edge& edge) const {
1332      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1333        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
1334    }
1335
1336    /// \brief Return the fractional matching arc (or edge) incident
1337    /// to the given node.
1338    ///
1339    /// This function returns one of the fractional matching arc (or
1340    /// edge) incident to the given node in the found matching or \c
1341    /// INVALID if the node is not covered by the matching or if the
1342    /// node is on an odd length cycle then it is the successor edge
1343    /// on the cycle.
1344    ///
1345    /// \pre Either run() or start() must be called before using this function.
1346    Arc matching(const Node& node) const {
1347      return (*_matching)[node];
1348    }
1349
1350    /// \brief Return a const reference to the matching map.
1351    ///
1352    /// This function returns a const reference to a node map that stores
1353    /// the matching arc (or edge) incident to each node.
1354    const MatchingMap& matchingMap() const {
1355      return *_matching;
1356    }
1357
1358    /// @}
1359
1360    /// \name Dual Solution
1361    /// Functions to get the dual solution.\n
1362    /// Either \ref run() or \ref start() function should be called before
1363    /// using them.
1364
1365    /// @{
1366
1367    /// \brief Return the value of the dual solution.
1368    ///
1369    /// This function returns the value of the dual solution.
1370    /// It should be equal to the primal value scaled by \ref dualScale
1371    /// "dual scale".
1372    ///
1373    /// \pre Either run() or start() must be called before using this function.
1374    Value dualValue() const {
1375      Value sum = 0;
1376      for (NodeIt n(_graph); n != INVALID; ++n) {
1377        sum += nodeValue(n);
1378      }
1379      return sum;
1380    }
1381
1382    /// \brief Return the dual value (potential) of the given node.
1383    ///
1384    /// This function returns the dual value (potential) of the given node.
1385    ///
1386    /// \pre Either run() or start() must be called before using this function.
1387    Value nodeValue(const Node& n) const {
1388      return (*_node_potential)[n];
1389    }
1390
1391    /// @}
1392
1393  };
1394
1395  /// \ingroup matching
1396  ///
1397  /// \brief Weighted fractional perfect matching in general graphs
1398  ///
1399  /// This class provides an efficient implementation of fractional
1400  /// matching algorithm. The implementation uses priority queues and
1401  /// provides \f$O(nm\log n)\f$ time complexity.
1402  ///
1403  /// The maximum weighted fractional perfect matching is a relaxation
1404  /// of the maximum weighted perfect matching problem where the odd
1405  /// set constraints are omitted.
1406  /// It can be formulated with the following linear program.
1407  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1408  /// \f[x_e \ge 0\quad \forall e\in E\f]
1409  /// \f[\max \sum_{e\in E}x_ew_e\f]
1410  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1411  /// \f$X\f$. The result must be the union of a matching with one
1412  /// value edges and a set of odd length cycles with half value edges.
1413  ///
1414  /// The algorithm calculates an optimal fractional matching and a
1415  /// proof of the optimality. The solution of the dual problem can be
1416  /// used to check the result of the algorithm. The dual linear
1417  /// problem is the following.
1418  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1419  /// \f[\min \sum_{u \in V}y_u \f]
1420  ///
1421  /// The algorithm can be executed with the run() function.
1422  /// After it the matching (the primal solution) and the dual solution
1423  /// can be obtained using the query functions.
1424  ///
1425  /// The primal solution is multiplied by
1426  /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
1427  /// If the value type is integer, then the dual
1428  /// solution is scaled by
1429  /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
1430  ///
1431  /// \tparam GR The undirected graph type the algorithm runs on.
1432  /// \tparam WM The type edge weight map. The default type is
1433  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1434#ifdef DOXYGEN
1435  template <typename GR, typename WM>
1436#else
1437  template <typename GR,
1438            typename WM = typename GR::template EdgeMap<int> >
1439#endif
1440  class MaxWeightedPerfectFractionalMatching {
1441  public:
1442
1443    /// The graph type of the algorithm
1444    typedef GR Graph;
1445    /// The type of the edge weight map
1446    typedef WM WeightMap;
1447    /// The value type of the edge weights
1448    typedef typename WeightMap::Value Value;
1449
1450    /// The type of the matching map
1451    typedef typename Graph::template NodeMap<typename Graph::Arc>
1452    MatchingMap;
1453
1454    /// \brief Scaling factor for primal solution
1455    ///
1456    /// Scaling factor for primal solution.
1457    static const int primalScale = 2;
1458
1459    /// \brief Scaling factor for dual solution
1460    ///
1461    /// Scaling factor for dual solution. It is equal to 4 or 1
1462    /// according to the value type.
1463    static const int dualScale =
1464      std::numeric_limits<Value>::is_integer ? 4 : 1;
1465
1466  private:
1467
1468    TEMPLATE_GRAPH_TYPEDEFS(Graph);
1469
1470    typedef typename Graph::template NodeMap<Value> NodePotential;
1471
1472    const Graph& _graph;
1473    const WeightMap& _weight;
1474
1475    MatchingMap* _matching;
1476    NodePotential* _node_potential;
1477
1478    int _node_num;
1479    bool _allow_loops;
1480
1481    enum Status {
1482      EVEN = -1, MATCHED = 0, ODD = 1
1483    };
1484
1485    typedef typename Graph::template NodeMap<Status> StatusMap;
1486    StatusMap* _status;
1487
1488    typedef typename Graph::template NodeMap<Arc> PredMap;
1489    PredMap* _pred;
1490
1491    typedef ExtendFindEnum<IntNodeMap> TreeSet;
1492
1493    IntNodeMap *_tree_set_index;
1494    TreeSet *_tree_set;
1495
1496    IntNodeMap *_delta2_index;
1497    BinHeap<Value, IntNodeMap> *_delta2;
1498
1499    IntEdgeMap *_delta3_index;
1500    BinHeap<Value, IntEdgeMap> *_delta3;
1501
1502    Value _delta_sum;
1503
1504    void createStructures() {
1505      _node_num = countNodes(_graph);
1506
1507      if (!_matching) {
1508        _matching = new MatchingMap(_graph);
1509      }
1510      if (!_node_potential) {
1511        _node_potential = new NodePotential(_graph);
1512      }
1513      if (!_status) {
1514        _status = new StatusMap(_graph);
1515      }
1516      if (!_pred) {
1517        _pred = new PredMap(_graph);
1518      }
1519      if (!_tree_set) {
1520        _tree_set_index = new IntNodeMap(_graph);
1521        _tree_set = new TreeSet(*_tree_set_index);
1522      }
1523      if (!_delta2) {
1524        _delta2_index = new IntNodeMap(_graph);
1525        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1526      }
1527      if (!_delta3) {
1528        _delta3_index = new IntEdgeMap(_graph);
1529        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1530      }
1531    }
1532
1533    void destroyStructures() {
1534      if (_matching) {
1535        delete _matching;
1536      }
1537      if (_node_potential) {
1538        delete _node_potential;
1539      }
1540      if (_status) {
1541        delete _status;
1542      }
1543      if (_pred) {
1544        delete _pred;
1545      }
1546      if (_tree_set) {
1547        delete _tree_set_index;
1548        delete _tree_set;
1549      }
1550      if (_delta2) {
1551        delete _delta2_index;
1552        delete _delta2;
1553      }
1554      if (_delta3) {
1555        delete _delta3_index;
1556        delete _delta3;
1557      }
1558    }
1559
1560    void matchedToEven(Node node, int tree) {
1561      _tree_set->insert(node, tree);
1562      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1563
1564      if (_delta2->state(node) == _delta2->IN_HEAP) {
1565        _delta2->erase(node);
1566      }
1567
1568      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1569        Node v = _graph.source(a);
1570        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1571          dualScale * _weight[a];
1572        if (node == v) {
1573          if (_allow_loops && _graph.direction(a)) {
1574            _delta3->push(a, rw / 2);
1575          }
1576        } else if ((*_status)[v] == EVEN) {
1577          _delta3->push(a, rw / 2);
1578        } else if ((*_status)[v] == MATCHED) {
1579          if (_delta2->state(v) != _delta2->IN_HEAP) {
1580            _pred->set(v, a);
1581            _delta2->push(v, rw);
1582          } else if ((*_delta2)[v] > rw) {
1583            _pred->set(v, a);
1584            _delta2->decrease(v, rw);
1585          }
1586        }
1587      }
1588    }
1589
1590    void matchedToOdd(Node node, int tree) {
1591      _tree_set->insert(node, tree);
1592      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1593
1594      if (_delta2->state(node) == _delta2->IN_HEAP) {
1595        _delta2->erase(node);
1596      }
1597    }
1598
1599    void evenToMatched(Node node, int tree) {
1600      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1601      Arc min = INVALID;
1602      Value minrw = std::numeric_limits<Value>::max();
1603      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1604        Node v = _graph.source(a);
1605        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1606          dualScale * _weight[a];
1607
1608        if (node == v) {
1609          if (_allow_loops && _graph.direction(a)) {
1610            _delta3->erase(a);
1611          }
1612        } else if ((*_status)[v] == EVEN) {
1613          _delta3->erase(a);
1614          if (minrw > rw) {
1615            min = _graph.oppositeArc(a);
1616            minrw = rw;
1617          }
1618        } else if ((*_status)[v]  == MATCHED) {
1619          if ((*_pred)[v] == a) {
1620            Arc mina = INVALID;
1621            Value minrwa = std::numeric_limits<Value>::max();
1622            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1623              Node va = _graph.target(aa);
1624              if ((*_status)[va] != EVEN ||
1625                  _tree_set->find(va) == tree) continue;
1626              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1627                dualScale * _weight[aa];
1628              if (minrwa > rwa) {
1629                minrwa = rwa;
1630                mina = aa;
1631              }
1632            }
1633            if (mina != INVALID) {
1634              _pred->set(v, mina);
1635              _delta2->increase(v, minrwa);
1636            } else {
1637              _pred->set(v, INVALID);
1638              _delta2->erase(v);
1639            }
1640          }
1641        }
1642      }
1643      if (min != INVALID) {
1644        _pred->set(node, min);
1645        _delta2->push(node, minrw);
1646      } else {
1647        _pred->set(node, INVALID);
1648      }
1649    }
1650
1651    void oddToMatched(Node node) {
1652      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1653      Arc min = INVALID;
1654      Value minrw = std::numeric_limits<Value>::max();
1655      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1656        Node v = _graph.source(a);
1657        if ((*_status)[v] != EVEN) continue;
1658        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1659          dualScale * _weight[a];
1660
1661        if (minrw > rw) {
1662          min = _graph.oppositeArc(a);
1663          minrw = rw;
1664        }
1665      }
1666      if (min != INVALID) {
1667        _pred->set(node, min);
1668        _delta2->push(node, minrw);
1669      } else {
1670        _pred->set(node, INVALID);
1671      }
1672    }
1673
1674    void alternatePath(Node even, int tree) {
1675      Node odd;
1676
1677      _status->set(even, MATCHED);
1678      evenToMatched(even, tree);
1679
1680      Arc prev = (*_matching)[even];
1681      while (prev != INVALID) {
1682        odd = _graph.target(prev);
1683        even = _graph.target((*_pred)[odd]);
1684        _matching->set(odd, (*_pred)[odd]);
1685        _status->set(odd, MATCHED);
1686        oddToMatched(odd);
1687
1688        prev = (*_matching)[even];
1689        _status->set(even, MATCHED);
1690        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1691        evenToMatched(even, tree);
1692      }
1693    }
1694
1695    void destroyTree(int tree) {
1696      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1697        if ((*_status)[n] == EVEN) {
1698          _status->set(n, MATCHED);
1699          evenToMatched(n, tree);
1700        } else if ((*_status)[n] == ODD) {
1701          _status->set(n, MATCHED);
1702          oddToMatched(n);
1703        }
1704      }
1705      _tree_set->eraseClass(tree);
1706    }
1707
1708    void augmentOnEdge(const Edge& edge) {
1709      Node left = _graph.u(edge);
1710      int left_tree = _tree_set->find(left);
1711
1712      alternatePath(left, left_tree);
1713      destroyTree(left_tree);
1714      _matching->set(left, _graph.direct(edge, true));
1715
1716      Node right = _graph.v(edge);
1717      int right_tree = _tree_set->find(right);
1718
1719      alternatePath(right, right_tree);
1720      destroyTree(right_tree);
1721      _matching->set(right, _graph.direct(edge, false));
1722    }
1723
1724    void augmentOnArc(const Arc& arc) {
1725      Node left = _graph.source(arc);
1726      _status->set(left, MATCHED);
1727      _matching->set(left, arc);
1728      _pred->set(left, arc);
1729
1730      Node right = _graph.target(arc);
1731      int right_tree = _tree_set->find(right);
1732
1733      alternatePath(right, right_tree);
1734      destroyTree(right_tree);
1735      _matching->set(right, _graph.oppositeArc(arc));
1736    }
1737
1738    void extendOnArc(const Arc& arc) {
1739      Node base = _graph.target(arc);
1740      int tree = _tree_set->find(base);
1741
1742      Node odd = _graph.source(arc);
1743      _tree_set->insert(odd, tree);
1744      _status->set(odd, ODD);
1745      matchedToOdd(odd, tree);
1746      _pred->set(odd, arc);
1747
1748      Node even = _graph.target((*_matching)[odd]);
1749      _tree_set->insert(even, tree);
1750      _status->set(even, EVEN);
1751      matchedToEven(even, tree);
1752    }
1753
1754    void cycleOnEdge(const Edge& edge, int tree) {
1755      Node nca = INVALID;
1756      std::vector<Node> left_path, right_path;
1757
1758      {
1759        std::set<Node> left_set, right_set;
1760        Node left = _graph.u(edge);
1761        left_path.push_back(left);
1762        left_set.insert(left);
1763
1764        Node right = _graph.v(edge);
1765        right_path.push_back(right);
1766        right_set.insert(right);
1767
1768        while (true) {
1769
1770          if (left_set.find(right) != left_set.end()) {
1771            nca = right;
1772            break;
1773          }
1774
1775          if ((*_matching)[left] == INVALID) break;
1776
1777          left = _graph.target((*_matching)[left]);
1778          left_path.push_back(left);
1779          left = _graph.target((*_pred)[left]);
1780          left_path.push_back(left);
1781
1782          left_set.insert(left);
1783
1784          if (right_set.find(left) != right_set.end()) {
1785            nca = left;
1786            break;
1787          }
1788
1789          if ((*_matching)[right] == INVALID) break;
1790
1791          right = _graph.target((*_matching)[right]);
1792          right_path.push_back(right);
1793          right = _graph.target((*_pred)[right]);
1794          right_path.push_back(right);
1795
1796          right_set.insert(right);
1797
1798        }
1799
1800        if (nca == INVALID) {
1801          if ((*_matching)[left] == INVALID) {
1802            nca = right;
1803            while (left_set.find(nca) == left_set.end()) {
1804              nca = _graph.target((*_matching)[nca]);
1805              right_path.push_back(nca);
1806              nca = _graph.target((*_pred)[nca]);
1807              right_path.push_back(nca);
1808            }
1809          } else {
1810            nca = left;
1811            while (right_set.find(nca) == right_set.end()) {
1812              nca = _graph.target((*_matching)[nca]);
1813              left_path.push_back(nca);
1814              nca = _graph.target((*_pred)[nca]);
1815              left_path.push_back(nca);
1816            }
1817          }
1818        }
1819      }
1820
1821      alternatePath(nca, tree);
1822      Arc prev;
1823
1824      prev = _graph.direct(edge, true);
1825      for (int i = 0; left_path[i] != nca; i += 2) {
1826        _matching->set(left_path[i], prev);
1827        _status->set(left_path[i], MATCHED);
1828        evenToMatched(left_path[i], tree);
1829
1830        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1831        _status->set(left_path[i + 1], MATCHED);
1832        oddToMatched(left_path[i + 1]);
1833      }
1834      _matching->set(nca, prev);
1835
1836      for (int i = 0; right_path[i] != nca; i += 2) {
1837        _status->set(right_path[i], MATCHED);
1838        evenToMatched(right_path[i], tree);
1839
1840        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1841        _status->set(right_path[i + 1], MATCHED);
1842        oddToMatched(right_path[i + 1]);
1843      }
1844
1845      destroyTree(tree);
1846    }
1847
1848    void extractCycle(const Arc &arc) {
1849      Node left = _graph.source(arc);
1850      Node odd = _graph.target((*_matching)[left]);
1851      Arc prev;
1852      while (odd != left) {
1853        Node even = _graph.target((*_matching)[odd]);
1854        prev = (*_matching)[odd];
1855        odd = _graph.target((*_matching)[even]);
1856        _matching->set(even, _graph.oppositeArc(prev));
1857      }
1858      _matching->set(left, arc);
1859
1860      Node right = _graph.target(arc);
1861      int right_tree = _tree_set->find(right);
1862      alternatePath(right, right_tree);
1863      destroyTree(right_tree);
1864      _matching->set(right, _graph.oppositeArc(arc));
1865    }
1866
1867  public:
1868
1869    /// \brief Constructor
1870    ///
1871    /// Constructor.
1872    MaxWeightedPerfectFractionalMatching(const Graph& graph,
1873                                         const WeightMap& weight,
1874                                         bool allow_loops = true)
1875      : _graph(graph), _weight(weight), _matching(0),
1876      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1877      _status(0),  _pred(0),
1878      _tree_set_index(0), _tree_set(0),
1879
1880      _delta2_index(0), _delta2(0),
1881      _delta3_index(0), _delta3(0),
1882
1883      _delta_sum() {}
1884
1885    ~MaxWeightedPerfectFractionalMatching() {
1886      destroyStructures();
1887    }
1888
1889    /// \name Execution Control
1890    /// The simplest way to execute the algorithm is to use the
1891    /// \ref run() member function.
1892
1893    ///@{
1894
1895    /// \brief Initialize the algorithm
1896    ///
1897    /// This function initializes the algorithm.
1898    void init() {
1899      createStructures();
1900
1901      for (NodeIt n(_graph); n != INVALID; ++n) {
1902        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1903      }
1904      for (EdgeIt e(_graph); e != INVALID; ++e) {
1905        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1906      }
1907
1908      for (NodeIt n(_graph); n != INVALID; ++n) {
1909        Value max = - std::numeric_limits<Value>::max();
1910        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1911          if (_graph.target(e) == n && !_allow_loops) continue;
1912          if ((dualScale * _weight[e]) / 2 > max) {
1913            max = (dualScale * _weight[e]) / 2;
1914          }
1915        }
1916        _node_potential->set(n, max);
1917
1918        _tree_set->insert(n);
1919
1920        _matching->set(n, INVALID);
1921        _status->set(n, EVEN);
1922      }
1923
1924      for (EdgeIt e(_graph); e != INVALID; ++e) {
1925        Node left = _graph.u(e);
1926        Node right = _graph.v(e);
1927        if (left == right && !_allow_loops) continue;
1928        _delta3->push(e, ((*_node_potential)[left] +
1929                          (*_node_potential)[right] -
1930                          dualScale * _weight[e]) / 2);
1931      }
1932    }
1933
1934    /// \brief Start the algorithm
1935    ///
1936    /// This function starts the algorithm.
1937    ///
1938    /// \pre \ref init() must be called before using this function.
1939    bool start() {
1940      enum OpType {
1941        D2, D3
1942      };
1943
1944      int unmatched = _node_num;
1945      while (unmatched > 0) {
1946        Value d2 = !_delta2->empty() ?
1947          _delta2->prio() : std::numeric_limits<Value>::max();
1948
1949        Value d3 = !_delta3->empty() ?
1950          _delta3->prio() : std::numeric_limits<Value>::max();
1951
1952        _delta_sum = d3; OpType ot = D3;
1953        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1954
1955        if (_delta_sum == std::numeric_limits<Value>::max()) {
1956          return false;
1957        }
1958
1959        switch (ot) {
1960        case D2:
1961          {
1962            Node n = _delta2->top();
1963            Arc a = (*_pred)[n];
1964            if ((*_matching)[n] == INVALID) {
1965              augmentOnArc(a);
1966              --unmatched;
1967            } else {
1968              Node v = _graph.target((*_matching)[n]);
1969              if ((*_matching)[n] !=
1970                  _graph.oppositeArc((*_matching)[v])) {
1971                extractCycle(a);
1972                --unmatched;
1973              } else {
1974                extendOnArc(a);
1975              }
1976            }
1977          } break;
1978        case D3:
1979          {
1980            Edge e = _delta3->top();
1981
1982            Node left = _graph.u(e);
1983            Node right = _graph.v(e);
1984
1985            int left_tree = _tree_set->find(left);
1986            int right_tree = _tree_set->find(right);
1987
1988            if (left_tree == right_tree) {
1989              cycleOnEdge(e, left_tree);
1990              --unmatched;
1991            } else {
1992              augmentOnEdge(e);
1993              unmatched -= 2;
1994            }
1995          } break;
1996        }
1997      }
1998      return true;
1999    }
2000
2001    /// \brief Run the algorithm.
2002    ///
2003    /// This method runs the \c %MaxWeightedPerfectFractionalMatching
2004    /// algorithm.
2005    ///
2006    /// \note mwfm.run() is just a shortcut of the following code.
2007    /// \code
2008    ///   mwpfm.init();
2009    ///   mwpfm.start();
2010    /// \endcode
2011    bool run() {
2012      init();
2013      return start();
2014    }
2015
2016    /// @}
2017
2018    /// \name Primal Solution
2019    /// Functions to get the primal solution, i.e. the maximum weighted
2020    /// matching.\n
2021    /// Either \ref run() or \ref start() function should be called before
2022    /// using them.
2023
2024    /// @{
2025
2026    /// \brief Return the weight of the matching.
2027    ///
2028    /// This function returns the weight of the found matching. This
2029    /// value is scaled by \ref primalScale "primal scale".
2030    ///
2031    /// \pre Either run() or start() must be called before using this function.
2032    Value matchingWeight() const {
2033      Value sum = 0;
2034      for (NodeIt n(_graph); n != INVALID; ++n) {
2035        if ((*_matching)[n] != INVALID) {
2036          sum += _weight[(*_matching)[n]];
2037        }
2038      }
2039      return sum * primalScale / 2;
2040    }
2041
2042    /// \brief Return the number of covered nodes in the matching.
2043    ///
2044    /// This function returns the number of covered nodes in the matching.
2045    ///
2046    /// \pre Either run() or start() must be called before using this function.
2047    int matchingSize() const {
2048      int num = 0;
2049      for (NodeIt n(_graph); n != INVALID; ++n) {
2050        if ((*_matching)[n] != INVALID) {
2051          ++num;
2052        }
2053      }
2054      return num;
2055    }
2056
2057    /// \brief Return \c true if the given edge is in the matching.
2058    ///
2059    /// This function returns \c true if the given edge is in the
2060    /// found matching. The result is scaled by \ref primalScale
2061    /// "primal scale".
2062    ///
2063    /// \pre Either run() or start() must be called before using this function.
2064    int matching(const Edge& edge) const {
2065      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
2066        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
2067    }
2068
2069    /// \brief Return the fractional matching arc (or edge) incident
2070    /// to the given node.
2071    ///
2072    /// This function returns one of the fractional matching arc (or
2073    /// edge) incident to the given node in the found matching or \c
2074    /// INVALID if the node is not covered by the matching or if the
2075    /// node is on an odd length cycle then it is the successor edge
2076    /// on the cycle.
2077    ///
2078    /// \pre Either run() or start() must be called before using this function.
2079    Arc matching(const Node& node) const {
2080      return (*_matching)[node];
2081    }
2082
2083    /// \brief Return a const reference to the matching map.
2084    ///
2085    /// This function returns a const reference to a node map that stores
2086    /// the matching arc (or edge) incident to each node.
2087    const MatchingMap& matchingMap() const {
2088      return *_matching;
2089    }
2090
2091    /// @}
2092
2093    /// \name Dual Solution
2094    /// Functions to get the dual solution.\n
2095    /// Either \ref run() or \ref start() function should be called before
2096    /// using them.
2097
2098    /// @{
2099
2100    /// \brief Return the value of the dual solution.
2101    ///
2102    /// This function returns the value of the dual solution.
2103    /// It should be equal to the primal value scaled by \ref dualScale
2104    /// "dual scale".
2105    ///
2106    /// \pre Either run() or start() must be called before using this function.
2107    Value dualValue() const {
2108      Value sum = 0;
2109      for (NodeIt n(_graph); n != INVALID; ++n) {
2110        sum += nodeValue(n);
2111      }
2112      return sum;
2113    }
2114
2115    /// \brief Return the dual value (potential) of the given node.
2116    ///
2117    /// This function returns the dual value (potential) of the given node.
2118    ///
2119    /// \pre Either run() or start() must be called before using this function.
2120    Value nodeValue(const Node& n) const {
2121      return (*_node_potential)[n];
2122    }
2123
2124    /// @}
2125
2126  };
2127
2128} //END OF NAMESPACE LEMON
2129
2130#endif //LEMON_FRACTIONAL_MATCHING_H
Note: See TracBrowser for help on using the repository browser.