#include <stdlib.h>
#include "../pnheaders/matrix.h"
#include "matrixmath.h"
#include "extendedmatrix.h"
#include "MemoryManager.h"
#include "spnbox.h"

static int CheckParams(matrix* Dm, matrix* Dp, int** X, int* XCount, matrix *Dcm, matrix *Dcp, int update);

static int isUnique(matrix* D, int* X, int XCount);

static actn_r BuildResult(matrix* Dm, matrix* Dp, int* cullPlace, int* cullTrans);

static int* cullPlaces(matrix* Dp, int* cullTransition);

static int* updateURT(matrix* Dm, matrix* Dp, matrix* Dcm, matrix* Dcp);

static int CheckParams(matrix* Dm, matrix* Dp, int** X, int* XCount, matrix *Dcm, matrix *Dcp, int update);

actn_r actn(matrix* Dm, matrix* Dp, int* X, int XCount, matrix *Dcm, matrix *Dcp, int update, int checkUnique)
{
	actn_r result;
	nltrans_r nltransResult;
	MemoryManager memory;
	matrix D; /*The incidence matrix*/
	int *cullTransition; /*An array of flags, one for each transition, set if it
	                     should be culled (is unraisable)*/
	int *cullPlace; /*An array of flags, one for each place, set if the place
									should be culled.*/	
	
	/*Initialize the error return value*/
	memset(&result, 0, sizeof(actn_r));
	
	/*Initialize the memory manager*/
	memory = CreateMemoryManager(5, 5, 0, 0);
	
	/*Check the parameters.*/
	if (! CheckParams(Dm, Dp, &X, &XCount, Dcm, Dcp, update)) return result;
	
	/*Get the incidence matrix - we'll need it later.*/
	D = SubtractMatrix(Dp, Dm, (matrix*) 2);
	ManageMatrix(&memory, &D);
	
	if (update)
	{
		/*If we are in update mode, use the updateURT function to get the
		unraisable transitions.*/
		cullTransition = updateURT(Dm, Dp, Dcm, Dcp);
	}
	else
	{
		/*Otherwise, use nltrans to get the list of unraisable transitions.*/
		nltransResult = nltrans(&D, X, XCount);
		if (! (nltransResult.dtrCount))
		{
			/*If the call failed, fail.*/
			FreeMemory(&memory);
			merror(0, "ACTN: nltrans call failed");
			return result;
		}
		else
		{
			cullTransition = nltransResult.dtr;			
		}
	}
	ManageMemory(&memory, cullTransition);
	
	/*Build the cullPlace flags, which will be cleared only for places that
	have nonzero input arcs from transitions that are not being culled.*/
	cullPlace = cullPlaces(Dp, cullTransition);
	ManageMemory(&memory, cullPlace);
	
	/*Build the final output.*/
	result = BuildResult(Dm, Dp, cullPlace, cullTransition);
	
	/*Do an isUnique test only if the checkUnique flag is set, since the test 
	does several linear programming problems and is thus rather slow.*/
	if (checkUnique)
	{
		result.unique = isUnique(&D, X, XCount);
	}
	
	/*Clean up.*/
	FreeMemory(&memory);
	
	return result;
}

/*******************************************************************************
isUnique examines the petri net and the given set of unraisable transitions to
determine if the live subnet that will be returned is unique. It is only called
when the actn caller requests it.*/
int isUnique(matrix* D, int* X, int XCount)
{
	/*We will be doing a linear programming problem for each raisable transition. 
	The constraint matrix will be the subset of D having only the	columns 
	corresponding to transitions not in X, with a row of 1s concatenated onto the 
	bottom, and the column corresponding to whichever transition is currently 
	being tested removed (or in this case, zeroed).*/
	matrix Dx, dx;
	int i, j, Transitions = 0;
	int* UnraisableT;
	double* B;
	short int *IntList;
	ipsolve_r result;	
	
	//If no transitions are being evaluated, return true.
	if ((Transitions = NumberOfColumns(*D) - XCount) <= 0) return 1;
	/*Build a flag array of the columns whose indices are in X.*/
	UnraisableT = tcalloc(NumberOfColumns(*D), sizeof(int));
	for (i = 0; i < XCount; i++)
	{
		UnraisableT[X[i]] = 1;
	}
	/*Build Dx. It should include only columns of D whose indices are not in X.
	It should be optimized for column ops.*/
	AllocateMatrixType(2, &Dx, Transitions, NumberOfRows(*D) + 1);	
	TransposeMatrix(&Dx);
	j = 0;
	for (i = 0; i < NumberOfColumns(*D); i++)
	{		
		if (! UnraisableT[i])
		{
			/*Fill in the last row.*/
			SetMatrixEl(&Dx, NumberOfRows(*D), j, 1);
			/*Fill in the other rows*/
			CopyBlock(NumberOfRows(*D), 1, D, 0, i, &Dx, 0, j++);
		}
	}
	//We are done with the flag array. Free memory.
	free(UnraisableT);
	
	/*Allocate space to record the values of the current column while it is
	zeroed. Make it a type-2 transpose to optimize the column-swap operation.*/
	AllocateMatrixType(1, &dx, 1, NumberOfRows(*D) + 1);
	TransposeMatrix(&dx);
	
	/*Build B. It is all zeroes except for a 1 in the last element.*/
	B = tcalloc(NumberOfRows(*D) + 1, sizeof(double));
	B[NumberOfRows(*D)] = 1;
	
	/*IntList needs to be all zeroes to override the default of all 1s.*/
	IntList = tcalloc(Transitions, sizeof(short int));
	
	/*Iterate through each transition*/
	for (i = 0; i < Transitions; i++)
	{
		/*If there is a previous transition to restore, restore it.*/
		if (i) SwapColumns(&Dx, &dx, i - 1, 0);		
		/*Swap out the current transition. Because dx is a column of zeros this
		will effectively zero the current transition and save a copy in dx for later
		restoration.*/
		SwapColumns(&Dx, &dx, i, 0);
		
		//ipsolve_r ipsolve(matrix* L, double* B, double* f, short int *IntList, double *ub, double *lb, short int *ctype);
		result = ipsolve(&Dx, B, 0, IntList, 0, 0, 0);
		if (result.res) free(result.res);
		
		/*If we get a valid solution, then free memory and return false.*/
		if (! strcmp(result.mhow, HOW_OK))
		{
			free(IntList);
			free(B);			
			DeallocateMatrix(&Dx);
			DeallocateMatrix(&dx);
			return 0;
		}
	}
	//If we get through all transition without returning false, return true.
	free(IntList);
	free(B);	
	DeallocateMatrix(&Dx);
	DeallocateMatrix(&dx);
	return 1;	
}

/*******************************************************************************
BuildResult takes the list of places and transitions to cull from the total net
and builds the active subnet. It returns the result structure that will be
returned by actn.
*/
actn_r BuildResult(matrix* Dm, matrix* Dp, int* cullPlace, int* cullTrans)
{
	int Places = 0, Transitions = 0, i, j, k, l;
	actn_r result;
	
	/*All matrices are zeroed initially.*/
	memset(&result, 0, sizeof(actn_r));
	
	/*Count the number of places and transitions that are going to be left after
	culling.*/
	for (i = 0; i < NumberOfRows(*Dm); i++)
	{
		if (! cullPlace[i]) Places++;
	}
	for (i = 0; i < NumberOfColumns(*Dm); i++)
	{
		if (! cullTrans[i]) Transitions++;
	}
	
	/*Allocate space. Dma/Dmp are the same size as the original matrices but
	contain only elements in non-culled rows and columns. Dmra/Dpra contain only
	non-culled rows and columns. TA holds the indices of kept transitions.*/
	if (Places && Transitions)
	{
		AllocateMatrixType(2, &result.Dmra, Places, Transitions);
		AllocateMatrixType(2, &result.Dpra, Places, Transitions);
	}
	AllocateMatrixType(2, &result.Dma, NumberOfRows(*Dm), NumberOfColumns(*Dm));
	AllocateMatrixType(2, &result.Dpa, NumberOfRows(*Dm), NumberOfColumns(*Dm));
	result.TA = tcalloc(Transitions, sizeof(int));
	result.TACount = Transitions;
	/*Do the copy. k and l indicate the row and column numbers in the Dmra/Dpra
	matrices. Fill in TA while we're at it.*/	
	l = 0;
	for (j = 0; j < NumberOfColumns(*Dm); j++)
	{
		if (cullTrans[j]) continue;
		k = 0;
		for (i = 0; i < NumberOfRows(*Dm); i++)
		{
			if (cullPlace[i]) continue;
			SetMatrixEl(&result.Dma, i, j, GetMatrixEl(Dm, i, j));
			SetMatrixEl(&result.Dpa, i, j, GetMatrixEl(Dp, i, j));
			SetMatrixEl(&result.Dmra, k, l, GetMatrixEl(Dm, i, j));
			SetMatrixEl(&result.Dpra, k, l, GetMatrixEl(Dp, i, j));
			k++;
		}
		result.TA[l] = j;
		l++;
	}	
	return result;
}
/*******************************************************************************
cullPlaces examines the Petri net and the list of transitions to cull and
returns an array of integer flags, one for each place, indicating whether that
place should be culled.
*/
int* cullPlaces(matrix* Dp, int* cullTransition)
{
	int i, j;
	int *cull;
	/*Allocate and initialize the return value cull. It will be clear only for
	places that have nonzero output arcs to transitions that are not being
	culled.*/
	cull = tcalloc(NumberOfRows(*Dp), sizeof(int));
	for (i = 0; i < NumberOfRows(*Dp); i++)
	{
		/*Look to see if the current place has any nonzero output arcs to
		transitions not being culled.*/
		for (j = 0; j < NumberOfColumns(*Dp); j++)
		{
			if (GetMatrixEl(Dp, i, j) && (! cullTransition[j]))
			{				
				break;
			}
		}
		/*If there were no nonzero output arcs set the place flag.*/
		if (j == NumberOfColumns(*Dp))
		{
			cull[i] = 1;
		}
	}
	return cull;
}

/*******************************************************************************
updateURT examines the Petri net and the predefined active subnet if given and
returns an array of transitions that are unraisable given the active subnet.
This is used to determine transitions to cull if in update mode.
*/
int* updateURT(matrix* Dm, matrix* Dp, matrix* Dcm, matrix* Dcp)
{
	/*This function builds a list of unraisable transitions in Dm/Dp by examining
	Dcm/Dcp.*/	
	int Places, Transitions, SubPlaces, SubTransitions, i, j, keepGoing;
	/*These are each flag arrays set to indicate live places or transitions at
	various times throughout the algorithm.*/
	int *URT, *Trans = 0, *Place = 0;	
	
	Places = NumberOfRows(*Dm);
	Transitions = NumberOfColumns(*Dm);
	SubPlaces = NumberOfRows(*Dcm);
	SubTransitions = NumberOfColumns(*Dcm);
	
	/*URT is the return value. It has a flag for each transition.*/
	URT = tcalloc(Transitions, sizeof(int));	
	/*Trans and Place both have flags only for each of the possible new
	transitions or places that can be added.*/
	if (Places - SubPlaces)
	{
		Place = tcalloc(Places - SubPlaces, sizeof(int));	
	}
	if (Transitions - SubTransitions)
	{
		Trans = tcalloc(Transitions - SubTransitions, sizeof(int));
	}
		
	/*URT is a flag vector with an element for each transition. It is initially
	cleared only for non-null columns of Dcm/Dcp, where column numbers in Dcm/Dcp 
	map	to the same column numbers in Dm/Dp.*/
	for (i = 0; i < SubTransitions; i++)
	{
		/*See if the current column of Dcm/Dcp is non-null.*/
		for (j = 0; j < SubPlaces; j++)
		{
			if (GetMatrixEl(Dcm, j, i) || GetMatrixEl(Dcp, j, i))
			{
				break;
			}
		}
		/*If it was null, set the URT flag.*/
		if (j == SubPlaces)
		{
			URT[i] = 1;
		}
	}
	/*Fill the additional elements of URT, the ones not taken care of above.*/
	for (i = SubTransitions; i < Transitions; i++)
	{
		URT[i] = 1;
	}
	
	/*Now, we loop. For the first iteration we have to do something special. Use
	keepGoing = -1 to indicate first loop iteration.*/
	keepGoing = -1;
	do
	{
		/*The place flag should be set for any place (in Dm/Dp beyond those mapped 
		to by Dcm/Dcp) such that the place has at least one nonzero output arc to 
		a transition that is flagged in the current transition list. It should
		be cleared for every other place. On the first loop iteration, rather
		than using the current transition list (which maps to the columns in the
		right of Dm/Dp that have no corresponding columns in Dcm/Dcp), we use logic 
		inverse of URT,	which maps to all columns in Dm/Dp. The keepGoing flag will 
		be set to -1 on	the first iteration.*/
		memset(Place, 0, sizeof(int) * (Places - SubPlaces));
		if (keepGoing > 0)
		{			
			for (i = SubPlaces; i < Places; i++)
			{
				/*Check to see if the current place has at least one nonzero output arc
				to a transition flagged in the current transition list.*/
				for (j = SubTransitions; j < Transitions; j++)
				{
					if (Trans[j - SubTransitions] && GetMatrixEl(Dm, i, j))
					{
						/*If it does, set the flag and break out of the test loop.*/
						Place[i - SubPlaces] = 1;
						break;
					}
				}
			}
		}
		else
		{
			for (i = SubPlaces; i < Places; i++)
			{
				/*Check to see if the current place has at least one nonzero output arc
				to a transition not flagged in URT.*/
				for (j = 0; j < SubTransitions; j++)
				{
					if ((! URT[j]) && GetMatrixEl(Dm, i, j))
					{
						/*If it does, set the flag and break out of the test loop.*/
						Place[i - SubPlaces] = 1;
						break;
					}
				}
			}
		}	
		/*Similarly, the transition flag should be set for any transition in Dm/Dp
		beyond those mapped to by Dcm/Dcp such that the transition has at least one
		nonzero input arc to a place that is flagged in the current place list.
		Additionally, we clear every flag of URT that corresponds to a
		transition whose flag is set in Trans. Finally, the overall loop will 
		continue only until Trans is entirely empty. Set the keepGoing flag
		if any element of Trans is set.*/
		keepGoing = 0;
		memset(Trans, 0, sizeof(int) * (Transitions - SubTransitions));
		for (j = SubTransitions; j < Transitions; j++)
		{
			for (i = SubPlaces; i < Places; i++)
			{
				if (Place[i - SubPlaces] && GetMatrixEl(Dp, i, j))
				{
					Trans[j - SubTransitions] = 1;
					/*Clear the appropriate URT flag.*/
					URT[j] = 0;
					/*Set the keepGoing flag.*/
					keepGoing = 1;
					break;
				}
			}
		}
	} while (keepGoing);
	/*Free memory*/
	if (Place) free(Place);
	if (Trans) free(Trans);
	/*Return*/
	return URT;
}

/*******************************************************************************
CheckParams checks the parameters to make sure required parameters are present,
checks some basic validity, and makes sure that all array pointers are null if 
their counts are zero and vice-versa.*/
int CheckParams(matrix* Dm, matrix* Dp, int** X, int* XCount, matrix *Dcm, matrix *Dcp, int update)
{
	int Rows0, Cols0, Rows1, Cols1;
	/*Dm and Dp are required and must be the same size.*/
	if (! (Dm && Dp))
	{
		merror(0, "ACTN: Dm or Dp was not given");
		return 0;
	}
	Rows0 = NumberOfRows(*Dm);
	Cols0 = NumberOfColumns(*Dm);
	if (Rows0 != NumberOfRows(*Dp) || Cols0 != NumberOfColumns(*Dp))
	{
		merror(0, "ACTN: The dimensions of Dm and Dp do not match");
		return 0;
	}
	
	/*X is not required. Make sure that XCount is zero if X is null and vice-
	versa.*/
	if (! *X)
	{
		*XCount = 0;
	}
	else if (! *XCount)
	{
		*X = 0;
	}
	
	/*If update is nonzero, Dcm and Dcp are required. They cannot be larger than
	Dm/Dp in either dimension and they must be the same size.*/
	if (update)
	{
		if (! (Dcm && Dcp))
		{
			merror(0, "ACTN: Dcm or Dcp was not given while update was nonzero");
			return 0;
		}
		Rows1 = NumberOfRows(*Dcm);
		Cols1 = NumberOfColumns(*Dcm);
		if (NumberOfRows(*Dcp) != Rows1 || NumberOfColumns(*Dcp) != Cols1)
		{
			merror(0, "ACTN: Dcm and Dcp are not the same size");
			return 0;
		}
		if (Rows1 > Rows0 || Cols1 > Cols0)
		{
			merror(0, "ACTN: Dcm and Dcp are larger than Dm and Dp");
			return 0;
		}
	}
	return 1;
}
