#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include "../pnheaders/general.h"
#include "../pnheaders/matrix.h"
#include "../pnheaders/pns.h"
#include "matrixmath.h"
#include "extendedmatrix.h"
#include "spnbox.h"

typedef struct dpState
{
	unsigned dptype   : 2; /*Type of supervision: 0 = deadlock prevention, 1 =
	                       liveness enforcement, 2 = T-liveness enforcement.*/
	unsigned unique   : 1; /*1 indicates a unique nonzero active subnet this
	                       iteration?*/
	unsigned perm     : 1; /*1 indicates successful siphon control.*/
	unsigned anetfail : 1; /*1 if the active subnet does not contain all
	                       transitions in the net (DP and LE case) or all 
											   transitions of T (T-LE case).*/
	unsigned incon    : 1; /*1 if initial constraints exist.*/
	unsigned uncontro : 1; /*1 if uncontrollable and/or unobservable transitions
	                       exist.*/
	unsigned usfl     : 1; /*1 if redundant constraints have been eliminated.*/
	unsigned nolive   : 1; /*1 if no liveness-enforcing supervisor exists.*/
	unsigned noTlive  : 1; /*1 if no T-liveness-enforcing supervisor exists.*/	
} dpState;

typedef struct dpData
{
	dpState state;
	/*Various matrices*/
	matrix Dm, Dp, DI, Dcm, Dcp, L, L0, M;
	/*Constraint-related vectors*/
	int *B, *B0, *G;
	/*Place lists:*/
	int *opl;	/*Original places*/
	int *ipl; /*Independent places*/
	int *dpl; /*Dependent places*/
	int *apl; /*Active places*/
	int *upl; /*Useful places*/
	int *cpl; /*New control places added during a third-level sub-iteration.*/
	int *ExcludeT, *IncludeT; /*Transitions that must be included in or excluded
	from use in subnet computation.*/
	int *TA; /*Active transitions.*/
	int *Tuc, *Tuo; /*Uncontrollable and unobservable transitions.*/
	
	/*Counts of the list*/
	int opCount, ipCount, dpCount, apCount, upCount, cpCount, ExcludeCount, IncludeCount, TACount, TucCount, TuoCount;
	
	/*The size of the target Petri net.*/
	int targetRows, targetCols;

	/*The log file pointers.*/
	FILE *log, *dat;
	
	/*Data related to primary processing runs*/
	int firstpass, restart, run, terminate;
	char *status;
	
	/*Data used during the inner iterations.*/
	int iteration; //Iteration number
	int uci, upd, consis; //flags
	matrix TSIPH; //All siphons found during a run.
	matrix S; //New siphons found during a run.
} dpData;

/*A set of marking constraints. The return type of UpdateMarkingConstraints,
which concatenates the initial marking constraints with the ordinary marking
constraints.*/
typedef struct MarkingConstraints
{	
	matrix L;
	int *B;
} MarkingConstraints;

/*Functions used by dp, implemented in this file. See individual functions for 
more information.*/
static int CheckParams(dp_p *p);
static void Sort(int* array, int length);
static void AddToArray(int** array, int length, int item);
static void AddToIndexArray(int** array, int* length, int item);
static void UpdateDPData(dpData *data, char* description, void *newData);
static MarkingConstraints UpdateMarkingConstraints(dpData *d);
static dpData CreateDpData(dp_p *p);
static void InitDpData(dp_p *p, dpData* d);
static void PrintMessage(dpData* d, char* message);
static void PrintState(dpData *d);
static void asbnmessage(dpData *d);
static void EnforceMarkingConstraints(dpData* d);
static void UpdateActiveNet(dpData* d);
static int ComputeSiphons(dpData* d);
static void PrintConstraints(dpData *d, matrix *l, int b, int *siphon, int cpl, int wn);
static void BuildSiphonConstraint(dpData* d, matrix* l, int* b, matrix* a, int *siphon);
static void AddSiphonExclusions(dpData* d, int* siphon);
static admcon_p BuildAdmconParams(dpData* d, int* siphon);
static int DoConstraintsApplyToAll(dpData *d, supervis_r* supervisRes, matrix* a, int* siphon);
static void CleanupSiphonChanges(dpData* d);
static void EnforceSiphonConstraints(dpData *d);
static void EnforceSiphonConstraints(dpData *d);
static void ReduceLB(dpData *d);
static char* ListTransitions(int* list, int listCount);
static void PrintConstraintSet(matrix* L, int* B, FILE* log);
static void FinalPrint(dpData* d);
static void FreeDpParts(dpData* d);

dp_r dp(dp_p* params)
{
	dp_r result;
	int i, j, flag;
	memset(&result, 0, sizeof(dp_r));
	/*Check the parameters.*/
	if (! CheckParams(params))
	{
		return result;
	}
	
	/*Initialize the data structure.*/
	dpData d = CreateDpData(params);
	/*Initial message printing.*/
	if (is_verbose() >= VRB_DP)
	{
		switch(d.state.dptype)
		{
		case 0:
			printf("Producing a deadlock prevention supervisor...\n");
			break;
		case 1:
			printf("Producing a liveness enforcing supervisor...\n");
			break;
		case 2:
			printf("Producing a T-liveness enforcing supervisor...\n");
			break;
		}
	}
	
	/*Make sure the constraints are consistent: Get a set of complete constraints
	and then check it using chkcons.*/	
	MarkingConstraints cons = UpdateMarkingConstraints(&d);
	if (cons.L.type)
	{
		chkcons_r chkconsRes = chkcons(&cons.L, cons.B, 0, 0);
		/*If resCount is 0, it means the constraints are redundant.*/
		if (! chkconsRes.resCount)
		{
			if (d.log)
			{
				fprintf(d.log, "\n\nInconsistent set of constraints. Relax (L, B) and (L0, B0).");
			}
			if (is_verbose() > VRB_DP)
			{
				printf("\n\nInconsistent set of constraints. Relax (L, B) and (L0, B0).\n");
			}
			/*Mark for termination.*/
			d.terminate = 1;
		}
		DeallocateMatrix(&cons.L);
		free(cons.B);
		DeallocateChkcons(&chkconsRes);
	}		
	
	/*Main processing begins here.*/
	while ((d.restart || d.firstpass) && ! d.terminate)
	{
		d.restart = 0;
		if (++d.run > 1 && is_verbose() >= VRB_DP)
		{
			printf("\n ====================== RUN %d ========================\n",d.run);
		}
		
		/*Setup the dp data for this run.*/
		InitDpData(params, &d);
		
		/*If there are marking constraints, enforce them.*/
		EnforceMarkingConstraints(&d);
		
		/*Build the useful places list: include all places initially.*/
		if (d.upl) free(d.upl);
		d.upCount = NumberOfRows(d.Dm);
		d.upl = tcalloc(d.upCount, sizeof(int));
		for (i = 0; i < d.upCount; i++)
		{
			d.upl[i] = i;
		}
		
		/*More log/verbose output output.*/
		if (d.log)
		{
			if (d.firstpass)
			{
				fprintf(d.dat, "\nStarting with the Petri net:");
			}
			else
			{
				fprintf(d.dat, "\n ================= RUN NO. %d ================\n", d.run);
			}
			fflush(d.dat);
			PrintState(&d);
		}
		if (is_verbose() >= VRB_DP)
		{
			asbnmessage(&d);
		}
		
		/*Do an msplit.*/
		msplit_r msplitRes = msplit(&d.Dm, &d.Dp, 0, &d.M, &d.L0, &d.L, d.ipl, d.ipCount, d.dpl, d.dpCount, 0, 0);
		UpdateDPData(&d, "free matrix Dm Dp M L0 L ipl ipCount free TD free intarray", &msplitRes);
		
		/*Update the listing of the active subnet.*/
		UpdateActiveNet(&d);
		
		/*Get ready for the inner iteration loop. Initialize TSIPH, the listing
		of active siphons found this primary iteration, set the iteration number to
		1, and clear the uci flag.*/
		/*TSIPH should be initialized for column ops, as a siphon is recorded as
		a column.*/
		if (d.TSIPH.type) DeallocateMatrix(&d.TSIPH);
		AllocateMatrixType(2, &d.TSIPH, 0, NumberOfRows(d.Dm));
		TransposeMatrix(&d.TSIPH);
		d.iteration = 1;
		d.uci = 0;
		
		/*The inner iteration loop.*/
		while (! d.terminate)
		{
			/*Log and verbosity stuff.*/
			d.upd = 1;
			if (d.log)
			{
				if (d.firstpass) fprintf(d.dat, "\nIteration %d", d.iteration);
				else fprintf(d.dat, "\nRun %d, iteration %d", d.run, d.iteration);
				fflush(d.dat);
			}
			
			/*Do the actual siphon computation. This function returns the number
			of new siphons. If there were no new siphons, terminate the loop.*/
			if (! ComputeSiphons(&d)) break;
			
			/*Build constraints from the siphons and enforce them.*/
			EnforceSiphonConstraints(&d);
			
			/*Cleanup and finish applying the changes made to this point.*/
			CleanupSiphonChanges(&d);
			d.iteration++;			
		}
		d.firstpass = 0;
		d.terminate = ! d.restart;	
	}
	/*Bring the marking constraints down to target size.*/
	if (NumberOfColumns(d.L) - d.targetRows > 0) 
	{
		RemoveColumns(&d.L, d.targetRows, NumberOfColumns(d.L) - d.targetRows, -1);
	}
	if (NumberOfColumns(d.L0) - d.targetRows > 0) 
	{
		RemoveColumns(&d.L0, d.targetRows, NumberOfColumns(d.L0) - d.targetRows, -1);
	}
	
	/*If there are any active transitions left, remove redundant constraints.*/
	if (d.TACount)
	{
		ReduceLB(&d);
	}
	
	/*Build the result.*/
	result.Lf = d.L;
	result.Bf = d.B;
	result.L0f = d.L0;
	result.B0f = d.B0;
	result.how = d.status;
	
	/*Do log printing and closeout.*/
	if (d.log)
	{
		FinalPrint(&d);
		fclose(d.log);
		fclose(d.dat);
		if (is_verbose() >= VRB_DP)
		{
			printf("\nOutput written to dp.log and dp.dat\n");
		}
	}
	
	/*Free up other memory.*/
	FreeDpParts(&d);
	
	return result;
}

/*******************************************************************************
FreeDpParts frees up parts of the dpData structure that are not being used as 
part of the return value.*/
void FreeDpParts(dpData* d)
{
	int i;
	DeallocateMatrix(&d->Dm);
	DeallocateMatrix(&d->Dp);
	DeallocateMatrix(&d->DI);
	DeallocateMatrix(&d->Dcm);
	DeallocateMatrix(&d->Dcp);
	DeallocateMatrix(&d->M);
	free(d->G);
	free(d->opl);
	free(d->ipl);
	free(d->dpl);
	free(d->apl);
	free(d->upl);
	free(d->cpl);
	free(d->ExcludeT);
	free(d->IncludeT);
	free(d->TA);
	free(d->Tuc);
	free(d->Tuo);
	
	DeallocateMatrix(&d->TSIPH);
	/*S may not ever have been allocated if no new siphons were found.*/
	if (d->S.type) DeallocateMatrix(&d->S);
	memset(d, 0, sizeof(dpData));
}

/*******************************************************************************
FinalPrint prints the final data to the log and the data file.*/
void FinalPrint(dpData* d)
{
	/*Generate a list that is the intersection of T and TA.*/
	int *TinTA;
	int TinTACount = 0, i, j;
	if (d->TACount && d->IncludeCount)
	{
		TinTA = tcalloc(d->TACount < d->IncludeCount ? d->TACount : d->IncludeCount, sizeof(int));
		for (i = 0, j = 0; i < d->TACount && j < d->IncludeCount; i++)
		{
			while(d->IncludeT[j] < d->TA[i]) j++;
			if (d->IncludeT[j] == d->TA[i])
			{
				TinTA[TinTACount++] = d->TA[i];
			}
		}
	}
	/*flags*/
	int failed = ! strcmp(d->status, HOW_FAILED);
	int impossible = ! strcmp(d->status, HOW_IMPOSSIBLE);
	int ok = ! strcmp(d->status, HOW_OK);
	char *achieve, *supname, *tmp;;
	
	fprintf(d->log, "\n\n===============================================================================\n\n");
	if (d->state.dptype == 0)
	{
		if (ok)
		{
			fprintf(d->log, "The generated supervisor is ");
			if (d->state.perm && d->state.unique)
			{
				fprintf(d->log, "the least restrictive ");
				if (d->state.anetfail)
				{
					fprintf(d->log, "T-");
					achieve = "T-live";
				}
				else
				{
					achieve = "live";
				}
				fprintf(d->log, "liveness enforcing supervisor.");
				if (d->state.anetfail)
				{
					fprintf(d->log, "\nThe set T is:\nT = %s", ListTransitions(d->TA, d->TACount));
				}
			}
			else if (d->state.perm)
			{
				fprintf(d->log, "guaranteed to prevent total deadlock");
				achieve = "deadlock-free";
				if (d->state.anetfail)
				{
					fprintf(d->log, " in the subsystem containing the transitions\nT = %s", ListTransitions(d->TA, d->TACount));
					fprintf(d->log, "The supervisor may not enforce T-liveness.\n");
					fprintf(d->log, "The supervisor is guaranteed to be at least as permissive as the least\n");
					fprintf(d->log, "restrictive T-liveness enforcing supervisor.");
				}
				else
				{
					fprintf(d->log, ".\nThe supervisor may not enforce liveness.\n");
					fprintf(d->log, "The supervisor is guaranteed to be at least as permissive as the least\n");
					fprintf(d->log, "restrictive liveness enforcing supervisor");
				}					
			}
			else
			{
				fprintf(d->log, "guaranteed to prevent total deadlock");
				achieve = "deadlock-free";
				if (d->state.anetfail)
				{
					fprintf(d->log, " in the subsystem containing the transitions\nT = %s", ListTransitions(d->TA, d->TACount));
					fprintf(d->log, "The supervisor may not enforce T-liveness.\n");
				}
				else
				{
					fprintf(d->log, ".\nThe supervisor may not enforce liveness.\n");
				}
			}
		}
		else if (impossible)
		{
			fprintf(d->log, "No deadlock prevention supervisors exist.");
			if (d->state.incon && d->state.uncontro)
			{
				fprintf(d->log, "\nDeadlock prevention supervisors may exist for other initial constraints.");
			}
			else if (d->state.incon)
			{
				fprintf(d->log, "\nDeadlock prevention supervisors exist for other initial constraints.");
			}
		}
		else if (failed)
		{
			fprintf(d->log, "The procedure could not generate a deadlock prevention supervisor.");
		}
	}
	else
	{		
		if (d->state.dptype == 1)
		{
			supname = "liveness";
			tmp = "live";			
		}
		else
		{
			supname = "T-liveness";
			tmp = "T-live";
		}
		if (ok)
		{
			fprintf(d->log, "The generated supervisor is ");
			if (d->state.unique && d->state.perm)
			{
				fprintf(d->log, "the least restrictive ");
				if (d->state.anetfail)
				{
					fprintf(d->log, "Tx-liveness");
					achieve = "Tx-live";
				}
				else
				{
					fprintf(d->log, "%s", supname);
					achieve = tmp;
				}
				fprintf(d->log, " enforcing supervisor");
				if (d->state.anetfail)
				{
					fprintf(d->log, "The set Tx is:\nTx = %s", ListTransitions(d->TA, d->TACount));
				}
			}
			else
			{
				fprintf(d->log, "guaranteed to enforce");
				if (d->state.anetfail)
				{
					fprintf(d->log, "Tx-liveness in the subsystem containing the transitions\n");
					fprintf(d->log, "Tx = %s", ListTransitions(TinTA, TinTACount));
					achieve = "Tx-live";
				}
				else
				{
					fprintf(d->log, "%s, ", supname);
					achieve = tmp;
				}
				fprintf(d->log, "and may not be the least restrictive");				
			}
			if (d->state.perm && ! d->state.unique)
			{
				fprintf(d->log, "\nThe supervisor is the least restrictive Ta-liveness ");
				fprintf(d->log, "enforcing supervisor for\nTa = %s", ListTransitions(d->TA, d->TACount));				
			}
			if (d->state.nolive)
			{
				fprintf(d->log, "\nNo liveness enforcing supervisor exists.");
			}
			if (d->state.noTlive)
			{
				fprintf(d->log, "\nNo T-livenss enforcing supervisor exists.");
			}
		}
		else if (failed)
		{
			fprintf(d->log, "The procedure could not generate a liveness or T-liveness enforcing supervisor.");		
		}
		else if (impossible)
		{
			fprintf(d->log, "No liveness/Tx-liveness enforcing supervisors exist, for any nonempty Tx.");
			fprintf(d->log, "\nEven deadlock prevention is impossible.");
		}
	}
	if (! ok) return;
	switch(d->state.dptype)
	{
	case 0:
		supname = "deadlock-free";
		break;
	case 1:
		supname = "live";
		break;
	case 2:
		supname = "T-live";
		break;
	}

	if (! (NumberOfRows(d->L) && NumberOfRows(d->L0)))
	{
		fprintf(d->log, "\nNo constraints are needed: the Petri net is %s for all initial\nmarkings.", supname);
	}
	else
	{
		if (NumberOfRows(d->L))
		{
			fprintf(d->log, "\nThe supervisor is defined by the following constraints:\n");
			PrintConstraintSet(&d->L, d->B, d->log);
		}
		if (NumberOfRows(d->L0))
		{
			if (NumberOfRows(d->L))
			{
				fprintf(d->log, "\n\nIn addition to the inequalities above, the initial marking \"i\" must satisfy:\n");
			}
			else
			{
				fprintf(d->log, "\n\nThe target Petri net is %s for all intiial markings \"i\" satisfying\n", achieve);
			}
			PrintConstraintSet(&d->L0, d->B0, d->log);
		}
	}
	if (d->state.usfl)
	{
		fprintf(d->log, "\n\n(The redundant constraints have been eliminated.)");
	}
	fflush(d->log);
	free(TinTA);
}

/*******************************************************************************
PrintConstraintSet prints the constraints represented by the matrix L and the
integer array B in pretty format to the file given.*/
void PrintConstraintSet(matrix* L, int* B, FILE* log)
{
	int i;
	matrix lrow;
	AllocateMatrix(&lrow, 1, NumberOfColumns(*L));
	for (i = 0; i < NumberOfRows(*L); i++)
	{
		CopyBlock(1, NumberOfColumns(*L), L, i, 0, &lrow, 0, 0);
		fprintf(log, "\n%s >= %d", avpr(&lrow, 0, DISPLAY_MATRIX | DISPLAY_OFFSET, "m", 0), B[i]);		
	}
	fflush(log);
	DeallocateMatrix(&lrow);
}

/*******************************************************************************
ListTransitions generates a pretty form of a transition index list.*/
char* ListTransitions(int* list, int listCount)
{
	static char* buffer = 0;
	static int bufferlength = 0;
	char element[32];
	int i = 0, charsPerElement = 3, newbufferlength;
	/*Compute the necessary buffer length.*/
	/*Take the log (to find the digit count) of the highest index.*/
	if (listCount) i = list[listCount - 1] + 1;
	while (i)
	{
		charsPerElement++;
		i /= 10;
	}
	if ((newbufferlength = charsPerElement * listCount + 4) > bufferlength)
	{
		free(buffer);
		bufferlength = newbufferlength;
		buffer = tcalloc(bufferlength, sizeof(char));
	}
	
	/*Fill the string.*/
	strcpy(buffer, "{");
	for (i = 0; i < listCount - 1; i++)
	{
		sprintf(element, "t%d, ", list[i] + 1);
		strcat(buffer, element);
	}
	if (listCount)
	{
		sprintf(element, "t%d", list[listCount - 1]);
		strcat(buffer, element);
	}
	strcat(buffer, "}\n");
	return buffer;
}

/*******************************************************************************
ReduceLB removes redundant constraints.*/
void ReduceLB(dpData *d)
{
	MarkingConstraints cons = UpdateMarkingConstraints(d);
	//If there are no constraints, return immediately.
	if (! cons.L.type) return;
	
	reduce_r rcons = reduce(&cons.L, cons.B);
	free(cons.B);
	DeallocateMatrix(&cons.L);
	int i, j, k, count = 0;
	/*We need to know which indices in the useful-places list returned by reduce
	correspond to L rows and which correspond to L0 rows. Find the break.*/
	Sort(rcons.indf, NumberOfRows(rcons.Lf));
	for (count = 0; count < NumberOfRows(rcons.Lf) && rcons.indf[count] < NumberOfRows(d->L); count++);
	/*The new marking constraint matrix will be all rows of the original that
	are present in the list of useful rows returned by reduce. Remove all other
	rows.*/	
	j = NumberOfRows(d->L) - 1;
	for (i = count - 1; i >= -1; i--)
	{
		while((i >= 0) ? j > rcons.indf[i] : j >= 0)
		{
			/*Remove the appropriate row from L.*/
			RemoveRows(&d->L, j, 1, -1);
			/*Remove the appropriate entry from B.*/
			for (k = j; k < NumberOfRows(d->L); k++)
			{
				d->B[k] = d->B[k + 1];
			}			
			j--;			
		}
		j--;
	}
	/*The new initial marking constraint will be all rows of the original that
	are present in the list of useful rows returned by reduce. Remove all other
	rows.*/
	j = NumberOfRows(d->L0) - 1;
	for (i = NumberOfRows(rcons.Lf) - 1; i >= count - 1; i--)
	{
		while((i >= count) ? j > rcons.indf[i] - NumberOfRows(d->L) : j >= 0)
		{
			RemoveRows(&d->L0, j, 1, -1);
			for (k = j; k < NumberOfRows(d->L0); k++)
			{
				d->B0[k] = d->B0[k + 1];
			}
			j--;
		}
		j--;
	}
	/*Set the useful flag.*/
	d->state.usfl = (NumberOfRows(d->L) + NumberOfRows(d->L0) - NumberOfRows(rcons.Lf)) ? 1 : 0;
	/*Free up memory used by the reduce-call result.*/
	DeallocateReduce(&rcons);	
}

/*******************************************************************************
EnforceSiphonConstraints creates and enforces new constraints based on the new
siphon lists.*/
void EnforceSiphonConstraints(dpData *d)
{
	/*We will need several vectors with elements for each independent place. 
	Allocate space now.*/	
	matrix l, l2, a;
	AllocateMatrix(&a, 0, 0);
	int siphon, i, j, k, b, x, flag, siphons;
	int* currentSiphon;
	/*Zero the (new) control place list.*/
	free(d->cpl);
	d->cpl = 0;
	d->cpCount = 0;
	siphons = NumberOfColumns(d->S);
	currentSiphon = tcalloc(NumberOfRows(d->S), sizeof(int));
	
	/*Iterate through each new siphon.*/
	for (siphon = 0; siphon < siphons; siphon++)
	{
		/*Get a copy of the current siphon as an integer array.*/
		for (i = 0; i < NumberOfRows(d->S); i++)
		{
			currentSiphon[i] = GetMatrixEl(&d->S, i, siphon) ? 1 : 0;			
		}
		AllocateMatrixType(2, &l, 1, d->ipCount);
		DeallocateMatrix(&a);
		AllocateMatrixType(1, &a, 1, NumberOfRows(d->Dm));
		
		char message[64];
		/*Message printing.*/
		if (is_verbose() >= VRB_DP)
		{
			sprintf(message, "\nWorking on siphon %d of %d...", siphon + 1, siphons);
			PrintMessage(d, message);
		}
		
		/*Build a constraint based on this siphon.*/
		BuildSiphonConstraint(d, &l, &b, 0, currentSiphon);
		
		/*Check to see if the constraint is redundant.*/
		MarkingConstraints cons = UpdateMarkingConstraints(d);
		chkcons_r chkconsRes;
		
		if (cons.L.type)
		{
			chkconsRes = chkcons(&cons.L, cons.B, &l, b);
		}
		else
		{
			/*If there are no constraints, then the proposed constraint is 
			nonredundant. Note this to be the case.*/
			memset(&chkconsRes, 0, sizeof(chkcons_r));
			chkconsRes.resCount = 1;
		}
		/*If res is not empty, the constraint is not redundant.*/
		if (chkconsRes.resCount)
		{
			DeallocateChkcons(&chkconsRes);
			/*Build the admcon parameters.*/
			admcon_p admconP = BuildAdmconParams(d, currentSiphon);
			/*Make the admcon call.*/
			admcon_r admconRes = admcon(&admconP);
			/*Free the memory that was allocated for the admcon parameters.*/
			free(admconP.ntr);			
			
			/*If admcon found an admissible constraint...*/
			if (admconRes.how)
			{
				/*Turn it into a constraint matrix and enforce it.*/
				for (i = 0; i < admconRes.lCount; i++)
				{
					SetMatrixEl(&a, 0, i, admconRes.l[i]);
				}
				supervis_r supervisRes = supervis(&d->Dp, &d->Dm, &a);
				
				/*Build l and b as before, but this time using a as the source.*/
				BuildSiphonConstraint(d, &l, &b, &a, 0);
				/*Check the constraints if applicable.*/
				if (d->consis)
				{
					/*We are checking constraints for consistency, not redundancy - add
					the proposed new constraint to the set of constraints. Use a copy
					so that the memory can be absorbed into L in an optimized operation.*/
					AddToArray(&cons.B, NumberOfRows(cons.L), b);
					AllocateMatrixType(2, &l2, 1, NumberOfColumns(l));
					CopyMatrix(&l, &l2);
					InsertRows(&cons.L, &l2, -1, -1);					
					chkconsRes = chkcons(&cons.L, cons.B, 0, 0);
				}
				else
				{
					memset(&chkconsRes, 0, sizeof(chkcons_r));
					chkconsRes.resCount = 1;
				}
				/*We want to find out if there exist any elements of the vector a such
				that their zero/nonzero state is the opposite of that of the current
				siphon.*/
				x = 0;
				for (i = 0; i < NumberOfRows(d->S); i++)
				{
					if ((! GetMatrixEl(&d->S, i, siphon)) != (! GetMatrixEl(&a, 0, i)))
					{
						/*If such a place is encountered, set the generic variable x to 2 for
						later use as a function parameter.*/
						x = 2;
						break;
					}
				}
				/*If no such place existed, clear perm.*/
				if (! x) d->state.perm = 0;
				
				/*If the constraints were consistent, get ready to add them either to
				L, B or L0, B0.*/
				if (chkconsRes.resCount)
				{
					DeallocateChkcons(&chkconsRes);
					/*Find out whether these constraints should apply to the marking at
					all times or just to the initial marking.*/
					if (DoConstraintsApplyToAll(d, &supervisRes, &a, currentSiphon))
					{
						/*Print constraint data.*/
						if (d->log) PrintConstraints(d, &l, b, currentSiphon, NumberOfRows(d->Dm), x);
					/*Add the new place index to the dependent place list, the useful
						places list, and the new control places list.*/
						AddToIndexArray(&d->dpl, &d->dpCount, NumberOfRows(supervisRes.Dfm) - 1);
						AddToIndexArray(&d->upl, &d->upCount, NumberOfRows(supervisRes.Dfm) - 1);
						AddToIndexArray(&d->cpl, &d->cpCount, NumberOfRows(supervisRes.Dfm) - 1);
						/*Dm and Dp should be set to the supervised net Dfm and Dfp.*/
						UpdateDPData(d, "Dm Dp", &supervisRes);						
						memset(&supervisRes, 0, sizeof(supervis_r));
						/*Add to the constraints. Make a copy of l first so that we can do
						optimized row-adds and not do too many allocations.*/						
						AllocateMatrixType(2, &l2, 1, NumberOfColumns(l));						
						CopyMatrix(&l, &l2);
						AddToArray(&d->B, NumberOfRows(d->L), b);
						AddToArray(&d->G, NumberOfRows(d->L), b);
						InsertRows(&d->L, &l, -1, -1);
						InsertRows(&d->M, &l2, -1, -1);
						/*Add a null row to the end of S and allocate enough space for the
						new siphon array.*/
						InsertNullRows(&d->S, -1, 1, -1);
						int *newCSiphon = tcalloc(NumberOfRows(d->S), sizeof(int));
						memcpy(newCSiphon, currentSiphon, sizeof(int) * (NumberOfRows(d->S) - 1));
						free(currentSiphon);
						currentSiphon = newCSiphon;						
					}
					else
					{						
						/*Add the built constraint to the initial marking constraints.*/
						if (d->log) PrintConstraints(d, &l, b, currentSiphon, 0, x);
						AddToArray(&d->B0, NumberOfRows(d->L0), b);
						InsertRows(&d->L0, &l, -1, -1);
					}
				}
				else
				{
					/*If the constraints were redundant, just do a constraint printing
				  with no new constraints. Also deallocate l.*/
					DeallocateMatrix(&l);
					if (d->log) PrintConstraints(d, 0, 0, currentSiphon, -1, x);
					d->upd = 0;
					/*Add some indices to X: those transitions such that they have nonzero
					output arcs from at least one of the places in the current siphon.*/
					AddSiphonExclusions(d, currentSiphon);
					d->state.perm = d->state.perm > 1 ? 2 : 0;
					d->restart = 1;
					d->terminate = 1;
				}
				/*Deallocate the supervis result if it still needs to be deallocated.*/
				DeallocateSupervis(&supervisRes);					
			}
			else
			{
				/*If admcon could not find admissible constraints, add to the excluded
				transition list in the same way we did earlier.*/
				AddSiphonExclusions(d, currentSiphon);
				d->uci = 1;
				d->state.perm = 0;
				d->upd = 0;
				if (d->log) PrintConstraints(d, 0, 0, currentSiphon, -2, 1);
			}
			DeallocateAdmcon(&admconRes);
		}
		else
		{
			if (d->log) PrintConstraints(d, 0, 0, currentSiphon, -1, 0);
		}
		/*If l has not been absorbed into something else, deallocate it here.*/
		if (l.type) DeallocateMatrix(&l);
		if (NumberOfRows(cons.L))
		{
			DeallocateMatrix(&cons.L);
			free(cons.B);
		}
	}		
	DeallocateMatrix(&a);
	free(currentSiphon);
}

/*******************************************************************************
CleanupSiphonChanges is a subroutine of EnforceSiphonConstraints. It performs a
transition split and an eac-net conversion on the net, then updates the active
subnet listing.*/
void CleanupSiphonChanges(dpData* d)
{
	int i, j, k;
	
	/*Transition split.*/
	msplit_r msplitRes = msplit(&d->Dm, &d->Dp, 0, &d->M, &d->L0, &d->L, d->ipl, d->ipCount, d->dpl, d->dpCount, 0, 0);
	UpdateDPData(d, "free matrix Dm Dp M L0 L ipl ipCount free TD free intarray", &msplitRes);
	
	/*If we are doing any kind of liveness enforcement, convert the net to an
	eac-net.*/
	if (d->state.dptype)
	{
		/*We need a flag array for the new control places.*/
		int *ncp = tcalloc(NumberOfRows(d->Dm), sizeof(int));
		for (i = 0; i < d->cpCount; i++)
		{
			ncp[d->cpl[i]] = 1;
		}
		pn2acpn_r pn2acpnRes = pn2acpn(&d->Dm, &d->Dp, ncp, &d->M, &d->L0, &d->L, d->ipl, d->ipCount, d->dpl, d->dpCount);
		UpdateDPData(d, "free matrix Dm Dp M L0 L ipl ipCount", &pn2acpnRes);
		
		free(ncp);
	}
	
	tactn_r tactnRes;
	actn_r actnRes;
	/*Update the active subnet.*/
	if (is_verbose() >= VRB_DP)
	{
		PrintMessage(d, "Updating the active subnet...");
	}
	if (d->state.dptype)
	{
		if (d->upd)
		{
			tactnRes = tactn(&d->Dm, &d->Dp, d->IncludeT, d->IncludeCount, d->ExcludeT, d->ExcludeCount, &d->Dcm, &d->Dcp, 0);
			UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &tactnRes);
		}
		else
		{
			tactnRes = tactn(&d->Dm, &d->Dp, d->IncludeT, d->IncludeCount, d->ExcludeT, d->ExcludeCount, 0, 0, 1);
			d->state.unique = tactnRes.unique;
			UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &tactnRes);			
		}
	}
	else
	{
		if (d->upd)
		{
			actnRes = actn(&d->Dm, &d->Dp, d->ExcludeT, d->ExcludeCount, &d->Dcm, &d->Dcp, 1, 0);
			UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &actnRes);
		}
		else
		{
			actnRes = actn(&d->Dm, &d->Dp, d->ExcludeT, d->ExcludeCount, &d->Dcm, &d->Dcp, 0, 1);
			d->state.unique = actnRes.unique;
			UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &actnRes);
		}
	}
	/*Update the exclusion list. It should include all transitions which are
	inactive. We can build this list from the active transition list already
	computed.*/
	free(d->ExcludeT);
	d->ExcludeCount = NumberOfColumns(d->Dm) - d->TACount;
	if (d->ExcludeCount)
	{
		d->ExcludeT = tcalloc(d->ExcludeCount, sizeof(int));
		j = 0;
		k = 0;
		for (i = 0; i < d->TACount; i++)
		{
			while (j < d->TA[i])
			{
				d->ExcludeT[k++] = j++;
			}
			j++;
		}
		while (k < d->ExcludeCount)
		{
			d->ExcludeT[k++] = j++;
		}
	}
	else
	{
		d->ExcludeT = 0;
	}
	/*If there are no active transitions, set failure and termination flags.*/
	if (! d->TACount)
	{
		d->terminate = 1;
		d->restart = 0;
		d->status = HOW_FAILED;
		if (d->firstpass && (! d->uci) && ((! d->state.dptype) || d->state.unique))
		{
			d->status = HOW_IMPOSSIBLE;
		}
	}
	/*Update the useful places list. It should be all independent places with
	indices valid for the target net and all dependent places.*/
	if (d->upd)
	{
		/*Find the number of independent places to use.*/		
		for (i = 0; i < d->ipCount; i++)
		{
			if (d->ipl[i] >= d->targetRows) break;
		}		
		/*i now contains the number of independent places that we want to use.
		Update the total count and allocate memory.*/
		free(d->upl);
		d->upCount = i + d->dpCount;
		d->upl = tcalloc(d->upCount, sizeof(int));
		/*Fill in the independent places.*/
		for (j = 0; j < i; j++)
		{
			d->upl[j] = d->ipl[j];
		}
		/*Fill in the dependent places*/
		for (j = 0; j < d->dpCount; j++)
		{
			d->upl[i + j] = d->dpl[j];
		}
		/*Sort the list.*/
		if (d->upCount) Sort(d->upl, d->upCount);
	}
}

/*******************************************************************************
DoConstraintsApplyToAll is a subroutine of EnforceSiphonConstraints. It 
determines if a set of constraints generated in an admcon call should be added
as constraints on marking at any time or only on initial marking. It returns 
nonzero to indicate the former and zero to indicate the latter.*/
int DoConstraintsApplyToAll(dpData *d, supervis_r* supervisRes, matrix* a, int *siphon)
{
	int flag = 0, i, j;
	/*Iterate through transitions*/
	for (i = 0; i < NumberOfColumns(supervisRes->Dfm); i++)
	{
		/*Is there an output arc from the new place and no input arc to the
		new place?*/
		if (GetMatrixEl(&supervisRes->Dfm, NumberOfRows(supervisRes->Dfm) - 1, i) && ! GetMatrixEl(&supervisRes->Dfp, NumberOfRows(supervisRes->Dfp) - 1, i))
		{
			/*Check to make sure there are no input arcs to any transitions in
			the siphon.*/						
			for (j = 0; j < NumberOfRows(d->S); j++)
			{
				if (GetMatrixEl(&supervisRes->Dfp, j, i) && siphon[j]) break;
			}
			/*If there were none, set the flag and break out of the loop.*/
			if (j == NumberOfRows(d->S))
			{
				flag = 1;
				break;
			}
		}
	}
	/*If the flag has no yet been set, it may be set if there are any
	elements of a such that they are not equal to the corresponding
	elements in the current siphon.*/
	if (! flag)
	{
		for (i = 0; i < NumberOfColumns(*a); i++)
		{
			if (GetMatrixEl(a, 0, i) != siphon[i])
			{
				flag = 1;
				break;
			}
		}
	}
	return flag;
}

/*******************************************************************************
BuildAdmconParams is a subroutine of EnforceSiphonConstraints. It fills an
admcon_p structure with the parameters necessary to make the admcon call.*/
admcon_p BuildAdmconParams(dpData* d, int* siphon)
{
	/*We'll need an updated copy of the active place list. Build it now.*/
	int *newApl = tcalloc(NumberOfRows(d->Dcm), sizeof(int));
	int newApCount = 0, i, j;
	for (i = 0; i < NumberOfRows(d->Dcm); i++)
	{
		for (j = 0; j < NumberOfColumns(d->Dcm); j++)
		{
			if (GetMatrixEl(&d->Dcm, i, j) || GetMatrixEl(&d->Dcp, i, j))
			{
				break;
			}
		}
		if (j != NumberOfColumns(d->Dcm)) newApl[newApCount++] = i;			
	}
	free(d->apl);
	d->apl = newApl;
	d->apCount = newApCount;
	
	/*Transform the siphon row to an admissible constraint using admcon.
	Start by building the admcon call parameters:*/
	admcon_p admconP;
	/*Most of the parameters come straight from the problem data.*/
	admconP.siphon = siphon;
	admconP.Dm = d->Dm;
	admconP.Dp = d->Dp;
	admconP.DI = d->DI;
	admconP.M  = d->M;
	admconP.ipl = d->ipl;
	admconP.dpl = d->dpl;
	admconP.opl = d->opl;
	admconP.apl = d->apl;
	admconP.Tuc = d->Tuc;
	admconP.Tuo = d->Tuo;
	admconP.ipCount = d->ipCount;				
	admconP.dpCount = d->dpCount;
	admconP.opCount = d->opCount;
	admconP.apCount = d->apCount;
	admconP.TucCount = d->TucCount;
	admconP.TuoCount = d->TuoCount;				
	/*For ntr, we build a list of all transitions not part of the original
	matrix.*/
	if (NumberOfColumns(d->Dm) - d->targetCols)
	{
		admconP.ntCount = NumberOfColumns(d->Dm) - d->targetCols;
		admconP.ntr = tcalloc(admconP.ntCount, sizeof(int));
		for (i = d->targetCols; i < NumberOfColumns(d->Dm); i++)
		{
			admconP.ntr[i - d->targetCols] = i;
		}
	}
	else
	{
		admconP.ntCount = 0;
		admconP.ntr = 0;
	}
	/*We know these parameters are well formed. Skip the (time-consuming)
	parameter check.*/
	admconP.skipParamCheck = 1;
	return admconP;
}

/*******************************************************************************
AddSiphonExclusions is called as a subroutine of EnforceSiphonConstraints. It
is used under certain circumstances when a siphon does not result in an
admissible constraint set. It adds to the excluded transition list those
transitions such that they have nonzero output arcs from at least one of the
places in the current siphon.*/
void AddSiphonExclusions(dpData* d, int* siphon)
{
	/*Build a flag array to indicate the transitions we'll have to add.*/
	int *ExcludeFlag = tcalloc(NumberOfColumns(d->Dm), sizeof(int));
	int i, j, k;
	for (i = 0; i < d->ExcludeCount; i++)
	{
		ExcludeFlag[d->ExcludeT[i]] = 1;
	}
	/*Find the new transitions.*/
	k = d->ExcludeCount;
	for (j = 0; j < NumberOfColumns(d->Dm); j++)
	{
		for (i = 0; i < NumberOfRows(d->S); i++)
		{
			if (GetMatrixEl(&d->Dm, i, j) && siphon[i] && ! ExcludeFlag[j])
			{
				ExcludeFlag[j] = 1;
				k++;
			}
		}
	}
	/*Now convert the flag array to an index list. Use the same memory;
	this is a waste of memory but avoids an allocation.*/
	free(d->ExcludeT);
	d->ExcludeCount = k;
	d->ExcludeT = ExcludeFlag;
	k = 0;
	for (i = 0; i < NumberOfColumns(d->Dm); i++)
	{
		if (d->ExcludeT[i]) d->ExcludeT[k++] = i;
	}
}

/*******************************************************************************
BuildSiphonConstraint is called as a subroutine of EnforceSiphonConstraints. It 
builds a single constraint (a row vector and a constant) based on either the row
vector in matrix form in a or the row vector in array form in siphon. It will
use a if it is present (pointer is non-null).*/
void BuildSiphonConstraint(dpData *d, matrix* l, int* b, matrix* a, int *siphon)
{	
	int i, j;
	/*Build a constraint row (l and b) ignoring M.*/
	for (i = 0; i < d->ipCount; i++)
	{
		SetMatrixEl(l, 0, i, a ? GetMatrixEl(a, 0, d->ipl[i]) : siphon[d->ipl[i]]);		
	}
	*b = 1;
	/*Take into account M if it exists.*/
	if (d->M.type && NumberOfRows(d->M) && NumberOfColumns(d->M))
	{
		/*Add to each element of l the matrix product of the rows of the siphon 
		corresponding to dependent places and the transformation matrix M. The 
		siphon matrix is all 1s and 0s, so we can skip actual multiplication and 
		just do a test for zero.*/
		for (i = 0; i < NumberOfColumns(d->M); i++)
		{			
			for (j = 0; j < d->dpCount; j++)
			{
				if (a ? GetMatrixEl(a, 0, d->dpl[j]) : siphon[d->dpl[j]])
				{
					SetMatrixEl(l, 0, i, GetMatrixEl(l, 0, i) + GetMatrixEl(&d->M, j, i));
				}
			}
		}
		/*Add to b the matrix product of the column vector G (transposed) and those
		elements of the current siphon which correspond to dependent places. Again,
		we don't have to actually multiply but merely test for nonzero.*/
		for (i = 0; i < d->dpCount; i++)
		{
			if (a ? GetMatrixEl(a, 0, d->dpl[i]) : siphon[d->dpl[i]]) *b += d->G[i];
		}
	}	
}
/*******************************************************************************
This functions prints constraints to the log files during the building and
enforcement of constraints based on siphons.*/
void PrintConstraints(dpData *d, matrix *l, int b, int *siphon, int cpl, int wn)
{
	int i;
	/*Print the siphon*/
	fprintf(d->log, "\nSiphon %s", fvpr(siphon, NumberOfRows(d->S), DISPLAY_INTS, 2));	
	switch(cpl)
	{
	case 0:
		/*Display.*/
		fprintf(d->log, "\nNo control places needed. Constraints below added to (L0, B0):");
		fprintf(d->log, "\n%s >= %d", avpr(l, NumberOfColumns(*l), DISPLAY_MATRIX | DISPLAY_OFFSET, 0, 0), b);
		break;
	case -1:
		fprintf(d->log, "\nNo constraints added.");
		break;
	case -2:
		fprintf(d->log, "\nInadmissible constraint.");
		break;
	default:
		fprintf(d->log, "\nPlace %d added to control the siphon and enforce:", cpl);
		fprintf(d->log, "\n%s >= %d", avpr(l, NumberOfColumns(*l), DISPLAY_MATRIX | DISPLAY_OFFSET, 0, 0), b);
		break;
	}
	if (wn == 1) fprintf(d->log, "\nSiphon control failure occurred.");
	else if (wn == 2) fprintf(d->log, "\nSuboptimal admissibility transformation.");
	fprintf(d->log, "\n");	
	fflush(d->log);
}

/*******************************************************************************
ComputeSiphons finds the minimal active siphons for an iteration. It returns the
number of new siphons.*/
int ComputeSiphons(dpData* d)
{
	/*A status message.*/
	if (is_verbose() > VRB_DP)
	{
		PrintMessage(d, "Computing the minimal active siphons...");
	}

	/*We'll needs lists of active transitions and places. TA will have just been
	updated before this function call. Update apl. To avoid having to iterate
	through the whole matrix twice, preallocate newApl to its maximum possible
	size.*/
	int *newApl = tcalloc(NumberOfRows(d->Dcm), sizeof(int));
	int i, j;
	d->apCount = 0;
	for (i = 0; i < NumberOfRows(d->Dcm); i++)
	{
		for (j = 0; j < NumberOfColumns(d->Dcm); j++)
		{
			if (GetMatrixEl(&d->Dcm, i, j) || GetMatrixEl(&d->Dcp, i, j))
			{
				newApl[d->apCount++] = i;
				break;
			}
		}
	}
	free(d->apl);
	d->apl = newApl;	
	
	/*Build a list of places that will be considered for siphon computation. This
	is any place such that it is present in both the useful and the active place
	lists.*/
	int* possibleSiphons = tcalloc(d->apCount < d->upCount ? d->apCount : d->upCount, sizeof(int));
	int possibleSiphonCount = 0;
	j = 0;
	for (i = 0; i < d->apCount && j < d->upCount; i++)
	{
		while (j < d->upCount && d->upl[j] < d->apl[i]) j++;
		if (j < d->upCount && d->upl[j] == d->apl[i])
		{
			possibleSiphons[possibleSiphonCount++] = d->apl[i];
		}
	}
	
	/*Make sure TSIPH is the right size.*/
	if (NumberOfRows(d->TSIPH) < NumberOfRows(d->Dm))
	{
		InsertNullRows(&d->TSIPH, -1, NumberOfRows(d->Dm) - NumberOfRows(d->TSIPH), -1);
	}
	/*Remove any columns from TSIPH such that they have no nonzero entries in rows
	corresponding to active places.*/
	j = 0;
	while (j < NumberOfColumns(d->TSIPH))
	{
		for (i = 0; i < d->apCount; i++)
		{
			if (GetMatrixEl(&d->TSIPH, d->apl[i], j)) break;
		}
		if (i == d->apCount)
		{
			RemoveColumns(&d->TSIPH, j, 1, -1);
		}
		else
		{
			j++;
		}
	}
	
	/*Minimal active siphon computation.*/
	asiph_r asiphRes = asiph(&d->Dm, &d->Dp, possibleSiphons, possibleSiphonCount, &d->TSIPH, &d->Dcm, &d->Dcp);
	UpdateDPData(d, "TSIPH S", &asiphRes);
	
	/*Remove from the new siphons list those siphons which contain no useful
	places.*/
	free(possibleSiphons);
	j = 0;
	while (j < NumberOfColumns(d->S))
	{
		for (i = 0; i < d->upCount; i++)
		{
			if (GetMatrixEl(&d->S, d->upl[i], j)) break;
		}
		if (d->upCount == i)
		{
			RemoveColumns(&d->S, j, 1, -1);
		}
		else
		{
			j++;
		}
	}
	
	/*Get the new siphon count.*/
	int newSiphons = NumberOfColumns(d->S);

	/*Make sure that ll the new siphons are in fact siphons.*/
	issiph_r issiphRes = issiph(&d->Dm, &d->Dp, &d->S, 1);
	
	if (issiphRes.flx)
	{
		DeallocateIssiph(&issiphRes);
		merror(0, "Siphon Computation Error");
	}
	
	DeallocateIssiph(&issiphRes);	
	
	/*Do printing.*/
	if (d->log)
	{
		if (d->firstpass)
		{
			fprintf(d->log, "\nIteration %d - %d new minimum siphons", d->iteration, newSiphons);
		}
		else
		{
			fprintf(d->log, "\nRun %d, iteration %d - %d new minimum siphons", d->run, d->iteration, newSiphons);
		}
		PrintState(d);
		fflush(d->log);
	}
	return newSiphons;
}

/*******************************************************************************
UpdateActiveNet updates the active subnet as appropriate for the type of
enforcement we are doing.*/
void UpdateActiveNet(dpData* d)
{
	int i, j, flag;
	if (d->state.dptype)
	{
		/*If we get here we are doing (T-)liveness enforcement.
		Transform the net to an EAC net.*/
		pn2acpn_r pn2acpnRes = pn2acpn(&d->Dm, &d->Dp, 0, &d->M, &d->L0, &d->L, d->ipl, d->ipCount, d->dpl, d->dpCount);
		UpdateDPData(d, "free matrix Dm Dp M L0 L ipl ipCount", &pn2acpnRes);
		
		/*Find the t-minimal active subnet.*/
		tactn_r tactnRes = tactn(&d->Dm, &d->Dp, d->IncludeT, d->IncludeCount, d->ExcludeT, d->ExcludeCount, 0, 0, 1);
		d->state.unique = tactnRes.unique;
		UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &tactnRes);
		
		/*Check to make sure that all required transitions are present in the
		subnet. We assume that TA and IncludeT entries have been properly sorted,
		from lowest to highest. Let i be the index in IncludeT and j be the index
		in TA.*/
		j = 0;
		flag = ! d->TACount;
		for (i = 0; i < d->IncludeCount && ! flag; i++)
		{
			while(j < d->TACount - 1 && d->TA[j] < d->IncludeT[i]) j++;
			if (d->TA[j] != d->IncludeT[i])
			{
				flag = 1;
				break;
			}
		}
		if (flag)
		{
			d->state.anetfail = 1;
			if (d->firstpass)
			{
				if (d->state.dptype == 2) d->state.noTlive = 1;
				else d->state.nolive = 1;
			}
		}			
	}
	else
	{
		/*If we get here, we are doing deadlock prevention, not (T-)liveness
		enforcement. Begin with active subnet finding.*/
		actn_r actnRes = actn(&d->Dm, &d->Dp, d->ExcludeT, d->ExcludeCount, 0, 0, 0, 1);
		d->state.unique = actnRes.unique;
		UpdateDPData(d, "Dcm Dcp free matrix free matrix TA TACount", &actnRes);
		/*Same test as above - make sure that all required transitions are
		present in the active subnet.*/
		j = 0;
		flag = ! d->TACount;
		for (i = 0; i < d->IncludeCount && ! flag; i++)
		{
			while(j < d->TACount - 1 && d->TA[j] < d->IncludeT[i]) j++;
			if (d->TA[j] != d->IncludeT[i])
			{
				flag = 1;
				break;
			}
		}
		if (flag)
		{
			d->state.anetfail = 1;
		}
	}
	/*Check to see if the task is impossible.*/
	if (! d->TACount)
	{
		d->terminate = 1;
		if (d->firstpass)
		{
			d->status = HOW_IMPOSSIBLE;
		}
		else
		{
			d->status = HOW_FAILED;
		}
		d->state.incon = 0;
	}
}

/*******************************************************************************
EnforceMarkingConstraints enforces the present marking constraints on the net.*/
void EnforceMarkingConstraints(dpData* d)
{
	if (NumberOfRows(d->L))
	{			
		supervis_r supervisRes = supervis(&d->Dp, &d->Dm, &d->L);
		UpdateDPData(d, "Dm Dp", &supervisRes);
		/*Create the dependent place list.*/
		d->dpCount = NumberOfRows(d->Dm) - d->targetRows;
		d->dpl = tcalloc(d->dpCount, sizeof(int));
		int i;
		for (i = d->targetRows; i < NumberOfRows(d->Dm); i++)
		{
			d->dpl[i - d->targetRows] = i;
		}		
		/*Setup M = L, G = B. Optimize for row operations.*/
		AllocateMatrixType(2, &d->M, NumberOfRows(d->L), NumberOfColumns(d->L));
		CopyMatrix(&d->L, &d->M);
		d->G = tcalloc(NumberOfRows(d->L), sizeof(int));
		memcpy(d->G, d->B, NumberOfRows(d->L) * sizeof(int));
	}
}

/*******************************************************************************
asbnmessage prints a per-major-iteration message.*/
void asbnmessage(dpData *d)
{
	if (d->run < 2)
	{
		printf("\nI");
	}
	else
	{
		printf("\nRun %d, i", d->run);
	}
	if (d->state.dptype < 2)
	{
		printf("teration 1: Computing the maximal active subnet...\n");
	}
	else
	{
		printf("teration 1: Computing a T-minimal active subnet...\n");
	}
}

/*******************************************************************************
This function prints the current state to the data and log files.*/
void PrintState(dpData *d)
{
	fprintf(d->log, "\nThe Petri net has %d places and %d transitions", NumberOfRows(d->Dm), NumberOfColumns(d->Dm));
	if (d->ExcludeT) fprintf(d->log, "\nThe current active subnet excludes the transitions %s", fvpr(d->ExcludeT, d->ExcludeCount, DISPLAY_INTS | DISPLAY_OFFSET, 2));
	
	int i, j;
	/*Build a list of the original places for vector display.*/
	if (d->targetRows)
	{
		int *originalPlaces = tcalloc(d->targetRows, sizeof(int));
		for (i = 1; i < d->targetRows; i++)
		{
			originalPlaces[i] = i;
		}
		/*Display it.*/		
		fprintf(d->log, "\n  %d original places: %s", d->targetRows, fvpr(originalPlaces, d->targetRows, DISPLAY_INTS | DISPLAY_OFFSET, 2));
		free(originalPlaces);
	}
	
	/*We need a list of all the dependent and independent places that have been
	added so far, that is, all those with indices greater than the maximum present
	in the original net.*/	
	int firstNewInd = -1, firstNewDep = -1;
	for (i = 0; i < d->ipCount; i++)
	{
		if (d->ipl[i] >= d->targetRows)
		{
			firstNewInd = i;
			break;
		}
	}
	for (i = 0; i < d->dpCount; i++)
	{
		if (d->dpl[i] >= d->targetRows)
		{
			firstNewDep = i;
			break;
		}
	}
	if (firstNewInd != -1)
	{	
		if (d->ipCount - firstNewInd == 1)
		{
			fprintf(d->log, "\n  One split place: %s", fvpr(d->ipl + firstNewInd, d->ipCount - firstNewInd, DISPLAY_INTS | DISPLAY_OFFSET, 2));
		}
		else if (d->ipCount - firstNewInd > 1)
		{
			fprintf(d->log, "\n  %d split places: %s", d->ipCount - firstNewInd, fvpr(d->ipl + firstNewInd, d->ipCount - firstNewInd, DISPLAY_INTS | DISPLAY_OFFSET, 2));
		}
	}
	if (firstNewDep != -1)
	{
		if (d->dpCount - firstNewDep == 1)
		{
			fprintf(d->log, "\n  One control place: %s", fvpr(d->dpl + firstNewDep, d->dpCount - firstNewDep, DISPLAY_INTS | DISPLAY_OFFSET, 2));
		}
		else if (d->dpCount - firstNewDep > 1)
		{
			fprintf(d->log, "\n  %d control places: %s", d->dpCount - firstNewDep, fvpr(d->dpl + firstNewDep, d->dpCount - firstNewDep, DISPLAY_INTS | DISPLAY_OFFSET, 2));
		}
	}
	fprintf(d->log, "\n");
	
	/*We want to display a modified version of M - one with null columns in place
	for each of the non-independent places, so that it has a column for each
	column of the net. Do this one row at a time. Allocate the row.*/
	int *MRow = tcalloc(NumberOfRows(d->Dm), sizeof(int));
	for (i = 0; i < NumberOfRows(d->M); i++)
	{
		memset(MRow, 0, sizeof(int) * NumberOfRows(d->Dm));
		for (j = 0; j < d->ipCount; j++)
		{
			MRow[d->ipl[j]] = GetMatrixEl(&d->M, i, j);
		}
		/*Display the row.*/
		fprintf(d->log, "\n  %d --> %s >= %d", d->dpl[i], avpr(MRow, NumberOfRows(d->Dm), DISPLAY_INTS | DISPLAY_OFFSET, 0, 0), d->G[i]);
	}
	free(MRow);
	
	/*Build a Petri net object we can use with displaypn to print the current net
	to the dat file.*/
	pns currentPn = createpn("mDm mDp", d->Dm, d->Dp);
	/*Do the display.*/
	fprintf(d->dat, "\n");
	displaypn(currentPn, d->dat);
	fprintf(d->dat, "\n");
	deallocpn(&currentPn);
	fflush(d->log);
	fflush(d->dat);
}

/*******************************************************************************
This function prints a status message with information about current run and
iteration.*/
void PrintMessage(dpData* d, char* message)
{
	if (d->iteration > 0)
	{
		if (d->run < 2)
		{
			printf("Iteration %d: %s\n", d->iteration, message);
		}
		else
		{
			printf("Run %d, iteration %d: %s\n", d->run, d->iteration, message);
		}
	}
	else
	{
		printf("%s\n", message);
	}
}

/*******************************************************************************
This function modifies the dpData structure passed to it to set it up for a
single run. It examines the firstpass bit to determine whether or not certain
matrices have already been set up correctly.*/
void InitDpData(dp_p *p, dpData* d)
{
	/*If this is the first pass, then the matrices have been properly initialized.
	Otherwise, make them new copies of the original parameter matrices.*/
	if (! d->firstpass)
	{
		DeallocateMatrix(&d->Dm);
		AllocateMatrixType(2, &d->Dm, NumberOfRows(p->Dm), NumberOfColumns(p->Dm));
		CopyMatrix(&p->Dm, &d->Dm);
		DeallocateMatrix(&d->Dp);
		AllocateMatrixType(2, &d->Dp, NumberOfRows(p->Dp), NumberOfColumns(p->Dp));
		CopyMatrix(&p->Dp, &d->Dp);
		DeallocateMatrix(&d->DI);
		d->DI = SubtractMatrix(&d->Dp, &d->Dm, (matrix*) 2);
		/*If either of the constraint matrices is not present, initialize it to have
		the right number of columns and zero rows so that it can be concatenated with
		the other to form the total constraint matrix.*/
		if (p->L.type)
		{
			DeallocateMatrix(&d->L);
			AllocateMatrixType(2, &d->L, NumberOfRows(p->L), NumberOfColumns(p->L));
			CopyMatrix(&p->L, &d->L);
		}
		else
		{
			DeallocateMatrix(&d->L);
			AllocateMatrixType(2, &d->L, 0, NumberOfRows(p->Dm));
		}
		if (p->L0.type)
		{
			DeallocateMatrix(&d->L0);
			AllocateMatrixType(2, &d->L0, NumberOfRows(p->L0), NumberOfColumns(p->L0));
			CopyMatrix(&p->L0, &d->L0);
		}
		else
		{
			DeallocateMatrix(&d->L0);
			AllocateMatrixType(2, &d->L0, 0, NumberOfRows(p->Dm));
		}
		if (p->B)
		{
			free(d->B);
			d->B = tcalloc(NumberOfRows(d->L), sizeof(int));
			memcpy(d->B, p->B, NumberOfRows(d->L) * sizeof(int));
		}
		if (p->B0)
		{
			free(d->B0);
			d->B0 = tcalloc(NumberOfRows(d->L0), sizeof(int));
			memcpy(d->B0, p->B0, NumberOfRows(d->L0) * sizeof(int));
		}
	}
	
	/*Initialize the original, independent, and dependent place lists.*/
	d->opCount = d->targetRows;
	d->ipCount = d->targetRows;
	d->dpCount = 0;
	if (d->opl) free(d->opl);
	if (d->ipl) free(d->ipl);
	if (d->dpl) free(d->dpl);
	d->opl = tcalloc(d->opCount, sizeof(int));
	d->ipl = tcalloc(d->ipCount, sizeof(int));
	d->dpl = 0;
	int i;
	for (i = 0; i < d->targetRows; i++)
	{
		d->opl[i] = i;
		d->ipl[i] = i;
	}

	/*Empty G and M.*/
	if (d->G) free(d->G);
	d->G = 0;
	if (d->M.type) DeallocateMatrix(&d->M);

	/*Initialize M to be the right size for an msplit.*/
	AllocateMatrix(&d->M, 0, d->targetRows);
	
	/*Remove from the exclusion list any transitions beyond those of the original
	matrices. This can be done by finding the index of the first list element that
	is too high and marking this as the list length.*/
	for (i = 0; i < d->ExcludeCount; i++)
	{
		if (d->ExcludeT[i] >= d->targetCols) break;
	}
	d->ExcludeCount = i;
	d->state.perm = 1;
	d->status = HOW_OK;
}
		
/*******************************************************************************
This function fills in the dpData structure containing information for a dp run
prior to the first run.*/
dpData CreateDpData(dp_p *p)
{
	dpData d;
	memset(&d, 0, sizeof(dpData));
	/*Copy matrices*/
	AllocateMatrixType(2, &d.Dm, NumberOfRows(p->Dm), NumberOfColumns(p->Dm));
	CopyMatrix(&p->Dm, &d.Dm);
	AllocateMatrixType(2, &d.Dp, NumberOfRows(p->Dp), NumberOfColumns(p->Dp));
	CopyMatrix(&p->Dp, &d.Dp);
	/*Get the target matrix DI.*/	
	d.DI = SubtractMatrix(&d.Dp, &d.Dm, (matrix*) 2);
	int i;
	
	/*If either of the constraint matrices is not present, initialize it to have
	the right number of columns and zero rows so that it can be concatenated with
	the other to form the total constraint matrix.*/
	if (p->L.type)
	{
		AllocateMatrixType(2, &d.L, NumberOfRows(p->L), NumberOfColumns(p->L));
		CopyMatrix(&p->L, &d.L);
	}
	else
	{
		AllocateMatrixType(2, &d.L, 0, NumberOfRows(p->Dm));
	}
	if (p->L0.type)
	{
		AllocateMatrixType(2, &d.L0, NumberOfRows(p->L0), NumberOfColumns(p->L0));
		CopyMatrix(&p->L0, &d.L0);
	}
	else
	{
		AllocateMatrixType(2, &d.L0, 0, NumberOfRows(p->Dm));
	}
	d.targetRows = NumberOfRows(p->Dm);
	d.targetCols = NumberOfColumns(p->Dm);
	/*Copy index lists (and fill in the corresponding flags).*/
	if (p->Tuc && p->TucCount)
	{
		d.Tuc = tcalloc(p->TucCount, sizeof(int));
		d.TucCount = p->TucCount;
		memcpy(d.Tuc, p->Tuc, sizeof(int) * p->TucCount);
		/*Make sure the index array is sorted. This will allow us to speed up
		certain processing later.*/
		Sort(d.Tuc, d.TucCount);
	}
	if (p->Tuo && p->TuoCount)
	{
		d.Tuo = tcalloc(p->TuoCount, sizeof(int));
		d.TuoCount = p->TuoCount;
		memcpy(d.Tuo, p->Tuo, sizeof(int) * p->TuoCount);
		/*Make sure the index array is sorted. This will allow us to speed up
		certain processing later.*/
		Sort(d.Tuo, d.TuoCount);
	}
	if (p->ExcludeT && p->ExcludeCount)
	{
		d.ExcludeT = tcalloc(p->ExcludeCount, sizeof(int));
		d.ExcludeCount = p->ExcludeCount;
		memcpy(d.ExcludeT, p->ExcludeT, sizeof(int) * p->ExcludeCount);
		/*Make sure the index array is sorted. This will allow us to speed up
		certain processing later.*/
		Sort(d.ExcludeT, d.ExcludeCount);
	}
	if (p->IncludeT && p->IncludeCount)
	{
		d.IncludeT = tcalloc(p->IncludeCount, sizeof(int));
		d.IncludeCount = p->IncludeCount;
		memcpy(d.IncludeT, p->IncludeT, sizeof(int) * p->IncludeCount);
		/*Make sure the index array is sorted. This will allow us to speed up
		certain processing later.*/
		Sort(d.IncludeT, d.IncludeCount);
	}
	else
	{
		/*If IncludeT is not present, it defaults to including all transitions.*/
		d.IncludeT = tcalloc(d.targetCols, sizeof(int));
		d.IncludeCount = d.targetCols;
		for (i = 0; i < d.targetCols; i++)
		{
			d.IncludeT[i] = i;
		}
	}
	/*Copy the constraint vectors*/
	if (p->B)
	{
		d.B = tcalloc(NumberOfRows(d.L), sizeof(int));
		memcpy(d.B, p->B, NumberOfRows(d.L) * sizeof(int));
	}
	if (p->B0)
	{
		d.B0 = tcalloc(NumberOfRows(d.L0), sizeof(int));
		memcpy(d.B0, p->B0, NumberOfRows(d.L0) * sizeof(int));
	}
	/*The only one of the state flags that does not default to zero is perm.*/
	d.state.perm = 1;
	/*If the DP_OPT_ENFORCELIVE option bit is set, we are doing liveness
	enforcement. If the IncludeT count is less than the total number of
	transitions, we are doing T-liveness enforcement.*/
	if (p->option & DP_OPT_ENFORCELIVE)
	{
		if (d.targetCols > d.IncludeCount)
		{
			d.state.dptype = 2;
		}
		else
		{
			d.state.dptype = 1;
		}
	}
	/*If Tuc or Tuo are present, set the uncontro state bit.*/
	if (d.Tuc || d.Tuo)
	{
		d.state.uncontro = 1;
	}
	/*If L or L0 are present (have at least one row), set the initial conditions 
	(incon) state bit.*/
	if (d.L.nr || d.L0.nr)
	{
		d.state.incon = 1;
		d.consis = 1;
	}
	/*If the DP_OPT_WRITELOG bit is set, attempt to open the log files. If either
	open fails, revert to no logging and issue a warning if verbosity is turned on.*/
	if (p->option & DP_OPT_WRITELOG)
	{
		char *mode = (p->option & DP_OPT_APPENDLOGS) ? "a" : "w";
		char *name = (p->option & DP_OPT_APPENDLOGS) ? "append" : "write";
		if (d.log = fopen("dp.log", mode))
		{
			if (! (d.dat = fopen("dp.dat", mode)))
			{
				if (is_verbose() >= VRB_DP) printf("Warning: Could not open file \"dp.dat\" for %s access. Logging disabled.\n", name);
				fclose(d.log);
				d.log = 0;				
			}
		}
		else if (is_verbose() >= VRB_DP) printf("Warning: Could not open file \"dp.log\" for %s access. Logging disabled.\n", name);
	}
	/*If logging is enabled go ahead and write the log headers.*/
	if (d.log)
	{
		char timeString[64];
		/*Matlab date format example: 24-Oct-2003 13:45:07*/
		time_t tm = time(0);
		strftime(timeString, 64, "%d-%b-%Y %H:%M:%S", localtime(&tm));
		fprintf(d.log, "\n%s -- DP.log -- Generated by DP\n", timeString);
		fprintf(d.dat, "\n%s -- DP.DAT -- Generated by DP\n", timeString);
		fflush(d.log);
		fflush(d.dat);
	}
	/*Set firstpass to 1 and leave restart as 0.*/
	d.firstpass = 1;
	return d;
}

/*******************************************************************************
This function checks the constraints given in L/B and L0/B0 of the dpData
structure and builds the set of complete marking constraints from them, returning
it in the MarkingConstraints structure.*/
MarkingConstraints UpdateMarkingConstraints(dpData *d)
{
	int count, cols;
	MarkingConstraints result;
	memset(&result, 0, sizeof(MarkingConstraints));
	if (NumberOfRows(d->L) || NumberOfRows(d->L0))
	{		
		count = NumberOfRows(d->L) + NumberOfRows(d->L0);
		cols = NumberOfColumns(d->L) ? NumberOfColumns(d->L) : NumberOfColumns(d->L0);
		AllocateMatrixType(2, &result.L, count, cols);
		CopyBlock(NumberOfRows(d->L), cols, &d->L, 0, 0, &result.L, 0, 0);
		CopyBlock(NumberOfRows(d->L0), cols, &d->L0, 0, 0, &result.L, NumberOfRows(d->L), 0);
	}
	if (count && (d->B || d->B0))
	{
		result.B = tcalloc(count, sizeof(int));
		if (d->B)
		{
			memcpy(result.B, d->B, sizeof(int) * NumberOfRows(d->L));
		}
		if (d->B0)
		{
			memcpy(result.B + NumberOfRows(d->L), d->B0, sizeof(int) * NumberOfRows(d->L0));
		}
	}
	return result;
}


/*******************************************************************************
This function is used to fill the appropriate elements of the dpData structure
with replacements from one of the structures returned by an spnbox function.
This is done by passing it a pointer to the dp data structure, a string 
describing where parts of the returned structure should be put, followed by a
pointer to the return structure itself.
Format of the description string is a series of whitespace-separated tokens,
each token representing an element of dpData. Tokens may the names of any element
in dpData. Additionally, tokens may be any of the following:
matrix         a dummy matrix, used to skip a piece of the newData structure the
               size of a matrix. It will be DEALLOCATED!
pointer        a dummy pointer. It will be DEALLOCATED!
staticpointer  a dummy pointer which will not be deallocated.
int            a dummy integer*/
void UpdateDPData(dpData *data, char* description, void *newData)
{
	void* next = newData;
	char* desc = description;
	char token[16];
	int i, transitions = NumberOfColumns(data->Dm);
	
	/*Go until we reach the terminating null.*/
	while (*desc)
	{
		/*Destination pointers.*/
		matrix* destMatrix = 0;
		int*    destInt = 0;
		int**   destIntArray = 0;
		int***  destTD = 0;		
		/*If we encounter an ignore instruction, we want to ignore the next
		paramater. If we encounter a deallocate instruction, we want to deallocate
		the space used by the next parameter.*/
		int ignore = 0, deallocate = 0;		
		do
		{
			/*Iterate past spaces.*/
			while (*desc == ' ' && *desc) desc++;
			/*Copy the current token.*/
			i = 0;
			while (*desc != ' ' && *desc)
			{
				token[i++] = *(desc++);
			}
			token[i] = '\0';
			if (! strcmp(token, "ignore")) ignore = 1;
			if (! strcmp(token, "free")) deallocate = 1;
		} while (!(strcmp(token, "ignore") && strcmp(token, "free")));
		
		/*Examine the token and take appropriate action.*/
		if (token[0])
		{
			if (! strcmp(token, "TD"))
			{
				if (ignore) destTD = (int***) 1;
				else if (deallocate) destTD = next;
			}
			else if (! strcmp(token, "Dm")) destMatrix = &data->Dm;
			else if (! strcmp(token, "Dp")) destMatrix = &data->Dp;
			else if (! strcmp(token, "DI")) destMatrix = &data->DI;
			else if (! strcmp(token, "Dcm")) destMatrix = &data->Dcm;
			else if (! strcmp(token, "Dcp")) destMatrix = &data->Dcp;
			else if (! strcmp(token, "L")) destMatrix = &data->L;
			else if (! strcmp(token, "L0")) 
			{
				destMatrix = &data->L0;
			}
			else if (! strcmp(token, "M")) destMatrix = &data->M;
			else if (! strcmp(token, "TSIPH")) destMatrix = &data->TSIPH;
			else if (! strcmp(token, "S")) destMatrix = &data->S;
			else if (! strcmp(token, "matrix"))
			{
				if (ignore) destMatrix = (matrix*) 1;
				else if (deallocate) destMatrix = next;
				else
				{
					merror(0, "DP: Internal Error: 'matrix' encounted without a 'free' or 'ignore' directive");
				}
			}			
			else if (! strcmp(token, "B")) destIntArray = &data->B;
			else if (! strcmp(token, "B0")) destIntArray = &data->B0;
			else if (! strcmp(token, "G")) destIntArray = &data->G;
			else if (! strcmp(token, "opl")) destIntArray = &data->opl;
			else if (! strcmp(token, "ipl")) destIntArray = &data->ipl;
			else if (! strcmp(token, "dpl")) destIntArray = &data->dpl;
			else if (! strcmp(token, "apl")) destIntArray = &data->apl;
			else if (! strcmp(token, "ExcludeT")) destIntArray = &data->ExcludeT;
			else if (! strcmp(token, "IncludeT")) destIntArray = &data->IncludeT;
			else if (! strcmp(token, "TA")) destIntArray = &data->TA;
			else if (! strcmp(token, "Tuc")) destIntArray = &data->Tuc;
			else if (! strcmp(token, "Tuo")) destIntArray = &data->Tuo;
			else if (! strcmp(token, "upl")) destIntArray = &data->upl;
			else if (! strcmp(token, "intarray"))
			{
				if (ignore) destIntArray = (int**) 1;
				else if (deallocate) destIntArray = next;
				else
				{
					merror(0, "DP: Internal Error: 'intarray' encountered without a 'free' or 'ignore' directive");
				}
			}
			else if (! strcmp(token, "opCount")) destInt = &data->opCount;
			else if (! strcmp(token, "ipCount")) destInt = &data->ipCount;
			else if (! strcmp(token, "dpCount")) destInt = &data->dpCount;
			else if (! strcmp(token, "apCount")) destInt = &data->apCount;
			else if (! strcmp(token, "upCount")) destInt = &data->upCount;			
			else if (! strcmp(token, "ExcludeCount")) destInt = &data->ExcludeCount;
			else if (! strcmp(token, "IncludeCount")) destInt = &data->IncludeCount;
			else if (! strcmp(token, "TACount")) destInt = &data->TACount;
			else if (! strcmp(token, "TucCount")) destInt = &data->TucCount;
			else if (! strcmp(token, "TuoCount")) destInt = &data->TuoCount;
			else if (! strcmp(token, "int")) destInt = (int*) 1;
			else
			{
				printf("DP: UpdateDPData: Unrecognized token '%s'.\n", token);
				printf("Update state: ignore = %d, free = %d.\n", ignore, deallocate);
				merror(0, "DP: Internal Error");
				return;
			}
			if (destMatrix)
			{				
				if ((! ignore) && destMatrix->type) DeallocateMatrix(destMatrix);
				if (! (deallocate || ignore)) *destMatrix = *((matrix*)next);
				destMatrix = next;
				next = ++destMatrix;
			}
			if (destInt)
			{
				if (! ignore) *destInt = *((int*)next);
				destInt = next;
				next = ++destInt;
			}
			if (destIntArray)
			{				
				if (! ignore) free(*destIntArray);
				if (! (deallocate || ignore)) *destIntArray = *((int**)next);
				destIntArray = next;
				next = ++destIntArray;
			}
			if (destTD)
			{				
				if (*destTD && ! ignore)
				{
					for (i = 0; i < transitions; i++)
					{
						free((*destTD)[i]);
					}
					free(*destTD);
				}
				if (! (deallocate || ignore)) *destTD = *((int***)next);
				destTD = next;
				next = ++destTD;
			}
		}
	}
}

/*******************************************************************************
This function adds an item to an index array, increments the count, and makes
sure the array is sorted properly. Note the Sort will make at only one pass
through an array that is already properly sorted, so this is low-cost.*/
void AddToIndexArray(int** array, int* length, int item)
{
	int *newArray = tcalloc(*length + 1, sizeof(int));
	if (length) memcpy(newArray, *array, sizeof(int) * (*length));
	newArray[(*length)++] = item;
	free(*array);
	*array = newArray;
	Sort(*array, *length);
}

/*******************************************************************************
This function adds an item to an integer array.*/
void AddToArray(int** array, int length, int item)
{
	int* newArray = tcalloc(length + 1, sizeof(int));
	if (length) memcpy(newArray, *array, sizeof(int) * length);
	newArray[length] = item;
	free(*array);
	*array = newArray;
}

/*******************************************************************************
This function does a bubble sort on the array, reordering it from lowest to
highest.*/
void Sort(int* array, int length)
{
	int i, swapped = 1, temp;
	while (swapped)
	{
		swapped = 0;
		for (i = 1; i < length; i++)
		{
			if (array[i] < array[i - 1])
			{
				temp = array[i];
				array[i] = array[i - 1];
				array[i - 1] = temp;
				swapped = 1;
			}
		}
	}	
}

/*******************************************************************************
This function makes sure the parameters are valid. Returns nonzero if valid and
zero otherwise.*/
int CheckParams(dp_p *p)
{
	/*Dm/Dp are required and must be the same size.*/
	if (! (p->Dm.type && p->Dp.type))
	{
		merror(0, "DP: One or both of Dm and Dp are missing");
		return 0;
	}
	if (NumberOfRows(p->Dm) != NumberOfRows(p->Dp) || NumberOfColumns(p->Dm) != NumberOfColumns(p->Dp))
	{
		merror(0, "DP: Dm and Dp are not the same size");
		return 0;
	}
	/*If Tuc and/or Tuo are present, make sure there are no out-of-range indices.*/
	int i, j;
	if (p->Tuc)
	{
		for (i = 0; i < p->TucCount; i++)
		{
			if (p->Tuc[i] < 0 || p->Tuc[i] >= NumberOfColumns(p->Dp))
			{
				merror(0, "DP: There is an out-of-range index in Tuc");
				return 0;
			}
		}
	}
	if (p->Tuo)
	{
		for (i = 0; i < p->TuoCount; i++)
		{
			if (p->Tuo[i] < 0 || p->Tuo[i] >= NumberOfColumns(p->Dp))
			{
				merror(0, "DP: There is an out-of-range index in Tuo");
				return 0;
			}
		}
	}
	/*Make sure there are no out-of-range indices in the exclusion or inclusion
	lists.*/
	if (p->ExcludeT)
	{
		for (i = 0; i < p->ExcludeCount; i++)
		{
			if (p->ExcludeT[i] < 0 || p->ExcludeT[i] >= NumberOfColumns(p->Dp))
			{
				merror(0, "DP: There is an out-of-range index in ExcludeT");
				return 0;
			}
		}
	}
	if (p->IncludeT)
	{
		for (i = 0; i < p->IncludeCount; i++)
		{
			if (p->IncludeT[i] < 0 || p->IncludeT[i] >= NumberOfColumns(p->Dp))
			{
				merror(0, "DP: There is an out-of-range index in IncludeT");
				return 0;
			}
		}
	}
	/*Make sure there are no duplicate entries in the exclusion/inclusion lists.*/
	if (p->IncludeT && p->ExcludeT)
	{
		for (i = 0; i < p->IncludeCount; i++)
		{
			for (j = 0; j < p->ExcludeCount; j++)
			{
				if (p->IncludeT[i] == p->ExcludeT[j])
				{
					merror(0, "DP: The inclusion and exclusion lists mark the same transition for both exclusion and inclusion");
					return 0;
				}
			}
		}
	}
	/*If constraint matrices are present, make sure that they have the proper
	number of columns.*/
	if (p->L.type)
	{
		if (NumberOfColumns(p->L) != NumberOfRows(p->Dm))
		{
			merror(0, "DP: The marking constraint matrix L does not have the proper size");
			return 0;
		}
	}
	if (p->L0.type)
	{
		if (NumberOfColumns(p->L0) != NumberOfRows(p->Dm))
		{
			merror(0, "DP: The initial marking constraint matrix L0 does not have the proper size");
			return 0;
		}
	}
	/*If constraint vectors are present, the constraint matrices must be present.*/
	if (p->B && ! p->L.type)
	{
		merror(0, "DP: The marking constraint vector B is present but the marking constraint matrix L is missing");
		return 0;
	}
	if (p->B0 && ! p->L0.type)
	{
		merror(0, "DP: The initial marking constraint vector B0 is present but the initial marking constraint matrix L0 is missing");
		return 0;
	}
	return 1;
}
