hungarian.c 11.2 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
/**
 * @file
 * @brief   Solving the Minimum Assignment Problem using the Hungarian Method.
 */
Matthias Braun's avatar
Matthias Braun committed
31
#include "config.h"
32

33
34
35
36
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

37
#include "util.h"
38
39
#include "xmalloc.h"
#include "debug.h"
40
#include "error.h"
41
#include "hungarian.h"
42
#include "raw_bitset.h"
43

44
DEBUG_ONLY(static firm_dbg_module_t *dbg;)
45

46
struct hungarian_problem_t {
47
48
49
50
51
52
53
54
55
	unsigned      num_rows;      /**< number of rows */
	unsigned      num_cols;      /**< number of columns */
	unsigned     *cost;          /**< the cost matrix */
	unsigned      max_cost;      /**< the maximal costs in the matrix */
	match_type_t  match_type;    /**< PERFECT or NORMAL matching */
	unsigned     *missing_left;  /**< bitset: left side nodes having no edge to
	                                  the right side */
	unsigned     *missing_right; /**< bitset: right side nodes having no edge to
	                              the left side */
56
57
};

58
59
static void hungarian_dump_f(FILE *f, const unsigned *cost,
                             unsigned num_rows, unsigned num_cols, int width)
60
{
61
	unsigned r, c;
62
63

	fprintf(f , "\n");
64
	for (r = 0; r < num_rows; r++) {
65
		fprintf(f, " [");
66
67
		for (c = 0; c < num_cols; c++) {
			fprintf(f, "%*u", width, cost[r*num_cols + c]);
68
69
70
71
72
73
		}
		fprintf(f, "]\n");
	}
	fprintf(f, "\n");
}

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

79
hungarian_problem_t *hungarian_new(unsigned num_rows, unsigned num_cols,
80
                                   match_type_t match_type)
81
{
82
	hungarian_problem_t *p = XMALLOCZ(hungarian_problem_t);
83

84
	FIRM_DBG_REGISTER(dbg, "firm.hungarian");
85
86
87
88
89

	/*
		Is the number of cols  not equal to number of rows ?
		If yes, expand with 0 - cols / 0 - cols
	*/
90
91
	num_rows = MAX(num_cols, num_rows);
	num_cols = num_rows;
92

93
94
	p->num_rows   = num_rows;
	p->num_cols   = num_cols;
95
96
97
98
99
100
101
102
	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) {
103
104
105
106
		p->missing_left  = rbitset_malloc(num_rows);
		p->missing_right = rbitset_malloc(num_cols);
		rbitset_set_all(p->missing_left,  num_rows);
		rbitset_set_all(p->missing_right, num_cols);
107
	}
108
109

	/* allocate space for cost matrix */
110
	p->cost = XMALLOCNZ(unsigned, num_rows * num_cols);
111
112
113
	return p;
}

114
115
void hungarian_prepare_cost_matrix(hungarian_problem_t *p,
                                   hungarian_mode_t mode)
116
{
117
	if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
118
119
120
121
		unsigned  r, c;
		unsigned  num_cols = p->num_cols;
		unsigned *cost     = p->cost;
		unsigned  max_cost = p->max_cost;
122
123
		for (r = 0; r < p->num_rows; r++) {
			for (c = 0; c < p->num_cols; c++) {
124
				cost[r*num_cols + c] = max_cost - cost[r*num_cols + c];
125
126
			}
		}
127
	} else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
128
		/* nothing to do */
129
	} else {
130
		panic("Unknown hungarian problem mode");
131
132
133
	}
}

134
void hungarian_add(hungarian_problem_t *p, unsigned left, unsigned right,
135
                   unsigned cost)
136
{
137
138
139
	assert(p->num_rows > left  && "Invalid row selected.");
	assert(p->num_cols > right && "Invalid column selected.");

140
141
	p->cost[left*p->num_cols + right] = cost;
	p->max_cost                       = MAX(p->max_cost, cost);
142
143

	if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
144
145
		rbitset_clear(p->missing_left, left);
		rbitset_clear(p->missing_right, right);
146
	}
147
148
}

149
void hungarian_remove(hungarian_problem_t *p, unsigned left, unsigned right)
150
{
151
152
153
	assert(p->num_rows > left  && "Invalid row selected.");
	assert(p->num_cols > right && "Invalid column selected.");

154
	p->cost[left*p->num_cols + right] = 0;
155
156

	if (p->match_type == HUNGARIAN_MATCH_NORMAL) {
157
158
		rbitset_set(p->missing_left, left);
		rbitset_set(p->missing_right, right);
159
	}
160
161
}

162
163
void hungarian_free(hungarian_problem_t* p)
{
164
165
166
	xfree(p->missing_left);
	xfree(p->missing_right);
	xfree(p->cost);
167
168
169
	xfree(p);
}

170
int hungarian_solve(hungarian_problem_t* p, unsigned *assignment,
171
                    unsigned *final_cost, unsigned cost_threshold)
172
{
173
174
175
176
177
178
179
180
181
182
183
184
	unsigned  res_cost     = 0;
	unsigned  num_rows     = p->num_rows;
	unsigned  num_cols     = p->num_cols;
	unsigned *cost         = p->cost;
	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);
185
186
187
188
189
190
	unsigned  r;
	unsigned  c;
	unsigned  t;
	unsigned  unmatched;

	memset(assignment, -1, num_rows * sizeof(assignment[0]));
191
192

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

195
	for (c = 0; c < num_cols; ++c) {
196
		unsigned col_mininum = cost[0*num_cols + c];
197

198
		for (r = 1; r < num_rows; ++r) {
199
200
			if (cost[r*num_cols + c] < col_mininum)
				col_mininum = cost[r*num_cols + c];
201
202
		}

203
		if (col_mininum == 0)
204
			continue;
205

206
		res_cost += col_mininum;
207
		for (r = 0; r < num_rows; ++r)
208
			cost[r*num_cols + c] -= col_mininum;
209
210
211
212
	}
	/* End subtract column minima in order to start with lots of zeros 12 */

	/* Begin initial state 16 */
213
	unmatched = 0;
214
215
216
217
218
	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;
219
220
	}

221
	for (r = 0; r < num_rows; ++r) {
222
		unsigned row_minimum = cost[r*num_cols + 0];
223

224
		for (c = 1; c < num_cols; ++c) {
225
226
			if (cost[r*num_cols + c] < row_minimum)
				row_minimum = cost[r*num_cols + c];
227
228
		}

229
		row_dec[r] = row_minimum;
230

231
		for (c = 0; c < num_cols; ++c) {
232
233
234
235
236
237
238
239
240
			if (cost[r*num_cols + c] != row_minimum)
				continue;
			if (row_mate[c] != (unsigned)-1)
				continue;

			col_mate[r] = c;
			row_mate[c] = r;
			DBG((dbg, LEVEL_1, "matching col %u == row %u\n", c, r));
			goto row_done;
241
242
		}

243
		col_mate[r] = (unsigned)-1;
244
245
		DBG((dbg, LEVEL_1, "node %u: unmatched row %u\n", unmatched, r));
		unchosen_row[unmatched++] = r;
246
247
248
249
250
row_done: ;
	}
	/* End initial state 16 */

	/* Begin Hungarian algorithm 18 */
251
	if (unmatched == 0)
252
253
		goto done;

254
	t = unmatched;
255
	for (;;) {
256
257
		unsigned q = 0;
		unsigned j;
258
		DBG((dbg, LEVEL_1, "Matched %u rows.\n", num_rows - t));
259

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

267
268
				for (c = 0; c < num_cols; ++c) {
					if (slack[c]) {
269
						int del = cost[r*num_cols + c] - s + col_inc[c];
270

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

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

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

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

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

335
			DBG((dbg, LEVEL_1, "rematching col %u == row %u\n", c, r));
336
			if (j == (unsigned)-1)
337
338
				break;

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

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

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

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

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

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

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

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

	/* End Hungarian algorithm 18 */

	/* collect the assigned values */
391
	for (r = 0; r < num_rows; ++r) {
392
393
		if (cost_threshold > 0
		    && cost[r*num_cols + col_mate[r]] >= cost_threshold)
394
			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
		for (r = 0; r < num_rows; ++r) {
402
403
			if (rbitset_is_set(p->missing_left, r)
			        || rbitset_is_set(p->missing_right, col_mate[r]))
404
				assignment[r] = -1;
405
406
407
		}
	}

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

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

417
	for (c = 0; c < num_cols; ++c)
418
		res_cost -= col_inc[c];
419

420
	DBG((dbg, LEVEL_1, "Cost is %d\n", res_cost));
421
422
423
424
425
426
427
428
429
430

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

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

	return 0;
435
}