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.
45 #include "hungarian.h"
47 struct hungarian_problem_t {
48 unsigned num_rows; /**< number of rows */
49 unsigned num_cols; /**< number of columns */
50 int **cost; /**< the cost matrix */
51 int max_cost; /**< the maximal costs in the matrix */
52 int match_type; /**< PERFECT or NORMAL matching */
53 bitset_t *missing_left; /**< left side nodes having no edge to the right side */
54 bitset_t *missing_right; /**< right side nodes having no edge to the left side */
56 DEBUG_ONLY(firm_dbg_module_t *dbg);
59 static void hungarian_dump_f(FILE *f, int **C, unsigned n_rows, unsigned n_cols,
65 for (r = 0; r < n_rows; r++) {
68 for (c = 0; c < n_cols; c++) {
69 fprintf(f, "%*d", width, C[r][c]);
76 void hungarian_print_cost_matrix(hungarian_problem_t *p, int width)
78 hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, width);
81 hungarian_problem_t *hungarian_new(unsigned n_rows, unsigned n_cols,
82 match_type_t match_type)
85 hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
87 FIRM_DBG_REGISTER(p->dbg, "firm.hungarian");
90 Is the number of cols not equal to number of rows ?
91 If yes, expand with 0 - cols / 0 - cols
93 n_rows = MAX(n_cols, n_rows);
96 obstack_init(&p->obst);
100 p->match_type = match_type;
103 In case of normal matching, we have to keep
104 track of nodes without edges to kill them in
105 the assignment later.
107 if (match_type == HUNGARIAN_MATCH_NORMAL) {
108 p->missing_left = bitset_obstack_alloc(&p->obst, n_rows);
109 p->missing_right = bitset_obstack_alloc(&p->obst, n_cols);
110 bitset_set_all(p->missing_left);
111 bitset_set_all(p->missing_right);
114 /* allocate space for cost matrix */
115 p->cost = OALLOCNZ(&p->obst, int*, n_rows);
116 for (r = 0; r < p->num_rows; r++)
117 p->cost[r] = OALLOCNZ(&p->obst, int, n_cols);
122 void hungarian_prepare_cost_matrix(hungarian_problem_t *p,
123 hungarian_mode_t mode)
127 if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
128 for (r = 0; r < p->num_rows; r++) {
129 for (c = 0; c < p->num_cols; c++) {
130 p->cost[r][c] = p->max_cost - p->cost[r][c];
133 } else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
136 panic("Unknown hungarian problem mode\n");
140 void hungarian_add(hungarian_problem_t *p, unsigned left, unsigned right,
143 assert(p->num_rows > left && "Invalid row selected.");
144 assert(p->num_cols > right && "Invalid column selected.");
147 p->cost[left][right] = cost;
148 p->max_cost = MAX(p->max_cost, cost);
150 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
151 bitset_clear(p->missing_left, left);
152 bitset_clear(p->missing_right, right);
156 void hungarian_remove(hungarian_problem_t *p, unsigned left, unsigned right)
158 assert(p->num_rows > left && "Invalid row selected.");
159 assert(p->num_cols > right && "Invalid column selected.");
161 /* Set cost[left][right] to 0. */
162 p->cost[left][right] = 0;
164 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
165 bitset_set(p->missing_left, left);
166 bitset_set(p->missing_right, right);
170 void hungarian_free(hungarian_problem_t* p)
172 obstack_free(&p->obst, NULL);
176 int hungarian_solve(hungarian_problem_t* p, unsigned *assignment,
177 int *final_cost, int cost_threshold)
180 unsigned num_rows = p->num_rows;
181 unsigned num_cols = p->num_cols;
182 unsigned *col_mate = XMALLOCNZ(unsigned, num_rows);
183 unsigned *row_mate = XMALLOCNZ(unsigned, num_cols);
184 unsigned *parent_row = XMALLOCNZ(unsigned, num_cols);
185 unsigned *unchosen_row = XMALLOCNZ(unsigned, num_rows);
186 int *row_dec = XMALLOCNZ(int, num_rows);
187 int *col_inc = XMALLOCNZ(int, num_cols);
188 int *slack = XMALLOCNZ(int, num_cols);
189 unsigned *slack_row = XMALLOCNZ(unsigned, num_rows);
195 memset(assignment, -1, num_rows * sizeof(assignment[0]));
197 /* Begin subtract column minima in order to start with lots of zeros 12 */
198 DBG((p->dbg, LEVEL_1, "Using heuristic\n"));
200 for (c = 0; c < num_cols; ++c) {
201 int s = p->cost[0][c];
203 for (r = 1; r < num_rows; ++r) {
204 if (p->cost[r][c] < s)
212 for (r = 0; r < num_rows; ++r)
215 /* End subtract column minima in order to start with lots of zeros 12 */
217 /* Begin initial state 16 */
219 for (c = 0; c < num_cols; ++c) {
220 row_mate[c] = (unsigned) -1;
221 parent_row[c] = (unsigned) -1;
226 for (r = 0; r < num_rows; ++r) {
227 int s = p->cost[r][0];
229 for (c = 1; c < num_cols; ++c) {
230 if (p->cost[r][c] < s)
236 for (c = 0; c < num_cols; ++c) {
237 if (s == p->cost[r][c] && row_mate[c] == (unsigned)-1) {
240 DBG((p->dbg, LEVEL_1, "matching col %d == row %d\n", c, r));
245 col_mate[r] = (unsigned)-1;
246 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, r));
247 unchosen_row[t++] = r;
250 /* End initial state 16 */
252 /* Begin Hungarian algorithm 18 */
260 DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", num_rows - t));
265 /* Begin explore node q of the forest 19 */
269 for (c = 0; c < num_cols; ++c) {
271 int del = p->cost[r][c] - s + col_inc[c];
273 if (del < slack[c]) {
275 if (row_mate[c] == (unsigned)-1)
280 DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[c], c, r));
281 unchosen_row[t++] = row_mate[c];
289 /* End explore node q of the forest 19 */
293 /* Begin introduce a new zero into the matrix 21 */
295 for (c = 0; c < num_cols; ++c) {
296 if (slack[c] && slack[c] < s)
300 for (q = 0; q < t; ++q)
301 row_dec[unchosen_row[q]] += s;
303 for (c = 0; c < num_cols; ++c) {
307 /* Begin look at a new zero 22 */
309 DBG((p->dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%d, %d]\n", s, r, c));
310 if (row_mate[c] == (unsigned)-1) {
311 for (j = c + 1; j < num_cols; ++j) {
318 DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[c], c, r));
319 unchosen_row[t++] = row_mate[c];
321 /* End look at a new zero 22 */
327 /* End introduce a new zero into the matrix 21 */
330 /* Begin update the matching 20 */
331 DBG((p->dbg, LEVEL_1, "Breakthrough at node %d of %d.\n", q, t));
337 DBG((p->dbg, LEVEL_1, "rematching col %d == row %d\n", c, r));
338 if (j == (unsigned)-1)
344 /* End update the matching 20 */
346 if (--unmatched == 0)
349 /* Begin get ready for another stage 17 */
351 for (c = 0; c < num_cols; ++c) {
356 for (r = 0; r < num_rows; ++r) {
357 if (col_mate[r] == (unsigned)-1) {
358 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, r));
359 unchosen_row[t++] = r;
362 /* End get ready for another stage 17 */
366 /* Begin double check the solution 23 */
367 for (r = 0; r < num_rows; ++r) {
368 for (c = 0; c < num_cols; ++c) {
369 if (p->cost[r][c] < row_dec[r] - col_inc[c])
374 for (r = 0; r < num_rows; ++r) {
376 if (c == (unsigned)-1 || p->cost[r][c] != 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 && p->cost[r][col_mate[r]] >= cost_threshold)
394 assignment[r] = -1; /* remove matching having cost > threshold */
396 assignment[r] = col_mate[r];
399 /* In case of normal matching: remove impossible ones */
400 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
401 for (r = 0; r < num_rows; ++r) {
402 if (bitset_is_set(p->missing_left, r)
403 || bitset_is_set(p->missing_right, col_mate[r]))
408 for (r = 0; r < num_rows; ++r) {
409 for (c = 0; c < num_cols; ++c) {
410 p->cost[r][c] = p->cost[r][c] - row_dec[r] + col_inc[c];
414 for (r = 0; r < num_rows; ++r)
417 for (c = 0; c < num_cols; ++c)
420 DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));
431 if (final_cost != NULL)