/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2006

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 #ifndef LEMON_BIPARTITE_MATCHING

 #define LEMON_BIPARTITE_MATCHING

 #include <functional>

 #include <lemon/bin_heap.h>

 #include <lemon/maps.h>

 #include <iostream>

 ///\ingroup matching

 ///\file

 ///\brief Maximum matching algorithms in bipartite graphs.

 namespace lemon {

   /// \ingroup matching

   ///

   /// \brief Bipartite Max Cardinality Matching algorithm

   ///

   /// Bipartite Max Cardinality Matching algorithm. This class implements

   /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time

   /// complexity.

   template <typename BpUGraph>

   class MaxBipartiteMatching {

   protected:

     typedef BpUGraph Graph;

     typedef typename Graph::Node Node;

     typedef typename Graph::ANodeIt ANodeIt;

     typedef typename Graph::BNodeIt BNodeIt;

     typedef typename Graph::UEdge UEdge;

     typedef typename Graph::UEdgeIt UEdgeIt;

     typedef typename Graph::IncEdgeIt IncEdgeIt;

     typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;

     typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;

   public:

     /// \brief Constructor.

     ///

     /// Constructor of the algorithm. 

     MaxBipartiteMatching(const BpUGraph& _graph) 

       : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}

     /// \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() {

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         anode_matching[it] = INVALID;

       }

       for (BNodeIt it(*graph); it != INVALID; ++it) {

         bnode_matching[it] = INVALID;

       }

       matching_size = 0;

     }

     /// \brief Initalize the data structures.

     ///

     /// It initalizes the data structures and creates a greedy

     /// matching.  From this matching sometimes it is faster to get

     /// the matching than from the initial empty matching.

     void greedyInit() {

       matching_size = 0;

       for (BNodeIt it(*graph); it != INVALID; ++it) {

         bnode_matching[it] = INVALID;

       }

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         anode_matching[it] = INVALID;

         for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {

           if (bnode_matching[graph->bNode(jt)] == INVALID) {

             anode_matching[it] = jt;

             bnode_matching[graph->bNode(jt)] = jt;

             ++matching_size;

             break;

           }

         }

       }

     }

     /// \brief Initalize the data structures with an initial matching.

     ///

     /// It initalizes the data structures with an initial matching.

     template <typename MatchingMap>

     void matchingInit(const MatchingMap& matching) {

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         anode_matching[it] = INVALID;

       }

       for (BNodeIt it(*graph); it != INVALID; ++it) {

         bnode_matching[it] = INVALID;

       }

       matching_size = 0;

       for (UEdgeIt it(*graph); it != INVALID; ++it) {

         if (matching[it]) {

           ++matching_size;

           anode_matching[graph->aNode(it)] = it;

           bnode_matching[graph->bNode(it)] = it;

         }

       }

     }

     /// \brief Initalize the data structures with an initial matching.

     ///

     /// It initalizes the data structures with an initial matching.

     /// \return %True when the given map contains really a matching.

     template <typename MatchingMap>

     void checkedMatchingInit(const MatchingMap& matching) {

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         anode_matching[it] = INVALID;

       }

       for (BNodeIt it(*graph); it != INVALID; ++it) {

         bnode_matching[it] = INVALID;

       }

       matching_size = 0;

       for (UEdgeIt it(*graph); it != INVALID; ++it) {

         if (matching[it]) {

           ++matching_size;

           if (anode_matching[graph->aNode(it)] != INVALID) {

             return false;

           }

           anode_matching[graph->aNode(it)] = it;

           if (bnode_matching[graph->aNode(it)] != INVALID) {

             return false;

           }

           bnode_matching[graph->bNode(it)] = it;

         }

       }

       return false;

     }

     /// \brief An augmenting phase of the Hopcroft-Karp algorithm

     ///

     /// It runs an augmenting phase of the Hopcroft-Karp

     /// algorithm. The phase finds maximum count of edge disjoint

     /// augmenting paths and augments on these paths. The algorithm

     /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is

     /// \f$ O(e) \f$ long.

     bool augment() {

       typename Graph::template ANodeMap<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);

       std::vector<Node> queue, bqueue;

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         if (anode_matching[it] == INVALID) {

           queue.push_back(it);

           areached[it] = true;

         }

       }

       bool success = false;

       while (!success && !queue.empty()) {

         std::vector<Node> newqueue;

         for (int i = 0; i < (int)queue.size(); ++i) {

           Node anode = queue[i];

           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {

             Node bnode = graph->bNode(jt);

             if (breached[bnode]) continue;

             breached[bnode] = true;

             bpred[bnode] = jt;

             if (bnode_matching[bnode] == INVALID) {

               bqueue.push_back(bnode);

               success = true;

             } else {           

               Node newanode = graph->aNode(bnode_matching[bnode]);

               if (!areached[newanode]) {

                 areached[newanode] = true;

                 newqueue.push_back(newanode);

               }

             }

           }

         }

         queue.swap(newqueue);

       }

       if (success) {

         typename Graph::template ANodeMap<bool> aused(*graph, false);

         for (int i = 0; i < (int)bqueue.size(); ++i) {

           Node bnode = bqueue[i];

           bool used = false;

           while (bnode != INVALID) {

             UEdge uedge = bpred[bnode];

             Node anode = graph->aNode(uedge);

             if (aused[anode]) {

               used = true;

               break;

             }

             bnode = anode_matching[anode] != INVALID ? 

               graph->bNode(anode_matching[anode]) : INVALID;

           }

           if (used) continue;

           bnode = bqueue[i];

           while (bnode != INVALID) {

             UEdge uedge = bpred[bnode];

             Node anode = graph->aNode(uedge);

             bnode_matching[bnode] = uedge;

             bnode = anode_matching[anode] != INVALID ? 

               graph->bNode(anode_matching[anode]) : INVALID;

             anode_matching[anode] = uedge;

             aused[anode] = true;

           }

           ++matching_size;

         }

       }

       return success;

     }

     /// \brief An augmenting phase of the Ford-Fulkerson algorithm

     ///

     /// It runs an augmenting phase of the Ford-Fulkerson

     /// algorithm. The phase finds only one augmenting path and 

     /// augments only on this paths. The algorithm consists at most 

     /// of \f$ O(n) \f$ simple phase and one phase is at most 

     /// \f$ O(e) \f$ long.

     bool simpleAugment() {

       typename Graph::template ANodeMap<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);

       std::vector<Node> queue;

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         if (anode_matching[it] == INVALID) {

           queue.push_back(it);

           areached[it] = true;

         }

       }

       while (!queue.empty()) {

         std::vector<Node> newqueue;

         for (int i = 0; i < (int)queue.size(); ++i) {

           Node anode = queue[i];

           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {

             Node bnode = graph->bNode(jt);

             if (breached[bnode]) continue;

             breached[bnode] = true;

             bpred[bnode] = jt;

             if (bnode_matching[bnode] == INVALID) {

               while (bnode != INVALID) {

                 UEdge uedge = bpred[bnode];

                 anode = graph->aNode(uedge);

                 bnode_matching[bnode] = uedge;

                 bnode = anode_matching[anode] != INVALID ? 

                   graph->bNode(anode_matching[anode]) : INVALID;

                 anode_matching[anode] = uedge;

               }

               ++matching_size;

               return true;

             } else {           

               Node newanode = graph->aNode(bnode_matching[bnode]);

               if (!areached[newanode]) {

                 areached[newanode] = true;

                 newqueue.push_back(newanode);

               }

             }

           }

         }

         queue.swap(newqueue);

       }

       return false;

     }

     /// \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() {

       greedyInit();

       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 Returns an minimum covering of the nodes.

     ///

     /// The minimum covering set problem is the dual solution of the

     /// maximum bipartite matching. It provides an solution for this

     /// problem what is proof of the optimality of the matching.

     /// \return The size of the cover set.

     template <typename CoverMap>

     int coverSet(CoverMap& covering) const {

       typename Graph::template ANodeMap<bool> areached(*graph, false);

       typename Graph::template BNodeMap<bool> breached(*graph, false);

       std::vector<Node> queue;

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         if (anode_matching[it] == INVALID) {

           queue.push_back(it);

         }

       }

       while (!queue.empty()) {

         std::vector<Node> newqueue;

         for (int i = 0; i < (int)queue.size(); ++i) {

           Node anode = queue[i];

           for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {

             Node bnode = graph->bNode(jt);

             if (breached[bnode]) continue;

             breached[bnode] = true;

             if (bnode_matching[bnode] != INVALID) {

               Node newanode = graph->aNode(bnode_matching[bnode]);

               if (!areached[newanode]) {

                 areached[newanode] = true;

                 newqueue.push_back(newanode);

               }

             }

           }

         }

         queue.swap(newqueue);

       }

       int size = 0;

       for (ANodeIt it(*graph); it != INVALID; ++it) {

         covering[it] = !areached[it] && anode_matching[it] != INVALID;

         if (!areached[it] && anode_matching[it] != INVALID) {

           ++size;

         }

       }

       for (BNodeIt it(*graph); it != INVALID; ++it) {

         covering[it] = breached[it];

         if (breached[it]) {

           ++size;

         }

       }

       return size;

     }

     /// \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) const {

       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) const {

       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) const {

       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) const {

       if (graph->aNode(node)) {

         return anode_matching[node];

       } else {

         return bnode_matching[node];

       }

     }

     /// \brief Gives back the number of the matching edges.

     ///

     /// Gives back the number of the matching edges.

     int matchingSize() const {

       return matching_size;

     }

     /// @}

   private:

     ANodeMatchingMap anode_matching;

     BNodeMatchingMap bnode_matching;

     const Graph *graph;

     int matching_size;

   };

   /// \ingroup matching

   ///

   /// \brief Maximum cardinality bipartite matching

   ///

   /// This function calculates the maximum cardinality matching

   /// in a bipartite graph. It gives back the matching in an undirected

   /// edge map.

   ///

   /// \param graph The bipartite graph.

   /// \retval matching The undirected edge map which will be set to 

   /// the matching.

   /// \return The size of the matching.

   template <typename BpUGraph, typename MatchingMap>

   int maxBipartiteMatching(const BpUGraph& graph, MatchingMap& matching) {

     MaxBipartiteMatching<BpUGraph> bpmatching(graph);

     bpmatching.run();

     bpmatching.matching(matching);

     return bpmatching.matchingSize();

   }

   /// \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* what() const throw() {

 	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  - (*weight)[jt] +

             anode_potential[anode] + bnode_potential[bnode];

           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 || 

                 bnode_potential[bnode] - bvalue > bestValue) {

               bestValue = bnode_potential[bnode] - bvalue;

               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 matching is optimal

     /// with the current number of edges 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) const {

       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) const {

       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) const {

       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) const {

       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) const {

       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;

   };

   /// \ingroup matching

   ///

   /// \brief Maximum weighted bipartite matching

   ///

   /// This function calculates the maximum weighted matching

   /// in a bipartite graph. It gives back the matching in an undirected

   /// edge map.

   ///

   /// \param graph The bipartite graph.

   /// \param weight The undirected edge map which contains the weights.

   /// \retval matching The undirected edge map which will be set to 

   /// the matching.

   /// \return The value of the matching.

   template <typename BpUGraph, typename WeightMap, typename MatchingMap>

   typename WeightMap::Value 

   maxWeightedBipartiteMatching(const BpUGraph& graph, const WeightMap& weight,

                                MatchingMap& matching) {

     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 

       bpmatching(graph, weight);

     bpmatching.run();

     bpmatching.matching(matching);

     return bpmatching.matchingValue();

   }

   /// \ingroup matching

   ///

   /// \brief Maximum weighted maximum cardinality bipartite matching

   ///

   /// This function calculates the maximum weighted of the maximum cardinality

   /// matchings of a bipartite graph. It gives back the matching in an 

   /// undirected edge map.

   ///

   /// \param graph The bipartite graph.

   /// \param weight The undirected edge map which contains the weights.

   /// \retval matching The undirected edge map which will be set to 

   /// the matching.

   /// \return The value of the matching.

   template <typename BpUGraph, typename WeightMap, typename MatchingMap>

   typename WeightMap::Value 

   maxWeightedMaxBipartiteMatching(const BpUGraph& graph, 

                                   const WeightMap& weight,

                                   MatchingMap& matching) {

     MaxWeightedBipartiteMatching<BpUGraph, WeightMap> 

       bpmatching(graph, weight);

     bpmatching.run(true);

     bpmatching.matching(matching);

     return bpmatching.matchingValue();

   }

   /// \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* what() const throw() {

 	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);

         }

       }

       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 + (*cost)[jt] + 

             anode_potential[anode] - bnode_potential[bnode];

           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) {

         return false;

       }

       matching_cost += 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.

     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 optimal with 

     /// the current number of edges 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) const {

       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) const {

       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) const {

       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) const {

       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) const {

       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;

   };

   /// \ingroup matching

   ///

   /// \brief Minimum cost maximum cardinality bipartite matching

   ///

   /// This function calculates the minimum cost matching of the maximum 

   /// cardinality matchings of a bipartite graph. It gives back the matching 

   /// in an undirected edge map.

   ///

   /// \param graph The bipartite graph.

   /// \param cost The undirected edge map which contains the costs.

   /// \retval matching The undirected edge map which will be set to 

   /// the matching.

   /// \return The cost of the matching.

   template <typename BpUGraph, typename CostMap, typename MatchingMap>

   typename CostMap::Value 

   minCostMaxBipartiteMatching(const BpUGraph& graph, 

                               const CostMap& cost,

                               MatchingMap& matching) {

     MinCostMaxBipartiteMatching<BpUGraph, CostMap> 

       bpmatching(graph, cost);

     bpmatching.run();

     bpmatching.matching(matching);

     return bpmatching.matchingCost();

   }

 }

 #endif