COIN-OR::LEMON - Graph Library

source: lemon/lemon/fractional_matching.h @ 1300:62dba6c90f35

Last change on this file since 1300:62dba6c90f35 was 1270:dceba191c00d, checked in by Alpar Juttner <alpar@…>, 11 years ago

Apply unify-sources.sh to the source tree

File size: 63.1 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-2013
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 lemon::MaxFractionalMatchingDefaultTraits
127    /// "traits 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      _delta1->clear();
1170      _delta2->clear();
1171      _delta3->clear();
1172      _tree_set->clear();
1173
1174      for (NodeIt n(_graph); n != INVALID; ++n) {
1175        Value max = 0;
1176        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1177          if (_graph.target(e) == n && !_allow_loops) continue;
1178          if ((dualScale * _weight[e]) / 2 > max) {
1179            max = (dualScale * _weight[e]) / 2;
1180          }
1181        }
1182        _node_potential->set(n, max);
1183        _delta1->push(n, max);
1184
1185        _tree_set->insert(n);
1186
1187        _matching->set(n, INVALID);
1188        _status->set(n, EVEN);
1189      }
1190
1191      for (EdgeIt e(_graph); e != INVALID; ++e) {
1192        Node left = _graph.u(e);
1193        Node right = _graph.v(e);
1194        if (left == right && !_allow_loops) continue;
1195        _delta3->push(e, ((*_node_potential)[left] +
1196                          (*_node_potential)[right] -
1197                          dualScale * _weight[e]) / 2);
1198      }
1199    }
1200
1201    /// \brief Start the algorithm
1202    ///
1203    /// This function starts the algorithm.
1204    ///
1205    /// \pre \ref init() must be called before using this function.
1206    void start() {
1207      enum OpType {
1208        D1, D2, D3
1209      };
1210
1211      int unmatched = _node_num;
1212      while (unmatched > 0) {
1213        Value d1 = !_delta1->empty() ?
1214          _delta1->prio() : std::numeric_limits<Value>::max();
1215
1216        Value d2 = !_delta2->empty() ?
1217          _delta2->prio() : std::numeric_limits<Value>::max();
1218
1219        Value d3 = !_delta3->empty() ?
1220          _delta3->prio() : std::numeric_limits<Value>::max();
1221
1222        _delta_sum = d3; OpType ot = D3;
1223        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1224        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1225
1226        switch (ot) {
1227        case D1:
1228          {
1229            Node n = _delta1->top();
1230            unmatchNode(n);
1231            --unmatched;
1232          }
1233          break;
1234        case D2:
1235          {
1236            Node n = _delta2->top();
1237            Arc a = (*_pred)[n];
1238            if ((*_matching)[n] == INVALID) {
1239              augmentOnArc(a);
1240              --unmatched;
1241            } else {
1242              Node v = _graph.target((*_matching)[n]);
1243              if ((*_matching)[n] !=
1244                  _graph.oppositeArc((*_matching)[v])) {
1245                extractCycle(a);
1246                --unmatched;
1247              } else {
1248                extendOnArc(a);
1249              }
1250            }
1251          } break;
1252        case D3:
1253          {
1254            Edge e = _delta3->top();
1255
1256            Node left = _graph.u(e);
1257            Node right = _graph.v(e);
1258
1259            int left_tree = _tree_set->find(left);
1260            int right_tree = _tree_set->find(right);
1261
1262            if (left_tree == right_tree) {
1263              cycleOnEdge(e, left_tree);
1264              --unmatched;
1265            } else {
1266              augmentOnEdge(e);
1267              unmatched -= 2;
1268            }
1269          } break;
1270        }
1271      }
1272    }
1273
1274    /// \brief Run the algorithm.
1275    ///
1276    /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
1277    ///
1278    /// \note mwfm.run() is just a shortcut of the following code.
1279    /// \code
1280    ///   mwfm.init();
1281    ///   mwfm.start();
1282    /// \endcode
1283    void run() {
1284      init();
1285      start();
1286    }
1287
1288    /// @}
1289
1290    /// \name Primal Solution
1291    /// Functions to get the primal solution, i.e. the maximum weighted
1292    /// matching.\n
1293    /// Either \ref run() or \ref start() function should be called before
1294    /// using them.
1295
1296    /// @{
1297
1298    /// \brief Return the weight of the matching.
1299    ///
1300    /// This function returns the weight of the found matching. This
1301    /// value is scaled by \ref primalScale "primal scale".
1302    ///
1303    /// \pre Either run() or start() must be called before using this function.
1304    Value matchingWeight() const {
1305      Value sum = 0;
1306      for (NodeIt n(_graph); n != INVALID; ++n) {
1307        if ((*_matching)[n] != INVALID) {
1308          sum += _weight[(*_matching)[n]];
1309        }
1310      }
1311      return sum * primalScale / 2;
1312    }
1313
1314    /// \brief Return the number of covered nodes in the matching.
1315    ///
1316    /// This function returns the number of covered nodes in the matching.
1317    ///
1318    /// \pre Either run() or start() must be called before using this function.
1319    int matchingSize() const {
1320      int num = 0;
1321      for (NodeIt n(_graph); n != INVALID; ++n) {
1322        if ((*_matching)[n] != INVALID) {
1323          ++num;
1324        }
1325      }
1326      return num;
1327    }
1328
1329    /// \brief Return \c true if the given edge is in the matching.
1330    ///
1331    /// This function returns \c true if the given edge is in the
1332    /// found matching. The result is scaled by \ref primalScale
1333    /// "primal scale".
1334    ///
1335    /// \pre Either run() or start() must be called before using this function.
1336    int matching(const Edge& edge) const {
1337      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1338        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
1339    }
1340
1341    /// \brief Return the fractional matching arc (or edge) incident
1342    /// to the given node.
1343    ///
1344    /// This function returns one of the fractional matching arc (or
1345    /// edge) incident to the given node in the found matching or \c
1346    /// INVALID if the node is not covered by the matching or if the
1347    /// node is on an odd length cycle then it is the successor edge
1348    /// on the cycle.
1349    ///
1350    /// \pre Either run() or start() must be called before using this function.
1351    Arc matching(const Node& node) const {
1352      return (*_matching)[node];
1353    }
1354
1355    /// \brief Return a const reference to the matching map.
1356    ///
1357    /// This function returns a const reference to a node map that stores
1358    /// the matching arc (or edge) incident to each node.
1359    const MatchingMap& matchingMap() const {
1360      return *_matching;
1361    }
1362
1363    /// @}
1364
1365    /// \name Dual Solution
1366    /// Functions to get the dual solution.\n
1367    /// Either \ref run() or \ref start() function should be called before
1368    /// using them.
1369
1370    /// @{
1371
1372    /// \brief Return the value of the dual solution.
1373    ///
1374    /// This function returns the value of the dual solution.
1375    /// It should be equal to the primal value scaled by \ref dualScale
1376    /// "dual scale".
1377    ///
1378    /// \pre Either run() or start() must be called before using this function.
1379    Value dualValue() const {
1380      Value sum = 0;
1381      for (NodeIt n(_graph); n != INVALID; ++n) {
1382        sum += nodeValue(n);
1383      }
1384      return sum;
1385    }
1386
1387    /// \brief Return the dual value (potential) of the given node.
1388    ///
1389    /// This function returns the dual value (potential) of the given node.
1390    ///
1391    /// \pre Either run() or start() must be called before using this function.
1392    Value nodeValue(const Node& n) const {
1393      return (*_node_potential)[n];
1394    }
1395
1396    /// @}
1397
1398  };
1399
1400  /// \ingroup matching
1401  ///
1402  /// \brief Weighted fractional perfect matching in general graphs
1403  ///
1404  /// This class provides an efficient implementation of fractional
1405  /// matching algorithm. The implementation uses priority queues and
1406  /// provides \f$O(nm\log n)\f$ time complexity.
1407  ///
1408  /// The maximum weighted fractional perfect matching is a relaxation
1409  /// of the maximum weighted perfect matching problem where the odd
1410  /// set constraints are omitted.
1411  /// It can be formulated with the following linear program.
1412  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1413  /// \f[x_e \ge 0\quad \forall e\in E\f]
1414  /// \f[\max \sum_{e\in E}x_ew_e\f]
1415  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1416  /// \f$X\f$. The result must be the union of a matching with one
1417  /// value edges and a set of odd length cycles with half value edges.
1418  ///
1419  /// The algorithm calculates an optimal fractional matching and a
1420  /// proof of the optimality. The solution of the dual problem can be
1421  /// used to check the result of the algorithm. The dual linear
1422  /// problem is the following.
1423  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1424  /// \f[\min \sum_{u \in V}y_u \f]
1425  ///
1426  /// The algorithm can be executed with the run() function.
1427  /// After it the matching (the primal solution) and the dual solution
1428  /// can be obtained using the query functions.
1429  ///
1430  /// The primal solution is multiplied by
1431  /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2".
1432  /// If the value type is integer, then the dual
1433  /// solution is scaled by
1434  /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4".
1435  ///
1436  /// \tparam GR The undirected graph type the algorithm runs on.
1437  /// \tparam WM The type edge weight map. The default type is
1438  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1439#ifdef DOXYGEN
1440  template <typename GR, typename WM>
1441#else
1442  template <typename GR,
1443            typename WM = typename GR::template EdgeMap<int> >
1444#endif
1445  class MaxWeightedPerfectFractionalMatching {
1446  public:
1447
1448    /// The graph type of the algorithm
1449    typedef GR Graph;
1450    /// The type of the edge weight map
1451    typedef WM WeightMap;
1452    /// The value type of the edge weights
1453    typedef typename WeightMap::Value Value;
1454
1455    /// The type of the matching map
1456    typedef typename Graph::template NodeMap<typename Graph::Arc>
1457    MatchingMap;
1458
1459    /// \brief Scaling factor for primal solution
1460    ///
1461    /// Scaling factor for primal solution.
1462    static const int primalScale = 2;
1463
1464    /// \brief Scaling factor for dual solution
1465    ///
1466    /// Scaling factor for dual solution. It is equal to 4 or 1
1467    /// according to the value type.
1468    static const int dualScale =
1469      std::numeric_limits<Value>::is_integer ? 4 : 1;
1470
1471  private:
1472
1473    TEMPLATE_GRAPH_TYPEDEFS(Graph);
1474
1475    typedef typename Graph::template NodeMap<Value> NodePotential;
1476
1477    const Graph& _graph;
1478    const WeightMap& _weight;
1479
1480    MatchingMap* _matching;
1481    NodePotential* _node_potential;
1482
1483    int _node_num;
1484    bool _allow_loops;
1485
1486    enum Status {
1487      EVEN = -1, MATCHED = 0, ODD = 1
1488    };
1489
1490    typedef typename Graph::template NodeMap<Status> StatusMap;
1491    StatusMap* _status;
1492
1493    typedef typename Graph::template NodeMap<Arc> PredMap;
1494    PredMap* _pred;
1495
1496    typedef ExtendFindEnum<IntNodeMap> TreeSet;
1497
1498    IntNodeMap *_tree_set_index;
1499    TreeSet *_tree_set;
1500
1501    IntNodeMap *_delta2_index;
1502    BinHeap<Value, IntNodeMap> *_delta2;
1503
1504    IntEdgeMap *_delta3_index;
1505    BinHeap<Value, IntEdgeMap> *_delta3;
1506
1507    Value _delta_sum;
1508
1509    void createStructures() {
1510      _node_num = countNodes(_graph);
1511
1512      if (!_matching) {
1513        _matching = new MatchingMap(_graph);
1514      }
1515      if (!_node_potential) {
1516        _node_potential = new NodePotential(_graph);
1517      }
1518      if (!_status) {
1519        _status = new StatusMap(_graph);
1520      }
1521      if (!_pred) {
1522        _pred = new PredMap(_graph);
1523      }
1524      if (!_tree_set) {
1525        _tree_set_index = new IntNodeMap(_graph);
1526        _tree_set = new TreeSet(*_tree_set_index);
1527      }
1528      if (!_delta2) {
1529        _delta2_index = new IntNodeMap(_graph);
1530        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1531      }
1532      if (!_delta3) {
1533        _delta3_index = new IntEdgeMap(_graph);
1534        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1535      }
1536    }
1537
1538    void destroyStructures() {
1539      if (_matching) {
1540        delete _matching;
1541      }
1542      if (_node_potential) {
1543        delete _node_potential;
1544      }
1545      if (_status) {
1546        delete _status;
1547      }
1548      if (_pred) {
1549        delete _pred;
1550      }
1551      if (_tree_set) {
1552        delete _tree_set_index;
1553        delete _tree_set;
1554      }
1555      if (_delta2) {
1556        delete _delta2_index;
1557        delete _delta2;
1558      }
1559      if (_delta3) {
1560        delete _delta3_index;
1561        delete _delta3;
1562      }
1563    }
1564
1565    void matchedToEven(Node node, int tree) {
1566      _tree_set->insert(node, tree);
1567      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1568
1569      if (_delta2->state(node) == _delta2->IN_HEAP) {
1570        _delta2->erase(node);
1571      }
1572
1573      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1574        Node v = _graph.source(a);
1575        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1576          dualScale * _weight[a];
1577        if (node == v) {
1578          if (_allow_loops && _graph.direction(a)) {
1579            _delta3->push(a, rw / 2);
1580          }
1581        } else if ((*_status)[v] == EVEN) {
1582          _delta3->push(a, rw / 2);
1583        } else if ((*_status)[v] == MATCHED) {
1584          if (_delta2->state(v) != _delta2->IN_HEAP) {
1585            _pred->set(v, a);
1586            _delta2->push(v, rw);
1587          } else if ((*_delta2)[v] > rw) {
1588            _pred->set(v, a);
1589            _delta2->decrease(v, rw);
1590          }
1591        }
1592      }
1593    }
1594
1595    void matchedToOdd(Node node, int tree) {
1596      _tree_set->insert(node, tree);
1597      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1598
1599      if (_delta2->state(node) == _delta2->IN_HEAP) {
1600        _delta2->erase(node);
1601      }
1602    }
1603
1604    void evenToMatched(Node node, int tree) {
1605      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1606      Arc min = INVALID;
1607      Value minrw = std::numeric_limits<Value>::max();
1608      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1609        Node v = _graph.source(a);
1610        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1611          dualScale * _weight[a];
1612
1613        if (node == v) {
1614          if (_allow_loops && _graph.direction(a)) {
1615            _delta3->erase(a);
1616          }
1617        } else if ((*_status)[v] == EVEN) {
1618          _delta3->erase(a);
1619          if (minrw > rw) {
1620            min = _graph.oppositeArc(a);
1621            minrw = rw;
1622          }
1623        } else if ((*_status)[v]  == MATCHED) {
1624          if ((*_pred)[v] == a) {
1625            Arc mina = INVALID;
1626            Value minrwa = std::numeric_limits<Value>::max();
1627            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1628              Node va = _graph.target(aa);
1629              if ((*_status)[va] != EVEN ||
1630                  _tree_set->find(va) == tree) continue;
1631              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1632                dualScale * _weight[aa];
1633              if (minrwa > rwa) {
1634                minrwa = rwa;
1635                mina = aa;
1636              }
1637            }
1638            if (mina != INVALID) {
1639              _pred->set(v, mina);
1640              _delta2->increase(v, minrwa);
1641            } else {
1642              _pred->set(v, INVALID);
1643              _delta2->erase(v);
1644            }
1645          }
1646        }
1647      }
1648      if (min != INVALID) {
1649        _pred->set(node, min);
1650        _delta2->push(node, minrw);
1651      } else {
1652        _pred->set(node, INVALID);
1653      }
1654    }
1655
1656    void oddToMatched(Node node) {
1657      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1658      Arc min = INVALID;
1659      Value minrw = std::numeric_limits<Value>::max();
1660      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1661        Node v = _graph.source(a);
1662        if ((*_status)[v] != EVEN) continue;
1663        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1664          dualScale * _weight[a];
1665
1666        if (minrw > rw) {
1667          min = _graph.oppositeArc(a);
1668          minrw = rw;
1669        }
1670      }
1671      if (min != INVALID) {
1672        _pred->set(node, min);
1673        _delta2->push(node, minrw);
1674      } else {
1675        _pred->set(node, INVALID);
1676      }
1677    }
1678
1679    void alternatePath(Node even, int tree) {
1680      Node odd;
1681
1682      _status->set(even, MATCHED);
1683      evenToMatched(even, tree);
1684
1685      Arc prev = (*_matching)[even];
1686      while (prev != INVALID) {
1687        odd = _graph.target(prev);
1688        even = _graph.target((*_pred)[odd]);
1689        _matching->set(odd, (*_pred)[odd]);
1690        _status->set(odd, MATCHED);
1691        oddToMatched(odd);
1692
1693        prev = (*_matching)[even];
1694        _status->set(even, MATCHED);
1695        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1696        evenToMatched(even, tree);
1697      }
1698    }
1699
1700    void destroyTree(int tree) {
1701      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1702        if ((*_status)[n] == EVEN) {
1703          _status->set(n, MATCHED);
1704          evenToMatched(n, tree);
1705        } else if ((*_status)[n] == ODD) {
1706          _status->set(n, MATCHED);
1707          oddToMatched(n);
1708        }
1709      }
1710      _tree_set->eraseClass(tree);
1711    }
1712
1713    void augmentOnEdge(const Edge& edge) {
1714      Node left = _graph.u(edge);
1715      int left_tree = _tree_set->find(left);
1716
1717      alternatePath(left, left_tree);
1718      destroyTree(left_tree);
1719      _matching->set(left, _graph.direct(edge, true));
1720
1721      Node right = _graph.v(edge);
1722      int right_tree = _tree_set->find(right);
1723
1724      alternatePath(right, right_tree);
1725      destroyTree(right_tree);
1726      _matching->set(right, _graph.direct(edge, false));
1727    }
1728
1729    void augmentOnArc(const Arc& arc) {
1730      Node left = _graph.source(arc);
1731      _status->set(left, MATCHED);
1732      _matching->set(left, arc);
1733      _pred->set(left, arc);
1734
1735      Node right = _graph.target(arc);
1736      int right_tree = _tree_set->find(right);
1737
1738      alternatePath(right, right_tree);
1739      destroyTree(right_tree);
1740      _matching->set(right, _graph.oppositeArc(arc));
1741    }
1742
1743    void extendOnArc(const Arc& arc) {
1744      Node base = _graph.target(arc);
1745      int tree = _tree_set->find(base);
1746
1747      Node odd = _graph.source(arc);
1748      _tree_set->insert(odd, tree);
1749      _status->set(odd, ODD);
1750      matchedToOdd(odd, tree);
1751      _pred->set(odd, arc);
1752
1753      Node even = _graph.target((*_matching)[odd]);
1754      _tree_set->insert(even, tree);
1755      _status->set(even, EVEN);
1756      matchedToEven(even, tree);
1757    }
1758
1759    void cycleOnEdge(const Edge& edge, int tree) {
1760      Node nca = INVALID;
1761      std::vector<Node> left_path, right_path;
1762
1763      {
1764        std::set<Node> left_set, right_set;
1765        Node left = _graph.u(edge);
1766        left_path.push_back(left);
1767        left_set.insert(left);
1768
1769        Node right = _graph.v(edge);
1770        right_path.push_back(right);
1771        right_set.insert(right);
1772
1773        while (true) {
1774
1775          if (left_set.find(right) != left_set.end()) {
1776            nca = right;
1777            break;
1778          }
1779
1780          if ((*_matching)[left] == INVALID) break;
1781
1782          left = _graph.target((*_matching)[left]);
1783          left_path.push_back(left);
1784          left = _graph.target((*_pred)[left]);
1785          left_path.push_back(left);
1786
1787          left_set.insert(left);
1788
1789          if (right_set.find(left) != right_set.end()) {
1790            nca = left;
1791            break;
1792          }
1793
1794          if ((*_matching)[right] == INVALID) break;
1795
1796          right = _graph.target((*_matching)[right]);
1797          right_path.push_back(right);
1798          right = _graph.target((*_pred)[right]);
1799          right_path.push_back(right);
1800
1801          right_set.insert(right);
1802
1803        }
1804
1805        if (nca == INVALID) {
1806          if ((*_matching)[left] == INVALID) {
1807            nca = right;
1808            while (left_set.find(nca) == left_set.end()) {
1809              nca = _graph.target((*_matching)[nca]);
1810              right_path.push_back(nca);
1811              nca = _graph.target((*_pred)[nca]);
1812              right_path.push_back(nca);
1813            }
1814          } else {
1815            nca = left;
1816            while (right_set.find(nca) == right_set.end()) {
1817              nca = _graph.target((*_matching)[nca]);
1818              left_path.push_back(nca);
1819              nca = _graph.target((*_pred)[nca]);
1820              left_path.push_back(nca);
1821            }
1822          }
1823        }
1824      }
1825
1826      alternatePath(nca, tree);
1827      Arc prev;
1828
1829      prev = _graph.direct(edge, true);
1830      for (int i = 0; left_path[i] != nca; i += 2) {
1831        _matching->set(left_path[i], prev);
1832        _status->set(left_path[i], MATCHED);
1833        evenToMatched(left_path[i], tree);
1834
1835        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1836        _status->set(left_path[i + 1], MATCHED);
1837        oddToMatched(left_path[i + 1]);
1838      }
1839      _matching->set(nca, prev);
1840
1841      for (int i = 0; right_path[i] != nca; i += 2) {
1842        _status->set(right_path[i], MATCHED);
1843        evenToMatched(right_path[i], tree);
1844
1845        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1846        _status->set(right_path[i + 1], MATCHED);
1847        oddToMatched(right_path[i + 1]);
1848      }
1849
1850      destroyTree(tree);
1851    }
1852
1853    void extractCycle(const Arc &arc) {
1854      Node left = _graph.source(arc);
1855      Node odd = _graph.target((*_matching)[left]);
1856      Arc prev;
1857      while (odd != left) {
1858        Node even = _graph.target((*_matching)[odd]);
1859        prev = (*_matching)[odd];
1860        odd = _graph.target((*_matching)[even]);
1861        _matching->set(even, _graph.oppositeArc(prev));
1862      }
1863      _matching->set(left, arc);
1864
1865      Node right = _graph.target(arc);
1866      int right_tree = _tree_set->find(right);
1867      alternatePath(right, right_tree);
1868      destroyTree(right_tree);
1869      _matching->set(right, _graph.oppositeArc(arc));
1870    }
1871
1872  public:
1873
1874    /// \brief Constructor
1875    ///
1876    /// Constructor.
1877    MaxWeightedPerfectFractionalMatching(const Graph& graph,
1878                                         const WeightMap& weight,
1879                                         bool allow_loops = true)
1880      : _graph(graph), _weight(weight), _matching(0),
1881      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1882      _status(0),  _pred(0),
1883      _tree_set_index(0), _tree_set(0),
1884
1885      _delta2_index(0), _delta2(0),
1886      _delta3_index(0), _delta3(0),
1887
1888      _delta_sum() {}
1889
1890    ~MaxWeightedPerfectFractionalMatching() {
1891      destroyStructures();
1892    }
1893
1894    /// \name Execution Control
1895    /// The simplest way to execute the algorithm is to use the
1896    /// \ref run() member function.
1897
1898    ///@{
1899
1900    /// \brief Initialize the algorithm
1901    ///
1902    /// This function initializes the algorithm.
1903    void init() {
1904      createStructures();
1905
1906      for (NodeIt n(_graph); n != INVALID; ++n) {
1907        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1908      }
1909      for (EdgeIt e(_graph); e != INVALID; ++e) {
1910        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1911      }
1912
1913      _delta2->clear();
1914      _delta3->clear();
1915      _tree_set->clear();
1916
1917      for (NodeIt n(_graph); n != INVALID; ++n) {
1918        Value max = - std::numeric_limits<Value>::max();
1919        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1920          if (_graph.target(e) == n && !_allow_loops) continue;
1921          if ((dualScale * _weight[e]) / 2 > max) {
1922            max = (dualScale * _weight[e]) / 2;
1923          }
1924        }
1925        _node_potential->set(n, max);
1926
1927        _tree_set->insert(n);
1928
1929        _matching->set(n, INVALID);
1930        _status->set(n, EVEN);
1931      }
1932
1933      for (EdgeIt e(_graph); e != INVALID; ++e) {
1934        Node left = _graph.u(e);
1935        Node right = _graph.v(e);
1936        if (left == right && !_allow_loops) continue;
1937        _delta3->push(e, ((*_node_potential)[left] +
1938                          (*_node_potential)[right] -
1939                          dualScale * _weight[e]) / 2);
1940      }
1941    }
1942
1943    /// \brief Start the algorithm
1944    ///
1945    /// This function starts the algorithm.
1946    ///
1947    /// \pre \ref init() must be called before using this function.
1948    bool start() {
1949      enum OpType {
1950        D2, D3
1951      };
1952
1953      int unmatched = _node_num;
1954      while (unmatched > 0) {
1955        Value d2 = !_delta2->empty() ?
1956          _delta2->prio() : std::numeric_limits<Value>::max();
1957
1958        Value d3 = !_delta3->empty() ?
1959          _delta3->prio() : std::numeric_limits<Value>::max();
1960
1961        _delta_sum = d3; OpType ot = D3;
1962        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1963
1964        if (_delta_sum == std::numeric_limits<Value>::max()) {
1965          return false;
1966        }
1967
1968        switch (ot) {
1969        case D2:
1970          {
1971            Node n = _delta2->top();
1972            Arc a = (*_pred)[n];
1973            if ((*_matching)[n] == INVALID) {
1974              augmentOnArc(a);
1975              --unmatched;
1976            } else {
1977              Node v = _graph.target((*_matching)[n]);
1978              if ((*_matching)[n] !=
1979                  _graph.oppositeArc((*_matching)[v])) {
1980                extractCycle(a);
1981                --unmatched;
1982              } else {
1983                extendOnArc(a);
1984              }
1985            }
1986          } break;
1987        case D3:
1988          {
1989            Edge e = _delta3->top();
1990
1991            Node left = _graph.u(e);
1992            Node right = _graph.v(e);
1993
1994            int left_tree = _tree_set->find(left);
1995            int right_tree = _tree_set->find(right);
1996
1997            if (left_tree == right_tree) {
1998              cycleOnEdge(e, left_tree);
1999              --unmatched;
2000            } else {
2001              augmentOnEdge(e);
2002              unmatched -= 2;
2003            }
2004          } break;
2005        }
2006      }
2007      return true;
2008    }
2009
2010    /// \brief Run the algorithm.
2011    ///
2012    /// This method runs the \c %MaxWeightedPerfectFractionalMatching
2013    /// algorithm.
2014    ///
2015    /// \note mwfm.run() is just a shortcut of the following code.
2016    /// \code
2017    ///   mwpfm.init();
2018    ///   mwpfm.start();
2019    /// \endcode
2020    bool run() {
2021      init();
2022      return start();
2023    }
2024
2025    /// @}
2026
2027    /// \name Primal Solution
2028    /// Functions to get the primal solution, i.e. the maximum weighted
2029    /// matching.\n
2030    /// Either \ref run() or \ref start() function should be called before
2031    /// using them.
2032
2033    /// @{
2034
2035    /// \brief Return the weight of the matching.
2036    ///
2037    /// This function returns the weight of the found matching. This
2038    /// value is scaled by \ref primalScale "primal scale".
2039    ///
2040    /// \pre Either run() or start() must be called before using this function.
2041    Value matchingWeight() const {
2042      Value sum = 0;
2043      for (NodeIt n(_graph); n != INVALID; ++n) {
2044        if ((*_matching)[n] != INVALID) {
2045          sum += _weight[(*_matching)[n]];
2046        }
2047      }
2048      return sum * primalScale / 2;
2049    }
2050
2051    /// \brief Return the number of covered nodes in the matching.
2052    ///
2053    /// This function returns the number of covered nodes in the matching.
2054    ///
2055    /// \pre Either run() or start() must be called before using this function.
2056    int matchingSize() const {
2057      int num = 0;
2058      for (NodeIt n(_graph); n != INVALID; ++n) {
2059        if ((*_matching)[n] != INVALID) {
2060          ++num;
2061        }
2062      }
2063      return num;
2064    }
2065
2066    /// \brief Return \c true if the given edge is in the matching.
2067    ///
2068    /// This function returns \c true if the given edge is in the
2069    /// found matching. The result is scaled by \ref primalScale
2070    /// "primal scale".
2071    ///
2072    /// \pre Either run() or start() must be called before using this function.
2073    int matching(const Edge& edge) const {
2074      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
2075        + (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
2076    }
2077
2078    /// \brief Return the fractional matching arc (or edge) incident
2079    /// to the given node.
2080    ///
2081    /// This function returns one of the fractional matching arc (or
2082    /// edge) incident to the given node in the found matching or \c
2083    /// INVALID if the node is not covered by the matching or if the
2084    /// node is on an odd length cycle then it is the successor edge
2085    /// on the cycle.
2086    ///
2087    /// \pre Either run() or start() must be called before using this function.
2088    Arc matching(const Node& node) const {
2089      return (*_matching)[node];
2090    }
2091
2092    /// \brief Return a const reference to the matching map.
2093    ///
2094    /// This function returns a const reference to a node map that stores
2095    /// the matching arc (or edge) incident to each node.
2096    const MatchingMap& matchingMap() const {
2097      return *_matching;
2098    }
2099
2100    /// @}
2101
2102    /// \name Dual Solution
2103    /// Functions to get the dual solution.\n
2104    /// Either \ref run() or \ref start() function should be called before
2105    /// using them.
2106
2107    /// @{
2108
2109    /// \brief Return the value of the dual solution.
2110    ///
2111    /// This function returns the value of the dual solution.
2112    /// It should be equal to the primal value scaled by \ref dualScale
2113    /// "dual scale".
2114    ///
2115    /// \pre Either run() or start() must be called before using this function.
2116    Value dualValue() const {
2117      Value sum = 0;
2118      for (NodeIt n(_graph); n != INVALID; ++n) {
2119        sum += nodeValue(n);
2120      }
2121      return sum;
2122    }
2123
2124    /// \brief Return the dual value (potential) of the given node.
2125    ///
2126    /// This function returns the dual value (potential) of the given node.
2127    ///
2128    /// \pre Either run() or start() must be called before using this function.
2129    Value nodeValue(const Node& n) const {
2130      return (*_node_potential)[n];
2131    }
2132
2133    /// @}
2134
2135  };
2136
2137} //END OF NAMESPACE LEMON
2138
2139#endif //LEMON_FRACTIONAL_MATCHING_H
Note: See TracBrowser for help on using the repository browser.