src/libprojectM/HungarianMethod.hpp

#ifndef HUNGARIAN_METHOD_HPP
#define HUNGARIAN_METHOD_HPP
//
//#include "Common.hpp"
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <limits>


/// A function object which calculates the maximum-weighted bipartite matching between
/// two sets via the hungarian method.
template <int N=20>
class HungarianMethod {
public :
static const int MAX_SIZE = N;

private:
int     n, max_match;       // n workers and n jobs
double  lx[N], ly[N];       // labels of X and Y parts
int     xy[N];              // xy[x] - vertex that is matched with x,
int     yx[N];              // yx[y] - vertex that is matched with y
bool    S[N], T[N];         // sets S and T in algorithm
double  slack[N];           // as in the algorithm description
double  slackx[N];          // slackx[y] such a vertex, that
                            //  l(slackx[y]) + l(y) - w(slackx[y],y) = slack[y]
int     prev[N];            // array for memorizing alternating paths

void init_labels(const double cost[N][N])
{
    memset(lx, 0, sizeof(lx));
    memset(ly, 0, sizeof(ly));
    for (int x = 0; x < n; x++)
        for (int y = 0; y < n; y++)
            lx[x] = std::max(lx[x], cost[x][y]);
}

void augment(const double cost[N][N]) //main function of the algorithm
{
    if (max_match == n) return;        //check wether matching is already perfect
    int x, y, root;                    //just counters and root vertex
    int q[N], wr = 0, rd = 0;          //q - queue for bfs, wr,rd - write and read
                                       //pos in queue
    memset(S, false, sizeof(S));       //init set S
    memset(T, false, sizeof(T));       //init set T
    memset(prev, -1, sizeof(prev));    //init set prev - for the alternating tree
    for (x = 0; x < n; x++)            //finding root of the tree
        if (xy[x] == -1)
        {
            q[wr++] = root = x;
            prev[x] = -2;
            S[x] = true;
            break;
        }

    for (y = 0; y < n; y++)            //initializing slack array
    {
        slack[y] = lx[root] + ly[y] - cost[root][y];
        slackx[y] = root;
    }
    while (true)                                                        //main cycle
    {
        while (rd < wr)                                                 //building tree with bfs cycle
        {
            x = q[rd++];                                                //current vertex from X part
            for (y = 0; y < n; y++)                                     //iterate through all edges in equality graph
                if (cost[x][y] == lx[x] + ly[y] &&  !T[y])
                {
                    if (yx[y] == -1) break;                             //an exposed vertex in Y found, so
                                                                        //augmenting path exists!
                    T[y] = true;                                        //else just add y to T,
                    q[wr++] = yx[y];                                    //add vertex yx[y], which is matched
                                                                        //with y, to the queue
                    add_to_tree(yx[y], x, cost);                              //add edges (x,y) and (y,yx[y]) to the tree
                }
            if (y < n) break;                                           //augmenting path found!
        }
        if (y < n) break;                                               //augmenting path found!

        update_labels();                                                //augmenting path not found, so improve labeling
        wr = rd = 0;
        for (y = 0; y < n; y++)
        //in this cycle we add edges that were added to the equality graph as a
        //result of improving the labeling, we add edge (slackx[y], y) to the tree if
        //and only if !T[y] &&  slack[y] == 0, also with this edge we add another one
        //(y, yx[y]) or augment the matching, if y was exposed
            if (!T[y] &&  slack[y] == 0)
            {
                if (yx[y] == -1)                                        //exposed vertex in Y found - augmenting path exists!
                {
                    x = slackx[y];
                    break;
                }
                else
                {
                    T[y] = true;                                        //else just add y to T,
                    if (!S[yx[y]])
                    {
                        q[wr++] = yx[y];                                //add vertex yx[y], which is matched with
                                                                        //y, to the queue
                        add_to_tree(yx[y], slackx[y],cost);                  //and add edges (x,y) and (y,
                                                                        //yx[y]) to the tree
                    }
                }
            }
        if (y < n) break;                                               //augmenting path found!
    }

    if (y < n)                                                          //we found augmenting path!
    {
        max_match++;                                                    //increment matching
        //in this cycle we inverse edges along augmenting path
        for (int cx = x, cy = y, ty; cx != -2; cx = prev[cx], cy = ty)
        {
            ty = xy[cx];
            yx[cy] = cx;
            xy[cx] = cy;
        }
        augment(cost);                                                      //recall function, go to step 1 of the algorithm
    }
}//end of augment() function

void update_labels()
{
    int x, y;
    double delta = std::numeric_limits<double>::max();
    for (y = 0; y < n; y++)            //calculate delta using slack
        if (!T[y])
            delta = std::min(delta, slack[y]);
    for (x = 0; x < n; x++)            //update X labels
        if (S[x]) lx[x] -= delta;
    for (y = 0; y < n; y++)            //update Y labels
        if (T[y]) ly[y] += delta;
    for (y = 0; y < n; y++)            //update slack array
        if (!T[y])
            slack[y] -= delta;
}

void add_to_tree(int x, int prevx, const double cost[N][N])
//x - current vertex,prevx - vertex from X before x in the alternating path,
//so we add edges (prevx, xy[x]), (xy[x], x)
{
    S[x] = true;                    //add x to S
    prev[x] = prevx;                //we need this when augmenting
    for (int y = 0; y < n; y++)    //update slacks, because we add new vertex to S
        if (lx[x] + ly[y] - cost[x][y] < slack[y])
        {
            slack[y] = lx[x] + ly[y] - cost[x][y];
            slackx[y] = x;
        }
}

public:
/// Computes the best matching of two sets given its cost matrix.
/// See the matching() method to get the computed match result.
/// \param cost a matrix of two sets I,J where cost[i][j] is the weight of edge i->j
/// \param logicalSize the number of elements in both I and J
/// \returns the total cost of the best matching
inline double operator()(const double cost[N][N], int logicalSize)
{

    n = logicalSize;
    assert(n <= N);
    double ret = 0;                      //weight of the optimal matching
    max_match = 0;                    //number of vertices in current matching
    memset(xy, -1, sizeof(xy));
    memset(yx, -1, sizeof(yx));
    init_labels(cost);                    //step 0
    augment(cost);                        //steps 1-3
    for (int x = 0; x < n; x++)       //forming answer there
        ret += cost[x][xy[x]];
    return ret;
}

/// Gets the matching element in 2nd set of the ith element in the first set
/// \param i the index of the ith element in the first set (passed in operator())
/// \returns an index j, denoting the matched jth element of the 2nd set
inline int matching(int i) const {
  return xy[i];
}


/// Gets the matching element in 1st set of the jth element in the 2nd set
/// \param j the index of the jth element in the 2nd set (passed in operator())
/// \returns an index i, denoting the matched ith element of the 1st set
/// \note inverseMatching(matching(i)) == i
inline int inverseMatching(int j) const {
  return yx[j];
}

};


#endif