linear-assignment.c

   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
  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
  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 }