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