#include "MemoryManager.h"
#include "spnbox.h"

static int CheckParams(matrix* D, int** URT, int* urtCount);
static int BuildIntegerProblem(MemoryManager* mem, matrix* D, int* URT, int urtCount, matrix* M, double** F, double** B, short int** IntList, short int** ctype);

nltrans_r nltrans(matrix* D, int* UnraisableTransitions, int urtCount)
{
	MemoryManager memory;
	matrix M;
	double *F, *B;
	short int *IntList, *ctype;
	int i, j, covCount, fContinue;
	int *cov;
	ipsolve_r solution;
	nltrans_r result;
	//Initialize the result to a value suitable for error return.
	result.dtr = 0;
	result.dtrCount = 0;
	
	//Check parameters
	if (! CheckParams(D, &UnraisableTransitions, &urtCount))
	{
		return result;
	}
	
	//Keep track of dynamically-allocated memory.
	memory = CreateMemoryManager(5, 1, 0, 0);
	
	/*Initialize the integer problem. The return value corresponds to the number
	of variables that should matched to cov values.*/	
	covCount = BuildIntegerProblem(&memory, D, UnraisableTransitions, urtCount, &M, &F, &B, &IntList, &ctype);
	
	/*Allocate space for the cov vector. cov(i) will be set to 1 if transition i
	belongs to a nonnegative invariant x. It should begin zeroed.*/
	cov = tcalloc(covCount, sizeof(int));
	
	/*fContinue indicates whether the loop should continue. Default to true.*/
	fContinue = 1;

	/*Main processing loop.*/
	do
	{
		//Build the last constraint row. It should be the logic inverse of cov.
		for (i = 0; i < covCount; i++)
		{			
			SetMatrixEl(&M, NumberOfRows(M) - 1, i, ! cov[i]);
		}	
		/*Solve the integer problem. Lower bounds default to zero and upper to
		infinity, which is just what we want.*/
		solution = ipsolve(&M, B, F, IntList, 0, 0, ctype);	
		
		/*Check the solution. If any constraint could not be solved for, break out 
		of the loop..*/
		fContinue = ! strcmp(solution.mhow, HOW_OK);
		
		if (fContinue)
		{
			/*If all constraints were solved, then take the logic or of each element
			of cov with a corresponding element from the solution. Ignore very small
			nonzero values in the solution.*/
			for (i = 0; i < covCount; i++)
			{
				if (solution.res[i] > HOW_MAX_ERROR || solution.res[i] < -HOW_MAX_ERROR)
				{
					cov[i] = 1;
				}
			}			
			/*Check. If all elements of cov are nonzero, we want to exit the loop.*/
			fContinue = 0;
			for (i = 0; i < covCount; i++)
			{
				if (! cov[i])
				{
					fContinue = 1;
					break;
				}
			}
		}
		DeallocateIpsolve(&solution);
	} while (fContinue);
	
	/*The return value is the logic inverse of cov.*/
	for (i = 0; i < covCount; i++)
	{
		cov[i] = ! cov[i];
	}
	
	result.dtrCount = covCount;
	result.dtr = cov;
	
	FreeMemory(&memory);
	
	return result;
}

/*******************************************************************************
BuildIntegerProblem builds a portion of the linear programming problem used to 
determine if a given transition is raisable. The portion it builds is that
portion which does not change from transition to transition. It takes incidence 
matrix, the list of predefined unraisable transitions, and pointers to the 
various parameters that will be used to defined the linear programming problem.*/
int BuildIntegerProblem(MemoryManager* mem, matrix* D, int* URT, int urtCount, matrix* M, double** F, double** B, short int** IntList, short int** ctype)
{
	int Constraints, Vars, i, j;
	/*Initialize the parameters that will be passed to ipsolve to describe an
	integer programming problem.*/
	
	/*The final constraint matrix will be D with a pair of rows beneath for each
	already unraisable transition and a minus identity tacked onto the right, with
	a single constraint row tacked on beneath it all. The last constraint row
	varies with loop iteration so it will not be calculated here.*/
	
	/*Find final sizes and allocate storage space.*/
	Constraints = NumberOfRows(*D) + (urtCount * 2) + 1;
	Vars = NumberOfColumns(*D) + Constraints - 1;
	MAllocateMatrixType(mem, 2, M, Constraints, Vars);
	
	/*Copy in D.*/
	CopyBlock(NumberOfRows(*D), NumberOfColumns(*D), D, 0, 0, M, 0, 0);
	
	/*One of the two rows inserted for each URT contains a single 1 in the 
	column corresponding to the transition. The other row contains a single -1 in
	the same position. This ensures that the constraints associated with this 
	transition can never be met.*/
	for (i = 0; i < urtCount; i++)
	{
		SetMatrixEl(M, NumberOfRows(*D) + i, URT[i], 1);
		SetMatrixEl(M, NumberOfRows(*D) + i + urtCount, URT[i], -1);
	}
	
	/*There should be a negative identity concatenated on to this on the 
	right-hand side.*/
	for (i = 0; i < Constraints - 1; i++)
	{
		SetMatrixEl(M, i, NumberOfColumns(*D) + i, -1);
	}
	
	/*The cost vector is all ones.*/
	(*F) = mcalloc(mem, Vars, sizeof(double));
	for (i = 0; i < Vars; i++)
	{
		(*F)[i] = 1.0;
	}
	/*The constraint vector is zeros except for the last one, which is 1.*/
	(*B) = mcalloc(mem, Constraints, sizeof(double));
	(*B)[Constraints - 1] = 1.0;
	
	/*The intlist should be zeros - nothing should be constrained as an integer.*/
	(*IntList) = mcalloc(mem, Vars, sizeof(short int));

	/*All constraints but the last are =. The last is >=*/
	(*ctype) = mcalloc(mem, Constraints, sizeof(short int));
	(*ctype)[Constraints - 1] = 1;

	/*The return value is the number of variables that will be used for cov.*/
	return NumberOfColumns(*D);
}

/*******************************************************************************
CheckParams checks the parameters for validity and makes sure that the URT
pointer is null if the urtCount is zero and vice-versa.*/
int CheckParams(matrix* D, int** URT, int* urtCount)
{
	int i;
	/*D is required*/
	if (! D)
	{
		merror(0, "NLTRANS: The pointer to the incidence matrix is null");
		return 0;
	}
	/*Make sure that urtCount is 0 if URT is null, and vice-versa.*/
	if (! *urtCount)
	{
		*URT = 0;
	}
	else if (! *URT)
	{
		*urtCount = 0;
	}
	/*Make sure that none of the transition numbers in the URT list are out of
	range.*/
	for (i = 0; i < *urtCount; i++)
	{
		if ((*URT)[i] >= NumberOfColumns(*D) || (*URT)[i] < 0)
		{
			printf("NLTRANS: Unraisable transition index %d is out of bounds.", (*URT)[i]);
			printf("It should satisfy 0 < index < %d.", NumberOfColumns(*D));
			merror(0, "NLTRANS: Unraisable Transition Index is out-of-bounds.");
			return 0;
		}
	}
	return 1;
}
	

