gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Add fractional matching algorithms (#314)
0 4 2
default
6 files changed with 2649 insertions and 1 deletions:
↑ Collapse diff ↑
Show white space 6 line context
1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2
 *
3
 * This file is a part of LEMON, a generic C++ optimization library.
4
 *
5
 * Copyright (C) 2003-2009
6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8
 *
9
 * Permission to use, modify and distribute this software is granted
10
 * provided that this copyright notice appears in all copies. For
11
 * precise terms see the accompanying LICENSE file.
12
 *
13
 * This software is provided "AS IS" with no warranty of any kind,
14
 * express or implied, and with no claim as to its suitability for any
15
 * purpose.
16
 *
17
 */
18

	
19
#ifndef LEMON_FRACTIONAL_MATCHING_H
20
#define LEMON_FRACTIONAL_MATCHING_H
21

	
22
#include <vector>
23
#include <queue>
24
#include <set>
25
#include <limits>
26

	
27
#include <lemon/core.h>
28
#include <lemon/unionfind.h>
29
#include <lemon/bin_heap.h>
30
#include <lemon/maps.h>
31
#include <lemon/assert.h>
32
#include <lemon/elevator.h>
33

	
34
///\ingroup matching
35
///\file
36
///\brief Fractional matching algorithms in general graphs.
37

	
38
namespace lemon {
39

	
40
  /// \brief Default traits class of MaxFractionalMatching class.
41
  ///
42
  /// Default traits class of MaxFractionalMatching class.
43
  /// \tparam GR Graph type.
44
  template <typename GR>
45
  struct MaxFractionalMatchingDefaultTraits {
46

	
47
    /// \brief The type of the graph the algorithm runs on.
48
    typedef GR Graph;
49

	
50
    /// \brief The type of the map that stores the matching.
51
    ///
52
    /// The type of the map that stores the matching arcs.
53
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
54
    typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
55

	
56
    /// \brief Instantiates a MatchingMap.
57
    ///
58
    /// This function instantiates a \ref MatchingMap.
59
    /// \param graph The graph for which we would like to define
60
    /// the matching map.
61
    static MatchingMap* createMatchingMap(const Graph& graph) {
62
      return new MatchingMap(graph);
63
    }
64

	
65
    /// \brief The elevator type used by MaxFractionalMatching algorithm.
66
    ///
67
    /// The elevator type used by MaxFractionalMatching algorithm.
68
    ///
69
    /// \sa Elevator
70
    /// \sa LinkedElevator
71
    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
72

	
73
    /// \brief Instantiates an Elevator.
74
    ///
75
    /// This function instantiates an \ref Elevator.
76
    /// \param graph The graph for which we would like to define
77
    /// the elevator.
78
    /// \param max_level The maximum level of the elevator.
79
    static Elevator* createElevator(const Graph& graph, int max_level) {
80
      return new Elevator(graph, max_level);
81
    }
82
  };
83

	
84
  /// \ingroup matching
85
  ///
86
  /// \brief Max cardinality fractional matching
87
  ///
88
  /// This class provides an implementation of fractional matching
89
  /// algorithm based on push-relabel principle.
90
  ///
91
  /// The maximum cardinality fractional matching is a relaxation of the
92
  /// maximum cardinality matching problem where the odd set constraints
93
  /// are omitted.
94
  /// It can be formulated with the following linear program.
95
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
96
  /// \f[x_e \ge 0\quad \forall e\in E\f]
97
  /// \f[\max \sum_{e\in E}x_e\f]
98
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
99
  /// \f$X\f$. The result can be represented as the union of a
100
  /// matching with one value edges and a set of odd length cycles
101
  /// with half value edges.
102
  ///
103
  /// The algorithm calculates an optimal fractional matching and a
104
  /// barrier. The number of adjacents of any node set minus the size
105
  /// of node set is a lower bound on the uncovered nodes in the
106
  /// graph. For maximum matching a barrier is computed which
107
  /// maximizes this difference.
108
  ///
109
  /// The algorithm can be executed with the run() function.  After it
110
  /// the matching (the primal solution) and the barrier (the dual
111
  /// solution) can be obtained using the query functions.
112
  ///
113
  /// The primal solution is multiplied by
114
  /// \ref MaxWeightedMatching::primalScale "2".
115
  ///
116
  /// \tparam GR The undirected graph type the algorithm runs on.
117
#ifdef DOXYGEN
118
  template <typename GR, typename TR>
119
#else
120
  template <typename GR,
121
            typename TR = MaxFractionalMatchingDefaultTraits<GR> >
122
#endif
123
  class MaxFractionalMatching {
124
  public:
125

	
126
    /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
127
    /// class" of the algorithm.
128
    typedef TR Traits;
129
    /// The type of the graph the algorithm runs on.
130
    typedef typename TR::Graph Graph;
131
    /// The type of the matching map.
132
    typedef typename TR::MatchingMap MatchingMap;
133
    /// The type of the elevator.
134
    typedef typename TR::Elevator Elevator;
135

	
136
    /// \brief Scaling factor for primal solution
137
    ///
138
    /// Scaling factor for primal solution.
139
    static const int primalScale = 2;
140

	
141
  private:
142

	
143
    const Graph &_graph;
144
    int _node_num;
145
    bool _allow_loops;
146
    int _empty_level;
147

	
148
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
149

	
150
    bool _local_matching;
151
    MatchingMap *_matching;
152

	
153
    bool _local_level;
154
    Elevator *_level;
155

	
156
    typedef typename Graph::template NodeMap<int> InDegMap;
157
    InDegMap *_indeg;
158

	
159
    void createStructures() {
160
      _node_num = countNodes(_graph);
161

	
162
      if (!_matching) {
163
        _local_matching = true;
164
        _matching = Traits::createMatchingMap(_graph);
165
      }
166
      if (!_level) {
167
        _local_level = true;
168
        _level = Traits::createElevator(_graph, _node_num);
169
      }
170
      if (!_indeg) {
171
        _indeg = new InDegMap(_graph);
172
      }
173
    }
174

	
175
    void destroyStructures() {
176
      if (_local_matching) {
177
        delete _matching;
178
      }
179
      if (_local_level) {
180
        delete _level;
181
      }
182
      if (_indeg) {
183
        delete _indeg;
184
      }
185
    }
186

	
187
    void postprocessing() {
188
      for (NodeIt n(_graph); n != INVALID; ++n) {
189
        if ((*_indeg)[n] != 0) continue;
190
        _indeg->set(n, -1);
191
        Node u = n;
192
        while ((*_matching)[u] != INVALID) {
193
          Node v = _graph.target((*_matching)[u]);
194
          _indeg->set(v, -1);
195
          Arc a = _graph.oppositeArc((*_matching)[u]);
196
          u = _graph.target((*_matching)[v]);
197
          _indeg->set(u, -1);
198
          _matching->set(v, a);
199
        }
200
      }
201

	
202
      for (NodeIt n(_graph); n != INVALID; ++n) {
203
        if ((*_indeg)[n] != 1) continue;
204
        _indeg->set(n, -1);
205

	
206
        int num = 1;
207
        Node u = _graph.target((*_matching)[n]);
208
        while (u != n) {
209
          _indeg->set(u, -1);
210
          u = _graph.target((*_matching)[u]);
211
          ++num;
212
        }
213
        if (num % 2 == 0 && num > 2) {
214
          Arc prev = _graph.oppositeArc((*_matching)[n]);
215
          Node v = _graph.target((*_matching)[n]);
216
          u = _graph.target((*_matching)[v]);
217
          _matching->set(v, prev);
218
          while (u != n) {
219
            prev = _graph.oppositeArc((*_matching)[u]);
220
            v = _graph.target((*_matching)[u]);
221
            u = _graph.target((*_matching)[v]);
222
            _matching->set(v, prev);
223
          }
224
        }
225
      }
226
    }
227

	
228
  public:
229

	
230
    typedef MaxFractionalMatching Create;
231

	
232
    ///\name Named Template Parameters
233

	
234
    ///@{
235

	
236
    template <typename T>
237
    struct SetMatchingMapTraits : public Traits {
238
      typedef T MatchingMap;
239
      static MatchingMap *createMatchingMap(const Graph&) {
240
        LEMON_ASSERT(false, "MatchingMap is not initialized");
241
        return 0; // ignore warnings
242
      }
243
    };
244

	
245
    /// \brief \ref named-templ-param "Named parameter" for setting
246
    /// MatchingMap type
247
    ///
248
    /// \ref named-templ-param "Named parameter" for setting MatchingMap
249
    /// type.
250
    template <typename T>
251
    struct SetMatchingMap
252
      : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
253
      typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
254
    };
255

	
256
    template <typename T>
257
    struct SetElevatorTraits : public Traits {
258
      typedef T Elevator;
259
      static Elevator *createElevator(const Graph&, int) {
260
        LEMON_ASSERT(false, "Elevator is not initialized");
261
        return 0; // ignore warnings
262
      }
263
    };
264

	
265
    /// \brief \ref named-templ-param "Named parameter" for setting
266
    /// Elevator type
267
    ///
268
    /// \ref named-templ-param "Named parameter" for setting Elevator
269
    /// type. If this named parameter is used, then an external
270
    /// elevator object must be passed to the algorithm using the
271
    /// \ref elevator(Elevator&) "elevator()" function before calling
272
    /// \ref run() or \ref init().
273
    /// \sa SetStandardElevator
274
    template <typename T>
275
    struct SetElevator
276
      : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
277
      typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
278
    };
279

	
280
    template <typename T>
281
    struct SetStandardElevatorTraits : public Traits {
282
      typedef T Elevator;
283
      static Elevator *createElevator(const Graph& graph, int max_level) {
284
        return new Elevator(graph, max_level);
285
      }
286
    };
287

	
288
    /// \brief \ref named-templ-param "Named parameter" for setting
289
    /// Elevator type with automatic allocation
290
    ///
291
    /// \ref named-templ-param "Named parameter" for setting Elevator
292
    /// type with automatic allocation.
293
    /// The Elevator should have standard constructor interface to be
294
    /// able to automatically created by the algorithm (i.e. the
295
    /// graph and the maximum level should be passed to it).
296
    /// However an external elevator object could also be passed to the
297
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
298
    /// before calling \ref run() or \ref init().
299
    /// \sa SetElevator
300
    template <typename T>
301
    struct SetStandardElevator
302
      : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
303
      typedef MaxFractionalMatching<Graph,
304
                                    SetStandardElevatorTraits<T> > Create;
305
    };
306

	
307
    /// @}
308

	
309
  protected:
310

	
311
    MaxFractionalMatching() {}
312

	
313
  public:
314

	
315
    /// \brief Constructor
316
    ///
317
    /// Constructor.
318
    ///
319
    MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
320
      : _graph(graph), _allow_loops(allow_loops),
321
        _local_matching(false), _matching(0),
322
        _local_level(false), _level(0),  _indeg(0)
323
    {}
324

	
325
    ~MaxFractionalMatching() {
326
      destroyStructures();
327
    }
328

	
329
    /// \brief Sets the matching map.
330
    ///
331
    /// Sets the matching map.
332
    /// If you don't use this function before calling \ref run() or
333
    /// \ref init(), an instance will be allocated automatically.
334
    /// The destructor deallocates this automatically allocated map,
335
    /// of course.
336
    /// \return <tt>(*this)</tt>
337
    MaxFractionalMatching& matchingMap(MatchingMap& map) {
338
      if (_local_matching) {
339
        delete _matching;
340
        _local_matching = false;
341
      }
342
      _matching = &map;
343
      return *this;
344
    }
345

	
346
    /// \brief Sets the elevator used by algorithm.
347
    ///
348
    /// Sets the elevator used by algorithm.
349
    /// If you don't use this function before calling \ref run() or
350
    /// \ref init(), an instance will be allocated automatically.
351
    /// The destructor deallocates this automatically allocated elevator,
352
    /// of course.
353
    /// \return <tt>(*this)</tt>
354
    MaxFractionalMatching& elevator(Elevator& elevator) {
355
      if (_local_level) {
356
        delete _level;
357
        _local_level = false;
358
      }
359
      _level = &elevator;
360
      return *this;
361
    }
362

	
363
    /// \brief Returns a const reference to the elevator.
364
    ///
365
    /// Returns a const reference to the elevator.
366
    ///
367
    /// \pre Either \ref run() or \ref init() must be called before
368
    /// using this function.
369
    const Elevator& elevator() const {
370
      return *_level;
371
    }
372

	
373
    /// \name Execution control
374
    /// The simplest way to execute the algorithm is to use one of the
375
    /// member functions called \c run(). \n
376
    /// If you need more control on the execution, first
377
    /// you must call \ref init() and then one variant of the start()
378
    /// member.
379

	
380
    /// @{
381

	
382
    /// \brief Initializes the internal data structures.
383
    ///
384
    /// Initializes the internal data structures and sets the initial
385
    /// matching.
386
    void init() {
387
      createStructures();
388

	
389
      _level->initStart();
390
      for (NodeIt n(_graph); n != INVALID; ++n) {
391
        _indeg->set(n, 0);
392
        _matching->set(n, INVALID);
393
        _level->initAddItem(n);
394
      }
395
      _level->initFinish();
396

	
397
      _empty_level = _node_num;
398
      for (NodeIt n(_graph); n != INVALID; ++n) {
399
        for (OutArcIt a(_graph, n); a != INVALID; ++a) {
400
          if (_graph.target(a) == n && !_allow_loops) continue;
401
          _matching->set(n, a);
402
          Node v = _graph.target((*_matching)[n]);
403
          _indeg->set(v, (*_indeg)[v] + 1);
404
          break;
405
        }
406
      }
407

	
408
      for (NodeIt n(_graph); n != INVALID; ++n) {
409
        if ((*_indeg)[n] == 0) {
410
          _level->activate(n);
411
        }
412
      }
413
    }
414

	
415
    /// \brief Starts the algorithm and computes a fractional matching
416
    ///
417
    /// The algorithm computes a maximum fractional matching.
418
    ///
419
    /// \param postprocess The algorithm computes first a matching
420
    /// which is a union of a matching with one value edges, cycles
421
    /// with half value edges and even length paths with half value
422
    /// edges. If the parameter is true, then after the push-relabel
423
    /// algorithm it postprocesses the matching to contain only
424
    /// matching edges and half value odd cycles.
425
    void start(bool postprocess = true) {
426
      Node n;
427
      while ((n = _level->highestActive()) != INVALID) {
428
        int level = _level->highestActiveLevel();
429
        int new_level = _level->maxLevel();
430
        for (InArcIt a(_graph, n); a != INVALID; ++a) {
431
          Node u = _graph.source(a);
432
          if (n == u && !_allow_loops) continue;
433
          Node v = _graph.target((*_matching)[u]);
434
          if ((*_level)[v] < level) {
435
            _indeg->set(v, (*_indeg)[v] - 1);
436
            if ((*_indeg)[v] == 0) {
437
              _level->activate(v);
438
            }
439
            _matching->set(u, a);
440
            _indeg->set(n, (*_indeg)[n] + 1);
441
            _level->deactivate(n);
442
            goto no_more_push;
443
          } else if (new_level > (*_level)[v]) {
444
            new_level = (*_level)[v];
445
          }
446
        }
447

	
448
        if (new_level + 1 < _level->maxLevel()) {
449
          _level->liftHighestActive(new_level + 1);
450
        } else {
451
          _level->liftHighestActiveToTop();
452
        }
453
        if (_level->emptyLevel(level)) {
454
          _level->liftToTop(level);
455
        }
456
      no_more_push:
457
        ;
458
      }
459
      for (NodeIt n(_graph); n != INVALID; ++n) {
460
        if ((*_matching)[n] == INVALID) continue;
461
        Node u = _graph.target((*_matching)[n]);
462
        if ((*_indeg)[u] > 1) {
463
          _indeg->set(u, (*_indeg)[u] - 1);
464
          _matching->set(n, INVALID);
465
        }
466
      }
467
      if (postprocess) {
468
        postprocessing();
469
      }
470
    }
471

	
472
    /// \brief Starts the algorithm and computes a perfect fractional
473
    /// matching
474
    ///
475
    /// The algorithm computes a perfect fractional matching. If it
476
    /// does not exists, then the algorithm returns false and the
477
    /// matching is undefined and the barrier.
478
    ///
479
    /// \param postprocess The algorithm computes first a matching
480
    /// which is a union of a matching with one value edges, cycles
481
    /// with half value edges and even length paths with half value
482
    /// edges. If the parameter is true, then after the push-relabel
483
    /// algorithm it postprocesses the matching to contain only
484
    /// matching edges and half value odd cycles.
485
    bool startPerfect(bool postprocess = true) {
486
      Node n;
487
      while ((n = _level->highestActive()) != INVALID) {
488
        int level = _level->highestActiveLevel();
489
        int new_level = _level->maxLevel();
490
        for (InArcIt a(_graph, n); a != INVALID; ++a) {
491
          Node u = _graph.source(a);
492
          if (n == u && !_allow_loops) continue;
493
          Node v = _graph.target((*_matching)[u]);
494
          if ((*_level)[v] < level) {
495
            _indeg->set(v, (*_indeg)[v] - 1);
496
            if ((*_indeg)[v] == 0) {
497
              _level->activate(v);
498
            }
499
            _matching->set(u, a);
500
            _indeg->set(n, (*_indeg)[n] + 1);
501
            _level->deactivate(n);
502
            goto no_more_push;
503
          } else if (new_level > (*_level)[v]) {
504
            new_level = (*_level)[v];
505
          }
506
        }
507

	
508
        if (new_level + 1 < _level->maxLevel()) {
509
          _level->liftHighestActive(new_level + 1);
510
        } else {
511
          _level->liftHighestActiveToTop();
512
          _empty_level = _level->maxLevel() - 1;
513
          return false;
514
        }
515
        if (_level->emptyLevel(level)) {
516
          _level->liftToTop(level);
517
          _empty_level = level;
518
          return false;
519
        }
520
      no_more_push:
521
        ;
522
      }
523
      if (postprocess) {
524
        postprocessing();
525
      }
526
      return true;
527
    }
528

	
529
    /// \brief Runs the algorithm
530
    ///
531
    /// Just a shortcut for the next code:
532
    ///\code
533
    /// init();
534
    /// start();
535
    ///\endcode
536
    void run(bool postprocess = true) {
537
      init();
538
      start(postprocess);
539
    }
540

	
541
    /// \brief Runs the algorithm to find a perfect fractional matching
542
    ///
543
    /// Just a shortcut for the next code:
544
    ///\code
545
    /// init();
546
    /// startPerfect();
547
    ///\endcode
548
    bool runPerfect(bool postprocess = true) {
549
      init();
550
      return startPerfect(postprocess);
551
    }
552

	
553
    ///@}
554

	
555
    /// \name Query Functions
556
    /// The result of the %Matching algorithm can be obtained using these
557
    /// functions.\n
558
    /// Before the use of these functions,
559
    /// either run() or start() must be called.
560
    ///@{
561

	
562

	
563
    /// \brief Return the number of covered nodes in the matching.
564
    ///
565
    /// This function returns the number of covered nodes in the matching.
566
    ///
567
    /// \pre Either run() or start() must be called before using this function.
568
    int matchingSize() const {
569
      int num = 0;
570
      for (NodeIt n(_graph); n != INVALID; ++n) {
571
        if ((*_matching)[n] != INVALID) {
572
          ++num;
573
        }
574
      }
575
      return num;
576
    }
577

	
578
    /// \brief Returns a const reference to the matching map.
579
    ///
580
    /// Returns a const reference to the node map storing the found
581
    /// fractional matching. This method can be called after
582
    /// running the algorithm.
583
    ///
584
    /// \pre Either \ref run() or \ref init() must be called before
585
    /// using this function.
586
    const MatchingMap& matchingMap() const {
587
      return *_matching;
588
    }
589

	
590
    /// \brief Return \c true if the given edge is in the matching.
591
    ///
592
    /// This function returns \c true if the given edge is in the
593
    /// found matching. The result is scaled by \ref primalScale
594
    /// "primal scale".
595
    ///
596
    /// \pre Either run() or start() must be called before using this function.
597
    int matching(const Edge& edge) const {
598
      return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
599
        (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
600
    }
601

	
602
    /// \brief Return the fractional matching arc (or edge) incident
603
    /// to the given node.
604
    ///
605
    /// This function returns one of the fractional matching arc (or
606
    /// edge) incident to the given node in the found matching or \c
607
    /// INVALID if the node is not covered by the matching or if the
608
    /// node is on an odd length cycle then it is the successor edge
609
    /// on the cycle.
610
    ///
611
    /// \pre Either run() or start() must be called before using this function.
612
    Arc matching(const Node& node) const {
613
      return (*_matching)[node];
614
    }
615

	
616
    /// \brief Returns true if the node is in the barrier
617
    ///
618
    /// The barrier is a subset of the nodes. If the nodes in the
619
    /// barrier have less adjacent nodes than the size of the barrier,
620
    /// then at least as much nodes cannot be covered as the
621
    /// difference of the two subsets.
622
    bool barrier(const Node& node) const {
623
      return (*_level)[node] >= _empty_level;
624
    }
625

	
626
    /// @}
627

	
628
  };
629

	
630
  /// \ingroup matching
631
  ///
632
  /// \brief Weighted fractional matching in general graphs
633
  ///
634
  /// This class provides an efficient implementation of fractional
635
  /// matching algorithm. The implementation is based on extensive use
636
  /// of priority queues and provides \f$O(nm\log n)\f$ time
637
  /// complexity.
638
  ///
639
  /// The maximum weighted fractional matching is a relaxation of the
640
  /// maximum weighted matching problem where the odd set constraints
641
  /// are omitted.
642
  /// It can be formulated with the following linear program.
643
  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
644
  /// \f[x_e \ge 0\quad \forall e\in E\f]
645
  /// \f[\max \sum_{e\in E}x_ew_e\f]
646
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
647
  /// \f$X\f$. The result must be the union of a matching with one
648
  /// value edges and a set of odd length cycles with half value edges.
649
  ///
650
  /// The algorithm calculates an optimal fractional matching and a
651
  /// proof of the optimality. The solution of the dual problem can be
652
  /// used to check the result of the algorithm. The dual linear
653
  /// problem is the following.
654
  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
655
  /// \f[y_u \ge 0 \quad \forall u \in V\f]
656
  /// \f[\min \sum_{u \in V}y_u \f] ///
657
  ///
658
  /// The algorithm can be executed with the run() function.
659
  /// After it the matching (the primal solution) and the dual solution
660
  /// can be obtained using the query functions.
661
  ///
662
  /// If the value type is integer, then the primal and the dual
663
  /// solutions are multiplied by
664
  /// \ref MaxWeightedMatching::primalScale "2" and
665
  /// \ref MaxWeightedMatching::dualScale "4" respectively.
666
  ///
667
  /// \tparam GR The undirected graph type the algorithm runs on.
668
  /// \tparam WM The type edge weight map. The default type is
669
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
670
#ifdef DOXYGEN
671
  template <typename GR, typename WM>
672
#else
673
  template <typename GR,
674
            typename WM = typename GR::template EdgeMap<int> >
675
#endif
676
  class MaxWeightedFractionalMatching {
677
  public:
678

	
679
    /// The graph type of the algorithm
680
    typedef GR Graph;
681
    /// The type of the edge weight map
682
    typedef WM WeightMap;
683
    /// The value type of the edge weights
684
    typedef typename WeightMap::Value Value;
685

	
686
    /// The type of the matching map
687
    typedef typename Graph::template NodeMap<typename Graph::Arc>
688
    MatchingMap;
689

	
690
    /// \brief Scaling factor for primal solution
691
    ///
692
    /// Scaling factor for primal solution. It is equal to 2 or 1
693
    /// according to the value type.
694
    static const int primalScale =
695
      std::numeric_limits<Value>::is_integer ? 2 : 1;
696

	
697
    /// \brief Scaling factor for dual solution
698
    ///
699
    /// Scaling factor for dual solution. It is equal to 4 or 1
700
    /// according to the value type.
701
    static const int dualScale =
702
      std::numeric_limits<Value>::is_integer ? 4 : 1;
703

	
704
  private:
705

	
706
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
707

	
708
    typedef typename Graph::template NodeMap<Value> NodePotential;
709

	
710
    const Graph& _graph;
711
    const WeightMap& _weight;
712

	
713
    MatchingMap* _matching;
714
    NodePotential* _node_potential;
715

	
716
    int _node_num;
717
    bool _allow_loops;
718

	
719
    enum Status {
720
      EVEN = -1, MATCHED = 0, ODD = 1
721
    };
722

	
723
    typedef typename Graph::template NodeMap<Status> StatusMap;
724
    StatusMap* _status;
725

	
726
    typedef typename Graph::template NodeMap<Arc> PredMap;
727
    PredMap* _pred;
728

	
729
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
730

	
731
    IntNodeMap *_tree_set_index;
732
    TreeSet *_tree_set;
733

	
734
    IntNodeMap *_delta1_index;
735
    BinHeap<Value, IntNodeMap> *_delta1;
736

	
737
    IntNodeMap *_delta2_index;
738
    BinHeap<Value, IntNodeMap> *_delta2;
739

	
740
    IntEdgeMap *_delta3_index;
741
    BinHeap<Value, IntEdgeMap> *_delta3;
742

	
743
    Value _delta_sum;
744

	
745
    void createStructures() {
746
      _node_num = countNodes(_graph);
747

	
748
      if (!_matching) {
749
        _matching = new MatchingMap(_graph);
750
      }
751
      if (!_node_potential) {
752
        _node_potential = new NodePotential(_graph);
753
      }
754
      if (!_status) {
755
        _status = new StatusMap(_graph);
756
      }
757
      if (!_pred) {
758
        _pred = new PredMap(_graph);
759
      }
760
      if (!_tree_set) {
761
        _tree_set_index = new IntNodeMap(_graph);
762
        _tree_set = new TreeSet(*_tree_set_index);
763
      }
764
      if (!_delta1) {
765
        _delta1_index = new IntNodeMap(_graph);
766
        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
767
      }
768
      if (!_delta2) {
769
        _delta2_index = new IntNodeMap(_graph);
770
        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
771
      }
772
      if (!_delta3) {
773
        _delta3_index = new IntEdgeMap(_graph);
774
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
775
      }
776
    }
777

	
778
    void destroyStructures() {
779
      if (_matching) {
780
        delete _matching;
781
      }
782
      if (_node_potential) {
783
        delete _node_potential;
784
      }
785
      if (_status) {
786
        delete _status;
787
      }
788
      if (_pred) {
789
        delete _pred;
790
      }
791
      if (_tree_set) {
792
        delete _tree_set_index;
793
        delete _tree_set;
794
      }
795
      if (_delta1) {
796
        delete _delta1_index;
797
        delete _delta1;
798
      }
799
      if (_delta2) {
800
        delete _delta2_index;
801
        delete _delta2;
802
      }
803
      if (_delta3) {
804
        delete _delta3_index;
805
        delete _delta3;
806
      }
807
    }
808

	
809
    void matchedToEven(Node node, int tree) {
810
      _tree_set->insert(node, tree);
811
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
812
      _delta1->push(node, (*_node_potential)[node]);
813

	
814
      if (_delta2->state(node) == _delta2->IN_HEAP) {
815
        _delta2->erase(node);
816
      }
817

	
818
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
819
        Node v = _graph.source(a);
820
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
821
          dualScale * _weight[a];
822
        if (node == v) {
823
          if (_allow_loops && _graph.direction(a)) {
824
            _delta3->push(a, rw / 2);
825
          }
826
        } else if ((*_status)[v] == EVEN) {
827
          _delta3->push(a, rw / 2);
828
        } else if ((*_status)[v] == MATCHED) {
829
          if (_delta2->state(v) != _delta2->IN_HEAP) {
830
            _pred->set(v, a);
831
            _delta2->push(v, rw);
832
          } else if ((*_delta2)[v] > rw) {
833
            _pred->set(v, a);
834
            _delta2->decrease(v, rw);
835
          }
836
        }
837
      }
838
    }
839

	
840
    void matchedToOdd(Node node, int tree) {
841
      _tree_set->insert(node, tree);
842
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
843

	
844
      if (_delta2->state(node) == _delta2->IN_HEAP) {
845
        _delta2->erase(node);
846
      }
847
    }
848

	
849
    void evenToMatched(Node node, int tree) {
850
      _delta1->erase(node);
851
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
852
      Arc min = INVALID;
853
      Value minrw = std::numeric_limits<Value>::max();
854
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
855
        Node v = _graph.source(a);
856
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
857
          dualScale * _weight[a];
858

	
859
        if (node == v) {
860
          if (_allow_loops && _graph.direction(a)) {
861
            _delta3->erase(a);
862
          }
863
        } else if ((*_status)[v] == EVEN) {
864
          _delta3->erase(a);
865
          if (minrw > rw) {
866
            min = _graph.oppositeArc(a);
867
            minrw = rw;
868
          }
869
        } else if ((*_status)[v]  == MATCHED) {
870
          if ((*_pred)[v] == a) {
871
            Arc mina = INVALID;
872
            Value minrwa = std::numeric_limits<Value>::max();
873
            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
874
              Node va = _graph.target(aa);
875
              if ((*_status)[va] != EVEN ||
876
                  _tree_set->find(va) == tree) continue;
877
              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
878
                dualScale * _weight[aa];
879
              if (minrwa > rwa) {
880
                minrwa = rwa;
881
                mina = aa;
882
              }
883
            }
884
            if (mina != INVALID) {
885
              _pred->set(v, mina);
886
              _delta2->increase(v, minrwa);
887
            } else {
888
              _pred->set(v, INVALID);
889
              _delta2->erase(v);
890
            }
891
          }
892
        }
893
      }
894
      if (min != INVALID) {
895
        _pred->set(node, min);
896
        _delta2->push(node, minrw);
897
      } else {
898
        _pred->set(node, INVALID);
899
      }
900
    }
901

	
902
    void oddToMatched(Node node) {
903
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
904
      Arc min = INVALID;
905
      Value minrw = std::numeric_limits<Value>::max();
906
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
907
        Node v = _graph.source(a);
908
        if ((*_status)[v] != EVEN) continue;
909
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
910
          dualScale * _weight[a];
911

	
912
        if (minrw > rw) {
913
          min = _graph.oppositeArc(a);
914
          minrw = rw;
915
        }
916
      }
917
      if (min != INVALID) {
918
        _pred->set(node, min);
919
        _delta2->push(node, minrw);
920
      } else {
921
        _pred->set(node, INVALID);
922
      }
923
    }
924

	
925
    void alternatePath(Node even, int tree) {
926
      Node odd;
927

	
928
      _status->set(even, MATCHED);
929
      evenToMatched(even, tree);
930

	
931
      Arc prev = (*_matching)[even];
932
      while (prev != INVALID) {
933
        odd = _graph.target(prev);
934
        even = _graph.target((*_pred)[odd]);
935
        _matching->set(odd, (*_pred)[odd]);
936
        _status->set(odd, MATCHED);
937
        oddToMatched(odd);
938

	
939
        prev = (*_matching)[even];
940
        _status->set(even, MATCHED);
941
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
942
        evenToMatched(even, tree);
943
      }
944
    }
945

	
946
    void destroyTree(int tree) {
947
      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
948
        if ((*_status)[n] == EVEN) {
949
          _status->set(n, MATCHED);
950
          evenToMatched(n, tree);
951
        } else if ((*_status)[n] == ODD) {
952
          _status->set(n, MATCHED);
953
          oddToMatched(n);
954
        }
955
      }
956
      _tree_set->eraseClass(tree);
957
    }
958

	
959

	
960
    void unmatchNode(const Node& node) {
961
      int tree = _tree_set->find(node);
962

	
963
      alternatePath(node, tree);
964
      destroyTree(tree);
965

	
966
      _matching->set(node, INVALID);
967
    }
968

	
969

	
970
    void augmentOnEdge(const Edge& edge) {
971
      Node left = _graph.u(edge);
972
      int left_tree = _tree_set->find(left);
973

	
974
      alternatePath(left, left_tree);
975
      destroyTree(left_tree);
976
      _matching->set(left, _graph.direct(edge, true));
977

	
978
      Node right = _graph.v(edge);
979
      int right_tree = _tree_set->find(right);
980

	
981
      alternatePath(right, right_tree);
982
      destroyTree(right_tree);
983
      _matching->set(right, _graph.direct(edge, false));
984
    }
985

	
986
    void augmentOnArc(const Arc& arc) {
987
      Node left = _graph.source(arc);
988
      _status->set(left, MATCHED);
989
      _matching->set(left, arc);
990
      _pred->set(left, arc);
991

	
992
      Node right = _graph.target(arc);
993
      int right_tree = _tree_set->find(right);
994

	
995
      alternatePath(right, right_tree);
996
      destroyTree(right_tree);
997
      _matching->set(right, _graph.oppositeArc(arc));
998
    }
999

	
1000
    void extendOnArc(const Arc& arc) {
1001
      Node base = _graph.target(arc);
1002
      int tree = _tree_set->find(base);
1003

	
1004
      Node odd = _graph.source(arc);
1005
      _tree_set->insert(odd, tree);
1006
      _status->set(odd, ODD);
1007
      matchedToOdd(odd, tree);
1008
      _pred->set(odd, arc);
1009

	
1010
      Node even = _graph.target((*_matching)[odd]);
1011
      _tree_set->insert(even, tree);
1012
      _status->set(even, EVEN);
1013
      matchedToEven(even, tree);
1014
    }
1015

	
1016
    void cycleOnEdge(const Edge& edge, int tree) {
1017
      Node nca = INVALID;
1018
      std::vector<Node> left_path, right_path;
1019

	
1020
      {
1021
        std::set<Node> left_set, right_set;
1022
        Node left = _graph.u(edge);
1023
        left_path.push_back(left);
1024
        left_set.insert(left);
1025

	
1026
        Node right = _graph.v(edge);
1027
        right_path.push_back(right);
1028
        right_set.insert(right);
1029

	
1030
        while (true) {
1031

	
1032
          if (left_set.find(right) != left_set.end()) {
1033
            nca = right;
1034
            break;
1035
          }
1036

	
1037
          if ((*_matching)[left] == INVALID) break;
1038

	
1039
          left = _graph.target((*_matching)[left]);
1040
          left_path.push_back(left);
1041
          left = _graph.target((*_pred)[left]);
1042
          left_path.push_back(left);
1043

	
1044
          left_set.insert(left);
1045

	
1046
          if (right_set.find(left) != right_set.end()) {
1047
            nca = left;
1048
            break;
1049
          }
1050

	
1051
          if ((*_matching)[right] == INVALID) break;
1052

	
1053
          right = _graph.target((*_matching)[right]);
1054
          right_path.push_back(right);
1055
          right = _graph.target((*_pred)[right]);
1056
          right_path.push_back(right);
1057

	
1058
          right_set.insert(right);
1059

	
1060
        }
1061

	
1062
        if (nca == INVALID) {
1063
          if ((*_matching)[left] == INVALID) {
1064
            nca = right;
1065
            while (left_set.find(nca) == left_set.end()) {
1066
              nca = _graph.target((*_matching)[nca]);
1067
              right_path.push_back(nca);
1068
              nca = _graph.target((*_pred)[nca]);
1069
              right_path.push_back(nca);
1070
            }
1071
          } else {
1072
            nca = left;
1073
            while (right_set.find(nca) == right_set.end()) {
1074
              nca = _graph.target((*_matching)[nca]);
1075
              left_path.push_back(nca);
1076
              nca = _graph.target((*_pred)[nca]);
1077
              left_path.push_back(nca);
1078
            }
1079
          }
1080
        }
1081
      }
1082

	
1083
      alternatePath(nca, tree);
1084
      Arc prev;
1085

	
1086
      prev = _graph.direct(edge, true);
1087
      for (int i = 0; left_path[i] != nca; i += 2) {
1088
        _matching->set(left_path[i], prev);
1089
        _status->set(left_path[i], MATCHED);
1090
        evenToMatched(left_path[i], tree);
1091

	
1092
        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1093
        _status->set(left_path[i + 1], MATCHED);
1094
        oddToMatched(left_path[i + 1]);
1095
      }
1096
      _matching->set(nca, prev);
1097

	
1098
      for (int i = 0; right_path[i] != nca; i += 2) {
1099
        _status->set(right_path[i], MATCHED);
1100
        evenToMatched(right_path[i], tree);
1101

	
1102
        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1103
        _status->set(right_path[i + 1], MATCHED);
1104
        oddToMatched(right_path[i + 1]);
1105
      }
1106

	
1107
      destroyTree(tree);
1108
    }
1109

	
1110
    void extractCycle(const Arc &arc) {
1111
      Node left = _graph.source(arc);
1112
      Node odd = _graph.target((*_matching)[left]);
1113
      Arc prev;
1114
      while (odd != left) {
1115
        Node even = _graph.target((*_matching)[odd]);
1116
        prev = (*_matching)[odd];
1117
        odd = _graph.target((*_matching)[even]);
1118
        _matching->set(even, _graph.oppositeArc(prev));
1119
      }
1120
      _matching->set(left, arc);
1121

	
1122
      Node right = _graph.target(arc);
1123
      int right_tree = _tree_set->find(right);
1124
      alternatePath(right, right_tree);
1125
      destroyTree(right_tree);
1126
      _matching->set(right, _graph.oppositeArc(arc));
1127
    }
1128

	
1129
  public:
1130

	
1131
    /// \brief Constructor
1132
    ///
1133
    /// Constructor.
1134
    MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
1135
                                  bool allow_loops = true)
1136
      : _graph(graph), _weight(weight), _matching(0),
1137
      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1138
      _status(0),  _pred(0),
1139
      _tree_set_index(0), _tree_set(0),
1140

	
1141
      _delta1_index(0), _delta1(0),
1142
      _delta2_index(0), _delta2(0),
1143
      _delta3_index(0), _delta3(0),
1144

	
1145
      _delta_sum() {}
1146

	
1147
    ~MaxWeightedFractionalMatching() {
1148
      destroyStructures();
1149
    }
1150

	
1151
    /// \name Execution Control
1152
    /// The simplest way to execute the algorithm is to use the
1153
    /// \ref run() member function.
1154

	
1155
    ///@{
1156

	
1157
    /// \brief Initialize the algorithm
1158
    ///
1159
    /// This function initializes the algorithm.
1160
    void init() {
1161
      createStructures();
1162

	
1163
      for (NodeIt n(_graph); n != INVALID; ++n) {
1164
        (*_delta1_index)[n] = _delta1->PRE_HEAP;
1165
        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1166
      }
1167
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1168
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1169
      }
1170

	
1171
      for (NodeIt n(_graph); n != INVALID; ++n) {
1172
        Value max = 0;
1173
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1174
          if (_graph.target(e) == n && !_allow_loops) continue;
1175
          if ((dualScale * _weight[e]) / 2 > max) {
1176
            max = (dualScale * _weight[e]) / 2;
1177
          }
1178
        }
1179
        _node_potential->set(n, max);
1180
        _delta1->push(n, max);
1181

	
1182
        _tree_set->insert(n);
1183

	
1184
        _matching->set(n, INVALID);
1185
        _status->set(n, EVEN);
1186
      }
1187

	
1188
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1189
        Node left = _graph.u(e);
1190
        Node right = _graph.v(e);
1191
        if (left == right && !_allow_loops) continue;
1192
        _delta3->push(e, ((*_node_potential)[left] +
1193
                          (*_node_potential)[right] -
1194
                          dualScale * _weight[e]) / 2);
1195
      }
1196
    }
1197

	
1198
    /// \brief Start the algorithm
1199
    ///
1200
    /// This function starts the algorithm.
1201
    ///
1202
    /// \pre \ref init() must be called before using this function.
1203
    void start() {
1204
      enum OpType {
1205
        D1, D2, D3
1206
      };
1207

	
1208
      int unmatched = _node_num;
1209
      while (unmatched > 0) {
1210
        Value d1 = !_delta1->empty() ?
1211
          _delta1->prio() : std::numeric_limits<Value>::max();
1212

	
1213
        Value d2 = !_delta2->empty() ?
1214
          _delta2->prio() : std::numeric_limits<Value>::max();
1215

	
1216
        Value d3 = !_delta3->empty() ?
1217
          _delta3->prio() : std::numeric_limits<Value>::max();
1218

	
1219
        _delta_sum = d3; OpType ot = D3;
1220
        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
1221
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1222

	
1223
        switch (ot) {
1224
        case D1:
1225
          {
1226
            Node n = _delta1->top();
1227
            unmatchNode(n);
1228
            --unmatched;
1229
          }
1230
          break;
1231
        case D2:
1232
          {
1233
            Node n = _delta2->top();
1234
            Arc a = (*_pred)[n];
1235
            if ((*_matching)[n] == INVALID) {
1236
              augmentOnArc(a);
1237
              --unmatched;
1238
            } else {
1239
              Node v = _graph.target((*_matching)[n]);
1240
              if ((*_matching)[n] !=
1241
                  _graph.oppositeArc((*_matching)[v])) {
1242
                extractCycle(a);
1243
                --unmatched;
1244
              } else {
1245
                extendOnArc(a);
1246
              }
1247
            }
1248
          } break;
1249
        case D3:
1250
          {
1251
            Edge e = _delta3->top();
1252

	
1253
            Node left = _graph.u(e);
1254
            Node right = _graph.v(e);
1255

	
1256
            int left_tree = _tree_set->find(left);
1257
            int right_tree = _tree_set->find(right);
1258

	
1259
            if (left_tree == right_tree) {
1260
              cycleOnEdge(e, left_tree);
1261
              --unmatched;
1262
            } else {
1263
              augmentOnEdge(e);
1264
              unmatched -= 2;
1265
            }
1266
          } break;
1267
        }
1268
      }
1269
    }
1270

	
1271
    /// \brief Run the algorithm.
1272
    ///
1273
    /// This method runs the \c %MaxWeightedMatching algorithm.
1274
    ///
1275
    /// \note mwfm.run() is just a shortcut of the following code.
1276
    /// \code
1277
    ///   mwfm.init();
1278
    ///   mwfm.start();
1279
    /// \endcode
1280
    void run() {
1281
      init();
1282
      start();
1283
    }
1284

	
1285
    /// @}
1286

	
1287
    /// \name Primal Solution
1288
    /// Functions to get the primal solution, i.e. the maximum weighted
1289
    /// matching.\n
1290
    /// Either \ref run() or \ref start() function should be called before
1291
    /// using them.
1292

	
1293
    /// @{
1294

	
1295
    /// \brief Return the weight of the matching.
1296
    ///
1297
    /// This function returns the weight of the found matching. This
1298
    /// value is scaled by \ref primalScale "primal scale".
1299
    ///
1300
    /// \pre Either run() or start() must be called before using this function.
1301
    Value matchingWeight() const {
1302
      Value sum = 0;
1303
      for (NodeIt n(_graph); n != INVALID; ++n) {
1304
        if ((*_matching)[n] != INVALID) {
1305
          sum += _weight[(*_matching)[n]];
1306
        }
1307
      }
1308
      return sum * primalScale / 2;
1309
    }
1310

	
1311
    /// \brief Return the number of covered nodes in the matching.
1312
    ///
1313
    /// This function returns the number of covered nodes in the matching.
1314
    ///
1315
    /// \pre Either run() or start() must be called before using this function.
1316
    int matchingSize() const {
1317
      int num = 0;
1318
      for (NodeIt n(_graph); n != INVALID; ++n) {
1319
        if ((*_matching)[n] != INVALID) {
1320
          ++num;
1321
        }
1322
      }
1323
      return num;
1324
    }
1325

	
1326
    /// \brief Return \c true if the given edge is in the matching.
1327
    ///
1328
    /// This function returns \c true if the given edge is in the
1329
    /// found matching. The result is scaled by \ref primalScale
1330
    /// "primal scale".
1331
    ///
1332
    /// \pre Either run() or start() must be called before using this function.
1333
    Value matching(const Edge& edge) const {
1334
      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
1335
        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
1336
        * primalScale / 2;
1337
    }
1338

	
1339
    /// \brief Return the fractional matching arc (or edge) incident
1340
    /// to the given node.
1341
    ///
1342
    /// This function returns one of the fractional matching arc (or
1343
    /// edge) incident to the given node in the found matching or \c
1344
    /// INVALID if the node is not covered by the matching or if the
1345
    /// node is on an odd length cycle then it is the successor edge
1346
    /// on the cycle.
1347
    ///
1348
    /// \pre Either run() or start() must be called before using this function.
1349
    Arc matching(const Node& node) const {
1350
      return (*_matching)[node];
1351
    }
1352

	
1353
    /// \brief Return a const reference to the matching map.
1354
    ///
1355
    /// This function returns a const reference to a node map that stores
1356
    /// the matching arc (or edge) incident to each node.
1357
    const MatchingMap& matchingMap() const {
1358
      return *_matching;
1359
    }
1360

	
1361
    /// @}
1362

	
1363
    /// \name Dual Solution
1364
    /// Functions to get the dual solution.\n
1365
    /// Either \ref run() or \ref start() function should be called before
1366
    /// using them.
1367

	
1368
    /// @{
1369

	
1370
    /// \brief Return the value of the dual solution.
1371
    ///
1372
    /// This function returns the value of the dual solution.
1373
    /// It should be equal to the primal value scaled by \ref dualScale
1374
    /// "dual scale".
1375
    ///
1376
    /// \pre Either run() or start() must be called before using this function.
1377
    Value dualValue() const {
1378
      Value sum = 0;
1379
      for (NodeIt n(_graph); n != INVALID; ++n) {
1380
        sum += nodeValue(n);
1381
      }
1382
      return sum;
1383
    }
1384

	
1385
    /// \brief Return the dual value (potential) of the given node.
1386
    ///
1387
    /// This function returns the dual value (potential) of the given node.
1388
    ///
1389
    /// \pre Either run() or start() must be called before using this function.
1390
    Value nodeValue(const Node& n) const {
1391
      return (*_node_potential)[n];
1392
    }
1393

	
1394
    /// @}
1395

	
1396
  };
1397

	
1398
  /// \ingroup matching
1399
  ///
1400
  /// \brief Weighted fractional perfect matching in general graphs
1401
  ///
1402
  /// This class provides an efficient implementation of fractional
1403
  /// matching algorithm. The implementation is based on extensive use
1404
  /// of priority queues and provides \f$O(nm\log n)\f$ time
1405
  /// complexity.
1406
  ///
1407
  /// The maximum weighted fractional perfect matching is a relaxation
1408
  /// of the maximum weighted perfect matching problem where the odd
1409
  /// set constraints are omitted.
1410
  /// It can be formulated with the following linear program.
1411
  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
1412
  /// \f[x_e \ge 0\quad \forall e\in E\f]
1413
  /// \f[\max \sum_{e\in E}x_ew_e\f]
1414
  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
1415
  /// \f$X\f$. The result must be the union of a matching with one
1416
  /// value edges and a set of odd length cycles with half value edges.
1417
  ///
1418
  /// The algorithm calculates an optimal fractional matching and a
1419
  /// proof of the optimality. The solution of the dual problem can be
1420
  /// used to check the result of the algorithm. The dual linear
1421
  /// problem is the following.
1422
  /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
1423
  /// \f[\min \sum_{u \in V}y_u \f] ///
1424
  ///
1425
  /// The algorithm can be executed with the run() function.
1426
  /// After it the matching (the primal solution) and the dual solution
1427
  /// can be obtained using the query functions.
1428

	
1429
  /// If the value type is integer, then the primal and the dual
1430
  /// solutions are multiplied by
1431
  /// \ref MaxWeightedMatching::primalScale "2" and
1432
  /// \ref MaxWeightedMatching::dualScale "4" respectively.
1433
  ///
1434
  /// \tparam GR The undirected graph type the algorithm runs on.
1435
  /// \tparam WM The type edge weight map. The default type is
1436
  /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1437
#ifdef DOXYGEN
1438
  template <typename GR, typename WM>
1439
#else
1440
  template <typename GR,
1441
            typename WM = typename GR::template EdgeMap<int> >
1442
#endif
1443
  class MaxWeightedPerfectFractionalMatching {
1444
  public:
1445

	
1446
    /// The graph type of the algorithm
1447
    typedef GR Graph;
1448
    /// The type of the edge weight map
1449
    typedef WM WeightMap;
1450
    /// The value type of the edge weights
1451
    typedef typename WeightMap::Value Value;
1452

	
1453
    /// The type of the matching map
1454
    typedef typename Graph::template NodeMap<typename Graph::Arc>
1455
    MatchingMap;
1456

	
1457
    /// \brief Scaling factor for primal solution
1458
    ///
1459
    /// Scaling factor for primal solution. It is equal to 2 or 1
1460
    /// according to the value type.
1461
    static const int primalScale =
1462
      std::numeric_limits<Value>::is_integer ? 2 : 1;
1463

	
1464
    /// \brief Scaling factor for dual solution
1465
    ///
1466
    /// Scaling factor for dual solution. It is equal to 4 or 1
1467
    /// according to the value type.
1468
    static const int dualScale =
1469
      std::numeric_limits<Value>::is_integer ? 4 : 1;
1470

	
1471
  private:
1472

	
1473
    TEMPLATE_GRAPH_TYPEDEFS(Graph);
1474

	
1475
    typedef typename Graph::template NodeMap<Value> NodePotential;
1476

	
1477
    const Graph& _graph;
1478
    const WeightMap& _weight;
1479

	
1480
    MatchingMap* _matching;
1481
    NodePotential* _node_potential;
1482

	
1483
    int _node_num;
1484
    bool _allow_loops;
1485

	
1486
    enum Status {
1487
      EVEN = -1, MATCHED = 0, ODD = 1
1488
    };
1489

	
1490
    typedef typename Graph::template NodeMap<Status> StatusMap;
1491
    StatusMap* _status;
1492

	
1493
    typedef typename Graph::template NodeMap<Arc> PredMap;
1494
    PredMap* _pred;
1495

	
1496
    typedef ExtendFindEnum<IntNodeMap> TreeSet;
1497

	
1498
    IntNodeMap *_tree_set_index;
1499
    TreeSet *_tree_set;
1500

	
1501
    IntNodeMap *_delta2_index;
1502
    BinHeap<Value, IntNodeMap> *_delta2;
1503

	
1504
    IntEdgeMap *_delta3_index;
1505
    BinHeap<Value, IntEdgeMap> *_delta3;
1506

	
1507
    Value _delta_sum;
1508

	
1509
    void createStructures() {
1510
      _node_num = countNodes(_graph);
1511

	
1512
      if (!_matching) {
1513
        _matching = new MatchingMap(_graph);
1514
      }
1515
      if (!_node_potential) {
1516
        _node_potential = new NodePotential(_graph);
1517
      }
1518
      if (!_status) {
1519
        _status = new StatusMap(_graph);
1520
      }
1521
      if (!_pred) {
1522
        _pred = new PredMap(_graph);
1523
      }
1524
      if (!_tree_set) {
1525
        _tree_set_index = new IntNodeMap(_graph);
1526
        _tree_set = new TreeSet(*_tree_set_index);
1527
      }
1528
      if (!_delta2) {
1529
        _delta2_index = new IntNodeMap(_graph);
1530
        _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
1531
      }
1532
      if (!_delta3) {
1533
        _delta3_index = new IntEdgeMap(_graph);
1534
        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
1535
      }
1536
    }
1537

	
1538
    void destroyStructures() {
1539
      if (_matching) {
1540
        delete _matching;
1541
      }
1542
      if (_node_potential) {
1543
        delete _node_potential;
1544
      }
1545
      if (_status) {
1546
        delete _status;
1547
      }
1548
      if (_pred) {
1549
        delete _pred;
1550
      }
1551
      if (_tree_set) {
1552
        delete _tree_set_index;
1553
        delete _tree_set;
1554
      }
1555
      if (_delta2) {
1556
        delete _delta2_index;
1557
        delete _delta2;
1558
      }
1559
      if (_delta3) {
1560
        delete _delta3_index;
1561
        delete _delta3;
1562
      }
1563
    }
1564

	
1565
    void matchedToEven(Node node, int tree) {
1566
      _tree_set->insert(node, tree);
1567
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1568

	
1569
      if (_delta2->state(node) == _delta2->IN_HEAP) {
1570
        _delta2->erase(node);
1571
      }
1572

	
1573
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1574
        Node v = _graph.source(a);
1575
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1576
          dualScale * _weight[a];
1577
        if (node == v) {
1578
          if (_allow_loops && _graph.direction(a)) {
1579
            _delta3->push(a, rw / 2);
1580
          }
1581
        } else if ((*_status)[v] == EVEN) {
1582
          _delta3->push(a, rw / 2);
1583
        } else if ((*_status)[v] == MATCHED) {
1584
          if (_delta2->state(v) != _delta2->IN_HEAP) {
1585
            _pred->set(v, a);
1586
            _delta2->push(v, rw);
1587
          } else if ((*_delta2)[v] > rw) {
1588
            _pred->set(v, a);
1589
            _delta2->decrease(v, rw);
1590
          }
1591
        }
1592
      }
1593
    }
1594

	
1595
    void matchedToOdd(Node node, int tree) {
1596
      _tree_set->insert(node, tree);
1597
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1598

	
1599
      if (_delta2->state(node) == _delta2->IN_HEAP) {
1600
        _delta2->erase(node);
1601
      }
1602
    }
1603

	
1604
    void evenToMatched(Node node, int tree) {
1605
      _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
1606
      Arc min = INVALID;
1607
      Value minrw = std::numeric_limits<Value>::max();
1608
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1609
        Node v = _graph.source(a);
1610
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1611
          dualScale * _weight[a];
1612

	
1613
        if (node == v) {
1614
          if (_allow_loops && _graph.direction(a)) {
1615
            _delta3->erase(a);
1616
          }
1617
        } else if ((*_status)[v] == EVEN) {
1618
          _delta3->erase(a);
1619
          if (minrw > rw) {
1620
            min = _graph.oppositeArc(a);
1621
            minrw = rw;
1622
          }
1623
        } else if ((*_status)[v]  == MATCHED) {
1624
          if ((*_pred)[v] == a) {
1625
            Arc mina = INVALID;
1626
            Value minrwa = std::numeric_limits<Value>::max();
1627
            for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
1628
              Node va = _graph.target(aa);
1629
              if ((*_status)[va] != EVEN ||
1630
                  _tree_set->find(va) == tree) continue;
1631
              Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
1632
                dualScale * _weight[aa];
1633
              if (minrwa > rwa) {
1634
                minrwa = rwa;
1635
                mina = aa;
1636
              }
1637
            }
1638
            if (mina != INVALID) {
1639
              _pred->set(v, mina);
1640
              _delta2->increase(v, minrwa);
1641
            } else {
1642
              _pred->set(v, INVALID);
1643
              _delta2->erase(v);
1644
            }
1645
          }
1646
        }
1647
      }
1648
      if (min != INVALID) {
1649
        _pred->set(node, min);
1650
        _delta2->push(node, minrw);
1651
      } else {
1652
        _pred->set(node, INVALID);
1653
      }
1654
    }
1655

	
1656
    void oddToMatched(Node node) {
1657
      _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
1658
      Arc min = INVALID;
1659
      Value minrw = std::numeric_limits<Value>::max();
1660
      for (InArcIt a(_graph, node); a != INVALID; ++a) {
1661
        Node v = _graph.source(a);
1662
        if ((*_status)[v] != EVEN) continue;
1663
        Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
1664
          dualScale * _weight[a];
1665

	
1666
        if (minrw > rw) {
1667
          min = _graph.oppositeArc(a);
1668
          minrw = rw;
1669
        }
1670
      }
1671
      if (min != INVALID) {
1672
        _pred->set(node, min);
1673
        _delta2->push(node, minrw);
1674
      } else {
1675
        _pred->set(node, INVALID);
1676
      }
1677
    }
1678

	
1679
    void alternatePath(Node even, int tree) {
1680
      Node odd;
1681

	
1682
      _status->set(even, MATCHED);
1683
      evenToMatched(even, tree);
1684

	
1685
      Arc prev = (*_matching)[even];
1686
      while (prev != INVALID) {
1687
        odd = _graph.target(prev);
1688
        even = _graph.target((*_pred)[odd]);
1689
        _matching->set(odd, (*_pred)[odd]);
1690
        _status->set(odd, MATCHED);
1691
        oddToMatched(odd);
1692

	
1693
        prev = (*_matching)[even];
1694
        _status->set(even, MATCHED);
1695
        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
1696
        evenToMatched(even, tree);
1697
      }
1698
    }
1699

	
1700
    void destroyTree(int tree) {
1701
      for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
1702
        if ((*_status)[n] == EVEN) {
1703
          _status->set(n, MATCHED);
1704
          evenToMatched(n, tree);
1705
        } else if ((*_status)[n] == ODD) {
1706
          _status->set(n, MATCHED);
1707
          oddToMatched(n);
1708
        }
1709
      }
1710
      _tree_set->eraseClass(tree);
1711
    }
1712

	
1713
    void augmentOnEdge(const Edge& edge) {
1714
      Node left = _graph.u(edge);
1715
      int left_tree = _tree_set->find(left);
1716

	
1717
      alternatePath(left, left_tree);
1718
      destroyTree(left_tree);
1719
      _matching->set(left, _graph.direct(edge, true));
1720

	
1721
      Node right = _graph.v(edge);
1722
      int right_tree = _tree_set->find(right);
1723

	
1724
      alternatePath(right, right_tree);
1725
      destroyTree(right_tree);
1726
      _matching->set(right, _graph.direct(edge, false));
1727
    }
1728

	
1729
    void augmentOnArc(const Arc& arc) {
1730
      Node left = _graph.source(arc);
1731
      _status->set(left, MATCHED);
1732
      _matching->set(left, arc);
1733
      _pred->set(left, arc);
1734

	
1735
      Node right = _graph.target(arc);
1736
      int right_tree = _tree_set->find(right);
1737

	
1738
      alternatePath(right, right_tree);
1739
      destroyTree(right_tree);
1740
      _matching->set(right, _graph.oppositeArc(arc));
1741
    }
1742

	
1743
    void extendOnArc(const Arc& arc) {
1744
      Node base = _graph.target(arc);
1745
      int tree = _tree_set->find(base);
1746

	
1747
      Node odd = _graph.source(arc);
1748
      _tree_set->insert(odd, tree);
1749
      _status->set(odd, ODD);
1750
      matchedToOdd(odd, tree);
1751
      _pred->set(odd, arc);
1752

	
1753
      Node even = _graph.target((*_matching)[odd]);
1754
      _tree_set->insert(even, tree);
1755
      _status->set(even, EVEN);
1756
      matchedToEven(even, tree);
1757
    }
1758

	
1759
    void cycleOnEdge(const Edge& edge, int tree) {
1760
      Node nca = INVALID;
1761
      std::vector<Node> left_path, right_path;
1762

	
1763
      {
1764
        std::set<Node> left_set, right_set;
1765
        Node left = _graph.u(edge);
1766
        left_path.push_back(left);
1767
        left_set.insert(left);
1768

	
1769
        Node right = _graph.v(edge);
1770
        right_path.push_back(right);
1771
        right_set.insert(right);
1772

	
1773
        while (true) {
1774

	
1775
          if (left_set.find(right) != left_set.end()) {
1776
            nca = right;
1777
            break;
1778
          }
1779

	
1780
          if ((*_matching)[left] == INVALID) break;
1781

	
1782
          left = _graph.target((*_matching)[left]);
1783
          left_path.push_back(left);
1784
          left = _graph.target((*_pred)[left]);
1785
          left_path.push_back(left);
1786

	
1787
          left_set.insert(left);
1788

	
1789
          if (right_set.find(left) != right_set.end()) {
1790
            nca = left;
1791
            break;
1792
          }
1793

	
1794
          if ((*_matching)[right] == INVALID) break;
1795

	
1796
          right = _graph.target((*_matching)[right]);
1797
          right_path.push_back(right);
1798
          right = _graph.target((*_pred)[right]);
1799
          right_path.push_back(right);
1800

	
1801
          right_set.insert(right);
1802

	
1803
        }
1804

	
1805
        if (nca == INVALID) {
1806
          if ((*_matching)[left] == INVALID) {
1807
            nca = right;
1808
            while (left_set.find(nca) == left_set.end()) {
1809
              nca = _graph.target((*_matching)[nca]);
1810
              right_path.push_back(nca);
1811
              nca = _graph.target((*_pred)[nca]);
1812
              right_path.push_back(nca);
1813
            }
1814
          } else {
1815
            nca = left;
1816
            while (right_set.find(nca) == right_set.end()) {
1817
              nca = _graph.target((*_matching)[nca]);
1818
              left_path.push_back(nca);
1819
              nca = _graph.target((*_pred)[nca]);
1820
              left_path.push_back(nca);
1821
            }
1822
          }
1823
        }
1824
      }
1825

	
1826
      alternatePath(nca, tree);
1827
      Arc prev;
1828

	
1829
      prev = _graph.direct(edge, true);
1830
      for (int i = 0; left_path[i] != nca; i += 2) {
1831
        _matching->set(left_path[i], prev);
1832
        _status->set(left_path[i], MATCHED);
1833
        evenToMatched(left_path[i], tree);
1834

	
1835
        prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
1836
        _status->set(left_path[i + 1], MATCHED);
1837
        oddToMatched(left_path[i + 1]);
1838
      }
1839
      _matching->set(nca, prev);
1840

	
1841
      for (int i = 0; right_path[i] != nca; i += 2) {
1842
        _status->set(right_path[i], MATCHED);
1843
        evenToMatched(right_path[i], tree);
1844

	
1845
        _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
1846
        _status->set(right_path[i + 1], MATCHED);
1847
        oddToMatched(right_path[i + 1]);
1848
      }
1849

	
1850
      destroyTree(tree);
1851
    }
1852

	
1853
    void extractCycle(const Arc &arc) {
1854
      Node left = _graph.source(arc);
1855
      Node odd = _graph.target((*_matching)[left]);
1856
      Arc prev;
1857
      while (odd != left) {
1858
        Node even = _graph.target((*_matching)[odd]);
1859
        prev = (*_matching)[odd];
1860
        odd = _graph.target((*_matching)[even]);
1861
        _matching->set(even, _graph.oppositeArc(prev));
1862
      }
1863
      _matching->set(left, arc);
1864

	
1865
      Node right = _graph.target(arc);
1866
      int right_tree = _tree_set->find(right);
1867
      alternatePath(right, right_tree);
1868
      destroyTree(right_tree);
1869
      _matching->set(right, _graph.oppositeArc(arc));
1870
    }
1871

	
1872
  public:
1873

	
1874
    /// \brief Constructor
1875
    ///
1876
    /// Constructor.
1877
    MaxWeightedPerfectFractionalMatching(const Graph& graph,
1878
                                         const WeightMap& weight,
1879
                                         bool allow_loops = true)
1880
      : _graph(graph), _weight(weight), _matching(0),
1881
      _node_potential(0), _node_num(0), _allow_loops(allow_loops),
1882
      _status(0),  _pred(0),
1883
      _tree_set_index(0), _tree_set(0),
1884

	
1885
      _delta2_index(0), _delta2(0),
1886
      _delta3_index(0), _delta3(0),
1887

	
1888
      _delta_sum() {}
1889

	
1890
    ~MaxWeightedPerfectFractionalMatching() {
1891
      destroyStructures();
1892
    }
1893

	
1894
    /// \name Execution Control
1895
    /// The simplest way to execute the algorithm is to use the
1896
    /// \ref run() member function.
1897

	
1898
    ///@{
1899

	
1900
    /// \brief Initialize the algorithm
1901
    ///
1902
    /// This function initializes the algorithm.
1903
    void init() {
1904
      createStructures();
1905

	
1906
      for (NodeIt n(_graph); n != INVALID; ++n) {
1907
        (*_delta2_index)[n] = _delta2->PRE_HEAP;
1908
      }
1909
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1910
        (*_delta3_index)[e] = _delta3->PRE_HEAP;
1911
      }
1912

	
1913
      for (NodeIt n(_graph); n != INVALID; ++n) {
1914
        Value max = - std::numeric_limits<Value>::max();
1915
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1916
          if (_graph.target(e) == n && !_allow_loops) continue;
1917
          if ((dualScale * _weight[e]) / 2 > max) {
1918
            max = (dualScale * _weight[e]) / 2;
1919
          }
1920
        }
1921
        _node_potential->set(n, max);
1922

	
1923
        _tree_set->insert(n);
1924

	
1925
        _matching->set(n, INVALID);
1926
        _status->set(n, EVEN);
1927
      }
1928

	
1929
      for (EdgeIt e(_graph); e != INVALID; ++e) {
1930
        Node left = _graph.u(e);
1931
        Node right = _graph.v(e);
1932
        if (left == right && !_allow_loops) continue;
1933
        _delta3->push(e, ((*_node_potential)[left] +
1934
                          (*_node_potential)[right] -
1935
                          dualScale * _weight[e]) / 2);
1936
      }
1937
    }
1938

	
1939
    /// \brief Start the algorithm
1940
    ///
1941
    /// This function starts the algorithm.
1942
    ///
1943
    /// \pre \ref init() must be called before using this function.
1944
    bool start() {
1945
      enum OpType {
1946
        D2, D3
1947
      };
1948

	
1949
      int unmatched = _node_num;
1950
      while (unmatched > 0) {
1951
        Value d2 = !_delta2->empty() ?
1952
          _delta2->prio() : std::numeric_limits<Value>::max();
1953

	
1954
        Value d3 = !_delta3->empty() ?
1955
          _delta3->prio() : std::numeric_limits<Value>::max();
1956

	
1957
        _delta_sum = d3; OpType ot = D3;
1958
        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
1959

	
1960
        if (_delta_sum == std::numeric_limits<Value>::max()) {
1961
          return false;
1962
        }
1963

	
1964
        switch (ot) {
1965
        case D2:
1966
          {
1967
            Node n = _delta2->top();
1968
            Arc a = (*_pred)[n];
1969
            if ((*_matching)[n] == INVALID) {
1970
              augmentOnArc(a);
1971
              --unmatched;
1972
            } else {
1973
              Node v = _graph.target((*_matching)[n]);
1974
              if ((*_matching)[n] !=
1975
                  _graph.oppositeArc((*_matching)[v])) {
1976
                extractCycle(a);
1977
                --unmatched;
1978
              } else {
1979
                extendOnArc(a);
1980
              }
1981
            }
1982
          } break;
1983
        case D3:
1984
          {
1985
            Edge e = _delta3->top();
1986

	
1987
            Node left = _graph.u(e);
1988
            Node right = _graph.v(e);
1989

	
1990
            int left_tree = _tree_set->find(left);
1991
            int right_tree = _tree_set->find(right);
1992

	
1993
            if (left_tree == right_tree) {
1994
              cycleOnEdge(e, left_tree);
1995
              --unmatched;
1996
            } else {
1997
              augmentOnEdge(e);
1998
              unmatched -= 2;
1999
            }
2000
          } break;
2001
        }
2002
      }
2003
      return true;
2004
    }
2005

	
2006
    /// \brief Run the algorithm.
2007
    ///
2008
    /// This method runs the \c %MaxWeightedMatching algorithm.
2009
    ///
2010
    /// \note mwfm.run() is just a shortcut of the following code.
2011
    /// \code
2012
    ///   mwpfm.init();
2013
    ///   mwpfm.start();
2014
    /// \endcode
2015
    bool run() {
2016
      init();
2017
      return start();
2018
    }
2019

	
2020
    /// @}
2021

	
2022
    /// \name Primal Solution
2023
    /// Functions to get the primal solution, i.e. the maximum weighted
2024
    /// matching.\n
2025
    /// Either \ref run() or \ref start() function should be called before
2026
    /// using them.
2027

	
2028
    /// @{
2029

	
2030
    /// \brief Return the weight of the matching.
2031
    ///
2032
    /// This function returns the weight of the found matching. This
2033
    /// value is scaled by \ref primalScale "primal scale".
2034
    ///
2035
    /// \pre Either run() or start() must be called before using this function.
2036
    Value matchingWeight() const {
2037
      Value sum = 0;
2038
      for (NodeIt n(_graph); n != INVALID; ++n) {
2039
        if ((*_matching)[n] != INVALID) {
2040
          sum += _weight[(*_matching)[n]];
2041
        }
2042
      }
2043
      return sum * primalScale / 2;
2044
    }
2045

	
2046
    /// \brief Return the number of covered nodes in the matching.
2047
    ///
2048
    /// This function returns the number of covered nodes in the matching.
2049
    ///
2050
    /// \pre Either run() or start() must be called before using this function.
2051
    int matchingSize() const {
2052
      int num = 0;
2053
      for (NodeIt n(_graph); n != INVALID; ++n) {
2054
        if ((*_matching)[n] != INVALID) {
2055
          ++num;
2056
        }
2057
      }
2058
      return num;
2059
    }
2060

	
2061
    /// \brief Return \c true if the given edge is in the matching.
2062
    ///
2063
    /// This function returns \c true if the given edge is in the
2064
    /// found matching. The result is scaled by \ref primalScale
2065
    /// "primal scale".
2066
    ///
2067
    /// \pre Either run() or start() must be called before using this function.
2068
    Value matching(const Edge& edge) const {
2069
      return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
2070
        * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
2071
        * primalScale / 2;
2072
    }
2073

	
2074
    /// \brief Return the fractional matching arc (or edge) incident
2075
    /// to the given node.
2076
    ///
2077
    /// This function returns one of the fractional matching arc (or
2078
    /// edge) incident to the given node in the found matching or \c
2079
    /// INVALID if the node is not covered by the matching or if the
2080
    /// node is on an odd length cycle then it is the successor edge
2081
    /// on the cycle.
2082
    ///
2083
    /// \pre Either run() or start() must be called before using this function.
2084
    Arc matching(const Node& node) const {
2085
      return (*_matching)[node];
2086
    }
2087

	
2088
    /// \brief Return a const reference to the matching map.
2089
    ///
2090
    /// This function returns a const reference to a node map that stores
2091
    /// the matching arc (or edge) incident to each node.
2092
    const MatchingMap& matchingMap() const {
2093
      return *_matching;
2094
    }
2095

	
2096
    /// @}
2097

	
2098
    /// \name Dual Solution
2099
    /// Functions to get the dual solution.\n
2100
    /// Either \ref run() or \ref start() function should be called before
2101
    /// using them.
2102

	
2103
    /// @{
2104

	
2105
    /// \brief Return the value of the dual solution.
2106
    ///
2107
    /// This function returns the value of the dual solution.
2108
    /// It should be equal to the primal value scaled by \ref dualScale
2109
    /// "dual scale".
2110
    ///
2111
    /// \pre Either run() or start() must be called before using this function.
2112
    Value dualValue() const {
2113
      Value sum = 0;
2114
      for (NodeIt n(_graph); n != INVALID; ++n) {
2115
        sum += nodeValue(n);
2116
      }
2117
      return sum;
2118
    }
2119

	
2120
    /// \brief Return the dual value (potential) of the given node.
2121
    ///
2122
    /// This function returns the dual value (potential) of the given node.
2123
    ///
2124
    /// \pre Either run() or start() must be called before using this function.
2125
    Value nodeValue(const Node& n) const {
2126
      return (*_node_potential)[n];
2127
    }
2128

	
2129
    /// @}
2130

	
2131
  };
2132

	
2133
} //END OF NAMESPACE LEMON
2134

	
2135
#endif //LEMON_FRACTIONAL_MATCHING_H
Ignore white space 6 line context
1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2
 *
3
 * This file is a part of LEMON, a generic C++ optimization library.
4
 *
5
 * Copyright (C) 2003-2009
6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8
 *
9
 * Permission to use, modify and distribute this software is granted
10
 * provided that this copyright notice appears in all copies. For
11
 * precise terms see the accompanying LICENSE file.
12
 *
13
 * This software is provided "AS IS" with no warranty of any kind,
14
 * express or implied, and with no claim as to its suitability for any
15
 * purpose.
16
 *
17
 */
18

	
19
#include <iostream>
20
#include <sstream>
21
#include <vector>
22
#include <queue>
23
#include <cstdlib>
24

	
25
#include <lemon/fractional_matching.h>
26
#include <lemon/smart_graph.h>
27
#include <lemon/concepts/graph.h>
28
#include <lemon/concepts/maps.h>
29
#include <lemon/lgf_reader.h>
30
#include <lemon/math.h>
31

	
32
#include "test_tools.h"
33

	
34
using namespace std;
35
using namespace lemon;
36

	
37
GRAPH_TYPEDEFS(SmartGraph);
38

	
39

	
40
const int lgfn = 4;
41
const std::string lgf[lgfn] = {
42
  "@nodes\n"
43
  "label\n"
44
  "0\n"
45
  "1\n"
46
  "2\n"
47
  "3\n"
48
  "4\n"
49
  "5\n"
50
  "6\n"
51
  "7\n"
52
  "@edges\n"
53
  "     label  weight\n"
54
  "7 4  0      984\n"
55
  "0 7  1      73\n"
56
  "7 1  2      204\n"
57
  "2 3  3      583\n"
58
  "2 7  4      565\n"
59
  "2 1  5      582\n"
60
  "0 4  6      551\n"
61
  "2 5  7      385\n"
62
  "1 5  8      561\n"
63
  "5 3  9      484\n"
64
  "7 5  10     904\n"
65
  "3 6  11     47\n"
66
  "7 6  12     888\n"
67
  "3 0  13     747\n"
68
  "6 1  14     310\n",
69

	
70
  "@nodes\n"
71
  "label\n"
72
  "0\n"
73
  "1\n"
74
  "2\n"
75
  "3\n"
76
  "4\n"
77
  "5\n"
78
  "6\n"
79
  "7\n"
80
  "@edges\n"
81
  "     label  weight\n"
82
  "2 5  0      710\n"
83
  "0 5  1      241\n"
84
  "2 4  2      856\n"
85
  "2 6  3      762\n"
86
  "4 1  4      747\n"
87
  "6 1  5      962\n"
88
  "4 7  6      723\n"
89
  "1 7  7      661\n"
90
  "2 3  8      376\n"
91
  "1 0  9      416\n"
92
  "6 7  10     391\n",
93

	
94
  "@nodes\n"
95
  "label\n"
96
  "0\n"
97
  "1\n"
98
  "2\n"
99
  "3\n"
100
  "4\n"
101
  "5\n"
102
  "6\n"
103
  "7\n"
104
  "@edges\n"
105
  "     label  weight\n"
106
  "6 2  0      553\n"
107
  "0 7  1      653\n"
108
  "6 3  2      22\n"
109
  "4 7  3      846\n"
110
  "7 2  4      981\n"
111
  "7 6  5      250\n"
112
  "5 2  6      539\n",
113

	
114
  "@nodes\n"
115
  "label\n"
116
  "0\n"
117
  "@edges\n"
118
  "     label  weight\n"
119
  "0 0  0      100\n"
120
};
121

	
122
void checkMaxFractionalMatchingCompile()
123
{
124
  typedef concepts::Graph Graph;
125
  typedef Graph::Node Node;
126
  typedef Graph::Edge Edge;
127

	
128
  Graph g;
129
  Node n;
130
  Edge e;
131

	
132
  MaxFractionalMatching<Graph> mat_test(g);
133
  const MaxFractionalMatching<Graph>&
134
    const_mat_test = mat_test;
135

	
136
  mat_test.init();
137
  mat_test.start();
138
  mat_test.start(true);
139
  mat_test.startPerfect();
140
  mat_test.startPerfect(true);
141
  mat_test.run();
142
  mat_test.run(true);
143
  mat_test.runPerfect();
144
  mat_test.runPerfect(true);
145

	
146
  const_mat_test.matchingSize();
147
  const_mat_test.matching(e);
148
  const_mat_test.matching(n);
149
  const MaxFractionalMatching<Graph>::MatchingMap& mmap =
150
    const_mat_test.matchingMap();
151
  e = mmap[n];
152

	
153
  const_mat_test.barrier(n);
154
}
155

	
156
void checkMaxWeightedFractionalMatchingCompile()
157
{
158
  typedef concepts::Graph Graph;
159
  typedef Graph::Node Node;
160
  typedef Graph::Edge Edge;
161
  typedef Graph::EdgeMap<int> WeightMap;
162

	
163
  Graph g;
164
  Node n;
165
  Edge e;
166
  WeightMap w(g);
167

	
168
  MaxWeightedFractionalMatching<Graph> mat_test(g, w);
169
  const MaxWeightedFractionalMatching<Graph>&
170
    const_mat_test = mat_test;
171

	
172
  mat_test.init();
173
  mat_test.start();
174
  mat_test.run();
175

	
176
  const_mat_test.matchingWeight();
177
  const_mat_test.matchingSize();
178
  const_mat_test.matching(e);
179
  const_mat_test.matching(n);
180
  const MaxWeightedFractionalMatching<Graph>::MatchingMap& mmap =
181
    const_mat_test.matchingMap();
182
  e = mmap[n];
183

	
184
  const_mat_test.dualValue();
185
  const_mat_test.nodeValue(n);
186
}
187

	
188
void checkMaxWeightedPerfectFractionalMatchingCompile()
189
{
190
  typedef concepts::Graph Graph;
191
  typedef Graph::Node Node;
192
  typedef Graph::Edge Edge;
193
  typedef Graph::EdgeMap<int> WeightMap;
194

	
195
  Graph g;
196
  Node n;
197
  Edge e;
198
  WeightMap w(g);
199

	
200
  MaxWeightedPerfectFractionalMatching<Graph> mat_test(g, w);
201
  const MaxWeightedPerfectFractionalMatching<Graph>&
202
    const_mat_test = mat_test;
203

	
204
  mat_test.init();
205
  mat_test.start();
206
  mat_test.run();
207

	
208
  const_mat_test.matchingWeight();
209
  const_mat_test.matching(e);
210
  const_mat_test.matching(n);
211
  const MaxWeightedPerfectFractionalMatching<Graph>::MatchingMap& mmap =
212
    const_mat_test.matchingMap();
213
  e = mmap[n];
214

	
215
  const_mat_test.dualValue();
216
  const_mat_test.nodeValue(n);
217
}
218

	
219
void checkFractionalMatching(const SmartGraph& graph,
220
                             const MaxFractionalMatching<SmartGraph>& mfm,
221
                             bool allow_loops = true) {
222
  int pv = 0;
223
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
224
    int indeg = 0;
225
    for (InArcIt a(graph, n); a != INVALID; ++a) {
226
      if (mfm.matching(graph.source(a)) == a) {
227
        ++indeg;
228
      }
229
    }
230
    if (mfm.matching(n) != INVALID) {
231
      check(indeg == 1, "Invalid matching");
232
      ++pv;
233
    } else {
234
      check(indeg == 0, "Invalid matching");
235
    }
236
  }
237
  check(pv == mfm.matchingSize(), "Wrong matching size");
238

	
239
  SmartGraph::NodeMap<bool> processed(graph, false);
240
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
241
    if (processed[n]) continue;
242
    processed[n] = true;
243
    if (mfm.matching(n) == INVALID) continue;
244
    int num = 1;
245
    Node v = graph.target(mfm.matching(n));
246
    while (v != n) {
247
      processed[v] = true;
248
      ++num;
249
      v = graph.target(mfm.matching(v));
250
    }
251
    check(num == 2 || num % 2 == 1, "Wrong cycle size");
252
    check(allow_loops || num != 1, "Wrong cycle size");
253
  }
254

	
255
  int anum = 0, bnum = 0;
256
  SmartGraph::NodeMap<bool> neighbours(graph, false);
257
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
258
    if (!mfm.barrier(n)) continue;
259
    ++anum;
260
    for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
261
      Node u = graph.source(a);
262
      if (!allow_loops && u == n) continue;
263
      if (!neighbours[u]) {
264
        neighbours[u] = true;
265
        ++bnum;
266
      }
267
    }
268
  }
269
  check(anum - bnum + mfm.matchingSize() == countNodes(graph),
270
        "Wrong barrier");
271
}
272

	
273
void checkPerfectFractionalMatching(const SmartGraph& graph,
274
                             const MaxFractionalMatching<SmartGraph>& mfm,
275
                             bool perfect, bool allow_loops = true) {
276
  if (perfect) {
277
    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
278
      int indeg = 0;
279
      for (InArcIt a(graph, n); a != INVALID; ++a) {
280
        if (mfm.matching(graph.source(a)) == a) {
281
          ++indeg;
282
        }
283
      }
284
      check(mfm.matching(n) != INVALID, "Invalid matching");
285
      check(indeg == 1, "Invalid matching");
286
    }
287
  } else {
288
    int anum = 0, bnum = 0;
289
    SmartGraph::NodeMap<bool> neighbours(graph, false);
290
    for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
291
      if (!mfm.barrier(n)) continue;
292
      ++anum;
293
      for (SmartGraph::InArcIt a(graph, n); a != INVALID; ++a) {
294
        Node u = graph.source(a);
295
        if (!allow_loops && u == n) continue;
296
        if (!neighbours[u]) {
297
          neighbours[u] = true;
298
          ++bnum;
299
        }
300
      }
301
    }
302
    check(anum - bnum > 0, "Wrong barrier");
303
  }
304
}
305

	
306
void checkWeightedFractionalMatching(const SmartGraph& graph,
307
                   const SmartGraph::EdgeMap<int>& weight,
308
                   const MaxWeightedFractionalMatching<SmartGraph>& mwfm,
309
                   bool allow_loops = true) {
310
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
311
    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
312
    int rw = mwfm.nodeValue(graph.u(e)) + mwfm.nodeValue(graph.v(e))
313
      - weight[e] * mwfm.dualScale;
314

	
315
    check(rw >= 0, "Negative reduced weight");
316
    check(rw == 0 || !mwfm.matching(e),
317
          "Non-zero reduced weight on matching edge");
318
  }
319

	
320
  int pv = 0;
321
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
322
    int indeg = 0;
323
    for (InArcIt a(graph, n); a != INVALID; ++a) {
324
      if (mwfm.matching(graph.source(a)) == a) {
325
        ++indeg;
326
      }
327
    }
328
    check(indeg <= 1, "Invalid matching");
329
    if (mwfm.matching(n) != INVALID) {
330
      check(mwfm.nodeValue(n) >= 0, "Invalid node value");
331
      check(indeg == 1, "Invalid matching");
332
      pv += weight[mwfm.matching(n)];
333
      SmartGraph::Node o = graph.target(mwfm.matching(n));
334
    } else {
335
      check(mwfm.nodeValue(n) == 0, "Invalid matching");
336
      check(indeg == 0, "Invalid matching");
337
    }
338
  }
339

	
340
  int dv = 0;
341
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
342
    dv += mwfm.nodeValue(n);
343
  }
344

	
345
  check(pv * mwfm.dualScale == dv * 2, "Wrong duality");
346

	
347
  SmartGraph::NodeMap<bool> processed(graph, false);
348
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
349
    if (processed[n]) continue;
350
    processed[n] = true;
351
    if (mwfm.matching(n) == INVALID) continue;
352
    int num = 1;
353
    Node v = graph.target(mwfm.matching(n));
354
    while (v != n) {
355
      processed[v] = true;
356
      ++num;
357
      v = graph.target(mwfm.matching(v));
358
    }
359
    check(num == 2 || num % 2 == 1, "Wrong cycle size");
360
    check(allow_loops || num != 1, "Wrong cycle size");
361
  }
362

	
363
  return;
364
}
365

	
366
void checkWeightedPerfectFractionalMatching(const SmartGraph& graph,
367
                const SmartGraph::EdgeMap<int>& weight,
368
                const MaxWeightedPerfectFractionalMatching<SmartGraph>& mwpfm,
369
                bool allow_loops = true) {
370
  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
371
    if (graph.u(e) == graph.v(e) && !allow_loops) continue;
372
    int rw = mwpfm.nodeValue(graph.u(e)) + mwpfm.nodeValue(graph.v(e))
373
      - weight[e] * mwpfm.dualScale;
374

	
375
    check(rw >= 0, "Negative reduced weight");
376
    check(rw == 0 || !mwpfm.matching(e),
377
          "Non-zero reduced weight on matching edge");
378
  }
379

	
380
  int pv = 0;
381
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
382
    int indeg = 0;
383
    for (InArcIt a(graph, n); a != INVALID; ++a) {
384
      if (mwpfm.matching(graph.source(a)) == a) {
385
        ++indeg;
386
      }
387
    }
388
    check(mwpfm.matching(n) != INVALID, "Invalid perfect matching");
389
    check(indeg == 1, "Invalid perfect matching");
390
    pv += weight[mwpfm.matching(n)];
391
    SmartGraph::Node o = graph.target(mwpfm.matching(n));
392
  }
393

	
394
  int dv = 0;
395
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
396
    dv += mwpfm.nodeValue(n);
397
  }
398

	
399
  check(pv * mwpfm.dualScale == dv * 2, "Wrong duality");
400

	
401
  SmartGraph::NodeMap<bool> processed(graph, false);
402
  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
403
    if (processed[n]) continue;
404
    processed[n] = true;
405
    if (mwpfm.matching(n) == INVALID) continue;
406
    int num = 1;
407
    Node v = graph.target(mwpfm.matching(n));
408
    while (v != n) {
409
      processed[v] = true;
410
      ++num;
411
      v = graph.target(mwpfm.matching(v));
412
    }
413
    check(num == 2 || num % 2 == 1, "Wrong cycle size");
414
    check(allow_loops || num != 1, "Wrong cycle size");
415
  }
416

	
417
  return;
418
}
419

	
420

	
421
int main() {
422

	
423
  for (int i = 0; i < lgfn; ++i) {
424
    SmartGraph graph;
425
    SmartGraph::EdgeMap<int> weight(graph);
426

	
427
    istringstream lgfs(lgf[i]);
428
    graphReader(graph, lgfs).
429
      edgeMap("weight", weight).run();
430

	
431
    bool perfect_with_loops;
432
    {
433
      MaxFractionalMatching<SmartGraph> mfm(graph, true);
434
      mfm.run();
435
      checkFractionalMatching(graph, mfm, true);
436
      perfect_with_loops = mfm.matchingSize() == countNodes(graph);
437
    }
438

	
439
    bool perfect_without_loops;
440
    {
441
      MaxFractionalMatching<SmartGraph> mfm(graph, false);
442
      mfm.run();
443
      checkFractionalMatching(graph, mfm, false);
444
      perfect_without_loops = mfm.matchingSize() == countNodes(graph);
445
    }
446

	
447
    {
448
      MaxFractionalMatching<SmartGraph> mfm(graph, true);
449
      bool result = mfm.runPerfect();
450
      checkPerfectFractionalMatching(graph, mfm, result, true);
451
      check(result == perfect_with_loops, "Wrong perfect matching");
452
    }
453

	
454
    {
455
      MaxFractionalMatching<SmartGraph> mfm(graph, false);
456
      bool result = mfm.runPerfect();
457
      checkPerfectFractionalMatching(graph, mfm, result, false);
458
      check(result == perfect_without_loops, "Wrong perfect matching");
459
    }
460

	
461
    {
462
      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, true);
463
      mwfm.run();
464
      checkWeightedFractionalMatching(graph, weight, mwfm, true);
465
    }
466

	
467
    {
468
      MaxWeightedFractionalMatching<SmartGraph> mwfm(graph, weight, false);
469
      mwfm.run();
470
      checkWeightedFractionalMatching(graph, weight, mwfm, false);
471
    }
472

	
473
    {
474
      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
475
                                                             true);
476
      bool perfect = mwpfm.run();
477
      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
478
            "Perfect matching found");
479
      check(perfect == perfect_with_loops, "Wrong perfect matching");
480

	
481
      if (perfect) {
482
        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, true);
483
      }
484
    }
485

	
486
    {
487
      MaxWeightedPerfectFractionalMatching<SmartGraph> mwpfm(graph, weight,
488
                                                             false);
489
      bool perfect = mwpfm.run();
490
      check(perfect == (mwpfm.matchingSize() == countNodes(graph)),
491
            "Perfect matching found");
492
      check(perfect == perfect_without_loops, "Wrong perfect matching");
493

	
494
      if (perfect) {
495
        checkWeightedPerfectFractionalMatching(graph, weight, mwpfm, false);
496
      }
497
    }
498

	
499
  }
500

	
501
  return 0;
502
}
Ignore white space 6 line context
... ...
@@ -351,3 +351,3 @@
351 351

	
352
\ref Circulation is a preflow push-relabel algorithm implemented directly 
352
\ref Circulation is a preflow push-relabel algorithm implemented directly
353 353
for finding feasible circulations, which is a somewhat different problem,
... ...
@@ -472,2 +472,9 @@
472 472
  perfect matching in general graphs.
473
- \ref MaxFractionalMatching Push-relabel algorithm for calculating
474
  maximum cardinality fractional matching in general graphs.
475
- \ref MaxWeightedFractionalMatching Augmenting path algorithm for calculating
476
  maximum weighted fractional matching in general graphs.
477
- \ref MaxWeightedPerfectFractionalMatching
478
  Augmenting path algorithm for calculating maximum weighted
479
  perfect fractional matching in general graphs.
473 480

	
Ignore white space 6 line context
... ...
@@ -83,2 +83,3 @@
83 83
	lemon/fourary_heap.h \
84
	lemon/fractional_matching.h \
84 85
	lemon/full_graph.h \
Ignore white space 6 line context
... ...
@@ -23,2 +23,3 @@
23 23
  euler_test
24
  fractional_matching_test
24 25
  gomory_hu_test
Ignore white space 6 line context
... ...
@@ -21,2 +21,3 @@
21 21
	test/euler_test \
22
	test/fractional_matching_test \
22 23
	test/gomory_hu_test \
... ...
@@ -67,2 +68,3 @@
67 68
test_euler_test_SOURCES = test/euler_test.cc
69
test_fractional_matching_test_SOURCES = test/fractional_matching_test.cc
68 70
test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
0 comments (0 inline)