diff options
| -rw-r--r-- | hungarian.c | 706 |
1 files changed, 706 insertions, 0 deletions
diff --git a/hungarian.c b/hungarian.c new file mode 100644 index 0000000..7171cf3 --- /dev/null +++ b/hungarian.c @@ -0,0 +1,706 @@ +#include <stdio.h> +#include <stdlib.h> + + +#define cell long +#define CELL_STR "%li" + +#define llong long long +#define byte char +#define boolean long +#define null 0 +#define true 1 +#define false 0 + + + +//new cell[X] +#define new_cells(X) new_longs(X) + +//new boolean[X] +#define new_booleans(X) new_longs(X) + +//new byte[X] +#define new_bytes(X) malloc(X) + +//new llong[X] +#define new_llongs(X) malloc((X) << 3) + +//new long[X] +#if !(defined __LP64__ || defined __LLP64__) + #define new_longs(X) malloc((X) << 2) /*32-bit*/ +#else + #define new_longs(X) malloc((X) << 3) /*64-bit*/ +#endif + +//new float[X] +#define new_floats(X) malloc((X) << 2) + +//new double[X] +#define new_doubles(X) malloc((X) << 3) + +//new ?[][X] +#define new_arrays(X) new_longs(X) + + +#ifdef DEBUG + #define debug(X) fprintf(stderr, "\e[31m%s\e[m\n", X) +#else + #define debug(X) +#endif + + + + +/** + * Cell marking: none + */ +#define UNMARKED 0L + +/** + * Cell marking: marked + */ +#define MARKED 1L + +/** + * Cell marking: prime + */ +#define PRIME 2L + + + +/** + * Bit set, a set of fixed number of bits/booleans + */ +typedef struct +{ + /** + * The set of all limbs, a limb consist of 64 bits + */ + cell* limbs; + + /** + * Singleton array with the index of the first non-zero limb + */ + long* first; + + /** + * Array the the index of the previous non-zero limb for echo limb + */ + long* prev; + + /** + * Array the the index of the next non-zero limb for echo limb + */ + long* next; + +} BitSet; + + + +long** kuhn_match(cell** table, long n, long m); +void kuhn_reduceRows(cell** t, long n, long m); +byte** kuhn_mark(cell** t, long n, long m); +boolean kuhn_isDone(byte** marks, boolean* colCovered, long n, long m); +long* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, long n, long m); +void kuhn_altMarks(byte** marks, long* altRow, long* altCol, long* colMarks, long* rowPrimes, long* prime, long n, long m); +void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, long n, long m); +long** kuhn_assign(byte** marks, long n, long m); + +BitSet new_BitSet(long size); +void BitSet_set(BitSet this, long i); +void BitSet_unset(BitSet this, long i); +long BitSet_any(BitSet this); + +long lb(llong x); + + + +void print(cell** t, long n, long m, long** assignment); + +int main(int argc, char** argv) +{ + unsigned a, d; + asm("cpuid"); + asm volatile("rdtsc" : "=a" (a), "=d" (d)); + srand(((llong)a) | (((llong)d) << 32LL)); + + + long n = 10, m = 15; + + long i; + cell** t = new_arrays(n); + cell** table = new_arrays(n); + long j; + if (argc < 2) + for (i = 0; i < n; i++) + { + *(t + i) = new_llongs(m); + *(table + i) = new_llongs(m); + for (j = 0; j < m; j++) + *(*(table + i) + j) = *(*(t + i) + j) = (cell)(random() & 63); + } + else + { + long x; + scanf("%li", &n); + scanf("%li", &m); + t = new_arrays(n); + table = new_arrays(n); + for (i = 0; i < n; i++) + { + *(t + i) = new_llongs(m); + *(table + i) = new_llongs(m); + for (j = 0; j < m; j++) + { + scanf(CELL_STR, &x); + *(*(table + i) + j) = *(*(t + i) + j) = x; + } + } + } + + //kuhn_match(table, n, m); + + printf("\nInput:\n\n"); + print(t, n, m, null); + + long** assignment = kuhn_match(table, n, m); + printf("\nOutput:\n\n"); + print(t, n, m, assignment); + + cell sum = 0; + for (i = 0; i < n; i++) + sum += *(*(t + *(*(assignment + i) + 0)) + *(*(assignment + i) + 1)); + printf("\n\nSum: %li\n\n", sum); + + return 0; +} + +void print(cell** t, long n, long m, long** assignment) +{ + long i, j; + + long** assigned = new_arrays(n); + for (i = 0; i < n; i++) + { + *(assigned + i) = new_longs(m); + for (j = 0; j < m; j++) + *(*(assigned + i) + j) = 0; + } + if (assignment != null) + for (i = 0; i < n; i++) + (*(*(assigned + **(assignment + i)) + *(*(assignment + i) + 1)))++; + + for (i = 0; i < n; i++) + { + printf(" "); + for (j = 0; j < m; j++) + { + if (*(*(assigned + i) + j)) + printf("\e[%lim", 30 + *(*(assigned + i) + j)); + printf("%5li%s\e[m ", (cell)(*(*(t + i) + j)), (*(*(assigned + i) + j) ? "^" : " ")); + } + printf("\n\n"); + } +} + + + +/** + * Calculates an optimal bipartite minimum weight matching using an + * O(n³)-time implementation of The Hungarian Algorithm, also known + * as Kuhn's Algorithm. This implemention is restricted to square + * tables. + * + * @param table The table in which to perform the matching + * @param n The dimension of the table + * @return The optimal assignment, an array of row–coloumn pairs + */ +long** kuhn_match(cell** table, long n, long m) +{ + long i; + + /* not copying table since it will only be used once */ + + kuhn_reduceRows(table, n, m); + byte** marks = kuhn_mark(table, n, m); + + boolean* rowCovered = new_booleans(n); + boolean* colCovered = new_booleans(m); + for (i = 0; i < n; i++) + { + *(rowCovered + i) = false; + *(colCovered + i) = false; + } + for (i = n; i < m; i++) + *(colCovered + i) = false; + + long* altRow = new_longs(n * m); + long* altCol = new_longs(n * m); + + long* rowPrimes = new_longs(n); + long* colMarks = new_longs(m); + + long* prime; + + for (;;) + { + if (kuhn_isDone(marks, colCovered, n, m)) + break; + + for (;;) + { + prime = kuhn_findPrime(table, marks, rowCovered, colCovered, n, m); + if (prime != null) + { + kuhn_altMarks(marks, altRow, altCol, colMarks, rowPrimes, prime, n, m); + for (i = 0; i < n; i++) + { + *(rowCovered + i) = false; + *(colCovered + i) = false; + } + for (i = n; i < m; i++) + *(colCovered + i) = false; + break; + } + kuhn_addAndSubtract(table, rowCovered, colCovered, n, m); + } + } + + return kuhn_assign(marks, n, m); +} + + +/** + * Reduces the values on each rows so that, for each row, the + * lowest cells value is zero, and all cells' values is decrease + * with the same value [the minium value in the row]. + * + * @param t The table in which to perform the reduction + * @param n The table's height + * @param m The table's width + */ +void kuhn_reduceRows(cell** t, long n, long m) +{ + long i, j; + cell min; + cell* ti; + for (i = 0; i < n; i++) + { + ti = *(t + i); + min = *ti; + for (j = 1; j < m; j++) + if (min > *(ti + j)) + min = *(ti + j); + + for (j = 0; j < m; j++) + *(ti + j) -= min; + } +} + + +/** + * Create a matrix with marking of cells in the table whose + * value is zero [minimal for the row]. Each marking will + * be on an unique row and an unique column. + * + * @param t The table in which to perform the reduction + * @param n The table's height + * @param m The table's width + * @return A matrix of markings as described in the summary + */ +byte** kuhn_mark(cell** t, long n, long m) +{ + long i, j; + byte** marks = new_arrays(n); + byte* marksi; + for (i = 0; i < n; i++) + { + marksi = *(marks + i) = new_bytes(m); + for (j = 0; j < m; j++) + *(marksi + j) = UNMARKED; + } + + boolean* rowCovered = new_booleans(n); + boolean* colCovered = new_booleans(m); + for (i = 0; i < n; i++) + { + *(rowCovered + i) = false; + *(colCovered + i) = false; + } + for (i = 0; i < m; i++) + *(colCovered + i) = false; + + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + if ((!*(rowCovered + i)) && (!*(colCovered + j)) && (*(*(t + i) + j) == 0)) + { + *(*(marks + i) + j) = MARKED; + *(rowCovered + i) = true; + *(colCovered + j) = true; + } + + return marks; +} + + +/** + * Determines whether the marking is complete, that is + * if each row has a marking which is on a unique column. + * + * @param marks The marking matrix + * @param colCovered An array which tells whether a column is covered + * @param n The table's height + * @param m The table's width + * @return Whether the marking is complete + */ +boolean kuhn_isDone(byte** marks, boolean* colCovered, long n, long m) +{ + long i, j; + for (j = 0; j < m; j++) + for (i = 0; i < n; i++) + if (*(*(marks + i) + j) == MARKED) + { + *(colCovered + j) = true; + break; + } + + long count = 0; + for (j = 0; j < m; j++) + if (*(colCovered + j)) + count++; + + return count == n; +} + + +/** + * Finds a prime + * + * @param t The table + * @param marks The marking matrix + * @param rowCovered Row cover array + * @param colCovered Column cover array + * @param n The table's height + * @param m The table's width + * @return The row and column of the found print, <code>null</code> will be returned if none can be found + */ +long* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, long n, long m) +{ + long i, j; + BitSet zeroes = new_BitSet(n * n); + + for (i = 0; i < n; i++) + if (!*(rowCovered + i)) + for (j = 0; j < m; j++) + if ((!*(colCovered + j)) && (*(*(t + i) + j) == 0)) + BitSet_set(zeroes, i * m + j); + + long p, row, col; + boolean markInRow; + + for (;;) + { + p = BitSet_any(zeroes); + if (p < 0) + return null; + + row = p / m; + col = p % m; + + *(*(marks + row) + col) = PRIME; + + markInRow = false; + for (j = 0; j < m; j++) + if (*(*(marks + row) + j) == MARKED) + { + markInRow = true; + col = j; + } + + if (markInRow) + { + *(rowCovered + row) = true; + *(colCovered + col) = false; + + for (i = 0; i < n; i++) + if ((*(*(t + i) + col) == 0) && (row != i)) + { + if ((!*(rowCovered + i)) && (!*(colCovered + col))) + BitSet_set(zeroes, i * m + col); + else + BitSet_unset(zeroes, i * m + col); + } + + for (j = 0; j < m; j++) + if ((*(*(t + row) + j) == 0) && (col != j)) + { + if ((!*(rowCovered + row)) && (!*(colCovered + j))) + BitSet_set(zeroes, row * m + j); + else + BitSet_unset(zeroes, row * m + j); + } + + if ((!*(rowCovered + row)) && (!*(colCovered + col))) + BitSet_set(zeroes, row * m + col); + else + BitSet_unset(zeroes, row * m + col); + } + else + { + long* rc = new_longs(2); + *rc = row; + *(rc + 1) = col; + return rc; + } + } +} + + +/** + * Removes all prime marks and modifies the marking + * + * @param marks The marking matrix + * @param altRow Marking modification path rows + * @param altCol Marking modification path columns + * @param colMarks Markings in the columns + * @param rowPrimes Primes in the rows + * @param prime The last found prime + * @param n The table's height + * @param m The table's width + */ +void kuhn_altMarks(byte** marks, long* altRow, long* altCol, long* colMarks, long* rowPrimes, long* prime, long n, long m) +{ + long index = 0, i, j; + *altRow = *prime; + *altCol = *(prime + 1); + + for (i = 0; i < n; i++) + { + *(rowPrimes + i) = -1; + *(colMarks + i) = -1; + } + for (i = n; i < m; i++) + *(colMarks + i) = -1; + + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + if (*(*(marks + i) + j) == MARKED) + *(colMarks + j) = i; + else if (*(*(marks + i) + j) == PRIME) + *(rowPrimes + i) = j; + + long row, col; + for (;;) + { + row = *(colMarks + *(altCol + index)); + if (row < 0) + break; + + index++; + *(altRow + index) = row; + *(altCol + index) = *(altCol + index - 1); + + col = *(rowPrimes + *(altRow + index)); + + index++; + *(altRow + index) = *(altRow + index - 1); + *(altCol + index) = col; + } + + byte* markx; + for (i = 0; i <= index; i++) + { + markx = *(marks + *(altRow + i)) + *(altCol + i); + if (*markx == MARKED) + *markx = UNMARKED; + else + *markx = MARKED; + } + + byte* marksi; + for (i = 0; i < n; i++) + { + marksi = *(marks + i); + for (j = 0; j < m; j++) + if (*(marksi + j) == PRIME) + *(marksi + j) = UNMARKED; + } +} + + +/** + * Depending on whether the cells' rows and columns are covered, + * the the minimum value in the table is added, subtracted or + * neither from the cells. + * + * @param t The table to manipulate + * @param rowCovered Array that tell whether the rows are covered + * @param colCovered Array that tell whether the columns are covered + * @param n The table's height + * @param m The table's width + */ +void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, long n, long m) +{ + long i, j; + cell min = 0x7FFFffffL; + for (i = 0; i < n; i++) + if (!*(rowCovered + i)) + for (j = 0; j < m; j++) + if ((!*(colCovered + j)) && (min > *(*(t + i) + j))) + min = *(*(t + i) + j); + + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + { + if (*(rowCovered + i)) + *(*(t + i) + j) += min; + if (*(colCovered + j) == false) + *(*(t + i) + j) -= min; + } +} + + +/** + * Creates a list of the assignment cells + * + * @param marks Matrix markings + * @param n The table's height + * @param m The table's width + * @return The assignment, an array of row–coloumn pairs + */ +long** kuhn_assign(byte** marks, long n, long m) +{ + long** assignment = new_arrays(n); + + long i, j; + for (i = 0; i < n; i++) + { + *(assignment + i) = new_longs(2); + for (j = 0; j < m; j++) + if (*(*(marks + i) + j) == MARKED) + { + **(assignment + i) = i; + *(*(assignment + i) + 1) = j; + } + } + + return assignment; +} + + +/** + * Constructor for BitSet + * + * @param size The (fixed) number of bits to bit set should contain + * @return The a unique BitSet instance with the specified size + */ +BitSet new_BitSet(long size) +{ + BitSet this; + + long c = size >> 6L; + if (size & 63L) + c++; + + this.limbs = new_llongs(c); + this.prev = new_longs(c + 1L); + this.next = new_longs(c + 1L); + *(this.first = new_longs(1)) = 0; + + long i; + for (i = 0; i < c; i++) + { + *(this.limbs + i) = 0LL; + *(this.prev + i) = *(this.next + i) = 0L; + } + *(this.prev + c) = *(this.next + c) = 0L; + + return this; +} + +/** + * Turns on a bit in a bit set + * + * @param this The bit set + * @param i The index of the bit to turn on + */ +void BitSet_set(BitSet this, long i) +{ + long j = i >> 6L; + llong old = *(this.limbs + j); + + *(this.limbs + j) |= 1LL << (llong)(i & 63L); + + if ((!*(this.limbs + j)) ^ (!old)) + { + j++; + *(this.prev + *(this.first)) = j; + *(this.prev + j) = 0; + *(this.next + j) = *(this.first); + *(this.first) = j; + } +} + +/** + * Turns off a bit in a bit set + * + * @param this The bit set + * @param i The index of the bit to turn off + */ +void BitSet_unset(BitSet this, long i) +{ + long j = i >> 6L; + llong old = *(this.limbs + j); + + *(this.limbs + j) &= ~(1LL << (llong)(i & 63L)); + + if ((!*(this.limbs + j)) ^ (!old)) + { + j++; + long p = *(this.prev + j); + long n = *(this.next + j); + *(this.prev + n) = p; + *(this.next + p) = n; + if (*(this.first) == j) + *(this.first) = *(this.next + j); + } +} + +/** + * Gets the index of any set bit in a bit set + * + * @param this The bit set + * @return The index of any set bit + */ +long BitSet_any(BitSet this) +{ + if (*(this.first) == 0L) + return -1; + + long i = *(this.first) - 1L; + return lb(*(this.limbs + i) & -*(this.limbs + i)) + (i << 6L); +} + + +/** + * Calculates the floored binary logarithm of a positive integer + * + * @param value The integer whose logarithm to calculate + * @return The floored binary logarithm of the integer + */ +long lb(llong value) +{ + long rc = 0L; + llong v = value; + + if (v & 0xFFFFFFFF00000000LL) { rc |= 32L; v >>= 32LL; } + if (v & 0x00000000FFFF0000LL) { rc |= 16L; v >>= 16LL; } + if (v & 0x000000000000FF00LL) { rc |= 8L; v >>= 8LL; } + if (v & 0x00000000000000F0LL) { rc |= 4L; v >>= 4LL; } + if (v & 0x000000000000000CLL) { rc |= 2L; v >>= 2LL; } + if (v & 0x0000000000000002LL) rc |= 1L; + + return rc; +} + |