1 /********************************************************************
2 ********************************************************************
4 ** libhungarian by Cyrill Stachniss, 2004
6 ** Added and adapted to libFirm by Christian Wuerdig, 2006
8 ** Solving the Minimum Assignment Problem using the
11 ** ** This file may be freely copied and distributed! **
13 ** Parts of the used code was originally provided by the
14 ** "Stanford GraphGase", but I made changes to this code.
15 ** As asked by the copyright node of the "Stanford GraphGase",
16 ** I hereby proclaim that this file are *NOT* part of the
17 ** "Stanford GraphGase" distrubition!
19 ** This file is distributed in the hope that it will be useful,
20 ** but WITHOUT ANY WARRANTY; without even the implied
21 ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
24 ********************************************************************
25 ********************************************************************/
29 * @brief Solving the Minimum Assignment Problem using the Hungarian Method.
43 #include "hungarian.h"
45 DEBUG_ONLY(static firm_dbg_module_t *dbg;)
47 struct hungarian_problem_t {
48 unsigned num_rows; /**< number of rows */
49 unsigned num_cols; /**< number of columns */
50 unsigned *cost; /**< the cost matrix */
51 unsigned max_cost; /**< the maximal costs in the matrix */
52 match_type_t match_type; /**< PERFECT or NORMAL matching */
53 unsigned *missing_left; /**< bitset: left side nodes having no edge to
55 unsigned *missing_right; /**< bitset: right side nodes having no edge to
59 static void hungarian_dump_f(FILE *f, const unsigned *cost,
60 unsigned num_rows, unsigned num_cols, int width)
65 for (r = 0; r < num_rows; r++) {
67 for (c = 0; c < num_cols; c++) {
68 fprintf(f, "%*u", width, cost[r*num_cols + c]);
75 void hungarian_print_cost_matrix(hungarian_problem_t *p, int width)
77 hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, width);
80 hungarian_problem_t *hungarian_new(unsigned num_rows, unsigned num_cols,
81 match_type_t match_type)
83 hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
85 FIRM_DBG_REGISTER(dbg, "firm.hungarian");
88 Is the number of cols not equal to number of rows ?
89 If yes, expand with 0 - cols / 0 - cols
91 num_rows = MAX(num_cols, num_rows);
94 p->num_rows = num_rows;
95 p->num_cols = num_cols;
96 p->match_type = match_type;
99 In case of normal matching, we have to keep
100 track of nodes without edges to kill them in
101 the assignment later.
103 if (match_type == HUNGARIAN_MATCH_NORMAL) {
104 p->missing_left = rbitset_malloc(num_rows);
105 p->missing_right = rbitset_malloc(num_cols);
106 rbitset_set_all(p->missing_left, num_rows);
107 rbitset_set_all(p->missing_right, num_cols);
110 /* allocate space for cost matrix */
111 p->cost = XMALLOCNZ(unsigned, num_rows * num_cols);
115 void hungarian_prepare_cost_matrix(hungarian_problem_t *p,
116 hungarian_mode_t mode)
118 if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
120 unsigned num_cols = p->num_cols;
121 unsigned *cost = p->cost;
122 unsigned max_cost = p->max_cost;
123 for (r = 0; r < p->num_rows; r++) {
124 for (c = 0; c < p->num_cols; c++) {
125 cost[r*num_cols + c] = max_cost - cost[r*num_cols + c];
128 } else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
131 panic("Unknown hungarian problem mode");
135 void hungarian_add(hungarian_problem_t *p, unsigned left, unsigned right,
138 assert(p->num_rows > left && "Invalid row selected.");
139 assert(p->num_cols > right && "Invalid column selected.");
141 p->cost[left*p->num_cols + right] = cost;
142 p->max_cost = MAX(p->max_cost, cost);
144 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
145 rbitset_clear(p->missing_left, left);
146 rbitset_clear(p->missing_right, right);
150 void hungarian_remove(hungarian_problem_t *p, unsigned left, unsigned right)
152 assert(p->num_rows > left && "Invalid row selected.");
153 assert(p->num_cols > right && "Invalid column selected.");
155 p->cost[left*p->num_cols + right] = 0;
157 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
158 rbitset_set(p->missing_left, left);
159 rbitset_set(p->missing_right, right);
163 void hungarian_free(hungarian_problem_t* p)
165 xfree(p->missing_left);
166 xfree(p->missing_right);
171 int hungarian_solve(hungarian_problem_t* p, unsigned *assignment,
172 unsigned *final_cost, unsigned cost_threshold)
174 unsigned res_cost = 0;
175 unsigned num_rows = p->num_rows;
176 unsigned num_cols = p->num_cols;
177 unsigned *cost = p->cost;
178 unsigned *col_mate = XMALLOCNZ(unsigned, num_rows);
179 unsigned *row_mate = XMALLOCNZ(unsigned, num_cols);
180 unsigned *parent_row = XMALLOCNZ(unsigned, num_cols);
181 unsigned *unchosen_row = XMALLOCNZ(unsigned, num_rows);
182 int *row_dec = XMALLOCNZ(int, num_rows);
183 int *col_inc = XMALLOCNZ(int, num_cols);
184 int *slack = XMALLOCNZ(int, num_cols);
185 unsigned *slack_row = XMALLOCNZ(unsigned, num_rows);
191 memset(assignment, -1, num_rows * sizeof(assignment[0]));
193 /* Begin subtract column minima in order to start with lots of zeros 12 */
194 DBG((dbg, LEVEL_1, "Using heuristic\n"));
196 for (c = 0; c < num_cols; ++c) {
197 unsigned col_mininum = cost[0*num_cols + c];
199 for (r = 1; r < num_rows; ++r) {
200 if (cost[r*num_cols + c] < col_mininum)
201 col_mininum = cost[r*num_cols + c];
204 if (col_mininum == 0)
207 res_cost += col_mininum;
208 for (r = 0; r < num_rows; ++r)
209 cost[r*num_cols + c] -= col_mininum;
211 /* End subtract column minima in order to start with lots of zeros 12 */
213 /* Begin initial state 16 */
215 for (c = 0; c < num_cols; ++c) {
216 row_mate[c] = (unsigned) -1;
217 parent_row[c] = (unsigned) -1;
222 for (r = 0; r < num_rows; ++r) {
223 unsigned row_minimum = cost[r*num_cols + 0];
225 for (c = 1; c < num_cols; ++c) {
226 if (cost[r*num_cols + c] < row_minimum)
227 row_minimum = cost[r*num_cols + c];
230 row_dec[r] = row_minimum;
232 for (c = 0; c < num_cols; ++c) {
233 if (cost[r*num_cols + c] != row_minimum)
235 if (row_mate[c] != (unsigned)-1)
240 DBG((dbg, LEVEL_1, "matching col %u == row %u\n", c, r));
244 col_mate[r] = (unsigned)-1;
245 DBG((dbg, LEVEL_1, "node %u: unmatched row %u\n", unmatched, r));
246 unchosen_row[unmatched++] = r;
249 /* End initial state 16 */
251 /* Begin Hungarian algorithm 18 */
259 DBG((dbg, LEVEL_1, "Matched %u rows.\n", num_rows - t));
264 /* Begin explore node q of the forest 19 */
268 for (c = 0; c < num_cols; ++c) {
270 int del = cost[r*num_cols + c] - s + col_inc[c];
272 if (del < slack[c]) {
274 if (row_mate[c] == (unsigned)-1)
279 DBG((dbg, LEVEL_1, "node %u: row %u == col %u -- row %u\n", t, row_mate[c], c, r));
280 unchosen_row[t++] = row_mate[c];
288 /* End explore node q of the forest 19 */
292 /* Begin introduce a new zero into the matrix 21 */
294 for (c = 0; c < num_cols; ++c) {
295 if (slack[c] && slack[c] < s)
299 for (q = 0; q < t; ++q)
300 row_dec[unchosen_row[q]] += s;
302 for (c = 0; c < num_cols; ++c) {
306 /* Begin look at a new zero 22 */
308 DBG((dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%u, %u]\n", s, r, c));
309 if (row_mate[c] == (unsigned)-1) {
310 for (j = c + 1; j < num_cols; ++j) {
317 DBG((dbg, LEVEL_1, "node %u: row %u == col %u -- row %u\n", t, row_mate[c], c, r));
318 unchosen_row[t++] = row_mate[c];
320 /* End look at a new zero 22 */
326 /* End introduce a new zero into the matrix 21 */
329 /* Begin update the matching 20 */
330 DBG((dbg, LEVEL_1, "Breakthrough at node %u of %u.\n", q, t));
336 DBG((dbg, LEVEL_1, "rematching col %u == row %u\n", c, r));
337 if (j == (unsigned)-1)
343 /* End update the matching 20 */
345 if (--unmatched == 0)
348 /* Begin get ready for another stage 17 */
350 for (c = 0; c < num_cols; ++c) {
351 parent_row[c] = (unsigned) -1;
355 for (r = 0; r < num_rows; ++r) {
356 if (col_mate[r] == (unsigned)-1) {
357 DBG((dbg, LEVEL_1, "node %u: unmatched row %u\n", t, r));
358 unchosen_row[t++] = r;
361 /* End get ready for another stage 17 */
365 /* Begin double check the solution 23 */
366 for (r = 0; r < num_rows; ++r) {
367 for (c = 0; c < num_cols; ++c) {
368 if ((int) cost[r*num_cols + c] < row_dec[r] - col_inc[c])
373 for (r = 0; r < num_rows; ++r) {
375 if (c == (unsigned)-1
376 || cost[r*num_cols + c] != (unsigned) (row_dec[r] - col_inc[c]))
380 for (r = c = 0; c < num_cols; ++c) {
387 /* End double check the solution 23 */
389 /* End Hungarian algorithm 18 */
391 /* collect the assigned values */
392 for (r = 0; r < num_rows; ++r) {
393 if (cost_threshold > 0
394 && cost[r*num_cols + col_mate[r]] >= cost_threshold)
395 assignment[r] = -1; /* remove matching having cost > threshold */
397 assignment[r] = col_mate[r];
400 /* In case of normal matching: remove impossible ones */
401 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
402 for (r = 0; r < num_rows; ++r) {
403 if (rbitset_is_set(p->missing_left, r)
404 || rbitset_is_set(p->missing_right, col_mate[r]))
409 for (r = 0; r < num_rows; ++r) {
410 for (c = 0; c < num_cols; ++c) {
411 cost[r*num_cols + c] = cost[r*num_cols + c] - row_dec[r] + col_inc[c];
415 for (r = 0; r < num_rows; ++r)
416 res_cost += row_dec[r];
418 for (c = 0; c < num_cols; ++c)
419 res_cost -= col_inc[c];
421 DBG((dbg, LEVEL_1, "Cost is %d\n", res_cost));
432 if (final_cost != NULL)
433 *final_cost = res_cost;