/* This is the main neural part of the project. This program can be
	 used to both train and execute the network.

	 To train, specify the following on the command line:
	 learning rate (should be VERY low; ~1E-9),
	 stopping MSE value,
	 weights file to load (if not specified, random weights are used).

	 Weights are saved to weights.dat only when the network finishes
	 training.

	 To execute, specify the following on the command line:
	 Bmp filename (64x64 uncompressed 24-bit bmp. Don't specify ".bmp"),
	 name of the weights file


	 Overall, there are many things in this code that aren't
	 pretty. This was meant to be written fast.
*/

#include <stdlib.h>
#include <time.h>
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <ctype.h>
#include <string.h>

#define ANGLESPER90							45
#define WEIGHTARRAYSIDELENGTH		(ANGLESPER90 + 1)
#define TANHSTRETCH							30

#include "honn.h"

typedef unsigned short int16;
typedef unsigned long  int32;

struct SBMPFileHeader
{
	int16					bfType;
	int32					bfSize;
	int32					bfReserved;
	int32					bfDataOffset;
} __attribute__ ((packed));

struct SBMPInfoHeader
{
	int32					biSize;
	int32					biWidth;
	int32					biHeight;
	int16					biPlanes;
	int16					biBitCount;
	int32					biCompression;
	int32					biSizeImage;
	int32					biXPelsPerMeter;
	int32					biYPelsPerMeter;
	int32					biClrUsed;
	int32					biClrImportant;
} __attribute__ ((packed));

inline double& max(double& a1, double& a2)
{
	return (a1 > a2) ? a1 : a2;
}

inline unsigned short round(double x)
{
	unsigned short out = (unsigned short)x;

	return (x-out < 0.5) ? out : out+1;
}

CDataVector::CDataVector(unsigned short numInputs)
	: m_pNext(NULL)
{
	m_pInputs = new char[numInputs];
}

CDataVector::~CDataVector(void)
{
	if(m_pInputs != NULL)
		delete m_pInputs;

	if(m_pNext != NULL)
		delete m_pNext;
}

CDataVector* CDataVector::Add(unsigned short numInputs)
{
	if(m_pNext != NULL)
		return m_pNext->Add(numInputs);

	m_pNext = new CDataVector(numInputs);
}

CDataVector* CDataVector::Add(char *pBitmapFile, char output)
{
	if(m_pNext != NULL)
		return m_pNext->Add(pBitmapFile, output);


	SBMPFileHeader BMPFileHeader;
	SBMPInfoHeader BMPInfoHeader;

	ifstream input(pBitmapFile);

	if(!input.is_open())
	{
		cerr << "Couldn't open input file";
		return NULL;
	}

	input.read(&BMPFileHeader, sizeof(BMPFileHeader));
	input.read(&BMPInfoHeader, sizeof(BMPInfoHeader));
	BMPInfoHeader.biSizeImage = BMPFileHeader.bfSize - BMPFileHeader.bfDataOffset;

	if(BMPInfoHeader.biBitCount != 24)
	{
		cerr << "Bit depths other than 24-bit aren't supported\n";
		return NULL;
	}

	input.seekg(BMPFileHeader.bfDataOffset);

	m_pNext = new CDataVector(BMPInfoHeader.biWidth*BMPInfoHeader.biHeight);

	for(unsigned long i=0; i < BMPInfoHeader.biWidth*BMPInfoHeader.biHeight; i++)
	{
		input.read(m_pNext->m_pInputs + i, 1);
		input.seekg(2, ios::cur);

		if(m_pNext->m_pInputs[i] == 0)
			m_pNext->m_pInputs[i] = '\x01';
		else
			m_pNext->m_pInputs[i] = '\x00';
	}

	//	m_pNext->m_width = BMPInfoHeader.biWidth;
	m_pNext->m_output = output;
}

CDataVector* CDataVector::getNthVector(unsigned short n)
{
	CDataVector *pVec = m_pNext;

	while(n--)
		pVec = pVec->m_pNext;

	return pVec;
}

CNode::CNode(unsigned short width,
						 double fracTrainingVectors,
						 char *pWeightsFile)
	:	m_width(width),
		m_numInputs(width * width),
		m_fracTrainingVectors(fracTrainingVectors),
		m_numEqConnections(WEIGHTARRAYSIDELENGTH*WEIGHTARRAYSIDELENGTH),
		m_numTuples( (width * width) * ((width * width)-1) * ((width * width)-2) / 6),
		m_numVectors(0),
		m_bConverged(false)
{
	m_pEqConnections = new SEqConnection[m_numEqConnections];
	m_pTuplePointers = new unsigned short[m_numTuples];

	if(pWeightsFile == NULL)
	{
		// Initialize the weights to a random double
		for(unsigned short weightIndex = 0; weightIndex < m_numEqConnections; weightIndex++)
		{
			m_pEqConnections[weightIndex].m_weight = 0.01 * (((double)rand() / RAND_MAX) - 0.5);
		}
	}
	else
	{
		ifstream weights(pWeightsFile);
		for(unsigned short weightIndex = 0;
				weightIndex < m_numEqConnections;
				weightIndex++)
			{
				weights.read(&m_pEqConnections[weightIndex].m_weight,
										 sizeof(m_pEqConnections[weightIndex].m_weight));
			}
	}

	// Initialize the tuple pointer array. Do this by cycling through the 
	unsigned long							tupleIndex;
	unsigned short						p1,p2,p3;
	unsigned char							x1, x2, x3, y1, y2, y3;
	double										length1, length2, length3;
	double										theta1, theta2;
	double										temp;
	unsigned short						weightIndex;

	ifstream mappingFile("mapping.dat");

	if(mappingFile.is_open())
	{
		mappingFile.read(m_pTuplePointers, m_numTuples * sizeof(*m_pTuplePointers));
	}
	else
	{
		tupleIndex = 0;

		for(p1=0; p1<m_numInputs; p1++)
		{
			y1 = p1 / m_width;
			x1 = p1 % m_width;

			for(p2=p1+1; p2<m_numInputs; p2++)
			{
				y2 = p2 / m_width;
				x2 = p2 % m_width;

				for(p3=p2+1; p3<m_numInputs; p3++)
				{
					y3 = p3 / m_width;
					x3 = p3 % m_width;

					length1 = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
					length2 = sqrt((x3-x1)*(x3-x1) + (y3-y1)*(y3-y1));
					length3 = sqrt((x3-x2)*(x3-x2) + (y3-y2)*(y3-y2));

					double& biggestLength = max(length1, max(length2, length3));
					if(biggestLength != length3)
					{
						temp = length3;
						length3 = biggestLength;
						biggestLength = temp;
					}
					if(length1 > length2)
					{
						temp = length1;
						length1 = length2;
						length2 = temp;
					}

					// some of the args are meant to streamline the code
					theta1 = acos((length2*length2 + length3*length3 - length1*length1)
												/ (2.0 * length2 *length3));

					theta2 = acos((length3*length3 + length1*length1 - length2*length2)
												/ (2.0 * length3 *length1));

					theta1 = theta1 * 2.0 * ANGLESPER90 / M_PI;
					theta2 = theta2 * 2.0 * ANGLESPER90 / M_PI;

					weightIndex = round(theta2) * WEIGHTARRAYSIDELENGTH + round(theta1);

					m_pTuplePointers[tupleIndex++] = weightIndex;
				}
			}
		}

		ofstream out("mapping.dat");
		out.write(m_pTuplePointers, m_numTuples * sizeof(*m_pTuplePointers));
	}
}

CNode::~CNode(void)
{
	delete[] m_pEqConnections;
	delete[] m_pTuplePointers;;
}

bool CNode::ReadDataVectors(char *pFileName)
{
	// EXECUTION
	if(pFileName)
	{
		char strstart[100] = "";

		strcat(strstart, pFileName);

		for(char i='0'; i<'4'; i++)
		{
			char strall[100] = "";
			char strnum[2] = "1";

			strnum[0] = i;
			strcat(strall, strstart);
			strcat(strall, strnum);
			strcat(strall, ".bmp");

			m_dataVectors.Add(strall,'\x01');
			m_numVectors++;
		}

	}
	else
	{
/*
		m_dataVectors.Add("bigsqcoarse0.bmp",'\xff');
		m_dataVectors.Add("bigsqcoarse1.bmp",'\xff');
		m_dataVectors.Add("bigsqcoarse2.bmp",'\xff');
		m_dataVectors.Add("bigsqcoarse3.bmp",'\xff');

		m_dataVectors.Add("bigtricoarse0.bmp",'\x01');
		m_dataVectors.Add("bigtricoarse1.bmp",'\x01');
		m_dataVectors.Add("bigtricoarse2.bmp",'\x01');
		m_dataVectors.Add("bigtricoarse3.bmp",'\x01');
*/
		m_dataVectors.Add("black.bmp",'\xff');
		m_dataVectors.Add("black.bmp",'\xff');
		m_dataVectors.Add("black.bmp",'\xff');
		m_dataVectors.Add("black.bmp",'\xff');

		m_dataVectors.Add("pat1coarse0.bmp",'\x01');
		m_dataVectors.Add("pat1coarse1.bmp",'\x01');
		m_dataVectors.Add("pat1coarse2.bmp",'\x01');
		m_dataVectors.Add("pat1coarse3.bmp",'\x01');
/*
		m_dataVectors.Add("pat3coarse0.bmp",'\x01');
		m_dataVectors.Add("pat3coarse1.bmp",'\x01');
		m_dataVectors.Add("pat3coarse2.bmp",'\x01');
		m_dataVectors.Add("pat3coarse3.bmp",'\x01');
*/

		m_numVectors += 8;
	}


/*	while(cin.peek() != EOF)
	{
		if(cin.peek() == '\n')
		{
			cin.ignore(1);
			continue;
		}

		CDataVector *pData = m_dataVectors.Add(m_numInputs);

		cin >> pData->m_output;

		for(unsigned short inputIndex = 0; inputIndex < m_numInputs; inputIndex++)
		{
			cin >> pData->m_pInputs[inputIndex];
		}

		m_numVectors++;
	}

	// Return with an error condition if no input vectors were specified
	if(m_numVectors == 0)
		return false;
*/

	m_lastTrainingVector = (unsigned short)(m_fracTrainingVectors * m_numVectors - 1);
	m_pLastTrainingVector = m_dataVectors.getNthVector(m_lastTrainingVector);

	return true;
}

bool CNode::ProcessEpoch(double MSEgoal)
{
	CDataVector *pData = m_dataVectors.m_pNext;

	m_MSE = 0;

	for(unsigned short vectorIndex = 0;
			vectorIndex <= m_lastTrainingVector;
			vectorIndex += 4)
	{
		double error = ProcessVector(pData);
		m_MSE += error * error;

		pData = pData->m_pNext->m_pNext->m_pNext->m_pNext;
	}

	m_MSE /= (m_lastTrainingVector + 1)/4;
	if(m_MSE <= MSEgoal)
	{
		m_bConverged = true;
		return false;
	}

	cout << m_MSE << '\n';

	return true;
}

double CNode::run(CDataVector *pData)
{
	unsigned long							tupleIndex;
	unsigned short						p1,p2,p3;
	double										sum = 0;

	if(pData == NULL)
	{
		pData = m_dataVectors.m_pNext;
	}
 
	unsigned char i;

	// cycle through the coarse-coded images
	for(i=0; i<4; i++)
	{
		tupleIndex = 0;
		for(p1=0; p1<m_numInputs; p1++)
		{
			for(p2=p1+1; p2<m_numInputs; p2++)
			{
				for(p3=p2+1; p3<m_numInputs; p3++)
				{
					//				if(1)
					if(m_pTuplePointers[tupleIndex] != 0)
					{
						sum += m_pEqConnections[m_pTuplePointers[tupleIndex]].m_weight * 
							pData->m_pInputs[p1] *
							pData->m_pInputs[p2] *
							pData->m_pInputs[p3];
					}

					tupleIndex++;
				}
			}
		}

		pData = pData->m_pNext;
	}

	return tanh(sum / TANHSTRETCH);
}

double CNode::ProcessVector(CDataVector *pData)
{
	double										output;
	double										error;

	double										correction;

	unsigned long							tupleIndex;
	unsigned short						p1,p2,p3;

	output = run(pData);

	error = pData->m_output - output;

	correction =
		CNode::learningRate *
		error *
		(1 - output*output);

	unsigned char i;

	// cycle through the coarse-coded images
	for(i=0; i<4; i++)
	{
		tupleIndex = 0;
		for(p1=0; p1<m_numInputs; p1++)
		{
			for(p2=p1+1; p2<m_numInputs; p2++)
			{
				for(p3=p2+1; p3<m_numInputs; p3++)
				{
					m_pEqConnections[m_pTuplePointers[tupleIndex]].m_weight +=
						correction *
						pData->m_pInputs[p1] *
						pData->m_pInputs[p2] *
						pData->m_pInputs[p3];

					tupleIndex++;
				}
			}
		}

		pData = pData->m_pNext;
	}

	return error;
}

void CNode::print(void)
{
	ofstream out ("weights.dat");

	for(unsigned short weightIndex = 0;
			weightIndex < m_numEqConnections;
			weightIndex++)
	{
		out.write(&m_pEqConnections[weightIndex].m_weight,
							sizeof(m_pEqConnections[weightIndex].m_weight));
	}
}

int main(int argc, char *argv[])
{
	CNode::learningRate = 0.1;
	double								MSEgoal = 0.01;
	unsigned long					maxEpochs = 1000;
	double								fracTrainingVectors = 1.0;

	unsigned short				numInputs, numOutputs, width;
	unsigned long					epochCount = 0;

	char *pWeightsFile = "weights.dat\0                                      ";

	numOutputs = 1;
	numInputs = 256;
	width = 16;

	if(numOutputs != 1)
	{
		cerr << "Multiple outputs not supported\n";
		return 1;
	}

	srand(time(NULL));

	if(argc > 1)
	{
		CNode::learningRate = atof(argv[1]);
		if(argc > 2)
		{
			MSEgoal = atof(argv[2]);
			if(argc > 3)
			{
				pWeightsFile = NULL;
			}
		}
	}

	if(isalpha(argv[1][0]))
	{
		CNode					node(width, fracTrainingVectors,argv[2]);

		node.ReadDataVectors(argv[1]);

		double output = node.run();

		cout << output << '\n';

		return 0;
	}

	CNode					node(width, fracTrainingVectors, pWeightsFile);

	if(!node.ReadDataVectors(NULL))
	{
		cerr << "No data vectors specified.\n";
		return 1;
	}

	while(node.ProcessEpoch(MSEgoal))
	{
		if(++epochCount == maxEpochs)
		{
			cout << "No convergence in " << epochCount << " epochs.\n";
			break;
		}
	}

	if(node.converged())
		cout << "Converged in " << epochCount << " epochs.\n";

	node.print();

	return 0;
}