[Lemon-commits] [lemon_svn] deba: r2694 - in hugo/trunk: lemon test
Lemon SVN
svn at lemon.cs.elte.hu
Mon Nov 6 20:54:32 CET 2006
Author: deba
Date: Fri Apr 14 20:07:33 2006
New Revision: 2694
Modified:
hugo/trunk/lemon/bipartite_matching.h
hugo/trunk/test/bipartite_matching_test.cc
Log:
MaxWeightedBipartiteMatching
MinCostMaxBipartiteMatching
Both algorithms are based on successive shortest
path algorithm with dijkstra shortest path
finding
Modified: hugo/trunk/lemon/bipartite_matching.h
==============================================================================
--- hugo/trunk/lemon/bipartite_matching.h (original)
+++ hugo/trunk/lemon/bipartite_matching.h Fri Apr 14 20:07:33 2006
@@ -19,8 +19,10 @@
#ifndef LEMON_BIPARTITE_MATCHING
#define LEMON_BIPARTITE_MATCHING
-#include <lemon/bpugraph_adaptor.h>
-#include <lemon/bfs.h>
+#include <functional>
+
+#include <lemon/bin_heap.h>
+#include <lemon/maps.h>
#include <iostream>
@@ -35,7 +37,7 @@
/// \brief Bipartite Max Cardinality Matching algorithm
///
/// Bipartite Max Cardinality Matching algorithm. This class implements
- /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
+ /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
/// complexity.
template <typename BpUGraph>
class MaxBipartiteMatching {
@@ -83,7 +85,7 @@
for (BNodeIt it(*graph); it != INVALID; ++it) {
bnode_matching[it] = INVALID;
}
- matching_value = 0;
+ matching_size = 0;
}
/// \brief Initalize the data structures.
@@ -92,7 +94,7 @@
/// matching. From this matching sometimes it is faster to get
/// the matching than from the initial empty matching.
void greedyInit() {
- matching_value = 0;
+ matching_size = 0;
for (BNodeIt it(*graph); it != INVALID; ++it) {
bnode_matching[it] = INVALID;
}
@@ -102,7 +104,7 @@
if (bnode_matching[graph->bNode(jt)] == INVALID) {
anode_matching[it] = jt;
bnode_matching[graph->bNode(jt)] = jt;
- ++matching_value;
+ ++matching_size;
break;
}
}
@@ -120,10 +122,10 @@
for (BNodeIt it(*graph); it != INVALID; ++it) {
bnode_matching[it] = INVALID;
}
- matching_value = 0;
+ matching_size = 0;
for (UEdgeIt it(*graph); it != INVALID; ++it) {
if (matching[it]) {
- ++matching_value;
+ ++matching_size;
anode_matching[graph->aNode(it)] = it;
bnode_matching[graph->bNode(it)] = it;
}
@@ -142,10 +144,10 @@
for (BNodeIt it(*graph); it != INVALID; ++it) {
bnode_matching[it] = INVALID;
}
- matching_value = 0;
+ matching_size = 0;
for (UEdgeIt it(*graph); it != INVALID; ++it) {
if (matching[it]) {
- ++matching_value;
+ ++matching_size;
if (anode_matching[graph->aNode(it)] != INVALID) {
return false;
}
@@ -246,7 +248,7 @@
aused[anode] = true;
}
- ++matching_value;
+ ++matching_size;
}
}
@@ -297,7 +299,7 @@
anode_matching[anode] = uedge;
}
- ++matching_value;
+ ++matching_size;
return true;
} else {
Node newanode = graph->aNode(bnode_matching[bnode]);
@@ -407,7 +409,7 @@
matching[anode_matching[it]] = true;
}
}
- return matching_value;
+ return matching_size;
}
/// \brief Set true all matching uedge in the map and the others to false.
@@ -419,7 +421,7 @@
for (UEdgeIt it(*graph); it != INVALID; ++it) {
matching[it] = it == anode_matching[graph->aNode(it)];
}
- return matching_value;
+ return matching_size;
}
@@ -445,8 +447,8 @@
/// \brief Gives back the number of the matching edges.
///
/// Gives back the number of the matching edges.
- int matchingValue() const {
- return matching_value;
+ int matchingSize() const {
+ return matching_size;
}
/// @}
@@ -457,7 +459,1022 @@
BNodeMatchingMap bnode_matching;
const Graph *graph;
- int matching_value;
+ int matching_size;
+
+ };
+
+ /// \brief Default traits class for weighted bipartite matching algoritms.
+ ///
+ /// Default traits class for weighted bipartite matching algoritms.
+ /// \param _BpUGraph The bipartite undirected graph type.
+ /// \param _WeightMap Type of weight map.
+ template <typename _BpUGraph, typename _WeightMap>
+ struct WeightedBipartiteMatchingDefaultTraits {
+ /// \brief The type of the weight of the undirected edges.
+ typedef typename _WeightMap::Value Value;
+
+ /// The undirected bipartite graph type the algorithm runs on.
+ typedef _BpUGraph BpUGraph;
+
+ /// The map of the edges weights
+ typedef _WeightMap WeightMap;
+
+ /// \brief The cross reference type used by heap.
+ ///
+ /// The cross reference type used by heap.
+ /// Usually it is \c Graph::NodeMap<int>.
+ typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
+
+ /// \brief Instantiates a HeapCrossRef.
+ ///
+ /// This function instantiates a \ref HeapCrossRef.
+ /// \param graph is the graph, to which we would like to define the
+ /// HeapCrossRef.
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+ return new HeapCrossRef(graph);
+ }
+
+ /// \brief The heap type used by weighted matching algorithms.
+ ///
+ /// The heap type used by weighted matching algorithms. It should
+ /// minimize the priorities and the heap's key type is the graph's
+ /// anode graph's node.
+ ///
+ /// \sa BinHeap
+ typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
+
+ /// \brief Instantiates a Heap.
+ ///
+ /// This function instantiates a \ref Heap.
+ /// \param crossref The cross reference of the heap.
+ static Heap *createHeap(HeapCrossRef& crossref) {
+ return new Heap(crossref);
+ }
+
+ };
+
+
+ /// \ingroup matching
+ ///
+ /// \brief Bipartite Max Weighted Matching algorithm
+ ///
+ /// This class implements the bipartite Max Weighted Matching
+ /// algorithm. It uses the successive shortest path algorithm to
+ /// calculate the maximum weighted matching in the bipartite
+ /// graph. The algorithm can be used also to calculate the maximum
+ /// cardinality maximum weighted matching. The time complexity
+ /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
+ /// heap implementation but this can be improved to
+ /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
+ ///
+ /// The algorithm also provides a potential function on the nodes
+ /// which a dual solution of the matching algorithm and it can be
+ /// used to proof the optimality of the given pimal solution.
+#ifdef DOXYGEN
+ template <typename _BpUGraph, typename _WeightMap, typename _Traits>
+#else
+ template <typename _BpUGraph,
+ typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
+ typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
+#endif
+ class MaxWeightedBipartiteMatching {
+ public:
+
+ typedef _Traits Traits;
+ typedef typename Traits::BpUGraph BpUGraph;
+ typedef typename Traits::WeightMap WeightMap;
+ typedef typename Traits::Value Value;
+
+ protected:
+
+ typedef typename Traits::HeapCrossRef HeapCrossRef;
+ typedef typename Traits::Heap Heap;
+
+
+ typedef typename BpUGraph::Node Node;
+ typedef typename BpUGraph::ANodeIt ANodeIt;
+ typedef typename BpUGraph::BNodeIt BNodeIt;
+ typedef typename BpUGraph::UEdge UEdge;
+ typedef typename BpUGraph::UEdgeIt UEdgeIt;
+ typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
+
+ typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
+ typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
+
+ typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
+ typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
+
+
+ public:
+
+ /// \brief \ref Exception for uninitialized parameters.
+ ///
+ /// This error represents problems in the initialization
+ /// of the parameters of the algorithms.
+ class UninitializedParameter : public lemon::UninitializedParameter {
+ public:
+ virtual const char* exceptionName() const {
+ return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
+ }
+ };
+
+ ///\name Named template parameters
+
+ ///@{
+
+ template <class H, class CR>
+ struct DefHeapTraits : public Traits {
+ typedef CR HeapCrossRef;
+ typedef H Heap;
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
+ throw UninitializedParameter();
+ }
+ static Heap *createHeap(HeapCrossRef &) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting heap
+ /// and cross reference type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting heap and cross
+ /// reference type
+ template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+ struct DefHeap
+ : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
+ DefHeapTraits<H, CR> > {
+ typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
+ DefHeapTraits<H, CR> > Create;
+ };
+
+ template <class H, class CR>
+ struct DefStandardHeapTraits : public Traits {
+ typedef CR HeapCrossRef;
+ typedef H Heap;
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+ return new HeapCrossRef(graph);
+ }
+ static Heap *createHeap(HeapCrossRef &crossref) {
+ return new Heap(crossref);
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting heap and
+ /// cross reference type with automatic allocation
+ ///
+ /// \ref named-templ-param "Named parameter" for setting heap and cross
+ /// reference type. It can allocate the heap and the cross reference
+ /// object if the cross reference's constructor waits for the graph as
+ /// parameter and the heap's constructor waits for the cross reference.
+ template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+ struct DefStandardHeap
+ : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
+ DefStandardHeapTraits<H, CR> > {
+ typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
+ DefStandardHeapTraits<H, CR> >
+ Create;
+ };
+
+ ///@}
+
+
+ /// \brief Constructor.
+ ///
+ /// Constructor of the algorithm.
+ MaxWeightedBipartiteMatching(const BpUGraph& _graph,
+ const WeightMap& _weight)
+ : graph(&_graph), weight(&_weight),
+ anode_matching(_graph), bnode_matching(_graph),
+ anode_potential(_graph), bnode_potential(_graph),
+ _heap_cross_ref(0), local_heap_cross_ref(false),
+ _heap(0), local_heap(0) {}
+
+ /// \brief Destructor.
+ ///
+ /// Destructor of the algorithm.
+ ~MaxWeightedBipartiteMatching() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the heap and the cross reference used by algorithm.
+ ///
+ /// Sets the heap and the cross reference used by algorithm.
+ /// If you don't use this function before calling \ref run(),
+ /// it will allocate one. The destuctor deallocates this
+ /// automatically allocated map, of course.
+ /// \return \c (*this)
+ MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
+ if(local_heap_cross_ref) {
+ delete _heap_cross_ref;
+ local_heap_cross_ref = false;
+ }
+ _heap_cross_ref = &crossRef;
+ if(local_heap) {
+ delete _heap;
+ local_heap = false;
+ }
+ _heap = &heap;
+ return *this;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to use
+ /// one of the member functions called \c run().
+ /// \n
+ /// If you need more control on the execution,
+ /// first you must call \ref init() or one alternative for it.
+ /// Finally \ref start() will perform the matching computation or
+ /// with step-by-step execution you can augment the solution.
+
+ /// @{
+
+ /// \brief Initalize the data structures.
+ ///
+ /// It initalizes the data structures and creates an empty matching.
+ void init() {
+ initStructures();
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ anode_matching[it] = INVALID;
+ anode_potential[it] = 0;
+ }
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ bnode_matching[it] = INVALID;
+ bnode_potential[it] = 0;
+ for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
+ if (-(*weight)[jt] <= bnode_potential[it]) {
+ bnode_potential[it] = -(*weight)[jt];
+ }
+ }
+ }
+ matching_value = 0;
+ matching_size = 0;
+ }
+
+
+ /// \brief An augmenting phase of the weighted matching algorithm
+ ///
+ /// It runs an augmenting phase of the weighted matching
+ /// algorithm. The phase finds the best augmenting path and
+ /// augments only on this paths.
+ ///
+ /// The algorithm consists at most
+ /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
+ /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
+ /// with binary heap.
+ /// \param decrease If the given parameter true the matching value
+ /// can be decreased in the augmenting phase. If we would like
+ /// to calculate the maximum cardinality maximum weighted matching
+ /// then we should let the algorithm to decrease the matching
+ /// value in order to increase the number of the matching edges.
+ bool augment(bool decrease = false) {
+
+ typename BpUGraph::template BNodeMap<Value> bdist(*graph);
+ typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
+
+ Node bestNode = INVALID;
+ Value bestValue = 0;
+
+ _heap->clear();
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
+ }
+
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] == INVALID) {
+ _heap->push(it, 0);
+ }
+ }
+
+ Value bdistMax = 0;
+ while (!_heap->empty()) {
+ Node anode = _heap->top();
+ Value avalue = _heap->prio();
+ _heap->pop();
+ for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
+ if (jt == anode_matching[anode]) continue;
+ Node bnode = graph->bNode(jt);
+ Value bvalue = avalue + anode_potential[anode] -
+ bnode_potential[bnode] - (*weight)[jt];
+ if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
+ bdist[bnode] = bvalue;
+ bpred[bnode] = jt;
+ }
+ if (bvalue > bdistMax) {
+ bdistMax = bvalue;
+ }
+ if (bnode_matching[bnode] != INVALID) {
+ Node newanode = graph->aNode(bnode_matching[bnode]);
+ switch (_heap->state(newanode)) {
+ case Heap::PRE_HEAP:
+ _heap->push(newanode, bvalue);
+ break;
+ case Heap::IN_HEAP:
+ if (bvalue < (*_heap)[newanode]) {
+ _heap->decrease(newanode, bvalue);
+ }
+ break;
+ case Heap::POST_HEAP:
+ break;
+ }
+ } else {
+ if (bestNode == INVALID ||
+ - bvalue - bnode_potential[bnode] > bestValue) {
+ bestValue = - bvalue - bnode_potential[bnode];
+ bestNode = bnode;
+ }
+ }
+ }
+ }
+
+ if (bestNode == INVALID || (!decrease && bestValue < 0)) {
+ return false;
+ }
+
+ matching_value += bestValue;
+ ++matching_size;
+
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ if (bpred[it] != INVALID) {
+ bnode_potential[it] += bdist[it];
+ } else {
+ bnode_potential[it] += bdistMax;
+ }
+ }
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] != INVALID) {
+ Node bnode = graph->bNode(anode_matching[it]);
+ if (bpred[bnode] != INVALID) {
+ anode_potential[it] += bdist[bnode];
+ } else {
+ anode_potential[it] += bdistMax;
+ }
+ }
+ }
+
+ while (bestNode != INVALID) {
+ UEdge uedge = bpred[bestNode];
+ Node anode = graph->aNode(uedge);
+
+ bnode_matching[bestNode] = uedge;
+ if (anode_matching[anode] != INVALID) {
+ bestNode = graph->bNode(anode_matching[anode]);
+ } else {
+ bestNode = INVALID;
+ }
+ anode_matching[anode] = uedge;
+ }
+
+
+ return true;
+ }
+
+ /// \brief Starts the algorithm.
+ ///
+ /// Starts the algorithm. It runs augmenting phases until the
+ /// optimal solution reached.
+ ///
+ /// \param maxCardinality If the given value is true it will
+ /// calculate the maximum cardinality maximum matching instead of
+ /// the maximum matching.
+ void start(bool maxCardinality = false) {
+ while (augment(maxCardinality)) {}
+ }
+
+ /// \brief Runs the algorithm.
+ ///
+ /// It just initalize the algorithm and then start it.
+ ///
+ /// \param maxCardinality If the given value is true it will
+ /// calculate the maximum cardinality maximum matching instead of
+ /// the maximum matching.
+ void run(bool maxCardinality = false) {
+ init();
+ start(maxCardinality);
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The result of the %Matching algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
+
+ /// \brief Gives back the potential in the NodeMap
+ ///
+ /// Gives back the potential in the NodeMap. The potential
+ /// is feasible if \f$ \pi(a) - \pi(b) - w(ab) = 0 \f$ for
+ /// each matching edges and \f$ \pi(a) - \pi(b) - w(ab) \ge 0 \f$
+ /// for each edges.
+ template <typename PotentialMap>
+ void potential(PotentialMap& potential) {
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ potential[it] = anode_potential[it];
+ }
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ potential[it] = bnode_potential[it];
+ }
+ }
+
+ /// \brief Set true all matching uedge in the map.
+ ///
+ /// Set true all matching uedge in the map. It does not change the
+ /// value mapped to the other uedges.
+ /// \return The number of the matching edges.
+ template <typename MatchingMap>
+ int quickMatching(MatchingMap& matching) {
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] != INVALID) {
+ matching[anode_matching[it]] = true;
+ }
+ }
+ return matching_size;
+ }
+
+ /// \brief Set true all matching uedge in the map and the others to false.
+ ///
+ /// Set true all matching uedge in the map and the others to false.
+ /// \return The number of the matching edges.
+ template <typename MatchingMap>
+ int matching(MatchingMap& matching) {
+ for (UEdgeIt it(*graph); it != INVALID; ++it) {
+ matching[it] = it == anode_matching[graph->aNode(it)];
+ }
+ return matching_size;
+ }
+
+
+ /// \brief Return true if the given uedge is in the matching.
+ ///
+ /// It returns true if the given uedge is in the matching.
+ bool matchingEdge(const UEdge& edge) {
+ return anode_matching[graph->aNode(edge)] == edge;
+ }
+
+ /// \brief Returns the matching edge from the node.
+ ///
+ /// Returns the matching edge from the node. If there is not such
+ /// edge it gives back \c INVALID.
+ UEdge matchingEdge(const Node& node) {
+ if (graph->aNode(node)) {
+ return anode_matching[node];
+ } else {
+ return bnode_matching[node];
+ }
+ }
+
+ /// \brief Gives back the sum of weights of the matching edges.
+ ///
+ /// Gives back the sum of weights of the matching edges.
+ Value matchingValue() const {
+ return matching_value;
+ }
+
+ /// \brief Gives back the number of the matching edges.
+ ///
+ /// Gives back the number of the matching edges.
+ int matchingSize() const {
+ return matching_size;
+ }
+
+ /// @}
+
+ private:
+
+ void initStructures() {
+ if (!_heap_cross_ref) {
+ local_heap_cross_ref = true;
+ _heap_cross_ref = Traits::createHeapCrossRef(*graph);
+ }
+ if (!_heap) {
+ local_heap = true;
+ _heap = Traits::createHeap(*_heap_cross_ref);
+ }
+ }
+
+ void destroyStructures() {
+ if (local_heap_cross_ref) delete _heap_cross_ref;
+ if (local_heap) delete _heap;
+ }
+
+
+ private:
+
+ const BpUGraph *graph;
+ const WeightMap* weight;
+
+ ANodeMatchingMap anode_matching;
+ BNodeMatchingMap bnode_matching;
+
+ ANodePotentialMap anode_potential;
+ BNodePotentialMap bnode_potential;
+
+ Value matching_value;
+ int matching_size;
+
+ HeapCrossRef *_heap_cross_ref;
+ bool local_heap_cross_ref;
+
+ Heap *_heap;
+ bool local_heap;
+
+ };
+
+ /// \brief Default traits class for minimum cost bipartite matching
+ /// algoritms.
+ ///
+ /// Default traits class for minimum cost bipartite matching
+ /// algoritms.
+ ///
+ /// \param _BpUGraph The bipartite undirected graph
+ /// type.
+ ///
+ /// \param _CostMap Type of cost map.
+ template <typename _BpUGraph, typename _CostMap>
+ struct MinCostMaxBipartiteMatchingDefaultTraits {
+ /// \brief The type of the cost of the undirected edges.
+ typedef typename _CostMap::Value Value;
+
+ /// The undirected bipartite graph type the algorithm runs on.
+ typedef _BpUGraph BpUGraph;
+
+ /// The map of the edges costs
+ typedef _CostMap CostMap;
+
+ /// \brief The cross reference type used by heap.
+ ///
+ /// The cross reference type used by heap.
+ /// Usually it is \c Graph::NodeMap<int>.
+ typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
+
+ /// \brief Instantiates a HeapCrossRef.
+ ///
+ /// This function instantiates a \ref HeapCrossRef.
+ /// \param graph is the graph, to which we would like to define the
+ /// HeapCrossRef.
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+ return new HeapCrossRef(graph);
+ }
+
+ /// \brief The heap type used by costed matching algorithms.
+ ///
+ /// The heap type used by costed matching algorithms. It should
+ /// minimize the priorities and the heap's key type is the graph's
+ /// anode graph's node.
+ ///
+ /// \sa BinHeap
+ typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
+
+ /// \brief Instantiates a Heap.
+ ///
+ /// This function instantiates a \ref Heap.
+ /// \param crossref The cross reference of the heap.
+ static Heap *createHeap(HeapCrossRef& crossref) {
+ return new Heap(crossref);
+ }
+
+ };
+
+
+ /// \ingroup matching
+ ///
+ /// \brief Bipartite Min Cost Matching algorithm
+ ///
+ /// This class implements the bipartite Min Cost Matching algorithm.
+ /// It uses the successive shortest path algorithm to calculate the
+ /// minimum cost maximum matching in the bipartite graph. The time
+ /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
+ /// default binary heap implementation but this can be improved to
+ /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
+ ///
+ /// The algorithm also provides a potential function on the nodes
+ /// which a dual solution of the matching algorithm and it can be
+ /// used to proof the optimality of the given pimal solution.
+#ifdef DOXYGEN
+ template <typename _BpUGraph, typename _CostMap, typename _Traits>
+#else
+ template <typename _BpUGraph,
+ typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
+ typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
+#endif
+ class MinCostMaxBipartiteMatching {
+ public:
+
+ typedef _Traits Traits;
+ typedef typename Traits::BpUGraph BpUGraph;
+ typedef typename Traits::CostMap CostMap;
+ typedef typename Traits::Value Value;
+
+ protected:
+
+ typedef typename Traits::HeapCrossRef HeapCrossRef;
+ typedef typename Traits::Heap Heap;
+
+
+ typedef typename BpUGraph::Node Node;
+ typedef typename BpUGraph::ANodeIt ANodeIt;
+ typedef typename BpUGraph::BNodeIt BNodeIt;
+ typedef typename BpUGraph::UEdge UEdge;
+ typedef typename BpUGraph::UEdgeIt UEdgeIt;
+ typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
+
+ typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
+ typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
+
+ typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
+ typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
+
+
+ public:
+
+ /// \brief \ref Exception for uninitialized parameters.
+ ///
+ /// This error represents problems in the initialization
+ /// of the parameters of the algorithms.
+ class UninitializedParameter : public lemon::UninitializedParameter {
+ public:
+ virtual const char* exceptionName() const {
+ return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
+ }
+ };
+
+ ///\name Named template parameters
+
+ ///@{
+
+ template <class H, class CR>
+ struct DefHeapTraits : public Traits {
+ typedef CR HeapCrossRef;
+ typedef H Heap;
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
+ throw UninitializedParameter();
+ }
+ static Heap *createHeap(HeapCrossRef &) {
+ throw UninitializedParameter();
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting heap
+ /// and cross reference type
+ ///
+ /// \ref named-templ-param "Named parameter" for setting heap and cross
+ /// reference type
+ template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+ struct DefHeap
+ : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
+ DefHeapTraits<H, CR> > {
+ typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
+ DefHeapTraits<H, CR> > Create;
+ };
+
+ template <class H, class CR>
+ struct DefStandardHeapTraits : public Traits {
+ typedef CR HeapCrossRef;
+ typedef H Heap;
+ static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
+ return new HeapCrossRef(graph);
+ }
+ static Heap *createHeap(HeapCrossRef &crossref) {
+ return new Heap(crossref);
+ }
+ };
+
+ /// \brief \ref named-templ-param "Named parameter" for setting heap and
+ /// cross reference type with automatic allocation
+ ///
+ /// \ref named-templ-param "Named parameter" for setting heap and cross
+ /// reference type. It can allocate the heap and the cross reference
+ /// object if the cross reference's constructor waits for the graph as
+ /// parameter and the heap's constructor waits for the cross reference.
+ template <class H, class CR = typename BpUGraph::template NodeMap<int> >
+ struct DefStandardHeap
+ : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
+ DefStandardHeapTraits<H, CR> > {
+ typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
+ DefStandardHeapTraits<H, CR> >
+ Create;
+ };
+
+ ///@}
+
+
+ /// \brief Constructor.
+ ///
+ /// Constructor of the algorithm.
+ MinCostMaxBipartiteMatching(const BpUGraph& _graph,
+ const CostMap& _cost)
+ : graph(&_graph), cost(&_cost),
+ anode_matching(_graph), bnode_matching(_graph),
+ anode_potential(_graph), bnode_potential(_graph),
+ _heap_cross_ref(0), local_heap_cross_ref(false),
+ _heap(0), local_heap(0) {}
+
+ /// \brief Destructor.
+ ///
+ /// Destructor of the algorithm.
+ ~MinCostMaxBipartiteMatching() {
+ destroyStructures();
+ }
+
+ /// \brief Sets the heap and the cross reference used by algorithm.
+ ///
+ /// Sets the heap and the cross reference used by algorithm.
+ /// If you don't use this function before calling \ref run(),
+ /// it will allocate one. The destuctor deallocates this
+ /// automatically allocated map, of course.
+ /// \return \c (*this)
+ MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
+ if(local_heap_cross_ref) {
+ delete _heap_cross_ref;
+ local_heap_cross_ref = false;
+ }
+ _heap_cross_ref = &crossRef;
+ if(local_heap) {
+ delete _heap;
+ local_heap = false;
+ }
+ _heap = &heap;
+ return *this;
+ }
+
+ /// \name Execution control
+ /// The simplest way to execute the algorithm is to use
+ /// one of the member functions called \c run().
+ /// \n
+ /// If you need more control on the execution,
+ /// first you must call \ref init() or one alternative for it.
+ /// Finally \ref start() will perform the matching computation or
+ /// with step-by-step execution you can augment the solution.
+
+ /// @{
+
+ /// \brief Initalize the data structures.
+ ///
+ /// It initalizes the data structures and creates an empty matching.
+ void init() {
+ initStructures();
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ anode_matching[it] = INVALID;
+ anode_potential[it] = 0;
+ }
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ bnode_matching[it] = INVALID;
+ bnode_potential[it] = 0;
+ }
+ matching_cost = 0;
+ matching_size = 0;
+ }
+
+
+ /// \brief An augmenting phase of the costed matching algorithm
+ ///
+ /// It runs an augmenting phase of the matching algorithm. The
+ /// phase finds the best augmenting path and augments only on this
+ /// paths.
+ ///
+ /// The algorithm consists at most
+ /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
+ /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
+ /// with binary heap.
+ bool augment() {
+
+ typename BpUGraph::template BNodeMap<Value> bdist(*graph);
+ typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
+
+ Node bestNode = INVALID;
+ Value bestValue = 0;
+
+ _heap->clear();
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
+ }
+
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] == INVALID) {
+ _heap->push(it, 0);
+ }
+ }
+
+ while (!_heap->empty()) {
+ Node anode = _heap->top();
+ Value avalue = _heap->prio();
+ _heap->pop();
+ for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
+ if (jt == anode_matching[anode]) continue;
+ Node bnode = graph->bNode(jt);
+ Value bvalue = avalue + (*cost)[jt] +
+ anode_potential[anode] - bnode_potential[bnode];
+ if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
+ bdist[bnode] = bvalue;
+ bpred[bnode] = jt;
+ }
+ if (bnode_matching[bnode] != INVALID) {
+ Node newanode = graph->aNode(bnode_matching[bnode]);
+ switch (_heap->state(newanode)) {
+ case Heap::PRE_HEAP:
+ _heap->push(newanode, bvalue);
+ break;
+ case Heap::IN_HEAP:
+ if (bvalue < (*_heap)[newanode]) {
+ _heap->decrease(newanode, bvalue);
+ }
+ break;
+ case Heap::POST_HEAP:
+ break;
+ }
+ } else {
+ if (bestNode == INVALID ||
+ bvalue + bnode_potential[bnode] < bestValue) {
+ bestValue = bvalue + bnode_potential[bnode];
+ bestNode = bnode;
+ }
+ }
+ }
+ }
+
+ if (bestNode == INVALID) {
+ return false;
+ }
+
+ matching_cost += bestValue;
+ ++matching_size;
+
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ if (bpred[it] != INVALID) {
+ bnode_potential[it] += bdist[it];
+ }
+ }
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] != INVALID) {
+ Node bnode = graph->bNode(anode_matching[it]);
+ if (bpred[bnode] != INVALID) {
+ anode_potential[it] += bdist[bnode];
+ }
+ }
+ }
+
+ while (bestNode != INVALID) {
+ UEdge uedge = bpred[bestNode];
+ Node anode = graph->aNode(uedge);
+
+ bnode_matching[bestNode] = uedge;
+ if (anode_matching[anode] != INVALID) {
+ bestNode = graph->bNode(anode_matching[anode]);
+ } else {
+ bestNode = INVALID;
+ }
+ anode_matching[anode] = uedge;
+ }
+
+
+ return true;
+ }
+
+ /// \brief Starts the algorithm.
+ ///
+ /// Starts the algorithm. It runs augmenting phases until the
+ /// optimal solution reached.
+ void start() {
+ while (augment()) {}
+ }
+
+ /// \brief Runs the algorithm.
+ ///
+ /// It just initalize the algorithm and then start it.
+ void run() {
+ init();
+ start();
+ }
+
+ /// @}
+
+ /// \name Query Functions
+ /// The result of the %Matching algorithm can be obtained using these
+ /// functions.\n
+ /// Before the use of these functions,
+ /// either run() or start() must be called.
+
+ ///@{
+
+ /// \brief Gives back the potential in the NodeMap
+ ///
+ /// Gives back the potential in the NodeMap. The potential
+ /// is feasible if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
+ /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
+ /// for each edges.
+ template <typename PotentialMap>
+ void potential(PotentialMap& potential) {
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ potential[it] = anode_potential[it];
+ }
+ for (BNodeIt it(*graph); it != INVALID; ++it) {
+ potential[it] = bnode_potential[it];
+ }
+ }
+
+ /// \brief Set true all matching uedge in the map.
+ ///
+ /// Set true all matching uedge in the map. It does not change the
+ /// value mapped to the other uedges.
+ /// \return The number of the matching edges.
+ template <typename MatchingMap>
+ int quickMatching(MatchingMap& matching) {
+ for (ANodeIt it(*graph); it != INVALID; ++it) {
+ if (anode_matching[it] != INVALID) {
+ matching[anode_matching[it]] = true;
+ }
+ }
+ return matching_size;
+ }
+
+ /// \brief Set true all matching uedge in the map and the others to false.
+ ///
+ /// Set true all matching uedge in the map and the others to false.
+ /// \return The number of the matching edges.
+ template <typename MatchingMap>
+ int matching(MatchingMap& matching) {
+ for (UEdgeIt it(*graph); it != INVALID; ++it) {
+ matching[it] = it == anode_matching[graph->aNode(it)];
+ }
+ return matching_size;
+ }
+
+
+ /// \brief Return true if the given uedge is in the matching.
+ ///
+ /// It returns true if the given uedge is in the matching.
+ bool matchingEdge(const UEdge& edge) {
+ return anode_matching[graph->aNode(edge)] == edge;
+ }
+
+ /// \brief Returns the matching edge from the node.
+ ///
+ /// Returns the matching edge from the node. If there is not such
+ /// edge it gives back \c INVALID.
+ UEdge matchingEdge(const Node& node) {
+ if (graph->aNode(node)) {
+ return anode_matching[node];
+ } else {
+ return bnode_matching[node];
+ }
+ }
+
+ /// \brief Gives back the sum of costs of the matching edges.
+ ///
+ /// Gives back the sum of costs of the matching edges.
+ Value matchingCost() const {
+ return matching_cost;
+ }
+
+ /// \brief Gives back the number of the matching edges.
+ ///
+ /// Gives back the number of the matching edges.
+ int matchingSize() const {
+ return matching_size;
+ }
+
+ /// @}
+
+ private:
+
+ void initStructures() {
+ if (!_heap_cross_ref) {
+ local_heap_cross_ref = true;
+ _heap_cross_ref = Traits::createHeapCrossRef(*graph);
+ }
+ if (!_heap) {
+ local_heap = true;
+ _heap = Traits::createHeap(*_heap_cross_ref);
+ }
+ }
+
+ void destroyStructures() {
+ if (local_heap_cross_ref) delete _heap_cross_ref;
+ if (local_heap) delete _heap;
+ }
+
+
+ private:
+
+ const BpUGraph *graph;
+ const CostMap* cost;
+
+ ANodeMatchingMap anode_matching;
+ BNodeMatchingMap bnode_matching;
+
+ ANodePotentialMap anode_potential;
+ BNodePotentialMap bnode_potential;
+
+ Value matching_cost;
+ int matching_size;
+
+ HeapCrossRef *_heap_cross_ref;
+ bool local_heap_cross_ref;
+
+ Heap *_heap;
+ bool local_heap;
};
Modified: hugo/trunk/test/bipartite_matching_test.cc
==============================================================================
--- hugo/trunk/test/bipartite_matching_test.cc (original)
+++ hugo/trunk/test/bipartite_matching_test.cc Fri Apr 14 20:07:33 2006
@@ -9,6 +9,8 @@
#include <lemon/graph_utils.h>
+#include <lemon/maps.h>
+
#include "test_tools.h"
using namespace std;
@@ -24,6 +26,13 @@
int n = argc > 1 ? atoi(argv[1]) : 100;
int m = argc > 2 ? atoi(argv[2]) : 100;
int e = argc > 3 ? atoi(argv[3]) : (int)((n+m) * log(n+m));
+ int c = argc > 4 ? atoi(argv[4]) : 100;
+
+ Graph::UEdgeMap<int> weight(graph);
+
+ int max_cardinality;
+ int max_weight;
+ int max_cardinality_max_weight;
for (int i = 0; i < n; ++i) {
Node node = graph.addANode();
@@ -36,7 +45,8 @@
for (int i = 0; i < e; ++i) {
Node aNode = aNodes[urandom(n)];
Node bNode = bNodes[urandom(m)];
- graph.addEdge(aNode, bNode);
+ UEdge uedge = graph.addEdge(aNode, bNode);
+ weight[uedge] = urandom(c);
}
{
@@ -57,9 +67,10 @@
int num = 0;
for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
if (mm[jt]) ++num;
- }
+ }
check(num <= 1, "INVALID PRIMAL");
}
+ max_cardinality = bpmatch.matchingSize();
}
{
@@ -111,20 +122,131 @@
}
{
- SwapBpUGraphAdaptor<Graph> swapgraph(graph);
- MaxBipartiteMatching<SwapBpUGraphAdaptor<Graph> > bpmatch(swapgraph);
+ MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
- bpmatch.run();
+ bpmatch.init();
+ max_weight = 0;
+ while (bpmatch.augment(true)) {
+
+ Graph::UEdgeMap<bool> mm(graph);
+ Graph::NodeMap<int> pm(graph);
+
+ bpmatch.matching(mm);
+ bpmatch.potential(pm);
+
+ for (UEdgeIt it(graph); it != INVALID; ++it) {
+ if (mm[it]) {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+ "INVALID DUAL");
+ } else {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+ "INVALID DUAL");
+ }
+ }
+
+ for (ANodeIt it(graph); it != INVALID; ++it) {
+ int num = 0;
+ for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+ if (mm[jt]) ++num;
+ }
+ check(num <= 1, "INVALID PRIMAL");
+ }
+ if (bpmatch.matchingValue() > max_weight) {
+ max_weight = bpmatch.matchingValue();
+ }
+ }
+ }
+
+ {
+ MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
+
+ bpmatch.run();
+
Graph::UEdgeMap<bool> mm(graph);
- Graph::NodeMap<bool> cs(graph);
+ Graph::NodeMap<int> pm(graph);
+
+ bpmatch.matching(mm);
+ bpmatch.potential(pm);
- check(bpmatch.coverSet(cs) == bpmatch.matching(mm), "INVALID PRIMAL-DUAL");
+ for (UEdgeIt it(graph); it != INVALID; ++it) {
+ if (mm[it]) {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+ "INVALID DUAL");
+ } else {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+ "INVALID DUAL");
+ }
+ }
+
+ for (ANodeIt it(graph); it != INVALID; ++it) {
+ int num = 0;
+ for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+ if (mm[jt]) ++num;
+ }
+ check(num <= 1, "INVALID PRIMAL");
+ }
+ check(max_weight == bpmatch.matchingValue(), "WRONG WEIGHT");
+ }
+
+ {
+ MaxWeightedBipartiteMatching<Graph> bpmatch(graph, weight);
+
+ bpmatch.run(true);
+
+ Graph::UEdgeMap<bool> mm(graph);
+ Graph::NodeMap<int> pm(graph);
+
+ bpmatch.matching(mm);
+ bpmatch.potential(pm);
for (UEdgeIt it(graph); it != INVALID; ++it) {
- check(cs[graph.aNode(it)] || cs[graph.bNode(it)], "INVALID DUAL");
+ if (mm[it]) {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] == 0,
+ "INVALID DUAL");
+ } else {
+ check(pm[graph.aNode(it)] - weight[it] - pm[graph.bNode(it)] >= 0,
+ "INVALID DUAL");
+ }
}
+
+ for (ANodeIt it(graph); it != INVALID; ++it) {
+ int num = 0;
+ for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+ if (mm[jt]) ++num;
+ }
+ check(num <= 1, "INVALID PRIMAL");
+ }
+ check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
+ max_cardinality_max_weight = bpmatch.matchingValue();
+
+ }
+
+ {
+ Graph::UEdgeMap<int> cost(graph);
+
+ cost = subMap(constMap<UEdge>(c), weight);
+
+ MinCostMaxBipartiteMatching<Graph> bpmatch(graph, cost);
+
+ bpmatch.run();
+
+ Graph::UEdgeMap<bool> mm(graph);
+ Graph::NodeMap<int> pm(graph);
+
+ bpmatch.matching(mm);
+ bpmatch.potential(pm);
+ for (UEdgeIt it(graph); it != INVALID; ++it) {
+ if (mm[it]) {
+ check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] == 0,
+ "INVALID DUAL");
+ } else {
+ check(pm[graph.aNode(it)] + cost[it] - pm[graph.bNode(it)] >= 0,
+ "INVALID DUAL");
+ }
+ }
+
for (ANodeIt it(graph); it != INVALID; ++it) {
int num = 0;
for (IncEdgeIt jt(graph, it); jt != INVALID; ++jt) {
@@ -132,6 +254,11 @@
}
check(num <= 1, "INVALID PRIMAL");
}
+
+ check(max_cardinality == bpmatch.matchingSize(), "WRONG SIZE");
+ check(max_cardinality * c - max_cardinality_max_weight
+ == bpmatch.matchingCost(), "WRONG SIZE");
+
}
return 0;
More information about the Lemon-commits
mailing list