/* -*- 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* exceptionName() const {
	return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
      }
    };

    ///\name Named template parameters

    ///@{

    template <class H, class CR>
    struct DefHeapTraits : public Traits {
      typedef CR HeapCrossRef;
      typedef H Heap;
      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
	throw UninitializedParameter();
      }
      static Heap *createHeap(HeapCrossRef &) {
	throw UninitializedParameter();
      }
    };

    /// \brief \ref named-templ-param "Named parameter" for setting heap 
    /// and cross reference type
    ///
    /// \ref named-templ-param "Named parameter" for setting heap and cross 
    /// reference type
    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
    struct DefHeap
      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
                                            DefHeapTraits<H, CR> > { 
      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
                                           DefHeapTraits<H, CR> > Create;
    };

    template <class H, class CR>
    struct DefStandardHeapTraits : public Traits {
      typedef CR HeapCrossRef;
      typedef H Heap;
      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
	return new HeapCrossRef(graph);
      }
      static Heap *createHeap(HeapCrossRef &crossref) {
	return new Heap(crossref);
      }
    };

    /// \brief \ref named-templ-param "Named parameter" for setting heap and 
    /// cross reference type with automatic allocation
    ///
    /// \ref named-templ-param "Named parameter" for setting heap and cross 
    /// reference type. It can allocate the heap and the cross reference 
    /// object if the cross reference's constructor waits for the graph as 
    /// parameter and the heap's constructor waits for the cross reference.
    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
    struct DefStandardHeap
      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
                                            DefStandardHeapTraits<H, CR> > { 
      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap, 
                                           DefStandardHeapTraits<H, CR> > 
      Create;
    };

    ///@}


    /// \brief Constructor.
    ///
    /// Constructor of the algorithm. 
    MaxWeightedBipartiteMatching(const BpUGraph& _graph, 
                                 const WeightMap& _weight) 
      : graph(&_graph), weight(&_weight),
        anode_matching(_graph), bnode_matching(_graph),
        anode_potential(_graph), bnode_potential(_graph),
        _heap_cross_ref(0), local_heap_cross_ref(false),
        _heap(0), local_heap(0) {}

    /// \brief Destructor.
    ///
    /// Destructor of the algorithm.
    ~MaxWeightedBipartiteMatching() {
      destroyStructures();
    }

    /// \brief Sets the heap and the cross reference used by algorithm.
    ///
    /// Sets the heap and the cross reference used by algorithm.
    /// If you don't use this function before calling \ref run(),
    /// it will allocate one. The destuctor deallocates this
    /// automatically allocated map, of course.
    /// \return \c (*this)
    MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
      if(local_heap_cross_ref) {
	delete _heap_cross_ref;
	local_heap_cross_ref = false;
      }
      _heap_cross_ref = &crossRef;
      if(local_heap) {
	delete _heap;
	local_heap = false;
      }
      _heap = &heap;
      return *this;
    }

    /// \name Execution control
    /// The simplest way to execute the algorithm is to use
    /// one of the member functions called \c run().
    /// \n
    /// If you need more control on the execution,
    /// first you must call \ref init() or one alternative for it.
    /// Finally \ref start() will perform the matching computation or
    /// with step-by-step execution you can augment the solution.

    /// @{

    /// \brief Initalize the data structures.
    ///
    /// It initalizes the data structures and creates an empty matching.
    void init() {
      initStructures();
      for (ANodeIt it(*graph); it != INVALID; ++it) {
        anode_matching[it] = INVALID;
        anode_potential[it] = 0;
      }
      for (BNodeIt it(*graph); it != INVALID; ++it) {
        bnode_matching[it] = INVALID;
        bnode_potential[it] = 0;
        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
          if ((*weight)[jt] > bnode_potential[it]) {
            bnode_potential[it] = (*weight)[jt];
          }
        }
      }
      matching_value = 0;
      matching_size = 0;
    }


    /// \brief An augmenting phase of the weighted matching algorithm
    ///
    /// It runs an augmenting phase of the weighted matching 
    /// algorithm. The phase finds the best augmenting path and 
    /// augments only on this paths. 
    ///
    /// The algorithm consists at most 
    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
    /// with binary heap.
    /// \param decrease If the given parameter true the matching value
    /// can be decreased in the augmenting phase. If we would like
    /// to calculate the maximum cardinality maximum weighted matching
    /// then we should let the algorithm to decrease the matching
    /// value in order to increase the number of the matching edges.
    bool augment(bool decrease = false) {

      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);

      Node bestNode = INVALID;
      Value bestValue = 0;

      _heap->clear();
      for (ANodeIt it(*graph); it != INVALID; ++it) {
        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
      }

      for (ANodeIt it(*graph); it != INVALID; ++it) {
        if (anode_matching[it] == INVALID) {
          _heap->push(it, 0);
        }
      }

      Value bdistMax = 0;
      while (!_heap->empty()) {
        Node anode = _heap->top();
        Value avalue = _heap->prio();
        _heap->pop();
        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
          if (jt == anode_matching[anode]) continue;
          Node bnode = graph->bNode(jt);
          Value bvalue = avalue  - (*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* exceptionName() const {
	return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
      }
    };

    ///\name Named template parameters

    ///@{

    template <class H, class CR>
    struct DefHeapTraits : public Traits {
      typedef CR HeapCrossRef;
      typedef H Heap;
      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
	throw UninitializedParameter();
      }
      static Heap *createHeap(HeapCrossRef &) {
	throw UninitializedParameter();
      }
    };

    /// \brief \ref named-templ-param "Named parameter" for setting heap 
    /// and cross reference type
    ///
    /// \ref named-templ-param "Named parameter" for setting heap and cross 
    /// reference type
    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
    struct DefHeap
      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
                                            DefHeapTraits<H, CR> > { 
      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
                                           DefHeapTraits<H, CR> > Create;
    };

    template <class H, class CR>
    struct DefStandardHeapTraits : public Traits {
      typedef CR HeapCrossRef;
      typedef H Heap;
      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
	return new HeapCrossRef(graph);
      }
      static Heap *createHeap(HeapCrossRef &crossref) {
	return new Heap(crossref);
      }
    };

    /// \brief \ref named-templ-param "Named parameter" for setting heap and 
    /// cross reference type with automatic allocation
    ///
    /// \ref named-templ-param "Named parameter" for setting heap and cross 
    /// reference type. It can allocate the heap and the cross reference 
    /// object if the cross reference's constructor waits for the graph as 
    /// parameter and the heap's constructor waits for the cross reference.
    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
    struct DefStandardHeap
      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
                                            DefStandardHeapTraits<H, CR> > { 
      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap, 
                                           DefStandardHeapTraits<H, CR> > 
      Create;
    };

    ///@}


    /// \brief Constructor.
    ///
    /// Constructor of the algorithm. 
    MinCostMaxBipartiteMatching(const BpUGraph& _graph, 
                                 const CostMap& _cost) 
      : graph(&_graph), cost(&_cost),
        anode_matching(_graph), bnode_matching(_graph),
        anode_potential(_graph), bnode_potential(_graph),
        _heap_cross_ref(0), local_heap_cross_ref(false),
        _heap(0), local_heap(0) {}

    /// \brief Destructor.
    ///
    /// Destructor of the algorithm.
    ~MinCostMaxBipartiteMatching() {
      destroyStructures();
    }

    /// \brief Sets the heap and the cross reference used by algorithm.
    ///
    /// Sets the heap and the cross reference used by algorithm.
    /// If you don't use this function before calling \ref run(),
    /// it will allocate one. The destuctor deallocates this
    /// automatically allocated map, of course.
    /// \return \c (*this)
    MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
      if(local_heap_cross_ref) {
	delete _heap_cross_ref;
	local_heap_cross_ref = false;
      }
      _heap_cross_ref = &crossRef;
      if(local_heap) {
	delete _heap;
	local_heap = false;
      }
      _heap = &heap;
      return *this;
    }

    /// \name Execution control
    /// The simplest way to execute the algorithm is to use
    /// one of the member functions called \c run().
    /// \n
    /// If you need more control on the execution,
    /// first you must call \ref init() or one alternative for it.
    /// Finally \ref start() will perform the matching computation or
    /// with step-by-step execution you can augment the solution.

    /// @{

    /// \brief Initalize the data structures.
    ///
    /// It initalizes the data structures and creates an empty matching.
    void init() {
      initStructures();
      for (ANodeIt it(*graph); it != INVALID; ++it) {
        anode_matching[it] = INVALID;
        anode_potential[it] = 0;
      }
      for (BNodeIt it(*graph); it != INVALID; ++it) {
        bnode_matching[it] = INVALID;
        bnode_potential[it] = 0;
      }
      matching_cost = 0;
      matching_size = 0;
    }


    /// \brief An augmenting phase of the costed matching algorithm
    ///
    /// It runs an augmenting phase of the matching algorithm. The
    /// phase finds the best augmenting path and augments only on this
    /// paths.
    ///
    /// The algorithm consists at most 
    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$ 
    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long 
    /// with binary heap.
    bool augment() {

      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);

      Node bestNode = INVALID;
      Value bestValue = 0;

      _heap->clear();
      for (ANodeIt it(*graph); it != INVALID; ++it) {
        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
      }

      for (ANodeIt it(*graph); it != INVALID; ++it) {
        if (anode_matching[it] == INVALID) {
          _heap->push(it, 0);
        }
      }

      while (!_heap->empty()) {
        Node anode = _heap->top();
        Value avalue = _heap->prio();
        _heap->pop();
        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
          if (jt == anode_matching[anode]) continue;
          Node bnode = graph->bNode(jt);
          Value bvalue = avalue + (*cost)[jt] + 
            anode_potential[anode] - bnode_potential[bnode];
          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
            bdist[bnode] = bvalue;
            bpred[bnode] = jt;
          }
          if (bnode_matching[bnode] != INVALID) {
            Node newanode = graph->aNode(bnode_matching[bnode]);
            switch (_heap->state(newanode)) {
            case Heap::PRE_HEAP:
              _heap->push(newanode, bvalue);
              break;
            case Heap::IN_HEAP:
              if (bvalue < (*_heap)[newanode]) {
                _heap->decrease(newanode, bvalue);
              }
              break;
            case Heap::POST_HEAP:
              break;
            }
          } else {
            if (bestNode == INVALID || 
                bvalue + bnode_potential[bnode] < bestValue) {
              bestValue = bvalue + bnode_potential[bnode];
              bestNode = bnode;
            }
          }
        }
      }

      if (bestNode == INVALID) {
        return false;
      }

      matching_cost += bestValue;
      ++matching_size;

      for (BNodeIt it(*graph); it != INVALID; ++it) {
        if (bpred[it] != INVALID) {
          bnode_potential[it] += bdist[it];
        }
      }
      for (ANodeIt it(*graph); it != INVALID; ++it) {
        if (anode_matching[it] != INVALID) {
          Node bnode = graph->bNode(anode_matching[it]);
          if (bpred[bnode] != INVALID) {
            anode_potential[it] += bdist[bnode];
          }
        }
      }

      while (bestNode != INVALID) {
        UEdge uedge = bpred[bestNode];
        Node anode = graph->aNode(uedge);
        
        bnode_matching[bestNode] = uedge;
        if (anode_matching[anode] != INVALID) {
          bestNode = graph->bNode(anode_matching[anode]);
        } else {
          bestNode = INVALID;
        }
        anode_matching[anode] = uedge;
      }


      return true;
    }

    /// \brief Starts the algorithm.
    ///
    /// Starts the algorithm. It runs augmenting phases until the
    /// optimal solution reached.
    void start() {
      while (augment()) {}
    }

    /// \brief Runs the algorithm.
    ///
    /// It just initalize the algorithm and then start it.
    void run() {
      init();
      start();
    }

    /// @}

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

    /// \brief Gives back the potential in the NodeMap
    ///
    /// Gives back the potential in the NodeMap. The potential is 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