blob: ecffc09be6ec5c1a9646ac152a899f3415f6321d [file] [log] [blame]
 /* * Based on: Jonker, R., & Volgenant, A. (1987). A shortest augmenting path * algorithm for dense and sparse linear assignment problems. Computing, * 38(4), 325-340. */ #include "cache.h" #include "linear-assignment.h" #define COST(column, row) cost[(column) + column_count * (row)] /* * The parameter `cost` is the cost matrix: the cost to assign column j to row * i is `cost[j + column_count * i]. */ void compute_assignment(int column_count, int row_count, int *cost, int *column2row, int *row2column) { int *v, *d; int *free_row, free_count = 0, saved_free_count, *pred, *col; int i, j, phase; if (column_count < 2) { memset(column2row, 0, sizeof(int) * column_count); memset(row2column, 0, sizeof(int) * row_count); return; } memset(column2row, -1, sizeof(int) * column_count); memset(row2column, -1, sizeof(int) * row_count); ALLOC_ARRAY(v, column_count); /* column reduction */ for (j = column_count - 1; j >= 0; j--) { int i1 = 0; for (i = 1; i < row_count; i++) if (COST(j, i1) > COST(j, i)) i1 = i; v[j] = COST(j, i1); if (row2column[i1] == -1) { /* row i1 unassigned */ row2column[i1] = j; column2row[j] = i1; } else { if (row2column[i1] >= 0) row2column[i1] = -2 - row2column[i1]; column2row[j] = -1; } } /* reduction transfer */ ALLOC_ARRAY(free_row, row_count); for (i = 0; i < row_count; i++) { int j1 = row2column[i]; if (j1 == -1) free_row[free_count++] = i; else if (j1 < -1) row2column[i] = -2 - j1; else { int min = COST(!j1, i) - v[!j1]; for (j = 1; j < column_count; j++) if (j != j1 && min > COST(j, i) - v[j]) min = COST(j, i) - v[j]; v[j1] -= min; } } if (free_count == (column_count < row_count ? row_count - column_count : 0)) { free(v); free(free_row); return; } /* augmenting row reduction */ for (phase = 0; phase < 2; phase++) { int k = 0; saved_free_count = free_count; free_count = 0; while (k < saved_free_count) { int u1, u2; int j1 = 0, j2, i0; i = free_row[k++]; u1 = COST(j1, i) - v[j1]; j2 = -1; u2 = INT_MAX; for (j = 1; j < column_count; j++) { int c = COST(j, i) - v[j]; if (u2 > c) { if (u1 < c) { u2 = c; j2 = j; } else { u2 = u1; u1 = c; j2 = j1; j1 = j; } } } if (j2 < 0) { j2 = j1; u2 = u1; } i0 = column2row[j1]; if (u1 < u2) v[j1] -= u2 - u1; else if (i0 >= 0) { j1 = j2; i0 = column2row[j1]; } if (i0 >= 0) { if (u1 < u2) free_row[--k] = i0; else free_row[free_count++] = i0; } row2column[i] = j1; column2row[j1] = i; } } /* augmentation */ saved_free_count = free_count; ALLOC_ARRAY(d, column_count); ALLOC_ARRAY(pred, column_count); ALLOC_ARRAY(col, column_count); for (free_count = 0; free_count < saved_free_count; free_count++) { int i1 = free_row[free_count], low = 0, up = 0, last, k; int min, c, u1; for (j = 0; j < column_count; j++) { d[j] = COST(j, i1) - v[j]; pred[j] = i1; col[j] = j; } j = -1; do { last = low; min = d[col[up++]]; for (k = up; k < column_count; k++) { j = col[k]; c = d[j]; if (c <= min) { if (c < min) { up = low; min = c; } col[k] = col[up]; col[up++] = j; } } for (k = low; k < up; k++) if (column2row[col[k]] == -1) goto update; /* scan a row */ do { int j1 = col[low++]; i = column2row[j1]; u1 = COST(j1, i) - v[j1] - min; for (k = up; k < column_count; k++) { j = col[k]; c = COST(j, i) - v[j] - u1; if (c < d[j]) { d[j] = c; pred[j] = i; if (c == min) { if (column2row[j] == -1) goto update; col[k] = col[up]; col[up++] = j; } } } } while (low != up); } while (low == up); update: /* updating of the column pieces */ for (k = 0; k < last; k++) { int j1 = col[k]; v[j1] += d[j1] - min; } /* augmentation */ do { if (j < 0) BUG("negative j: %d", j); i = pred[j]; column2row[j] = i; SWAP(j, row2column[i]); } while (i1 != i); } free(col); free(pred); free(d); free(v); free(free_row); }