Merge branch 'ep/worktree-quiet-option'
[git] / linear-assignment.c
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                         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 }