code/postprocess/hungarian.c

   1 /********************************************************************
   2  ********************************************************************
   3  **
   4  ** libhungarian by Cyrill Stachniss, 2004
   5  ** -- modified (very) slightly by Benjamin Auder, 2010
   6  ** -- (verbose and printing routines deleted, because used in a R package)
   7  **
   8  **
   9  ** Solving the Minimum Assignment Problem using the
  10  ** Hungarian Method.
  11  **
  12  ** ** This file may be freely copied and distributed! **
  13  **
  14  ** Parts of the used code was originally provided by the
  15  ** "Stanford GraphGase", but I made changes to this code.
  16  ** As asked by  the copyright node of the "Stanford GraphGase",
  17  ** I hereby proclaim that this file are *NOT* part of the
  18  ** "Stanford GraphGase" distribution!
  19  **
  20  ** This file is distributed in the hope that it will be useful,
  21  ** but WITHOUT ANY WARRANTY; without even the implied
  22  ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  23  ** PURPOSE.
  24  **
  25  ********************************************************************
  26  ********************************************************************/
  27
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30
  31 #define HUNGARIAN_NOT_ASSIGNED 0
  32 #define HUNGARIAN_ASSIGNED 1
  33
  34 #define HUNGARIAN_MODE_MINIMIZE_COST   0
  35 #define HUNGARIAN_MODE_MAXIMIZE_UTIL 1
  36
  37 typedef struct {
  38     int num_rows;
  39     int num_cols;
  40     int** cost;
  41     int** assignment;
  42 } hungarian_problem_t;
  43
  44 #define INF (0x7FFFFFFF)
  45 #define hungarian_test_alloc(X) do {if ((void *)(X) == NULL) fprintf(stderr, "Out of memory in %s, (%s, line %d).\n", __FUNCTION__, __FILE__, __LINE__); } while (0)
  46
  47 /** This method initialize the hungarian_problem structure and init
  48  *  the  cost matrices (missing lines or columns are filled with 0).
  49  *  It returns the size of the quadratic(!) assignment matrix. **/
  50 void hungarian_init(hungarian_problem_t* p, int* cost_matrix, int rows, int cols, int mode) {
  51
  52     int i,j, org_cols, org_rows;
  53     int max_cost;
  54     max_cost = 0;
  55
  56     org_cols = cols;
  57     org_rows = rows;
  58
  59     // is the number of cols  not equal to number of rows ?
  60     // if yes, expand with 0-cols / 0-cols
  61     rows = (cols < rows) ? rows : cols;
  62     cols = rows;
  63
  64     p->num_rows = rows;
  65     p->num_cols = cols;
  66
  67     p->cost = (int**)calloc(rows,sizeof(int*));
  68     hungarian_test_alloc(p->cost);
  69     p->assignment = (int**)calloc(rows,sizeof(int*));
  70     hungarian_test_alloc(p->assignment);
  71
  72     for(i=0; i<p->num_rows; i++) {
  73         p->cost[i] = (int*)calloc(cols,sizeof(int));
  74         hungarian_test_alloc(p->cost[i]);
  75         p->assignment[i] = (int*)calloc(cols,sizeof(int));
  76         hungarian_test_alloc(p->assignment[i]);
  77         for(j=0; j<p->num_cols; j++) {
  78         p->cost[i][j] =  (i < org_rows && j < org_cols) ? cost_matrix[i*(p->num_cols)+j] : 0;
  79         p->assignment[i][j] = 0;
  80
  81         if (max_cost < p->cost[i][j])
  82         max_cost = p->cost[i][j];
  83         }
  84     }
  85
  86     if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
  87         for(i=0; i<p->num_rows; i++) {
  88             for(j=0; j<p->num_cols; j++) {
  89                 p->cost[i][j] =  max_cost - p->cost[i][j];
  90             }
  91         }
  92     }
  93 }
  94
  95 /** Free the memory allocated by init. **/
  96 void hungarian_free(hungarian_problem_t* p) {
  97     int i;
  98     for(i=0; i<p->num_rows; i++) {
  99         free(p->cost[i]);
 100         free(p->assignment[i]);
 101     }
 102     free(p->cost);
 103     free(p->assignment);
 104     p->cost = NULL;
 105     p->assignment = NULL;
 106 }
 107
 108 /** This method computes the optimal assignment. **/
 109 void hungarian_solve(hungarian_problem_t* p) {
 110
 111     int i, j, m, n, k, l, s, t, q, unmatched, cost;
 112     int* col_mate;
 113     int* row_mate;
 114     int* parent_row;
 115     int* unchosen_row;
 116     int* row_dec;
 117     int* col_inc;
 118     int* slack;
 119     int* slack_row;
 120
 121     cost=0;
 122     m =p->num_rows;
 123     n =p->num_cols;
 124
 125     col_mate = (int*)calloc(p->num_rows,sizeof(int));
 126     hungarian_test_alloc(col_mate);
 127     unchosen_row = (int*)calloc(p->num_rows,sizeof(int));
 128     hungarian_test_alloc(unchosen_row);
 129     row_dec  = (int*)calloc(p->num_rows,sizeof(int));
 130     hungarian_test_alloc(row_dec);
 131     slack_row  = (int*)calloc(p->num_rows,sizeof(int));
 132     hungarian_test_alloc(slack_row);
 133
 134     row_mate = (int*)calloc(p->num_cols,sizeof(int));
 135     hungarian_test_alloc(row_mate);
 136     parent_row = (int*)calloc(p->num_cols,sizeof(int));
 137     hungarian_test_alloc(parent_row);
 138     col_inc = (int*)calloc(p->num_cols,sizeof(int));
 139     hungarian_test_alloc(col_inc);
 140     slack = (int*)calloc(p->num_cols,sizeof(int));
 141     hungarian_test_alloc(slack);
 142
 143     for (i=0;i<p->num_rows;i++) {
 144         col_mate[i]=0;
 145         unchosen_row[i]=0;
 146         row_dec[i]=0;
 147         slack_row[i]=0;
 148     }
 149     for (j=0;j<p->num_cols;j++) {
 150         row_mate[j]=0;
 151         parent_row[j] = 0;
 152         col_inc[j]=0;
 153         slack[j]=0;
 154     }
 155
 156     for (i=0;i<p->num_rows;++i)
 157         for (j=0;j<p->num_cols;++j)
 158             p->assignment[i][j]=HUNGARIAN_NOT_ASSIGNED;
 159
 160     // Begin subtract column minima in order to start with lots of zeroes 12
 161     for (l=0;l<n;l++) {
 162         s=p->cost[0][l];
 163         for (k=1;k<m;k++)
 164         if (p->cost[k][l]<s)
 165         s=p->cost[k][l];
 166         cost+=s;
 167         if (s!=0)
 168         for (k=0;k<m;k++)
 169         p->cost[k][l]-=s;
 170     }
 171     // End subtract column minima in order to start with lots of zeroes 12
 172
 173     // Begin initial state 16
 174     t=0;
 175     for (l=0;l<n;l++) {
 176         row_mate[l]= -1;
 177         parent_row[l]= -1;
 178         col_inc[l]=0;
 179         slack[l]=INF;
 180     }
 181     for (k=0;k<m;k++) {
 182         s=p->cost[k][0];
 183         for (l=1;l<n;l++)
 184         if (p->cost[k][l]<s)
 185         s=p->cost[k][l];
 186         row_dec[k]=s;
 187         for (l=0;l<n;l++)
 188         if (s==p->cost[k][l] && row_mate[l]<0) {
 189             col_mate[k]=l;
 190             row_mate[l]=k;
 191             goto row_done;
 192         }
 193         col_mate[k]= -1;
 194         unchosen_row[t++]=k;
 195         row_done:
 196         ;
 197     }
 198     // End initial state 16
 199
 200     // Begin Hungarian algorithm 18
 201     if (t==0)
 202         goto done;
 203     unmatched=t;
 204     while (1) {
 205         q=0;
 206         while (1) {
 207             while (q<t) {
 208                 // Begin explore node q of the forest 19
 209                 {
 210                 k=unchosen_row[q];
 211                 s=row_dec[k];
 212                 for (l=0;l<n;l++)
 213                     if (slack[l]) {
 214                         int del;
 215                         del=p->cost[k][l]-s+col_inc[l];
 216                         if (del<slack[l]) {
 217                             if (del==0) {
 218                                 if (row_mate[l]<0)
 219                                     goto breakthru;
 220                                 slack[l]=0;
 221                                 parent_row[l]=k;
 222                                 unchosen_row[t++]=row_mate[l];
 223                             }
 224                             else {
 225                                 slack[l]=del;
 226                                 slack_row[l]=k;
 227                             }
 228                         }
 229                     }
 230                 }
 231                 // End explore node q of the forest 19
 232                 q++;
 233             }
 234
 235             // Begin introduce a new zero into the matrix 21
 236             s=INF;
 237             for (l=0;l<n;l++)
 238                 if (slack[l] && slack[l]<s)
 239                     s=slack[l];
 240             for (q=0;q<t;q++)
 241                 row_dec[unchosen_row[q]]+=s;
 242             for (l=0;l<n;l++) {
 243                 if (slack[l]) {
 244                     slack[l]-=s;
 245                     if (slack[l]==0) {
 246                         // Begin look at a new zero 22
 247                         k=slack_row[l];
 248                         if (row_mate[l]<0) {
 249                             for (j=l+1;j<n;j++)
 250                                 if (slack[j]==0)
 251                                     col_inc[j]+=s;
 252                             goto breakthru;
 253                         }
 254                         else {
 255                             parent_row[l]=k;
 256                             unchosen_row[t++]=row_mate[l];
 257                         }
 258                         // End look at a new zero 22
 259                     }
 260                 }
 261                 else
 262                     col_inc[l]+=s;
 263                     // End introduce a new zero into the matrix 21
 264             }
 265         }
 266         breakthru:
 267         // Begin update the matching 20
 268         while (1) {
 269             j=col_mate[k];
 270             col_mate[k]=l;
 271             row_mate[l]=k;
 272             if (j<0)
 273             break;
 274             k=parent_row[j];
 275             l=j;
 276         }
 277         // End update the matching 20
 278         if (--unmatched==0)
 279             goto done;
 280         // Begin get ready for another stage 17
 281         t=0;
 282         for (l=0;l<n;l++) {
 283             parent_row[l]= -1;
 284             slack[l]=INF;
 285         }
 286         for (k=0;k<m;k++)
 287             if (col_mate[k]<0)
 288                 unchosen_row[t++]=k;
 289         // End get ready for another stage 17
 290     }
 291
 292     done:
 293
 294     // Begin doublecheck the solution 23
 295     for (k=0;k<m;k++)
 296         for (l=0;l<n;l++)
 297         if (p->cost[k][l]<row_dec[k]-col_inc[l])
 298         exit(0);
 299     for (k=0;k<m;k++) {
 300         l=col_mate[k];
 301         if (l<0 || p->cost[k][l]!=row_dec[k]-col_inc[l])
 302         exit(0);
 303     }
 304     k=0;
 305     for (l=0;l<n;l++)
 306         if (col_inc[l])
 307             k++;
 308     if (k>m)
 309         exit(0);
 310     // End doublecheck the solution 23
 311     // End Hungarian algorithm 18
 312
 313     for (i=0;i<m;++i) {
 314         p->assignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED;
 315         /*TRACE("%d - %d\n", i, col_mate[i]);*/
 316     }
 317     for (k=0;k<m;++k) {
 318         for (l=0;l<n;++l) {
 319             p->cost[k][l]=p->cost[k][l]-row_dec[k]+col_inc[l];
 320         }
 321     }
 322     for (i=0;i<m;i++)
 323         cost+=row_dec[i];
 324     for (i=0;i<n;i++)
 325         cost-=col_inc[i];
 326
 327     free(slack);
 328     free(col_inc);
 329     free(parent_row);
 330     free(row_mate);
 331     free(slack_row);
 332     free(row_dec);
 333     free(unchosen_row);
 334     free(col_mate);
 335 }
 336
 337 /* get the optimal assignment, by calling hungarian_solve above */
 338 void computeCoefSimil(int* P1, int* P2, int* maxInd, int* n, double* coefSimil) {
 339
 340     /* first, determine weights by computing intersections
 341      * reads like: weights[i*(*n)+j] == gain when cluster i in P1
 342      * is assigned to cluster j in P2
 343      */
 344     int i,j,k,inter;
 345     int* utils = malloc((*maxInd) * (*maxInd) * sizeof(int));
 346     for (i=1; i <= *maxInd; i++) {
 347         for (j=1; j <= *maxInd; j++) {
 348             //try label i with label j
 349             inter=0;
 350             for (k=0; k < *n; k++) {
 351                 if (P1[k]==i && P2[k]==j) inter++;
 352             }
 353             utils[(i-1)*(*maxInd)+(j-1)] = inter;
 354         }
 355     }
 356
 357     //then solve problem :
 358     hungarian_problem_t p;
 359     hungarian_init(&p, utils, *maxInd, *maxInd, HUNGARIAN_MODE_MAXIMIZE_UTIL);
 360     hungarian_solve(&p);
 361
 362     //and now get optimal assignment :
 363     *coefSimil = 0.0;
 364     for (i=0; i < *maxInd; i++) {
 365         for (j=0; j < *maxInd; j++) {
 366             if (p.assignment[i][j]) *coefSimil += utils[i*(*maxInd)+j];
 367         }
 368     }
 369     *coefSimil /= (double)(*n);
 370     hungarian_free(&p);
 371 }