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

static int CheckParams(matrix* Dm, matrix* Dp, int **IncludeT, int *IncludeTCount, int ** ExcludeT, int *ExcludeTCount, matrix *Dcm, matrix *Dcp, MemoryManager* mem);

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

static int* GetKeepT0(matrix *D, int *ExcludeFlag, int *IncludeFlag, int ExcludedTCount);

static void minactn(int *keepTransition, matrix *D, int *IncludeFlag);

static int* maxactn(matrix *D, int *IncludeFlag, int *ExcludeFlag);

static int* keepPlaces(matrix* Dp, int* keepTransition);

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

static int isUnique(matrix* D, int* keepTransition, int *IncludeFlag, int *ExcludeFlag);

tactn_r tactn(matrix* Dm, matrix* Dp, int* IncludeT, int IncludeTCount, int* ExcludeT, int ExcludeTCount, matrix *Dcm, matrix *Dcp, int checkUnique)
{
	tactn_r result;
	MemoryManager memory;
	matrix D; /*The incidence matrix*/
	int *keepTransition; /*An array of flags, one for each transition, set if it
	                     should be culled (is unraisable)*/
	int *keepPlace; /*An array of flags, one for each place, set if the place
									should be culled.*/	
	int *IncludeFlag, *ExcludeFlag;
	int i;	
	
	/*Initialize the error return value*/
	memset(&result, 0, sizeof(tactn_r));
	
	/*Initialize the memory manager*/
	memory = CreateMemoryManager(10, 5, 0, 0);
	
	/*Check the parameters.*/
	if (! CheckParams(Dm, Dp, &IncludeT, &IncludeTCount, &ExcludeT, &ExcludeTCount, Dcm, Dcp, &memory)) return result;
	
	/*Get the incidence matrix - we'll need it later.*/
	D = SubtractMatrix(Dp, Dm, (matrix*) 2);
	ManageMatrix(&memory, &D);
	/*Get flag arrays corresponding to the IncludeT and ExcludeT index arrays.*/
	IncludeFlag = mcalloc(&memory, NumberOfColumns(D), sizeof(int));
	ExcludeFlag = mcalloc(&memory, NumberOfColumns(D), sizeof(int));
	for (i = 0; i < IncludeTCount; i++)
	{
		IncludeFlag[IncludeT[i]] = 1;
	}
	for (i = 0; i < ExcludeTCount; i++)
	{
		ExcludeFlag[ExcludeT[i]] = 1;
	}
	
	if (Dcm && Dcp)
	{
		/*If we are in update mode, use the updateActn function to get the
		keepable transitions.*/
		keepTransition = updateActn(Dm, Dp, Dcm, Dcp);
	}
	else
	{
		/*Otherwise, determine if we need to use minactn or maxactn to get the
		transitions in the active subnet. GetKeepT0 returns the initial keepable-
		transition flag vector needed by minactn if it should be run and a null
		pointer if maxactn should be run.*/
		if (keepTransition = GetKeepT0(&D, ExcludeFlag, IncludeFlag, ExcludeTCount))
		{
			minactn(keepTransition, &D, IncludeFlag);
		}
		else
		{
			keepTransition = maxactn(&D, IncludeFlag, ExcludeFlag);
		}		
	}
	ManageMemory(&memory, keepTransition);
	
	/*Build the keepPlace flags, which will be set only for places that
	have nonzero input arcs from transitions that are being kept.*/
	keepPlace = keepPlaces(Dp, keepTransition);
	ManageMemory(&memory, keepPlace);
	
	/*Build the final output.*/
	result = BuildResult(Dm, Dp, keepPlace, keepTransition);
	
	/*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, keepTransition, IncludeFlag, ExcludeFlag);
	}
	
	/*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* keepTransition, int *IncludeFlag, int *ExcludeFlag)
{
	/*We will be doing a linear programming problem for each transition such that 
	it is flagged in keepTransition, not flagged in IncludeFlag, and not flagged 
	in ExcludeFlag.
	The constraint matrix will be the subset of D having only the	columns 
	corresponding to transitions not flagged in ExcludeFlag, 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, StoredT;
	double* B;
	short int *IntList;
	ipsolve_r result;
	
	for (i = 0; i < NumberOfColumns(*D); i++)
	{
		if (! ExcludeFlag[i]) Transitions++;
	}
	//If no transitions are being evaluated, return true.
	if (! Transitions) return 1;
	
	/*Build Dx. It should include only columns of D not flagged in ExcludeFlag.*/
	AllocateMatrixType(2, &Dx, NumberOfRows(*D) + 1, Transitions);	
	j = 0;
	for (i = 0; i < NumberOfColumns(*D); i++)
	{
		
		if (! ExcludeFlag[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++);
		}
	}
	
	/*Allocate space to record the values of the current column while it is
	zeroed.*/
	AllocateMatrixType(1, &dx, NumberOfRows(*D) + 1, 1);	
	
	/*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));
	
	/*We will be zeroing and restoring certain columns. Store the index of the 
	currently-stored column in StoredT. Defaults to -1 to indicate no stored
	column.*/
	StoredT = -1;
	
	/*Iterate through the transitions*/
	for (i = 0; i < NumberOfColumns(*D); i++)
	{
		/*Skip each transition such that it is not flagged in keepTransition,
		flagged in IncludeFlag, or flagged in ExcludeFlag.*/
		if ((! keepTransition[i]) || IncludeFlag[i] || ExcludeFlag[i]) continue;
		
		/*Zero a column of Dx and store it temporarily in dx for later restoration*/
		for (j = 0; j <= NumberOfRows(*D); j++)
		{
			//Restore the previous transition (if there is one to restore)
			if (StoredT >= 0) SetMatrixEl(&Dx, j, i - 1, GetMatrixEl(&dx, j, 0));
			//Copy the current transition
			SetMatrixEl(&dx, j, 0, GetMatrixEl(&Dx, j, i));
			//Zero the current transition.
			SetMatrixEl(&Dx, j, i, 0);
		}
		/*Set the index of the newly-stored column.*/
		StoredT = i;
		
		result = ipsolve(&Dx, B, 0, IntList, 0, 0, 0);		
		/*If we get a valid solution, then free memory and return false.*/
		if (! strcmp(result.mhow, HOW_OK))
		{
			DeallocateIpsolve(&result);
			free(IntList);
			free(B);			
			DeallocateMatrix(&Dx);
			DeallocateMatrix(&dx);
			return 0;
		}
		else
		{
			DeallocateIpsolve(&result);
		}
	}
	//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 keep from the total net
and builds the active subnet. It returns the result structure that will be
returned by tactn.*/
tactn_r BuildResult(matrix* Dm, matrix* Dp, int* keepPlace, int* keepTrans)
{
	int Places = 0, Transitions = 0, i, j, k, l;
	tactn_r result;
	
	/*All matrices are zeroed initially.*/
	memset(&result, 0, sizeof(tactn_r));
	
	/*Count the number of places and transitions that are going to be left after
	culling.*/
	for (i = 0; i < NumberOfRows(*Dm); i++)
	{
		if (keepPlace[i]) Places++;
	}
	for (i = 0; i < NumberOfColumns(*Dm); i++)
	{
		if (keepTrans[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 (! keepTrans[j]) continue;
		k = 0;
		for (i = 0; i < NumberOfRows(*Dm); i++)
		{
			if (! keepPlace[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;
}
/*******************************************************************************
keepPlaces examines the Petri net and the list of transitions to keep and
returns an array of integer flags, one for each place, indicating whether that
place should be kept.*/
int* keepPlaces(matrix* Dp, int* keepTransition)
{
	int i, j;
	int *keep;
	/*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.*/
	keep = 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 being kept.*/
		for (j = 0; j < NumberOfColumns(*Dp); j++)
		{
			if (GetMatrixEl(Dp, i, j) && (keepTransition[j]))
			{				
				keep[i] = 1;
				break;
			}
		}
	}
	return keep;
}

/*******************************************************************************
maxactn finds the maximal active subnet. It functions by examining D, the
include and exclude transition flags, and building an initial keepable
transition flag array, and then passing it to minactn for processing. It returns
a pointer to the keepable transition flag array.*/
int *maxactn(matrix *D, int *IncludeFlag, int *ExcludeFlag)
{
	/*We will be building the initial keepable transition flag array. This will
	be done via a series of linear programming problems, using as constraint
	matrix a copy of the Petri net incidence matrix with certain columns nulled
	out and with a single constraint row of all ones (except for the null 
	columns).*/
	/*The flags used to indicate which transitions are to be non-nulled for the
	linear programming problem of the current iteration of the processing loop.
	We keep a copy of the flags for the previous iteration so that we can avoid
	recopying or renulling a column that was already copied or nulled.*/
	int *oldFlag, *newFlag, *keepTransition;
	/*Other parts of the linear programming problem.*/
	short *IntList;
	double *B;
	matrix L;
	ipsolve_r lpResult;
	
	/*Other variables*/
	MemoryManager mem;
	int i, j, keepGoing, DRows, DCols;	
	DRows = NumberOfRows(*D);
	DCols = NumberOfColumns(*D);
	
	/*Allocate space*/
	mem = CreateMemoryManager(5, 1, 0, 0);
	MAllocateMatrixType(&mem, 2, &L, DRows + 1, DCols);
	IntList = mcalloc(&mem, DCols, sizeof(short));
	B = mcalloc(&mem, DRows + 1, sizeof(double));
	oldFlag = mcalloc(&mem, DCols, sizeof(int));
	newFlag = mcalloc(&mem, DCols, sizeof(int));
	keepTransition = tcalloc(DCols, sizeof(int));
	
	/*Initialize the flags, B, and the preprocessing keepGoing flag. The 
	preprocessing loop runs only until all elements of the Flag arrays are 0.*/
	keepGoing = 0;
	B[DRows] = 1.0;
	for (i = 0; i < DCols; i++)
	{
		if (IncludeFlag[i] && (! ExcludeFlag[i]))
		{
			newFlag[i] = 1;
			keepGoing = 1;
		}
	}
	
	/*The preprocessing loop*/
	while (keepGoing)
	{
		/*Build the constraint matrix. Null out any columns that were filled last
		time and fill in any columns that were nulled last time.*/
		for (i = 0; i < DCols; i++)
		{
			if (newFlag[i] && (! oldFlag[i]))
			{
				CopyBlock(DRows, 1, D, 0, i, &L, 0, i);
				SetMatrixEl(&L, DRows, i, 1);
			}
			else if ((! newFlag[i]) && oldFlag[i])
			{
				for (j = 0; j <= DRows; j++)
				{
					SetMatrixEl(&L, j, i, 0);
				}
			}
		}
		
		/*Solve the linear programmig problem. Cost, bounds, and constraint types
		can all be left as defaults.*/
		lpResult = ipsolve(&L, B, 0, IntList, 0, 0, 0);
		
		/*If the solve failed, terminate preprocessing immediately.*/
		if (strcmp(lpResult.mhow, HOW_OK))
		{			
			if (lpResult.res) free(lpResult.res);
			break;
		}
		else
		/*Otherwise, set all keepable transition flags with nonzero elements in the
		solution vector and clear all transition flags corresponding to nonzero
		elements in the solution vector. Copy the new flag array to the old one
		first. Also use the loop to set the keepGoing flag only if there are any
		set flags in the new flag array.*/
		{
			keepGoing = 0;
			memcpy(oldFlag, newFlag, sizeof(int) * DCols);
			for (i = 0; i < DCols; i++)
			{
				/*Ignore solution elements from null columns and very small nonzero
				solution elements.*/
				if (newFlag[i] && (lpResult.res[i] < -HOW_MAX_ERROR || lpResult.res[i] > HOW_MAX_ERROR))
				{
					keepTransition[i] = 1;
					newFlag[i] = 0;
				}
				/*Set the keepGoung flag if a nonzero flag element is found.*/
				if (newFlag[i]) keepGoing = 1;
			}
			free(lpResult.res);
		}
	}
	/*Build a new transition-to-include flag array for the minactn call. It should
	be the logic-and of the actual transition-to-include flag array and the
	initial keepable transitions flag array. Use the newFlag array since the
	storage space is already allocated.*/
	for (i = 0; i < DCols; i++)
	{
		if (IncludeFlag[i] && keepTransition[i])
		{
			newFlag[i] = 1;
		}
		else
		{
			newFlag[i] = 0;
		}
	}
	/*Do the minactn call. This will make appropriate modifications to the
	keepTransition flags.*/
	minactn(keepTransition, D, newFlag);
	
	/*Free memory and return.*/
	FreeMemory(&mem);
	return keepTransition;	
}

/*******************************************************************************
minactn finds the minimal active subnet of the given Petri net. It takes a
pointer to the array of transition which are to be kept, initialized by the
GetKeepT0 function, which it modifies in an iterative fashion. It also takes
a pointer to the incidence matrix and the flags indicating which transitions
must be included in the minimal active subnet. It modifes the keepTransition 
flag array passed to it.*/
void minactn(int *keepTransition, matrix *D, int *IncludeFlag)
{
	matrix L;
	int i, j, k, Cols, Rows;
	int *LFlags;
	short *IntList;
	double *B;
	double *LoBounds;
	ipsolve_r lpResult;
	Cols = NumberOfColumns(*D);
	Rows = NumberOfRows(*D);
	MemoryManager mem;
	/*We will be iterating through each transition in D and doing a linear
	programming problem. The constraint matrix will be a copy of D with
	constraint columns for each zeroed element in keepTransition and for
	the current column removed or zeroed.*/
	/*Allocate space and initialize.*/
	mem = CreateMemoryManager(5, 1, 0, 0);
	MAllocateMatrixType(&mem, 2, &L, Rows, Cols);	
	IntList = mcalloc(&mem, Cols, sizeof(short));
	B = mcalloc(&mem, Rows, sizeof(double));
	LoBounds = mcalloc(&mem, Cols, sizeof(double));	
	/*Lo bounds are 0 for columns not flagged in IncludeFlag and 1 for the rest.*/
	for (i = 0; i < Cols; i++)
	{
		if (IncludeFlag[i]) LoBounds[i] = 1.0;
	}
	/*This flag array indicates whether a particular column of L has been filled
	(flag set) or nulled (flag cleared).*/
	LFlags = mcalloc(&mem, Cols, sizeof(int));
	
	for (i = 0; i < Cols; i++)
	{
		/*We only process a subset of transitions - those which have not already
		been marked as unkeepable and those which are NOT flagged in IncludeFlag.*/
		if (keepTransition[i] && (! IncludeFlag[i]))
		{
			/*Zero or fill columns of the constraint matrix as appropriate.*/
			for (j = 0; j < Cols; j++)
			{
				/*We zero the current column and any column not flagged in
				keepTransition.*/
				if (((! keepTransition[j]) || i == j))
				{
					/*Zero the column only if it has not already been zeroed.*/
					if (LFlags[j])
					{
						for (k = 0; k < Rows; k++)
						{
							SetMatrixEl(&L, k, j, 0);
						}
						LFlags[j] = 0;
					}					
				}
				/*Otherwise, the column should be filled if it is not already full.*/
				else if (! LFlags[j])
				{
					CopyBlock(Rows, 1, D, 0, j, &L, 0, j);
					LFlags[j] = 1;
				}
			}
			/*Solve the linear programming problem. Cost defaults to all 1s, upper
			bounds default to infinity, and constraint type defaults to >=. No 
			variables should be restricted to being integers.*/
			lpResult = ipsolve(&L, B, 0, IntList, 0, LoBounds, 0);
			
			/*If the linear programming problem was infeasible, break, leave the 
			keepTransition flag array as-is, and return immediately.*/
			if (strcmp(lpResult.mhow, HOW_OK))
			{
				if (lpResult.res) free(lpResult.res);
				break;
			}
			/*Otherwise, set the keepTransition flag array to the logic equivalent of
			the lp solution.*/
			else
			{
				memset(keepTransition, 0, sizeof(int) * Cols);
				for (j = 0; j < Cols; j++)
				{
					/*Ignore very small nonzero values in the solution vector. Also
					examine only solution elements corresponding to non-nulled columns,
					as noted in LFlags.*/
					if ((lpResult.res[j] < -HOW_MAX_ERROR || lpResult.res[j] > HOW_MAX_ERROR) && LFlags[j])
					{
						keepTransition[j] = 1;
					}
				}
				free(lpResult.res);
			}
		}
	}
	FreeMemory(&mem);
}

/*******************************************************************************
GetKeepT0 is called in non-update mode. It examines the incidence matrix and
excluded transition list and determines whether to call MaxActn or MinActn to
get the list of transitions in the active subnet. If minactn is to be called,
it returns a pointer to the initial keepable transition flag array that must be 
given as a parameter to minactn. Otherwise it returns a null pointer.*/
int* GetKeepT0(matrix *D, int *ExcludeFlag, int *IncludeFlag, int ExcludedTCount)
{
	/*We will be doing a linear programming problem.*/
	matrix L; // The constraint matrix	
	double *LoBound; // Lower bounds on variables - the upper bounds will default.
	double *B; // Constraint vector
	short *IntList; //The list of which variables in the LP problem must be integral
	int VarCount, ConCount, i, j;
	MemoryManager mem;
	ipsolve_r lpresult;
	int* result;
	/*We are using as constraints only those columns of D such that they are not
	excluded by ExcludedT. That will be equal to the total transition count less
	the excluded transition count.*/
	VarCount = NumberOfColumns(*D) - ExcludedTCount;
	ConCount = NumberOfRows(*D);
	
	/*Allocate space. B and IntList will be filled with zeroes by the allocation
	and since this is the final value they should have we won't have to do any
	further processing.*/
	mem = CreateMemoryManager(3, 1, 0, 0);
	MAllocateMatrixType(&mem, 2, &L, ConCount, VarCount);	
	LoBound = mcalloc(&mem, VarCount, sizeof(double));
	B = mcalloc(&mem, ConCount, sizeof(double));
	IntList = mcalloc(&mem, VarCount, sizeof(short));
	
	/*Initialize variable boundaries and build the constraint matrix. The upper 
	bounds will default to infinity. The lower bounds should be 1 for each 
	transition such that it is included in IncludedT and 0 for the rest. Iterate 
	through the flags in the loop. There is	an element in LoBound for each zeroed 
	element in ExcludedFlag. We copy to L each column of D such that it is not
	excluded by ExcludeFlag.*/
	j = 0;
	for (i = 0; i < NumberOfColumns(*D); i++)
	{
		if (! ExcludeFlag[i])
		{
			/*Copy the column*/
			CopyBlock(ConCount, 1, D, 0, i, &L, 0, j);
			/*Set the lower boundary if necessary */
			if (IncludeFlag[i])
			{
				LoBound[j] = 1.0;
			}
			/*Increment the constraint index.*/
			j++;
		}
	}
	
	/*Do the linear programming problem. The original matlab implementation used
	linprog, which accepts only <= constraints, and negated the constraint vector
	and matrix. ipsolve defaults to >= constraints, so the constraints vector and
	matrix will not need to be negated. Cost vector, upper bounds, and constraint
	type should all be left as default.*/
	lpresult = ipsolve(&L, B, 0, IntList, 0, LoBound, 0);

	if (lpresult.res) ManageMemory(&mem, lpresult.res);
	/*If the problem was not solved, we'll be doing a maxactn to find the active
	subnet. Return a null pointer since the X0 vector will not be needed.*/
	if (strcmp(lpresult.mhow, HOW_OK))
	{
		FreeMemory(&mem);
		return 0;
	}
	else
	{
		/*Otherwise, prepare the initial keepTransition vector. This should have an
		element for each transition. Each value corresponding to a non-excluded
		transition should be set to one if the corresponding lp solution is nonzero.
		All the rest should be zero.*/
		result = tcalloc(NumberOfColumns(*D), sizeof(int));
		j = 0;
		for (i = 0; i < NumberOfColumns(*D); i++)
		{
			if (! ExcludeFlag[i])
			{
				/*Ignore very small nonzero solution values.*/
				if (lpresult.res[j] < -HOW_MAX_ERROR || lpresult.res[j] > HOW_MAX_ERROR)
				{
					result[i] = 1;
				}
				j++;
			}
		}
		FreeMemory(&mem);
		return result;
	}
}

/*******************************************************************************
updateActn examines the Petri net and the predefined active subnet if given and
returns an array of transitions that are to be kept given the active subnet.
This is used to determine transitions to keep if in update mode.
*/
int* updateActn(matrix* Dm, matrix* Dp, matrix* Dcm, matrix* 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 *Actn, *Trans = 0, *Place = 0;	
	
	Places = NumberOfRows(*Dm);
	Transitions = NumberOfColumns(*Dm);
	SubPlaces = NumberOfRows(*Dcm);
	SubTransitions = NumberOfColumns(*Dcm);
	
	/*Actn is the return value. It has a flag for each transition.*/
	Actn = 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));
	}
		
	/*Actn is a flag vector with an element for each transition. It is initially
	set 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 not null, set the Actn flag.*/
		if (j != SubPlaces)
		{
			Actn[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 flagged in Actn.*/
				for (j = 0; j < SubTransitions; j++)
				{
					if (Actn[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 set every flag of Actn 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;
					/*Set the appropriate Actn flag.*/
					Actn[j] = 1;
					/*Set the keepGoing flag.*/
					keepGoing = 1;
					break;
				}
			}
		}
	} while (keepGoing);
	/*Free memory*/
	if (Place) free(Place);
	if (Trans) free(Trans);
	/*Return*/
	return Actn;
}

/*******************************************************************************
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 **IncludeT, int *IncludeTCount, int ** ExcludeT, int *ExcludeTCount, matrix *Dcm, matrix *Dcp, MemoryManager* mem)
{
	int Rows0, Cols0, Rows1, Cols1, i;
	/*Dm and Dp are required and must be the same size.*/
	if (! (Dm && Dp))
	{
		merror(0, "TACTN: Dm or Dp was not given");
		return 0;
	}
	Rows0 = NumberOfRows(*Dm);
	Cols0 = NumberOfColumns(*Dm);
	if (Rows0 != NumberOfRows(*Dp) || Cols0 != NumberOfColumns(*Dp))
	{
		merror(0, "TACTN: Dm and Dp are not the same size");
		return 0;
	}
	
	/*ExcludeT is not required. Make sure that the count is zero if the pointer is 
	null and vice-versa.*/
	if (! *ExcludeT)
	{
		*ExcludeTCount = 0;
	}
	else if (! *ExcludeTCount)
	{
		*ExcludeT = 0;
	}
	
	/*IncludeT is not required. Make sure the count is zero if the pointer is null
	and vice-versa.*/
	if (! *IncludeT)
	{
		*IncludeTCount = 0;
	}
	else if (! *IncludeTCount)
	{
		*IncludeT = 0;
	}
	/*If IncludeT was not given it defaults to including the indices of every
	transitions.*/
	if (! *IncludeT)
	{
		*IncludeT = mcalloc(mem, Cols0, sizeof(int));
		*IncludeTCount = Cols0;
		for (i = 1; i < Cols0; i++)
		{
			(*IncludeT)[i] = i;
		}
	}
	/*Update mode will be entered if Dcm and Dcp are present. If just one is
	present it will be ignored. If both are present they must be the same size.*/
	if (Dcm && Dcp)
	{
		Rows1 = NumberOfRows(*Dcm);
		Cols1 = NumberOfColumns(*Dcm);
		if (NumberOfRows(*Dcp) != Rows1 || NumberOfColumns(*Dcp) != Cols1)
		{
			merror(0, "TACTN: Dcm and Dcp are not the same size");
			return 0;
		}
		if (Rows1 > Rows0 || Cols1 > Cols0)
		{
			merror(0, "TACTN: Dcm and Dcp are larger than Dm and Dp");
			return 0;
		}
	}
	return 1;
}
