nsz Git - libfirm/blob - ir/adt/hungarian.c

   1 /********************************************************************
   2  ********************************************************************
   3  **
   4  ** libhungarian by Cyrill Stachniss, 2004
   5  **
   6  ** Added and adapted to libFirm by Christian Wuerdig, 2006
   7  **
   8  ** Solving the Minimum Assignment Problem using the
   9  ** Hungarian Method.
  10  **
  11  ** ** This file may be freely copied and distributed! **
  12  **
  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!
  18  **
  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
  22  ** PURPOSE.
  23  **
  24  ********************************************************************
  25  ********************************************************************/
  26
  27 /* $Id$ */
  28
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <assert.h>
  32
  33 #include "irtools.h"
  34 #include "xmalloc.h"
  35 #include "debug.h"
  36 #include "obst.h"
  37
  38 #include "hungarian.h"
  39
  40 #define INF (0x7FFFFFFF)
  41
  42 struct _hungarian_problem_t {
  43         int num_rows;
  44         int num_cols;
  45         int **cost;
  46         int width;
  47         int max_cost;
  48         struct obstack obst;
  49         DEBUG_ONLY(firm_dbg_module_t *dbg);
  50 };
  51
  52 static INLINE void *get_init_mem(struct obstack *obst, long sz) {
  53         void *p = obstack_alloc(obst, sz);
  54         memset(p, 0, sz);
  55         return p;
  56 }
  57
  58 static void hungarian_dump_f(FILE *f, int **C, int rows, int cols, int width) {
  59         int i, j;
  60
  61         fprintf(f , "\n");
  62         for (i = 0; i < rows; i++) {
  63                 fprintf(f, " [");
  64                 for (j = 0; j < cols; j++) {
  65                         fprintf(f, "%*d", width, C[i][j]);
  66                 }
  67                 fprintf(f, "]\n");
  68         }
  69         fprintf(f, "\n");
  70 }
  71
  72 void hungarian_print_costmatrix(hungarian_problem_t *p) {
  73         hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, p->width);
  74 }
  75
  76 /**
  77  * Create the object and allocate memory for the data structures.
  78  */
  79 hungarian_problem_t *hungarian_new(int rows, int cols, int width) {
  80         int i;
  81         int max_cost = 0;
  82         hungarian_problem_t *p = xmalloc(sizeof(*p));
  83
  84         memset(p, 0, sizeof(p));
  85
  86         FIRM_DBG_REGISTER(p->dbg, "firm.hungarian");
  87
  88         /*
  89                 Is the number of cols  not equal to number of rows ?
  90                 If yes, expand with 0 - cols / 0 - cols
  91         */
  92         rows = MAX(cols, rows);
  93         cols = rows;
  94
  95         obstack_init(&p->obst);
  96
  97         p->num_rows = rows;
  98         p->num_cols = cols;
  99         p->width    = width;
 100         p->cost     = (int **)get_init_mem(&p->obst, rows * sizeof(p->cost[0]));
 101
 102         /* allocate space for cost matrix */
 103         for (i = 0; i < p->num_rows; i++)
 104                 p->cost[i] = (int *)get_init_mem(&p->obst, cols * sizeof(p->cost[0][0]));
 105
 106         return p;
 107 }
 108
 109 /**
 110  * Prepare the cost matrix.
 111  */
 112 void hungarian_prepare_cost_matrix(hungarian_problem_t *p, int mode) {
 113         int i, j;
 114
 115         if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
 116                 for (i = 0; i < p->num_rows; i++) {
 117                         for (j = 0; j < p->num_cols; j++) {
 118                                 p->cost[i][j] = p->max_cost - p->cost[i][j];
 119                         }
 120                 }
 121         }
 122         else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
 123                 /* nothing to do */
 124         }
 125         else
 126                 fprintf(stderr, "Unknown mode. Mode was set to HUNGARIAN_MODE_MINIMIZE_COST.\n");
 127 }
 128
 129 /**
 130  * Set cost[left][right] to cost.
 131  */
 132 void hungarian_add(hungarian_problem_t *p, int left, int right, int cost) {
 133         assert(p->num_rows > left  && "Invalid row selected.");
 134         assert(p->num_cols > right && "Invalid column selected.");
 135
 136         p->cost[left][right] = cost;
 137         p->max_cost          = MAX(p->max_cost, cost);
 138 }
 139
 140 /**
 141  * Set cost[left][right] to 0.
 142  */
 143 void hungarian_remv(hungarian_problem_t *p, int left, int right) {
 144         assert(p->num_rows > left  && "Invalid row selected.");
 145         assert(p->num_cols > right && "Invalid column selected.");
 146
 147         p->cost[left][right] = 0;
 148 }
 149
 150 /**
 151  * Frees all allocated memory.
 152  */
 153 void hungarian_free(hungarian_problem_t* p) {
 154         obstack_free(&p->obst, NULL);
 155         xfree(p);
 156 }
 157
 158 /**
 159  * Do the assignment.
 160  */
 161 int hungarian_solve(hungarian_problem_t* p, int *assignment) {
 162         int i, j, m, n, k, l, s, t, q, unmatched, cost;
 163         int *col_mate;
 164         int *row_mate;
 165         int *parent_row;
 166         int *unchosen_row;
 167         int *row_dec;
 168         int *col_inc;
 169         int *slack;
 170         int *slack_row;
 171
 172         cost = 0;
 173         m    = p->num_rows;
 174         n    = p->num_cols;
 175
 176         col_mate     = xcalloc(p->num_rows, sizeof(col_mate[0]));
 177         unchosen_row = xcalloc(p->num_rows, sizeof(unchosen_row[0]));
 178         row_dec      = xcalloc(p->num_rows, sizeof(row_dec[0]));
 179         slack_row    = xcalloc(p->num_rows, sizeof(slack_row[0]));
 180
 181         row_mate     = xcalloc(p->num_cols, sizeof(row_mate[0]));
 182         parent_row   = xcalloc(p->num_cols, sizeof(parent_row[0]));
 183         col_inc      = xcalloc(p->num_cols, sizeof(col_inc[0]));
 184         slack        = xcalloc(p->num_cols, sizeof(slack[0]));
 185
 186 #if 0
 187         for (i = 0; i < p->num_rows; ++i) {
 188                 col_mate[i]     = 0;
 189                 unchosen_row[i] = 0;
 190                 row_dec[i]      = 0;
 191                 slack_row[i]=0;
 192         }
 193         for (j=0;j<p->num_cols;j++) {
 194                 row_mate[j]=0;
 195                 parent_row[j] = 0;
 196                 col_inc[j]=0;
 197                 slack[j]=0;
 198         }
 199 #endif
 200
 201         memset(assignment, -1, m * sizeof(assignment[0]));
 202
 203         /* Begin subtract column minima in order to start with lots of zeros 12 */
 204         DBG((p->dbg, LEVEL_1, "Using heuristic\n"));
 205
 206         for (l = 0; l < n; ++l) {
 207                 s = p->cost[0][l];
 208
 209                 for (k = 1; k < m; ++k) {
 210                         if (p->cost[k][l] < s)
 211                                 s = p->cost[k][l];
 212                 }
 213
 214                 cost += s;
 215
 216                 if (s != 0) {
 217                         for (k = 0; k < m; ++k)
 218                                 p->cost[k][l] -= s;
 219                 }
 220         }
 221         /* End subtract column minima in order to start with lots of zeros 12 */
 222
 223         /* Begin initial state 16 */
 224         t = 0;
 225         for (l = 0; l < n; ++l) {
 226                 row_mate[l]   = -1;
 227                 parent_row[l] = -1;
 228                 col_inc[l]    = 0;
 229                 slack[l]      = INF;
 230         }
 231
 232         for (k = 0; k < m; ++k) {
 233                 s = p->cost[k][0];
 234
 235                 for (l = 1; l < n; ++l) {
 236                         if (p->cost[k][l] < s)
 237                                 s = p->cost[k][l];
 238                 }
 239
 240                 row_dec[k] = s;
 241
 242                 for (l = 0; l < n; ++l) {
 243                         if (s == p->cost[k][l] && row_mate[l] < 0) {
 244                                 col_mate[k] = l;
 245                                 row_mate[l] = k;
 246                                 DBG((p->dbg, LEVEL_1, "matching col %d == row %d\n", l, k));
 247                                 goto row_done;
 248                         }
 249                 }
 250
 251                 col_mate[k] = -1;
 252                 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
 253                 unchosen_row[t++] = k;
 254 row_done: ;
 255         }
 256         /* End initial state 16 */
 257
 258         /* Begin Hungarian algorithm 18 */
 259         if (t == 0)
 260                 goto done;
 261
 262         unmatched = t;
 263         while (1) {
 264                 DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", m - t));
 265                 q = 0;
 266
 267                 while (1) {
 268                         while (q < t) {
 269                                 /* Begin explore node q of the forest 19 */
 270                                 k = unchosen_row[q];
 271                                 s = row_dec[k];
 272
 273                                 for (l = 0; l < n; ++l) {
 274                                         if (slack[l]) {
 275                                                 int del = p->cost[k][l] - s + col_inc[l];
 276
 277                                                 if (del < slack[l]) {
 278                                                         if (del == 0) {
 279                                                                 if (row_mate[l] < 0)
 280                                                                         goto breakthru;
 281
 282                                                                 slack[l]      = 0;
 283                                                                 parent_row[l] = k;
 284                                                                 DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
 285                                                                 unchosen_row[t++] = row_mate[l];
 286                                                         }
 287                                                         else {
 288                                                                 slack[l]     = del;
 289                                                                 slack_row[l] = k;
 290                                                         }
 291                                                 }
 292                                         }
 293                                 }
 294                                 /* End explore node q of the forest 19 */
 295                                 q++;
 296                         }
 297
 298                         /* Begin introduce a new zero into the matrix 21 */
 299                         s = INF;
 300                         for (l = 0; l < n; ++l) {
 301                                 if (slack[l] && slack[l] < s)
 302                                         s = slack[l];
 303                         }
 304
 305                         for (q = 0; q < t; ++q)
 306                                 row_dec[unchosen_row[q]] += s;
 307
 308                         for (l = 0; l < n; ++l) {
 309                                 if (slack[l]) {
 310                                         slack[l] -= s;
 311                                         if (slack[l] == 0) {
 312                                                 /* Begin look at a new zero 22 */
 313                                                 k = slack_row[l];
 314                                                 DBG((p->dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%d, %d]\n", s, k, l));
 315                                                 if (row_mate[l] < 0) {
 316                                                         for (j = l + 1; j < n; ++j) {
 317                                                                 if (slack[j] == 0)
 318                                                                         col_inc[j] += s;
 319                                                         }
 320                                                         goto breakthru;
 321                                                 }
 322                                                 else {
 323                                                         parent_row[l] = k;
 324                                                         DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[l], l, k));
 325                                                         unchosen_row[t++] = row_mate[l];
 326                                                 }
 327                                                 /* End look at a new zero 22 */
 328                                         }
 329                                 }
 330                                 else {
 331                                         col_inc[l] += s;
 332                                 }
 333                         }
 334                         /* End introduce a new zero into the matrix 21 */
 335                 }
 336 breakthru:
 337                 /* Begin update the matching 20 */
 338                 DBG((p->dbg, LEVEL_1, "Breakthrough at node %d of %d.\n", q, t));
 339                 while (1) {
 340                         j           = col_mate[k];
 341                         col_mate[k] = l;
 342                         row_mate[l] = k;
 343
 344                         DBG((p->dbg, LEVEL_1, "rematching col %d == row %d\n", l, k));
 345                         if (j < 0)
 346                                 break;
 347
 348                         k = parent_row[j];
 349                         l = j;
 350                 }
 351                 /* End update the matching 20 */
 352
 353                 if (--unmatched == 0)
 354                         goto done;
 355
 356                 /* Begin get ready for another stage 17 */
 357                 t = 0;
 358                 for (l = 0; l < n; ++l) {
 359                         parent_row[l] = -1;
 360                         slack[l]      = INF;
 361                 }
 362
 363                 for (k = 0; k < m; ++k) {
 364                         if (col_mate[k] < 0) {
 365                                 DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, k));
 366                                 unchosen_row[t++] = k;
 367                         }
 368                 }
 369                 /* End get ready for another stage 17 */
 370         }
 371 done:
 372
 373         /* Begin double check the solution 23 */
 374         for (k = 0; k < m; ++k) {
 375                 for (l = 0; l < n; ++l) {
 376                         if (p->cost[k][l] < row_dec[k] - col_inc[l])
 377                                 return -1;
 378                 }
 379         }
 380
 381         for (k = 0; k < m; ++k) {
 382                 l = col_mate[k];
 383                 if (l < 0 || p->cost[k][l] != row_dec[k] - col_inc[l])
 384                         return -2;
 385         }
 386
 387         for (k = l = 0; l < n; ++l) {
 388                 if (col_inc[l])
 389                         k++;
 390         }
 391
 392         if (k > m)
 393                 return -3;
 394         /* End double check the solution 23 */
 395
 396         /* End Hungarian algorithm 18 */
 397
 398         /* collect the assigned values */
 399         for (i = 0; i < m; ++i) {
 400                 assignment[i] = col_mate[i];
 401         }
 402
 403         for (k = 0; k < m; ++k) {
 404                 for (l = 0; l < n; ++l) {
 405                         p->cost[k][l] = p->cost[k][l] - row_dec[k] + col_inc[l];
 406                 }
 407         }
 408
 409         for (i = 0; i < m; ++i)
 410                 cost += row_dec[i];
 411
 412         for (i = 0; i < n; ++i)
 413                 cost -= col_inc[i];
 414
 415         DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));
 416
 417         xfree(slack);
 418         xfree(col_inc);
 419         xfree(parent_row);
 420         xfree(row_mate);
 421         xfree(slack_row);
 422         xfree(row_dec);
 423         xfree(unchosen_row);
 424         xfree(col_mate);
 425
 426         return 0;
 427 }