THIS IS A TEST INSTANCE ONLY! REPOSITORIES CAN BE DELETED AT ANY TIME!

Git Source Code Mirror - This is a publish-only repository and all pull requests are ignored. Please follow Documentation/SubmittingPatches procedure for any of your improvements.
git
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

208 lines
4.1KB

  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. #define COST(column, row) cost[(column) + column_count * (row)]
  9. /*
  10. * The parameter `cost` is the cost matrix: the cost to assign column j to row
  11. * i is `cost[j + column_count * i].
  12. */
  13. void compute_assignment(int column_count, int row_count, int *cost,
  14. int *column2row, int *row2column)
  15. {
  16. int *v, *d;
  17. int *free_row, free_count = 0, saved_free_count, *pred, *col;
  18. int i, j, phase;
  19. if (column_count < 2) {
  20. memset(column2row, 0, sizeof(int) * column_count);
  21. memset(row2column, 0, sizeof(int) * row_count);
  22. return;
  23. }
  24. memset(column2row, -1, sizeof(int) * column_count);
  25. memset(row2column, -1, sizeof(int) * row_count);
  26. ALLOC_ARRAY(v, column_count);
  27. /* column reduction */
  28. for (j = column_count - 1; j >= 0; j--) {
  29. int i1 = 0;
  30. for (i = 1; i < row_count; i++)
  31. if (COST(j, i1) > COST(j, i))
  32. i1 = i;
  33. v[j] = COST(j, i1);
  34. if (row2column[i1] == -1) {
  35. /* row i1 unassigned */
  36. row2column[i1] = j;
  37. column2row[j] = i1;
  38. } else {
  39. if (row2column[i1] >= 0)
  40. row2column[i1] = -2 - row2column[i1];
  41. column2row[j] = -1;
  42. }
  43. }
  44. /* reduction transfer */
  45. ALLOC_ARRAY(free_row, row_count);
  46. for (i = 0; i < row_count; i++) {
  47. int j1 = row2column[i];
  48. if (j1 == -1)
  49. free_row[free_count++] = i;
  50. else if (j1 < -1)
  51. row2column[i] = -2 - j1;
  52. else {
  53. int min = COST(!j1, i) - v[!j1];
  54. for (j = 1; j < column_count; j++)
  55. if (j != j1 && min > COST(j, i) - v[j])
  56. min = COST(j, i) - v[j];
  57. v[j1] -= min;
  58. }
  59. }
  60. if (free_count ==
  61. (column_count < row_count ? row_count - column_count : 0)) {
  62. free(v);
  63. free(free_row);
  64. return;
  65. }
  66. /* augmenting row reduction */
  67. for (phase = 0; phase < 2; phase++) {
  68. int k = 0;
  69. saved_free_count = free_count;
  70. free_count = 0;
  71. while (k < saved_free_count) {
  72. int u1, u2;
  73. int j1 = 0, j2, i0;
  74. i = free_row[k++];
  75. u1 = COST(j1, i) - v[j1];
  76. j2 = -1;
  77. u2 = INT_MAX;
  78. for (j = 1; j < column_count; j++) {
  79. int c = COST(j, i) - v[j];
  80. if (u2 > c) {
  81. if (u1 < c) {
  82. u2 = c;
  83. j2 = j;
  84. } else {
  85. u2 = u1;
  86. u1 = c;
  87. j2 = j1;
  88. j1 = j;
  89. }
  90. }
  91. }
  92. if (j2 < 0) {
  93. j2 = j1;
  94. u2 = u1;
  95. }
  96. i0 = column2row[j1];
  97. if (u1 < u2)
  98. v[j1] -= u2 - u1;
  99. else if (i0 >= 0) {
  100. j1 = j2;
  101. i0 = column2row[j1];
  102. }
  103. if (i0 >= 0) {
  104. if (u1 < u2)
  105. free_row[--k] = i0;
  106. else
  107. free_row[free_count++] = i0;
  108. }
  109. row2column[i] = j1;
  110. column2row[j1] = i;
  111. }
  112. }
  113. /* augmentation */
  114. saved_free_count = free_count;
  115. ALLOC_ARRAY(d, column_count);
  116. ALLOC_ARRAY(pred, column_count);
  117. ALLOC_ARRAY(col, column_count);
  118. for (free_count = 0; free_count < saved_free_count; free_count++) {
  119. int i1 = free_row[free_count], low = 0, up = 0, last, k;
  120. int min, c, u1;
  121. for (j = 0; j < column_count; j++) {
  122. d[j] = COST(j, i1) - v[j];
  123. pred[j] = i1;
  124. col[j] = j;
  125. }
  126. j = -1;
  127. do {
  128. last = low;
  129. min = d[col[up++]];
  130. for (k = up; k < column_count; k++) {
  131. j = col[k];
  132. c = d[j];
  133. if (c <= min) {
  134. if (c < min) {
  135. up = low;
  136. min = c;
  137. }
  138. col[k] = col[up];
  139. col[up++] = j;
  140. }
  141. }
  142. for (k = low; k < up; k++)
  143. if (column2row[col[k]] == -1)
  144. goto update;
  145. /* scan a row */
  146. do {
  147. int j1 = col[low++];
  148. i = column2row[j1];
  149. u1 = COST(j1, i) - v[j1] - min;
  150. for (k = up; k < column_count; k++) {
  151. j = col[k];
  152. c = COST(j, i) - v[j] - u1;
  153. if (c < d[j]) {
  154. d[j] = c;
  155. pred[j] = i;
  156. if (c == min) {
  157. if (column2row[j] == -1)
  158. goto update;
  159. col[k] = col[up];
  160. col[up++] = j;
  161. }
  162. }
  163. }
  164. } while (low != up);
  165. } while (low == up);
  166. update:
  167. /* updating of the column pieces */
  168. for (k = 0; k < last; k++) {
  169. int j1 = col[k];
  170. v[j1] += d[j1] - min;
  171. }
  172. /* augmentation */
  173. do {
  174. if (j < 0)
  175. BUG("negative j: %d", j);
  176. i = pred[j];
  177. column2row[j] = i;
  178. SWAP(j, row2column[i]);
  179. } while (i1 != i);
  180. }
  181. free(col);
  182. free(pred);
  183. free(d);
  184. free(v);
  185. free(free_row);
  186. }