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