blob: 81800371b70ee8bec2314a15c1e7832dc2aad6a8 [file] [log] [blame]
/* AUTORIGHTS
Copyright (C) 2007 Princeton University
This file is part of Ferret Toolkit.
Ferret Toolkit is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
This program 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. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <values.h>
#include <cass.h>
#define M_STAR 1
#define M_PRIME 2
#define M_FLIP 4
/* http://www.netlib.org/utk/lsi/pcwLSI/text/node222.html */
void pas (int n, char **M)
{
int i, j;
printf("---\n");
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
if (M[i][j] == M_STAR) printf("0*");
else if (M[i][j] == M_PRIME) printf("0'");
else printf("--");
}
printf("\n");
}
}
float as_solve (int n, const float **_cost, char **M)
{
int i, j, k;
float result;
float min;
char *r_cov, *c_cov;
float **cost;
int cnt;
/* preprocess */
cost = (float **)type_matrix_alloc(float, n, n);
for (i = 0; i < n; i++)
{
min = _cost[i][0];
for (j = 1; j < n; j++) if (_cost[i][j] < min) min = _cost[i][j];
for (j = 0; j < n; j++) cost[i][j] = _cost[i][j] - min;
}
for (j = 0; j < n; j++)
{
min = cost[0][j];
for (i = 1; i < n; i++) if (cost[i][j] < min) min = cost[i][j];
for (i = 0; i < n; i++) cost[i][j] -= min;
}
assert(M != NULL);
/* stage 2 */
r_cov = (char *)alloca(n);
c_cov = (char *)alloca(n);
step1:
//printf("S1\n");
memset(r_cov, 0, n);
memset(c_cov, 0, n);
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
if (cost[i][j] == 0 && !r_cov[i] && !c_cov[j])
{
r_cov[i] = c_cov[j] = 1;
M[i][j] = M_STAR;
}
else M[i][j] = 0;
}
//pas(n, M);
step2:
//printf("S2\n");
cnt = 0;
memset(r_cov, 0, n);
memset(c_cov, 0, n);
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
if (M[i][j] == M_STAR)
{
c_cov[j] = 1;
cnt++;
break;
}
//pas(n, M);
if (cnt == n) goto done;
step3:
//printf("S3\n");
//pas(n, M);
for (j = 0; j < n; j++) if (!c_cov[j])
for (i = 0; i < n; i++) if (!r_cov[i])
if (cost[i][j] == 0)
{
M[i][j] = M_PRIME;
for (k = 0; k < n; k++) if (M[i][k] == M_STAR) break;
if (k >= n) goto step4;
c_cov[k] = 0;
r_cov[i] = 1;
goto step3;
}
goto step5;
step4:
//printf("S4 (%d,%d)\n", i, j);
//pas(n, M);
/* i, j from step 3 */
M[i][j] |= M_FLIP;
for (;;)
{
for (i = 0; i < n; i++) if (M[i][j] == M_STAR) break;
if (i >= n) break;
M[i][j] |= M_FLIP;
for (j = 0; j < n; j++) if (M[i][j] == M_PRIME) break;
assert(j < n);
M[i][j] |= M_FLIP;
}
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
if ((M[i][j] == (M_FLIP | M_PRIME)) || (M[i][j] == M_STAR))
{
M[i][j] = M_STAR;
}
else M[i][j] = 0;
}
goto step2;
step5:
// printf("S5\n");
min = MAXFLOAT;
for (i = 0; i < n; i++) if (!r_cov[i])
for (j = 0; j < n; j++) if (!c_cov[j])
{
if (cost[i][j] < min) min = cost[i][j];
}
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
if (r_cov[i]) cost[i][j] += min;
if (!c_cov[j]) cost[i][j] -= min;
}
goto step3;
done:
matrix_free(cost);
result = 0;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
if (M[i][j]) result += _cost[i][j];
return result;
}