linear-assignment: a function to solve least-cost assignment problems
[git/git.git] / linear-assignment.c
CommitLineData
22d87333
JS
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 */
15void 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
179update:
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 SWAP(j, row2column[i]);
193 } while (i1 != i);
194 }
195
196 free(col);
197 free(pred);
198 free(d);
199 free(v);
200 free(free_row);
201}