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

/*ModifyD makes modifications to the Petri net I/O 
matrices for the given transition, based on the
current number of transitions and places and the number
of transitions that are being added.*/
static void ModifyD(msplit_r* Ret, int Places, int Transitions, int NewTransitions, int Transition);

/*ModifyTD fills out the elements of array of dead-transition arrays for the 
current transition.*/
static void ModifyTD(msplit_r* Ret, int Transition, int NewTransitions, int TDCount);

/*ModifyML0 modifies ML0 for the current transition.*/
static void ModifyML0(msplit_r* Ret, int* ipl, int ipCount, int* dpl, int MColCount, int L0ColCount, int Transition, int NewTransitions);

/*AllocateTransition space allocates space for and initializes some of the 
data that will be returned.*/
static void AllocateTransitionSpace(msplit_r* Ret, matrix* Dm, matrix* Dp, int* mask, matrix* M, matrix* L0, matrix* L, int** TD, int* TDCount, int* ipl, int ipCount);

/*CountNewTransitions returns an integer corresponding to the number of places
and transitions that will need to be added to split the given transition.*/
static int CountNewTransitions(matrix* Dm, int Transition);

/*CheckParams checks the validity of the parameters and initializes some of them
to default values.*/
static int CheckParams(MemoryManager* mem, matrix* Dm, matrix* Dp, int** mask, matrix* M, matrix* L0, matrix* L, int** ipl, int* ipCount, int** dpl, int* dpCount, int*** TD, int** TDCount);

msplit_r msplit(matrix* Dm, matrix* Dp, int* mask, matrix* M, matrix* L0, matrix* L, int* ipl, int ipCount, int* dpl, int dpCount, int** TD, int* TDCount)
{
	msplit_r Ret;
	MemoryManager memory;
	int Places, Transitions, NewTransitions, j, MCols, L0Cols, Count;
	int OldTransitions;

	/*Initialize the returned structure to an error value.*/
	memset(&Ret, 0, sizeof(msplit_r));
	
	/*Initialize the memory manager.*/
	memory = CreateMemoryManager(1, 3, 0, 0);

	/*Check/set to default the parameters and return if parameters were bad*/
	if (! CheckParams(&memory, Dm, Dp, &mask, M, L0, L, &ipl, &ipCount, &dpl, &dpCount, &TD, &TDCount))
	{
		return Ret;
		FreeMemory(&memory);
		merror(0, "MSPLIT: The parameters were in error.");
	}
	
	Places = NumberOfRows(*Dm);
	Transitions = NumberOfColumns(*Dm);
	OldTransitions = NumberOfColumns(*Dm);
	/*Allocate storage space and initialize the return structure.*/
	AllocateTransitionSpace(&Ret, Dm, Dp, mask, M, L0, L, TD, TDCount, ipl, ipCount);
	
	/*Get some initial matrix and vector sizes*/
	if (M)
	{
		MCols = NumberOfColumns(*M);
	}
	else
	{
		MCols = -1;
	}
	if (L0)
	{
		L0Cols = NumberOfColumns(*L0);
	}
	else
	{
		L0Cols = -1;
	}	
	
	/*Main processing loop.*/
	for (j = 0; j < OldTransitions; j++)
	{
		/*If we are not supposed to process this transition skip to the next one.*/
		if (! mask[j]) continue;
		
		/*Make sure we know how many transitions are being added for this
		transition.*/
		NewTransitions = CountNewTransitions(Dm, j);
		
		/*Modify Mf and/or L0f if necessary.*/
		if (M || L0)
		{
			ModifyML0(&Ret, ipl, ipCount, dpl, MCols, L0Cols, j, NewTransitions);
		}
		/*Modify the counts if necessary*/
		if (M)
		{			
			MCols += NewTransitions;
		}
		if (L0)
		{
			L0Cols += NewTransitions;
		}
		if (TD)
		{
			ModifyTD(&Ret, j, NewTransitions, Transitions);			
		}
		/*Modify the i/o matrices*/
		ModifyD(&Ret, Places, Transitions, NewTransitions, j);
		/*Update place and transition counts*/
		Transitions += NewTransitions;
		Places += NewTransitions;
	}
	/*Compute the incidence matrix.*/
	Ret.Df = SubtractMatrix(&Ret.Dfp, &Ret.Dfm, (matrix*) 2);
	/*Free memory and return*/
	FreeMemory(&memory);
	return Ret;
}

void ModifyD(msplit_r* Ret, int Places, int Transitions, int NewTransitions, int Transition)
{
	int i, j, k;
	/*In the original Matlab implementation, a copy of Dm was made
	beforehand and used as the source for all GetMatrixEls. However,
	the use of the copy was not strictly necessary if some rearranging of
	of operations was done.*/
	for (j = 0; j < NewTransitions; j++)
	{
		/*For each output arc of weight n feeding the transition being split,
		there should be n transitions with output arcs of weight 1 from the 
		original arc's place.*/
		for (i = 0; i < Places; i++)
		{
			k = GetMatrixEl(&Ret->Dfm, i, Transition);			
			SetMatrixEl(&Ret->Dfm, i, Transitions + j, (k > j + 1) ? 1 : 0);
		}
	}
	for (j = 0; j < NewTransitions; j++)
	{
		/*Create an input arc from the new transition to the new place.*/
		SetMatrixEl(&Ret->Dfp, Places + j, Transitions + j, 1);
		if (j == 0)
		{
			/*Each output arc of weight greater than 1 feeding the transition being
			split should be changed to have a weight of 1.*/
			for (i = 0; i < Places; i++)
			{
				k = GetMatrixEl(&Ret->Dfm, i, Transition);
				SetMatrixEl(&Ret->Dfm, i, Transition, (k > 0) ? 1 : 0);
			}
			/*The newly-added places (modified by this line and the next non-trivial
			line) should each have an output arc of weight 1 feeding (for the first
			one) the transition being split, the last original transition, or or
			all but the last of the newly-added transitions*/
			SetMatrixEl(&Ret->Dfm, Places, Transition, 1);
		}
		else
		{
			SetMatrixEl(&Ret->Dfm, Places + j, Transitions + j - 1, 1);
		}
	}
}
			
void ModifyTD(msplit_r* Ret, int Transition, int NewTransitions, int TDCount)
{
	int j, i;
	for (j = 0; j < NewTransitions; j++)
	{
		Ret->TDf[TDCount + j] = tcalloc(Ret->TDfCount[Transition] + j + 1, sizeof(int));
		Ret->TDfCount[TDCount + j] = Ret->TDfCount[Transition] + j + 1;
		for (i = 0; i < Ret->TDfCount[Transition]; i++)
		{
			Ret->TDf[TDCount + j][i] = Ret->TDf[Transition][i];
		}
		for (i = 0; i <= j; i++)
		{
			Ret->TDf[TDCount + j][i + Ret->TDfCount[Transition]] = i + TDCount;
		}
	}
}

void ModifyML0(msplit_r* Ret, int* ipl, int ipCount, int* dpl, int MColCount, int L0ColCount, int Transition, int NewTransitions)
{
	int i, j, k, dm, m, l0;
	/*We iterate through all the transitions currently being added.*/		
	for (j = 0; j < NewTransitions; j++)
	{
		/*The column of M corresponding to the transition being added is filled 
		according to this expression:
		M(:,f+j) = M*((Om(ipl,i) > j).*(Om(ipl,i) - j)) - (Om(dpl,i) > j).*
		(Om(dpl,i) - j)
		The column of L0 corresponding to the transition being added is filed
		according to this expression:
		L0(:,r+j) = L0*((Om(ipl,i) > j).*(Om(ipl,i) - j));
		*/		
		if (MColCount > 0)
		{
			for (i = 0; i < NumberOfRows(Ret->Mf); i++)
			{
				m = 0;			
				/*Get the products M*Dm(ipl,i) and L0*Dm(ipl,i) for the element of 
				Dm > j. Add one to j throughout, since in the original matlab
				implementation uses 1-starting indexes, and thus j's that begin at 1, 
				and	although in C we are using zero-indexes and thus j's that begin at
				zero, the arc weights they are being compared to still begin at 1 as in
				Matlab.*/
				for (k = 0; k < ipCount; k++)
				{
					dm = GetMatrixEl(&Ret->Dfm, ipl[k], Transition);
					if (dm > j + 1)
					{
						m += GetMatrixEl(&Ret->Mf, i, ipl[k]) * (dm - j - 1);					
					}
				}
				/*Now subtract from that (Om(dpl,i) > j).*(Om(dpl,i) - j). Add 1 to j 
				throughout as above.*/
				dm = GetMatrixEl(&Ret->Dfm, dpl[i], Transition);
				if (dm > j + 1)
				{
					m -= (dm - j - 1);
				}
				/*Set the appropriate element of Mf.*/
				SetMatrixEl(&Ret->Mf, i, MColCount + j, m);
			}
		}
		if (L0ColCount > 0)
		{
			for (i = 0; i < NumberOfRows(Ret->L0f); i++)
			{
				l0 = 0;
				/*Get the products L0*Dm(ipl,i) and L0*Dm(ipl,i) for the element of 
				Dm > j*/
				for (k = 0; k < ipCount; k++)
				{
					dm = GetMatrixEl(&Ret->Dfm, ipl[k], Transition);
					if (dm > j + 1)
					{
						l0 += GetMatrixEl(&Ret->L0f, i, ipl[k]) * (dm - j - 1);
					}
				}
				/*Set the appropriate element of L0f*/
				SetMatrixEl(&Ret->L0f, i, L0ColCount + j, l0);
			}
		}	
	}	
}

void AllocateTransitionSpace(msplit_r* Ret, matrix* Dm, matrix* Dp, int* mask, matrix* M, matrix* L0, matrix* L, int** TD, int* TDCount, int* ipl, int ipCount)
{
	int ToAdd = 0, Places, Trans, i, j;
	
	
	Places = NumberOfRows(*Dm);
	Trans = NumberOfColumns(*Dm);
	
	/*Count how many transitions will be added overall.*/	
	for (j = 0; j < Trans; j++)
	{		
		ToAdd += CountNewTransitions(Dm, j);
	}
	
	/*For each new transition we add a row and a column to the incidence
	matrix. The matrix contents default to original.*/	
	AllocateMatrixType(2, &Ret->Dfm, Places + ToAdd, Trans + ToAdd);
	CopyBlock(Places, Trans, Dm, 0, 0, &Ret->Dfm, 0, 0);
	
	AllocateMatrixType(2, &Ret->Dfp, Places + ToAdd, Trans + ToAdd);
	CopyBlock(Places, Trans, Dp, 0, 0, &Ret->Dfp, 0, 0);
	
	/*For each new transition we add a column to M and L, and L0.*/
	if (M)
	{
		AllocateMatrixType(2, &Ret->Mf, NumberOfRows(*M), NumberOfColumns(*M) + ToAdd);
		CopyBlock(NumberOfRows(*M), NumberOfColumns(*M), M, 0, 0, &Ret->Mf, 0, 0);
	}
	
	if (L)
	{
		AllocateMatrixType(2, &Ret->Lf, NumberOfRows(*L), NumberOfColumns(*L) + ToAdd);
		CopyBlock(NumberOfRows(*L), NumberOfColumns(*L), L, 0, 0, &Ret->Lf, 0, 0);
	}

	if (L0)
	{
		AllocateMatrixType(2, &Ret->L0f, NumberOfRows(*L0), NumberOfColumns(*L0) + ToAdd);
		CopyBlock(NumberOfRows(*L0), NumberOfColumns(*L0), L0, 0, 0, &Ret->L0f, 0, 0);
	}
	
	/*Each new transition gets its own element in the TD array. Prexisting
	elements in the TD array will remain unchanged.*/	
	if (TD)
	{
		Ret->TDf = tcalloc(Trans + ToAdd, sizeof(int*));
		Ret->TDfCount = tcalloc(Trans + ToAdd, sizeof(int));
		for (i = 0; i < Trans; i++)
		{
			Ret->TDfCount[i] = TDCount[i];
			Ret->TDf[i] = tcalloc(TDCount[i], sizeof(int));
			for (j = 0; j < TDCount[i]; j++)
			{
				Ret->TDf[i][j] = TD[i][j];
			}
		}
	}
	/*The new list of independent places will include all the original places 
	and all the places that will be added.*/
	Ret->iplf = tcalloc(ipCount + ToAdd, sizeof(int));
	Ret->ipCountf = ipCount + ToAdd;
	for (i = 0; i < ipCount; i++)
	{
		Ret->iplf[i] = ipl[i];
	}
	for (i = 0; i < ToAdd; i++)
	{
		Ret->iplf[i + ipCount] = Places + i;
	}	
}

/*Count the number of new transitions that need to be added for the given 
transition number. If the transition number is negative, count the number
of transitions that will be added to reach that transition.*/
int CountNewTransitions(matrix* Dm, int Transition)
{
	int i, total, c;
	total = 0;
	/*For each transition with which we work, the number of new transitions
	added is equal to the largest arc weight that passes from to that
	transition. Find the max and add it to the total.*/
	for (i = 0; i < NumberOfRows(*Dm); i++)
	{
		c = GetMatrixEl(Dm, i, Transition);
		if (c > total) total = c;
	}
	if (--total < 0) total = 0;
	return total;
}
			

int CheckParams(MemoryManager* mem, matrix* Dm, matrix* Dp, int** mask, matrix* M, matrix* L0, matrix* L, int** ipl, int* ipCount, int** dpl, int* dpCount, int*** TD, int** TDCount)
{
	/*Parameter requirements:
	Dm and Dp are required and must be of identical size.
	mask is not required but has a default which will need to be computed.
	M is not required but if present has required dimensions. These are not
		checked by the original matlab function but probably ought to be checked.
	L0 is not required but if present has a required number of columns.
	L is not required but if present has a required number of columns.
	ipl and dpl are not required but have defaults that will need to be
	computed.
	TD is not required; TDCount should be nulled if TD is null and vice-versa.*/
	int Places, Transitions, i, j;
	/*Make sure the required parameters have been passed*/
	if (! (Dm && Dp))
	{
		merror(0, "MSPLIT: A required parameter was not given");
		return 0;
	}
	/*Get a record of place/transition numbers and make sure that i/o matrix
	dimensions agree*/
	Places = NumberOfRows(*Dm);
	Transitions = NumberOfColumns(*Dm);
	
	if (Places != NumberOfRows(*Dp) || Transitions != NumberOfColumns(*Dp))
	{
		merror(0, "MSPLIT: Dm and Dp do not agree in size");
		return 0;
	}
	/*Setup the default mask if necessary.*/
	if (! (*mask))
	{
		(*mask) = mcalloc(mem, Transitions, sizeof(int));
		for (j = 0; j < Transitions; j++)
		{
			for (i = 0; i < Places; i++)
			{
				if (GetMatrixEl(Dm, i, j) > 1)
				{
					(*mask)[j] = 1;
					break;
				}
			}
		}
	}
	/*Make sure place counts and array pointers agree if one or the other
	is null.*/
	if (! *ipCount) 
	{
		(*ipl) = 0;
	}
	else if (! *ipl)
	{
		(*ipCount) = 0;
	}
	if (! *dpCount) 
	{
		(*dpl) = 0;
	}
	else if (! *dpl)
	{
		(*dpCount) = 0;
	}
	/*Setup the default independent/dependent place lists if necessary.*/
	if (! (*ipl))
	{
		(*ipl) = mcalloc(mem, Places, sizeof(int));
		(*ipCount) = Places;
		(*dpl) = 0;
		(*dpCount) = 0;
		for (i = 1; i < Places; i++)
		{
			(*ipl)[i] = i;
		}
	}
	/*Make sure the independent/dependent place list sizes agree.*/
	if ((*ipCount) + (*dpCount) != Places)
	{
		merror(0, "The lengths of lists of dependent and independent places do not agree with the Petri net dimensions");
		return 0;
	}
	
	/*If M has been given, make sure it has the right dimensions.*/
	if (M)
	{
		if (NumberOfRows(*M) != (*dpCount) || NumberOfColumns(*M) != (*ipCount))
		{
			merror(0, "The dimensions of the matrix M and the number of dependent and independent places do not agree");
			return 0;
		}
	}
	
	/*If L0 has been given make sure it is of the proper size. It must
	have the same number of columns as there are independent places.*/
	if (L0)
	{
		if (NumberOfColumns(*L0) != *ipCount)
		{
			merror(0, "The dimensions of L0 and the number of independent places do not agree");
			return 0;
		}
	}
	
	/*If L has been given make sure it is of the proper size. It must have the 
	same number of columns as there are independent places.*/
	if (L)
	{
		if (NumberOfColumns(*L0) != *ipCount)
		{
			merror(0, "The dimensions of L and the number of independent places do not agree");
			return 0;
		}
	}
	/*Make sure the nullity of TD and TDCount agree*/
	if (! TD) TDCount = 0;
	if (! TDCount) TD = 0;
	return 1;
}
