/* 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;

}