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.
46 #include "hungarian.h"
48 #define INF (0x7FFFFFFF)
50 struct _hungarian_problem_t {
51 int num_rows; /**< number of rows */
52 int num_cols; /**< number of columns */
53 int **cost; /**< the cost matrix */
54 int width; /**< the width for cost matrix dumper */
55 int max_cost; /**< the maximal costs in the matrix */
56 int match_type; /**< PERFECT or NORMAL matching */
57 bitset_t *missing_left; /**< left side nodes having no edge to the right side */
58 bitset_t *missing_right; /**< right side nodes having no edge to the left side */
60 DEBUG_ONLY(firm_dbg_module_t *dbg);
63 static INLINE void *get_init_mem(struct obstack *obst, long sz) {
64 void *p = obstack_alloc(obst, sz);
69 static void hungarian_dump_f(FILE *f, int **C, int rows, int cols, int width) {
73 for (i = 0; i < rows; i++) {
75 for (j = 0; j < cols; j++) {
76 fprintf(f, "%*d", width, C[i][j]);
83 void hungarian_print_costmatrix(hungarian_problem_t *p) {
84 hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, p->width);
88 * Create the object and allocate memory for the data structures.
90 hungarian_problem_t *hungarian_new(int rows, int cols, int width, int match_type) {
92 hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
94 FIRM_DBG_REGISTER(p->dbg, "firm.hungarian");
97 Is the number of cols not equal to number of rows ?
98 If yes, expand with 0 - cols / 0 - cols
100 rows = MAX(cols, rows);
103 obstack_init(&p->obst);
108 p->match_type = match_type;
111 In case of normal matching, we have to keep
112 track of nodes without edges to kill them in
113 the assignment later.
115 if (match_type == HUNGARIAN_MATCH_NORMAL) {
116 p->missing_left = bitset_obstack_alloc(&p->obst, rows);
117 p->missing_right = bitset_obstack_alloc(&p->obst, cols);
118 bitset_set_all(p->missing_left);
119 bitset_set_all(p->missing_right);
122 /* allocate space for cost matrix */
123 p->cost = (int **)get_init_mem(&p->obst, rows * sizeof(p->cost[0]));
124 for (i = 0; i < p->num_rows; i++)
125 p->cost[i] = (int *)get_init_mem(&p->obst, cols * sizeof(p->cost[0][0]));
131 * Prepare the cost matrix.
133 void hungarian_prepare_cost_matrix(hungarian_problem_t *p, int mode) {
136 if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
137 for (i = 0; i < p->num_rows; i++) {
138 for (j = 0; j < p->num_cols; j++) {
139 p->cost[i][j] = p->max_cost - p->cost[i][j];
143 else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
147 fprintf(stderr, "Unknown mode. Mode was set to HUNGARIAN_MODE_MINIMIZE_COST.\n");
151 * Set cost[left][right] to cost.
153 void hungarian_add(hungarian_problem_t *p, int left, int right, int cost) {
154 assert(p->num_rows > left && "Invalid row selected.");
155 assert(p->num_cols > right && "Invalid column selected.");
158 p->cost[left][right] = cost;
159 p->max_cost = MAX(p->max_cost, cost);
161 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
162 bitset_clear(p->missing_left, left);
163 bitset_clear(p->missing_right, right);
168 * Set cost[left][right] to 0.
170 void hungarian_remv(hungarian_problem_t *p, int left, int right) {
171 assert(p->num_rows > left && "Invalid row selected.");
172 assert(p->num_cols > right && "Invalid column selected.");
174 p->cost[left][right] = 0;
176 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
177 bitset_set(p->missing_left, left);
178 bitset_set(p->missing_right, right);
183 * Frees all allocated memory.
185 void hungarian_free(hungarian_problem_t* p) {
186 obstack_free(&p->obst, NULL);
193 int hungarian_solve(hungarian_problem_t* p, int *assignment, int *final_cost, int cost_threshold) {
194 int i, j, m, n, k, l, s, t, q, unmatched, cost;
208 col_mate = XMALLOCNZ(int, p->num_rows);
209 unchosen_row = XMALLOCNZ(int, p->num_rows);
210 row_dec = XMALLOCNZ(int, p->num_rows);
211 slack_row = XMALLOCNZ(int, p->num_rows);
213 row_mate = XMALLOCNZ(int, p->num_cols);
214 parent_row = XMALLOCNZ(int, p->num_cols);
215 col_inc = XMALLOCNZ(int, p->num_cols);
216 slack = XMALLOCNZ(int, p->num_cols);
218 memset(assignment, -1, m * sizeof(assignment[0]));
220 /* Begin subtract column minima in order to start with lots of zeros 12 */
221 DBG((p->dbg, LEVEL_1, "Using heuristic\n"));
223 for (l = 0; l < n; ++l) {
226 for (k = 1; k < m; ++k) {
227 if (p->cost[k][l] < s)
234 for (k = 0; k < m; ++k)
238 /* End subtract column minima in order to start with lots of zeros 12 */
240 /* Begin initial state 16 */
242 for (l = 0; l < n; ++l) {
249 for (k = 0; k < m; ++k) {
252 for (l = 1; l < n; ++l) {
253 if (p->cost[k][l] < s)
259 for (l = 0; l < n; ++l) {
260 if (s == p->cost[k][l] && row_mate[l] < 0) {
263 DBG((p->dbg, LEVEL_1, "matching col %d == row %d\n", l, k));
269 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
270 unchosen_row[t++] = k;
273 /* End initial state 16 */
275 /* Begin Hungarian algorithm 18 */
281 DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", m - t));
286 /* Begin explore node q of the forest 19 */
290 for (l = 0; l < n; ++l) {
292 int del = p->cost[k][l] - s + col_inc[l];
294 if (del < slack[l]) {
301 DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
302 unchosen_row[t++] = row_mate[l];
311 /* End explore node q of the forest 19 */
315 /* Begin introduce a new zero into the matrix 21 */
317 for (l = 0; l < n; ++l) {
318 if (slack[l] && slack[l] < s)
322 for (q = 0; q < t; ++q)
323 row_dec[unchosen_row[q]] += s;
325 for (l = 0; l < n; ++l) {
329 /* Begin look at a new zero 22 */
331 DBG((p->dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%d, %d]\n", s, k, l));
332 if (row_mate[l] < 0) {
333 for (j = l + 1; j < n; ++j) {
341 DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
342 unchosen_row[t++] = row_mate[l];
344 /* End look at a new zero 22 */
351 /* End introduce a new zero into the matrix 21 */
354 /* Begin update the matching 20 */
355 DBG((p->dbg, LEVEL_1, "Breakthrough at node %d of %d.\n", q, t));
361 DBG((p->dbg, LEVEL_1, "rematching col %d == row %d\n", l, k));
368 /* End update the matching 20 */
370 if (--unmatched == 0)
373 /* Begin get ready for another stage 17 */
375 for (l = 0; l < n; ++l) {
380 for (k = 0; k < m; ++k) {
381 if (col_mate[k] < 0) {
382 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
383 unchosen_row[t++] = k;
386 /* End get ready for another stage 17 */
390 /* Begin double check the solution 23 */
391 for (k = 0; k < m; ++k) {
392 for (l = 0; l < n; ++l) {
393 if (p->cost[k][l] < row_dec[k] - col_inc[l])
398 for (k = 0; k < m; ++k) {
400 if (l < 0 || p->cost[k][l] != row_dec[k] - col_inc[l])
404 for (k = l = 0; l < n; ++l) {
411 /* End double check the solution 23 */
413 /* End Hungarian algorithm 18 */
415 /* collect the assigned values */
416 for (i = 0; i < m; ++i) {
417 if (cost_threshold > 0 && p->cost[i][col_mate[i]] >= cost_threshold)
418 assignment[i] = -1; /* remove matching having cost > threshold */
420 assignment[i] = col_mate[i];
423 /* In case of normal matching: remove impossible ones */
424 if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
425 for (i = 0; i < m; ++i) {
426 if (bitset_is_set(p->missing_left, i) || bitset_is_set(p->missing_right, col_mate[i]))
431 for (k = 0; k < m; ++k) {
432 for (l = 0; l < n; ++l) {
433 p->cost[k][l] = p->cost[k][l] - row_dec[k] + col_inc[l];
437 for (i = 0; i < m; ++i)
440 for (i = 0; i < n; ++i)
443 DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));