#include <stdio.h>
#include <omp.h>
#include <float.h>
#include <malloc.h>
#include <stdlib.h>
// #include <time.h>
#define maxThreads 30
// #define expectedVertices 12
#define maxWidth 400
#define maxHeight 300
#define currentMaxDegree 5
#define sizeOfGraph 73
#define maxFaceSize 10
#define maxNumberOfFaces 20


#define allInformation tableau,vertices,anglesUpTo,vertexDegrees,faceArray,numberOfFaces,heightOfTableau,widthOfTableau,&currentRow,&currentSlackVariable
#define addInequalityToAnglesInfo tableau,currentRow,currentSlackVariable,heightOfTableau,widthOfTableau,inequalityValue


#define ALL_VERTEX 100
#define ALL_FACE 101
#define ADJACENT_FACE 102
#define OPPOSITE_FACE 103


#define INDIVIDUAL_ANGLE 300
#define SUM_ANGLE 301
#define SUBTRACT_ANGLE 302

#define GEQ 400
#define LEQ 401


#define separatorASCIICode 32

int expectedVertices=0;
int totalCorrectGraphs=0;
int totalProcessedGraphs = 0;
int outputAsASCII = 0;
int feasableSimplex(float tableau[maxHeight][maxWidth],int heightOfTableau,int widthOfTableau)
{
	
	/*
	Input arguments: heightOfTableau is the height of the tableau not including the objective row
	widthOfTableau is the number of variables in the tableau
	Runs the simplex algorithm on a tableau and returns false
	if the objective function is not sufficiently small at the end of execution*/
	int pivotColumn,pivotRow; //provides the location of the pivot

	int continueSimplex = 1;
	int k=0;
	float currentMaximum; //used to calculate the maximum from an array
	float currentMinimum; //used to calculate the minimum from an array
	float valueAtPivotColumn; //used to calculate the pivot row
	float pivotValue;
	int i,j; //inner and outerloop variables
	float objectiveValue;
	float subtractionFactor = 0; //used when subtracting one row from another
	float subtractionAmount = 0;
	while(1)
	{
		//Step 1 - Determine which column to pivot on (or whether we are optimal)
		currentMaximum = 0;
		pivotColumn = -1;
		for(i = 0; i < widthOfTableau; i++)
		{
			if(tableau[heightOfTableau][i] > currentMaximum)
			{
				currentMaximum = tableau[heightOfTableau][i];
				pivotColumn = i;
			}
		}
		if(pivotColumn == -1)
		{
			break; 
		}//we cannot optimize any further
		/*Step 2 - Determine which row to pivot on (or whether the problem is unbounded)*/
		currentMinimum = FLT_MAX;
		pivotRow = -1;
		//test each row for minimum ratio
		for(i = 0; i < heightOfTableau; i++)
		{
			valueAtPivotColumn = tableau[i][pivotColumn]; 
			/*valueAtPivotColumn stores the value of the tableau on row i, column pivotColumn*/
			if(valueAtPivotColumn>0 && (((tableau[i][widthOfTableau])/valueAtPivotColumn)<currentMinimum))
			{
				currentMinimum = (tableau[i][widthOfTableau])/valueAtPivotColumn;
				pivotRow = i;
			}
		}
		if(pivotRow ==-1)
		{
			#pragma omp critical(printerror)
			{
				fprintf(stderr,"An unbounded problem was found\n");
			}
			break;
		}

		/*Step 3 - Divide every entry in the pivot row by the pivot number*/
		pivotValue = tableau[pivotRow][pivotColumn];

		for(i = 0; i < (widthOfTableau + 1); i++)
		{
			tableau[pivotRow][i] = (tableau[pivotRow][i])/pivotValue;
		}

		/*Step 4 - Subtract multiple of new pivot row from each other row*/

		for(i = 0; i < (heightOfTableau + 1); i++)
		{
			subtractionFactor = tableau[i][pivotColumn];
			if(i!=pivotRow && (subtractionFactor !=0))
			{
				for(j=0;j<(widthOfTableau+1);j++)
				{
					subtractionAmount = subtractionFactor * (tableau[pivotRow][j]);
				
					tableau[i][j]  -=  subtractionAmount ;
				}
			}
		}
				
	}
	objectiveValue = tableau[heightOfTableau][widthOfTableau];
	if(objectiveValue>-0.001 && objectiveValue<0.001)
	{
		return 1;
	}
	return 0;
	//check whether or not we should allow the graph

}
void splitGraphToVertices(int graph[sizeOfGraph],int lengthOfGraph,int vertexDegrees[expectedVertices],int vertices[expectedVertices][currentMaxDegree])
{

	int currentEdge = 0; //used to temporarily store the current edge or a 0 if the end of the vertex is specified
	int currentVertexPosition = 0; //the current degree of the vertex
	int currentGraphEntry=0; //the current entry in the graph
	int vertexNumber = 0; //how many vertices have we currently processed
	int edgeProcessed = 0; //how many edges have we processed for the current vertex
	for(currentGraphEntry = 1;currentGraphEntry < lengthOfGraph; currentGraphEntry++)
	{
		currentEdge = graph[currentGraphEntry];
		if(currentEdge != 0)
		{
			vertices[vertexNumber][currentVertexPosition] = currentEdge - 1; //plantri outputs vertices starting from vertex 1, but we would prefer to label them from 0
			currentVertexPosition++;
		}
		else
		{
			//we have finished processing this vertex
			vertexDegrees[vertexNumber] = currentVertexPosition;
			currentVertexPosition = 0;
			vertexNumber++; 
		}
	}
	if(vertexNumber != expectedVertices)
	{
		fprintf(stderr,"A graph was passed to splitGraphToVertices with the incorrect number of vertices. %i vertices were passed and %i vertices were expected\n",vertexNumber,expectedVertices);
	}
}
int findNextVertex(int vertices[expectedVertices][currentMaxDegree],int previousVertex,int nextVertex,int currentVertexDegree)
{

	/* We consider the list of vertices connected to the vertex we are currently at. 
	Once we have found the edge that we just took to reach this vertex (by finding the previous vertex in the aforementioned list)
	we take the next edge in the list to get the next vertex (clockwise)
	*/
	int i = 0;
	int returnValue = -1;
	for(i = 0; i < (currentVertexDegree - 1); i++)
	{
		if((vertices[nextVertex][i]) == previousVertex)
		{
			//We have found the edge that we just took
			returnValue = vertices[nextVertex][i+1]; //this is the next edge clockwise
			break;
		}
	}
	if(i == (currentVertexDegree - 1))
	{
		/*the final edge of the current vertex connects to the previous vertex: 
		thus the next edge to traverse clockwise is the first edge listed*/
		
		returnValue = vertices[nextVertex][0];
	}
	return returnValue;
}
void getAnglesUpTo(int vertexDegrees[expectedVertices],int anglesUpTo[expectedVertices])
{
	int i;
	anglesUpTo[0]=0;
	for(i=1;i<expectedVertices;i++)
	{
		anglesUpTo[i]=anglesUpTo[i-1]+vertexDegrees[i-1];
	}
}
int getAngleID(int vertices[expectedVertices][currentMaxDegree],int vertexDegrees[expectedVertices],int anglesUpTo[expectedVertices],int centreVertex,int anticlockwiseVertex)
{
	//gets the angleID for the angle firstVertex centreVertex clockwiseVertex
	int i = 0; //used to iterate through all vertices connected to centreVertex
	int degreeOfCentreVertex = vertexDegrees[centreVertex];
	int angleId;
	for(i=0;i<degreeOfCentreVertex;i++)
	{
		if(vertices[centreVertex][i]==anticlockwiseVertex)
		{
			angleId = (i + anglesUpTo[centreVertex]);
			return angleId;
		}
	}
	
	fprintf(stderr,"There was an error trying to setup an inequality: the required angle could not be found");
	return -1;

}
void addInequalityToAngles(float tableau[maxHeight][maxWidth],int *currentRow,int *currentSlackVariable,int heightOfTableau,int widthOfTableau
	,float inequalityValue,int totalAngles,int anglesToApplyTo[],int coefficiants[])
{
	/*Adds an inequality for a collection of angles (anglesToApplyTo) to the table. The inequality given is of the form sum a_i c_i + s_i = inequalityValue, where a_i is the angle and c_i is the coefficiant
	Accepts only positive floats for inequalityValue
	Modifies the objective row if an artificial variable is necessary.
	*/
	//iterator variable: used to clear the current row and add angles
	float multiplier = 1; //if inequalityValue < 0 we will simply multiply the entire inequality by -1 and include an arbitrary variable.
	float valueToInsert; //used as temporary storage for the value that needs to be inserted in the tableau after multiplication 
	int i;
	if(inequalityValue < 0) multiplier = -1;
	for(i = 0; i < widthOfTableau;i++)
	{
		tableau[*currentRow][i] = 0; //set the entire row to 0
	}
	for(i = 0; i < totalAngles;i++)
	{
		valueToInsert = multiplier * coefficiants[i];
		tableau[*currentRow][anglesToApplyTo[i]] = valueToInsert; //insert each angle
	}
	tableau[*currentRow][*currentSlackVariable] = multiplier; //the slack variable is positive or negative depending on whether or not we are adding a row with a negative inequalityValue
	valueToInsert = inequalityValue * multiplier;
	tableau[*currentRow][widthOfTableau] = valueToInsert;
	if(multiplier == -1)
	{
		//we now need to add this entire row to the objective row to cancel any effects of artificial variables
		for(i = 0; i<totalAngles; i++)
		{
			valueToInsert = multiplier * coefficiants[i];
			tableau[heightOfTableau][anglesToApplyTo[i]] += valueToInsert; //insert each angle
		}
		tableau[heightOfTableau][*currentSlackVariable] = multiplier;
		
		valueToInsert = inequalityValue * multiplier;
		tableau[heightOfTableau][widthOfTableau] += valueToInsert;
	}
	//make sure that this row is not overwritten
	*currentSlackVariable = (*currentSlackVariable) + 1;
	*currentRow = (*currentRow) + 1;
}
int addInequality(float tableau[maxHeight][maxWidth], int vertices[expectedVertices][currentMaxDegree], int anglesUpTo[expectedVertices], int vertexDegrees[expectedVertices],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],
	int numberOfFaces[maxFaceSize],int heightOfTableau, int widthOfTableau,  int *currentRow, int *currentSlackVariable,
	int applicableAngles, int operationOnAngles, int faceSize, float inequalityValue, int typeOfInequality)
{
	/*Adds the inequality specified by the parameters to tableau. 
	applicableAngles is either ALL_VERTEX,ALL_FACE,ADJACENT_FACE or OPPOSITE_FACE depending on which angles the inequality applies to
	operation is either INDIVIDUAL_ANGLE SUM_ANGLE OR SUBTRACT_ANGLE depending on whether the left hand side of the inequality refers to individual angles, the sum of the angles or one angle subtract the other
	faceSize is the faces the inequality applies to: 3 refers to triangles, 4 refers to quadrilaterals e.t.c
	currentRow is the currentRow of the tableau: it is updated by this proceedure
	objectiveRow is the position of the objective row in the tableau: it is used so that no slack variables are left behind
	inequality value is the right hand side of the inequality

	Returns -1 if the inequality was not successfully added, and 0 if it was.
	*/
	
	int i,j; //iterator variables
	int *anglesToApplyTo; //the angles to send to addInequalityToAngles
	int *coefficiants; //the coefficiants to send to addInequalityToAngles
	int totalAngles = anglesUpTo[(expectedVertices - 1)] + vertexDegrees[(expectedVertices - 1)]; 
	int numberOfAnglesInInequality;
	int halfFaceSize = faceSize / 2; //integer division here is deliberate: we will use if required later to find out whether or not the face has an even number of sides (used for opposite angle code)
	int previousVertex,nextVertex; //used for adding an inequality involving adjacent angles (prevents accessing the array more than necessary) 
	int previousOppositeVertex,nextOppositeVertex; //used for adding an inequality involving adjacent angles
	int multiplier = 1;
	if(typeOfInequality == GEQ)
	{
		multiplier = -1;
	}
	inequalityValue *= multiplier;
	if(applicableAngles == ALL_VERTEX)
	{
		//applies to all angles tableau,currentRow,currentArtificialVariable,currentSlackVariable,heightOfTableau,widthOfTableau,inequalityValue
		if(operationOnAngles == INDIVIDUAL_ANGLE)
		{
			//each individual angle has a separate inequality that applies to it
			anglesToApplyTo = (int *)malloc(sizeof(int));
			coefficiants = (int *)malloc(sizeof(int));
			for(i = 0; i < totalAngles; i++)
			{
				anglesToApplyTo[0]=i;
				coefficiants[0]=multiplier;
				addInequalityToAngles(addInequalityToAnglesInfo,1,anglesToApplyTo,coefficiants);
			}
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUM_ANGLE)
		{
			//the sum of all angles is on the LHS of the inequality
			
			anglesToApplyTo = (int *)malloc(currentMaxDegree* sizeof(int));
			coefficiants = (int *)malloc(currentMaxDegree* sizeof(int));
			for(i = 0; i < expectedVertices;i++)
			{
				numberOfAnglesInInequality = vertexDegrees[i];
				for(j = 0; j < numberOfAnglesInInequality;j++)
				{
					anglesToApplyTo[j] = anglesUpTo[i] + j; 
					coefficiants[j] = multiplier;
				}
				
				addInequalityToAngles(addInequalityToAnglesInfo,numberOfAnglesInInequality,anglesToApplyTo,coefficiants);
			}
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUBTRACT_ANGLE)
		{
			#pragma omp critical(printerror)
			{

				fprintf(stderr,"The inequality operation SUBTRACT_ANGLE was requested, but more than two angles were used.");
			}
			return -1;//it does not make sense to refer to SUBTRACT_ANGLE if there are more than 2 angles
		}
		else
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"An incorrect inequality operation was requested: the operation requested was %i",operationOnAngles);
			}
			return -1; 
		}

	}
	else if(applicableAngles == ALL_FACE)
	{
		if(faceSize < 3)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"You cannot add inequalities involving faces of less than 3 sides");
			}
			return -1;
		}
		if(operationOnAngles == INDIVIDUAL_ANGLE)
		{
			
			anglesToApplyTo = (int *)malloc(sizeof(int));
			coefficiants = (int *)malloc(sizeof(int));
			for(i = 0;i < numberOfFaces[(faceSize - 1)];i++)
			{
				//get the interior angle about each vertex in the face
				//the first angle is tricky, so do that one separately
				anglesToApplyTo[0]=getAngleID(vertices,vertexDegrees,anglesUpTo,faceArray[(faceSize - 1)][i][0],faceArray[(faceSize - 1)][i][(faceSize - 1)]);
				coefficiants[0]=multiplier;
				addInequalityToAngles(addInequalityToAnglesInfo,1,anglesToApplyTo,coefficiants);
				for(j = 1; j < faceSize;j++)
				{
					anglesToApplyTo[0]=getAngleID(vertices,vertexDegrees,anglesUpTo,faceArray[(faceSize - 1)][i][j],faceArray[(faceSize - 1)][i][(j - 1)]);
					coefficiants[0]=multiplier;
					
					addInequalityToAngles(addInequalityToAnglesInfo,1,anglesToApplyTo,coefficiants);
					
				}
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUM_ANGLE)
		{
			anglesToApplyTo = (int *)malloc(faceSize * sizeof(int));
			coefficiants = (int *)malloc(faceSize * sizeof(int));
			for(i = 0; i<numberOfFaces[(faceSize - 1)]; i++)
			{
				//get the interior angle about each vertex in the face
				anglesToApplyTo[0] = getAngleID(vertices,vertexDegrees,anglesUpTo,faceArray[(faceSize - 1)][i][0],faceArray[(faceSize - 1)][i][(faceSize - 1)]);
				coefficiants[0] = multiplier;
				for(j = 1; j < faceSize;j++)
				{
					anglesToApplyTo[j] = getAngleID(vertices,vertexDegrees,anglesUpTo,faceArray[(faceSize - 1)][i][j],faceArray[(faceSize - 1)][i][(j - 1)]);
					coefficiants[j] = multiplier;
					
				}
				addInequalityToAngles(addInequalityToAnglesInfo,faceSize,anglesToApplyTo,coefficiants);
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUBTRACT_ANGLE)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"The inequality operation SUBTRACT_ANGLE was requested, but more than two angles were used.");
			}
			return -1;

		}
		else
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"An incorrect inequality operation was requested: the operation requested was %i",operationOnAngles);
			}
			return -1; //it does not make sense to refer to SUBTRACT_ANGLE if there are more than 2 angles
		}
	}
	else if(applicableAngles == ADJACENT_FACE)
	{
		if(faceSize <3)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"You cannot add inequalities involving faces of less than 3 sides");
			}
			return -1;
		}
		if(operationOnAngles == INDIVIDUAL_ANGLE)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"The inequality operation INDIVIDUAL_ANGLE was requested, but exactly two angles will be in each inequality");
			}
			return -1;
		}
		else if(operationOnAngles == SUM_ANGLE)
		{
			anglesToApplyTo = (int *)malloc(2 * sizeof(int));
			coefficiants = (int *)malloc(2 * sizeof(int));
			//since this requires summing exactly two angles, we might as well set the coefficiants here
			coefficiants[0] = multiplier;
			coefficiants[1] = multiplier;
			for(i = 0; i<numberOfFaces[(faceSize - 1)]; i++)
			{
				/*we need two vertices to know what angle is being referred to: we cycle through the ordered pairs of vertices in order to generate all inequalities involving adjacent angles. 
				In doing so, it is best to store the previous angle (and one of the vertices) to save looking it up again for the next inequality*/
				previousVertex = faceArray[(faceSize - 1)][i][(faceSize - 1)];
				nextVertex = faceArray[(faceSize - 1)][i][0];
				anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
				for(j = 0; j < faceSize;j++)
				{
					//add the first angle to the inequality
					anglesToApplyTo[0] = anglesToApplyTo[1];
					
					//add the second angle to the inequality
					previousVertex = nextVertex;
					if(j == (faceSize - 1))
					{
						nextVertex = faceArray[(faceSize - 1)][i][0];
					}
					else
					{
						nextVertex = faceArray[(faceSize -1)][i][(j+1)]; //get the next vertex if we haven't reached the end of the list of vertices.
					}
					anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
				}
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUBTRACT_ANGLE)
		{
			anglesToApplyTo = (int *)malloc(2 * sizeof(int));
			coefficiants = (int *)malloc(2 * sizeof(int));
			for(i = 0; i<numberOfFaces[(faceSize - 1)]; i++)
			{
				/*we need two vertices to know what angle is being referred to: we cycle through the ordered pairs of vertices in order to generate all inequalities involving adjacent angles. 
				In doing so, it is best to store the previous angle (and one of the vertices) to save looking it up again for the next inequality*/
				previousVertex = faceArray[(faceSize - 1)][i][(faceSize - 1)];
				nextVertex = faceArray[(faceSize - 1)][i][0];
				anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
				for(j = 0; j < faceSize;j++)
				{
					//add the first angle to the inequality
					anglesToApplyTo[0] = anglesToApplyTo[1];
					
					//add the second angle to the inequality
					previousVertex = nextVertex;
					if(j == (faceSize - 1))
					{
						nextVertex = faceArray[(faceSize - 1)][i][0];
					}
					else
					{
						nextVertex = faceArray[(faceSize -1)][i][(j+1)]; //get the next vertex if we haven't reached the end of the list of vertices.
					}
					anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
					//Now we know which angles to apply each inequality to. 
					coefficiants[0] = multiplier;
					coefficiants[1] = (-1 * multiplier); 
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
					//To do the second inequality that comes with this pair of angles, simply flip the coefficiants
					coefficiants[0] = (-1 * multiplier);
					coefficiants[1] = multiplier;
					
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
				}
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
	}
	else if(applicableAngles == OPPOSITE_FACE)
	{
		if(faceSize <3)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"You cannot add inequalities involving faces of less than 3 sides");
			}
			return -1;
		}
		if((halfFaceSize * 2) != faceSize)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"You cannot add inequalities for opposite angles if the face does not have an even number of sides");
			}
			return -1;
		}
		if(operationOnAngles == INDIVIDUAL_ANGLE)
		{

			#pragma omp critical(printerror)
			{
				fprintf(stderr,"The inequality operation INDIVIDUAL_ANGLE was requested, but exactly two angles will be in each inequality");
			}
			return -1;
		}
		else if(operationOnAngles == SUM_ANGLE)
		{
			anglesToApplyTo = (int *)malloc(2 * sizeof(int));
			coefficiants = (int *)malloc(2 * sizeof(int));
			//since this requires summing exactly two angles, we might as well set the coefficiants here
			coefficiants[0] = multiplier;
			coefficiants[1] = multiplier;
			for(i = 0; i<numberOfFaces[(faceSize - 1)]; i++)
			{
				previousVertex = faceArray[(faceSize - 1)][i][(faceSize - 1)];
				previousOppositeVertex = faceArray[(faceSize - 1)][i][(halfFaceSize - 1)]; //the opposite angle should be half the facesize away from the original angle
				for(j = 0; j < halfFaceSize;j++)
				{
					//work out what the next vertex is
					nextVertex = faceArray[(faceSize - 1)][i][j];
					anglesToApplyTo[0] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
					previousVertex = nextVertex;
					//now we need to get the vertex opposite
					nextOppositeVertex = faceArray[(faceSize - 1)][i][(j+halfFaceSize)];

					anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextOppositeVertex,previousOppositeVertex);
					previousOppositeVertex = nextOppositeVertex;
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
				}
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
		else if(operationOnAngles == SUBTRACT_ANGLE)
		{
			anglesToApplyTo = (int *)malloc(2 * sizeof(int));
			coefficiants = (int *)malloc(2 * sizeof(int));
			

			
			for(i = 0; i<numberOfFaces[(faceSize - 1)]; i++)
			{
				previousVertex = faceArray[(faceSize - 1)][i][(faceSize - 1)];
				previousOppositeVertex = faceArray[(faceSize - 1)][i][(halfFaceSize - 1)]; //the opposite angle should be half the facesize away from the original angle
				for(j = 0; j < halfFaceSize;j++)
				{
					
					//work out what the next vertex is
					nextVertex = faceArray[(faceSize - 1)][i][j];
					anglesToApplyTo[0] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextVertex,previousVertex);
					previousVertex = nextVertex;
					//now we need to get the vertex opposite
					nextOppositeVertex = faceArray[(faceSize - 1)][i][(j+halfFaceSize)];

					anglesToApplyTo[1] = getAngleID(vertices,vertexDegrees,anglesUpTo,nextOppositeVertex,previousOppositeVertex);
					previousOppositeVertex = nextOppositeVertex;
					coefficiants[0] = multiplier;
					coefficiants[1] = (-1 * multiplier);
					//set the coefficiants for one of the two inequalities added by this option
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
					coefficiants[0] = (-1 * multiplier);
					coefficiants[1] = multiplier;
					addInequalityToAngles(addInequalityToAnglesInfo,2,anglesToApplyTo,coefficiants);
				}
			}
			
			free(anglesToApplyTo);
			free(coefficiants);
		}
	}
	return 0;
}

void retrieveFaces(int vertexDegrees[expectedVertices],int vertices[expectedVertices][currentMaxDegree],int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize],int numberOfFaces[maxFaceSize])
{
	/*Input is vertexDegrees (the degrees of each vertex)
	and vertices (the degrees of each vertex and its neighbours in clockwise order).
	This procedure modifies the values of faceArray and numberOfFaces
	numberOfFaces stores the number of faces in each entry e.g. numberOfFaces[2] is the number of triangles found
	faceArray stores information about types of faces in the graph 
	e.g. faceArray[3][1][2] is the third vertex in the second face from the collection of quadrilateral faces that the graph contains
	*/
	int vertexNumber = 0; //which vertex are we currently processing 
	int nextVertex; //the next vertex in the face
	int previousVertex; //the previous vertex in the face
	int faceSides = 0;//stores the number of sides of the face
	int edgeProcessed; //the current edge that we are processing
	int currentFace[maxFaceSize];
	int i; //used for iteration
	int alreadyCountedFace=0;// if the face has already been counted then this flag will be triggered
	// Begin by setting the numberOfFaces to 0 : this way we can increment them
	for(i = 0; i<maxFaceSize;i++)
	{
		numberOfFaces[i] = 0;
	}
	/*The strategy for processing faceArray and faceArray lengths is to go through each vertex (except the last two - no face can contain just two vertices)
	and try to work out which face it is connected to. We must use the assumption of the clockwise ordering of the edges in vertices */
	for(vertexNumber = 0; vertexNumber < (expectedVertices - 2); vertexNumber++) 
	{
		for(edgeProcessed = 0; edgeProcessed < vertexDegrees[vertexNumber]; edgeProcessed++)
		{
			//We traverse the face until we return to the starting vertex
			previousVertex = vertexNumber;
			nextVertex = vertices[previousVertex][edgeProcessed];
			alreadyCountedFace = 0;
			/*Avoid faces which have already been counted*/
			
			if(!(nextVertex < previousVertex)) // if nextVertex < previousVertex then all faces attached to this vertex will have already been added
			{
				currentFace[0] = vertexNumber;

				faceSides = 1;

				while(nextVertex != vertexNumber)
				{

					currentFace[faceSides] = nextVertex; //add this vertex to the list of vertices traversed
					nextVertex = findNextVertex(vertices,previousVertex,nextVertex,vertexDegrees[nextVertex]);
					previousVertex = currentFace[faceSides];
					faceSides++;

					/*Avoid faces which have already been counted*/
					if(nextVertex<vertexNumber)
					{
						alreadyCountedFace = 1; 
						break;
					}
				}
				if(!alreadyCountedFace && faceSides > 2)
				{
					// we need to add this face to the list of faces 
					for(i = 0; i < faceSides; i++)
					{
						faceArray[(faceSides - 1)][numberOfFaces[(faceSides - 1)]][i] = currentFace[i];

					}
					numberOfFaces[(faceSides - 1)]++;
				}
			}
			
		}
	}
}

int makeTableau(int graph[sizeOfGraph],int lengthOfGraph)
{
	//takes a graph and produces a tableau ready for the simplex algorithm.
	int vertexDegrees[expectedVertices];
	int anglesUpTo[expectedVertices];
	float tableau[maxHeight][maxWidth];
	int vertices[expectedVertices][currentMaxDegree];
	int faceArray[maxFaceSize][maxNumberOfFaces][maxFaceSize];
	int numberOfFaces[maxFaceSize];
	int heightOfTableau;
	int widthOfTableau;
	int currentRow = 0;
	int currentSlackVariable = 0;
	int totalInequalities,totalSlackVariables,totalAngles;
	int i; //used to iterate through the tableau
	

	splitGraphToVertices(graph,lengthOfGraph,vertexDegrees,vertices);
	retrieveFaces(vertexDegrees,vertices,faceArray,numberOfFaces);
	getAnglesUpTo(vertexDegrees,anglesUpTo);
	
	
	totalAngles = anglesUpTo[(expectedVertices - 1)] + vertexDegrees[(expectedVertices - 1)]; //the number of angles 


	/*MODIFY THIS TO SUIT YOUR OWN SIMPLEX PROBLEM 
	*/
	totalInequalities = totalAngles; //each angle induces an inequality
	totalInequalities += (2 * expectedVertices); //there are 2 inequalities for each vertex: the sum of all angles around a point <= 2pi, >=2pi
	totalInequalities += 3 * (numberOfFaces[2]); //there are 3 inequalities for each triangle : each angle in a triangle is less than or equal to arccos(1/3)
	totalInequalities += 4 * (numberOfFaces[3]); //four inequalities for each quadrilateral : each angle in a quadrilateral is less than or equal to 2arccos(1/3)
	totalInequalities += 4 * (numberOfFaces[3]); //an additional 4 inequalities for each quadrilateral: there are exactly 2 pairs of opposite angles in a quadrilateral and 2 inequalities for each pair
	totalInequalities += 8 * (numberOfFaces[3]); //there are 4 adjacent pairs of angles in a quadrilateral and each pair has exactly two inequalities
	totalInequalities += 5 * (numberOfFaces[4]); //there is one inequality for each angle in a pentagon: there are 5 angles in a pentagon

	totalSlackVariables = totalInequalities; //each inequality will introduce a slack variable


	//DO NOT EDIT THE NEXT FEW LINES

	widthOfTableau = totalAngles + totalSlackVariables;
	heightOfTableau = totalInequalities;

	currentSlackVariable = totalAngles;//start from the final angle
	for(i=0;i<(widthOfTableau + 1);i++)
	{
		tableau[heightOfTableau][i] = 0;
	}
	/*MODIFY THIS TO SUIT YOUR OWN SIMPLEX PROBLEM 
	*/
	addInequality(allInformation,ALL_VERTEX,INDIVIDUAL_ANGLE,0,1.2309f,GEQ);
	addInequality(allInformation,ALL_VERTEX,SUM_ANGLE,0,6.2832f,LEQ);
	addInequality(allInformation,ALL_VERTEX,SUM_ANGLE,0,6.2831f,GEQ);
	addInequality(allInformation,ALL_FACE,INDIVIDUAL_ANGLE,3,1.231f,LEQ);
	addInequality(allInformation,ALL_FACE,INDIVIDUAL_ANGLE,4,2.462f,LEQ);
	addInequality(allInformation,OPPOSITE_FACE,SUBTRACT_ANGLE,4,0,LEQ);
	addInequality(allInformation,ADJACENT_FACE,SUM_ANGLE,4,3.8213f,LEQ);
	addInequality(allInformation,ADJACENT_FACE,SUM_ANGLE,4,3.6928f,GEQ);
	addInequality(allInformation,ALL_FACE,INDIVIDUAL_ANGLE,5, 3.6929f,LEQ);
	if(feasableSimplex(tableau,heightOfTableau,widthOfTableau)) return 1;


	return 0;
}
void convertGraphToASCII(int graph[sizeOfGraph],int lengthOfGraph)
{
	int i; //used to iterate through the graph
	for(i=1;i<lengthOfGraph;i++) //ignore the first character
	{
		if(graph[i]==0) //the end of a vertex
		{
			graph[i] = separatorASCIICode;
		}
		else
		{
			graph[i] += 96;
		}
	}
}
void passGraph(int graph[sizeOfGraph], int length)
{
	int i;
	if(makeTableau(graph,length)) //this graph was permitted
	{
		if(outputAsASCII)
		{
				convertGraphToASCII(graph,length);
		}
		#pragma omp critical(writefile) //avoids OMP race conditions
		{
			if(outputAsASCII)
			{

				printf("%i ",expectedVertices);
			}
			else
			{
				printf("%c",expectedVertices);
			}
			for(i = 1; i<length;i++)
			{
				printf("%c",graph[i]); //print the graph if it is accepted
			}
			if(outputAsASCII)
			{

				printf("\n");
			}
			totalCorrectGraphs++;
		}
	}
}
int main()
{
	int vertexCount; //which vertex are we currently processing
	int graph[sizeOfGraph]; //the information contained within each graph (per thread)
	int lengthOfGraph = 0; //the length of the information contained within each graph (per thread)
	char a; //used to read the file
	int exitLoop = 0;
	int threadsToUse = maxThreads;
	#pragma omp parallel
	{ 
		//this section of code automatically determines the optimal number of threads, constrained by maxThreads
		if(omp_get_thread_num() == 0) //only the first thread should execute this piece of code
		{
			if(threadsToUse>omp_get_num_threads())
			{
				threadsToUse = omp_get_num_threads(); //if the user has asked for too many threads then the number of threads will be reset to the optimal number
			}
			fprintf(stderr,"%i threads will be used to execute this program \n",threadsToUse);
		}
	}
	#pragma omp barrier
	#pragma omp parallel private(graph,lengthOfGraph,a) shared(exitLoop) num_threads(threadsToUse) //each thread has it's own graph, length of graph and input character(a). The flag to exit the loop is shared between all threads
	{
		a = 0;
		while(!exitLoop)
		{
			//if exitloop is set then all threads will eventually exit this loop
			#pragma omp critical(readGraph) //reading a graph needs to be done in serial, not parallel.
			{
				//a thread may enter this section of code even though there are no graphs left to process. Consequently we need to check whether or not another thread has set exitloop
				lengthOfGraph = 0;
				if(!exitLoop) //we have not necessarily reached the end of file
				{

					a = getchar(); //read the first character
					if(expectedVertices==0) {expectedVertices = a;}
					if(a!=expectedVertices)
					{
						if(a!=EOF)
						{

							fprintf(stderr,"A graph started with %i vertices when %i vertices were expected\n",a,expectedVertices);
						}
						else
						{
							exitLoop = 1; //now we have reached the end of the file. Signal to all other threads to stop processing new graphs
						}

					}
					else
					{
						//read in the graph: there are some vertices left
						graph[0] = a;
						for(vertexCount=0;vertexCount<expectedVertices;vertexCount++)
						{
							a = getchar();
							while(a!=0)
							{
								lengthOfGraph++;
								graph[lengthOfGraph] = a; //write this edge to the current graph

								a=getchar();
							}

							lengthOfGraph++;
							graph[lengthOfGraph]=a;
						}
					}
				}
			}
			if(lengthOfGraph != 0) //we have no graph to process
			{

				passGraph(graph,(lengthOfGraph + 1)); //test this graph
				#pragma omp atomic
				totalProcessedGraphs++; //this operation also needs to be done in serial 
				lengthOfGraph = 0;
			}
		}
	}
	fprintf(stderr,"%i graphs were processed, and %i graphs remain\n",totalProcessedGraphs,totalCorrectGraphs);

}