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