bipartite matching based on hungarian method
authorChristian Würdig <chriswue@ipd.info.uni-karlsruhe.de>
Wed, 13 Sep 2006 12:00:16 +0000 (12:00 +0000)
committerChristian Würdig <chriswue@ipd.info.uni-karlsruhe.de>
Wed, 13 Sep 2006 12:00:16 +0000 (12:00 +0000)
[r8238]

ir/adt/hungarian.c [new file with mode: 0644]
ir/adt/hungarian.h [new file with mode: 0644]

diff --git a/ir/adt/hungarian.c b/ir/adt/hungarian.c
new file mode 100644 (file)
index 0000000..d3e9507
--- /dev/null
@@ -0,0 +1,427 @@
+/********************************************************************
+ ********************************************************************
+ **
+ ** libhungarian by Cyrill Stachniss, 2004
+ **
+ ** Added and adapted to libFirm by Christian Wuerdig, 2006
+ **
+ ** Solving the Minimum Assignment Problem using the
+ ** Hungarian Method.
+ **
+ ** ** This file may be freely copied and distributed! **
+ **
+ ** Parts of the used code was originally provided by the
+ ** "Stanford GraphGase", but I made changes to this code.
+ ** As asked by  the copyright node of the "Stanford GraphGase",
+ ** I hereby proclaim that this file are *NOT* part of the
+ ** "Stanford GraphGase" distrubition!
+ **
+ ** This file is distributed in the hope that it will be useful,
+ ** but WITHOUT ANY WARRANTY; without even the implied
+ ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ ** PURPOSE.
+ **
+ ********************************************************************
+ ********************************************************************/
+
+/* $Id$ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#include "irtools.h"
+#include "xmalloc.h"
+#include "debug.h"
+#include "obst.h"
+
+#include "hungarian.h"
+
+#define INF (0x7FFFFFFF)
+
+struct _hungarian_problem_t {
+       int num_rows;
+       int num_cols;
+       int **cost;
+       int width;
+       int max_cost;
+       struct obstack obst;
+       DEBUG_ONLY(firm_dbg_module_t *dbg);
+};
+
+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;
+
+       fprintf(f , "\n");
+       for (i = 0; i < rows; i++) {
+               fprintf(f, " [");
+               for (j = 0; j < cols; j++) {
+                       fprintf(f, "%*d", width, C[i][j]);
+               }
+               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);
+}
+
+/**
+ * Create the object and allocate memory for the data structures.
+ */
+hungarian_problem_t *hungarian_new(int rows, int cols, int width) {
+       int i;
+       int max_cost = 0;
+       hungarian_problem_t *p = xmalloc(sizeof(*p));
+
+       memset(p, 0, sizeof(p));
+
+       FIRM_DBG_REGISTER(p->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);
+
+       p->num_rows = rows;
+       p->num_cols = cols;
+       p->width    = width;
+       p->cost     = (int **)get_init_mem(&p->obst, rows * sizeof(p->cost[0]));
+
+       /* allocate space for cost matrix */
+       for (i = 0; i < p->num_rows; i++)
+               p->cost[i] = (int *)get_init_mem(&p->obst, cols * sizeof(p->cost[0][0]));
+
+       return p;
+}
+
+/**
+ * Prepare the cost matrix.
+ */
+void hungarian_prepare_cost_matrix(hungarian_problem_t *p, int mode) {
+       int i, j;
+
+       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];
+                       }
+               }
+       }
+       else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
+               /* nothing to do */
+       }
+       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) {
+       assert(p->num_rows > left  && "Invalid row selected.");
+       assert(p->num_cols > right && "Invalid column selected.");
+
+       p->cost[left][right] = cost;
+       p->max_cost          = MAX(p->max_cost, cost);
+}
+
+/**
+ * Set cost[left][right] to 0.
+ */
+void hungarian_remv(hungarian_problem_t *p, int left, int right) {
+       assert(p->num_rows > left  && "Invalid row selected.");
+       assert(p->num_cols > right && "Invalid column selected.");
+
+       p->cost[left][right] = 0;
+}
+
+/**
+ * Frees all allocated memory.
+ */
+void hungarian_free(hungarian_problem_t* p) {
+       obstack_free(&p->obst, NULL);
+       xfree(p);
+}
+
+/**
+ * Do the assignment.
+ */
+int hungarian_solve(hungarian_problem_t* p, int *assignment) {
+       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     = xcalloc(p->num_rows, sizeof(col_mate[0]));
+       unchosen_row = xcalloc(p->num_rows, sizeof(unchosen_row[0]));
+       row_dec      = xcalloc(p->num_rows, sizeof(row_dec[0]));
+       slack_row    = xcalloc(p->num_rows, sizeof(slack_row[0]));
+
+       row_mate     = xcalloc(p->num_cols, sizeof(row_mate[0]));
+       parent_row   = xcalloc(p->num_cols, sizeof(parent_row[0]));
+       col_inc      = xcalloc(p->num_cols, sizeof(col_inc[0]));
+       slack        = xcalloc(p->num_cols, sizeof(slack[0]));
+
+#if 0
+       for (i = 0; i < p->num_rows; ++i) {
+               col_mate[i]     = 0;
+               unchosen_row[i] = 0;
+               row_dec[i]      = 0;
+               slack_row[i]=0;
+       }
+       for (j=0;j<p->num_cols;j++) {
+               row_mate[j]=0;
+               parent_row[j] = 0;
+               col_inc[j]=0;
+               slack[j]=0;
+       }
+#endif
+
+       memset(assignment, -1, m * sizeof(assignment[0]));
+
+       /* Begin subtract column minima in order to start with lots of zeros 12 */
+       DBG((p->dbg, LEVEL_1, "Using heuristic\n"));
+
+       for (l = 0; l < n; ++l) {
+               s = p->cost[0][l];
+
+               for (k = 1; k < m; ++k) {
+                       if (p->cost[k][l] < s)
+                               s = p->cost[k][l];
+               }
+
+               cost += s;
+
+               if (s != 0) {
+                       for (k = 0; k < m; ++k)
+                               p->cost[k][l] -= s;
+               }
+       }
+       /* 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;
+       }
+
+       for (k = 0; k < m; ++k) {
+               s = p->cost[k][0];
+
+               for (l = 1; l < n; ++l) {
+                       if (p->cost[k][l] < s)
+                               s = p->cost[k][l];
+               }
+
+               row_dec[k] = s;
+
+               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;
+                       }
+               }
+
+               col_mate[k] = -1;
+               DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
+               unchosen_row[t++] = k;
+row_done: ;
+       }
+       /* End initial state 16 */
+
+       /* Begin Hungarian algorithm 18 */
+       if (t == 0)
+               goto done;
+
+       unmatched = t;
+       while (1) {
+               DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", m - t));
+               q = 0;
+
+               while (1) {
+                       while (q < t) {
+                               /* Begin explore node q of the forest 19 */
+                               k = unchosen_row[q];
+                               s = row_dec[k];
+
+                               for (l = 0; l < n; ++l) {
+                                       if (slack[l]) {
+                                               int del = p->cost[k][l] - s + col_inc[l];
+
+                                               if (del < slack[l]) {
+                                                       if (del == 0) {
+                                                               if (row_mate[l] < 0)
+                                                                       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;
+                                                       }
+                                               }
+                                       }
+                               }
+                               /* End explore node q of the forest 19 */
+                               q++;
+                       }
+
+                       /* 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];
+                       }
+
+                       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) {
+                                               /* 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) {
+                                                               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];
+                                               }
+                                               /* End look at a new zero 22 */
+                                       }
+                               }
+                               else {
+                                       col_inc[l] += 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)
+                               break;
+
+                       k = parent_row[j];
+                       l = j;
+               }
+               /* End update the matching 20 */
+
+               if (--unmatched == 0)
+                       goto done;
+
+               /* Begin get ready for another stage 17 */
+               t = 0;
+               for (l = 0; l < n; ++l) {
+                       parent_row[l] = -1;
+                       slack[l]      = INF;
+               }
+
+               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;
+                       }
+               }
+               /* 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])
+                               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])
+                       return -2;
+       }
+
+       for (k = l = 0; l < n; ++l) {
+               if (col_inc[l])
+                       k++;
+       }
+
+       if (k > m)
+               return -3;
+       /* End double check the solution 23 */
+
+       /* End Hungarian algorithm 18 */
+
+       /* collect the assigned values */
+       for (i = 0; i < m; ++i) {
+               assignment[i] = col_mate[i];
+       }
+
+       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 (i = 0; i < m; ++i)
+               cost += row_dec[i];
+
+       for (i = 0; i < n; ++i)
+               cost -= col_inc[i];
+
+       DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));
+
+       xfree(slack);
+       xfree(col_inc);
+       xfree(parent_row);
+       xfree(row_mate);
+       xfree(slack_row);
+       xfree(row_dec);
+       xfree(unchosen_row);
+       xfree(col_mate);
+
+       return 0;
+}
diff --git a/ir/adt/hungarian.h b/ir/adt/hungarian.h
new file mode 100644 (file)
index 0000000..cc302e5
--- /dev/null
@@ -0,0 +1,82 @@
+/********************************************************************
+ ********************************************************************
+ **
+ ** libhungarian by Cyrill Stachniss, 2004
+ **
+ **
+ ** Solving the Minimum Assignment Problem using the
+ ** Hungarian Method.
+ **
+ ** ** This file may be freely copied and distributed! **
+ **
+ ** Parts of the used code was originally provided by the
+ ** "Stanford GraphGase", but I made changes to this code.
+ ** As asked by  the copyright node of the "Stanford GraphGase",
+ ** I hereby proclaim that this file are *NOT* part of the
+ ** "Stanford GraphGase" distrubition!
+ **
+ ** This file is distributed in the hope that it will be useful,
+ ** but WITHOUT ANY WARRANTY; without even the implied
+ ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ ** PURPOSE.
+ **
+ ********************************************************************
+ ********************************************************************/
+
+#ifndef _HUNGARIAN_H_
+#define _HUNGARIAN_H_
+
+#define HUNGARIAN_MODE_MINIMIZE_COST 0
+#define HUNGARIAN_MODE_MAXIMIZE_UTIL 1
+
+typedef struct _hungarian_problem_t hungarian_problem_t;
+
+/**
+ * This method initialize the hungarian_problem structure and init
+ * the cost matrix (missing lines or columns are filled with 0).
+ *
+ * @param rows  Number of rows in the given matrix
+ * @param cols  Number of cols in the given matrix
+ * @param width Element width for matrix dumping
+ * @return The problem object.
+ */
+hungarian_problem_t *hungarian_new(int rows, int cols, int width);
+
+/**
+ * Adds an edge from left to right with some costs.
+ */
+void hungarian_add(hungarian_problem_t *p, int left, int right, int cost);
+
+/**
+* Removes the edge from left to right.
+*/
+void hungarian_remv(hungarian_problem_t *p, int left, int right);
+
+/**
+ * Prepares the cost matrix, dependend on the given mode.
+ *
+ * @param p     The hungarian object
+ * @param mode  HUNGARIAN_MODE_MAXIMIZE_UTIL or HUNGARIAN_MODE_MINIMIZE_COST (default)
+ */
+void hungarian_prepare_cost_matrix(hungarian_problem_t *p, int mode);
+
+/**
+ * Destroys the hungarian object.
+ */
+void hungarian_free(hungarian_problem_t *p);
+
+/**
+ * This method computes the optimal assignment.
+ * @param p          The hungarian object
+ * @param assignment The final assignment
+ * @return Negative value if solution is invalid, 0 otherwise
+ */
+int hungarian_solve(hungarian_problem_t *p, int *assignment);
+
+/**
+ * Print the cost matrix.
+ * @param p The hungarian object
+ */
+void hungarian_print_costmatrix(hungarian_problem_t *p);
+
+#endif /* _HUNGARIAN_H_ */