| 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 | } |