git

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 "git-compat-util.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
	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