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