1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2009
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
19 #ifndef LEMON_FRACTIONAL_MATCHING_H
20 #define LEMON_FRACTIONAL_MATCHING_H
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>
36 ///\brief Fractional matching algorithms in general graphs.
40 /// \brief Default traits class of MaxFractionalMatching class.
42 /// Default traits class of MaxFractionalMatching class.
43 /// \tparam GR Graph type.
44 template <typename GR>
45 struct MaxFractionalMatchingDefaultTraits {
47 /// \brief The type of the graph the algorithm runs on.
50 /// \brief The type of the map that stores the matching.
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;
56 /// \brief Instantiates a MatchingMap.
58 /// This function instantiates a \ref MatchingMap.
59 /// \param graph The graph for which we would like to define
61 static MatchingMap* createMatchingMap(const Graph& graph) {
62 return new MatchingMap(graph);
65 /// \brief The elevator type used by MaxFractionalMatching algorithm.
67 /// The elevator type used by MaxFractionalMatching algorithm.
70 /// \sa LinkedElevator
71 typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
73 /// \brief Instantiates an Elevator.
75 /// This function instantiates an \ref Elevator.
76 /// \param graph The graph for which we would like to define
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);
86 /// \brief Max cardinality fractional matching
88 /// This class provides an implementation of fractional matching
89 /// algorithm based on push-relabel principle.
91 /// The maximum cardinality fractional matching is a relaxation of the
92 /// maximum cardinality matching problem where the odd set constraints
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.
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.
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.
113 /// The primal solution is multiplied by
114 /// \ref MaxFractionalMatching::primalScale "2".
116 /// \tparam GR The undirected graph type the algorithm runs on.
118 template <typename GR, typename TR>
120 template <typename GR,
121 typename TR = MaxFractionalMatchingDefaultTraits<GR> >
123 class MaxFractionalMatching {
126 /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
127 /// class" of the algorithm.
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;
136 /// \brief Scaling factor for primal solution
138 /// Scaling factor for primal solution.
139 static const int primalScale = 2;
148 TEMPLATE_GRAPH_TYPEDEFS(Graph);
150 bool _local_matching;
151 MatchingMap *_matching;
156 typedef typename Graph::template NodeMap<int> InDegMap;
159 void createStructures() {
160 _node_num = countNodes(_graph);
163 _local_matching = true;
164 _matching = Traits::createMatchingMap(_graph);
168 _level = Traits::createElevator(_graph, _node_num);
171 _indeg = new InDegMap(_graph);
175 void destroyStructures() {
176 if (_local_matching) {
187 void postprocessing() {
188 for (NodeIt n(_graph); n != INVALID; ++n) {
189 if ((*_indeg)[n] != 0) continue;
192 while ((*_matching)[u] != INVALID) {
193 Node v = _graph.target((*_matching)[u]);
195 Arc a = _graph.oppositeArc((*_matching)[u]);
196 u = _graph.target((*_matching)[v]);
198 _matching->set(v, a);
202 for (NodeIt n(_graph); n != INVALID; ++n) {
203 if ((*_indeg)[n] != 1) continue;
207 Node u = _graph.target((*_matching)[n]);
210 u = _graph.target((*_matching)[u]);
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);
219 prev = _graph.oppositeArc((*_matching)[u]);
220 v = _graph.target((*_matching)[u]);
221 u = _graph.target((*_matching)[v]);
222 _matching->set(v, prev);
230 typedef MaxFractionalMatching Create;
232 ///\name Named Template Parameters
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
245 /// \brief \ref named-templ-param "Named parameter" for setting
248 /// \ref named-templ-param "Named parameter" for setting MatchingMap
250 template <typename T>
251 struct SetMatchingMap
252 : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
253 typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
256 template <typename T>
257 struct SetElevatorTraits : public Traits {
259 static Elevator *createElevator(const Graph&, int) {
260 LEMON_ASSERT(false, "Elevator is not initialized");
261 return 0; // ignore warnings
265 /// \brief \ref named-templ-param "Named parameter" for setting
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>
276 : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
277 typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
280 template <typename T>
281 struct SetStandardElevatorTraits : public Traits {
283 static Elevator *createElevator(const Graph& graph, int max_level) {
284 return new Elevator(graph, max_level);
288 /// \brief \ref named-templ-param "Named parameter" for setting
289 /// Elevator type with automatic allocation
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().
300 template <typename T>
301 struct SetStandardElevator
302 : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
303 typedef MaxFractionalMatching<Graph,
304 SetStandardElevatorTraits<T> > Create;
311 MaxFractionalMatching() {}
315 /// \brief Constructor
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)
325 ~MaxFractionalMatching() {
329 /// \brief Sets the matching map.
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,
336 /// \return <tt>(*this)</tt>
337 MaxFractionalMatching& matchingMap(MatchingMap& map) {
338 if (_local_matching) {
340 _local_matching = false;
346 /// \brief Sets the elevator used by algorithm.
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,
353 /// \return <tt>(*this)</tt>
354 MaxFractionalMatching& elevator(Elevator& elevator) {
357 _local_level = false;
363 /// \brief Returns a const reference to the elevator.
365 /// Returns a const reference to the elevator.
367 /// \pre Either \ref run() or \ref init() must be called before
368 /// using this function.
369 const Elevator& elevator() const {
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()
382 /// \brief Initializes the internal data structures.
384 /// Initializes the internal data structures and sets the initial
390 for (NodeIt n(_graph); n != INVALID; ++n) {
392 _matching->set(n, INVALID);
393 _level->initAddItem(n);
395 _level->initFinish();
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);
408 for (NodeIt n(_graph); n != INVALID; ++n) {
409 if ((*_indeg)[n] == 0) {
415 /// \brief Starts the algorithm and computes a fractional matching
417 /// The algorithm computes a maximum fractional matching.
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) {
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) {
439 _matching->set(u, a);
440 _indeg->set(n, (*_indeg)[n] + 1);
441 _level->deactivate(n);
443 } else if (new_level > (*_level)[v]) {
444 new_level = (*_level)[v];
448 if (new_level + 1 < _level->maxLevel()) {
449 _level->liftHighestActive(new_level + 1);
451 _level->liftHighestActiveToTop();
453 if (_level->emptyLevel(level)) {
454 _level->liftToTop(level);
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);
472 /// \brief Starts the algorithm and computes a perfect fractional
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.
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) {
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) {
499 _matching->set(u, a);
500 _indeg->set(n, (*_indeg)[n] + 1);
501 _level->deactivate(n);
503 } else if (new_level > (*_level)[v]) {
504 new_level = (*_level)[v];
508 if (new_level + 1 < _level->maxLevel()) {
509 _level->liftHighestActive(new_level + 1);
511 _level->liftHighestActiveToTop();
512 _empty_level = _level->maxLevel() - 1;
515 if (_level->emptyLevel(level)) {
516 _level->liftToTop(level);
517 _empty_level = level;
529 /// \brief Runs the algorithm
531 /// Just a shortcut for the next code:
536 void run(bool postprocess = true) {
541 /// \brief Runs the algorithm to find a perfect fractional matching
543 /// Just a shortcut for the next code:
548 bool runPerfect(bool postprocess = true) {
550 return startPerfect(postprocess);
555 /// \name Query Functions
556 /// The result of the %Matching algorithm can be obtained using these
558 /// Before the use of these functions,
559 /// either run() or start() must be called.
563 /// \brief Return the number of covered nodes in the matching.
565 /// This function returns the number of covered nodes in the matching.
567 /// \pre Either run() or start() must be called before using this function.
568 int matchingSize() const {
570 for (NodeIt n(_graph); n != INVALID; ++n) {
571 if ((*_matching)[n] != INVALID) {
578 /// \brief Returns a const reference to the matching map.
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.
584 /// \pre Either \ref run() or \ref init() must be called before
585 /// using this function.
586 const MatchingMap& matchingMap() const {
590 /// \brief Return \c true if the given edge is in the matching.
592 /// This function returns \c true if the given edge is in the
593 /// found matching. The result is scaled by \ref primalScale
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);
602 /// \brief Return the fractional matching arc (or edge) incident
603 /// to the given node.
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
611 /// \pre Either run() or start() must be called before using this function.
612 Arc matching(const Node& node) const {
613 return (*_matching)[node];
616 /// \brief Returns true if the node is in the barrier
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;
630 /// \ingroup matching
632 /// \brief Weighted fractional matching in general graphs
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.
638 /// The maximum weighted fractional matching is a relaxation of the
639 /// maximum weighted matching problem where the odd set constraints
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.
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]
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.
661 /// If the value type is integer, then the primal and the dual
662 /// solutions are multiplied by
663 /// \ref MaxWeightedFractionalMatching::primalScale "2" and
664 /// \ref MaxWeightedFractionalMatching::dualScale "4" respectively.
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>".
670 template <typename GR, typename WM>
672 template <typename GR,
673 typename WM = typename GR::template EdgeMap<int> >
675 class MaxWeightedFractionalMatching {
678 /// The graph type of the algorithm
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;
685 /// The type of the matching map
686 typedef typename Graph::template NodeMap<typename Graph::Arc>
689 /// \brief Scaling factor for primal solution
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;
696 /// \brief Scaling factor for dual solution
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;
705 TEMPLATE_GRAPH_TYPEDEFS(Graph);
707 typedef typename Graph::template NodeMap<Value> NodePotential;
710 const WeightMap& _weight;
712 MatchingMap* _matching;
713 NodePotential* _node_potential;
719 EVEN = -1, MATCHED = 0, ODD = 1
722 typedef typename Graph::template NodeMap<Status> StatusMap;
725 typedef typename Graph::template NodeMap<Arc> PredMap;
728 typedef ExtendFindEnum<IntNodeMap> TreeSet;
730 IntNodeMap *_tree_set_index;
733 IntNodeMap *_delta1_index;
734 BinHeap<Value, IntNodeMap> *_delta1;
736 IntNodeMap *_delta2_index;
737 BinHeap<Value, IntNodeMap> *_delta2;
739 IntEdgeMap *_delta3_index;
740 BinHeap<Value, IntEdgeMap> *_delta3;
744 void createStructures() {
745 _node_num = countNodes(_graph);
748 _matching = new MatchingMap(_graph);
750 if (!_node_potential) {
751 _node_potential = new NodePotential(_graph);
754 _status = new StatusMap(_graph);
757 _pred = new PredMap(_graph);
760 _tree_set_index = new IntNodeMap(_graph);
761 _tree_set = new TreeSet(*_tree_set_index);
764 _delta1_index = new IntNodeMap(_graph);
765 _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
768 _delta2_index = new IntNodeMap(_graph);
769 _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
772 _delta3_index = new IntEdgeMap(_graph);
773 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
777 void destroyStructures() {
781 if (_node_potential) {
782 delete _node_potential;
791 delete _tree_set_index;
795 delete _delta1_index;
799 delete _delta2_index;
803 delete _delta3_index;
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]);
813 if (_delta2->state(node) == _delta2->IN_HEAP) {
814 _delta2->erase(node);
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];
822 if (_allow_loops && _graph.direction(a)) {
823 _delta3->push(a, rw / 2);
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) {
830 _delta2->push(v, rw);
831 } else if ((*_delta2)[v] > rw) {
833 _delta2->decrease(v, rw);
839 void matchedToOdd(Node node, int tree) {
840 _tree_set->insert(node, tree);
841 _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
843 if (_delta2->state(node) == _delta2->IN_HEAP) {
844 _delta2->erase(node);
848 void evenToMatched(Node node, int tree) {
849 _delta1->erase(node);
850 _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
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];
859 if (_allow_loops && _graph.direction(a)) {
862 } else if ((*_status)[v] == EVEN) {
865 min = _graph.oppositeArc(a);
868 } else if ((*_status)[v] == MATCHED) {
869 if ((*_pred)[v] == a) {
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];
883 if (mina != INVALID) {
885 _delta2->increase(v, minrwa);
887 _pred->set(v, INVALID);
893 if (min != INVALID) {
894 _pred->set(node, min);
895 _delta2->push(node, minrw);
897 _pred->set(node, INVALID);
901 void oddToMatched(Node node) {
902 _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
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];
912 min = _graph.oppositeArc(a);
916 if (min != INVALID) {
917 _pred->set(node, min);
918 _delta2->push(node, minrw);
920 _pred->set(node, INVALID);
924 void alternatePath(Node even, int tree) {
927 _status->set(even, MATCHED);
928 evenToMatched(even, tree);
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);
938 prev = (*_matching)[even];
939 _status->set(even, MATCHED);
940 _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
941 evenToMatched(even, tree);
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);
955 _tree_set->eraseClass(tree);
959 void unmatchNode(const Node& node) {
960 int tree = _tree_set->find(node);
962 alternatePath(node, tree);
965 _matching->set(node, INVALID);
969 void augmentOnEdge(const Edge& edge) {
970 Node left = _graph.u(edge);
971 int left_tree = _tree_set->find(left);
973 alternatePath(left, left_tree);
974 destroyTree(left_tree);
975 _matching->set(left, _graph.direct(edge, true));
977 Node right = _graph.v(edge);
978 int right_tree = _tree_set->find(right);
980 alternatePath(right, right_tree);
981 destroyTree(right_tree);
982 _matching->set(right, _graph.direct(edge, false));
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);
991 Node right = _graph.target(arc);
992 int right_tree = _tree_set->find(right);
994 alternatePath(right, right_tree);
995 destroyTree(right_tree);
996 _matching->set(right, _graph.oppositeArc(arc));
999 void extendOnArc(const Arc& arc) {
1000 Node base = _graph.target(arc);
1001 int tree = _tree_set->find(base);
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);
1009 Node even = _graph.target((*_matching)[odd]);
1010 _tree_set->insert(even, tree);
1011 _status->set(even, EVEN);
1012 matchedToEven(even, tree);
1015 void cycleOnEdge(const Edge& edge, int tree) {
1017 std::vector<Node> left_path, right_path;
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);
1025 Node right = _graph.v(edge);
1026 right_path.push_back(right);
1027 right_set.insert(right);
1031 if (left_set.find(right) != left_set.end()) {
1036 if ((*_matching)[left] == INVALID) break;
1038 left = _graph.target((*_matching)[left]);
1039 left_path.push_back(left);
1040 left = _graph.target((*_pred)[left]);
1041 left_path.push_back(left);
1043 left_set.insert(left);
1045 if (right_set.find(left) != right_set.end()) {
1050 if ((*_matching)[right] == INVALID) break;
1052 right = _graph.target((*_matching)[right]);
1053 right_path.push_back(right);
1054 right = _graph.target((*_pred)[right]);
1055 right_path.push_back(right);
1057 right_set.insert(right);
1061 if (nca == INVALID) {
1062 if ((*_matching)[left] == INVALID) {
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);
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);
1082 alternatePath(nca, tree);
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);
1091 prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1092 _status->set(left_path[i + 1], MATCHED);
1093 oddToMatched(left_path[i + 1]);
1095 _matching->set(nca, prev);
1097 for (int i = 0; right_path[i] != nca; i += 2) {
1098 _status->set(right_path[i], MATCHED);
1099 evenToMatched(right_path[i], tree);
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]);
1109 void extractCycle(const Arc &arc) {
1110 Node left = _graph.source(arc);
1111 Node odd = _graph.target((*_matching)[left]);
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));
1119 _matching->set(left, arc);
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));
1130 /// \brief 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),
1140 _delta1_index(0), _delta1(0),
1141 _delta2_index(0), _delta2(0),
1142 _delta3_index(0), _delta3(0),
1146 ~MaxWeightedFractionalMatching() {
1147 destroyStructures();
1150 /// \name Execution Control
1151 /// The simplest way to execute the algorithm is to use the
1152 /// \ref run() member function.
1156 /// \brief Initialize the algorithm
1158 /// This function initializes the algorithm.
1162 for (NodeIt n(_graph); n != INVALID; ++n) {
1163 (*_delta1_index)[n] = _delta1->PRE_HEAP;
1164 (*_delta2_index)[n] = _delta2->PRE_HEAP;
1166 for (EdgeIt e(_graph); e != INVALID; ++e) {
1167 (*_delta3_index)[e] = _delta3->PRE_HEAP;
1170 for (NodeIt n(_graph); n != INVALID; ++n) {
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;
1178 _node_potential->set(n, max);
1179 _delta1->push(n, max);
1181 _tree_set->insert(n);
1183 _matching->set(n, INVALID);
1184 _status->set(n, EVEN);
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);
1197 /// \brief Start the algorithm
1199 /// This function starts the algorithm.
1201 /// \pre \ref init() must be called before using this function.
1207 int unmatched = _node_num;
1208 while (unmatched > 0) {
1209 Value d1 = !_delta1->empty() ?
1210 _delta1->prio() : std::numeric_limits<Value>::max();
1212 Value d2 = !_delta2->empty() ?
1213 _delta2->prio() : std::numeric_limits<Value>::max();
1215 Value d3 = !_delta3->empty() ?
1216 _delta3->prio() : std::numeric_limits<Value>::max();
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; }
1225 Node n = _delta1->top();
1232 Node n = _delta2->top();
1233 Arc a = (*_pred)[n];
1234 if ((*_matching)[n] == INVALID) {
1238 Node v = _graph.target((*_matching)[n]);
1239 if ((*_matching)[n] !=
1240 _graph.oppositeArc((*_matching)[v])) {
1250 Edge e = _delta3->top();
1252 Node left = _graph.u(e);
1253 Node right = _graph.v(e);
1255 int left_tree = _tree_set->find(left);
1256 int right_tree = _tree_set->find(right);
1258 if (left_tree == right_tree) {
1259 cycleOnEdge(e, left_tree);
1270 /// \brief Run the algorithm.
1272 /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
1274 /// \note mwfm.run() is just a shortcut of the following code.
1286 /// \name Primal Solution
1287 /// Functions to get the primal solution, i.e. the maximum weighted
1289 /// Either \ref run() or \ref start() function should be called before
1294 /// \brief Return the weight of the matching.
1296 /// This function returns the weight of the found matching. This
1297 /// value is scaled by \ref primalScale "primal scale".
1299 /// \pre Either run() or start() must be called before using this function.
1300 Value matchingWeight() const {
1302 for (NodeIt n(_graph); n != INVALID; ++n) {
1303 if ((*_matching)[n] != INVALID) {
1304 sum += _weight[(*_matching)[n]];
1307 return sum * primalScale / 2;
1310 /// \brief Return the number of covered nodes in the matching.
1312 /// This function returns the number of covered nodes in the matching.
1314 /// \pre Either run() or start() must be called before using this function.
1315 int matchingSize() const {
1317 for (NodeIt n(_graph); n != INVALID; ++n) {
1318 if ((*_matching)[n] != INVALID) {
1325 /// \brief Return \c true if the given edge is in the matching.
1327 /// This function returns \c true if the given edge is in the
1328 /// found matching. The result is scaled by \ref primalScale
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)
1338 /// \brief Return the fractional matching arc (or edge) incident
1339 /// to the given node.
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
1347 /// \pre Either run() or start() must be called before using this function.
1348 Arc matching(const Node& node) const {
1349 return (*_matching)[node];
1352 /// \brief Return a const reference to the matching map.
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 {
1362 /// \name Dual Solution
1363 /// Functions to get the dual solution.\n
1364 /// Either \ref run() or \ref start() function should be called before
1369 /// \brief Return the value of the dual solution.
1371 /// This function returns the value of the dual solution.
1372 /// It should be equal to the primal value scaled by \ref dualScale
1375 /// \pre Either run() or start() must be called before using this function.
1376 Value dualValue() const {
1378 for (NodeIt n(_graph); n != INVALID; ++n) {
1379 sum += nodeValue(n);
1384 /// \brief Return the dual value (potential) of the given node.
1386 /// This function returns the dual value (potential) of the given node.
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];
1397 /// \ingroup matching
1399 /// \brief Weighted fractional perfect matching in general graphs
1401 /// This class provides an efficient implementation of fractional
1402 /// matching algorithm. The implementation uses priority queues and
1403 /// provides \f$O(nm\log n)\f$ time complexity.
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.
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]
1421 /// \f[\min \sum_{u \in V}y_u \f]
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.
1427 /// If the value type is integer, then the primal and the dual
1428 /// solutions are multiplied by
1429 /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2" and
1430 /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4" respectively.
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>".
1436 template <typename GR, typename WM>
1438 template <typename GR,
1439 typename WM = typename GR::template EdgeMap<int> >
1441 class MaxWeightedPerfectFractionalMatching {
1444 /// The graph type of the algorithm
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;
1451 /// The type of the matching map
1452 typedef typename Graph::template NodeMap<typename Graph::Arc>
1455 /// \brief Scaling factor for primal solution
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;
1462 /// \brief Scaling factor for dual solution
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;
1471 TEMPLATE_GRAPH_TYPEDEFS(Graph);
1473 typedef typename Graph::template NodeMap<Value> NodePotential;
1475 const Graph& _graph;
1476 const WeightMap& _weight;
1478 MatchingMap* _matching;
1479 NodePotential* _node_potential;
1485 EVEN = -1, MATCHED = 0, ODD = 1
1488 typedef typename Graph::template NodeMap<Status> StatusMap;
1491 typedef typename Graph::template NodeMap<Arc> PredMap;
1494 typedef ExtendFindEnum<IntNodeMap> TreeSet;
1496 IntNodeMap *_tree_set_index;
1499 IntNodeMap *_delta2_index;
1500 BinHeap<Value, IntNodeMap> *_delta2;
1502 IntEdgeMap *_delta3_index;
1503 BinHeap<Value, IntEdgeMap> *_delta3;
1507 void createStructures() {
1508 _node_num = countNodes(_graph);
1511 _matching = new MatchingMap(_graph);
1513 if (!_node_potential) {
1514 _node_potential = new NodePotential(_graph);
1517 _status = new StatusMap(_graph);
1520 _pred = new PredMap(_graph);
1523 _tree_set_index = new IntNodeMap(_graph);
1524 _tree_set = new TreeSet(*_tree_set_index);
1527 _delta2_index = new IntNodeMap(_graph);
1528 _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1531 _delta3_index = new IntEdgeMap(_graph);
1532 _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1536 void destroyStructures() {
1540 if (_node_potential) {
1541 delete _node_potential;
1550 delete _tree_set_index;
1554 delete _delta2_index;
1558 delete _delta3_index;
1563 void matchedToEven(Node node, int tree) {
1564 _tree_set->insert(node, tree);
1565 _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1567 if (_delta2->state(node) == _delta2->IN_HEAP) {
1568 _delta2->erase(node);
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];
1576 if (_allow_loops && _graph.direction(a)) {
1577 _delta3->push(a, rw / 2);
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) {
1584 _delta2->push(v, rw);
1585 } else if ((*_delta2)[v] > rw) {
1587 _delta2->decrease(v, rw);
1593 void matchedToOdd(Node node, int tree) {
1594 _tree_set->insert(node, tree);
1595 _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1597 if (_delta2->state(node) == _delta2->IN_HEAP) {
1598 _delta2->erase(node);
1602 void evenToMatched(Node node, int tree) {
1603 _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
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];
1612 if (_allow_loops && _graph.direction(a)) {
1615 } else if ((*_status)[v] == EVEN) {
1618 min = _graph.oppositeArc(a);
1621 } else if ((*_status)[v] == MATCHED) {
1622 if ((*_pred)[v] == a) {
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];
1636 if (mina != INVALID) {
1637 _pred->set(v, mina);
1638 _delta2->increase(v, minrwa);
1640 _pred->set(v, INVALID);
1646 if (min != INVALID) {
1647 _pred->set(node, min);
1648 _delta2->push(node, minrw);
1650 _pred->set(node, INVALID);
1654 void oddToMatched(Node node) {
1655 _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
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];
1665 min = _graph.oppositeArc(a);
1669 if (min != INVALID) {
1670 _pred->set(node, min);
1671 _delta2->push(node, minrw);
1673 _pred->set(node, INVALID);
1677 void alternatePath(Node even, int tree) {
1680 _status->set(even, MATCHED);
1681 evenToMatched(even, tree);
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);
1691 prev = (*_matching)[even];
1692 _status->set(even, MATCHED);
1693 _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1694 evenToMatched(even, tree);
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);
1708 _tree_set->eraseClass(tree);
1711 void augmentOnEdge(const Edge& edge) {
1712 Node left = _graph.u(edge);
1713 int left_tree = _tree_set->find(left);
1715 alternatePath(left, left_tree);
1716 destroyTree(left_tree);
1717 _matching->set(left, _graph.direct(edge, true));
1719 Node right = _graph.v(edge);
1720 int right_tree = _tree_set->find(right);
1722 alternatePath(right, right_tree);
1723 destroyTree(right_tree);
1724 _matching->set(right, _graph.direct(edge, false));
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);
1733 Node right = _graph.target(arc);
1734 int right_tree = _tree_set->find(right);
1736 alternatePath(right, right_tree);
1737 destroyTree(right_tree);
1738 _matching->set(right, _graph.oppositeArc(arc));
1741 void extendOnArc(const Arc& arc) {
1742 Node base = _graph.target(arc);
1743 int tree = _tree_set->find(base);
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);
1751 Node even = _graph.target((*_matching)[odd]);
1752 _tree_set->insert(even, tree);
1753 _status->set(even, EVEN);
1754 matchedToEven(even, tree);
1757 void cycleOnEdge(const Edge& edge, int tree) {
1759 std::vector<Node> left_path, right_path;
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);
1767 Node right = _graph.v(edge);
1768 right_path.push_back(right);
1769 right_set.insert(right);
1773 if (left_set.find(right) != left_set.end()) {
1778 if ((*_matching)[left] == INVALID) break;
1780 left = _graph.target((*_matching)[left]);
1781 left_path.push_back(left);
1782 left = _graph.target((*_pred)[left]);
1783 left_path.push_back(left);
1785 left_set.insert(left);
1787 if (right_set.find(left) != right_set.end()) {
1792 if ((*_matching)[right] == INVALID) break;
1794 right = _graph.target((*_matching)[right]);
1795 right_path.push_back(right);
1796 right = _graph.target((*_pred)[right]);
1797 right_path.push_back(right);
1799 right_set.insert(right);
1803 if (nca == INVALID) {
1804 if ((*_matching)[left] == INVALID) {
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);
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);
1824 alternatePath(nca, tree);
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);
1833 prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1834 _status->set(left_path[i + 1], MATCHED);
1835 oddToMatched(left_path[i + 1]);
1837 _matching->set(nca, prev);
1839 for (int i = 0; right_path[i] != nca; i += 2) {
1840 _status->set(right_path[i], MATCHED);
1841 evenToMatched(right_path[i], tree);
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]);
1851 void extractCycle(const Arc &arc) {
1852 Node left = _graph.source(arc);
1853 Node odd = _graph.target((*_matching)[left]);
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));
1861 _matching->set(left, arc);
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));
1872 /// \brief 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),
1883 _delta2_index(0), _delta2(0),
1884 _delta3_index(0), _delta3(0),
1888 ~MaxWeightedPerfectFractionalMatching() {
1889 destroyStructures();
1892 /// \name Execution Control
1893 /// The simplest way to execute the algorithm is to use the
1894 /// \ref run() member function.
1898 /// \brief Initialize the algorithm
1900 /// This function initializes the algorithm.
1904 for (NodeIt n(_graph); n != INVALID; ++n) {
1905 (*_delta2_index)[n] = _delta2->PRE_HEAP;
1907 for (EdgeIt e(_graph); e != INVALID; ++e) {
1908 (*_delta3_index)[e] = _delta3->PRE_HEAP;
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;
1919 _node_potential->set(n, max);
1921 _tree_set->insert(n);
1923 _matching->set(n, INVALID);
1924 _status->set(n, EVEN);
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);
1937 /// \brief Start the algorithm
1939 /// This function starts the algorithm.
1941 /// \pre \ref init() must be called before using this function.
1947 int unmatched = _node_num;
1948 while (unmatched > 0) {
1949 Value d2 = !_delta2->empty() ?
1950 _delta2->prio() : std::numeric_limits<Value>::max();
1952 Value d3 = !_delta3->empty() ?
1953 _delta3->prio() : std::numeric_limits<Value>::max();
1955 _delta_sum = d3; OpType ot = D3;
1956 if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1958 if (_delta_sum == std::numeric_limits<Value>::max()) {
1965 Node n = _delta2->top();
1966 Arc a = (*_pred)[n];
1967 if ((*_matching)[n] == INVALID) {
1971 Node v = _graph.target((*_matching)[n]);
1972 if ((*_matching)[n] !=
1973 _graph.oppositeArc((*_matching)[v])) {
1983 Edge e = _delta3->top();
1985 Node left = _graph.u(e);
1986 Node right = _graph.v(e);
1988 int left_tree = _tree_set->find(left);
1989 int right_tree = _tree_set->find(right);
1991 if (left_tree == right_tree) {
1992 cycleOnEdge(e, left_tree);
2004 /// \brief Run the algorithm.
2006 /// This method runs the \c %MaxWeightedPerfectFractionalMatching
2009 /// \note mwfm.run() is just a shortcut of the following code.
2021 /// \name Primal Solution
2022 /// Functions to get the primal solution, i.e. the maximum weighted
2024 /// Either \ref run() or \ref start() function should be called before
2029 /// \brief Return the weight of the matching.
2031 /// This function returns the weight of the found matching. This
2032 /// value is scaled by \ref primalScale "primal scale".
2034 /// \pre Either run() or start() must be called before using this function.
2035 Value matchingWeight() const {
2037 for (NodeIt n(_graph); n != INVALID; ++n) {
2038 if ((*_matching)[n] != INVALID) {
2039 sum += _weight[(*_matching)[n]];
2042 return sum * primalScale / 2;
2045 /// \brief Return the number of covered nodes in the matching.
2047 /// This function returns the number of covered nodes in the matching.
2049 /// \pre Either run() or start() must be called before using this function.
2050 int matchingSize() const {
2052 for (NodeIt n(_graph); n != INVALID; ++n) {
2053 if ((*_matching)[n] != INVALID) {
2060 /// \brief Return \c true if the given edge is in the matching.
2062 /// This function returns \c true if the given edge is in the
2063 /// found matching. The result is scaled by \ref primalScale
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)
2073 /// \brief Return the fractional matching arc (or edge) incident
2074 /// to the given node.
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
2082 /// \pre Either run() or start() must be called before using this function.
2083 Arc matching(const Node& node) const {
2084 return (*_matching)[node];
2087 /// \brief Return a const reference to the matching map.
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 {
2097 /// \name Dual Solution
2098 /// Functions to get the dual solution.\n
2099 /// Either \ref run() or \ref start() function should be called before
2104 /// \brief Return the value of the dual solution.
2106 /// This function returns the value of the dual solution.
2107 /// It should be equal to the primal value scaled by \ref dualScale
2110 /// \pre Either run() or start() must be called before using this function.
2111 Value dualValue() const {
2113 for (NodeIt n(_graph); n != INVALID; ++n) {
2114 sum += nodeValue(n);
2119 /// \brief Return the dual value (potential) of the given node.
2121 /// This function returns the dual value (potential) of the given node.
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];
2132 } //END OF NAMESPACE LEMON
2134 #endif //LEMON_FRACTIONAL_MATCHING_H