lemon/bipartite_matching.h
author deba
Fri, 02 Mar 2007 18:04:28 +0000
changeset 2386 81b47fc5c444
parent 2352 5e273e0bd5e2
child 2391 14a343be7a5a
permissions -rw-r--r--
Hard Warning checking
- based on the remark of the ZIB user
- we do not use -Winline
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2006
     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_BIPARTITE_MATCHING
    20 #define LEMON_BIPARTITE_MATCHING
    21 
    22 #include <functional>
    23 
    24 #include <lemon/bin_heap.h>
    25 #include <lemon/maps.h>
    26 
    27 #include <iostream>
    28 
    29 ///\ingroup matching
    30 ///\file
    31 ///\brief Maximum matching algorithms in bipartite graphs.
    32 
    33 namespace lemon {
    34 
    35   /// \ingroup matching
    36   ///
    37   /// \brief Bipartite Max Cardinality Matching algorithm
    38   ///
    39   /// Bipartite Max Cardinality Matching algorithm. This class implements
    40   /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
    41   /// complexity.
    42   template <typename BpUGraph>
    43   class MaxBipartiteMatching {
    44   protected:
    45 
    46     typedef BpUGraph Graph;
    47 
    48     typedef typename Graph::Node Node;
    49     typedef typename Graph::ANodeIt ANodeIt;
    50     typedef typename Graph::BNodeIt BNodeIt;
    51     typedef typename Graph::UEdge UEdge;
    52     typedef typename Graph::UEdgeIt UEdgeIt;
    53     typedef typename Graph::IncEdgeIt IncEdgeIt;
    54 
    55     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
    56     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
    57 
    58 
    59   public:
    60 
    61     /// \brief Constructor.
    62     ///
    63     /// Constructor of the algorithm. 
    64     MaxBipartiteMatching(const BpUGraph& _graph) 
    65       : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
    66 
    67     /// \name Execution control
    68     /// The simplest way to execute the algorithm is to use
    69     /// one of the member functions called \c run().
    70     /// \n
    71     /// If you need more control on the execution,
    72     /// first you must call \ref init() or one alternative for it.
    73     /// Finally \ref start() will perform the matching computation or
    74     /// with step-by-step execution you can augment the solution.
    75 
    76     /// @{
    77 
    78     /// \brief Initalize the data structures.
    79     ///
    80     /// It initalizes the data structures and creates an empty matching.
    81     void init() {
    82       for (ANodeIt it(*graph); it != INVALID; ++it) {
    83         anode_matching[it] = INVALID;
    84       }
    85       for (BNodeIt it(*graph); it != INVALID; ++it) {
    86         bnode_matching[it] = INVALID;
    87       }
    88       matching_size = 0;
    89     }
    90 
    91     /// \brief Initalize the data structures.
    92     ///
    93     /// It initalizes the data structures and creates a greedy
    94     /// matching.  From this matching sometimes it is faster to get
    95     /// the matching than from the initial empty matching.
    96     void greedyInit() {
    97       matching_size = 0;
    98       for (BNodeIt it(*graph); it != INVALID; ++it) {
    99         bnode_matching[it] = INVALID;
   100       }
   101       for (ANodeIt it(*graph); it != INVALID; ++it) {
   102         anode_matching[it] = INVALID;
   103         for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
   104           if (bnode_matching[graph->bNode(jt)] == INVALID) {
   105             anode_matching[it] = jt;
   106             bnode_matching[graph->bNode(jt)] = jt;
   107             ++matching_size;
   108             break;
   109           }
   110         }
   111       }
   112     }
   113 
   114     /// \brief Initalize the data structures with an initial matching.
   115     ///
   116     /// It initalizes the data structures with an initial matching.
   117     template <typename MatchingMap>
   118     void matchingInit(const MatchingMap& mm) {
   119       for (ANodeIt it(*graph); it != INVALID; ++it) {
   120         anode_matching[it] = INVALID;
   121       }
   122       for (BNodeIt it(*graph); it != INVALID; ++it) {
   123         bnode_matching[it] = INVALID;
   124       }
   125       matching_size = 0;
   126       for (UEdgeIt it(*graph); it != INVALID; ++it) {
   127         if (mm[it]) {
   128           ++matching_size;
   129           anode_matching[graph->aNode(it)] = it;
   130           bnode_matching[graph->bNode(it)] = it;
   131         }
   132       }
   133     }
   134 
   135     /// \brief Initalize the data structures with an initial matching.
   136     ///
   137     /// It initalizes the data structures with an initial matching.
   138     /// \return %True when the given map contains really a matching.
   139     template <typename MatchingMap>
   140     void checkedMatchingInit(const MatchingMap& mm) {
   141       for (ANodeIt it(*graph); it != INVALID; ++it) {
   142         anode_matching[it] = INVALID;
   143       }
   144       for (BNodeIt it(*graph); it != INVALID; ++it) {
   145         bnode_matching[it] = INVALID;
   146       }
   147       matching_size = 0;
   148       for (UEdgeIt it(*graph); it != INVALID; ++it) {
   149         if (mm[it]) {
   150           ++matching_size;
   151           if (anode_matching[graph->aNode(it)] != INVALID) {
   152             return false;
   153           }
   154           anode_matching[graph->aNode(it)] = it;
   155           if (bnode_matching[graph->aNode(it)] != INVALID) {
   156             return false;
   157           }
   158           bnode_matching[graph->bNode(it)] = it;
   159         }
   160       }
   161       return false;
   162     }
   163 
   164     /// \brief An augmenting phase of the Hopcroft-Karp algorithm
   165     ///
   166     /// It runs an augmenting phase of the Hopcroft-Karp
   167     /// algorithm. This phase finds maximum count of edge disjoint
   168     /// augmenting paths and augments on these paths. The algorithm
   169     /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
   170     /// \f$ O(e) \f$ long.
   171     bool augment() {
   172 
   173       typename Graph::template ANodeMap<bool> areached(*graph, false);
   174       typename Graph::template BNodeMap<bool> breached(*graph, false);
   175 
   176       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
   177 
   178       std::vector<Node> queue, bqueue;
   179       for (ANodeIt it(*graph); it != INVALID; ++it) {
   180         if (anode_matching[it] == INVALID) {
   181           queue.push_back(it);
   182           areached[it] = true;
   183         }
   184       }
   185 
   186       bool success = false;
   187 
   188       while (!success && !queue.empty()) {
   189         std::vector<Node> newqueue;
   190         for (int i = 0; i < int(queue.size()); ++i) {
   191           Node anode = queue[i];
   192           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   193             Node bnode = graph->bNode(jt);
   194             if (breached[bnode]) continue;
   195             breached[bnode] = true;
   196             bpred[bnode] = jt;
   197             if (bnode_matching[bnode] == INVALID) {
   198               bqueue.push_back(bnode);
   199               success = true;
   200             } else {           
   201               Node newanode = graph->aNode(bnode_matching[bnode]);
   202               if (!areached[newanode]) {
   203                 areached[newanode] = true;
   204                 newqueue.push_back(newanode);
   205               }
   206             }
   207           }
   208         }
   209         queue.swap(newqueue);
   210       }
   211 
   212       if (success) {
   213 
   214         typename Graph::template ANodeMap<bool> aused(*graph, false);
   215         
   216         for (int i = 0; i < int(bqueue.size()); ++i) {
   217           Node bnode = bqueue[i];
   218 
   219           bool used = false;
   220 
   221           while (bnode != INVALID) {
   222             UEdge uedge = bpred[bnode];
   223             Node anode = graph->aNode(uedge);
   224             
   225             if (aused[anode]) {
   226               used = true;
   227               break;
   228             }
   229             
   230             bnode = anode_matching[anode] != INVALID ? 
   231               graph->bNode(anode_matching[anode]) : INVALID;
   232             
   233           }
   234           
   235           if (used) continue;
   236 
   237           bnode = bqueue[i];
   238           while (bnode != INVALID) {
   239             UEdge uedge = bpred[bnode];
   240             Node anode = graph->aNode(uedge);
   241             
   242             bnode_matching[bnode] = uedge;
   243 
   244             bnode = anode_matching[anode] != INVALID ? 
   245               graph->bNode(anode_matching[anode]) : INVALID;
   246             
   247             anode_matching[anode] = uedge;
   248 
   249             aused[anode] = true;
   250           }
   251           ++matching_size;
   252 
   253         }
   254       }
   255       return success;
   256     }
   257 
   258     /// \brief An augmenting phase of the Ford-Fulkerson algorithm
   259     ///
   260     /// It runs an augmenting phase of the Ford-Fulkerson
   261     /// algorithm. This phase finds only one augmenting path and 
   262     /// augments only on this paths. The algorithm consists at most 
   263     /// of \f$ O(n) \f$ simple phase and one phase is at most 
   264     /// \f$ O(e) \f$ long.
   265     bool simpleAugment() {
   266 
   267       typename Graph::template ANodeMap<bool> areached(*graph, false);
   268       typename Graph::template BNodeMap<bool> breached(*graph, false);
   269 
   270       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
   271 
   272       std::vector<Node> queue;
   273       for (ANodeIt it(*graph); it != INVALID; ++it) {
   274         if (anode_matching[it] == INVALID) {
   275           queue.push_back(it);
   276           areached[it] = true;
   277         }
   278       }
   279 
   280       while (!queue.empty()) {
   281         std::vector<Node> newqueue;
   282         for (int i = 0; i < int(queue.size()); ++i) {
   283           Node anode = queue[i];
   284           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   285             Node bnode = graph->bNode(jt);
   286             if (breached[bnode]) continue;
   287             breached[bnode] = true;
   288             bpred[bnode] = jt;
   289             if (bnode_matching[bnode] == INVALID) {
   290               while (bnode != INVALID) {
   291                 UEdge uedge = bpred[bnode];
   292                 anode = graph->aNode(uedge);
   293                 
   294                 bnode_matching[bnode] = uedge;
   295                 
   296                 bnode = anode_matching[anode] != INVALID ? 
   297                   graph->bNode(anode_matching[anode]) : INVALID;
   298                 
   299                 anode_matching[anode] = uedge;
   300                 
   301               }
   302               ++matching_size;
   303               return true;
   304             } else {           
   305               Node newanode = graph->aNode(bnode_matching[bnode]);
   306               if (!areached[newanode]) {
   307                 areached[newanode] = true;
   308                 newqueue.push_back(newanode);
   309               }
   310             }
   311           }
   312         }
   313         queue.swap(newqueue);
   314       }
   315       
   316       return false;
   317     }
   318 
   319     /// \brief Starts the algorithm.
   320     ///
   321     /// Starts the algorithm. It runs augmenting phases until the optimal
   322     /// solution reached.
   323     void start() {
   324       while (augment()) {}
   325     }
   326 
   327     /// \brief Runs the algorithm.
   328     ///
   329     /// It just initalize the algorithm and then start it.
   330     void run() {
   331       greedyInit();
   332       start();
   333     }
   334 
   335     /// @}
   336 
   337     /// \name Query Functions
   338     /// The result of the %Matching algorithm can be obtained using these
   339     /// functions.\n
   340     /// Before the use of these functions,
   341     /// either run() or start() must be called.
   342     
   343     ///@{
   344 
   345     /// \brief Returns an minimum covering of the nodes.
   346     ///
   347     /// The minimum covering set problem is the dual solution of the
   348     /// maximum bipartite matching. It provides an solution for this
   349     /// problem what is proof of the optimality of the matching.
   350     /// \return The size of the cover set.
   351     template <typename CoverMap>
   352     int coverSet(CoverMap& covering) const {
   353 
   354       typename Graph::template ANodeMap<bool> areached(*graph, false);
   355       typename Graph::template BNodeMap<bool> breached(*graph, false);
   356       
   357       std::vector<Node> queue;
   358       for (ANodeIt it(*graph); it != INVALID; ++it) {
   359         if (anode_matching[it] == INVALID) {
   360           queue.push_back(it);
   361         }
   362       }
   363 
   364       while (!queue.empty()) {
   365         std::vector<Node> newqueue;
   366         for (int i = 0; i < int(queue.size()); ++i) {
   367           Node anode = queue[i];
   368           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   369             Node bnode = graph->bNode(jt);
   370             if (breached[bnode]) continue;
   371             breached[bnode] = true;
   372             if (bnode_matching[bnode] != INVALID) {
   373               Node newanode = graph->aNode(bnode_matching[bnode]);
   374               if (!areached[newanode]) {
   375                 areached[newanode] = true;
   376                 newqueue.push_back(newanode);
   377               }
   378             }
   379           }
   380         }
   381         queue.swap(newqueue);
   382       }
   383 
   384       int size = 0;
   385       for (ANodeIt it(*graph); it != INVALID; ++it) {
   386         covering[it] = !areached[it] && anode_matching[it] != INVALID;
   387         if (!areached[it] && anode_matching[it] != INVALID) {
   388           ++size;
   389         }
   390       }
   391       for (BNodeIt it(*graph); it != INVALID; ++it) {
   392         covering[it] = breached[it];
   393         if (breached[it]) {
   394           ++size;
   395         }
   396       }
   397       return size;
   398     }
   399 
   400     /// \brief Set true all matching uedge in the map.
   401     /// 
   402     /// Set true all matching uedge in the map. It does not change the
   403     /// value mapped to the other uedges.
   404     /// \return The number of the matching edges.
   405     template <typename MatchingMap>
   406     int quickMatching(MatchingMap& mm) const {
   407       for (ANodeIt it(*graph); it != INVALID; ++it) {
   408         if (anode_matching[it] != INVALID) {
   409           mm[anode_matching[it]] = true;
   410         }
   411       }
   412       return matching_size;
   413     }
   414 
   415     /// \brief Set true all matching uedge in the map and the others to false.
   416     /// 
   417     /// Set true all matching uedge in the map and the others to false.
   418     /// \return The number of the matching edges.
   419     template <typename MatchingMap>
   420     int matching(MatchingMap& mm) const {
   421       for (UEdgeIt it(*graph); it != INVALID; ++it) {
   422         mm[it] = it == anode_matching[graph->aNode(it)];
   423       }
   424       return matching_size;
   425     }
   426 
   427 
   428     /// \brief Return true if the given uedge is in the matching.
   429     /// 
   430     /// It returns true if the given uedge is in the matching.
   431     bool matchingEdge(const UEdge& edge) const {
   432       return anode_matching[graph->aNode(edge)] == edge;
   433     }
   434 
   435     /// \brief Returns the matching edge from the node.
   436     /// 
   437     /// Returns the matching edge from the node. If there is not such
   438     /// edge it gives back \c INVALID.
   439     UEdge matchingEdge(const Node& node) const {
   440       if (graph->aNode(node)) {
   441         return anode_matching[node];
   442       } else {
   443         return bnode_matching[node];
   444       }
   445     }
   446 
   447     /// \brief Gives back the number of the matching edges.
   448     ///
   449     /// Gives back the number of the matching edges.
   450     int matchingSize() const {
   451       return matching_size;
   452     }
   453 
   454     /// @}
   455 
   456   private:
   457 
   458     ANodeMatchingMap anode_matching;
   459     BNodeMatchingMap bnode_matching;
   460     const Graph *graph;
   461 
   462     int matching_size;
   463   
   464   };
   465 
   466   /// \ingroup matching
   467   ///
   468   /// \brief Maximum cardinality bipartite matching
   469   ///
   470   /// This function calculates the maximum cardinality matching
   471   /// in a bipartite graph. It gives back the matching in an undirected
   472   /// edge map.
   473   ///
   474   /// \param graph The bipartite graph.
   475   /// \retval matching The undirected edge map which will be set to 
   476   /// the matching.
   477   /// \return The size of the matching.
   478   template <typename BpUGraph, typename MatchingMap>
   479   int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {
   480     MaxBipartiteMatching<BpUGraph> bpmatching(graph);
   481     bpmatching.run();
   482     bpmatching.matching(matching);
   483     return bpmatching.matchingSize();
   484   }
   485 
   486   /// \brief Default traits class for weighted bipartite matching algoritms.
   487   ///
   488   /// Default traits class for weighted bipartite matching algoritms.
   489   /// \param _BpUGraph The bipartite undirected graph type.
   490   /// \param _WeightMap Type of weight map.
   491   template <typename _BpUGraph, typename _WeightMap>
   492   struct WeightedBipartiteMatchingDefaultTraits {
   493     /// \brief The type of the weight of the undirected edges.
   494     typedef typename _WeightMap::Value Value;
   495 
   496     /// The undirected bipartite graph type the algorithm runs on. 
   497     typedef _BpUGraph BpUGraph;
   498 
   499     /// The map of the edges weights
   500     typedef _WeightMap WeightMap;
   501 
   502     /// \brief The cross reference type used by heap.
   503     ///
   504     /// The cross reference type used by heap.
   505     /// Usually it is \c Graph::NodeMap<int>.
   506     typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
   507 
   508     /// \brief Instantiates a HeapCrossRef.
   509     ///
   510     /// This function instantiates a \ref HeapCrossRef. 
   511     /// \param graph is the graph, to which we would like to define the 
   512     /// HeapCrossRef.
   513     static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
   514       return new HeapCrossRef(graph);
   515     }
   516     
   517     /// \brief The heap type used by weighted matching algorithms.
   518     ///
   519     /// The heap type used by weighted matching algorithms. It should
   520     /// minimize the priorities and the heap's key type is the graph's
   521     /// anode graph's node.
   522     ///
   523     /// \sa BinHeap
   524     typedef BinHeap<Value, HeapCrossRef> Heap;
   525     
   526     /// \brief Instantiates a Heap.
   527     ///
   528     /// This function instantiates a \ref Heap. 
   529     /// \param crossref The cross reference of the heap.
   530     static Heap *createHeap(HeapCrossRef& crossref) {
   531       return new Heap(crossref);
   532     }
   533 
   534   };
   535 
   536 
   537   /// \ingroup matching
   538   ///
   539   /// \brief Bipartite Max Weighted Matching algorithm
   540   ///
   541   /// This class implements the bipartite Max Weighted Matching
   542   /// algorithm.  It uses the successive shortest path algorithm to
   543   /// calculate the maximum weighted matching in the bipartite
   544   /// graph. The algorithm can be used also to calculate the maximum
   545   /// cardinality maximum weighted matching. The time complexity
   546   /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
   547   /// heap implementation but this can be improved to 
   548   /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
   549   ///
   550   /// The algorithm also provides a potential function on the nodes
   551   /// which a dual solution of the matching algorithm and it can be
   552   /// used to proof the optimality of the given pimal solution.
   553 #ifdef DOXYGEN
   554   template <typename _BpUGraph, typename _WeightMap, typename _Traits>
   555 #else
   556   template <typename _BpUGraph, 
   557             typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
   558             typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
   559 #endif
   560   class MaxWeightedBipartiteMatching {
   561   public:
   562 
   563     typedef _Traits Traits;
   564     typedef typename Traits::BpUGraph BpUGraph;
   565     typedef typename Traits::WeightMap WeightMap;
   566     typedef typename Traits::Value Value;
   567 
   568   protected:
   569 
   570     typedef typename Traits::HeapCrossRef HeapCrossRef;
   571     typedef typename Traits::Heap Heap; 
   572 
   573     
   574     typedef typename BpUGraph::Node Node;
   575     typedef typename BpUGraph::ANodeIt ANodeIt;
   576     typedef typename BpUGraph::BNodeIt BNodeIt;
   577     typedef typename BpUGraph::UEdge UEdge;
   578     typedef typename BpUGraph::UEdgeIt UEdgeIt;
   579     typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
   580 
   581     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
   582     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
   583 
   584     typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
   585     typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
   586 
   587 
   588   public:
   589 
   590     /// \brief \ref Exception for uninitialized parameters.
   591     ///
   592     /// This error represents problems in the initialization
   593     /// of the parameters of the algorithms.
   594     class UninitializedParameter : public lemon::UninitializedParameter {
   595     public:
   596       virtual const char* what() const throw() {
   597 	return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
   598       }
   599     };
   600 
   601     ///\name Named template parameters
   602 
   603     ///@{
   604 
   605     template <class H, class CR>
   606     struct DefHeapTraits : public Traits {
   607       typedef CR HeapCrossRef;
   608       typedef H Heap;
   609       static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
   610 	throw UninitializedParameter();
   611       }
   612       static Heap *createHeap(HeapCrossRef &) {
   613 	throw UninitializedParameter();
   614       }
   615     };
   616 
   617     /// \brief \ref named-templ-param "Named parameter" for setting heap 
   618     /// and cross reference type
   619     ///
   620     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   621     /// reference type
   622     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
   623     struct DefHeap
   624       : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
   625                                             DefHeapTraits<H, CR> > { 
   626       typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
   627                                            DefHeapTraits<H, CR> > Create;
   628     };
   629 
   630     template <class H, class CR>
   631     struct DefStandardHeapTraits : public Traits {
   632       typedef CR HeapCrossRef;
   633       typedef H Heap;
   634       static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
   635 	return new HeapCrossRef(graph);
   636       }
   637       static Heap *createHeap(HeapCrossRef &crossref) {
   638 	return new Heap(crossref);
   639       }
   640     };
   641 
   642     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
   643     /// cross reference type with automatic allocation
   644     ///
   645     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   646     /// reference type. It can allocate the heap and the cross reference 
   647     /// object if the cross reference's constructor waits for the graph as 
   648     /// parameter and the heap's constructor waits for the cross reference.
   649     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
   650     struct DefStandardHeap
   651       : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
   652                                             DefStandardHeapTraits<H, CR> > { 
   653       typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
   654                                            DefStandardHeapTraits<H, CR> > 
   655       Create;
   656     };
   657 
   658     ///@}
   659 
   660 
   661     /// \brief Constructor.
   662     ///
   663     /// Constructor of the algorithm. 
   664     MaxWeightedBipartiteMatching(const BpUGraph& _graph, 
   665                                  const WeightMap& _weight) 
   666       : graph(&_graph), weight(&_weight),
   667         anode_matching(_graph), bnode_matching(_graph),
   668         anode_potential(_graph), bnode_potential(_graph),
   669         _heap_cross_ref(0), local_heap_cross_ref(false),
   670         _heap(0), local_heap(0) {}
   671 
   672     /// \brief Destructor.
   673     ///
   674     /// Destructor of the algorithm.
   675     ~MaxWeightedBipartiteMatching() {
   676       destroyStructures();
   677     }
   678 
   679     /// \brief Sets the heap and the cross reference used by algorithm.
   680     ///
   681     /// Sets the heap and the cross reference used by algorithm.
   682     /// If you don't use this function before calling \ref run(),
   683     /// it will allocate one. The destuctor deallocates this
   684     /// automatically allocated map, of course.
   685     /// \return \c (*this)
   686     MaxWeightedBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
   687       if(local_heap_cross_ref) {
   688 	delete _heap_cross_ref;
   689 	local_heap_cross_ref = false;
   690       }
   691       _heap_cross_ref = &cr;
   692       if(local_heap) {
   693 	delete _heap;
   694 	local_heap = false;
   695       }
   696       _heap = &hp;
   697       return *this;
   698     }
   699 
   700     /// \name Execution control
   701     /// The simplest way to execute the algorithm is to use
   702     /// one of the member functions called \c run().
   703     /// \n
   704     /// If you need more control on the execution,
   705     /// first you must call \ref init() or one alternative for it.
   706     /// Finally \ref start() will perform the matching computation or
   707     /// with step-by-step execution you can augment the solution.
   708 
   709     /// @{
   710 
   711     /// \brief Initalize the data structures.
   712     ///
   713     /// It initalizes the data structures and creates an empty matching.
   714     void init() {
   715       initStructures();
   716       for (ANodeIt it(*graph); it != INVALID; ++it) {
   717         anode_matching[it] = INVALID;
   718         anode_potential[it] = 0;
   719       }
   720       for (BNodeIt it(*graph); it != INVALID; ++it) {
   721         bnode_matching[it] = INVALID;
   722         bnode_potential[it] = 0;
   723         for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
   724           if ((*weight)[jt] > bnode_potential[it]) {
   725             bnode_potential[it] = (*weight)[jt];
   726           }
   727         }
   728       }
   729       matching_value = 0;
   730       matching_size = 0;
   731     }
   732 
   733 
   734     /// \brief An augmenting phase of the weighted matching algorithm
   735     ///
   736     /// It runs an augmenting phase of the weighted matching 
   737     /// algorithm. This phase finds the best augmenting path and 
   738     /// augments only on this paths. 
   739     ///
   740     /// The algorithm consists at most 
   741     /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
   742     /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
   743     /// with binary heap.
   744     /// \param decrease If the given parameter true the matching value
   745     /// can be decreased in the augmenting phase. If we would like
   746     /// to calculate the maximum cardinality maximum weighted matching
   747     /// then we should let the algorithm to decrease the matching
   748     /// value in order to increase the number of the matching edges.
   749     bool augment(bool decrease = false) {
   750 
   751       typename BpUGraph::template BNodeMap<Value> bdist(*graph);
   752       typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
   753 
   754       Node bestNode = INVALID;
   755       Value bestValue = 0;
   756 
   757       _heap->clear();
   758       for (ANodeIt it(*graph); it != INVALID; ++it) {
   759         (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
   760       }
   761 
   762       for (ANodeIt it(*graph); it != INVALID; ++it) {
   763         if (anode_matching[it] == INVALID) {
   764           _heap->push(it, 0);
   765         }
   766       }
   767 
   768       Value bdistMax = 0;
   769       while (!_heap->empty()) {
   770         Node anode = _heap->top();
   771         Value avalue = _heap->prio();
   772         _heap->pop();
   773         for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   774           if (jt == anode_matching[anode]) continue;
   775           Node bnode = graph->bNode(jt);
   776           Value bvalue = avalue  - (*weight)[jt] +
   777             anode_potential[anode] + bnode_potential[bnode];
   778           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
   779             bdist[bnode] = bvalue;
   780             bpred[bnode] = jt;
   781           }
   782           if (bvalue > bdistMax) {
   783             bdistMax = bvalue;
   784           }
   785           if (bnode_matching[bnode] != INVALID) {
   786             Node newanode = graph->aNode(bnode_matching[bnode]);
   787             switch (_heap->state(newanode)) {
   788             case Heap::PRE_HEAP:
   789               _heap->push(newanode, bvalue);
   790               break;
   791             case Heap::IN_HEAP:
   792               if (bvalue < (*_heap)[newanode]) {
   793                 _heap->decrease(newanode, bvalue);
   794               }
   795               break;
   796             case Heap::POST_HEAP:
   797               break;
   798             }
   799           } else {
   800             if (bestNode == INVALID || 
   801                 bnode_potential[bnode] - bvalue > bestValue) {
   802               bestValue = bnode_potential[bnode] - bvalue;
   803               bestNode = bnode;
   804             }
   805           }
   806         }
   807       }
   808 
   809       if (bestNode == INVALID || (!decrease && bestValue < 0)) {
   810         return false;
   811       }
   812 
   813       matching_value += bestValue;
   814       ++matching_size;
   815 
   816       for (BNodeIt it(*graph); it != INVALID; ++it) {
   817         if (bpred[it] != INVALID) {
   818           bnode_potential[it] -= bdist[it];
   819         } else {
   820           bnode_potential[it] -= bdistMax;
   821         }
   822       }
   823       for (ANodeIt it(*graph); it != INVALID; ++it) {
   824         if (anode_matching[it] != INVALID) {
   825           Node bnode = graph->bNode(anode_matching[it]);
   826           if (bpred[bnode] != INVALID) {
   827             anode_potential[it] += bdist[bnode];
   828           } else {
   829             anode_potential[it] += bdistMax;
   830           }
   831         }
   832       }
   833 
   834       while (bestNode != INVALID) {
   835         UEdge uedge = bpred[bestNode];
   836         Node anode = graph->aNode(uedge);
   837         
   838         bnode_matching[bestNode] = uedge;
   839         if (anode_matching[anode] != INVALID) {
   840           bestNode = graph->bNode(anode_matching[anode]);
   841         } else {
   842           bestNode = INVALID;
   843         }
   844         anode_matching[anode] = uedge;
   845       }
   846 
   847 
   848       return true;
   849     }
   850 
   851     /// \brief Starts the algorithm.
   852     ///
   853     /// Starts the algorithm. It runs augmenting phases until the
   854     /// optimal solution reached.
   855     ///
   856     /// \param maxCardinality If the given value is true it will
   857     /// calculate the maximum cardinality maximum matching instead of
   858     /// the maximum matching.
   859     void start(bool maxCardinality = false) {
   860       while (augment(maxCardinality)) {}
   861     }
   862 
   863     /// \brief Runs the algorithm.
   864     ///
   865     /// It just initalize the algorithm and then start it.
   866     ///
   867     /// \param maxCardinality If the given value is true it will
   868     /// calculate the maximum cardinality maximum matching instead of
   869     /// the maximum matching.
   870     void run(bool maxCardinality = false) {
   871       init();
   872       start(maxCardinality);
   873     }
   874 
   875     /// @}
   876 
   877     /// \name Query Functions
   878     /// The result of the %Matching algorithm can be obtained using these
   879     /// functions.\n
   880     /// Before the use of these functions,
   881     /// either run() or start() must be called.
   882     
   883     ///@{
   884 
   885     /// \brief Gives back the potential in the NodeMap
   886     ///
   887     /// Gives back the potential in the NodeMap. The matching is optimal
   888     /// with the current number of edges if \f$ \pi(a) + \pi(b) - w(ab) = 0 \f$
   889     /// for each matching edges and \f$ \pi(a) + \pi(b) - w(ab) \ge 0 \f$
   890     /// for each edges. 
   891     template <typename PotentialMap>
   892     void potential(PotentialMap& pt) const {
   893       for (ANodeIt it(*graph); it != INVALID; ++it) {
   894         pt[it] = anode_potential[it];
   895       }
   896       for (BNodeIt it(*graph); it != INVALID; ++it) {
   897         pt[it] = bnode_potential[it];
   898       }
   899     }
   900 
   901     /// \brief Set true all matching uedge in the map.
   902     /// 
   903     /// Set true all matching uedge in the map. It does not change the
   904     /// value mapped to the other uedges.
   905     /// \return The number of the matching edges.
   906     template <typename MatchingMap>
   907     int quickMatching(MatchingMap& mm) const {
   908       for (ANodeIt it(*graph); it != INVALID; ++it) {
   909         if (anode_matching[it] != INVALID) {
   910           mm[anode_matching[it]] = true;
   911         }
   912       }
   913       return matching_size;
   914     }
   915 
   916     /// \brief Set true all matching uedge in the map and the others to false.
   917     /// 
   918     /// Set true all matching uedge in the map and the others to false.
   919     /// \return The number of the matching edges.
   920     template <typename MatchingMap>
   921     int matching(MatchingMap& mm) const {
   922       for (UEdgeIt it(*graph); it != INVALID; ++it) {
   923         mm[it] = it == anode_matching[graph->aNode(it)];
   924       }
   925       return matching_size;
   926     }
   927 
   928 
   929     /// \brief Return true if the given uedge is in the matching.
   930     /// 
   931     /// It returns true if the given uedge is in the matching.
   932     bool matchingEdge(const UEdge& edge) const {
   933       return anode_matching[graph->aNode(edge)] == edge;
   934     }
   935 
   936     /// \brief Returns the matching edge from the node.
   937     /// 
   938     /// Returns the matching edge from the node. If there is not such
   939     /// edge it gives back \c INVALID.
   940     UEdge matchingEdge(const Node& node) const {
   941       if (graph->aNode(node)) {
   942         return anode_matching[node];
   943       } else {
   944         return bnode_matching[node];
   945       }
   946     }
   947 
   948     /// \brief Gives back the sum of weights of the matching edges.
   949     ///
   950     /// Gives back the sum of weights of the matching edges.
   951     Value matchingValue() const {
   952       return matching_value;
   953     }
   954 
   955     /// \brief Gives back the number of the matching edges.
   956     ///
   957     /// Gives back the number of the matching edges.
   958     int matchingSize() const {
   959       return matching_size;
   960     }
   961 
   962     /// @}
   963 
   964   private:
   965 
   966     void initStructures() {
   967       if (!_heap_cross_ref) {
   968 	local_heap_cross_ref = true;
   969 	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
   970       }
   971       if (!_heap) {
   972 	local_heap = true;
   973 	_heap = Traits::createHeap(*_heap_cross_ref);
   974       }
   975     }
   976 
   977     void destroyStructures() {
   978       if (local_heap_cross_ref) delete _heap_cross_ref;
   979       if (local_heap) delete _heap;
   980     }
   981 
   982 
   983   private:
   984     
   985     const BpUGraph *graph;
   986     const WeightMap* weight;
   987 
   988     ANodeMatchingMap anode_matching;
   989     BNodeMatchingMap bnode_matching;
   990 
   991     ANodePotentialMap anode_potential;
   992     BNodePotentialMap bnode_potential;
   993 
   994     Value matching_value;
   995     int matching_size;
   996 
   997     HeapCrossRef *_heap_cross_ref;
   998     bool local_heap_cross_ref;
   999 
  1000     Heap *_heap;
  1001     bool local_heap;
  1002   
  1003   };
  1004 
  1005   /// \ingroup matching
  1006   ///
  1007   /// \brief Maximum weighted bipartite matching
  1008   ///
  1009   /// This function calculates the maximum weighted matching
  1010   /// in a bipartite graph. It gives back the matching in an undirected
  1011   /// edge map.
  1012   ///
  1013   /// \param graph The bipartite graph.
  1014   /// \param weight The undirected edge map which contains the weights.
  1015   /// \retval matching The undirected edge map which will be set to 
  1016   /// the matching.
  1017   /// \return The value of the matching.
  1018   template <typename BpUGraph, typename WeightMap, typename MatchingMap>
  1019   typename WeightMap::Value 
  1020   maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,
  1021                                MatchingMap& matching) {
  1022     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
  1023       bpmatching(graph, weight);
  1024     bpmatching.run();
  1025     bpmatching.matching(matching);
  1026     return bpmatching.matchingValue();
  1027   }
  1028 
  1029   /// \ingroup matching
  1030   ///
  1031   /// \brief Maximum weighted maximum cardinality bipartite matching
  1032   ///
  1033   /// This function calculates the maximum weighted of the maximum cardinality
  1034   /// matchings of a bipartite graph. It gives back the matching in an 
  1035   /// undirected edge map.
  1036   ///
  1037   /// \param graph The bipartite graph.
  1038   /// \param weight The undirected edge map which contains the weights.
  1039   /// \retval matching The undirected edge map which will be set to 
  1040   /// the matching.
  1041   /// \return The value of the matching.
  1042   template <typename BpUGraph, typename WeightMap, typename MatchingMap>
  1043   typename WeightMap::Value 
  1044   maxWeightedMaxBipartiteMatching(const BpUGraph& graph, 
  1045                                   const WeightMap& weight,
  1046                                   MatchingMap& matching) {
  1047     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 
  1048       bpmatching(graph, weight);
  1049     bpmatching.run(true);
  1050     bpmatching.matching(matching);
  1051     return bpmatching.matchingValue();
  1052   }
  1053 
  1054   /// \brief Default traits class for minimum cost bipartite matching
  1055   /// algoritms.
  1056   ///
  1057   /// Default traits class for minimum cost bipartite matching
  1058   /// algoritms.  
  1059   ///
  1060   /// \param _BpUGraph The bipartite undirected graph
  1061   /// type.  
  1062   ///
  1063   /// \param _CostMap Type of cost map.
  1064   template <typename _BpUGraph, typename _CostMap>
  1065   struct MinCostMaxBipartiteMatchingDefaultTraits {
  1066     /// \brief The type of the cost of the undirected edges.
  1067     typedef typename _CostMap::Value Value;
  1068 
  1069     /// The undirected bipartite graph type the algorithm runs on. 
  1070     typedef _BpUGraph BpUGraph;
  1071 
  1072     /// The map of the edges costs
  1073     typedef _CostMap CostMap;
  1074 
  1075     /// \brief The cross reference type used by heap.
  1076     ///
  1077     /// The cross reference type used by heap.
  1078     /// Usually it is \c Graph::NodeMap<int>.
  1079     typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
  1080 
  1081     /// \brief Instantiates a HeapCrossRef.
  1082     ///
  1083     /// This function instantiates a \ref HeapCrossRef. 
  1084     /// \param graph is the graph, to which we would like to define the 
  1085     /// HeapCrossRef.
  1086     static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
  1087       return new HeapCrossRef(graph);
  1088     }
  1089     
  1090     /// \brief The heap type used by costed matching algorithms.
  1091     ///
  1092     /// The heap type used by costed matching algorithms. It should
  1093     /// minimize the priorities and the heap's key type is the graph's
  1094     /// anode graph's node.
  1095     ///
  1096     /// \sa BinHeap
  1097     typedef BinHeap<Value, HeapCrossRef> Heap;
  1098     
  1099     /// \brief Instantiates a Heap.
  1100     ///
  1101     /// This function instantiates a \ref Heap. 
  1102     /// \param crossref The cross reference of the heap.
  1103     static Heap *createHeap(HeapCrossRef& crossref) {
  1104       return new Heap(crossref);
  1105     }
  1106 
  1107   };
  1108 
  1109 
  1110   /// \ingroup matching
  1111   ///
  1112   /// \brief Bipartite Min Cost Matching algorithm
  1113   ///
  1114   /// This class implements the bipartite Min Cost Matching algorithm.
  1115   /// It uses the successive shortest path algorithm to calculate the
  1116   /// minimum cost maximum matching in the bipartite graph. The time
  1117   /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
  1118   /// default binary heap implementation but this can be improved to
  1119   /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
  1120   ///
  1121   /// The algorithm also provides a potential function on the nodes
  1122   /// which a dual solution of the matching algorithm and it can be
  1123   /// used to proof the optimality of the given pimal solution.
  1124 #ifdef DOXYGEN
  1125   template <typename _BpUGraph, typename _CostMap, typename _Traits>
  1126 #else
  1127   template <typename _BpUGraph, 
  1128             typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
  1129             typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
  1130 #endif
  1131   class MinCostMaxBipartiteMatching {
  1132   public:
  1133 
  1134     typedef _Traits Traits;
  1135     typedef typename Traits::BpUGraph BpUGraph;
  1136     typedef typename Traits::CostMap CostMap;
  1137     typedef typename Traits::Value Value;
  1138 
  1139   protected:
  1140 
  1141     typedef typename Traits::HeapCrossRef HeapCrossRef;
  1142     typedef typename Traits::Heap Heap; 
  1143 
  1144     
  1145     typedef typename BpUGraph::Node Node;
  1146     typedef typename BpUGraph::ANodeIt ANodeIt;
  1147     typedef typename BpUGraph::BNodeIt BNodeIt;
  1148     typedef typename BpUGraph::UEdge UEdge;
  1149     typedef typename BpUGraph::UEdgeIt UEdgeIt;
  1150     typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
  1151 
  1152     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
  1153     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
  1154 
  1155     typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
  1156     typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
  1157 
  1158 
  1159   public:
  1160 
  1161     /// \brief \ref Exception for uninitialized parameters.
  1162     ///
  1163     /// This error represents problems in the initialization
  1164     /// of the parameters of the algorithms.
  1165     class UninitializedParameter : public lemon::UninitializedParameter {
  1166     public:
  1167       virtual const char* what() const throw() {
  1168 	return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
  1169       }
  1170     };
  1171 
  1172     ///\name Named template parameters
  1173 
  1174     ///@{
  1175 
  1176     template <class H, class CR>
  1177     struct DefHeapTraits : public Traits {
  1178       typedef CR HeapCrossRef;
  1179       typedef H Heap;
  1180       static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
  1181 	throw UninitializedParameter();
  1182       }
  1183       static Heap *createHeap(HeapCrossRef &) {
  1184 	throw UninitializedParameter();
  1185       }
  1186     };
  1187 
  1188     /// \brief \ref named-templ-param "Named parameter" for setting heap 
  1189     /// and cross reference type
  1190     ///
  1191     /// \ref named-templ-param "Named parameter" for setting heap and cross 
  1192     /// reference type
  1193     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
  1194     struct DefHeap
  1195       : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
  1196                                             DefHeapTraits<H, CR> > { 
  1197       typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
  1198                                            DefHeapTraits<H, CR> > Create;
  1199     };
  1200 
  1201     template <class H, class CR>
  1202     struct DefStandardHeapTraits : public Traits {
  1203       typedef CR HeapCrossRef;
  1204       typedef H Heap;
  1205       static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
  1206 	return new HeapCrossRef(graph);
  1207       }
  1208       static Heap *createHeap(HeapCrossRef &crossref) {
  1209 	return new Heap(crossref);
  1210       }
  1211     };
  1212 
  1213     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
  1214     /// cross reference type with automatic allocation
  1215     ///
  1216     /// \ref named-templ-param "Named parameter" for setting heap and cross 
  1217     /// reference type. It can allocate the heap and the cross reference 
  1218     /// object if the cross reference's constructor waits for the graph as 
  1219     /// parameter and the heap's constructor waits for the cross reference.
  1220     template <class H, class CR = typename BpUGraph::template NodeMap<int> >
  1221     struct DefStandardHeap
  1222       : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
  1223                                             DefStandardHeapTraits<H, CR> > { 
  1224       typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
  1225                                            DefStandardHeapTraits<H, CR> > 
  1226       Create;
  1227     };
  1228 
  1229     ///@}
  1230 
  1231 
  1232     /// \brief Constructor.
  1233     ///
  1234     /// Constructor of the algorithm. 
  1235     MinCostMaxBipartiteMatching(const BpUGraph& _graph, 
  1236                                  const CostMap& _cost) 
  1237       : graph(&_graph), cost(&_cost),
  1238         anode_matching(_graph), bnode_matching(_graph),
  1239         anode_potential(_graph), bnode_potential(_graph),
  1240         _heap_cross_ref(0), local_heap_cross_ref(false),
  1241         _heap(0), local_heap(0) {}
  1242 
  1243     /// \brief Destructor.
  1244     ///
  1245     /// Destructor of the algorithm.
  1246     ~MinCostMaxBipartiteMatching() {
  1247       destroyStructures();
  1248     }
  1249 
  1250     /// \brief Sets the heap and the cross reference used by algorithm.
  1251     ///
  1252     /// Sets the heap and the cross reference used by algorithm.
  1253     /// If you don't use this function before calling \ref run(),
  1254     /// it will allocate one. The destuctor deallocates this
  1255     /// automatically allocated map, of course.
  1256     /// \return \c (*this)
  1257     MinCostMaxBipartiteMatching& heap(Heap& hp, HeapCrossRef &cr) {
  1258       if(local_heap_cross_ref) {
  1259 	delete _heap_cross_ref;
  1260 	local_heap_cross_ref = false;
  1261       }
  1262       _heap_cross_ref = &cr;
  1263       if(local_heap) {
  1264 	delete _heap;
  1265 	local_heap = false;
  1266       }
  1267       _heap = &hp;
  1268       return *this;
  1269     }
  1270 
  1271     /// \name Execution control
  1272     /// The simplest way to execute the algorithm is to use
  1273     /// one of the member functions called \c run().
  1274     /// \n
  1275     /// If you need more control on the execution,
  1276     /// first you must call \ref init() or one alternative for it.
  1277     /// Finally \ref start() will perform the matching computation or
  1278     /// with step-by-step execution you can augment the solution.
  1279 
  1280     /// @{
  1281 
  1282     /// \brief Initalize the data structures.
  1283     ///
  1284     /// It initalizes the data structures and creates an empty matching.
  1285     void init() {
  1286       initStructures();
  1287       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1288         anode_matching[it] = INVALID;
  1289         anode_potential[it] = 0;
  1290       }
  1291       for (BNodeIt it(*graph); it != INVALID; ++it) {
  1292         bnode_matching[it] = INVALID;
  1293         bnode_potential[it] = 0;
  1294       }
  1295       matching_cost = 0;
  1296       matching_size = 0;
  1297     }
  1298 
  1299 
  1300     /// \brief An augmenting phase of the costed matching algorithm
  1301     ///
  1302     /// It runs an augmenting phase of the matching algorithm. The
  1303     /// phase finds the best augmenting path and augments only on this
  1304     /// paths.
  1305     ///
  1306     /// The algorithm consists at most 
  1307     /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
  1308     /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
  1309     /// with binary heap.
  1310     bool augment() {
  1311 
  1312       typename BpUGraph::template BNodeMap<Value> bdist(*graph);
  1313       typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
  1314 
  1315       Node bestNode = INVALID;
  1316       Value bestValue = 0;
  1317 
  1318       _heap->clear();
  1319       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1320         (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
  1321       }
  1322 
  1323       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1324         if (anode_matching[it] == INVALID) {
  1325           _heap->push(it, 0);
  1326         }
  1327       }
  1328       Value bdistMax = 0;
  1329 
  1330       while (!_heap->empty()) {
  1331         Node anode = _heap->top();
  1332         Value avalue = _heap->prio();
  1333         _heap->pop();
  1334         for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
  1335           if (jt == anode_matching[anode]) continue;
  1336           Node bnode = graph->bNode(jt);
  1337           Value bvalue = avalue + (*cost)[jt] + 
  1338             anode_potential[anode] - bnode_potential[bnode];
  1339           if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
  1340             bdist[bnode] = bvalue;
  1341             bpred[bnode] = jt;
  1342           }
  1343           if (bvalue > bdistMax) {
  1344             bdistMax = bvalue;
  1345           }
  1346           if (bnode_matching[bnode] != INVALID) {
  1347             Node newanode = graph->aNode(bnode_matching[bnode]);
  1348             switch (_heap->state(newanode)) {
  1349             case Heap::PRE_HEAP:
  1350               _heap->push(newanode, bvalue);
  1351               break;
  1352             case Heap::IN_HEAP:
  1353               if (bvalue < (*_heap)[newanode]) {
  1354                 _heap->decrease(newanode, bvalue);
  1355               }
  1356               break;
  1357             case Heap::POST_HEAP:
  1358               break;
  1359             }
  1360           } else {
  1361             if (bestNode == INVALID || 
  1362                 bvalue + bnode_potential[bnode] < bestValue) {
  1363               bestValue = bvalue + bnode_potential[bnode];
  1364               bestNode = bnode;
  1365             }
  1366           }
  1367         }
  1368       }
  1369 
  1370       if (bestNode == INVALID) {
  1371         return false;
  1372       }
  1373 
  1374       matching_cost += bestValue;
  1375       ++matching_size;
  1376 
  1377       for (BNodeIt it(*graph); it != INVALID; ++it) {
  1378         if (bpred[it] != INVALID) {
  1379           bnode_potential[it] += bdist[it];
  1380         } else {
  1381           bnode_potential[it] += bdistMax;
  1382         }
  1383       }
  1384       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1385         if (anode_matching[it] != INVALID) {
  1386           Node bnode = graph->bNode(anode_matching[it]);
  1387           if (bpred[bnode] != INVALID) {
  1388             anode_potential[it] += bdist[bnode];
  1389           } else {
  1390             anode_potential[it] += bdistMax;
  1391           }
  1392         }
  1393       }
  1394 
  1395       while (bestNode != INVALID) {
  1396         UEdge uedge = bpred[bestNode];
  1397         Node anode = graph->aNode(uedge);
  1398         
  1399         bnode_matching[bestNode] = uedge;
  1400         if (anode_matching[anode] != INVALID) {
  1401           bestNode = graph->bNode(anode_matching[anode]);
  1402         } else {
  1403           bestNode = INVALID;
  1404         }
  1405         anode_matching[anode] = uedge;
  1406       }
  1407 
  1408 
  1409       return true;
  1410     }
  1411 
  1412     /// \brief Starts the algorithm.
  1413     ///
  1414     /// Starts the algorithm. It runs augmenting phases until the
  1415     /// optimal solution reached.
  1416     void start() {
  1417       while (augment()) {}
  1418     }
  1419 
  1420     /// \brief Runs the algorithm.
  1421     ///
  1422     /// It just initalize the algorithm and then start it.
  1423     void run() {
  1424       init();
  1425       start();
  1426     }
  1427 
  1428     /// @}
  1429 
  1430     /// \name Query Functions
  1431     /// The result of the %Matching algorithm can be obtained using these
  1432     /// functions.\n
  1433     /// Before the use of these functions,
  1434     /// either run() or start() must be called.
  1435     
  1436     ///@{
  1437 
  1438     /// \brief Gives back the potential in the NodeMap
  1439     ///
  1440     /// Gives back the potential in the NodeMap. The potential is optimal with 
  1441     /// the current number of edges if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
  1442     /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
  1443     /// for each edges.
  1444     template <typename PotentialMap>
  1445     void potential(PotentialMap& pt) const {
  1446       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1447         pt[it] = anode_potential[it];
  1448       }
  1449       for (BNodeIt it(*graph); it != INVALID; ++it) {
  1450         pt[it] = bnode_potential[it];
  1451       }
  1452     }
  1453 
  1454     /// \brief Set true all matching uedge in the map.
  1455     /// 
  1456     /// Set true all matching uedge in the map. It does not change the
  1457     /// value mapped to the other uedges.
  1458     /// \return The number of the matching edges.
  1459     template <typename MatchingMap>
  1460     int quickMatching(MatchingMap& mm) const {
  1461       for (ANodeIt it(*graph); it != INVALID; ++it) {
  1462         if (anode_matching[it] != INVALID) {
  1463           mm[anode_matching[it]] = true;
  1464         }
  1465       }
  1466       return matching_size;
  1467     }
  1468 
  1469     /// \brief Set true all matching uedge in the map and the others to false.
  1470     /// 
  1471     /// Set true all matching uedge in the map and the others to false.
  1472     /// \return The number of the matching edges.
  1473     template <typename MatchingMap>
  1474     int matching(MatchingMap& mm) const {
  1475       for (UEdgeIt it(*graph); it != INVALID; ++it) {
  1476         mm[it] = it == anode_matching[graph->aNode(it)];
  1477       }
  1478       return matching_size;
  1479     }
  1480 
  1481 
  1482     /// \brief Return true if the given uedge is in the matching.
  1483     /// 
  1484     /// It returns true if the given uedge is in the matching.
  1485     bool matchingEdge(const UEdge& edge) const {
  1486       return anode_matching[graph->aNode(edge)] == edge;
  1487     }
  1488 
  1489     /// \brief Returns the matching edge from the node.
  1490     /// 
  1491     /// Returns the matching edge from the node. If there is not such
  1492     /// edge it gives back \c INVALID.
  1493     UEdge matchingEdge(const Node& node) const {
  1494       if (graph->aNode(node)) {
  1495         return anode_matching[node];
  1496       } else {
  1497         return bnode_matching[node];
  1498       }
  1499     }
  1500 
  1501     /// \brief Gives back the sum of costs of the matching edges.
  1502     ///
  1503     /// Gives back the sum of costs of the matching edges.
  1504     Value matchingCost() const {
  1505       return matching_cost;
  1506     }
  1507 
  1508     /// \brief Gives back the number of the matching edges.
  1509     ///
  1510     /// Gives back the number of the matching edges.
  1511     int matchingSize() const {
  1512       return matching_size;
  1513     }
  1514 
  1515     /// @}
  1516 
  1517   private:
  1518 
  1519     void initStructures() {
  1520       if (!_heap_cross_ref) {
  1521 	local_heap_cross_ref = true;
  1522 	_heap_cross_ref = Traits::createHeapCrossRef(*graph);
  1523       }
  1524       if (!_heap) {
  1525 	local_heap = true;
  1526 	_heap = Traits::createHeap(*_heap_cross_ref);
  1527       }
  1528     }
  1529 
  1530     void destroyStructures() {
  1531       if (local_heap_cross_ref) delete _heap_cross_ref;
  1532       if (local_heap) delete _heap;
  1533     }
  1534 
  1535 
  1536   private:
  1537     
  1538     const BpUGraph *graph;
  1539     const CostMap* cost;
  1540 
  1541     ANodeMatchingMap anode_matching;
  1542     BNodeMatchingMap bnode_matching;
  1543 
  1544     ANodePotentialMap anode_potential;
  1545     BNodePotentialMap bnode_potential;
  1546 
  1547     Value matching_cost;
  1548     int matching_size;
  1549 
  1550     HeapCrossRef *_heap_cross_ref;
  1551     bool local_heap_cross_ref;
  1552 
  1553     Heap *_heap;
  1554     bool local_heap;
  1555   
  1556   };
  1557 
  1558   /// \ingroup matching
  1559   ///
  1560   /// \brief Minimum cost maximum cardinality bipartite matching
  1561   ///
  1562   /// This function calculates the minimum cost matching of the maximum 
  1563   /// cardinality matchings of a bipartite graph. It gives back the matching 
  1564   /// in an undirected edge map.
  1565   ///
  1566   /// \param graph The bipartite graph.
  1567   /// \param cost The undirected edge map which contains the costs.
  1568   /// \retval matching The undirected edge map which will be set to 
  1569   /// the matching.
  1570   /// \return The cost of the matching.
  1571   template <typename BpUGraph, typename CostMap, typename MatchingMap>
  1572   typename CostMap::Value 
  1573   minCostMaxBipartiteMatching(const BpUGraph& graph, 
  1574                               const CostMap& cost,
  1575                               MatchingMap& matching) {
  1576     MinCostMaxBipartiteMatching<BpUGraph, CostMap> 
  1577       bpmatching(graph, cost);
  1578     bpmatching.run();
  1579     bpmatching.matching(matching);
  1580     return bpmatching.matchingCost();
  1581   }
  1582 
  1583 }
  1584 
  1585 #endif