[thirdparty/git.git] / linear-assignment.c

/*
 * Based on: Jonker, R., & Volgenant, A. (1987). <i>A shortest augmenting path
 * algorithm for dense and sparse linear assignment problems</i>. Computing,
 * 38(4), 325-340.
 */
#include "cache.h"
#include "linear-assignment.h"

#define COST(column, row) cost[(column) + column_count * (row)]

/*
 * The parameter `cost` is the cost matrix: the cost to assign column j to row
 * i is `cost[j + column_count * i].
 */
void compute_assignment(int column_count, int row_count, int *cost,
			int *column2row, int *row2column)
{
	int *v, *d;
	int *free_row, free_count = 0, saved_free_count, *pred, *col;
	int i, j, phase;

	if (column_count < 2) {
		memset(column2row, 0, sizeof(int) * column_count);
		memset(row2column, 0, sizeof(int) * row_count);
		return;
	}

	memset(column2row, -1, sizeof(int) * column_count);
	memset(row2column, -1, sizeof(int) * row_count);
	ALLOC_ARRAY(v, column_count);

	/* column reduction */
	for (j = column_count - 1; j >= 0; j--) {
		int i1 = 0;

		for (i = 1; i < row_count; i++)
			if (COST(j, i1) > COST(j, i))
				i1 = i;
		v[j] = COST(j, i1);
		if (row2column[i1] == -1) {
			/* row i1 unassigned */
			row2column[i1] = j;
			column2row[j] = i1;
		} else {
			if (row2column[i1] >= 0)
				row2column[i1] = -2 - row2column[i1];
			column2row[j] = -1;
		}
	}

	/* reduction transfer */
	ALLOC_ARRAY(free_row, row_count);
	for (i = 0; i < row_count; i++) {
		int j1 = row2column[i];
		if (j1 == -1)
			free_row[free_count++] = i;
		else if (j1 < -1)
			row2column[i] = -2 - j1;
		else {
			int min = COST(!j1, i) - v[!j1];
			for (j = 1; j < column_count; j++)
				if (j != j1 && min > COST(j, i) - v[j])
					min = COST(j, i) - v[j];
			v[j1] -= min;
		}
	}

	if (free_count ==
	    (column_count < row_count ? row_count - column_count : 0)) {
		free(v);
		free(free_row);
		return;
	}

	/* augmenting row reduction */
	for (phase = 0; phase < 2; phase++) {
		int k = 0;

		saved_free_count = free_count;
		free_count = 0;
		while (k < saved_free_count) {
			int u1, u2;
			int j1 = 0, j2, i0;

			i = free_row[k++];
			u1 = COST(j1, i) - v[j1];
			j2 = -1;
			u2 = INT_MAX;
			for (j = 1; j < column_count; j++) {
				int c = COST(j, i) - v[j];
				if (u2 > c) {
					if (u1 < c) {
						u2 = c;
						j2 = j;
					} else {
						u2 = u1;
						u1 = c;
						j2 = j1;
						j1 = j;
					}
				}
			}
			if (j2 < 0) {
				j2 = j1;
				u2 = u1;
			}

			i0 = column2row[j1];
			if (u1 < u2)
				v[j1] -= u2 - u1;
			else if (i0 >= 0) {
				j1 = j2;
				i0 = column2row[j1];
			}

			if (i0 >= 0) {
				if (u1 < u2)
					free_row[--k] = i0;
				else
					free_row[free_count++] = i0;
			}
			row2column[i] = j1;
			column2row[j1] = i;
		}
	}

	/* augmentation */
	saved_free_count = free_count;
	ALLOC_ARRAY(d, column_count);
	ALLOC_ARRAY(pred, column_count);
	ALLOC_ARRAY(col, column_count);
	for (free_count = 0; free_count < saved_free_count; free_count++) {
		int i1 = free_row[free_count], low = 0, up = 0, last, k;
		int min, c, u1;

		for (j = 0; j < column_count; j++) {
			d[j] = COST(j, i1) - v[j];
			pred[j] = i1;
			col[j] = j;
		}

		j = -1;
		do {
			last = low;
			min = d[col[up++]];
			for (k = up; k < column_count; k++) {
				j = col[k];
				c = d[j];
				if (c <= min) {
					if (c < min) {
						up = low;
						min = c;
					}
					col[k] = col[up];
					col[up++] = j;
				}
			}
			for (k = low; k < up; k++)
				if (column2row[col[k]] == -1)
					goto update;

			/* scan a row */
			do {
				int j1 = col[low++];

				i = column2row[j1];
				u1 = COST(j1, i) - v[j1] - min;
				for (k = up; k < column_count; k++) {
					j = col[k];
					c = COST(j, i) - v[j] - u1;
					if (c < d[j]) {
						d[j] = c;
						pred[j] = i;
						if (c == min) {
							if (column2row[j] == -1)
								goto update;
							col[k] = col[up];
							col[up++] = j;
						}
					}
				}
			} while (low != up);
		} while (low == up);

update:
		/* updating of the column pieces */
		for (k = 0; k < last; k++) {
			int j1 = col[k];
			v[j1] += d[j1] - min;
		}

		/* augmentation */
		do {
			if (j < 0)
				BUG("negative j: %d", j);
			i = pred[j];
			column2row[j] = i;
			SWAP(j, row2column[i]);
		} while (i1 != i);
	}

	free(col);
	free(pred);
	free(d);
	free(v);
	free(free_row);
}
Commit	Line	Data
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	*/
	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
e467a90c TG	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
22d87333 JS	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	}