hungarian.c 10.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
/********************************************************************
 ********************************************************************
 **
 ** libhungarian by Cyrill Stachniss, 2004
 **
 ** Added and adapted to libFirm by Christian Wuerdig, 2006
 **
 ** Solving the Minimum Assignment Problem using the
 ** Hungarian Method.
 **
 ** ** This file may be freely copied and distributed! **
 **
 ** Parts of the used code was originally provided by the
 ** "Stanford GraphGase", but I made changes to this code.
 ** As asked by  the copyright node of the "Stanford GraphGase",
 ** I hereby proclaim that this file are *NOT* part of the
 ** "Stanford GraphGase" distrubition!
 **
 ** This file is distributed in the hope that it will be useful,
 ** but WITHOUT ANY WARRANTY; without even the implied
 ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 ** PURPOSE.
 **
 ********************************************************************
 ********************************************************************/

27
28
29
30
31
/**
 * @file
 * @brief   Solving the Minimum Assignment Problem using the Hungarian Method.
 * @version $Id$
 */
Matthias Braun's avatar
Matthias Braun committed
32
#include "config.h"
33

34
35
36
37
38
39
40
41
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#include "irtools.h"
#include "xmalloc.h"
#include "debug.h"
#include "obst.h"
42
#include "bitset.h"
43
#include "error.h"
44
45
46

#include "hungarian.h"

47
48
49
struct hungarian_problem_t {
	unsigned num_rows;          /**< number of rows */
	unsigned num_cols;          /**< number of columns */
50
51
52
53
54
	int      **cost;            /**< the cost matrix */
	int      max_cost;          /**< the maximal costs in the matrix */
	int      match_type;        /**< PERFECT or NORMAL matching */
	bitset_t *missing_left;     /**< left side nodes having no edge to the right side */
	bitset_t *missing_right;    /**< right side nodes having no edge to the left side */
55
56
57
58
	struct obstack obst;
	DEBUG_ONLY(firm_dbg_module_t *dbg);
};

59
60
static void hungarian_dump_f(FILE *f, int **C, unsigned n_rows, unsigned n_cols,
                             int width)
61
{
62
	unsigned r;
63
64

	fprintf(f , "\n");
65
66
	for (r = 0; r < n_rows; r++) {
		unsigned c;
67
		fprintf(f, " [");
68
69
		for (c = 0; c < n_cols; c++) {
			fprintf(f, "%*d", width, C[r][c]);
70
71
72
73
74
75
		}
		fprintf(f, "]\n");
	}
	fprintf(f, "\n");
}

76
77
void hungarian_print_cost_matrix(hungarian_problem_t *p, int width)
{
Matthias Braun's avatar
Matthias Braun committed
78
	hungarian_dump_f(stderr, p->cost, p->num_rows, p->num_cols, width);
79
80
}

81
82
hungarian_problem_t *hungarian_new(unsigned n_rows, unsigned n_cols,
                                   match_type_t match_type)
83
{
84
	unsigned r;
85
	hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
86
87
88
89
90
91
92

	FIRM_DBG_REGISTER(p->dbg, "firm.hungarian");

	/*
		Is the number of cols  not equal to number of rows ?
		If yes, expand with 0 - cols / 0 - cols
	*/
93
94
	n_rows = MAX(n_cols, n_rows);
	n_cols = n_rows;
95
96
97

	obstack_init(&p->obst);

98
99
	p->num_rows   = n_rows;
	p->num_cols   = n_cols;
100
101
102
103
104
105
106
107
	p->match_type = match_type;

	/*
		In case of normal matching, we have to keep
		track of nodes without edges to kill them in
		the assignment later.
	*/
	if (match_type == HUNGARIAN_MATCH_NORMAL) {
108
109
		p->missing_left  = bitset_obstack_alloc(&p->obst, n_rows);
		p->missing_right = bitset_obstack_alloc(&p->obst, n_cols);
110
111
112
		bitset_set_all(p->missing_left);
		bitset_set_all(p->missing_right);
	}
113
114

	/* allocate space for cost matrix */
115
116
117
	p->cost = OALLOCNZ(&p->obst, int*, n_rows);
	for (r = 0; r < p->num_rows; r++)
		p->cost[r] = OALLOCNZ(&p->obst, int, n_cols);
118
119
120
121

	return p;
}

122
123
void hungarian_prepare_cost_matrix(hungarian_problem_t *p,
                                   hungarian_mode_t mode)
124
{
125
	unsigned r, c;
126
127

	if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
128
129
130
		for (r = 0; r < p->num_rows; r++) {
			for (c = 0; c < p->num_cols; c++) {
				p->cost[r][c] = p->max_cost - p->cost[r][c];
131
132
			}
		}
133
	} else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
134
		/* nothing to do */
135
136
	} else {
		panic("Unknown hungarian problem mode\n");
137
138
139
	}
}

140
141
void hungarian_add(hungarian_problem_t *p, unsigned left, unsigned right,
                   int cost)
142
{
143
144
	assert(p->num_rows > left  && "Invalid row selected.");
	assert(p->num_cols > right && "Invalid column selected.");
145
	assert(cost >= 0);
146
147
148

	p->cost[left][right] = cost;
	p->max_cost          = MAX(p->max_cost, cost);
149
150
151
152
153

	if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
		bitset_clear(p->missing_left, left);
		bitset_clear(p->missing_right, right);
	}
154
155
}

156
void hungarian_remove(hungarian_problem_t *p, unsigned left, unsigned right)
157
{
158
159
160
	assert(p->num_rows > left  && "Invalid row selected.");
	assert(p->num_cols > right && "Invalid column selected.");

161
	/* Set cost[left][right] to 0. */
162
	p->cost[left][right] = 0;
163
164
165
166
167

	if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
		bitset_set(p->missing_left, left);
		bitset_set(p->missing_right, right);
	}
168
169
}

170
171
void hungarian_free(hungarian_problem_t* p)
{
172
173
174
175
	obstack_free(&p->obst, NULL);
	xfree(p);
}

176
177
int hungarian_solve(hungarian_problem_t* p, unsigned *assignment,
                    int *final_cost, int cost_threshold)
178
{
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
	int       cost          = 0;
	unsigned  num_rows      = p->num_rows;
	unsigned  num_cols      = p->num_cols;
	unsigned *col_mate      = XMALLOCNZ(unsigned, num_rows);
	unsigned *row_mate      = XMALLOCNZ(unsigned, num_cols);
	unsigned *parent_row    = XMALLOCNZ(unsigned, num_cols);
	unsigned *unchosen_row  = XMALLOCNZ(unsigned, num_rows);
	int      *row_dec       = XMALLOCNZ(int, num_rows);
	int      *col_inc       = XMALLOCNZ(int, num_cols);
	int      *slack         = XMALLOCNZ(int, num_cols);
	unsigned *slack_row     = XMALLOCNZ(unsigned, num_rows);
	unsigned  r;
	unsigned  c;
	unsigned  t;
	unsigned  unmatched;

	memset(assignment, -1, num_rows * sizeof(assignment[0]));
196
197
198
199

	/* Begin subtract column minima in order to start with lots of zeros 12 */
	DBG((p->dbg, LEVEL_1, "Using heuristic\n"));

200
201
	for (c = 0; c < num_cols; ++c) {
		int s = p->cost[0][c];
202

203
204
205
		for (r = 1; r < num_rows; ++r) {
			if (p->cost[r][c] < s)
				s = p->cost[r][c];
206
207
208
		}

		cost += s;
209
210
		if (s == 0)
			continue;
211

212
213
		for (r = 0; r < num_rows; ++r)
			p->cost[r][c] -= s;
214
215
216
217
218
	}
	/* End subtract column minima in order to start with lots of zeros 12 */

	/* Begin initial state 16 */
	t = 0;
219
220
221
222
223
	for (c = 0; c < num_cols; ++c) {
		row_mate[c]   = (unsigned) -1;
		parent_row[c] = (unsigned) -1;
		col_inc[c]    = 0;
		slack[c]      = INT_MAX;
224
225
	}

226
227
	for (r = 0; r < num_rows; ++r) {
		int s = p->cost[r][0];
228

229
230
231
		for (c = 1; c < num_cols; ++c) {
			if (p->cost[r][c] < s)
				s = p->cost[r][c];
232
233
		}

234
		row_dec[r] = s;
235

236
237
238
239
240
		for (c = 0; c < num_cols; ++c) {
			if (s == p->cost[r][c] && row_mate[c] == (unsigned)-1) {
				col_mate[r] = c;
				row_mate[c] = r;
				DBG((p->dbg, LEVEL_1, "matching col %d == row %d\n", c, r));
241
242
243
244
				goto row_done;
			}
		}

245
246
247
		col_mate[r] = (unsigned)-1;
		DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, r));
		unchosen_row[t++] = r;
248
249
250
251
252
253
254
255
256
row_done: ;
	}
	/* End initial state 16 */

	/* Begin Hungarian algorithm 18 */
	if (t == 0)
		goto done;

	unmatched = t;
257
	for (;;) {
258
259
260
		unsigned q = 0;
		unsigned j;
		DBG((p->dbg, LEVEL_1, "Matched %d rows.\n", num_rows - t));
261

262
		for (;;) {
263
			int s;
264
265
			while (q < t) {
				/* Begin explore node q of the forest 19 */
266
267
				r = unchosen_row[q];
				s = row_dec[r];
268

269
270
271
				for (c = 0; c < num_cols; ++c) {
					if (slack[c]) {
						int del = p->cost[r][c] - s + col_inc[c];
272

273
						if (del < slack[c]) {
274
							if (del == 0) {
275
								if (row_mate[c] == (unsigned)-1)
276
277
									goto breakthru;

278
279
280
281
282
283
284
								slack[c]      = 0;
								parent_row[c] = r;
								DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[c], c, r));
								unchosen_row[t++] = row_mate[c];
							} else {
								slack[c]     = del;
								slack_row[c] = r;
285
286
287
288
289
290
291
292
293
							}
						}
					}
				}
				/* End explore node q of the forest 19 */
				q++;
			}

			/* Begin introduce a new zero into the matrix 21 */
294
295
296
297
			s = INT_MAX;
			for (c = 0; c < num_cols; ++c) {
				if (slack[c] && slack[c] < s)
					s = slack[c];
298
299
300
301
302
			}

			for (q = 0; q < t; ++q)
				row_dec[unchosen_row[q]] += s;

303
304
305
306
			for (c = 0; c < num_cols; ++c) {
				if (slack[c]) {
					slack[c] -= s;
					if (slack[c] == 0) {
307
						/* Begin look at a new zero 22 */
308
309
310
311
						r = slack_row[c];
						DBG((p->dbg, LEVEL_1, "Decreasing uncovered elements by %d produces zero at [%d, %d]\n", s, r, c));
						if (row_mate[c] == (unsigned)-1) {
							for (j = c + 1; j < num_cols; ++j) {
312
313
314
315
								if (slack[j] == 0)
									col_inc[j] += s;
							}
							goto breakthru;
316
317
318
319
						} else {
							parent_row[c] = r;
							DBG((p->dbg, LEVEL_1, "node %d: row %d == col %d -- row %d\n", t, row_mate[c], c, r));
							unchosen_row[t++] = row_mate[c];
320
321
322
						}
						/* End look at a new zero 22 */
					}
323
324
				} else {
					col_inc[c] += s;
325
326
327
328
329
330
331
				}
			}
			/* End introduce a new zero into the matrix 21 */
		}
breakthru:
		/* Begin update the matching 20 */
		DBG((p->dbg, LEVEL_1, "Breakthrough at node %d of %d.\n", q, t));
332
		for (;;) {
333
334
335
			j           = col_mate[r];
			col_mate[r] = c;
			row_mate[c] = r;
336

337
338
			DBG((p->dbg, LEVEL_1, "rematching col %d == row %d\n", c, r));
			if (j == (unsigned)-1)
339
340
				break;

341
342
			r = parent_row[j];
			c = j;
343
344
345
346
347
348
349
350
		}
		/* End update the matching 20 */

		if (--unmatched == 0)
			goto done;

		/* Begin get ready for another stage 17 */
		t = 0;
351
352
353
		for (c = 0; c < num_cols; ++c) {
			parent_row[c] = -1;
			slack[c]      = INT_MAX;
354
355
		}

356
357
358
359
		for (r = 0; r < num_rows; ++r) {
			if (col_mate[r] == (unsigned)-1) {
				DBG((p->dbg, LEVEL_1, "node %d: unmatched row %d\n", t, r));
				unchosen_row[t++] = r;
360
361
362
363
364
365
366
			}
		}
		/* End get ready for another stage 17 */
	}
done:

	/* Begin double check the solution 23 */
367
368
369
	for (r = 0; r < num_rows; ++r) {
		for (c = 0; c < num_cols; ++c) {
			if (p->cost[r][c] < row_dec[r] - col_inc[c])
370
371
372
373
				return -1;
		}
	}

374
375
376
	for (r = 0; r < num_rows; ++r) {
		c = col_mate[r];
		if (c == (unsigned)-1 || p->cost[r][c] != row_dec[r] - col_inc[c])
377
378
379
			return -2;
	}

380
381
382
	for (r = c = 0; c < num_cols; ++c) {
		if (col_inc[c])
			r++;
383
384
	}

385
	if (r > num_rows)
386
387
388
389
390
391
		return -3;
	/* End double check the solution 23 */

	/* End Hungarian algorithm 18 */

	/* collect the assigned values */
392
393
394
	for (r = 0; r < num_rows; ++r) {
		if (cost_threshold > 0 && p->cost[r][col_mate[r]] >= cost_threshold)
			assignment[r] = -1; /* remove matching having cost > threshold */
Christian Würdig's avatar
Christian Würdig committed
395
		else
396
			assignment[r] = col_mate[r];
397
398
	}

399
400
	/* In case of normal matching: remove impossible ones */
	if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
401
402
403
404
		for (r = 0; r < num_rows; ++r) {
			if (bitset_is_set(p->missing_left, r)
			        || bitset_is_set(p->missing_right, col_mate[r]))
				assignment[r] = -1;
405
406
407
		}
	}

408
409
410
	for (r = 0; r < num_rows; ++r) {
		for (c = 0; c < num_cols; ++c) {
			p->cost[r][c] = p->cost[r][c] - row_dec[r] + col_inc[c];
411
412
413
		}
	}

414
415
	for (r = 0; r < num_rows; ++r)
		cost += row_dec[r];
416

417
418
	for (c = 0; c < num_cols; ++c)
		cost -= col_inc[c];
419
420
421
422
423
424
425
426
427
428
429
430

	DBG((p->dbg, LEVEL_1, "Cost is %d\n", cost));

	xfree(slack);
	xfree(col_inc);
	xfree(parent_row);
	xfree(row_mate);
	xfree(slack_row);
	xfree(row_dec);
	xfree(unchosen_row);
	xfree(col_mate);

431
432
	if (final_cost != NULL)
		*final_cost = cost;
433
434

	return 0;
435
}