#include <stdio.h>
#include "../pnheaders/general.h"
#include "../pnheaders/matrix.h"
#include "spnbox.h"
#include "MemoryManager.h"
#include "extendedmatrix.h"

typedef struct asiph_data
{
	int *numInputPlaces;
	int *curTrans;
	int *counters;
	matrix M;
	matrix N;
	matrix MT;
	matrix NT;
	matrix newRow;
	matrix newPlace;	
	matrix MNS;
	matrix NS;
} asiph_data;

static int CheckParams(matrix *Dm, matrix *Dp, int **list, int *listCount, matrix *PS, matrix **Dma, matrix **Dpa, MemoryManager *mem);
static matrix SXInit(matrix *Dm, matrix *Dp, matrix *Dma, matrix* Dpa, matrix *PS, int **list, int* listCount, asiph_data* data);
static int* insert(asiph_data* data, matrix* vector, int* index);
static int *DataInit(matrix *Dm, matrix *Dp, matrix *PS, matrix *sx, asiph_data* data);
static int GetNewTrans(matrix *Dm, matrix *Dp, matrix *Places, int *Flags);
static void OuterLoopInit(asiph_data* data, matrix *Dm, matrix *Dp, MemoryManager *mem);
static void MiddleLoopInit(asiph_data* data, asiph_r* result, int siphon);
static void InnerLoopInit(asiph_data* data, matrix *Dm, matrix *Dp);
static int Check1(matrix *MNS, matrix *newRow);
static int Check2(matrix* M, matrix* newRow);
static int* DoInnerLoop(asiph_data* data, matrix* Dm, matrix *Dp, int* index);
static void PrintPercent(int num, int denom);
static asiph_r BuildResult(matrix *Dm, matrix *Dp, asiph_data* data, int *index);

asiph_r asiph(matrix *Dm, matrix *Dp, int* list, int listCount, matrix *PS, matrix *Dma, matrix *Dpa)
{
	MemoryManager mem;
	mem = CreateMemoryManager(10, 10, 0, 0);
	asiph_r result;
	memset(&result, 0, sizeof(asiph_r));
	/*Check parameters.*/
	if (! CheckParams(Dm, Dp, &list, &listCount, PS, &Dma, &Dpa, &mem))
	{
		return result;
	}
	asiph_data data;
	memset(&data, 0, sizeof(asiph_data));
	
	/*Initialize sx (and a few other things)*/
	matrix sx = SXInit(Dm, Dp, Dma, Dpa, PS, &list, &listCount, &data);
	ManageMatrix(&mem, &sx);
	/*Initialize some other parts of the loop data. Build the newIndex array.*/
	int *newIndex = DataInit(Dm, Dp, PS, &sx, &data);
	
	/*Allocate and initialize various quantities that do not change or the space
	for which does not change over the iterations of the outer loop.*/
	OuterLoopInit(&data, Dm, Dp, &mem);
	
	/*Iterate through each of the rows of NS. Each row corresponds to a siphon
	that needs to be processed. If verbosity is high enough we'll be printing a
	progress report as we go. Prepare.*/
	int siphon, originalSiphons;
	originalSiphons = NumberOfRows(data.NS);
	if (is_verbose() >= VRB_ASIPH)
	{
		printf("ASIPH: Primary Processing Loop Progress: ");
		PrintPercent(0, 1);
	}
	for (siphon = 0; siphon < originalSiphons; siphon++)
	{
		/*Initialize the stuff that doesn't change within the middle loop.*/
		MiddleLoopInit(&data, &result, siphon);
		/*The middle loop runs until we have removed all rows from M/N*/
		while (NumberOfRows(data.M))
		{
			/*Initialize the stuff for the innermost loop run.*/
			InnerLoopInit(&data, Dm, Dp);
			/*Run the inner loop. This has been moved to its own function.*/
			newIndex = DoInnerLoop(&data, Dm, Dp, newIndex);
			/*Now concatenate MT onto M and NT onto N.*/
			if (NumberOfRows(data.MT)) InsertRows(&data.M, &data.MT, -1, -1);
			if (NumberOfRows(data.NT)) InsertRows(&data.N, &data.NT, -1, -1);
			/*Remove the first row of M and N.*/
			RemoveRows(&data.M, 0, 1, -1);
			RemoveRows(&data.N, 0, 1, -1);
		}
		/*Deallocate whatever still needs to be deallocated - M and N.*/
		DeallocateMatrix(&data.M);
		DeallocateMatrix(&data.N);
		/*If verbosity is high enough, print a progress percentage.*/
		if (is_verbose() >= VRB_ASIPH)
		{
			PrintPercent(siphon + 1, originalSiphons);
		}
	}
	if (is_verbose() >= VRB_ASIPH) 
	{
		PrintPercent(1, 1);
		printf("\n");
	}
	/*If MNS is empty, build it based on if there are any source transitions.*/
	result = BuildResult(Dm, Dp, &data, newIndex);
	FreeMemory(&mem);
	return result;
}

/*******************************************************************************
This function builds the result structure.*/
asiph_r BuildResult(matrix *Dm, matrix *Dp, asiph_data* data, int *index)
{
	int i, j, oFlag, iFlag, flag = 0;
	int *sourceF;
	sourceF = tcalloc(NumberOfColumns(*Dm), sizeof(int));
	
	/*Build a set of source transition flags.*/
	
	/*If MNS is empty, build a new MNS algorithmically.*/
	if (! NumberOfRows(data->MNS))
	{
		/*Find out if there are any source transitions - transitions such that they
		have both input and output arcs.*/
		int oFlag, iFlag;
		for (j = 0; j < NumberOfColumns(*Dm); j++)
		{
			iFlag = 0;
			oFlag = 0;
			for (i - 0; i < NumberOfRows(*Dm); i++)
			{
				if (GetMatrixEl(Dm, i, j)) oFlag = 1;
				if (GetMatrixEl(Dp, i, j)) iFlag = 1;				
			}
			if (iFlag && oFlag)
			{
				sourceF[j] = 1;
				flag = 1;
			}
		}
		/*If there exist source transitions, MNS should be set to a zero-by-m matrix,
		m the number of places. Otherwise, it should be set to a 1-m matrix of 1s.*/
		if (flag)
		{
			if (data->MNS.type) DeallocateMatrix(&data->MNS);
			AllocateMatrixType(2, &data->MNS, 0, NumberOfRows(*Dm));
		}
		else
		{
			if (data->MNS.type) DeallocateMatrix(&data->MNS);
			AllocateMatrixType(2, &data->MNS, 1, NumberOfRows(*Dm));
			for (i = 0; i < NumberOfRows(*Dm); i++)
			{
				SetMatrixEl(&data->MNS, 0, i, 1);
			}
		}
		/*Build index*/
		if (index) free(index);
		index = tcalloc(NumberOfRows(*Dm), sizeof(int));
		for (i = 0; i < NumberOfRows(*Dm); i++)
		{
			index[i] = 1;
		}
	}
	/*The result NMS should be a transpose copy of MNS with columns not flagged
	in newindex removed. The result MS should be a transpose copy of MNS. Count
	how many columns that is. Both of the result matrices should be optimized
	for column operations.*/
	j = 0;
	for (i = 0; i < NumberOfRows(data->MNS); i++)
	{
		if (index[i]) j++;
	}
	
	asiph_r result;
	AllocateMatrixType(2, &result.NMS, j, NumberOfColumns(data->MNS));
	AllocateMatrixType(2, &result.MS, NumberOfRows(data->MNS), NumberOfColumns(data->MNS));
	TransposeMatrix(&result.NMS);
	TransposeMatrix(&result.MS);
	
	TransposeMatrix(&data->MNS);
	CopyMatrix(&data->MNS, &result.MS);
	j = 0;
	for (i = 0; i < NumberOfColumns(data->MNS); i++)
	{
		if (index[i])
		{
			CopyBlock(NumberOfRows(data->MNS), 1, &data->MNS, 0, i, &result.NMS, 0, j++);
		}
	}
	if (data->MNS.type) DeallocateMatrix(&data->MNS);
	if (data->NS.type) DeallocateMatrix(&data->NS);
	free(sourceF);
	free(index);
	return result;
}
				
/*******************************************************************************
This function displays a real-time (changing) progress percentage. Start by
calling it with a numerator of zero to pring the initial zero-percent text.
Every call for which the division num / denom is nonzero overwrites the previously
displayed percentage - or the last four characters if no percentage has been
displayed yet.*/
void PrintPercent(int num, int denom)
{
	static int lastPercent = -1;
	int percent;
	percent = (num * 100) / denom;
	if (percent != lastPercent)
	{
		lastPercent = percent;
		if (percent) printf("\b\b\b\b");
		printf("%3d%%", percent);
		fflush(stdout);
	}	
}

/*******************************************************************************
This function is the innermost processing loop. This adds rows to MT and NT, the
tentative siphons to add to the M and N temporary siphon tables.*/
int* DoInnerLoop(asiph_data* data, matrix* Dm, matrix *Dp, int *index)
{
	int flag = 1, place, trans, nullPlace, places;
	while (flag)
	{
		/*Allocate the newRow matrix.*/
		AllocateMatrixType(2, &data->newRow, 1, NumberOfRows(*Dm));
		
		/*This next section replaces the building of plc and the filling of tst in
		the original matlab algorithm. Essently, we want to iterate through all
		transitions that have been flagged in curTrans. For each one, we want to
		find the index of the counter[i]-th place that has an output arc to the 
		transition, where i is the transition index. If there is no such place, we
		set the nullPlace flag. Otherwise, we set the newRow[index] flag.*/
		nullPlace = 0;
		for (trans = 0; trans < NumberOfColumns(*Dm); trans++)
		{
			if (data->curTrans[trans])
			{
				places = 0;
				for (place = 0; place < NumberOfRows(*Dm); place++)
				{
					if (GetMatrixEl(Dm, place, trans))
					{
						if (++places == data->counters[trans]) break;
					}
				}
				/*If place is less than the upper limit, then it is the index we want.
				Otherwise the index we want does not exist.*/
				if (place < NumberOfRows(*Dm))
				{
					SetMatrixEl(&data->newRow, 0, place, 1);
				}
				else
				{
					nullPlace = 1;
				}
			}
		}
		
		/*If the nullPlace flag is set, we found some ith transition without a
		counter[i]-th output arc. Free up memory used by the newRow matrix and the
		newPlace matrix.*/
		if (nullPlace)
		{
			if (data->newRow.type) DeallocateMatrix(&data->newRow);
			if (data->newPlace.type) DeallocateMatrix(&data->newPlace);
		}
		/*Otherwise, we can go ahead and add to MT and NT.*/
		else
		{
			/*Get the newPlaces matrix. There is an element of this for each place,
			and it is set to the logic-and of newRow and (logic-not 1st row M).
			While we're at it, logic-or newRow with the first row of M.*/
			AllocateMatrixType(2, &data->newPlace, 1, NumberOfRows(*Dm));
			for (place = 0; place < NumberOfRows(*Dm); place++)
			{
				/*Set newPlace to newRow && ! M(1, :)*/
				SetMatrixEl(&data->newPlace, 0, place, (GetMatrixEl(&data->newRow, 0, place) && ! GetMatrixEl(&data->M, 0, place)) ? 1 : 0);
				/*Set newRow to newRow | M(1, :)*/
				if (GetMatrixEl(&data->M, 0, place))
				{
					SetMatrixEl(&data->newRow, 0, place, 1);
				}
			}
			
			/*If check1 returns false, we won't be modifying anything.*/
			if (Check1(&data->MNS, &data->newRow))
			{
				/*If Check1 returns true, how we modify depends on if there are any
				transitions that meet the GetNewTrans criteria, given newRow as the
				place flags. If there are such transitions, run Check2 and modify MT
				and NT. If not, modify MNS and newIndex.*/
				if (GetNewTrans(Dm, Dp, &data->newRow, 0))
				{
					if (Check2(&data->M, &data->newRow))
					{
						InsertRows(&data->MT, &data->newRow, -1, -1);
						InsertRows(&data->NT, &data->newPlace, -1, -1);
					}
					else
					{
						DeallocateMatrix(&data->newRow);
						DeallocateMatrix(&data->newPlace);
					}
				}			
				else
				{
					/*Clear the loop flag. It may be set again later, but it may not be.*/
					flag = 0;
					/*Modify MNS and newIndex via insert. The row vector to use is newRow*/
					index = insert(data, &data->newRow, index);
					/*Deallocate newPlace.*/
					DeallocateMatrix(&data->newPlace);
				}
			}
			else
			{
				/*Deallocate the unused matrices.*/
				DeallocateMatrix(&data->newRow);
				DeallocateMatrix(&data->newPlace);
			}
		}
		/*Increment the counters and adjust the flag. We are examining only counters
		for which curTrans is set.*/
		int i;
		for (i = 0; i < NumberOfColumns(*Dm); i++)
		{
			if (data->curTrans[i])
			{
				/*The limit on each counter is the number of input places found earlier.*/
				if (data->counters[i] < data->numInputPlaces[i])
				{
					data->counters[i]++;
					flag = 1;
					break;
				}
				else
				{
					data->counters[i] = 0;
					flag = 0;
				}
			}
		}	
	}
	return index;
}

/*******************************************************************************
Check type 2. This runs a check on M and newRow and returns a logical result.*/
int Check2(matrix* M, matrix* newRow)
{
	/*This check is the same as Check1 except that the second half, examining 
	rows of M, will return false if there exists any row of M such that it
	equals newRow exactly.*/
	int i, j;
	/*If all elements of newRow are 1, return false.*/
	for (i = 0; i < NumberOfColumns(*M); i++)
	{
		if (GetMatrixEl(newRow, 0, i) != 1) break;
	}
	if (i == NumberOfColumns(*M)) return 0;
	
	/*Check each row of M. If there exists a row of M that is equal to newRow,
	return false.*/
	for (i = 0; i < NumberOfRows(*M); i++)
	{
		for (j = 0; j < NumberOfColumns(*M); j++)
		{
			if (GetMatrixEl(M, i, j) != GetMatrixEl(newRow, 0, j)) break;
		}
		if (j == NumberOfColumns(*M))
		{
			return 0;
		}
	}
	return 1;
}
	


/*******************************************************************************
Check type 1. This runs a check on MNS and newRow and returns a logical result.*/
int Check1(matrix *MNS, matrix *newRow)
{
	int i, j;
	/*If all elements of newRow are 1, return false.*/
	for (i = 0; i < NumberOfColumns(*MNS); i++)
	{
		if (GetMatrixEl(newRow, 0, i) != 1) break;
	}
	if (i == NumberOfColumns(*MNS)) return 0;
	
	/*Check each row of MNS. If there exists a row of MNS such that it has no
	element that is true when the corresponding element in newRow is false, return
	false.*/
	for (i = 0; i < NumberOfRows(*MNS); i++)
	{
		for (j = 0; j < NumberOfColumns(*MNS); j++)
		{
			if (GetMatrixEl(MNS, i, j) && ! GetMatrixEl(newRow, 0, j)) break;
		}
		if (j == NumberOfColumns(*MNS))
		{
			return 0;
		}
	}
	return 1;
}

/*******************************************************************************
This function sets up all parts of the loop data structure that need to be 
initialized for the innermost loop.*/
void InnerLoopInit(asiph_data* data, matrix *Dm, matrix *Dp)
{
	/*We'll need to get a transition flag array and store it in curTrans in
	the loop data structure. Use the algorithm in GetNewTrans, with the first
	row of M as the place flags.*/
	GetNewTrans(Dm, Dp, &data->M, data->curTrans);
	
	/*MT and NT get initialized to matrices with no rows and as many columns as
	there are places. They should be set up for fast row ops.*/
	AllocateMatrixType(2, &data->MT, 0, NumberOfRows(*Dm));
	AllocateMatrixType(2, &data->NT, 0, NumberOfRows(*Dm));
	
	/*Initialize the counters to 1.*/
	int i;
	for (i = 0; i < NumberOfColumns(*Dm); i++)
	{
		data->counters[i] = 1;
	}
}
	
/*******************************************************************************
This function sets up all parts of the loop data structure that need to be
initialized for the middle loop. This is M and N.*/
void MiddleLoopInit(asiph_data* data, asiph_r* result, int siphon)
{
	/*Initialize M and N each to be a copy of the current row of NS, that is, the
	current siphon. M and N should be type-2 untransposed for optimized row ops.*/
	AllocateMatrixType(2, &data->M, 1, NumberOfColumns(data->NS));
	AllocateMatrixType(2, &data->N, 1, NumberOfColumns(data->NS));
	CopyBlock(1, NumberOfColumns(data->NS), &data->NS, siphon, 0, &data->M, 0, 0);
	CopyBlock(1, NumberOfColumns(data->NS), &data->NS, siphon, 0, &data->N, 0, 0);
}

/*******************************************************************************
This function sets up all parts of the loop data structure that do not change
ever. This includes two allocations and one allocation and data setup.*/
void OuterLoopInit(asiph_data* data, matrix *Dm, matrix *Dp, MemoryManager *mem)
{
	/*There is one quantity to be allocated and inialized and two more to be 
	allocated. All are integer arrays with an element for each transition.
	Allocate:*/
	data->numInputPlaces = mcalloc(mem, NumberOfColumns(*Dp), sizeof(int));
	data->curTrans = mcalloc(mem, NumberOfColumns(*Dp), sizeof(int));
	data->counters = mcalloc(mem, NumberOfColumns(*Dp), sizeof(int));
	
	/*numInputPlaces has an element for each transition, showing the number of
	output arcs that go to that transition.*/
	int i, j, num;
	for (j = 0; j < NumberOfColumns(*Dm); j++)
	{
		for (i = 0; i < NumberOfRows(*Dm); i++)
		{
			if (GetMatrixEl(Dm, i, j)) data->numInputPlaces[j]++;
		}
	}
}

/*******************************************************************************
This function performs a logic operation to develop a set of new transition
flags. It examines places given in the matrix Places and fills in flags in the
integer array Flags. It returns the number of flags set in Flags.
The matrix Places should be a matrix with as many columns as there are places.
The first row will be interpreted as a flag array. 
The array Flags should have an element for each transition. If Flags is a null 
pointer it is ignored and the function optimizes by returning as soon as any 
transition meeting the criteria is encountered.*/
int GetNewTrans(matrix *Dm, matrix *Dp, matrix *Places, int *Flags)
{
	/*We want to set the flag for each transition such that it has at least one
	input arc going to one of the flagged places and no output arcs coming from any
	of the flagged places.*/
	int transition, place, iFlag, oFlag, flagged;
	flagged = 0;
	for (transition = 0; transition < NumberOfColumns(*Dm); transition++)
	{
		iFlag = 0;
		oFlag = 0;
		for (place = 0; place < NumberOfRows(*Dm); place++)
		{
			/*Examine only if the current place is flagged.*/
			if (GetMatrixEl(Places, 0, place))
			{				
				/*If there is an input arc, set the iFlag.*/
				if (GetMatrixEl(Dp, place, transition))
				{
					iFlag = 1;
				}
				/*If there is an output arc, set the oFlag and break.*/
				if (GetMatrixEl(Dm, place, transition))
				{
					oFlag = 1;
					break;
				}
			}
		}
		/*Does the transition meet the criteria?*/
		if (iFlag && ! oFlag)
		{
			/*If Flags is present, set it and increment the counter. If not, we are
			in fast mode. Merely return 1 now.*/
			if (Flags)
			{
				Flags[transition] = 1;
				flagged++;
			}
			else
			{
				return 1;
			}
		}
		else if (Flags)
		{
			Flags[transition] = 0;
		}
	}
	return flagged;
}

/*******************************************************************************
This function does initial building on MNS and NS. It also builds the initial
newIndex vector and returns it to the main function for later use.*/
int *DataInit(matrix *Dm, matrix *Dp, matrix *PS, matrix *sx, asiph_data* data)
{
	/*newIndex defaults to an array of zeros with one element for each column of
	PS. If PS is empty (null pointer), newIndex becomes a null pointer as well.*/
	int *newIndex = PS ? tcalloc(NumberOfColumns(*PS), sizeof(int)) : 0;
	
	/*We are going to iterate through the columns of sx. Each column is an active
	siphon. For each siphon, we determine if there are any transitions such that
	they have at least one input arc to a place in the siphon but no output arcs
	from any places in the siphon. If there are, we modify NS. If not, we modify 
	MNS and newIndex.*/
	int siphon, place, transition, iFlag, oFlag;
	for (siphon = 0; siphon < NumberOfColumns(*sx); siphon++)
	{
		/*We need a row vector containing a copy of the current siphon column.
		Build it.*/
		matrix rowSX;
		AllocateMatrixType(2, &rowSX, 1, NumberOfRows(*sx));
		for (place = 0; place < NumberOfRows(*sx); place++)
		{
			SetMatrixEl(&rowSX, 0, place, GetMatrixEl(sx, place, siphon));
		}
		
		/*Find out whether the transitions meet the criteria via GetNewTrans. It will
		return 1 if it finds any transitions that meet the criteria.*/
		if (GetNewTrans(Dm, Dp, &rowSX, 0))
		{
			/*We are going to add a row to NS consisting of the current column of sx.*/
			/*Build a matrix (type-2 untransposed for optimization) and do an 
			insertion. This row insertion will be optimized)*/
			InsertRows(&data->NS, &rowSX, -1, -1);
		}
		/*If no siphons met the criteria, use ins to modify MNS and newIndex.*/
		else
		{
			newIndex = insert(data, &rowSX, newIndex);
		}
	}
	return newIndex;
}

/*******************************************************************************
This function modifies MNS and newindex based on a given row vector.*/
int* insert(asiph_data* data, matrix* vector, int* index)
{
	/*We need an array of flags, one for each row of the MNS vector. The flags
	will be set if that row should be removed.*/
	int *removeMNSRow = tcalloc(NumberOfRows(data->MNS), sizeof(int));
	int removeCount = 0;
	/*Iterate through each row of MNS*/
	int i, j, k, flag = 1;
	for (i = 0; i < NumberOfRows(data->MNS); i++)
	{
		/*Take the difference MNS - vector. Fill sumDeltas with the sum of the 
		differences at each element. Set the minusDeltas flag if any of the differences
		are negative.*/
		int sumDeltas = 0, minusDeltas = 0, j, delta;
		for (j = 0; j < NumberOfColumns(data->MNS); j++)
		{
			delta = GetMatrixEl(&data->MNS, i, j) - GetMatrixEl(vector, 0, j);
			if (delta < 0) minusDeltas = 1;
			sumDeltas += delta;
		}
		/*If any of the differences were negative, do nothing.*/
		if (! minusDeltas)
		{
			/*If the flag is still set...*/
			if (flag)
			{
				/*If the sum of the deltas is nonzero, set the current row of MNS to vector
				and set the index flag.*/
				if (sumDeltas)
				{
					CopyBlock(1, NumberOfColumns(data->MNS), vector, 0, 0, &data->MNS, i, 0);
					index[i] = 1;			
				}
				/*Clear the flag*/
				flag = 0;
			}
			/*If the flag isn't still set, mark the current MNS row for removal and
			clear the index flag.*/
			else
			{
				removeCount++;
				removeMNSRow[i] = 1;
				index[i] = 0;
			}
		}
	}
	/*Initialize space for the new index. We'll be removing all the elements
	flagged in removeMNSRow. We'll be adding an element iff flag is still set.*/
	int* newIndex = tcalloc(NumberOfRows(data->MNS) - removeCount + flag, sizeof(int));
	
	/*Fill in the new index and remove rows from MNS.*/
	j = NumberOfRows(data->MNS) - removeCount;
	/*if (! mwIsSafeAddr(index, NumberOfRows(data->MNS) * sizeof(int)))
	{
		merror(0, "BAD MEMORY!");
	}*/
	for (i = NumberOfRows(data->MNS) - 1; i >= 0; i--)
	{
		if (removeMNSRow[i])
		{
			RemoveRows(&data->MNS, i, 1, -1);
		}
		else
		{
			newIndex[--j] = index[i];
		}
	}
	/*Fill in the last element of index and the last row of MNS if necessary.*/
	if (flag)
	{
		InsertRows(&data->MNS, vector, -1, -1);
		newIndex[NumberOfRows(data->MNS) - 1] = 1;
	}
	else
	{
		/*If vector's memory was not merged into MNS in the row insertion, deallocate
		vector here for consistency's sake.*/
		DeallocateMatrix(vector);
	}
	/*Deallocate the flag array.*/
	free(removeMNSRow);
	/*Free the old index.*/
	free(index);
	/*Return the new index.*/
	return newIndex;	
}

/*******************************************************************************
This function finds the precomputed set of active siphons that are present in
the active subnet. It also performs some initialization of result.*/
matrix SXInit(matrix *Dm, matrix *Dp, matrix *Dma, matrix* Dpa, matrix *PS, int **list, int* listCount, asiph_data* data)
{
	matrix sx;
	int i, j, k, newCount;
	int *newList;
	MemoryManager memory = CreateMemoryManager(2, 2, 0, 0);
	asiph_r recursive;
	/*We need to build a list of the indices of places that are in the active
	subnet (have a non-null row in Dma/Dpa). Count them.*/
	newCount = NumberOfRows(*Dma);
	for (i = 0; i < NumberOfRows(*Dma); i++)
	{
		for (j = 0; j < NumberOfColumns(*Dma); j++)
		{
			if (GetMatrixEl(Dma, i, j) || GetMatrixEl(Dpa, i, j)) break;
		}
		/*Was the row null?*/
		if (j == NumberOfColumns(*Dma))
		{
			newCount--;
		}
	}
	/*If the entire net is inactive, then asiph should return empty matrices.
	To signify this, return an empty SXInit matrix.*/
	if (! newCount)
	{
		memset(&sx, 0, sizeof(matrix));
		FreeMemory(&memory);
		return sx;
	}
	/*Otherwise, build the actual list of non-null Dma/Dpa row indices.*/
	k = 0;
	newList = mcalloc(&memory, newCount, sizeof(int));
	for (i = 0; i < NumberOfRows(*Dma); i++)
	{
		for (j = 0; j < NumberOfColumns(*Dma); j++)
		{
			if (GetMatrixEl(Dma, i, j) || GetMatrixEl(Dpa, i, j)) break;
		}
		/*Was the row non-null?*/
		if (j < NumberOfColumns(*Dma))
		{
			newList[k++] = i;
		}
	}	

	/*If the entire net is active, then sx is simply a modified identity.*/
	if (MatrixEqualMatrix(Dm, Dma) && MatrixEqualMatrix(Dp, Dpa))
	{
		/*There will be a column of sx for each place the index of which is present 
		in the list. Each of these will have exactly one 1, the row number of which
		is the same as the original place number.*/
		AllocateMatrixType(2, &sx, NumberOfRows(*Dm), *listCount);
		for (i = 0; i < *listCount; i++)
		{
			SetMatrixEl(&sx, (*list)[i], i, 1);
		}
	}
	/*Otherwise, we use the new list and build a new known PS and call asiph 
	recursively.*/
	else
	{		
		matrix NewPS;
		if (PS)
		{
			/*NewPS should be a copy of PS, with rows corresponding to inactive places
			(null rows in Dma/Dpa) zeroed and with duplicate columns removed.
			Allocate as a type-2 transpose to optimize for rapid column removal (see
			RemoveColumns).*/
			MAllocateMatrixType(&memory, 2, &NewPS, NumberOfColumns(*PS), NumberOfRows(*PS));
			TransposeMatrix(&NewPS);
			CopyMatrix(PS, &NewPS);
			/*Zero rows corresponding to zero rows of both Dma and Dpa.*/
			for (i = 0; i < NumberOfRows(*Dma); i++)
			{
				for (j = 0; j < NumberOfColumns(*Dma); j++)
				{
					if (GetMatrixEl(Dma, i, j) || GetMatrixEl(Dpa, i, j))
					{
						break;
					}
				}
				/*Was the row zero?*/
				if (j == NumberOfColumns(*Dma))
				{
					/*If so, zero the corresponding row of NewPS.*/
					MakeZeroRow(&NewPS, i);
				}
			}
			/*Remove duplicate columns*/
			for (i = 0; i < NumberOfColumns(NewPS) - 1; i++)
			{
				for (j = i + 1; j < NumberOfColumns(NewPS); j++)
				{
					/*Check the columns i and j.?*/
					for (k = 0; k < NumberOfRows(NewPS); k++)
					{
						if (GetMatrixEl(&NewPS, k, i) != GetMatrixEl(&NewPS, k, j))
						{
							break;
						}
					}
					/*Are they identical?*/
					if (k == NumberOfRows(NewPS))
					{
						/*If so, remove the jth column.*/
						RemoveColumns(&NewPS, j, 1, -1);
					}
				}
			}
		}
		
		/*Call asiph recursively on Dma and Dpa.*/
		recursive = asiph(Dma, Dpa, newList, newCount, PS ? &NewPS : 0, 0, 0);
		sx = recursive.MS;
		
		/*Free memory if necessary.*/
		if (recursive.NMS.type) DeallocateMatrix(&recursive.NMS);
	}
	/*Setup the data matrices appropriately. NMS should be the transpose copy
	of PS. If PS is empty, then setup NMS with zero rows and the same number
	of columns as there are places so that row additions will work.
	NS should be an empty matrix setup with zero rows and the same number of
	columns as there are places.
	Both matrices should be type-2 untransposed to optimize for row ops.*/	
	if (PS)
	{
		AllocateMatrixType(2, &data->MNS, NumberOfColumns(*PS), NumberOfRows(*PS));
		TransposeMatrix(PS);
		CopyMatrix(PS, &data->MNS);
		TransposeMatrix(PS);
	}
	else
	{
		AllocateMatrixType(2, &data->MNS, 0, NumberOfRows(*Dm));
	}
	AllocateMatrixType(2, &data->NS, 0, NumberOfRows(*Dm));
	FreeMemory(&memory);
	return sx;
}

/*******************************************************************************
This function checks the parameters to verify that the right ones are present,
sizes match, indices are within range, and so forth.*/
int CheckParams(matrix *Dm, matrix *Dp, int **list, int *listCount, matrix *PS, matrix **Dma, matrix **Dpa, MemoryManager *mem)
{
	int i, j;
	/*Dm and Dp are required and must be the same size.*/
	if (! (Dm && Dp))
	{
		merror(0, "ASIPH: A required parameter was a null pointer");
		return 0;
	}
	if (NumberOfRows(*Dm) != NumberOfRows(*Dp) || NumberOfColumns(*Dm) != NumberOfColumns(*Dp))
	{
		merror(0, "ASIPH: The input and output matrices are not the same size");
		return 0;
	}
	
	/*List is not required. If not present, it defaults to including the
	indices of each place. It is checked to be sure it includes no duplicate or
	spurious indices.*/
	if (! *list) *listCount = 0;
	if (! *listCount)
	{
		*listCount = NumberOfRows(*Dm);
		*list = mcalloc(mem, *listCount, sizeof(int));
		for (i = 1; i < *listCount; i++)
		{
			(*list)[i] = i;
		}
	}
	else
	{
		for (i = 0; i < *listCount; i++)
		{
			/*Make sure no entries are out-of-range.*/
			if ((*list)[i] < 0 || (*list)[i] >= NumberOfRows(*Dm))
			{
				FreeMemory(mem);
				merror(0, "ASIPH: One of the indices in the place list is out of range");
				return 0;
			}
			/*Make sure no entries are duplicate.*/
			for (j = i + 1; j < *listCount; j++)
			{
				if ((*list)[i] == (*list)[j])
				{
					FreeMemory(mem);
					merror(0, "ASIPH: The place list contains duplicate entries");
					return 0;
				}
			}
		}
	}
	/*PS is not required but if present must have the same number of rows as 
	Dm/Dp.*/
	if (PS)
	{
		if (NumberOfRows(*PS) != NumberOfRows(*Dm))
		{	
			FreeMemory(mem);
			merror(0, "ASIPH: PS does not have the same number of rows as the Petri net matrices");
			return 0;
		}
	}
	/*Dma and Dpa are not required. If present they must be the same size as Dp & 
	Dp. If not present they will be found via call to actn.*/
	if (*Dma && *Dpa)
	{
		if (NumberOfRows(**Dma) != NumberOfRows(**Dpa) || NumberOfColumns(**Dma) != NumberOfColumns(**Dpa))
		{
			FreeMemory(mem);
			merror(0, "ASIPH: Dma and/or Dpa is not the same size as Dm and Dp");
			return 0;
		}
	}
	else
	{
		actn_r result;
		result = actn(Dm, Dp, 0, 0, 0, 0, 0, 0);
		if (result.Dmra.type) DeallocateMatrix(&result.Dmra);
		if (result.Dpra.type) DeallocateMatrix(&result.Dpra);

		*Dma = mmalloc(mem, sizeof(matrix));
		**Dma = result.Dma;
		*Dpa = mmalloc(mem, sizeof(matrix));
		**Dpa = result.Dpa;
		ManageMatrix(mem, *Dma);
		ManageMatrix(mem, *Dpa);

		if (result.TA) free(result.TA);
	}
	return 1;
}
