COIN-OR::LEMON - Graph Library

Changeset 2051:08652c1763f6 in lemon-0.x for lemon


Ignore:
Timestamp:
04/14/06 20:07:33 (14 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@2694
Message:

MaxWeightedBipartiteMatching?
MinCostMaxBipartiteMatching?

Both algorithms are based on successive shortest
path algorithm with dijkstra shortest path
finding

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/bipartite_matching.h

    r2040 r2051  
    2020#define LEMON_BIPARTITE_MATCHING
    2121
    22 #include <lemon/bpugraph_adaptor.h>
    23 #include <lemon/bfs.h>
     22#include <functional>
     23
     24#include <lemon/bin_heap.h>
     25#include <lemon/maps.h>
    2426
    2527#include <iostream>
     
    3638  ///
    3739  /// Bipartite Max Cardinality Matching algorithm. This class implements
    38   /// the Hopcroft-Karp algorithm wich has \f$ O(e\sqrt{n}) \f$ time
     40  /// the Hopcroft-Karp algorithm which has \f$ O(e\sqrt{n}) \f$ time
    3941  /// complexity.
    4042  template <typename BpUGraph>
     
    8486        bnode_matching[it] = INVALID;
    8587      }
    86       matching_value = 0;
     88      matching_size = 0;
    8789    }
    8890
     
    9395    /// the matching than from the initial empty matching.
    9496    void greedyInit() {
    95       matching_value = 0;
     97      matching_size = 0;
    9698      for (BNodeIt it(*graph); it != INVALID; ++it) {
    9799        bnode_matching[it] = INVALID;
     
    103105            anode_matching[it] = jt;
    104106            bnode_matching[graph->bNode(jt)] = jt;
    105             ++matching_value;
     107            ++matching_size;
    106108            break;
    107109          }
     
    121123        bnode_matching[it] = INVALID;
    122124      }
    123       matching_value = 0;
     125      matching_size = 0;
    124126      for (UEdgeIt it(*graph); it != INVALID; ++it) {
    125127        if (matching[it]) {
    126           ++matching_value;
     128          ++matching_size;
    127129          anode_matching[graph->aNode(it)] = it;
    128130          bnode_matching[graph->bNode(it)] = it;
     
    143145        bnode_matching[it] = INVALID;
    144146      }
    145       matching_value = 0;
     147      matching_size = 0;
    146148      for (UEdgeIt it(*graph); it != INVALID; ++it) {
    147149        if (matching[it]) {
    148           ++matching_value;
     150          ++matching_size;
    149151          if (anode_matching[graph->aNode(it)] != INVALID) {
    150152            return false;
     
    247249            aused[anode] = true;
    248250          }
    249           ++matching_value;
     251          ++matching_size;
    250252
    251253        }
     
    298300               
    299301              }
    300               ++matching_value;
     302              ++matching_size;
    301303              return true;
    302304            } else {           
     
    408410        }
    409411      }
    410       return matching_value;
     412      return matching_size;
    411413    }
    412414
     
    420422        matching[it] = it == anode_matching[graph->aNode(it)];
    421423      }
    422       return matching_value;
     424      return matching_size;
    423425    }
    424426
     
    446448    ///
    447449    /// Gives back the number of the matching edges.
    448     int matchingValue() const {
    449       return matching_value;
     450    int matchingSize() const {
     451      return matching_size;
    450452    }
    451453
     
    458460    const Graph *graph;
    459461
    460     int matching_value;
     462    int matching_size;
    461463 
    462464  };
    463465
     466  /// \brief Default traits class for weighted bipartite matching algoritms.
     467  ///
     468  /// Default traits class for weighted bipartite matching algoritms.
     469  /// \param _BpUGraph The bipartite undirected graph type.
     470  /// \param _WeightMap Type of weight map.
     471  template <typename _BpUGraph, typename _WeightMap>
     472  struct WeightedBipartiteMatchingDefaultTraits {
     473    /// \brief The type of the weight of the undirected edges.
     474    typedef typename _WeightMap::Value Value;
     475
     476    /// The undirected bipartite graph type the algorithm runs on.
     477    typedef _BpUGraph BpUGraph;
     478
     479    /// The map of the edges weights
     480    typedef _WeightMap WeightMap;
     481
     482    /// \brief The cross reference type used by heap.
     483    ///
     484    /// The cross reference type used by heap.
     485    /// Usually it is \c Graph::NodeMap<int>.
     486    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
     487
     488    /// \brief Instantiates a HeapCrossRef.
     489    ///
     490    /// This function instantiates a \ref HeapCrossRef.
     491    /// \param graph is the graph, to which we would like to define the
     492    /// HeapCrossRef.
     493    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
     494      return new HeapCrossRef(graph);
     495    }
     496   
     497    /// \brief The heap type used by weighted matching algorithms.
     498    ///
     499    /// The heap type used by weighted matching algorithms. It should
     500    /// minimize the priorities and the heap's key type is the graph's
     501    /// anode graph's node.
     502    ///
     503    /// \sa BinHeap
     504    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
     505   
     506    /// \brief Instantiates a Heap.
     507    ///
     508    /// This function instantiates a \ref Heap.
     509    /// \param crossref The cross reference of the heap.
     510    static Heap *createHeap(HeapCrossRef& crossref) {
     511      return new Heap(crossref);
     512    }
     513
     514  };
     515
     516
     517  /// \ingroup matching
     518  ///
     519  /// \brief Bipartite Max Weighted Matching algorithm
     520  ///
     521  /// This class implements the bipartite Max Weighted Matching
     522  /// algorithm.  It uses the successive shortest path algorithm to
     523  /// calculate the maximum weighted matching in the bipartite
     524  /// graph. The algorithm can be used also to calculate the maximum
     525  /// cardinality maximum weighted matching. The time complexity
     526  /// of the algorithm is \f$ O(ne\log(n)) \f$ with the default binary
     527  /// heap implementation but this can be improved to
     528  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
     529  ///
     530  /// The algorithm also provides a potential function on the nodes
     531  /// which a dual solution of the matching algorithm and it can be
     532  /// used to proof the optimality of the given pimal solution.
     533#ifdef DOXYGEN
     534  template <typename _BpUGraph, typename _WeightMap, typename _Traits>
     535#else
     536  template <typename _BpUGraph,
     537            typename _WeightMap = typename _BpUGraph::template UEdgeMap<int>,
     538            typename _Traits = WeightedBipartiteMatchingDefaultTraits<_BpUGraph, _WeightMap> >
     539#endif
     540  class MaxWeightedBipartiteMatching {
     541  public:
     542
     543    typedef _Traits Traits;
     544    typedef typename Traits::BpUGraph BpUGraph;
     545    typedef typename Traits::WeightMap WeightMap;
     546    typedef typename Traits::Value Value;
     547
     548  protected:
     549
     550    typedef typename Traits::HeapCrossRef HeapCrossRef;
     551    typedef typename Traits::Heap Heap;
     552
     553   
     554    typedef typename BpUGraph::Node Node;
     555    typedef typename BpUGraph::ANodeIt ANodeIt;
     556    typedef typename BpUGraph::BNodeIt BNodeIt;
     557    typedef typename BpUGraph::UEdge UEdge;
     558    typedef typename BpUGraph::UEdgeIt UEdgeIt;
     559    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
     560
     561    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
     562    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
     563
     564    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
     565    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
     566
     567
     568  public:
     569
     570    /// \brief \ref Exception for uninitialized parameters.
     571    ///
     572    /// This error represents problems in the initialization
     573    /// of the parameters of the algorithms.
     574    class UninitializedParameter : public lemon::UninitializedParameter {
     575    public:
     576      virtual const char* exceptionName() const {
     577        return "lemon::MaxWeightedBipartiteMatching::UninitializedParameter";
     578      }
     579    };
     580
     581    ///\name Named template parameters
     582
     583    ///@{
     584
     585    template <class H, class CR>
     586    struct DefHeapTraits : public Traits {
     587      typedef CR HeapCrossRef;
     588      typedef H Heap;
     589      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
     590        throw UninitializedParameter();
     591      }
     592      static Heap *createHeap(HeapCrossRef &) {
     593        throw UninitializedParameter();
     594      }
     595    };
     596
     597    /// \brief \ref named-templ-param "Named parameter" for setting heap
     598    /// and cross reference type
     599    ///
     600    /// \ref named-templ-param "Named parameter" for setting heap and cross
     601    /// reference type
     602    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
     603    struct DefHeap
     604      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
     605                                            DefHeapTraits<H, CR> > {
     606      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
     607                                           DefHeapTraits<H, CR> > Create;
     608    };
     609
     610    template <class H, class CR>
     611    struct DefStandardHeapTraits : public Traits {
     612      typedef CR HeapCrossRef;
     613      typedef H Heap;
     614      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
     615        return new HeapCrossRef(graph);
     616      }
     617      static Heap *createHeap(HeapCrossRef &crossref) {
     618        return new Heap(crossref);
     619      }
     620    };
     621
     622    /// \brief \ref named-templ-param "Named parameter" for setting heap and
     623    /// cross reference type with automatic allocation
     624    ///
     625    /// \ref named-templ-param "Named parameter" for setting heap and cross
     626    /// reference type. It can allocate the heap and the cross reference
     627    /// object if the cross reference's constructor waits for the graph as
     628    /// parameter and the heap's constructor waits for the cross reference.
     629    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
     630    struct DefStandardHeap
     631      : public MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
     632                                            DefStandardHeapTraits<H, CR> > {
     633      typedef MaxWeightedBipartiteMatching<BpUGraph, WeightMap,
     634                                           DefStandardHeapTraits<H, CR> >
     635      Create;
     636    };
     637
     638    ///@}
     639
     640
     641    /// \brief Constructor.
     642    ///
     643    /// Constructor of the algorithm.
     644    MaxWeightedBipartiteMatching(const BpUGraph& _graph,
     645                                 const WeightMap& _weight)
     646      : graph(&_graph), weight(&_weight),
     647        anode_matching(_graph), bnode_matching(_graph),
     648        anode_potential(_graph), bnode_potential(_graph),
     649        _heap_cross_ref(0), local_heap_cross_ref(false),
     650        _heap(0), local_heap(0) {}
     651
     652    /// \brief Destructor.
     653    ///
     654    /// Destructor of the algorithm.
     655    ~MaxWeightedBipartiteMatching() {
     656      destroyStructures();
     657    }
     658
     659    /// \brief Sets the heap and the cross reference used by algorithm.
     660    ///
     661    /// Sets the heap and the cross reference used by algorithm.
     662    /// If you don't use this function before calling \ref run(),
     663    /// it will allocate one. The destuctor deallocates this
     664    /// automatically allocated map, of course.
     665    /// \return \c (*this)
     666    MaxWeightedBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
     667      if(local_heap_cross_ref) {
     668        delete _heap_cross_ref;
     669        local_heap_cross_ref = false;
     670      }
     671      _heap_cross_ref = &crossRef;
     672      if(local_heap) {
     673        delete _heap;
     674        local_heap = false;
     675      }
     676      _heap = &heap;
     677      return *this;
     678    }
     679
     680    /// \name Execution control
     681    /// The simplest way to execute the algorithm is to use
     682    /// one of the member functions called \c run().
     683    /// \n
     684    /// If you need more control on the execution,
     685    /// first you must call \ref init() or one alternative for it.
     686    /// Finally \ref start() will perform the matching computation or
     687    /// with step-by-step execution you can augment the solution.
     688
     689    /// @{
     690
     691    /// \brief Initalize the data structures.
     692    ///
     693    /// It initalizes the data structures and creates an empty matching.
     694    void init() {
     695      initStructures();
     696      for (ANodeIt it(*graph); it != INVALID; ++it) {
     697        anode_matching[it] = INVALID;
     698        anode_potential[it] = 0;
     699      }
     700      for (BNodeIt it(*graph); it != INVALID; ++it) {
     701        bnode_matching[it] = INVALID;
     702        bnode_potential[it] = 0;
     703        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
     704          if (-(*weight)[jt] <= bnode_potential[it]) {
     705            bnode_potential[it] = -(*weight)[jt];
     706          }
     707        }
     708      }
     709      matching_value = 0;
     710      matching_size = 0;
     711    }
     712
     713
     714    /// \brief An augmenting phase of the weighted matching algorithm
     715    ///
     716    /// It runs an augmenting phase of the weighted matching
     717    /// algorithm. The phase finds the best augmenting path and
     718    /// augments only on this paths.
     719    ///
     720    /// The algorithm consists at most
     721    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
     722    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
     723    /// with binary heap.
     724    /// \param decrease If the given parameter true the matching value
     725    /// can be decreased in the augmenting phase. If we would like
     726    /// to calculate the maximum cardinality maximum weighted matching
     727    /// then we should let the algorithm to decrease the matching
     728    /// value in order to increase the number of the matching edges.
     729    bool augment(bool decrease = false) {
     730
     731      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
     732      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
     733
     734      Node bestNode = INVALID;
     735      Value bestValue = 0;
     736
     737      _heap->clear();
     738      for (ANodeIt it(*graph); it != INVALID; ++it) {
     739        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
     740      }
     741
     742      for (ANodeIt it(*graph); it != INVALID; ++it) {
     743        if (anode_matching[it] == INVALID) {
     744          _heap->push(it, 0);
     745        }
     746      }
     747
     748      Value bdistMax = 0;
     749      while (!_heap->empty()) {
     750        Node anode = _heap->top();
     751        Value avalue = _heap->prio();
     752        _heap->pop();
     753        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
     754          if (jt == anode_matching[anode]) continue;
     755          Node bnode = graph->bNode(jt);
     756          Value bvalue = avalue + anode_potential[anode] -
     757            bnode_potential[bnode] - (*weight)[jt];
     758          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
     759            bdist[bnode] = bvalue;
     760            bpred[bnode] = jt;
     761          }
     762          if (bvalue > bdistMax) {
     763            bdistMax = bvalue;
     764          }
     765          if (bnode_matching[bnode] != INVALID) {
     766            Node newanode = graph->aNode(bnode_matching[bnode]);
     767            switch (_heap->state(newanode)) {
     768            case Heap::PRE_HEAP:
     769              _heap->push(newanode, bvalue);
     770              break;
     771            case Heap::IN_HEAP:
     772              if (bvalue < (*_heap)[newanode]) {
     773                _heap->decrease(newanode, bvalue);
     774              }
     775              break;
     776            case Heap::POST_HEAP:
     777              break;
     778            }
     779          } else {
     780            if (bestNode == INVALID ||
     781                - bvalue - bnode_potential[bnode] > bestValue) {
     782              bestValue = - bvalue - bnode_potential[bnode];
     783              bestNode = bnode;
     784            }
     785          }
     786        }
     787      }
     788
     789      if (bestNode == INVALID || (!decrease && bestValue < 0)) {
     790        return false;
     791      }
     792
     793      matching_value += bestValue;
     794      ++matching_size;
     795
     796      for (BNodeIt it(*graph); it != INVALID; ++it) {
     797        if (bpred[it] != INVALID) {
     798          bnode_potential[it] += bdist[it];
     799        } else {
     800          bnode_potential[it] += bdistMax;
     801        }
     802      }
     803      for (ANodeIt it(*graph); it != INVALID; ++it) {
     804        if (anode_matching[it] != INVALID) {
     805          Node bnode = graph->bNode(anode_matching[it]);
     806          if (bpred[bnode] != INVALID) {
     807            anode_potential[it] += bdist[bnode];
     808          } else {
     809            anode_potential[it] += bdistMax;
     810          }
     811        }
     812      }
     813
     814      while (bestNode != INVALID) {
     815        UEdge uedge = bpred[bestNode];
     816        Node anode = graph->aNode(uedge);
     817       
     818        bnode_matching[bestNode] = uedge;
     819        if (anode_matching[anode] != INVALID) {
     820          bestNode = graph->bNode(anode_matching[anode]);
     821        } else {
     822          bestNode = INVALID;
     823        }
     824        anode_matching[anode] = uedge;
     825      }
     826
     827
     828      return true;
     829    }
     830
     831    /// \brief Starts the algorithm.
     832    ///
     833    /// Starts the algorithm. It runs augmenting phases until the
     834    /// optimal solution reached.
     835    ///
     836    /// \param maxCardinality If the given value is true it will
     837    /// calculate the maximum cardinality maximum matching instead of
     838    /// the maximum matching.
     839    void start(bool maxCardinality = false) {
     840      while (augment(maxCardinality)) {}
     841    }
     842
     843    /// \brief Runs the algorithm.
     844    ///
     845    /// It just initalize the algorithm and then start it.
     846    ///
     847    /// \param maxCardinality If the given value is true it will
     848    /// calculate the maximum cardinality maximum matching instead of
     849    /// the maximum matching.
     850    void run(bool maxCardinality = false) {
     851      init();
     852      start(maxCardinality);
     853    }
     854
     855    /// @}
     856
     857    /// \name Query Functions
     858    /// The result of the %Matching algorithm can be obtained using these
     859    /// functions.\n
     860    /// Before the use of these functions,
     861    /// either run() or start() must be called.
     862   
     863    ///@{
     864
     865    /// \brief Gives back the potential in the NodeMap
     866    ///
     867    /// Gives back the potential in the NodeMap. The potential
     868    /// is feasible if \f$ \pi(a) - \pi(b) - w(ab) = 0 \f$ for
     869    /// each matching edges and \f$ \pi(a) - \pi(b) - w(ab) \ge 0 \f$
     870    /// for each edges.
     871    template <typename PotentialMap>
     872    void potential(PotentialMap& potential) {
     873      for (ANodeIt it(*graph); it != INVALID; ++it) {
     874        potential[it] = anode_potential[it];
     875      }
     876      for (BNodeIt it(*graph); it != INVALID; ++it) {
     877        potential[it] = bnode_potential[it];
     878      }
     879    }
     880
     881    /// \brief Set true all matching uedge in the map.
     882    ///
     883    /// Set true all matching uedge in the map. It does not change the
     884    /// value mapped to the other uedges.
     885    /// \return The number of the matching edges.
     886    template <typename MatchingMap>
     887    int quickMatching(MatchingMap& matching) {
     888      for (ANodeIt it(*graph); it != INVALID; ++it) {
     889        if (anode_matching[it] != INVALID) {
     890          matching[anode_matching[it]] = true;
     891        }
     892      }
     893      return matching_size;
     894    }
     895
     896    /// \brief Set true all matching uedge in the map and the others to false.
     897    ///
     898    /// Set true all matching uedge in the map and the others to false.
     899    /// \return The number of the matching edges.
     900    template <typename MatchingMap>
     901    int matching(MatchingMap& matching) {
     902      for (UEdgeIt it(*graph); it != INVALID; ++it) {
     903        matching[it] = it == anode_matching[graph->aNode(it)];
     904      }
     905      return matching_size;
     906    }
     907
     908
     909    /// \brief Return true if the given uedge is in the matching.
     910    ///
     911    /// It returns true if the given uedge is in the matching.
     912    bool matchingEdge(const UEdge& edge) {
     913      return anode_matching[graph->aNode(edge)] == edge;
     914    }
     915
     916    /// \brief Returns the matching edge from the node.
     917    ///
     918    /// Returns the matching edge from the node. If there is not such
     919    /// edge it gives back \c INVALID.
     920    UEdge matchingEdge(const Node& node) {
     921      if (graph->aNode(node)) {
     922        return anode_matching[node];
     923      } else {
     924        return bnode_matching[node];
     925      }
     926    }
     927
     928    /// \brief Gives back the sum of weights of the matching edges.
     929    ///
     930    /// Gives back the sum of weights of the matching edges.
     931    Value matchingValue() const {
     932      return matching_value;
     933    }
     934
     935    /// \brief Gives back the number of the matching edges.
     936    ///
     937    /// Gives back the number of the matching edges.
     938    int matchingSize() const {
     939      return matching_size;
     940    }
     941
     942    /// @}
     943
     944  private:
     945
     946    void initStructures() {
     947      if (!_heap_cross_ref) {
     948        local_heap_cross_ref = true;
     949        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
     950      }
     951      if (!_heap) {
     952        local_heap = true;
     953        _heap = Traits::createHeap(*_heap_cross_ref);
     954      }
     955    }
     956
     957    void destroyStructures() {
     958      if (local_heap_cross_ref) delete _heap_cross_ref;
     959      if (local_heap) delete _heap;
     960    }
     961
     962
     963  private:
     964   
     965    const BpUGraph *graph;
     966    const WeightMap* weight;
     967
     968    ANodeMatchingMap anode_matching;
     969    BNodeMatchingMap bnode_matching;
     970
     971    ANodePotentialMap anode_potential;
     972    BNodePotentialMap bnode_potential;
     973
     974    Value matching_value;
     975    int matching_size;
     976
     977    HeapCrossRef *_heap_cross_ref;
     978    bool local_heap_cross_ref;
     979
     980    Heap *_heap;
     981    bool local_heap;
     982 
     983  };
     984
     985  /// \brief Default traits class for minimum cost bipartite matching
     986  /// algoritms.
     987  ///
     988  /// Default traits class for minimum cost bipartite matching
     989  /// algoritms. 
     990  ///
     991  /// \param _BpUGraph The bipartite undirected graph
     992  /// type. 
     993  ///
     994  /// \param _CostMap Type of cost map.
     995  template <typename _BpUGraph, typename _CostMap>
     996  struct MinCostMaxBipartiteMatchingDefaultTraits {
     997    /// \brief The type of the cost of the undirected edges.
     998    typedef typename _CostMap::Value Value;
     999
     1000    /// The undirected bipartite graph type the algorithm runs on.
     1001    typedef _BpUGraph BpUGraph;
     1002
     1003    /// The map of the edges costs
     1004    typedef _CostMap CostMap;
     1005
     1006    /// \brief The cross reference type used by heap.
     1007    ///
     1008    /// The cross reference type used by heap.
     1009    /// Usually it is \c Graph::NodeMap<int>.
     1010    typedef typename BpUGraph::template NodeMap<int> HeapCrossRef;
     1011
     1012    /// \brief Instantiates a HeapCrossRef.
     1013    ///
     1014    /// This function instantiates a \ref HeapCrossRef.
     1015    /// \param graph is the graph, to which we would like to define the
     1016    /// HeapCrossRef.
     1017    static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
     1018      return new HeapCrossRef(graph);
     1019    }
     1020   
     1021    /// \brief The heap type used by costed matching algorithms.
     1022    ///
     1023    /// The heap type used by costed matching algorithms. It should
     1024    /// minimize the priorities and the heap's key type is the graph's
     1025    /// anode graph's node.
     1026    ///
     1027    /// \sa BinHeap
     1028    typedef BinHeap<typename BpUGraph::Node, Value, HeapCrossRef> Heap;
     1029   
     1030    /// \brief Instantiates a Heap.
     1031    ///
     1032    /// This function instantiates a \ref Heap.
     1033    /// \param crossref The cross reference of the heap.
     1034    static Heap *createHeap(HeapCrossRef& crossref) {
     1035      return new Heap(crossref);
     1036    }
     1037
     1038  };
     1039
     1040
     1041  /// \ingroup matching
     1042  ///
     1043  /// \brief Bipartite Min Cost Matching algorithm
     1044  ///
     1045  /// This class implements the bipartite Min Cost Matching algorithm.
     1046  /// It uses the successive shortest path algorithm to calculate the
     1047  /// minimum cost maximum matching in the bipartite graph. The time
     1048  /// complexity of the algorithm is \f$ O(ne\log(n)) \f$ with the
     1049  /// default binary heap implementation but this can be improved to
     1050  /// \f$ O(n^2\log(n)+ne) \f$ if we use fibonacci heaps.
     1051  ///
     1052  /// The algorithm also provides a potential function on the nodes
     1053  /// which a dual solution of the matching algorithm and it can be
     1054  /// used to proof the optimality of the given pimal solution.
     1055#ifdef DOXYGEN
     1056  template <typename _BpUGraph, typename _CostMap, typename _Traits>
     1057#else
     1058  template <typename _BpUGraph,
     1059            typename _CostMap = typename _BpUGraph::template UEdgeMap<int>,
     1060            typename _Traits = MinCostMaxBipartiteMatchingDefaultTraits<_BpUGraph, _CostMap> >
     1061#endif
     1062  class MinCostMaxBipartiteMatching {
     1063  public:
     1064
     1065    typedef _Traits Traits;
     1066    typedef typename Traits::BpUGraph BpUGraph;
     1067    typedef typename Traits::CostMap CostMap;
     1068    typedef typename Traits::Value Value;
     1069
     1070  protected:
     1071
     1072    typedef typename Traits::HeapCrossRef HeapCrossRef;
     1073    typedef typename Traits::Heap Heap;
     1074
     1075   
     1076    typedef typename BpUGraph::Node Node;
     1077    typedef typename BpUGraph::ANodeIt ANodeIt;
     1078    typedef typename BpUGraph::BNodeIt BNodeIt;
     1079    typedef typename BpUGraph::UEdge UEdge;
     1080    typedef typename BpUGraph::UEdgeIt UEdgeIt;
     1081    typedef typename BpUGraph::IncEdgeIt IncEdgeIt;
     1082
     1083    typedef typename BpUGraph::template ANodeMap<UEdge> ANodeMatchingMap;
     1084    typedef typename BpUGraph::template BNodeMap<UEdge> BNodeMatchingMap;
     1085
     1086    typedef typename BpUGraph::template ANodeMap<Value> ANodePotentialMap;
     1087    typedef typename BpUGraph::template BNodeMap<Value> BNodePotentialMap;
     1088
     1089
     1090  public:
     1091
     1092    /// \brief \ref Exception for uninitialized parameters.
     1093    ///
     1094    /// This error represents problems in the initialization
     1095    /// of the parameters of the algorithms.
     1096    class UninitializedParameter : public lemon::UninitializedParameter {
     1097    public:
     1098      virtual const char* exceptionName() const {
     1099        return "lemon::MinCostMaxBipartiteMatching::UninitializedParameter";
     1100      }
     1101    };
     1102
     1103    ///\name Named template parameters
     1104
     1105    ///@{
     1106
     1107    template <class H, class CR>
     1108    struct DefHeapTraits : public Traits {
     1109      typedef CR HeapCrossRef;
     1110      typedef H Heap;
     1111      static HeapCrossRef *createHeapCrossRef(const BpUGraph &) {
     1112        throw UninitializedParameter();
     1113      }
     1114      static Heap *createHeap(HeapCrossRef &) {
     1115        throw UninitializedParameter();
     1116      }
     1117    };
     1118
     1119    /// \brief \ref named-templ-param "Named parameter" for setting heap
     1120    /// and cross reference type
     1121    ///
     1122    /// \ref named-templ-param "Named parameter" for setting heap and cross
     1123    /// reference type
     1124    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
     1125    struct DefHeap
     1126      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
     1127                                            DefHeapTraits<H, CR> > {
     1128      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
     1129                                           DefHeapTraits<H, CR> > Create;
     1130    };
     1131
     1132    template <class H, class CR>
     1133    struct DefStandardHeapTraits : public Traits {
     1134      typedef CR HeapCrossRef;
     1135      typedef H Heap;
     1136      static HeapCrossRef *createHeapCrossRef(const BpUGraph &graph) {
     1137        return new HeapCrossRef(graph);
     1138      }
     1139      static Heap *createHeap(HeapCrossRef &crossref) {
     1140        return new Heap(crossref);
     1141      }
     1142    };
     1143
     1144    /// \brief \ref named-templ-param "Named parameter" for setting heap and
     1145    /// cross reference type with automatic allocation
     1146    ///
     1147    /// \ref named-templ-param "Named parameter" for setting heap and cross
     1148    /// reference type. It can allocate the heap and the cross reference
     1149    /// object if the cross reference's constructor waits for the graph as
     1150    /// parameter and the heap's constructor waits for the cross reference.
     1151    template <class H, class CR = typename BpUGraph::template NodeMap<int> >
     1152    struct DefStandardHeap
     1153      : public MinCostMaxBipartiteMatching<BpUGraph, CostMap,
     1154                                            DefStandardHeapTraits<H, CR> > {
     1155      typedef MinCostMaxBipartiteMatching<BpUGraph, CostMap,
     1156                                           DefStandardHeapTraits<H, CR> >
     1157      Create;
     1158    };
     1159
     1160    ///@}
     1161
     1162
     1163    /// \brief Constructor.
     1164    ///
     1165    /// Constructor of the algorithm.
     1166    MinCostMaxBipartiteMatching(const BpUGraph& _graph,
     1167                                 const CostMap& _cost)
     1168      : graph(&_graph), cost(&_cost),
     1169        anode_matching(_graph), bnode_matching(_graph),
     1170        anode_potential(_graph), bnode_potential(_graph),
     1171        _heap_cross_ref(0), local_heap_cross_ref(false),
     1172        _heap(0), local_heap(0) {}
     1173
     1174    /// \brief Destructor.
     1175    ///
     1176    /// Destructor of the algorithm.
     1177    ~MinCostMaxBipartiteMatching() {
     1178      destroyStructures();
     1179    }
     1180
     1181    /// \brief Sets the heap and the cross reference used by algorithm.
     1182    ///
     1183    /// Sets the heap and the cross reference used by algorithm.
     1184    /// If you don't use this function before calling \ref run(),
     1185    /// it will allocate one. The destuctor deallocates this
     1186    /// automatically allocated map, of course.
     1187    /// \return \c (*this)
     1188    MinCostMaxBipartiteMatching& heap(Heap& heap, HeapCrossRef &crossRef) {
     1189      if(local_heap_cross_ref) {
     1190        delete _heap_cross_ref;
     1191        local_heap_cross_ref = false;
     1192      }
     1193      _heap_cross_ref = &crossRef;
     1194      if(local_heap) {
     1195        delete _heap;
     1196        local_heap = false;
     1197      }
     1198      _heap = &heap;
     1199      return *this;
     1200    }
     1201
     1202    /// \name Execution control
     1203    /// The simplest way to execute the algorithm is to use
     1204    /// one of the member functions called \c run().
     1205    /// \n
     1206    /// If you need more control on the execution,
     1207    /// first you must call \ref init() or one alternative for it.
     1208    /// Finally \ref start() will perform the matching computation or
     1209    /// with step-by-step execution you can augment the solution.
     1210
     1211    /// @{
     1212
     1213    /// \brief Initalize the data structures.
     1214    ///
     1215    /// It initalizes the data structures and creates an empty matching.
     1216    void init() {
     1217      initStructures();
     1218      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1219        anode_matching[it] = INVALID;
     1220        anode_potential[it] = 0;
     1221      }
     1222      for (BNodeIt it(*graph); it != INVALID; ++it) {
     1223        bnode_matching[it] = INVALID;
     1224        bnode_potential[it] = 0;
     1225      }
     1226      matching_cost = 0;
     1227      matching_size = 0;
     1228    }
     1229
     1230
     1231    /// \brief An augmenting phase of the costed matching algorithm
     1232    ///
     1233    /// It runs an augmenting phase of the matching algorithm. The
     1234    /// phase finds the best augmenting path and augments only on this
     1235    /// paths.
     1236    ///
     1237    /// The algorithm consists at most
     1238    /// of \f$ O(n) \f$ phase and one phase is \f$ O(n\log(n)+e) \f$
     1239    /// long with Fibonacci heap or \f$ O((n+e)\log(n)) \f$ long
     1240    /// with binary heap.
     1241    bool augment() {
     1242
     1243      typename BpUGraph::template BNodeMap<Value> bdist(*graph);
     1244      typename BpUGraph::template BNodeMap<UEdge> bpred(*graph, INVALID);
     1245
     1246      Node bestNode = INVALID;
     1247      Value bestValue = 0;
     1248
     1249      _heap->clear();
     1250      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1251        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
     1252      }
     1253
     1254      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1255        if (anode_matching[it] == INVALID) {
     1256          _heap->push(it, 0);
     1257        }
     1258      }
     1259
     1260      while (!_heap->empty()) {
     1261        Node anode = _heap->top();
     1262        Value avalue = _heap->prio();
     1263        _heap->pop();
     1264        for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
     1265          if (jt == anode_matching[anode]) continue;
     1266          Node bnode = graph->bNode(jt);
     1267          Value bvalue = avalue + (*cost)[jt] +
     1268            anode_potential[anode] - bnode_potential[bnode];
     1269          if (bpred[bnode] == INVALID || bvalue < bdist[bnode]) {
     1270            bdist[bnode] = bvalue;
     1271            bpred[bnode] = jt;
     1272          }
     1273          if (bnode_matching[bnode] != INVALID) {
     1274            Node newanode = graph->aNode(bnode_matching[bnode]);
     1275            switch (_heap->state(newanode)) {
     1276            case Heap::PRE_HEAP:
     1277              _heap->push(newanode, bvalue);
     1278              break;
     1279            case Heap::IN_HEAP:
     1280              if (bvalue < (*_heap)[newanode]) {
     1281                _heap->decrease(newanode, bvalue);
     1282              }
     1283              break;
     1284            case Heap::POST_HEAP:
     1285              break;
     1286            }
     1287          } else {
     1288            if (bestNode == INVALID ||
     1289                bvalue + bnode_potential[bnode] < bestValue) {
     1290              bestValue = bvalue + bnode_potential[bnode];
     1291              bestNode = bnode;
     1292            }
     1293          }
     1294        }
     1295      }
     1296
     1297      if (bestNode == INVALID) {
     1298        return false;
     1299      }
     1300
     1301      matching_cost += bestValue;
     1302      ++matching_size;
     1303
     1304      for (BNodeIt it(*graph); it != INVALID; ++it) {
     1305        if (bpred[it] != INVALID) {
     1306          bnode_potential[it] += bdist[it];
     1307        }
     1308      }
     1309      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1310        if (anode_matching[it] != INVALID) {
     1311          Node bnode = graph->bNode(anode_matching[it]);
     1312          if (bpred[bnode] != INVALID) {
     1313            anode_potential[it] += bdist[bnode];
     1314          }
     1315        }
     1316      }
     1317
     1318      while (bestNode != INVALID) {
     1319        UEdge uedge = bpred[bestNode];
     1320        Node anode = graph->aNode(uedge);
     1321       
     1322        bnode_matching[bestNode] = uedge;
     1323        if (anode_matching[anode] != INVALID) {
     1324          bestNode = graph->bNode(anode_matching[anode]);
     1325        } else {
     1326          bestNode = INVALID;
     1327        }
     1328        anode_matching[anode] = uedge;
     1329      }
     1330
     1331
     1332      return true;
     1333    }
     1334
     1335    /// \brief Starts the algorithm.
     1336    ///
     1337    /// Starts the algorithm. It runs augmenting phases until the
     1338    /// optimal solution reached.
     1339    void start() {
     1340      while (augment()) {}
     1341    }
     1342
     1343    /// \brief Runs the algorithm.
     1344    ///
     1345    /// It just initalize the algorithm and then start it.
     1346    void run() {
     1347      init();
     1348      start();
     1349    }
     1350
     1351    /// @}
     1352
     1353    /// \name Query Functions
     1354    /// The result of the %Matching algorithm can be obtained using these
     1355    /// functions.\n
     1356    /// Before the use of these functions,
     1357    /// either run() or start() must be called.
     1358   
     1359    ///@{
     1360
     1361    /// \brief Gives back the potential in the NodeMap
     1362    ///
     1363    /// Gives back the potential in the NodeMap. The potential
     1364    /// is feasible if \f$ \pi(a) - \pi(b) + w(ab) = 0 \f$ for
     1365    /// each matching edges and \f$ \pi(a) - \pi(b) + w(ab) \ge 0 \f$
     1366    /// for each edges.
     1367    template <typename PotentialMap>
     1368    void potential(PotentialMap& potential) {
     1369      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1370        potential[it] = anode_potential[it];
     1371      }
     1372      for (BNodeIt it(*graph); it != INVALID; ++it) {
     1373        potential[it] = bnode_potential[it];
     1374      }
     1375    }
     1376
     1377    /// \brief Set true all matching uedge in the map.
     1378    ///
     1379    /// Set true all matching uedge in the map. It does not change the
     1380    /// value mapped to the other uedges.
     1381    /// \return The number of the matching edges.
     1382    template <typename MatchingMap>
     1383    int quickMatching(MatchingMap& matching) {
     1384      for (ANodeIt it(*graph); it != INVALID; ++it) {
     1385        if (anode_matching[it] != INVALID) {
     1386          matching[anode_matching[it]] = true;
     1387        }
     1388      }
     1389      return matching_size;
     1390    }
     1391
     1392    /// \brief Set true all matching uedge in the map and the others to false.
     1393    ///
     1394    /// Set true all matching uedge in the map and the others to false.
     1395    /// \return The number of the matching edges.
     1396    template <typename MatchingMap>
     1397    int matching(MatchingMap& matching) {
     1398      for (UEdgeIt it(*graph); it != INVALID; ++it) {
     1399        matching[it] = it == anode_matching[graph->aNode(it)];
     1400      }
     1401      return matching_size;
     1402    }
     1403
     1404
     1405    /// \brief Return true if the given uedge is in the matching.
     1406    ///
     1407    /// It returns true if the given uedge is in the matching.
     1408    bool matchingEdge(const UEdge& edge) {
     1409      return anode_matching[graph->aNode(edge)] == edge;
     1410    }
     1411
     1412    /// \brief Returns the matching edge from the node.
     1413    ///
     1414    /// Returns the matching edge from the node. If there is not such
     1415    /// edge it gives back \c INVALID.
     1416    UEdge matchingEdge(const Node& node) {
     1417      if (graph->aNode(node)) {
     1418        return anode_matching[node];
     1419      } else {
     1420        return bnode_matching[node];
     1421      }
     1422    }
     1423
     1424    /// \brief Gives back the sum of costs of the matching edges.
     1425    ///
     1426    /// Gives back the sum of costs of the matching edges.
     1427    Value matchingCost() const {
     1428      return matching_cost;
     1429    }
     1430
     1431    /// \brief Gives back the number of the matching edges.
     1432    ///
     1433    /// Gives back the number of the matching edges.
     1434    int matchingSize() const {
     1435      return matching_size;
     1436    }
     1437
     1438    /// @}
     1439
     1440  private:
     1441
     1442    void initStructures() {
     1443      if (!_heap_cross_ref) {
     1444        local_heap_cross_ref = true;
     1445        _heap_cross_ref = Traits::createHeapCrossRef(*graph);
     1446      }
     1447      if (!_heap) {
     1448        local_heap = true;
     1449        _heap = Traits::createHeap(*_heap_cross_ref);
     1450      }
     1451    }
     1452
     1453    void destroyStructures() {
     1454      if (local_heap_cross_ref) delete _heap_cross_ref;
     1455      if (local_heap) delete _heap;
     1456    }
     1457
     1458
     1459  private:
     1460   
     1461    const BpUGraph *graph;
     1462    const CostMap* cost;
     1463
     1464    ANodeMatchingMap anode_matching;
     1465    BNodeMatchingMap bnode_matching;
     1466
     1467    ANodePotentialMap anode_potential;
     1468    BNodePotentialMap bnode_potential;
     1469
     1470    Value matching_cost;
     1471    int matching_size;
     1472
     1473    HeapCrossRef *_heap_cross_ref;
     1474    bool local_heap_cross_ref;
     1475
     1476    Heap *_heap;
     1477    bool local_heap;
     1478 
     1479  };
     1480
    4641481}
    4651482
Note: See TracChangeset for help on using the changeset viewer.