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