|
| 1 | +/* |
| 2 | + * Based on: Jonker, R., & Volgenant, A. (1987). <i>A shortest augmenting path |
| 3 | + * algorithm for dense and sparse linear assignment problems</i>. Computing, |
| 4 | + * 38(4), 325-340. |
| 5 | + */ |
| 6 | +#include "cache.h" |
| 7 | +#include "linear-assignment.h" |
| 8 | + |
| 9 | +#define COST(column, row) cost[(column) + column_count * (row)] |
| 10 | + |
| 11 | +/* |
| 12 | + * The parameter `cost` is the cost matrix: the cost to assign column j to row |
| 13 | + * i is `cost[j + column_count * i]. |
| 14 | + */ |
| 15 | +void compute_assignment(int column_count, int row_count, int *cost, |
| 16 | + int *column2row, int *row2column) |
| 17 | +{ |
| 18 | + int *v, *d; |
| 19 | + int *free_row, free_count = 0, saved_free_count, *pred, *col; |
| 20 | + int i, j, phase; |
| 21 | + |
| 22 | + memset(column2row, -1, sizeof(int) * column_count); |
| 23 | + memset(row2column, -1, sizeof(int) * row_count); |
| 24 | + ALLOC_ARRAY(v, column_count); |
| 25 | + |
| 26 | + /* column reduction */ |
| 27 | + for (j = column_count - 1; j >= 0; j--) { |
| 28 | + int i1 = 0; |
| 29 | + |
| 30 | + for (i = 1; i < row_count; i++) |
| 31 | + if (COST(j, i1) > COST(j, i)) |
| 32 | + i1 = i; |
| 33 | + v[j] = COST(j, i1); |
| 34 | + if (row2column[i1] == -1) { |
| 35 | + /* row i1 unassigned */ |
| 36 | + row2column[i1] = j; |
| 37 | + column2row[j] = i1; |
| 38 | + } else { |
| 39 | + if (row2column[i1] >= 0) |
| 40 | + row2column[i1] = -2 - row2column[i1]; |
| 41 | + column2row[j] = -1; |
| 42 | + } |
| 43 | + } |
| 44 | + |
| 45 | + /* reduction transfer */ |
| 46 | + ALLOC_ARRAY(free_row, row_count); |
| 47 | + for (i = 0; i < row_count; i++) { |
| 48 | + int j1 = row2column[i]; |
| 49 | + if (j1 == -1) |
| 50 | + free_row[free_count++] = i; |
| 51 | + else if (j1 < -1) |
| 52 | + row2column[i] = -2 - j1; |
| 53 | + else { |
| 54 | + int min = COST(!j1, i) - v[!j1]; |
| 55 | + for (j = 1; j < column_count; j++) |
| 56 | + if (j != j1 && min > COST(j, i) - v[j]) |
| 57 | + min = COST(j, i) - v[j]; |
| 58 | + v[j1] -= min; |
| 59 | + } |
| 60 | + } |
| 61 | + |
| 62 | + if (free_count == |
| 63 | + (column_count < row_count ? row_count - column_count : 0)) { |
| 64 | + free(v); |
| 65 | + free(free_row); |
| 66 | + return; |
| 67 | + } |
| 68 | + |
| 69 | + /* augmenting row reduction */ |
| 70 | + for (phase = 0; phase < 2; phase++) { |
| 71 | + int k = 0; |
| 72 | + |
| 73 | + saved_free_count = free_count; |
| 74 | + free_count = 0; |
| 75 | + while (k < saved_free_count) { |
| 76 | + int u1, u2; |
| 77 | + int j1 = 0, j2, i0; |
| 78 | + |
| 79 | + i = free_row[k++]; |
| 80 | + u1 = COST(j1, i) - v[j1]; |
| 81 | + j2 = -1; |
| 82 | + u2 = INT_MAX; |
| 83 | + for (j = 1; j < column_count; j++) { |
| 84 | + int c = COST(j, i) - v[j]; |
| 85 | + if (u2 > c) { |
| 86 | + if (u1 < c) { |
| 87 | + u2 = c; |
| 88 | + j2 = j; |
| 89 | + } else { |
| 90 | + u2 = u1; |
| 91 | + u1 = c; |
| 92 | + j2 = j1; |
| 93 | + j1 = j; |
| 94 | + } |
| 95 | + } |
| 96 | + } |
| 97 | + if (j2 < 0) { |
| 98 | + j2 = j1; |
| 99 | + u2 = u1; |
| 100 | + } |
| 101 | + |
| 102 | + i0 = column2row[j1]; |
| 103 | + if (u1 < u2) |
| 104 | + v[j1] -= u2 - u1; |
| 105 | + else if (i0 >= 0) { |
| 106 | + j1 = j2; |
| 107 | + i0 = column2row[j1]; |
| 108 | + } |
| 109 | + |
| 110 | + if (i0 >= 0) { |
| 111 | + if (u1 < u2) |
| 112 | + free_row[--k] = i0; |
| 113 | + else |
| 114 | + free_row[free_count++] = i0; |
| 115 | + } |
| 116 | + row2column[i] = j1; |
| 117 | + column2row[j1] = i; |
| 118 | + } |
| 119 | + } |
| 120 | + |
| 121 | + /* augmentation */ |
| 122 | + saved_free_count = free_count; |
| 123 | + ALLOC_ARRAY(d, column_count); |
| 124 | + ALLOC_ARRAY(pred, column_count); |
| 125 | + ALLOC_ARRAY(col, column_count); |
| 126 | + for (free_count = 0; free_count < saved_free_count; free_count++) { |
| 127 | + int i1 = free_row[free_count], low = 0, up = 0, last, k; |
| 128 | + int min, c, u1; |
| 129 | + |
| 130 | + for (j = 0; j < column_count; j++) { |
| 131 | + d[j] = COST(j, i1) - v[j]; |
| 132 | + pred[j] = i1; |
| 133 | + col[j] = j; |
| 134 | + } |
| 135 | + |
| 136 | + j = -1; |
| 137 | + do { |
| 138 | + last = low; |
| 139 | + min = d[col[up++]]; |
| 140 | + for (k = up; k < column_count; k++) { |
| 141 | + j = col[k]; |
| 142 | + c = d[j]; |
| 143 | + if (c <= min) { |
| 144 | + if (c < min) { |
| 145 | + up = low; |
| 146 | + min = c; |
| 147 | + } |
| 148 | + col[k] = col[up]; |
| 149 | + col[up++] = j; |
| 150 | + } |
| 151 | + } |
| 152 | + for (k = low; k < up; k++) |
| 153 | + if (column2row[col[k]] == -1) |
| 154 | + goto update; |
| 155 | + |
| 156 | + /* scan a row */ |
| 157 | + do { |
| 158 | + int j1 = col[low++]; |
| 159 | + |
| 160 | + i = column2row[j1]; |
| 161 | + u1 = COST(j1, i) - v[j1] - min; |
| 162 | + for (k = up; k < column_count; k++) { |
| 163 | + j = col[k]; |
| 164 | + c = COST(j, i) - v[j] - u1; |
| 165 | + if (c < d[j]) { |
| 166 | + d[j] = c; |
| 167 | + pred[j] = i; |
| 168 | + if (c == min) { |
| 169 | + if (column2row[j] == -1) |
| 170 | + goto update; |
| 171 | + col[k] = col[up]; |
| 172 | + col[up++] = j; |
| 173 | + } |
| 174 | + } |
| 175 | + } |
| 176 | + } while (low != up); |
| 177 | + } while (low == up); |
| 178 | + |
| 179 | +update: |
| 180 | + /* updating of the column pieces */ |
| 181 | + for (k = 0; k < last; k++) { |
| 182 | + int j1 = col[k]; |
| 183 | + v[j1] += d[j1] - min; |
| 184 | + } |
| 185 | + |
| 186 | + /* augmentation */ |
| 187 | + do { |
| 188 | + if (j < 0) |
| 189 | + BUG("negative j: %d", j); |
| 190 | + i = pred[j]; |
| 191 | + column2row[j] = i; |
| 192 | + k = j; |
| 193 | + j = row2column[i]; |
| 194 | + row2column[i] = k; |
| 195 | + } while (i1 != i); |
| 196 | + } |
| 197 | + |
| 198 | + free(col); |
| 199 | + free(pred); |
| 200 | + free(d); |
| 201 | + free(v); |
| 202 | + free(free_row); |
| 203 | +} |
0 commit comments