/**
* @file
* @brief Solving the Minimum Assignment Problem using the Hungarian Method.
- * @version $Id$
*/
#include "config.h"
#include <stdlib.h>
#include <assert.h>
-#include "irtools.h"
+#include "util.h"
#include "xmalloc.h"
#include "debug.h"
-#include "obst.h"
#include "bitset.h"
+#include "error.h"
#include "hungarian.h"
-#define INF (0x7FFFFFFF)
-
-struct _hungarian_problem_t {
- int num_rows; /**< number of rows */
- int num_cols; /**< number of columns */
- int **cost; /**< the cost matrix */
- int width; /**< the width for cost matrix dumper */
- int max_cost; /**< the maximal costs in the matrix */
- int match_type; /**< PERFECT or NORMAL matching */
- bitset_t *missing_left; /**< left side nodes having no edge to the right side */
- bitset_t *missing_right; /**< right side nodes having no edge to the left side */
- struct obstack obst;
- DEBUG_ONLY(firm_dbg_module_t *dbg);
+DEBUG_ONLY(static firm_dbg_module_t *dbg;)
+
+struct hungarian_problem_t {
+ unsigned num_rows; /**< number of rows */
+ unsigned num_cols; /**< number of columns */
+ unsigned *cost; /**< the cost matrix */
+ unsigned max_cost; /**< the maximal costs in the matrix */
+ match_type_t match_type; /**< PERFECT or NORMAL matching */
+ unsigned *missing_left; /**< bitset: left side nodes having no edge to
+ the right side */
+ unsigned *missing_right; /**< bitset: right side nodes having no edge to
+ the left side */
};
-static inline void *get_init_mem(struct obstack *obst, long sz) {
- void *p = obstack_alloc(obst, sz);
- memset(p, 0, sz);
- return p;
-}
-
-static void hungarian_dump_f(FILE *f, int **C, int rows, int cols, int width) {
- int i, j;
+static void hungarian_dump_f(FILE *f, const unsigned *cost,
+ unsigned num_rows, unsigned num_cols, int width)
+{
+ unsigned r, c;
fprintf(f , "\n");
- for (i = 0; i < rows; i++) {
+ for (r = 0; r < num_rows; r++) {
fprintf(f, " [");
- for (j = 0; j < cols; j++) {
- fprintf(f, "%*d", width, C[i][j]);
+ for (c = 0; c < num_cols; c++) {
+ fprintf(f, "%*u", width, cost[r*num_cols + c]);
}
fprintf(f, "]\n");
}
fprintf(f, "\n");
}
-void hungarian_print_costmatrix(hungarian_problem_t *p) {
- hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, p->width);
+void hungarian_print_cost_matrix(hungarian_problem_t *p, int width)
+{
+ hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, width);
}
-/**
- * Create the object and allocate memory for the data structures.
- */
-hungarian_problem_t *hungarian_new(int rows, int cols, int width, int match_type) {
- int i;
+hungarian_problem_t *hungarian_new(unsigned num_rows, unsigned num_cols,
+ match_type_t match_type)
+{
hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
- FIRM_DBG_REGISTER(p->dbg, "firm.hungarian");
+ FIRM_DBG_REGISTER(dbg, "firm.hungarian");
/*
Is the number of cols not equal to number of rows ?
If yes, expand with 0 - cols / 0 - cols
*/
- rows = MAX(cols, rows);
- cols = rows;
-
- obstack_init(&p->obst);
+ num_rows = MAX(num_cols, num_rows);
+ num_cols = num_rows;
- p->num_rows = rows;
- p->num_cols = cols;
- p->width = width;
+ p->num_rows = num_rows;
+ p->num_cols = num_cols;
p->match_type = match_type;
/*
the assignment later.
*/
if (match_type == HUNGARIAN_MATCH_NORMAL) {
- p->missing_left = bitset_obstack_alloc(&p->obst, rows);
- p->missing_right = bitset_obstack_alloc(&p->obst, cols);
- bitset_set_all(p->missing_left);
- bitset_set_all(p->missing_right);
+ p->missing_left = rbitset_malloc(num_rows);
+ p->missing_right = rbitset_malloc(num_cols);
+ rbitset_set_all(p->missing_left, num_rows);
+ rbitset_set_all(p->missing_right, num_cols);
}
/* allocate space for cost matrix */
- p->cost = (int **)get_init_mem(&p->obst, rows * sizeof(p->cost[0]));
- for (i = 0; i < p->num_rows; i++)
- p->cost[i] = (int *)get_init_mem(&p->obst, cols * sizeof(p->cost[0][0]));
-
+ p->cost = XMALLOCNZ(unsigned, num_rows * num_cols);
return p;
}
-/**
- * Prepare the cost matrix.
- */
-void hungarian_prepare_cost_matrix(hungarian_problem_t *p, int mode) {
- int i, j;
-
+void hungarian_prepare_cost_matrix(hungarian_problem_t *p,
+ hungarian_mode_t mode)
+{
if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
- for (i = 0; i < p->num_rows; i++) {
- for (j = 0; j < p->num_cols; j++) {
- p->cost[i][j] = p->max_cost - p->cost[i][j];
+ unsigned r, c;
+ unsigned num_cols = p->num_cols;
+ unsigned *cost = p->cost;
+ unsigned max_cost = p->max_cost;
+ for (r = 0; r < p->num_rows; r++) {
+ for (c = 0; c < p->num_cols; c++) {
+ cost[r*num_cols + c] = max_cost - cost[r*num_cols + c];
}
}
- }
- else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
+ } else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
/* nothing to do */
+ } else {
+ panic("Unknown hungarian problem mode");
}
- else
- fprintf(stderr, "Unknown mode. Mode was set to HUNGARIAN_MODE_MINIMIZE_COST.\n");
}
-/**
- * Set cost[left][right] to cost.
- */
-void hungarian_add(hungarian_problem_t *p, int left, int right, int cost) {
+void hungarian_add(hungarian_problem_t *p, unsigned left, unsigned right,
+ unsigned cost)
+{
assert(p->num_rows > left && "Invalid row selected.");
assert(p->num_cols > right && "Invalid column selected.");
- assert(cost >= 0);
- p->cost[left][right] = cost;
- p->max_cost = MAX(p->max_cost, cost);
+ p->cost[left*p->num_cols + right] = cost;
+ p->max_cost = MAX(p->max_cost, cost);
if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
- bitset_clear(p->missing_left, left);
- bitset_clear(p->missing_right, right);
+ rbitset_clear(p->missing_left, left);
+ rbitset_clear(p->missing_right, right);
}
}
-/**
- * Set cost[left][right] to 0.
- */
-void hungarian_remv(hungarian_problem_t *p, int left, int right) {
+void hungarian_remove(hungarian_problem_t *p, unsigned left, unsigned right)
+{
assert(p->num_rows > left && "Invalid row selected.");
assert(p->num_cols > right && "Invalid column selected.");
- p->cost[left][right] = 0;
+ p->cost[left*p->num_cols + right] = 0;
if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
- bitset_set(p->missing_left, left);
- bitset_set(p->missing_right, right);
+ rbitset_set(p->missing_left, left);
+ rbitset_set(p->missing_right, right);
}
}
-/**
- * Frees all allocated memory.
- */
-void hungarian_free(hungarian_problem_t* p) {
- obstack_free(&p->obst, NULL);
+void hungarian_free(hungarian_problem_t* p)
+{
+ xfree(p->missing_left);
+ xfree(p->missing_right);
+ xfree(p->cost);
xfree(p);
}
-/**
- * Do the assignment.
- */
-int hungarian_solve(hungarian_problem_t* p, int *assignment, int *final_cost, int cost_threshold) {
- int i, j, m, n, k, l, s, t, q, unmatched, cost;
- int *col_mate;
- int *row_mate;
- int *parent_row;
- int *unchosen_row;
- int *row_dec;
- int *col_inc;
- int *slack;
- int *slack_row;
-
- cost = 0;
- m = p->num_rows;
- n = p->num_cols;
-
- col_mate = XMALLOCNZ(int, p->num_rows);
- unchosen_row = XMALLOCNZ(int, p->num_rows);
- row_dec = XMALLOCNZ(int, p->num_rows);
- slack_row = XMALLOCNZ(int, p->num_rows);
-
- row_mate = XMALLOCNZ(int, p->num_cols);
- parent_row = XMALLOCNZ(int, p->num_cols);
- col_inc = XMALLOCNZ(int, p->num_cols);
- slack = XMALLOCNZ(int, p->num_cols);
-
- memset(assignment, -1, m * sizeof(assignment[0]));
+int hungarian_solve(hungarian_problem_t* p, unsigned *assignment,
+ unsigned *final_cost, unsigned cost_threshold)
+{
+ unsigned res_cost = 0;
+ unsigned num_rows = p->num_rows;
+ unsigned num_cols = p->num_cols;
+ unsigned *cost = p->cost;
+ unsigned *col_mate = XMALLOCNZ(unsigned, num_rows);
+ unsigned *row_mate = XMALLOCNZ(unsigned, num_cols);
+ unsigned *parent_row = XMALLOCNZ(unsigned, num_cols);
+ unsigned *unchosen_row = XMALLOCNZ(unsigned, num_rows);
+ int *row_dec = XMALLOCNZ(int, num_rows);
+ int *col_inc = XMALLOCNZ(int, num_cols);
+ int *slack = XMALLOCNZ(int, num_cols);
+ unsigned *slack_row = XMALLOCNZ(unsigned, num_rows);
+ unsigned r;
+ unsigned c;
+ unsigned t;
+ unsigned unmatched;
+
+ memset(assignment, -1, num_rows * sizeof(assignment[0]));
/* Begin subtract column minima in order to start with lots of zeros 12 */
- DBG((p->dbg, LEVEL_1, "Using heuristic\n"));
+ DBG((dbg, LEVEL_1, "Using heuristic\n"));
- for (l = 0; l < n; ++l) {
- s = p->cost[0][l];
+ for (c = 0; c < num_cols; ++c) {
+ unsigned col_mininum = cost[0*num_cols + c];
- for (k = 1; k < m; ++k) {
- if (p->cost[k][l] < s)
- s = p->cost[k][l];
+ for (r = 1; r < num_rows; ++r) {
+ if (cost[r*num_cols + c] < col_mininum)
+ col_mininum = cost[r*num_cols + c];
}
- cost += s;
+ if (col_mininum == 0)
+ continue;
- if (s != 0) {
- for (k = 0; k < m; ++k)
- p->cost[k][l] -= s;
- }
+ res_cost += col_mininum;
+ for (r = 0; r < num_rows; ++r)
+ cost[r*num_cols + c] -= col_mininum;
}
/* End subtract column minima in order to start with lots of zeros 12 */
/* Begin initial state 16 */
- t = 0;
- for (l = 0; l < n; ++l) {
- row_mate[l] = -1;
- parent_row[l] = -1;
- col_inc[l] = 0;
- slack[l] = INF;
+ unmatched = 0;
+ for (c = 0; c < num_cols; ++c) {
+ row_mate[c] = (unsigned) -1;
+ parent_row[c] = (unsigned) -1;
+ col_inc[c] = 0;
+ slack[c] = INT_MAX;
}
- for (k = 0; k < m; ++k) {
- s = p->cost[k][0];
+ for (r = 0; r < num_rows; ++r) {
+ unsigned row_minimum = cost[r*num_cols + 0];
- for (l = 1; l < n; ++l) {
- if (p->cost[k][l] < s)
- s = p->cost[k][l];
+ for (c = 1; c < num_cols; ++c) {
+ if (cost[r*num_cols + c] < row_minimum)
+ row_minimum = cost[r*num_cols + c];
}
- row_dec[k] = s;
+ row_dec[r] = row_minimum;
- for (l = 0; l < n; ++l) {
- if (s == p->cost[k][l] && row_mate[l] < 0) {
- col_mate[k] = l;
- row_mate[l] = k;
- DBG((p->dbg, LEVEL_1, "matching col %d == row %d\n", l, k));
- goto row_done;
- }
+ for (c = 0; c < num_cols; ++c) {
+ if (cost[r*num_cols + c] != row_minimum)
+ continue;
+ if (row_mate[c] != (unsigned)-1)
+ continue;
+
+ col_mate[r] = c;
+ row_mate[c] = r;
+ DBG((dbg, LEVEL_1, "matching col %u == row %u\n", c, r));
+ goto row_done;
}
- col_mate[k] = -1;
- DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
- unchosen_row[t++] = k;
+ col_mate[r] = (unsigned)-1;
+ DBG((dbg, LEVEL_1, "node %u: unmatched row %u\n", unmatched, r));
+ unchosen_row[unmatched++] = r;
row_done: ;
}
/* End initial state 16 */
/* Begin Hungarian algorithm 18 */
- if (t == 0)
+ if (unmatched == 0)
goto done;
- unmatched = t;
- while (1) {
- DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", m - t));
- q = 0;
+ t = unmatched;
+ for (;;) {
+ unsigned q = 0;
+ unsigned j;
+ DBG((dbg, LEVEL_1, "Matched %u rows.\n", num_rows - t));
- while (1) {
+ for (;;) {
+ int s;
while (q < t) {
/* Begin explore node q of the forest 19 */
- k = unchosen_row[q];
- s = row_dec[k];
+ r = unchosen_row[q];
+ s = row_dec[r];
- for (l = 0; l < n; ++l) {
- if (slack[l]) {
- int del = p->cost[k][l] - s + col_inc[l];
+ for (c = 0; c < num_cols; ++c) {
+ if (slack[c]) {
+ int del = cost[r*num_cols + c] - s + col_inc[c];
- if (del < slack[l]) {
+ if (del < slack[c]) {
if (del == 0) {
- if (row_mate[l] < 0)
+ if (row_mate[c] == (unsigned)-1)
goto breakthru;
- slack[l] = 0;
- parent_row[l] = k;
- DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
- unchosen_row[t++] = row_mate[l];
- }
- else {
- slack[l] = del;
- slack_row[l] = k;
+ slack[c] = 0;
+ parent_row[c] = r;
+ DBG((dbg, LEVEL_1, "node %u: row %u == col %u -- row %u\n", t, row_mate[c], c, r));
+ unchosen_row[t++] = row_mate[c];
+ } else {
+ slack[c] = del;
+ slack_row[c] = r;
}
}
}
}
/* Begin introduce a new zero into the matrix 21 */
- s = INF;
- for (l = 0; l < n; ++l) {
- if (slack[l] && slack[l] < s)
- s = slack[l];
+ s = INT_MAX;
+ for (c = 0; c < num_cols; ++c) {
+ if (slack[c] && slack[c] < s)
+ s = slack[c];
}
for (q = 0; q < t; ++q)
row_dec[unchosen_row[q]] += s;
- for (l = 0; l < n; ++l) {
- if (slack[l]) {
- slack[l] -= s;
- if (slack[l] == 0) {
+ for (c = 0; c < num_cols; ++c) {
+ if (slack[c]) {
+ slack[c] -= s;
+ if (slack[c] == 0) {
/* Begin look at a new zero 22 */
- k = slack_row[l];
- DBG((p->dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%d, %d]\n", s, k, l));
- if (row_mate[l] < 0) {
- for (j = l + 1; j < n; ++j) {
+ r = slack_row[c];
+ DBG((dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%u, %u]\n", s, r, c));
+ if (row_mate[c] == (unsigned)-1) {
+ for (j = c + 1; j < num_cols; ++j) {
if (slack[j] == 0)
col_inc[j] += s;
}
goto breakthru;
- }
- else {
- parent_row[l] = k;
- DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
- unchosen_row[t++] = row_mate[l];
+ } else {
+ parent_row[c] = r;
+ DBG((dbg, LEVEL_1, "node %u: row %u == col %u -- row %u\n", t, row_mate[c], c, r));
+ unchosen_row[t++] = row_mate[c];
}
/* End look at a new zero 22 */
}
- }
- else {
- col_inc[l] += s;
+ } else {
+ col_inc[c] += s;
}
}
/* End introduce a new zero into the matrix 21 */
}
breakthru:
/* Begin update the matching 20 */
- DBG((p->dbg, LEVEL_1, "Breakthrough at node %d of %d.\n", q, t));
- while (1) {
- j = col_mate[k];
- col_mate[k] = l;
- row_mate[l] = k;
-
- DBG((p->dbg, LEVEL_1, "rematching col %d == row %d\n", l, k));
- if (j < 0)
+ DBG((dbg, LEVEL_1, "Breakthrough at node %u of %u.\n", q, t));
+ for (;;) {
+ j = col_mate[r];
+ col_mate[r] = c;
+ row_mate[c] = r;
+
+ DBG((dbg, LEVEL_1, "rematching col %u == row %u\n", c, r));
+ if (j == (unsigned)-1)
break;
- k = parent_row[j];
- l = j;
+ r = parent_row[j];
+ c = j;
}
/* End update the matching 20 */
/* Begin get ready for another stage 17 */
t = 0;
- for (l = 0; l < n; ++l) {
- parent_row[l] = -1;
- slack[l] = INF;
+ for (c = 0; c < num_cols; ++c) {
+ parent_row[c] = (unsigned) -1;
+ slack[c] = INT_MAX;
}
- for (k = 0; k < m; ++k) {
- if (col_mate[k] < 0) {
- DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
- unchosen_row[t++] = k;
+ for (r = 0; r < num_rows; ++r) {
+ if (col_mate[r] == (unsigned)-1) {
+ DBG((dbg, LEVEL_1, "node %u: unmatched row %u\n", t, r));
+ unchosen_row[t++] = r;
}
}
/* End get ready for another stage 17 */
done:
/* Begin double check the solution 23 */
- for (k = 0; k < m; ++k) {
- for (l = 0; l < n; ++l) {
- if (p->cost[k][l] < row_dec[k] - col_inc[l])
+ for (r = 0; r < num_rows; ++r) {
+ for (c = 0; c < num_cols; ++c) {
+ if ((int) cost[r*num_cols + c] < row_dec[r] - col_inc[c])
return -1;
}
}
- for (k = 0; k < m; ++k) {
- l = col_mate[k];
- if (l < 0 || p->cost[k][l] != row_dec[k] - col_inc[l])
+ for (r = 0; r < num_rows; ++r) {
+ c = col_mate[r];
+ if (c == (unsigned)-1
+ || cost[r*num_cols + c] != (unsigned) (row_dec[r] - col_inc[c]))
return -2;
}
- for (k = l = 0; l < n; ++l) {
- if (col_inc[l])
- k++;
+ for (r = c = 0; c < num_cols; ++c) {
+ if (col_inc[c])
+ r++;
}
- if (k > m)
+ if (r > num_rows)
return -3;
/* End double check the solution 23 */
/* End Hungarian algorithm 18 */
/* collect the assigned values */
- for (i = 0; i < m; ++i) {
- if (cost_threshold > 0 && p->cost[i][col_mate[i]] >= cost_threshold)
- assignment[i] = -1; /* remove matching having cost > threshold */
+ for (r = 0; r < num_rows; ++r) {
+ if (cost_threshold > 0
+ && cost[r*num_cols + col_mate[r]] >= cost_threshold)
+ assignment[r] = -1; /* remove matching having cost > threshold */
else
- assignment[i] = col_mate[i];
+ assignment[r] = col_mate[r];
}
/* In case of normal matching: remove impossible ones */
if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
- for (i = 0; i < m; ++i) {
- if (bitset_is_set(p->missing_left, i) || bitset_is_set(p->missing_right, col_mate[i]))
- assignment[i] = -1;
+ for (r = 0; r < num_rows; ++r) {
+ if (rbitset_is_set(p->missing_left, r)
+ || rbitset_is_set(p->missing_right, col_mate[r]))
+ assignment[r] = -1;
}
}
- for (k = 0; k < m; ++k) {
- for (l = 0; l < n; ++l) {
- p->cost[k][l] = p->cost[k][l] - row_dec[k] + col_inc[l];
+ for (r = 0; r < num_rows; ++r) {
+ for (c = 0; c < num_cols; ++c) {
+ cost[r*num_cols + c] = cost[r*num_cols + c] - row_dec[r] + col_inc[c];
}
}
- for (i = 0; i < m; ++i)
- cost += row_dec[i];
+ for (r = 0; r < num_rows; ++r)
+ res_cost += row_dec[r];
- for (i = 0; i < n; ++i)
- cost -= col_inc[i];
+ for (c = 0; c < num_cols; ++c)
+ res_cost -= col_inc[c];
- DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));
+ DBG((dbg, LEVEL_1, "Cost is %d\n", res_cost));
xfree(slack);
xfree(col_inc);
xfree(unchosen_row);
xfree(col_mate);
- *final_cost = cost;
+ if (final_cost != NULL)
+ *final_cost = res_cost;
return 0;
}
projects / libfirm / blobdiff