#include "../pnheaders/general.h"
#include "../pnheaders/matrix.h"
#include "MemoryManager.h"
#include "matrixmath.h"
#include "extendedmatrix.h"

typedef struct invarData
{
	matrix B, C;
	int *Ip, *Im, *I;
	MemoryGrower gmem;
} invarData;

static int InitMatrices(matrix* A, invarData* d);
static void AddGcdRows(invarData* d, int j);
static void FindNonMinimalRows(invarData* d);
static void EliminateNullSpace(invarData* d, matrix* A);

matrix invar(matrix* A)
{
	invarData d;
	if (! InitMatrices(A, &d))
	{
		return d.B;
	}
	
	/*We are going to be using several flag lists that will need to change
	size with A and C. Use a memory grower (a fast memory-area growing
	routine set in MemoryManager.c)*/
	d.gmem = CreateMemoryGrower(3, NumberOfRows(d.B) * sizeof(int));
	int i, j, flag;
	d.Ip = growCalloc(&d.gmem, 0, NumberOfRows(d.B), sizeof(int));
	d.Im = growCalloc(&d.gmem, 0, NumberOfRows(d.B), sizeof(int));
	d.I  = growCalloc(&d.gmem, 0, NumberOfRows(d.B), sizeof(int));
	
	/*Iterate through columns of B.*/
	for (j = 0; j < NumberOfColumns(d.B); j++)
	{
		/*Flag the places with input and output arcs from the current
		transition.*/
		flag = 0;
		memset(d.Im, 0, sizeof(int) * NumberOfRows(d.B));
		memset(d.Ip, 0, sizeof(int) * NumberOfRows(d.B));
		for (i = 0; i < NumberOfRows(d.B); i++)
		{
			if (GetMatrixEl(&d.B, i, j) < 0)
			{
				d.Im[i] = 1;
				flag |= 1;
			}
			if (GetMatrixEl(&d.B, i, j) > 0) 
			{
				d.Ip[i] = 1;
				flag |= 2;
			}
		}
		/*Do the gcd row adds only if there is at least one flagged place each in
		Ip and Im..*/
		if (flag == 3)
		{
			AddGcdRows(&d, j);
			/*Eliminate all rows not flagged in I.*/
			for (i = NumberOfRows(d.B) - 1; i >= 0; i--)
			{
				if (! d.I[i])
				{
					RemoveRows(&d.B, i, 1, -1);
					RemoveRows(&d.C, i, 1, -1);
				}
			}
			if (! NumberOfRows(d.B)) break;
			
			/*Flag the rows of C whose support is not minimal with respect to
			other rows of C.*/
			FindNonMinimalRows(&d);
			/*Remove those rows.*/
			for (i = NumberOfRows(d.C) - 1; i >= 0; i--)
			{
				if (d.I[i])
				{
					RemoveRows(&d.B, i, 1, -1);
					RemoveRows(&d.C, i, 1, -1);
				}
			}
			/*Make sure I, Ip, and Im are as large as they need to be for the next
			loop iteration.*/
			d.I = growCalloc(&d.gmem, d.I, sizeof(int), NumberOfRows(d.C));
			d.Ip = growCalloc(&d.gmem, d.Ip, sizeof(int), NumberOfRows(d.C));
			d.Im = growCalloc(&d.gmem, d.Im, sizeof(int), NumberOfRows(d.C));
			if (! NumberOfRows(d.B)) break;
		}		
	}
	/*Eliminate rows of C which are not part of null space.*/
	TransposeMatrix(A);
	EliminateNullSpace(&d, A);
	TransposeMatrix(A);
	TransposeMatrix(&d.C);
	/*Clean up C by reducing all the arc weights.*/
	for (j = 0; j < NumberOfColumns(d.C); j++)
	{
		int div;
		if ((div = gcdv(&d.C, -1, j)) > 1)
		{
			for (i = 0; i < NumberOfRows(d.C); i++)
			{
				SetMatrixEl(&d.C, i, j, GetMatrixEl(&d.C, i, j) / div);
			}
		}
	}
	/*Free memory.*/
	DeallocateMatrix(&d.B);
	FreeMemoryGrower(&d.gmem);
	
	return d.C;
}

/*******************************************************************************
EliminateNonNullSpace eliminates rows which are not part of null space.*/
void EliminateNullSpace(invarData* d, matrix* A)
{
	/*Get the matrix product of C and A.*/
	matrix X = MultiplyMatrix(&d->C, A, (matrix*) 2);
	int i, j;
	/*We'll be using I to flag rows that should be kept. Clear it now.*/
	memset(d->I, 0, sizeof(int) * NumberOfRows(d->C));
	/*Iterate through rows of C.*/
	for (i = 0; i < NumberOfRows(d->C); i++)
	{
		/*Find out if the corresponding row in X has any non-null entries.*/
		for (j = 0; j < NumberOfColumns(X); j++)
		{
			if (GetMatrixEl(&X, i, j)) break;
		}
		/*If there were no non-null entries, the row should be flagged for
		keeping. Otherwise leave the flag blank to request removal.*/
		if (j == NumberOfColumns(X)) d->I[i] = 1;
	}
	/*Remove any rows of C not flagged in I.*/
	for (i = NumberOfRows(d->C) - 1; i >= 0; i--)
	{
		if (! d->I[i]) RemoveRows(&d->C, i, 1, -1);
	}
	DeallocateMatrix(&X);
}

/*******************************************************************************
FindNonMinimimalRows finds rows of C whose support is not minimal with respect
to other rows of C and flags them for removal in I.*/
void FindNonMinimalRows(invarData* d)
{
	int i, j, k, less, more;
	/*Make sure I is big enough to hold all the new flags.*/
	d->I = growCalloc(&d->gmem, d->I, sizeof(int), NumberOfRows(d->C));
	/*Clear I.*/
	memset(d->I, 0, sizeof(int) * NumberOfRows(d->C));
	for (i = 0; i < NumberOfRows(d->C); i++)
	{
		if (d->I[i]) continue;
		for (j = i + 1; j < NumberOfRows(d->C); j++)
		{
			if (d->I[j]) continue;
			/*Check to see if logic(C(i,:)) >= logic(C(j,:)) or if logic(C(i,:)) <= 
			logic(C(j,:))*/
			less = 1;
			more = 1;
			for (k = 0; k < NumberOfColumns(d->C) && (less || more); k++)
			{
				if (more && GetMatrixEl(&d->C, j, k) && ! GetMatrixEl(&d->C, i, k))
				{
					more = 0;
				}
				if (less && GetMatrixEl(&d->C, i, k) && ! GetMatrixEl(&d->C, j, k))
				{
					less = 0;
				}
			}
			if (more)
			{
				d->I[i] = 1;
			}
			else if (less)
			{
				d->I[j] = 1;
			}
		}
	}
}

/*******************************************************************************
AddGcdRows adds some rows to the B/C pair based on the gcd of particular
elements.*/
void AddGcdRows(invarData* d, int j)
{
	int i, k, l, originalrows;
	/*Make sure that*/
	originalrows = NumberOfRows(d->B);
	memset(d->I, 0, originalrows * sizeof(int));
	for (i = 0; i < originalrows; i++)
	{
		/*Flag in I any row of B null in the current column.*/
		if (! GetMatrixEl(&d->B, i, j)) d->I[i] = 1;
		
		/*Skip from here any places without input arcs.*/
		if (! d->Ip[i]) continue;
		/*Iterate through all the places again. We will examine input places in the
		outer loop and output places in the inner loop.*/
		for (k = 0; k < originalrows; k++)
		{
			if (! d->Im[k]) continue;
			int z = gcdv(0, GetMatrixEl(&d->B, i, j), GetMatrixEl(&d->B, k, j));
			int x = GetMatrixEl(&d->B, k, j) / z;
			int y = GetMatrixEl(&d->B, i, j) / z;
			/*Add a new row to BC.*/
			InsertNullRows(&d->B, -1, 1, -1);
			InsertNullRows(&d->C, -1, 1, -1);
			/*The new row should be found algorithmically.*/
			for (l = 0; l < NumberOfColumns(d->B); l++)
			{
				int c = (-x * GetMatrixEl(&d->B, i, l)) + (y * GetMatrixEl(&d->B, k, l));
				SetMatrixEl(&d->B, NumberOfRows(d->B) - 1, l, c);
			}
			for (l = 0; l < NumberOfColumns(d->C); l++)
			{
				int c = (-x * GetMatrixEl(&d->C, i, l)) + (y * GetMatrixEl(&d->C, k, l));
				SetMatrixEl(&d->C, NumberOfRows(d->C) - 1, l, c);
			}
			/*Flag the new row in I.*/
			d->I = growCalloc(&d->gmem, d->I, sizeof(int), NumberOfRows(d->C));
			d->I[NumberOfRows(d->C) - 1] = 1;
		}
	}
}

/*******************************************************************************
InitMatrices initializes the matrices B and C.*/
int InitMatrices(matrix* A, invarData* d)
{
	memset(d, 0, sizeof(invarData));
	if (! A)
	{
		merror(0, "INVAR: Null matrix pointer");
		return 0;
	}
	/*B should be a transpose copy of A optimized for row ops.*/
	AllocateMatrixType(2, &d->B, NumberOfColumns(*A), NumberOfRows(*A));
	TransposeMatrix(&d->B);
	CopyMatrix(A, &d->B);
	TransposeMatrix(&d->B);

	/*C should be an identity, number of rows of A, optimized for row ops.*/
	AllocateMatrixType(2, &d->C, NumberOfRows(d->B), NumberOfRows(d->B));
	int i;
	for (i = 0; i < NumberOfRows(d->B); i++)
	{
		SetMatrixEl(&d->C, i, i, 1);
	}
	return 1;
}
