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