COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/fractional_matching.h @ 871:86613aa28a0c

Last change on this file since 871:86613aa28a0c was 871:86613aa28a0c, checked in by Balazs Dezso <deba@…>, 10 years ago

Fix documentation issues (#314)

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