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 ↑
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
#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 128 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
... ...
@@ -288,249 +288,256 @@
288 288
implemented in LEMON.
289 289
*/
290 290

	
291 291
/**
292 292
@defgroup search Graph Search
293 293
@ingroup algs
294 294
\brief Common graph search algorithms.
295 295

	
296 296
This group contains the common graph search algorithms, namely
297 297
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
298 298
*/
299 299

	
300 300
/**
301 301
@defgroup shortest_path Shortest Path Algorithms
302 302
@ingroup algs
303 303
\brief Algorithms for finding shortest paths.
304 304

	
305 305
This group contains the algorithms for finding shortest paths in digraphs.
306 306

	
307 307
 - \ref Dijkstra algorithm for finding shortest paths from a source node
308 308
   when all arc lengths are non-negative.
309 309
 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
310 310
   from a source node when arc lenghts can be either positive or negative,
311 311
   but the digraph should not contain directed cycles with negative total
312 312
   length.
313 313
 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
314 314
   for solving the \e all-pairs \e shortest \e paths \e problem when arc
315 315
   lenghts can be either positive or negative, but the digraph should
316 316
   not contain directed cycles with negative total length.
317 317
 - \ref Suurballe A successive shortest path algorithm for finding
318 318
   arc-disjoint paths between two nodes having minimum total length.
319 319
*/
320 320

	
321 321
/**
322 322
@defgroup max_flow Maximum Flow Algorithms
323 323
@ingroup algs
324 324
\brief Algorithms for finding maximum flows.
325 325

	
326 326
This group contains the algorithms for finding maximum flows and
327 327
feasible circulations.
328 328

	
329 329
The \e maximum \e flow \e problem is to find a flow of maximum value between
330 330
a single source and a single target. Formally, there is a \f$G=(V,A)\f$
331 331
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
332 332
\f$s, t \in V\f$ source and target nodes.
333 333
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
334 334
following optimization problem.
335 335

	
336 336
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
337 337
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
338 338
    \quad \forall u\in V\setminus\{s,t\} \f]
339 339
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
340 340

	
341 341
LEMON contains several algorithms for solving maximum flow problems:
342 342
- \ref EdmondsKarp Edmonds-Karp algorithm.
343 343
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
344 344
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
345 345
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
346 346

	
347 347
In most cases the \ref Preflow "Preflow" algorithm provides the
348 348
fastest method for computing a maximum flow. All implementations
349 349
also provide functions to query the minimum cut, which is the dual
350 350
problem of maximum flow.
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,
354 354
but it is strongly related to maximum flow.
355 355
For more information, see \ref Circulation.
356 356
*/
357 357

	
358 358
/**
359 359
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
360 360
@ingroup algs
361 361

	
362 362
\brief Algorithms for finding minimum cost flows and circulations.
363 363

	
364 364
This group contains the algorithms for finding minimum cost flows and
365 365
circulations. For more information about this problem and its dual
366 366
solution see \ref min_cost_flow "Minimum Cost Flow Problem".
367 367

	
368 368
LEMON contains several algorithms for this problem.
369 369
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
370 370
   pivot strategies.
371 371
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
372 372
   cost scaling.
373 373
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
374 374
   capacity scaling.
375 375
 - \ref CancelAndTighten The Cancel and Tighten algorithm.
376 376
 - \ref CycleCanceling Cycle-Canceling algorithms.
377 377

	
378 378
In general NetworkSimplex is the most efficient implementation,
379 379
but in special cases other algorithms could be faster.
380 380
For example, if the total supply and/or capacities are rather small,
381 381
CapacityScaling is usually the fastest algorithm (without effective scaling).
382 382
*/
383 383

	
384 384
/**
385 385
@defgroup min_cut Minimum Cut Algorithms
386 386
@ingroup algs
387 387

	
388 388
\brief Algorithms for finding minimum cut in graphs.
389 389

	
390 390
This group contains the algorithms for finding minimum cut in graphs.
391 391

	
392 392
The \e minimum \e cut \e problem is to find a non-empty and non-complete
393 393
\f$X\f$ subset of the nodes with minimum overall capacity on
394 394
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
395 395
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
396 396
cut is the \f$X\f$ solution of the next optimization problem:
397 397

	
398 398
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
399 399
    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]
400 400

	
401 401
LEMON contains several algorithms related to minimum cut problems:
402 402

	
403 403
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
404 404
  in directed graphs.
405 405
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
406 406
  calculating minimum cut in undirected graphs.
407 407
- \ref GomoryHu "Gomory-Hu tree computation" for calculating
408 408
  all-pairs minimum cut in undirected graphs.
409 409

	
410 410
If you want to find minimum cut just between two distinict nodes,
411 411
see the \ref max_flow "maximum flow problem".
412 412
*/
413 413

	
414 414
/**
415 415
@defgroup graph_properties Connectivity and Other Graph Properties
416 416
@ingroup algs
417 417
\brief Algorithms for discovering the graph properties
418 418

	
419 419
This group contains the algorithms for discovering the graph properties
420 420
like connectivity, bipartiteness, euler property, simplicity etc.
421 421

	
422 422
\image html edge_biconnected_components.png
423 423
\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth
424 424
*/
425 425

	
426 426
/**
427 427
@defgroup planar Planarity Embedding and Drawing
428 428
@ingroup algs
429 429
\brief Algorithms for planarity checking, embedding and drawing
430 430

	
431 431
This group contains the algorithms for planarity checking,
432 432
embedding and drawing.
433 433

	
434 434
\image html planar.png
435 435
\image latex planar.eps "Plane graph" width=\textwidth
436 436
*/
437 437

	
438 438
/**
439 439
@defgroup matching Matching Algorithms
440 440
@ingroup algs
441 441
\brief Algorithms for finding matchings in graphs and bipartite graphs.
442 442

	
443 443
This group contains the algorithms for calculating
444 444
matchings in graphs and bipartite graphs. The general matching problem is
445 445
finding a subset of the edges for which each node has at most one incident
446 446
edge.
447 447

	
448 448
There are several different algorithms for calculate matchings in
449 449
graphs.  The matching problems in bipartite graphs are generally
450 450
easier than in general graphs. The goal of the matching optimization
451 451
can be finding maximum cardinality, maximum weight or minimum cost
452 452
matching. The search can be constrained to find perfect or
453 453
maximum cardinality matching.
454 454

	
455 455
The matching algorithms implemented in LEMON:
456 456
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
457 457
  for calculating maximum cardinality matching in bipartite graphs.
458 458
- \ref PrBipartiteMatching Push-relabel algorithm
459 459
  for calculating maximum cardinality matching in bipartite graphs.
460 460
- \ref MaxWeightedBipartiteMatching
461 461
  Successive shortest path algorithm for calculating maximum weighted
462 462
  matching and maximum weighted bipartite matching in bipartite graphs.
463 463
- \ref MinCostMaxBipartiteMatching
464 464
  Successive shortest path algorithm for calculating minimum cost maximum
465 465
  matching in bipartite graphs.
466 466
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
467 467
  maximum cardinality matching in general graphs.
468 468
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
469 469
  maximum weighted matching in general graphs.
470 470
- \ref MaxWeightedPerfectMatching
471 471
  Edmond's blossom shrinking algorithm for calculating maximum weighted
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

	
474 481
\image html bipartite_matching.png
475 482
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
476 483
*/
477 484

	
478 485
/**
479 486
@defgroup spantree Minimum Spanning Tree Algorithms
480 487
@ingroup algs
481 488
\brief Algorithms for finding minimum cost spanning trees and arborescences.
482 489

	
483 490
This group contains the algorithms for finding minimum cost spanning
484 491
trees and arborescences.
485 492
*/
486 493

	
487 494
/**
488 495
@defgroup auxalg Auxiliary Algorithms
489 496
@ingroup algs
490 497
\brief Auxiliary algorithms implemented in LEMON.
491 498

	
492 499
This group contains some algorithms implemented in LEMON
493 500
in order to make it easier to implement complex algorithms.
494 501
*/
495 502

	
496 503
/**
497 504
@defgroup approx Approximation Algorithms
498 505
@ingroup algs
499 506
\brief Approximation algorithms.
500 507

	
501 508
This group contains the approximation and heuristic algorithms
502 509
implemented in LEMON.
503 510
*/
504 511

	
505 512
/**
506 513
@defgroup gen_opt_group General Optimization Tools
507 514
\brief This group contains some general optimization frameworks
508 515
implemented in LEMON.
509 516

	
510 517
This group contains some general optimization frameworks
511 518
implemented in LEMON.
512 519
*/
513 520

	
514 521
/**
515 522
@defgroup lp_group Lp and Mip Solvers
516 523
@ingroup gen_opt_group
517 524
\brief Lp and Mip solver interfaces for LEMON.
518 525

	
519 526
This group contains Lp and Mip solver interfaces for LEMON. The
520 527
various LP solvers could be used in the same manner with this
521 528
interface.
522 529
*/
523 530

	
524 531
/**
525 532
@defgroup lp_utils Tools for Lp and Mip Solvers
526 533
@ingroup lp_group
527 534
\brief Helper tools to the Lp and Mip solvers.
528 535

	
529 536
This group adds some helper tools to general optimization framework
530 537
implemented in LEMON.
531 538
*/
532 539

	
533 540
/**
534 541
@defgroup metah Metaheuristics
535 542
@ingroup gen_opt_group
536 543
\brief Metaheuristics for LEMON library.
Ignore white space 6 line context
... ...
@@ -20,122 +20,123 @@
20 20

	
21 21
lemon_libemon_la_CXXFLAGS = \
22 22
	$(AM_CXXFLAGS) \
23 23
	$(GLPK_CFLAGS) \
24 24
	$(CPLEX_CFLAGS) \
25 25
	$(SOPLEX_CXXFLAGS) \
26 26
	$(CLP_CXXFLAGS) \
27 27
	$(CBC_CXXFLAGS)
28 28

	
29 29
lemon_libemon_la_LDFLAGS = \
30 30
	$(GLPK_LIBS) \
31 31
	$(CPLEX_LIBS) \
32 32
	$(SOPLEX_LIBS) \
33 33
	$(CLP_LIBS) \
34 34
	$(CBC_LIBS)
35 35

	
36 36
if HAVE_GLPK
37 37
lemon_libemon_la_SOURCES += lemon/glpk.cc
38 38
endif
39 39

	
40 40
if HAVE_CPLEX
41 41
lemon_libemon_la_SOURCES += lemon/cplex.cc
42 42
endif
43 43

	
44 44
if HAVE_SOPLEX
45 45
lemon_libemon_la_SOURCES += lemon/soplex.cc
46 46
endif
47 47

	
48 48
if HAVE_CLP
49 49
lemon_libemon_la_SOURCES += lemon/clp.cc
50 50
endif
51 51

	
52 52
if HAVE_CBC
53 53
lemon_libemon_la_SOURCES += lemon/cbc.cc
54 54
endif
55 55

	
56 56
lemon_HEADERS += \
57 57
	lemon/adaptors.h \
58 58
	lemon/arg_parser.h \
59 59
	lemon/assert.h \
60 60
	lemon/bellman_ford.h \
61 61
	lemon/bfs.h \
62 62
	lemon/bin_heap.h \
63 63
	lemon/binom_heap.h \
64 64
	lemon/bucket_heap.h \
65 65
	lemon/cbc.h \
66 66
	lemon/circulation.h \
67 67
	lemon/clp.h \
68 68
	lemon/color.h \
69 69
	lemon/concept_check.h \
70 70
	lemon/connectivity.h \
71 71
	lemon/counter.h \
72 72
	lemon/core.h \
73 73
	lemon/cplex.h \
74 74
	lemon/dfs.h \
75 75
	lemon/dijkstra.h \
76 76
	lemon/dim2.h \
77 77
	lemon/dimacs.h \
78 78
	lemon/edge_set.h \
79 79
	lemon/elevator.h \
80 80
	lemon/error.h \
81 81
	lemon/euler.h \
82 82
	lemon/fib_heap.h \
83 83
	lemon/fourary_heap.h \
84
	lemon/fractional_matching.h \
84 85
	lemon/full_graph.h \
85 86
	lemon/glpk.h \
86 87
	lemon/gomory_hu.h \
87 88
	lemon/graph_to_eps.h \
88 89
	lemon/grid_graph.h \
89 90
	lemon/hypercube_graph.h \
90 91
	lemon/kary_heap.h \
91 92
	lemon/kruskal.h \
92 93
	lemon/hao_orlin.h \
93 94
	lemon/lgf_reader.h \
94 95
	lemon/lgf_writer.h \
95 96
	lemon/list_graph.h \
96 97
	lemon/lp.h \
97 98
	lemon/lp_base.h \
98 99
	lemon/lp_skeleton.h \
99 100
	lemon/maps.h \
100 101
	lemon/matching.h \
101 102
	lemon/math.h \
102 103
	lemon/min_cost_arborescence.h \
103 104
	lemon/nauty_reader.h \
104 105
	lemon/network_simplex.h \
105 106
	lemon/pairing_heap.h \
106 107
	lemon/path.h \
107 108
	lemon/preflow.h \
108 109
	lemon/radix_heap.h \
109 110
	lemon/radix_sort.h \
110 111
	lemon/random.h \
111 112
	lemon/smart_graph.h \
112 113
	lemon/soplex.h \
113 114
	lemon/suurballe.h \
114 115
	lemon/time_measure.h \
115 116
	lemon/tolerance.h \
116 117
	lemon/unionfind.h \
117 118
	lemon/bits/windows.h
118 119

	
119 120
bits_HEADERS += \
120 121
	lemon/bits/alteration_notifier.h \
121 122
	lemon/bits/array_map.h \
122 123
	lemon/bits/bezier.h \
123 124
	lemon/bits/default_map.h \
124 125
	lemon/bits/edge_set_extender.h \
125 126
	lemon/bits/enable_if.h \
126 127
	lemon/bits/graph_adaptor_extender.h \
127 128
	lemon/bits/graph_extender.h \
128 129
	lemon/bits/map_extender.h \
129 130
	lemon/bits/path_dump.h \
130 131
	lemon/bits/solver_bits.h \
131 132
	lemon/bits/traits.h \
132 133
	lemon/bits/variant.h \
133 134
	lemon/bits/vector_map.h
134 135

	
135 136
concept_HEADERS += \
136 137
	lemon/concepts/digraph.h \
137 138
	lemon/concepts/graph.h \
138 139
	lemon/concepts/graph_components.h \
139 140
	lemon/concepts/heap.h \
140 141
	lemon/concepts/maps.h \
141 142
	lemon/concepts/path.h
Ignore white space 6 line context
1 1
INCLUDE_DIRECTORIES(
2 2
  ${PROJECT_SOURCE_DIR}
3 3
  ${PROJECT_BINARY_DIR}
4 4
)
5 5

	
6 6
LINK_DIRECTORIES(
7 7
  ${PROJECT_BINARY_DIR}/lemon
8 8
)
9 9

	
10 10
SET(TESTS
11 11
  adaptors_test
12 12
  bellman_ford_test
13 13
  bfs_test
14 14
  circulation_test
15 15
  connectivity_test
16 16
  counter_test
17 17
  dfs_test
18 18
  digraph_test
19 19
  dijkstra_test
20 20
  dim_test
21 21
  edge_set_test
22 22
  error_test
23 23
  euler_test
24
  fractional_matching_test
24 25
  gomory_hu_test
25 26
  graph_copy_test
26 27
  graph_test
27 28
  graph_utils_test
28 29
  hao_orlin_test
29 30
  heap_test
30 31
  kruskal_test
31 32
  maps_test
32 33
  matching_test
33 34
  min_cost_arborescence_test
34 35
  min_cost_flow_test
35 36
  path_test
36 37
  preflow_test
37 38
  radix_sort_test
38 39
  random_test
39 40
  suurballe_test
40 41
  time_measure_test
41 42
  unionfind_test
42 43
)
43 44

	
44 45
IF(LEMON_HAVE_LP)
45 46
  ADD_EXECUTABLE(lp_test lp_test.cc)
46 47
  SET(LP_TEST_LIBS lemon)
47 48

	
48 49
  IF(LEMON_HAVE_GLPK)
49 50
    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${GLPK_LIBRARIES})
50 51
  ENDIF()
51 52
  IF(LEMON_HAVE_CPLEX)
52 53
    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${CPLEX_LIBRARIES})
53 54
  ENDIF()
54 55
  IF(LEMON_HAVE_CLP)
55 56
    SET(LP_TEST_LIBS ${LP_TEST_LIBS} ${COIN_CLP_LIBRARIES})
56 57
  ENDIF()
57 58

	
58 59
  TARGET_LINK_LIBRARIES(lp_test ${LP_TEST_LIBS})
59 60
  ADD_TEST(lp_test lp_test)
60 61

	
61 62
  IF(WIN32 AND LEMON_HAVE_GLPK)
62 63
    GET_TARGET_PROPERTY(TARGET_LOC lp_test LOCATION)
63 64
    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
64 65
    ADD_CUSTOM_COMMAND(TARGET lp_test POST_BUILD
65 66
      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/glpk.dll ${TARGET_PATH}
66 67
      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/libltdl3.dll ${TARGET_PATH}
67 68
      COMMAND ${CMAKE_COMMAND} -E copy ${GLPK_BIN_DIR}/zlib1.dll ${TARGET_PATH}
68 69
    )
69 70
  ENDIF()
70 71

	
71 72
  IF(WIN32 AND LEMON_HAVE_CPLEX)
72 73
    GET_TARGET_PROPERTY(TARGET_LOC lp_test LOCATION)
73 74
    GET_FILENAME_COMPONENT(TARGET_PATH ${TARGET_LOC} PATH)
74 75
    ADD_CUSTOM_COMMAND(TARGET lp_test POST_BUILD
75 76
      COMMAND ${CMAKE_COMMAND} -E copy ${CPLEX_BIN_DIR}/cplex91.dll ${TARGET_PATH}
76 77
    )
77 78
  ENDIF()
78 79
ENDIF()
79 80

	
80 81
IF(LEMON_HAVE_MIP)
81 82
  ADD_EXECUTABLE(mip_test mip_test.cc)
82 83
  SET(MIP_TEST_LIBS lemon)
83 84

	
84 85
  IF(LEMON_HAVE_GLPK)
85 86
    SET(MIP_TEST_LIBS ${MIP_TEST_LIBS} ${GLPK_LIBRARIES})
86 87
  ENDIF()
87 88
  IF(LEMON_HAVE_CPLEX)
Ignore white space 6 line context
1 1
EXTRA_DIST += \
2 2
	test/CMakeLists.txt
3 3

	
4 4
noinst_HEADERS += \
5 5
	test/graph_test.h \
6 6
	test/test_tools.h
7 7

	
8 8
check_PROGRAMS += \
9 9
	test/adaptors_test \
10 10
	test/bellman_ford_test \
11 11
	test/bfs_test \
12 12
	test/circulation_test \
13 13
	test/connectivity_test \
14 14
	test/counter_test \
15 15
	test/dfs_test \
16 16
	test/digraph_test \
17 17
	test/dijkstra_test \
18 18
	test/dim_test \
19 19
	test/edge_set_test \
20 20
	test/error_test \
21 21
	test/euler_test \
22
	test/fractional_matching_test \
22 23
	test/gomory_hu_test \
23 24
	test/graph_copy_test \
24 25
	test/graph_test \
25 26
	test/graph_utils_test \
26 27
	test/hao_orlin_test \
27 28
	test/heap_test \
28 29
	test/kruskal_test \
29 30
	test/maps_test \
30 31
	test/matching_test \
31 32
	test/min_cost_arborescence_test \
32 33
	test/min_cost_flow_test \
33 34
	test/path_test \
34 35
	test/preflow_test \
35 36
	test/radix_sort_test \
36 37
	test/random_test \
37 38
	test/suurballe_test \
38 39
	test/test_tools_fail \
39 40
	test/test_tools_pass \
40 41
	test/time_measure_test \
41 42
	test/unionfind_test
42 43

	
43 44
test_test_tools_pass_DEPENDENCIES = demo
44 45

	
45 46
if HAVE_LP
46 47
check_PROGRAMS += test/lp_test
47 48
endif HAVE_LP
48 49
if HAVE_MIP
49 50
check_PROGRAMS += test/mip_test
50 51
endif HAVE_MIP
51 52

	
52 53
TESTS += $(check_PROGRAMS)
53 54
XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
54 55

	
55 56
test_adaptors_test_SOURCES = test/adaptors_test.cc
56 57
test_bellman_ford_test_SOURCES = test/bellman_ford_test.cc
57 58
test_bfs_test_SOURCES = test/bfs_test.cc
58 59
test_circulation_test_SOURCES = test/circulation_test.cc
59 60
test_counter_test_SOURCES = test/counter_test.cc
60 61
test_connectivity_test_SOURCES = test/connectivity_test.cc
61 62
test_dfs_test_SOURCES = test/dfs_test.cc
62 63
test_digraph_test_SOURCES = test/digraph_test.cc
63 64
test_dijkstra_test_SOURCES = test/dijkstra_test.cc
64 65
test_dim_test_SOURCES = test/dim_test.cc
65 66
test_edge_set_test_SOURCES = test/edge_set_test.cc
66 67
test_error_test_SOURCES = test/error_test.cc
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
69 71
test_graph_copy_test_SOURCES = test/graph_copy_test.cc
70 72
test_graph_test_SOURCES = test/graph_test.cc
71 73
test_graph_utils_test_SOURCES = test/graph_utils_test.cc
72 74
test_heap_test_SOURCES = test/heap_test.cc
73 75
test_kruskal_test_SOURCES = test/kruskal_test.cc
74 76
test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
75 77
test_lp_test_SOURCES = test/lp_test.cc
76 78
test_maps_test_SOURCES = test/maps_test.cc
77 79
test_mip_test_SOURCES = test/mip_test.cc
78 80
test_matching_test_SOURCES = test/matching_test.cc
79 81
test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
80 82
test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
81 83
test_path_test_SOURCES = test/path_test.cc
82 84
test_preflow_test_SOURCES = test/preflow_test.cc
83 85
test_radix_sort_test_SOURCES = test/radix_sort_test.cc
84 86
test_suurballe_test_SOURCES = test/suurballe_test.cc
85 87
test_random_test_SOURCES = test/random_test.cc
86 88
test_test_tools_fail_SOURCES = test/test_tools_fail.cc
87 89
test_test_tools_pass_SOURCES = test/test_tools_pass.cc
88 90
test_time_measure_test_SOURCES = test/time_measure_test.cc
89 91
test_unionfind_test_SOURCES = test/unionfind_test.cc
0 comments (0 inline)