float Neuron::Calc(float *inputs)
{
	float out = 0;

	//add up all (weight * input) pairs
	for (int x=0;x<num_inputs;x++)
		out += weights[x] * inputs[x];

	//sigmoid reponse curve (broken)
	//return ( 1 / (sig_response + (float)exp( - (out / 10) )) );

	//sinus response curve
	return sinf( out / sig_response);
} 

/*=======================================

	cMyError::Header File
	
	By Matt Barnett;

	Description: Simple error reporting.

	Usage:

		cMyError cMyError1;

		//quiet error printf()
		cMyError1.Error("vertex number %i is %s",vertex,stringproblem);

		//message box and printf()
		//same as FatalError(); except it doesnt look so intimidating
		cMyError1.AlertError("feature not supported %s",stringfeature);

		//call this to report a fatal error
		//just before you exit(-1);
		cMyError1.FatalError("Couldnt find file %s",stringlostfile);
		exit(-1);
		

=======================================*/

#ifndef		_MYERROR_CLASS_H
#define		_MYERROR_CLASS_H


#include	<windows.h>							//comment for non Win32 platforms


#include	<cstdio>							//for printf

//#define VERBOSE_CREATION						//class debug


class cMyError
{
	public:							
			
		cMyError();											//constructor
		~cMyError();										//destructon

		void FatalError(const char *strString, ...);		//MessageBox and printf()
		void AlertError(const char *strString, ...);		//exactly the same as FatalError (except title)
		void Error(const char *strString, ...);				//printf()		
		void Alert(const char *strString, ...);
		void Log(const char *strString, ...);

		static void isFlowHere(const char *filename, int linenumber)
		{
			printf("Program flow reached %s %i\n",filename,linenumber);
		}

		//static bool logall;
};

//macros are shouted and they use a ; so the usage is ...
//
//	CMYERROR_ISFLOWHERE
//  and not CMYERROR_ISFLOWHERE;
//
//which helps make it clear that its a macro we're calling
//

#define CMYERROR_ISFLOWHERE cMyError::isFlowHere(__FILE__,__LINE__);

#endif 

/*=======================================

	cMyError::Implementation

	*See header file for more information
	
	By Matt Barnett;

=======================================*/

#include "myerror.h"
#include "leakdetect.h"

/*=======================================
cMyError::Constructor
=======================================*/

cMyError::cMyError()
{
#ifdef VERBOSE_CREATION
	printf("cMyError::Constructor::Empty\n");
#endif
}

/*=======================================
cMyError::Destructor
=======================================*/

cMyError::~cMyError()
{

#ifdef VERBOSE_CREATION
	printf("cMyError::Destructor::Empty\n");
#endif

}

/*=========================================================
void cMyError::FatalError()
==========================================================*/

void cMyError::FatalError(const char *strString, ...)
{
	char		strText[256];				// This will hold our text to display
	va_list		argumentPtr;				// This will hold the pointer to the argument list
	float		unitLength=0.0f;			// This will store the length of the 3D Font in unit length
		
	if (strString == NULL)					// Check if a string was given
		return;								// Don't render anything then

	va_start(argumentPtr, strString);			// Parse the arguments out of the string
	vsprintf(strText, strString, argumentPtr);	// Now add the arguments into the full string
	va_end(argumentPtr);						// This resets and frees the pointer to the argument list.

	MessageBox(NULL,strText,"Fatal Error!!!",MB_OK|MB_ICONSTOP);
	printf("FATAL ERROR: %s\n",strText);
}

/*=========================================================
void cMyError::AlertError()
==========================================================*/

void cMyError::AlertError(const char *strString, ...)
{
	char		strText[256];				// This will hold our text to display
	va_list		argumentPtr;				// This will hold the pointer to the argument list
	float		unitLength=0.0f;			// This will store the length of the 3D Font in unit length
		
	if (strString == NULL)					// Check if a string was given
		return;								// Don't render anything then

	va_start(argumentPtr, strString);			// Parse the arguments out of the string
	vsprintf(strText, strString, argumentPtr);	// Now add the arguments into the full string
	va_end(argumentPtr);						// This resets and frees the pointer to the argument list.

	MessageBox(NULL,strText,"Non fatal error!",MB_OK|MB_ICONEXCLAMATION);
	printf("ERROR(loud): %s\n",strText);
}

/*=========================================================
void cMyError::Error()
==========================================================*/

void cMyError::Error(const char *strString, ...)
{
	char		strText[256];				// This will hold our text to display
	va_list		argumentPtr;				// This will hold the pointer to the argument list
	float		unitLength=0.0f;			// This will store the length of the 3D Font in unit length
		
	if (strString == NULL)					// Check if a string was given
		return;								// Don't render anything then

	va_start(argumentPtr, strString);			// Parse the arguments out of the string
	vsprintf(strText, strString, argumentPtr);	// Now add the arguments into the full string
	va_end(argumentPtr);						// This resets and frees the pointer to the argument list.

	printf("ERROR(quiet): %s\n",strText);
}

/*=========================================================
void cMyError::Alert()
==========================================================*/

void cMyError::Alert(const char *strString, ...)
{
	char		strText[256];				// This will hold our text to display
	va_list		argumentPtr;				// This will hold the pointer to the argument list
	float		unitLength=0.0f;			// This will store the length of the 3D Font in unit length
		
	if (strString == NULL)					// Check if a string was given
		return;								// Don't render anything then

	va_start(argumentPtr, strString);			// Parse the arguments out of the string
	vsprintf(strText, strString, argumentPtr);	// Now add the arguments into the full string
	va_end(argumentPtr);						// This resets and frees the pointer to the argument list.

	MessageBox(NULL,strText,"Information",MB_OK);
	printf("Alert(loud): %s\n",strText);
}
/*=========================================================
void cMyError::Log()
==========================================================*/

void cMyError::Log(const char *strString, ...)
{
	char		strText[256];				// This will hold our text to display
	va_list		argumentPtr;				// This will hold the pointer to the argument list
	float		unitLength=0.0f;			// This will store the length of the 3D Font in unit length
		
	if (strString == NULL)					// Check if a string was given
		return;								// Don't render anything then

	va_start(argumentPtr, strString);			// Parse the arguments out of the string
	vsprintf(strText, strString, argumentPtr);	// Now add the arguments into the full string
	va_end(argumentPtr);						// This resets and frees the pointer to the argument list.

	printf("Log: %s\n",strText);
} 

/*=======================================

	Memory leak detection header file

=======================================*/

//#define LEAKDBG

#ifdef LEAKDBG
	#define _CRTDBG_MAP_ALLOC
		#include <stdlib.h>
		#include <crtdbg.h>
#endif

#ifdef LEAKDBG
  #define meMalloc(s)       _malloc_dbg(s, _NORMAL_BLOCK, __FILE__, __LINE__)
  #define meCalloc(c, s)    _calloc_dbg(c, s, _NORMAL_BLOCK, __FILE__, __LINE__)
  #define meRealloc(p, s)   _realloc_dbg(p, s, _NORMAL_BLOCK, __FILE__, __LINE__)
  #define meExpand(p, s)    _expand_dbg(p, s, _NORMAL_BLOCK, __FILE__, __LINE__)
  #define meFree(p)         _free_dbg(p, _NORMAL_BLOCK)
  #define meMemSize(p)      _msize_dbg(p, _NORMAL_BLOCK)

  #define meNew new(_NORMAL_BLOCK, __FILE__, __LINE__)
  #define meDelete delete

  // Set to dump leaks at the program exit.
  #define meInitMemoryCheck() \
    _CrtSetDbgFlag( _CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF)

#else
  #define meMalloc malloc
  #define meCalloc calloc
  #define meRealloc realloc
  #define meExpand _expand
  #define meFree free
  #define meMemSize _msize

  #define meNew new
  #define meDelete delete

  #define meInitMemoryCheck()

#endif 

/*---------------------------------------------
		
	MatNN library : Test App 1 (cpp)

		* This is not an terrain poly pushing demo.

		This shows a silly brainwave of the neural
		net. The nets are totally untrained (random weights).

		The neural net (in this example) has 3 input neurons
		that take the x / y position of a 2d grid and the time/ticks
		returned by GetTickCount(); It has 3 outputs = RGB color which
		is also used as a height.
		 
		We/me/i graph the response of the the net over a small 2d plane
		of x/y values, at a specific time (using GetTickCount() ) to create
		a pretty swirly jaggy heightmap thingy (tm)

		keys : r = random neural net
			   a = decrease response
			   d = increase response

	By Matt Barnett

		(see References)	* See matnn.h for more information
	
---------------------------------------------*/

//basic glut and math stuff for camera etc
 
#include <ctime>
#include <cstdio>
#include <windows.h>
#include <GL/glut.h>

#include "../../../imath/imath.h"
#include "../../../matnn/matnn.h"

GLfloat grey[] = {0.2,0.2,0.2,1.0};
GLfloat above_origin[] = {0.0,0.0,3.0,1.0};
GLfloat white[] = {1.0,1.0,1.0,1.0};
GLfloat	shiny[] = {50.0};

#define WINDOW_SIZE			1024		//
#define	SCALE_MOTION		100			//scale down by

int x,y,mx,my,xs,ys;


int oldmx = 0;
int oldmy = 0;

enum modetype { MODE_ROTATE,MODE_MOVE,MODE_CAM };
modetype mode;
iquaternion4 camquat;

//globals for testing here

using MatLayeredNeuralNetLib::NeuralNet;
float resp = 0.45f;
float inputs[3];
float outputs[3];


//change the number of hidden layers
//(number of inputs, number of outputs, number of neurons in hidden layers, number of hidden layers)
NeuralNet n1(3,3,4,3,inputs,outputs);

 
//------------------------
	
void drawaxes()
{
  glDisable(GL_LIGHTING);

  glBegin(GL_LINES);
  glColor3f(1.0,0.0,0.0);glVertex3f(-0.5,0.0,0.0);glVertex3f(1.2,0.0,0.0);
  glColor3f(0.0,1.0,0.0);glVertex3f(0.0,-0.5,0.0);glVertex3f(0.0,1.2,0.0);
  glColor3f(0.0,0.0,1.0);glVertex3f(0.0,0.0,-0.5);glVertex3f(0.0,0.0,1.2);
  glEnd();
}

void display(void)
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	imatrix16 rm(camquat);
	glTranslatef(0.0,0.0,-0.8);
	glMultMatrixf(rm.m);
	
	//instert drawing here
	
	float floorsize = 1.5;

	//if you want more tri-strips then decrease this value
	float floorstep = 0.025f;

	float i,j;

	float fftime = (float)GetTickCount() / 1500;

	glPushMatrix();
	glTranslatef((float)-floorsize/2,(float)-floorsize/2,0.0);

	for (i=0;i<floorsize;i+=floorstep)
	{
		glBegin(GL_TRIANGLE_STRIP);
		for (j=0;j<floorsize;j+=floorstep)
		{
			
				inputs[0] = j;
				inputs[1] = i;
				inputs[2] = fftime;	
				n1.RunNN();

				glColor3f(outputs[0],outputs[1],outputs[2]);			//rem this out and try one of the below
				//glVertex3f(i,j,outputs[2] / 9.0 )				
				//glVertex3f(i,j,0.0f);
				glVertex3f(i,j,(outputs[0]+outputs[1]+outputs[2]) / 35.0f);

				inputs[0] = j+floorstep;
				inputs[1] = i+floorstep;
				inputs[2] = fftime;	
				n1.RunNN();

				glColor3f(outputs[0],outputs[1],outputs[2]);			//rem this out and try one of the below (chose the same as above)
				//glVertex3f(i+floorstep,j+floorstep,outputs[2] / 9.0 );
				//glVertex3f(i+floorstep,j+floorstep,0.0f);
				glVertex3f(i+floorstep,j+floorstep,(outputs[0]+outputs[1]+outputs[2])  / 35.0f);

		}
		glEnd();
	}


	drawaxes();
	glPopMatrix();


	glPopMatrix();
	glutSwapBuffers();
	glutPostRedisplay();
}



void reshape(int w, int h)
{
  glViewport(0,0,(GLsizei) w,(GLsizei) h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(90,(float)w/h,0.01,100.0);
}



void motion(int x, int y)
{
	mx = x;my = y;

	xs = mx - oldmx;
	ys = oldmy - my;

	oldmx = mx;
	oldmy = my;

	iquaternion4 multquat;

	if (mode == MODE_ROTATE)
	{	
	
	}
	else if (mode == MODE_MOVE)
	{

	}
	else if (mode == MODE_CAM)
	{
		multquat.makeEulerRotation( ivector3( (float)ys/SCALE_MOTION,-(float)xs/SCALE_MOTION,0.0 ) );
		camquat.multiply(multquat);
	}
}


void mouse(int button,int state,int x,int y)
{
	oldmx = x;
	oldmy = y;

	if (button == GLUT_LEFT_BUTTON)
		mode = MODE_ROTATE;

	if (button == GLUT_RIGHT_BUTTON)
		mode = MODE_MOVE;

	
	if (button == GLUT_MIDDLE_BUTTON)
		mode = MODE_CAM;
}

void init(void)
{

	srand( (unsigned)time( NULL ) );	


	glEnable(GL_DEPTH_TEST);
	//glEnable(GL_LIGHTING);
	//glEnable(GL_LIGHT0);
	//glLightfv(GL_LIGHT0,GL_POSITION,above_origin);
	
	//glEnable(GL_BLEND);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	//glMaterialfv(GL_FRONT,GL_DIFFUSE,grey);
	//glMaterialfv(GL_FRONT,GL_SPECULAR,white);
	//glMaterialfv(GL_FRONT,GL_SHININESS,shiny);

	//glEnable(GL_POLYGON_OFFSET_FILL);
	//glPolygonOffset(0.1,0.1);
	//glDepthFunc(GL_LEQUAL);
	//glEnable(GL_DEPTH_TEST);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
	//glEnable(GL_BLEND);
	//glEnable(GL_POLYGON_SMOOTH);
	//glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	glPointSize(2.0f);
	

	inputs[0] = 0;
	inputs[1] = 0;

	n1.SetResponse(resp);
}

void key(GLubyte key, int x, int y)
{
  switch (key)
    {
			case 'a':resp -= 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);
			break;

			case 'd':resp += 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);
			
			break;

			case 'r':im_random_seed_time();n1.Randomize();

			break;

			case '1':

			break;

			case '!':

			break;


			case 'l':

			break;


			case 'q':exit(0);
			break;
    }
	
}

int main(int argc, char **argv)
{

	meInitMemoryCheck();
	
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
	glutInitWindowSize(WINDOW_SIZE,WINDOW_SIZE);
	glutInitWindowPosition(100,10);
	glutCreateWindow("Matts test");

	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(key);
	glutMotionFunc(motion);
	glutMouseFunc(mouse);

	init();
	glutMainLoop();

	return 0;
} 

/*---------------------------------------------
		
	MatNN library : Test App 2 (cpp)


		This one uses a neural net to control the position and color
		of some opengl points (large pointsize).

		I guess this could be used to control a bunch of blended particles
		in an interesting way.


		keys : r = random neural net
			   a = decrease response
			   d = increase response

	By Matt Barnett

		(see References)	* See matnn.h for more information
	
---------------------------------------------*/

//basic glut and math stuff for camera etc

#include <ctime>
#include <cstdio>
#include <windows.h>
#include <GL/glut.h>

#include "../../../imath/imath.h"
#include "../../../matnn/matnn.h"

GLfloat grey[] = {0.2,0.2,0.2,1.0};
GLfloat above_origin[] = {0.0,0.0,3.0,1.0};
GLfloat white[] = {1.0,1.0,1.0,1.0};
GLfloat	shiny[] = {50.0};

#define WINDOW_SIZE			1024		//
#define	SCALE_MOTION		100			//scale down by

int x,y,mx,my,xs,ys;
int oldmx = 0;
int oldmy = 0;

enum modetype { MODE_ROTATE,MODE_MOVE,MODE_CAM };
modetype mode;
iquaternion4 camquat;

//globals for testing here

using MatLayeredNeuralNetLib::NeuralNet;
float resp = 1.0f;
float inputs[6];
float outputs[6];
NeuralNet n1(6,6,6,2,inputs,outputs);


//------------------------
	
void drawaxes()
{
  glDisable(GL_LIGHTING);

  glBegin(GL_LINES);
  glColor3f(1.0,0.0,0.0);glVertex3f(-0.5,0.0,0.0);glVertex3f(1.2,0.0,0.0);
  glColor3f(0.0,1.0,0.0);glVertex3f(0.0,-0.5,0.0);glVertex3f(0.0,1.2,0.0);
  glColor3f(0.0,0.0,1.0);glVertex3f(0.0,0.0,-0.5);glVertex3f(0.0,0.0,1.2);
  glEnd();
}

void display(void)
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	imatrix16 rm(camquat);
	glTranslatef(0.0,0.0,-1.7);
	glMultMatrixf(rm.m);
	
	//instert drawing here
	
	float floorsize = 1.5;
	float floorstep = 0.25f;
	float i,j,k,l,m;

	float fftime = (float)GetTickCount() / 1000;
	//float bp1 = 0.0;
	//float bp2 = 1.0;
	//float bp3 = 0.0;

	glPushMatrix();
	glTranslatef((float)-floorsize/2,(float)-floorsize/2,0.0);
	glBegin(GL_POINTS);
	for (i=0;i<floorsize;i+=floorstep)
	{
		for (j=0;j<floorsize;j+=floorstep)
		{	
			for (k=0;k<floorsize;k+=floorstep)
			{	
				for (l=0;l<floorsize;l+=floorstep)
				{
				//	for (m=0;m<floorsize;m+=floorstep)
				//	{
						inputs[0] = j;
						inputs[1] = i;
						inputs[2] = k;	
						inputs[3] = l;
						inputs[4] = 0.0;
						inputs[5] = fftime;
						n1.RunNN();

						glColor3f(outputs[0],outputs[1],outputs[2]);
						glVertex3f(outputs[3],outputs[4],outputs[5]);
				//	}
				}
			}
		}
	
	}
	glEnd();

	drawaxes();
	glPopMatrix();


	glPopMatrix();
	glutSwapBuffers();
	glutPostRedisplay();
}



void reshape(int w, int h)
{
  glViewport(0,0,(GLsizei) w,(GLsizei) h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(90,(float)w/h,0.01,100.0);
}



void motion(int x, int y)
{
	mx = x;my = y;

	xs = mx - oldmx;
	ys = oldmy - my;

	oldmx = mx;
	oldmy = my;

	iquaternion4 multquat;

	if (mode == MODE_ROTATE)
	{	
	
	}
	else if (mode == MODE_MOVE)
	{

	}
	else if (mode == MODE_CAM)
	{
		multquat.makeEulerRotation( ivector3( (float)ys/SCALE_MOTION,-(float)xs/SCALE_MOTION,0.0 ) );
		camquat.multiply(multquat);
	}
}


void mouse(int button,int state,int x,int y)
{
	oldmx = x;
	oldmy = y;

	if (button == GLUT_LEFT_BUTTON)
		mode = MODE_ROTATE;

	if (button == GLUT_RIGHT_BUTTON)
		mode = MODE_MOVE;

	
	if (button == GLUT_MIDDLE_BUTTON)
		mode = MODE_CAM;
}

void init(void)
{

	srand( (unsigned)time( NULL ) );	


	glEnable(GL_DEPTH_TEST);
	//glEnable(GL_LIGHTING);
	//glEnable(GL_LIGHT0);
	//glLightfv(GL_LIGHT0,GL_POSITION,above_origin);
	
	//glEnable(GL_BLEND);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	//glMaterialfv(GL_FRONT,GL_DIFFUSE,grey);
	//glMaterialfv(GL_FRONT,GL_SPECULAR,white);
	//glMaterialfv(GL_FRONT,GL_SHININESS,shiny);

	//glEnable(GL_POLYGON_OFFSET_FILL);
	//glPolygonOffset(0.1,0.1);
	//glDepthFunc(GL_LEQUAL);
	//glEnable(GL_DEPTH_TEST);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
	//glEnable(GL_BLEND);
	//glEnable(GL_POLYGON_SMOOTH);
	//glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	glPointSize(15.0f);
	printf("%f",log(100) );

	inputs[0] = 0;
	inputs[1] = 0;

	n1.SetResponse(resp);
}

void key(GLubyte key, int x, int y)
{
  switch (key)
    {
			case 'a':resp -= 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);
			break;

			case 'd':resp += 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);
			
			break;

			case 'r':im_random_seed_time();n1.Randomize();

			break;

			case '1':

			break;

			case '!':

			break;


			case 'l':

			break;


			case 'q':exit(0);
			break;
    }
	
}

int main(int argc, char **argv)
{
	meInitMemoryCheck();

	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
	glutInitWindowSize(WINDOW_SIZE,WINDOW_SIZE);
	glutInitWindowPosition(100,10);
	glutCreateWindow("Matts test");

	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(key);
	glutMotionFunc(motion);
	glutMouseFunc(mouse);

	init();
	glutMainLoop();

	return 0;
} 

/*---------------------------------------------
		
	MatNN library : Test App 4 (cpp)

		This shows the copy constructor working.
		(4 different instances of a net are drawn,
		that were created by copy constructors)

		It will show the assignment operator
		working aswell, when/if i decide its right
		to assign neural nets.

		keys : r = random neural nets
			   a = decrease response
			   d = increase response

	By Matt Barnett

		(see References)	* See matnn.h for more information
	
---------------------------------------------*/

//basic glut and math stuff for camera etc

#include <ctime>
#include <cstdio>
#include <windows.h>
#include <GL/glut.h>

#include "../../../imath/imath.h"
#include "../../../matnn/matnn.h"

GLfloat grey[] = {0.2,0.2,0.2,1.0};
GLfloat above_origin[] = {0.0,0.0,3.0,1.0};
GLfloat white[] = {1.0,1.0,1.0,1.0};
GLfloat	shiny[] = {50.0};

#define WINDOW_SIZE			1024		//
#define	SCALE_MOTION		100			//scale down by

int x,y,mx,my,xs,ys;
int oldmx = 0;
int oldmy = 0;

enum modetype { MODE_ROTATE,MODE_MOVE,MODE_CAM };
modetype mode;
iquaternion4 camquat;

//globals for testing here

using MatLayeredNeuralNetLib::NeuralNet;
float resp = 1.0f;
float inputs[6];
float outputs[3];

NeuralNet *n0 = meNew NeuralNet(6,3,4,5,inputs,outputs);	//on the heap, so that we can delete it showing that n1 is copyied properly (no crosslinks)
NeuralNet n1 = *n0;											//calls the copy ctor , (allowed)

//------------------------
	
void drawaxes()
{
  glDisable(GL_LIGHTING);

  glBegin(GL_LINES);
  glColor3f(1.0,0.0,0.0);glVertex3f(-0.5,0.0,0.0);glVertex3f(1.2,0.0,0.0);
  glColor3f(0.0,1.0,0.0);glVertex3f(0.0,-0.5,0.0);glVertex3f(0.0,1.2,0.0);
  glColor3f(0.0,0.0,1.0);glVertex3f(0.0,0.0,-0.5);glVertex3f(0.0,0.0,1.2);
  glEnd();
}

void DrawBrainWave(NeuralNet &in,float floorsize,float floorstep,int x,int y)
{

	float i,j;

	float fftime = (float)GetTickCount() / 1500;
	float bp1 = 0.0;
	float bp2 = 1.0;
	float bp3 = 0.0;

	glPushMatrix();	
	glTranslatef((float)x,(float)y,0.0);


	for (i=0;i<floorsize;i+=floorstep)
	{
		glBegin(GL_TRIANGLE_STRIP);
		for (j=0;j<floorsize;j+=floorstep)
		{
			
				inputs[0] = j;
				inputs[1] = i;
				inputs[2] = fftime;	
				inputs[3] = bp1;
				inputs[4] = bp2;
				inputs[5] = bp3;
				in.RunNN();

				glColor3f(outputs[0],outputs[1],outputs[2]);
				glVertex3f(i,j,outputs[2]);	
				
				inputs[0] = j+floorstep;
				inputs[1] = i+floorstep;
				inputs[2] = fftime;	
				
				inputs[3] = bp1;
				inputs[4] = bp2;
				inputs[5] = bp3;
				in.RunNN();

				glColor3f(outputs[0],outputs[1],outputs[2]);
				glVertex3f(i+floorstep,j+floorstep,outputs[2]);
		}
		glEnd();
	}


	drawaxes();
	glPopMatrix();
}

void display(void)
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	imatrix16 rm(camquat);
	glTranslatef(0.0,0.0,-3.0);
	glMultMatrixf(rm.m);

	glTranslatef((float)-1,(float)-1,0.0);

	//instert drawing here 
	DrawBrainWave(n1,1.5f,0.08f,0,0);
	DrawBrainWave(n2,1.5f,0.08f,2,2);
	DrawBrainWave(n3,1.5f,0.08f,0,2);
	DrawBrainWave(n4,1.5f,0.08f,2,0);

	glutSwapBuffers();
	glutPostRedisplay();
}



void reshape(int w, int h)
{
  glViewport(0,0,(GLsizei) w,(GLsizei) h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(90,(float)w/h,0.01,100.0);
}



void motion(int x, int y)
{
	mx = x;my = y;

	xs = mx - oldmx;
	ys = oldmy - my;

	oldmx = mx;
	oldmy = my;

	iquaternion4 multquat;

	if (mode == MODE_ROTATE)
	{	
	
	}
	else if (mode == MODE_MOVE)
	{

	}
	else if (mode == MODE_CAM)
	{
		multquat.makeEulerRotation( ivector3( (float)ys/SCALE_MOTION,-(float)xs/SCALE_MOTION,0.0 ) );
		camquat.multiply(multquat);
	}
}


void mouse(int button,int state,int x,int y)
{
	oldmx = x;
	oldmy = y;

	if (button == GLUT_LEFT_BUTTON)
		mode = MODE_ROTATE;

	if (button == GLUT_RIGHT_BUTTON)
		mode = MODE_MOVE;

	
	if (button == GLUT_MIDDLE_BUTTON)
		mode = MODE_CAM;
}

void init(void)
{

	srand( (unsigned)time( NULL ) );	


	glEnable(GL_DEPTH_TEST);
	//glEnable(GL_LIGHTING);
	//glEnable(GL_LIGHT0);
	//glLightfv(GL_LIGHT0,GL_POSITION,above_origin);
	
	//glEnable(GL_BLEND);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	//glMaterialfv(GL_FRONT,GL_DIFFUSE,grey);
	//glMaterialfv(GL_FRONT,GL_SPECULAR,white);
	//glMaterialfv(GL_FRONT,GL_SHININESS,shiny);

	//glEnable(GL_POLYGON_OFFSET_FILL);
	//glPolygonOffset(0.1,0.1);
	//glDepthFunc(GL_LEQUAL);
	//glEnable(GL_DEPTH_TEST);
	//glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
	//glEnable(GL_BLEND);
	//glEnable(GL_POLYGON_SMOOTH);
	//glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	
	//here we delete n2

	meDelete n0;

	glPointSize(4.0f);

	//wont compile (not allowed)
	//n1 = n2;

	inputs[0] = 0;
	inputs[1] = 0;

	//n1.SetResponse(resp);
}

void key(GLubyte key, int x, int y)
{
  switch (key)
    {
			case 'a':resp -= 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);n2.SetResponse(resp);n3.SetResponse(resp);n4.SetResponse(resp);
			break;

			case 'd':resp += 0.05f;printf("response %f\n",resp);n1.SetResponse(resp);n2.SetResponse(resp);n3.SetResponse(resp);n4.SetResponse(resp);
			
			break;

			case 'r':im_random_seed_time();n1.Randomize();n2.Randomize();n3.Randomize();n4.Randomize();

			break;

			case '1':

			break;

			case '!':

			break;


			case 'l':

			break;


			case 'q':exit(0);
			break;
    }
	
}

int main(int argc, char **argv)
{
	meInitMemoryCheck();

	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
	glutInitWindowSize(WINDOW_SIZE,WINDOW_SIZE);
	glutInitWindowPosition(100,10);
	glutCreateWindow("Matts test");

	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(key);
	glutMotionFunc(motion);
	glutMouseFunc(mouse);

	init();
	glutMainLoop();

	return 0;
} 

/*---------------------------------------------
		
	MatNN library 

	By Matt Barnett
		(see References)

	-----------
	description
	-----------

		This file just keeps some ugly private NeuralNet
		classes out of the main header

	------------
	classes
	------------

		Neuron
		NeuronLayer	

	------------
	class design
	------------

	-----------
	References:	
	-----------

		* See matnn.h.
	
---------------------------------------------*/


//-------- headers etc -------
//---------------------------- 

#ifndef MATNN_NEURON_H_HEADER_132424623254
#define MATNN_NEURON_H_HEADER_132424623254

//----------------------------
//----------------------------

namespace MatLayeredNeuralNetLib
{
	//private class for neurons (clients cant mess with neurons)
	class Neuron
	{
		friend class NeuronLayer;

		private:

			Neuron(Neuron ©);				//copy
			Neuron &operator=(Neuron &rhs);		//assign

			Neuron(int in_inputs,float &in_sig_response);
			~Neuron();

			inline float Calc(float *inputs);
			void RandomizeWeights();
			inline void Mutate(float chance,float extent);

			float *weights;						//array of weights, owned by this neuron
			int num_inputs;
			float &sig_response;
	};	

	//private class for neuron layers (clients cant mess with layers)
	class NeuronLayer
	{
		friend class NeuralNet;
		friend class TrainableNeuralNet;

		private:

			NeuronLayer(NeuronLayer ©);				//copy
			NeuronLayer &operator=(NeuronLayer &rhs);	//assign

			NeuronLayer(int in_neurons,int in_num_inputs,float *finputs);
			~NeuronLayer();

			inline void RunLayer();
			void Randomize();
			inline int Mutate(float chance,float extent);

			inline float GetOutput(int neuron_index)			 { return output[neuron_index]; }
			inline void  SetResponse(float in_sig_response)		 { sig_response = in_sig_response; }

			int num_neurons;
			float *input;					//float array pointer (points to parent layers outputs);
			float *output;					//float array pointer (owned by this layer)
			Neuron **neurons;				//array of pointers to neurons
			float sig_response;
	};	
}

#endif 

/*---------------------------------------------
		
	MatNN library 

	By Matt Barnett
		(see References)

  
	-----------
	usage
	-----------

		I find that this is the section that i most often
		look for in other peoples docs, and neglect in my own.

		Look at the examples for usage comments, or examin
		the public interface.

	-----------
	description
	-----------

		This is not really a library yet, but its designed so it
		can be extended (hence the namespace) , and is as general as possible.
	
	  
		The Neural Net Specification
		(this is the same as the this library's
		pre-cursor, my other nn classes)

		* Neurons are made in layers

			input  layer  (any number of neurons)
			hidden layers (each has the same number of neurons)
			output layer  (any number of neurons)

		* Each layer (starting with the input layer) uses its neurons
		  to calculate the outputs for all the neurons in this layer.
		  This is stored in an array which the subsequent layer "knows".

		* The inital layer has no neuon inputs, this is where the input
		  from the "environment" enters the net. The input layer
		  has a pointer to inputs that are specified by the client

	-----------
	portability
	-----------

		* Depends on imath.h (see imath.h) (needs OpenGL, this is bad for a class a general as this)
		* Error messages need windows, hack the error class to support non windows
		* Neural nets are written to text files for easy editing and portability

	------------
	class goals
	------------

		NeuralNet : Serves as a base class for more specific neural nets (basic neural net function)
					Will not contain any code for training/adjustment of the net.

	------------
	class design
	------------

		The design is such that, the net class constructs itself
		appropriatly for the given parameters and is then ready to
		run.  This helps keep the rest of the interface small.

		Neural nets are written to a formatted text file for storage.
		This makes it easy to understand / edit and should make it
		more portable. (ie no struct size depenence)

		* No consideration given to const/non const neural nets/neurons/layers	( YET ! )
		* There is no default constructor.

		* Class design leaves some bounds checking to the client !!!!!!!
		  
		  (s/he simply has to stay within the bounds of his own arrays and be able
		   to remeber/work out the size of the array, (only once) )


	-----------
	References:	
	-----------

		* See imath.h, myerror.h
		
	
	
---------------------------------------------*/



//-------- headers etc -------
//---------------------------- 

#ifndef MATNN_MAIN_H_HEADER_012023254
#define MATNN_MAIN_H_HEADER_012023254

#include "../leakdetect.h"

#include <iostream>

//using std::cout;
//using std::endl;

#include <math.h>						//fpu maths
#include <cassert>						//assert
//#include <windows.h>					//opengl needs windows (comment for non WIN32 platforms)
//#include <gl/gl.h>					

#include "matnn_neuron.h"
#include "../imath/imath.h"
#include "../myerror.h"

//----------------------------
//----------------------------


//these limits are quite arbitrary
//just to prevent the user making an insane
//net just to test the memory handling (which isnt that hot)

#define MATNN_MAX_HIDDEN_LAYERS		256				//256 seems a reasonable limit for these
#define MATNN_MAX_IN_NEURONS		256				//
#define MATNN_MAX_OUT_NEURONS		256				//
#define MATNN_MAX_HIDDEN_NEURONS	256				//


namespace MatLayeredNeuralNetLib
{
	class NeuralNet
	{
		//neural net public interface
		
		public:
			
			//the client is responsible for storage of input & output arrays
			//but they only need to be set once, and then the you change the input
			//just be fiddling the array.

			NeuralNet(NeuralNet ©);				//copy
			NeuralNet &operator=(NeuralNet &rhs);	//assign

			NeuralNet(int in_neurons,int out_neurons,int hidden_neurons,int hidden_layers,float *in_finput,float *in_foutput);
			~NeuralNet();

			void RunNN();												//runs the neural net
			void SetResponse(float in_sig_response);					//alters the response constant for all neurons and layers in the net
			void SetInputOutput(float *in_finput,float *in_foutput);	//alters input/output arrays (null checking)
			void Randomize();											//dangerous to call on a trained net (it totally randomises weights to -1.0 + 1.0)



		protected:
		
			int num_input_neurons;					
			int num_output_neurons;					
			int num_hidden_neurons;					
			int num_hidden_layers;					

			NeuronLayer	*input_layer;				//will point to the first layer
			NeuronLayer	*output_layer;				//will point to the last layer
			NeuronLayer	**parray_hidden_layers;		//an array of pointers for the hidden layers
		
		private:
			float *input;
			float *output;

			float sig_response;
	};
	
	
	class TrainableNeuralNet : public NeuralNet
	{
		public:

			TrainableNeuralNet(int in_neurons,int out_neurons,int hidden_neurons,int hidden_layers,float *in_finput,float *in_foutput) : NeuralNet(in_neurons,out_neurons,hidden_neurons,hidden_layers,in_finput,in_foutput) {}
			~TrainableNeuralNet() {}
		
			int Mutate(float mutation_chance,float mutation_extent);		//returns the number of mutations that occured
			//ref or copy?
			//TrainableNeuralNet Clone()
			//TrainableNeuralNet Cross(TrainableNeuralNet &mate)
			//void Merge()
	};
	
}

#endif 

/*---------------------------------------------
		
	MatNN library : Implementation (cpp)

		* See Header file for more information

	By Matt Barnett
		(see References)
	
---------------------------------------------*/

#include "matnn.h"

//----------------------------
//----------------------------
using MatLayeredNeuralNetLib::Neuron;
using MatLayeredNeuralNetLib::NeuronLayer;
using MatLayeredNeuralNetLib::NeuralNet;
using MatLayeredNeuralNetLib::TrainableNeuralNet;

//the constructors look a bit nasty because a reference to
//the response variable/constant is passed down by the NeuralNet
//owner class (this allows one copy per NeuralNet , and when it changes
//there is no need to update each neuron/layer copy)

//
//	MatNN::Neuron
//

//do not copy neurons
//Neuron::Neuron(Neuron ©)
//{	 
//}
//Neuron &Neuron::operator=(Neuron &rhs)
//{
//
//	return *this;
//}


Neuron::Neuron(int in_inputs,float &in_sig_response) : sig_response(in_sig_response)
{
	num_inputs = in_inputs;

	//inputs = 0;
	//inputs = meNew Neuron * [num_inputs];
	
	weights = 0;
	weights = meNew float [num_inputs];

	RandomizeWeights();
}

void Neuron::RandomizeWeights()
{
	//randomize each weight 0.0f - 1.0f
	for (int i=0;i<num_inputs;i++)
	{
		weights[i] = RANDOM_CLAMP;
	}
}

void Neuron::Mutate(float chance,float extent)
{
	//mutate one of my weights by +/- extent
	//int i = RANDOM_NUM * num_inputs;
	//weights[i] += RANDOM_CLAMP * extent;
	
	//mutate any/some/all neuron weights
	for (int i=0;i<num_inputs;i++)
		if (RANDOM_NUM < chance) weights[i] += RANDOM_CLAMP * extent;
}

float Neuron::Calc(float *inputs)
{
	float out = 0;
	for (int x=0;x<num_inputs;x++)
	{
		out += weights[x] * inputs[x];
	
	}

	//sigmoid reponse curve
	//return ( 1 / (sig_response + (float)exp( - (out / 10) )) );
	
	//sinus response curve
	return im_sin( out / sig_response);
}

Neuron::~Neuron()
{
	//cleanup
	meDelete [] weights;weights = 0;
}

//
//	MatNN::NeuronLayer
//

NeuronLayer::NeuronLayer(NeuronLayer ©)
{
	num_neurons = copy.num_neurons;

	//cannot be assigned yet!!
	//this must be fixed by the caller
	input = 0;

	output = 0;
	output = meNew float [num_neurons];
	neurons = 0;
	neurons = meNew Neuron * [num_neurons];	
	
	//same number of neurons as previous layer
	for (int i=0;i<num_neurons;i++)
	{
		//create neurons
		neurons[i] = 0;	
		neurons[i] = meNew Neuron(copy.neurons[0]->num_inputs,sig_response);		//one float input for each neuron in the environmental/input layer
	
		//copy over the weights from copy
		memcpy( neurons[i]->weights,copy.neurons[i]->weights,sizeof(float) * copy.neurons[0]->num_inputs );
	}
}

//No =
//NeuronLayer &NeuronLayer::operator=(NeuronLayer &rhs)
//{
//
//	return *this;
//}

NeuronLayer::NeuronLayer(int in_neurons,int in_num_inputs,float *finputs)
{
	num_neurons = in_neurons;

	input = 0;
	input = finputs;

	output = 0;
	output = meNew float [num_neurons];
	neurons = 0;
	neurons = meNew Neuron * [num_neurons];	
	
	for (int i=0;i<num_neurons;i++)
	{
		neurons[i] = 0;	
		neurons[i] = meNew Neuron(in_num_inputs,sig_response);		//one float input for each neuron in the environmental/input layer
	}
	
	if (!input)
	{
		printf("NeuronLayer::NeuronLayer() no input array specified!\n");
		exit(0);
	}	
}

NeuronLayer::~NeuronLayer()
{	
	meDelete [] output;
	output = 0;

	for (int i=0;i<num_neurons;i++)
	{
		meDelete neurons[i];
		neurons[i] = 0;	
	}
	
	meDelete [] neurons;
	neurons = 0;
}	


void NeuronLayer::RunLayer()
{
	for (int x = 0;x<num_neurons;x++)
	{
		output[x] = neurons[x]->Calc(input);
	}
}

void NeuronLayer::Randomize()
{
	for (int x = 0;x<num_neurons;x++)
		neurons[x]->RandomizeWeights();
} 

int NeuronLayer::Mutate(float chance,float extent)
{
	int mutated = 0;

	//if chance, then mutate 
	for (int x = 0;x<num_neurons;x++)
	{
		if (RANDOM_NUM < chance)
		{
			neurons[x]->Mutate(chance,extent);
			mutated++;
		}
	}

	return mutated;
} 

//
//	MatNN::NeuralNet
//

NeuralNet::NeuralNet(int in_neurons,int out_neurons,int hidden_neurons,int hidden_layers,float *in_finput,float *in_foutput)
{
	Sleep(30);
	im_random_seed_time();

	//sanity checks first
	if (num_hidden_layers < 0 || num_hidden_layers > MATNN_MAX_HIDDEN_LAYERS)
	{
		cMyError e;e.FatalError("hidden_layers must be  0 => < = %i!",MATNN_MAX_HIDDEN_LAYERS);
		exit(EXIT_FAILURE);
	}

	if (in_neurons < 0 || in_neurons > MATNN_MAX_IN_NEURONS)
	{
		cMyError e;e.FatalError("in_neurons must be  0 => < = %i!",MATNN_MAX_IN_NEURONS);
		exit(EXIT_FAILURE);
	}

	if (out_neurons < 0 || out_neurons > MATNN_MAX_OUT_NEURONS)
	{
		cMyError e;e.FatalError("out_neurons must be  0 => < = %i!",MATNN_MAX_OUT_NEURONS);
		exit(EXIT_FAILURE);
	}

	if (hidden_neurons < 0 || hidden_neurons > MATNN_MAX_HIDDEN_NEURONS)
	{
		cMyError e;e.FatalError("hidden_neurons must be  0 => < = %i!",MATNN_MAX_HIDDEN_NEURONS);
		exit(EXIT_FAILURE);
	}


	//get input & output pointers
	input = in_finput;
	output = in_foutput;

	//grap those parameters
	num_input_neurons = in_neurons;
	num_output_neurons = out_neurons;
	num_hidden_neurons = hidden_neurons;
	num_hidden_layers = hidden_layers;
	
	//null em
	input_layer = 0;		
	output_layer = 0;
	parray_hidden_layers = 0;

	//create input layer
	input_layer = meNew NeuronLayer(num_input_neurons,num_input_neurons,input);		

	//an array of pointers to hidden layers	
	if (num_hidden_layers) 
	{
		parray_hidden_layers = meNew NeuronLayer * [num_hidden_layers];
		if (!parray_hidden_layers)
		{
			cMyError e;e.FatalError("Could not allocate %i hidden layers!",num_hidden_layers);
			
			meDelete input_layer;

			exit(EXIT_FAILURE);
		}
	}

	//null each one and create it
	for (int i=0;i<num_hidden_layers;i++)
	{
		parray_hidden_layers[i] = 0;

		if (i == 0)
			parray_hidden_layers[i] = meNew NeuronLayer(num_hidden_neurons, input_layer->num_neurons,input_layer->output);
		else
			parray_hidden_layers[i] = meNew NeuronLayer(num_hidden_neurons, parray_hidden_layers[i - 1]->num_neurons , parray_hidden_layers[i - 1]->output);
	}

	//output layer attached to last hidden layer
	//or to input layer (special case)
	if (num_hidden_layers == 0)
		output_layer = meNew NeuronLayer(num_output_neurons , input_layer->num_neurons , input_layer->output);
	else
		output_layer = meNew NeuronLayer(num_output_neurons , parray_hidden_layers[num_hidden_layers - 1]->num_neurons , parray_hidden_layers[num_hidden_layers - 1]->output);



	//anything reasonable is ok here
	//as long as the net has been constructed
	SetResponse(1.0f);

}

//No =
//NeuralNet &NeuralNet::operator=(NeuralNet &rhs)
//{
//	return *this;
//}

NeuralNet::NeuralNet(NeuralNet ©)
{
	//starting with the simple stuff
	
	sig_response = copy.sig_response;
	num_input_neurons = copy.num_input_neurons;
	num_output_neurons = copy.num_input_neurons;
	num_hidden_neurons = copy.num_hidden_neurons;
	num_hidden_layers = copy.num_hidden_layers;

	//null layers
	input_layer = 0;		
	output_layer = 0;
	parray_hidden_layers = 0;
	
	//get input & output pointers
	//its assumed that users want copied
	//nets to share input & output pointer
	//(this is the most obvious behaviour i think)

	input = copy.input;
	output = copy.output;

	//after copying layers
	//we must fix the input pointer to the new net

	//create input layer as copy
	input_layer = meNew NeuronLayer(*copy.input_layer);
	input_layer->input = input;

	//an array of pointers to hidden layers
	parray_hidden_layers = meNew NeuronLayer * [num_hidden_layers];
	//null each one and create it as copy
	for (int i=0;i<num_hidden_layers;i++) 
	{
		parray_hidden_layers[i] = 0;
		parray_hidden_layers[i] = meNew NeuronLayer(*copy.parray_hidden_layers[i]);

		if (i == 0)
			parray_hidden_layers[i]->input = input_layer->output;
		else
		{
			parray_hidden_layers[i]->input = parray_hidden_layers[i - 1]->output;
		}
	}
	//output layer attached to last hidden layer(?? or input layer ??) as copy
	output_layer = meNew NeuronLayer(*copy.output_layer);

	if (num_hidden_layers == 0)
		output_layer->input = input_layer->output;
	else
		output_layer->input = parray_hidden_layers[num_hidden_layers - 1]->output;


	SetResponse(copy.sig_response);
	
	//that should do it
}

NeuralNet::~NeuralNet()
{
	//cleanup
	meDelete input_layer;input_layer = 0;
	meDelete output_layer;output_layer = 0;

	//cleanup an array of pointers to allocated layers

	for (int i=0;i<num_hidden_layers;i++)
	{
		meDelete parray_hidden_layers[i];
		parray_hidden_layers[i] = 0;
	}

	meDelete [] parray_hidden_layers;
	parray_hidden_layers = 0;
}

void NeuralNet::RunNN()
{
	//run each layer
	input_layer->RunLayer();
	for (int i=0;i<num_hidden_layers;i++)
		parray_hidden_layers[i]->RunLayer();

	output_layer->RunLayer();

	//copy output from last layer
	for (i=0;i<num_output_neurons;i++)
		output[i] = output_layer->GetOutput(i);
}

void NeuralNet::SetResponse(float in_sig_response)
{
	sig_response = in_sig_response;
	input_layer->SetResponse(in_sig_response);
	for (int i=0;i<num_hidden_layers;i++)
		parray_hidden_layers[i]->SetResponse(in_sig_response);
	output_layer->SetResponse(in_sig_response);
}

void NeuralNet::SetInputOutput(float *in_finput,float *in_foutput)
{
	if (!in_finput || !in_foutput)
	{ 
		cMyError e;e.AlertError("SetInputOutput(), dont pass null pointers!!");return;
	}
	input = in_finput;output = in_foutput; 
}

void NeuralNet::Randomize()
{
	input_layer->Randomize();
	for (int i=0;i<num_hidden_layers;i++)
	{
		parray_hidden_layers[i]->Randomize();
	}
	output_layer->Randomize();
}

//
//	MatNN::TrainableNeuralNet
//

int TrainableNeuralNet::Mutate(float chance,float extent)
{
	int mutated = 0;

	mutated += input_layer->Mutate(chance,extent);
	for (int i=0;i<num_hidden_layers;i++)
	{
		mutated += parray_hidden_layers[i]->Mutate(chance,extent);
	}
	mutated += output_layer->Mutate(chance,extent);

	return mutated;
} 

/*-----------------------------------------------------------------------
		
	Inline Math library : Matrix and Quaternion implementation
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/


#ifndef _IMATH_MQ_IMP_HEADER_4778097634
#define _IMATH_MQ_IMP_HEADER_4778097634

//----------------------------
//----------------------------

class iquaternion4
{
	public:
		
		//quaternion data (w,x,y,z)
		union
		{
			float q[4];	
			struct { float w, x, y, z; };
		};
		
		//-------------- operators ---
		//---------------------------- 
		inline float operator[](int subscript) 
		{
			assert( subscript < 4 && subscript >= 0 );
			return q[subscript];
		}

		//assignment
		inline const iquaternion4 &operator=(const iquaternion4 &right)
		{
			if (&right != this) {				//check for self assign
				w = right.w;
				x = right.x;					//copy right to left					
				y = right.y;	
				z = right.z;	
			}
			return *this;						//enables x = y = z						
		}
	
		//----------------------------
		//----------------------------
		
		//------- construction etc ---
		//----------------------------

		//default ctor
		inline iquaternion4()		
		{
			makeMultiplicationIdentity(); 
		}
		
		//ctor											
		inline iquaternion4(float setw,float setx,float sety,float setz)
		{
			set(setw,setx,sety,setz);
		}												
	
		//ctor											
		inline iquaternion4(float angle,ivector3 &axis)
		{
			makeAngleAxisRotation(angle,axis);
		}
		
		//ctor											
		explicit inline iquaternion4(ivector3 &euler)
		{
			makeEulerRotation(euler);
		}	

		//dtor(unnessesary)
		inline ~iquaternion4()	
		{
		}						
	
		//----------------------------
		//---------------------------- 

		//--------- void functions ---
		//----------------------------
		inline void set(float setw,float setx,float sety,float setz)
		{
			w = setw;
			x =	setx;
			y = sety;
			z = setz;
		}
		
		//dumps the quaternion componets
		inline const void cdump(const char *caption)  
		{
			cout << caption << w << " " << x << " " << y << " " << z << endl;
		}	

		//theroy (identity quaternion)
		//
		//There are two identity quaternions.
		//In this context, we only use the one that multiplies with another quaternion 
		//to give the same quaternion ( the multiplication identity )

		inline void makeMultiplicationIdentity()
		{
			q[0] = 1.0;		//w
			q[1] = 0.0;		//x
			q[2] = 0.0;		//y
			q[3] = 0.0;		//z
		}

		//from Ref 5. (* See imath.h)
		//(http://gamedev.net/reference/articles/article1691.asp#Q53)
		//Qr = Q1.Q2 = ( w1.w2 - v1.v2, w1.v2 + w2.v1 + v1 x v2 )
		//
		//where v1 = (x,y,z) of Q1
		//	w1 = (w)     of Q1
		//	v2 = (x,y,z) of Q2
		//	w2 = (w)     of Q2
		//
	    //and both . and x are the standard vector dot and cross products.
		//its easy(ish) to follow that the multiply function 
		//doesnt chane the quaternion if multiply is the  identity
		inline void multiply(const iquaternion4 &multiply)
		{
			iquaternion4 temp;
			temp.x = ((w * multiply.x) + (x * multiply.w) + (y * multiply.z) - (z * multiply.y));
			
			//for multiply == identity quaternion
			//x = 1 * x + 0 * w + 0 * y + 0 * z;
			//x = x

			temp.y = ((w * multiply.y) - (x * multiply.z) + (y * multiply.w) + (z * multiply.x));
			temp.z = ((w * multiply.z) + (x * multiply.y) - (y * multiply.x) + (z * multiply.w));
			temp.w = ((w * multiply.w) - (x * multiply.x) - (y * multiply.y) - (z * multiply.z));
			*this = temp;
		}

		
		inline void makeAngleAxisRotation(float angle,ivector3 &axis)
		{
			const float hangle(angle * 0.5f);
			const float sina(im_sin(hangle) );
			x = (axis.x * sina);
			y = (axis.y * sina);
			z = (axis.z * sina);
			w = cos(hangle);
		}

		//from Ref 5. (* See imath.h)
		//(more or less unchanged from what that articles has)
		//(http://gamedev.net/reference/articles/article1691.asp#Q60)
			//" The following more or less comes from:
			//  http://vered.rose.utoronto.ca/people/david_dir/GEMS/GEMS.html "

		inline void makeEulerRotation(ivector3 &angles)
		{
			//Pitch->X axis, Yaw->Y axis, Roll->Z axis
		
			const float sinPitch(im_sin(angles.y*0.5f));
			const float cosPitch(im_cos(angles.y*0.5f));

			const float sinYaw(im_sin(angles.z*0.5f));
			const float cosYaw(im_cos(angles.z*0.5f));

			const float sinRoll(im_sin(angles.x*0.5f));
			const float cosRoll(im_cos(angles.x*0.5f));
			const float cosPitch_cosYaw(cosPitch*cosYaw);
			const float sinPitch_sinYaw(sinPitch*sinYaw);

			x = sinRoll * cosPitch_cosYaw   - cosRoll * sinPitch_sinYaw;
			y = cosRoll * sinPitch * cosYaw + sinRoll * cosPitch * sinYaw;
			z = cosRoll * cosPitch * sinYaw - sinRoll * sinPitch * cosYaw;
			w = cosRoll * cosPitch_cosYaw   + sinRoll * sinPitch_sinYaw;
		}

		//this is the alternative
		//multiplying euler angles in the correct order
		//i _think_ that the above is simply a re-arrangment
		//but i think the ordering here is wrong or somthing
		//(whatever, the above works the way i want and this doesnt)

		//
		//inline void makeEulerRotation(ivector3 &angles)
		//{
		//	iquaternion4 qx( im_cos(angles.x/2),im_sin(angles.x/2),0,0 );
		//	iquaternion4 qy( im_cos(angles.y/2),0,im_sin(angles.y/2),0 );	
		//	iquaternion4 qz( im_cos(angles.z/2),0,0,im_sin(angles.z/2) );
		//	qx.multiply ( qy );
		//	qx.multiply ( qz );
		//	*this = qx;
		//}
		//
	
		//---------------------------- 
		//---------------------------- 

		//-------------- functions ---
		//----------------------------

		//well this seems to draw the expected vector
		//but its easier to extract the basis from the
		//matrix, since the quat is converted to one in the
		//end
		const ivector3 &direction()
		{
			ivector3 out;
			out.x = 2.0f * ( x * z - w * y );
			out.y = 2.0f * ( y * z + w * x );
			out.z = 1.0f - 2.0f * ( x * x - y * y );

			out.normalise();
			return out;
		}

		inline float dot(const iquaternion4 &right)
		{
			return ( (x * right.x) + (y * right.y) +
					 (z * right.z) + (w * right.w) );
		}
		inline float lengthsqr()
		{
			return dot(*this);
		}
		inline float length()
		{
			return sqrt( lengthsqr() );
		}

		//just like a vector ( UNTESTED )
		inline void normalise()
		{
			const float leng(length() );
			x /= leng;
			y /= leng;
			z /= leng;
			w /= leng;
		}

		//?? inline iquaternion4 &lerp(iquaternion4 &other,float t);
		//?? inline iquaternion4 &slerp(iquaternion4 &other,float t);	
		//?? inline void invert(iquaternion4 &multiply);

		//---------------------------- 
		//---------------------------- 
};

//-------------- imatrix16 -----------------------------------------------------------
//
//  ---------------------------------
//  Row and clomn major format issues
//  ---------------------------------
//
//	Currently, this is an OpenGL freindly matrix
//  class.  Direct3D stores matricies in a different
//  order to OpenGL, and this would need to be accounted
//  for if Direct3D was used.  Also, an engine might need
//  to use both formats if it supports both APIs.  This is
//  why i have differed implementing this feature untill i
//  know how to get it right.
//
//

//matrix in colomn major format (consistant with opengl)
// 
//		eg { m[0] m[4] m[8]  m[12] 
//		 	 m[1] m[5] m[9]  m[13]
//			 m[2] m[6] m[10] m[14]
//			 m[3] m[7] m[11] m[15] }
//

//note: it would be cool to have a union for the matrix data allowing sensible names
//      for different parts of the matrix.  For instance, the first 3x3 block represents
//      the rotation in the form of the orthonormal basis axes. The matrix also stores scaling etc.
//
//	    somthing like...
//
//							
//		union
//		{
//			float m[16];
//			struct { float col0[4],col1[4],col2[4],col3[4]; };
//			
//		};		
//
//---------------------------------------------------------------------------------

class imatrix16
{
	public:

		//------- construction etc ---
		//----------------------------
		float m[16];																						//matrix data
		inline imatrix16()										 { makeIdentity(); }						//ctor											
		inline imatrix16(iplane4 &cast_floor,ivector3 lightpos)  { makeProjection(cast_floor,lightpos); }	//ctor
		explicit inline imatrix16(iquaternion4 &createfrom)		 { makefromQuaternion(createfrom); }		//ctor

		inline ~imatrix16()	{}																			//dtor(unnessesary)
		//---------------------------- 
		//---------------------------- 

		//--------- void functions ---
		//----------------------------	
		
		//dumps the matrix componets
		inline const void cdump(const char *caption)  
		{
			cout << caption << " row major\n";

			cout << caption << " row 1 " << m[0] << " " << m[1] << " " << m[2] << " " << m[3] << endl;
			cout << caption << " row 2 " << m[4] << " " << m[5] << " " << m[6] << " " << m[7] << endl;
			cout << caption << " row 3 " << m[8] << " " << m[9] << " " << m[10] << " " << m[11] << endl;
			cout << caption << " row 4 " << m[2] << " " << m[13] << " " << m[14] << " " << m[15] << endl;
			
			cout << caption << " col major\n";

			cout << caption << " col 1 " << m[0] << " " << m[4] << " " << m[8] << " " << m[12] << endl;
			cout << caption << " col 2 " << m[1] << " " << m[5] << " " << m[9] << " " << m[13] << endl;
			cout << caption << " col 3 " << m[2] << " " << m[6] << " " << m[10] << " " << m[14] << endl;
			cout << caption << " col 4 " << m[3] << " " << m[7] << " " << m[11] << " " << m[15] << endl << endl;
		}	

		//inline void invert
		//inline void transpose
		inline void makeIdentity() 
		{		
			m[0] = 1.0;
			m[1] = 0.0;
			m[2] = 0.0;
			m[3] = 0.0;

			m[4] = 0.0;
			m[5] = 1.0; 
			m[6] = 0.0;
			m[7] = 0.0;

			m[8] = 0.0; 
			m[9] = 0.0; 
			m[10] = 1.0;
			m[11] = 0.0; 

			//last right hand colomn
			m[12] = 0.0; 
			m[13] = 0.0; 
			m[14] = 0.0; 
			m[15] = 1.0; 
		}

		//this makes a projection matrix
		//which can be fed into a graphics API
		//
		//is by no means perfect as a shadow/projection
		//matrix (still a work in progress)

		//from http://www.sgi.com/software/opengl/advanced96/node46.html
		inline void makeProjection(iplane4 &cast_floor,ivector3 &light_pos) 
		{
			float light_posw = 0.0;							//0.0?

			float dot = cast_floor.norm.x * light_pos.x +
					    cast_floor.norm.y * light_pos.y +
						cast_floor.norm.z * light_pos.z +
						cast_floor.d * light_posw;				

			m[0] = dot - light_pos.x * cast_floor.norm.x;
			m[1] = 0.0 - light_pos.y * cast_floor.norm.x;
			m[2] = 0.0 - light_pos.z * cast_floor.norm.x;
			m[3] = 0.0 - light_posw  * cast_floor.norm.x;
    
			m[4] = 0.0 - light_pos.x * cast_floor.norm.y;
			m[5] = dot - light_pos.y * cast_floor.norm.y;
			m[6] = 0.0 - light_pos.z * cast_floor.norm.y;
			m[7] = 0.0 - light_posw  * cast_floor.norm.y;
    
			m[8 ] = 0.0 - light_pos.x * cast_floor.norm.z;
			m[9 ] = 0.0 - light_pos.y * cast_floor.norm.z;
			m[10] = dot - light_pos.z * cast_floor.norm.z;
			m[11] = 0.0 - light_posw  * cast_floor.norm.z;
    
			m[12] = 0.0 - light_pos.x * cast_floor.d;
			m[13] = 0.0 - light_pos.y * cast_floor.d;
			m[14] = 0.0 - light_pos.z * cast_floor.d;
			m[15] = dot - light_posw  * cast_floor.d;
		} 
	
		
		//from gametutorials.net (see ref 6)
		inline void makefromQuaternion(iquaternion4 &q)		 
		{
			
			// First row
			m[ 0] = 1.0f - 2.0f * ( q.y * q.y + q.z * q.z );  
			m[ 1] = 2.0f * ( q.x * q.y - q.w * q.z );  
			m[ 2] = 2.0f * ( q.x * q.z + q.w * q.y );  
			m[ 3] = 0.0f;  

			// Second row
			m[ 4] = 2.0f * ( q.x * q.y + q.w * q.z );  
			m[ 5] = 1.0f - 2.0f * ( q.x * q.x + q.z * q.z );  
			m[ 6] = 2.0f * ( q.y * q.z - q.w * q.x );  
			m[ 7] = 0.0f;  

			// Third row
			m[ 8] = 2.0f * ( q.x * q.z - q.w * q.y );  
			m[ 9] = 2.0f * ( q.y * q.z + q.w * q.x );  
			m[10] = 1.0f - 2.0f * ( q.x * q.x + q.y * q.y );  
			m[11] = 0.0f;  

			// Fourth row
			m[12] = 0;  
			m[13] = 0;  
			m[14] = 0;  
			m[15] = 1.0f;
		}

		//this functions could be avoided by
		//providing better access to the matrix's
		//internal values
		//based on theroy (Q5 of matrix FAQ) (see ref 5)
		inline ibasis9 &getBasis()
		{
			//prepare the basis based on the matrix
			ivector3 xa( m[0],m[1],m[2] );
			ivector3 ya( m[4],m[5],m[6] );
			ivector3 za( m[8],m[9],m[10] );
			//make it and return it
			ibasis9 out(xa,ya,za);
	
			return out;
		}

		//
		void multiply_vectors(ivector3 *vector_array,int num_vectors)
		{
			assert(vector_array);
			assert(num_vectors == 1);

			ivector3 out;
			out.x = m[0] * vector_array->x + m[4] * vector_array->y + m[8] * vector_array->z;
			out.y = m[1] * vector_array->x + m[5] * vector_array->y + m[9] * vector_array->z;
			out.z = m[2] * vector_array->x + m[6] * vector_array->y + m[10] * vector_array->z;
			*vector_array = out;
		}	

		//??? inline void makePerspectiveProjection() {} 
		//??? inline void makeOrthographicProjection() {}
		//??? inline void makeEulerRotation(float angle, ivector3 &axis) {}
		//??? inline void makeAngleAxisRotation(float angle, ivector3 &axis) {}

		//---------------------------- 
		//---------------------------- 
};

#endif
 

/*-----------------------------------------------------------------------
		
	Inline Math library : Lookup table implementation
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

#ifndef _IMATH_LUT_IMP_HEADER_223457642
#define _IMATH_LUT_IMP_HEADER_223457642

//------- float_lookuptable ------------------------------------
//
//	Before i begin.. Yes, i know that the c++ 
//  standard library has a lookup table template class.
//
//
//  This is a simple generic lookup table(lut) class
//  for functions that return a single float
//	and take a single float parameter
//	(GPG2 see ref 2) has a better lookuptable)
//	
//	this works well for sinf etc
//	
//.	With bounds checking its quite a bit slower.
//	Without it, i think a lookup maybe amount to
//	only a few cycles (WRONG!, in retrospect, see asm anal),
//  because looking up consists of
//	an inline (no fuction call) function
//	which performs a single divide, and 
//  uses the clever real2int to get an int
//	to index the array with
//
//  ------ design notes -----
//
//  i decided to forgo inheritance
//
//  (it would be nice to have a base class
//   lookuptable and inhert specific
//	 table types)
//  
//  templates would probably be a good idea
//
//  (but an int indexing table would be
//   quite different than a float "indexing"
//   one, and could be faster, so the template
//   would have to account for that)
//
//	In the end decided to make this a simple class
//  with inlined functions specific to functions of
//	the form
//
//  float somefunction(float someparam);
//
//	becomes...
// 
//  float somelut.lookup(someparam);
//
//  The other lut class, is specifically
//  for generating random numbers and is
//  pretty much hard-coded to that effect.
//
//  ---------- performance notes ---------------
//  --------------------------------------------
//
//	*test yourself using test_speed(int);
//	*test in both Debug and Release modes
//    (only take Release mode seriously, but
//		its usefull to see the difference between
//			debug & release)
//
//
//	Debug mode is still faster than parent functions
//	which are them selves slower than in Release mode.
//	
//  faster than parent tri functions
//  generally slower for sqrtf
//  faster than on-the fly random number generation
//		(not yet suited to rng, which really
//		 deserves its own lut class)
//
//  --------------------------------------------

//  ---------- mem usage notes ---------------
//  ------------------------------------------
//
//100 floats	400bytes			//maybe for some stuff
//1000 floats = 4k					//vagly accurate
//10000 floats = 40k etc			//pretty accurate
//	
//	Need more accuracy? 
//  (fix class before using more samples)
//  
//	1)
//
//	I've made that the upper limit for the number of 
//	samples. Not due to memory usage but processor
//  cache usage. (im guessing here, from accumlated knowledge)
//	We want to be sure that the whole table fits in the cache
//  so that its really fast. (this afik occurs automagically)
//	If the table is to big, the chances are that the processor
//  will thrash the table. Much like a gfx card will thrash textures
//  between system/texture memory if the application tries to use to
//  much texture memory.					

//  ------------------------------------------

//if you want a quick report when a lookup table constructs then define FLUT_VERBOSE
//if you want to remove the construction message, rem FLUT_VERBOSE out

#define FLUT_VERBOSE

class float_lookuptable
{
	public:

	//constructor provides everything
	float_lookuptable(int innumsamples,float instart,float inend,float (*function)(float) )
	{
		//store the pointer to the function for testing later

		parentfunction = function;
		num_samples = innumsamples;
		
		if (innumsamples > 10000)
		{
			printf("float_lookuptable::float_lookuptable() num_samples > 10000!! rtfm\n");
			exit(0);
		}


		start = instart;
		end = inend;
		samples = 0;
		samples = meNew float [num_samples];
		assert(samples);

		//may as well fabs() to be sure we get a positive
		//delta (as this is precalc
		float diff = fabs(end - start);
		delta = diff / (float)(num_samples);
	  
		float x = start;

		#ifdef FLUT_VERBOSE
			printf("float_lookuptable( start(%f),end(%f),delta(%f),num(%i) )\n",start,end,delta,num_samples);
		#endif

		for (int i = 0; i < num_samples; i++)
		{
			//assert(index < num_samples);
		
			x = start + (diff * (float)i / (num_samples) );
			samples[i] = (*parentfunction)(x);
		}
	}	

	//looks up like the function
	inline float lookup(float lookup_what) const
	{
		lookup_what /= delta; 
		int il = im_real2int(lookup_what);		//fast
		//int il = (int)lookup_what ;			//doesnt seem any slower
		
		//we only need to check the index is safe
		//to access the array.  These are faster
		//than comparing floats
		// (we could use a smaller precision index
		//  and possibly gain more speed)



		//was previously this line (which caused the assertion(il < num_samples) to fail (+3v!L />U|\|7i/\/\3 3/>/>0z)
		//while (il >= num_samples) il -= num_samples;
		while (il >= num_samples) il -= num_samples;
		while (il < 0) il += num_samples;
		
		assert(il >= 0);
		assert(il < num_samples);

		return samples[il];
	}

	
	void test_speed(int test_samples) const
	{	
		printf("\n\ntesting speed\n\n");

		//i think as long as we do enough samples
		//it will be an ok test

		long before,after,deltatime;
		
		float some_number = two_pi/1.23456789;		//get an angle
		float out;									//temp to store output of function
		int x = 0;									//counter

		before = GetTickCount();					//time before
	
		for (x=0;x<test_samples;x++)				//do the test
			out = lookup(some_number);
		
		after = GetTickCount();						//time after
		deltatime = after - before;					//delta

		printf("lookup(%f) x %i = (%f) in %i msec\n",some_number,samples,out,deltatime);

		before = GetTickCount();					//time before

		for (x=0;x<test_samples;x++)				//do the test
			out = parentfunction(some_number);

		after = GetTickCount();						//time after
		deltatime = after - before;					//delta


		printf("parentfunction(%f) x %i = (%f) in %i msec\n",some_number,samples,out,deltatime);
	}

	~float_lookuptable()
	{
		meDelete [] samples;
	}


	private:
			float *samples;
			int num_samples;
			float start,end;
			float delta;
			float (*parentfunction)(float);


	//no assignment or copy kids
	operator=(float_lookuptable &);
	float_lookuptable(float_lookuptable &);
};

//------------------------------------- 
//hard coded random number lookup table
//-------------------------------------

//note: this is slower in debug mode than
//      rand();
//
//      but in release mode it beats rand()
//      however, because rand()  is so fast
//      its probably not worth the memory

//quick prototype
void im_random_seed_time();

class float_random_lookuptable
{
	public:
		float_random_lookuptable(int numsamples)
		{
			//seed the generator
			im_random_seed_time();

			//class general
			num_samples = numsamples;
			current_index = 0;

			//alloc memory
			samples = 0;
			samples = meNew float [num_samples];
			assert(samples);

			//fill the array
			for (int x=0;x<num_samples;x++)
			{
				samples[x] = ((float)rand()/(RAND_MAX+1));
			}

		}
		~float_random_lookuptable()
		{
			meDelete [] samples;
		}

		inline float lookup_random_num()
		{
			return samples[current_index];
		
			current_index++;

			//loop around again
			if (current_index > num_samples) current_index = 0;
		
		}
			
		//not const becauase it calls lookup_random_num
		void test_speed(int test_samples)
		{	
			printf("\n\ntesting speed\n\n");

			//i think as long as we do enough samples
			//it will be an ok test

			long before,after,deltatime;

			float out;									//temp to store output of function
			int x = 0;									//counter

			before = GetTickCount();					//time before
		
			for (x=0;x<test_samples;x++)				//do the test
				out = lookup_random_num();
			
			after = GetTickCount();						//time after
			deltatime = after - before;					//delta

			printf("lookup_random_num() x %i = %i msec\n",test_samples,deltatime);

			before = GetTickCount();					//time before

			for (x=0;x<test_samples;x++)				//do the test
				out = ((float)rand()/(RAND_MAX+1));

			after = GetTickCount();						//time after
			deltatime = after - before;					//delta

			printf("((float)rand()/(RAND_MAX+1)) x %i = %i msec\n",test_samples,deltatime);
		}

	private:
		int num_samples,current_index;
		float *samples;
};

#endif 

/*---------------------------------------------
		
	Inline Math library 
	(imath , im_)

	By Matt Barnett
		(see References)

	-----------
	portability
	-----------

		exclusivly uses the internal float datatype throughout
		and precision and portability are dependent on that.

	------------
	classes
	------------
		
		ivector3
		ibasis9			
		imatrix16		
		iquaternion4
		iray6
		iplane4
		ibox3
		

  	----------------
	other components
	----------------

		* float constants collection

		* lookup table class (for functions such as sin/cos/tan/acos etc)
		* random number lookup table class

		* interfaces to generalise common math functions
		
			im_random_num()			//uses lookup tables
			im_random_clamp()
			im_sin()				//uses lookup tables
			im_cos()				//uses lookup tables
			im_acos()
			im_tan()				//uses lookup tables
			im_sqrt()			


  	-------------------------------
	external components [optional]
	-------------------------------
		
		ibox3
		(closly intergrated with imath)

	------------
	class design
	------------

		imath.h					//include this file
		imath.cpp				//add this file to project
		imath.lib?				//or link to this?

		this file, imath.h sould include the files

		imath_vmq.h				//vector mathrix and quats all together
		imath_trace.h			//rays and planes together
		imath_box.h				//bounding boxes
		imath_lut.h				//lookup tables

		* All functions are inlined (the i stands for inline) ( this should be optimal, ie: the compilier chooses )

		* Operators are generally not implemented using member functions calls.  This is deliberate.  Member function
		  calls quite specifically operate on an object, and operators act between them.
		
		* Integral types are passed/returned by value
		  math objects are passed/returned by reference where appropriate

		* All members are public (for core classes). Classes are
		  designed to mimic integral types(hmm, dont that too seriously).

	-----------
	References:		( INCOMPLETE )
	-----------
		
		 1. PI constants from windows calculator :)
			Allthough you cant store such an accurate PI very easily
			its worth noting that wincalc will give pi to 21 decimal places.

		 2. Game Programming Gems 2
		 3. Elemetry Linear Algebra	  
		 4. Real2Int function By Sree (* See imath_real2int.h )
		 5. www.gametutorials.com	  (specifically the quaternion example, great site)
		 6. matrix & quaternion FAQ	  (everything about m's and q's)
		 7. Effective C++ 2e		  (great book)
		 8. flipcode http://www.flipcode.com/cgi-bin/msg.cgi?showThread=00002476&forum=general&id=-1
		    (angle between 2 vectors thread) Paul Hope's angle code
	
---------------------------------------------*/



//-------- headers etc -------
//---------------------------- 

#ifndef _IMATH_MAIN_H_HEADER_069389570
#define _IMATH_MAIN_H_HEADER_069389570

#include "../leakdetect.h"

#include <iostream>

using std::cout;
using std::endl;

#include <math.h>						//fpu maths
#include <cassert>						//assert
#include <windows.h>					//opengl needs windows (comment for non WIN32 platforms)
#include <gl/gl.h>						//opengl for drawing vectors etc (minimal, not important)

#include "../myerror.h"

//----------------------------
//----------------------------



//----------------------------
//----------------------------

//alias/type/macro your favorites
#define		RANDOM_NUM							im_random_num()
#define		RANDOM_CLAMP						im_random_clamp()

//macroed settings
#define		_IVECTOR3_VERBOSE_NORMALISE_FAIL_

//----------------------------
//----------------------------



//----------------------------
//----------------------------

class ivector3;
class float_lookuptable;

//----------------------------
//----------------------------



//---- constants -------------------
//
// see imath.cpp for
// constant definitions and comments
//----------------------------------

extern const float  cnf_epsilon;			//"small" value (considered 0)
extern const float			 pi;			
extern const float		half_pi;			
extern const float		 two_pi;
extern const float _pi_over_180;
extern const float _180_over_pi;

//extern const ivector3 xaxis;
//extern const ivector3 yaxis;
//extern const ivector3 zaxis;
//extern const ivector3 zero;

//extern const ivector3 red;  
//extern const ivector3 green;
//extern const ivector3 blue;
//extern const ivector3 black;

//extern const ivector3 cyan;
//extern const ivector3 magenta;
//extern const ivector3 yellow;

extern const float_lookuptable sintable;
extern const float_lookuptable costable;
extern const float_lookuptable tantable;

//----------------------------
//----------------------------



//-- imath component headers --
//-----------------------------

#include "imath_real2int.h"		//need first by lut
#include "imath_lut.h"			//lookup tables
#include "imath_core.h"			//trig, sqrt,  etc

#include "imath_vector.h"		//vector is needed by trace
#include "imath_trace.h"		//trace is needed by matrix
#include "imath_mq.h"			//matrix and quaternion

#include "imath_box.h"			//bounding boxess

//---------------------------- 
//---------------------------- 

#endif 

/*---------------------------------------------
		
	Inline Math library : Constants implementation (cpp)
	(imath , im_)

	* See Header file for more information

	By Matt Barnett
		(see References)

---------------------------------------------*/

#include "imath.h"

//---- float constants --------
//----------------------------- 

//data ranges (according to msvc help)
//a 4 byte float can hold 7 digits
//an 8 byte double can hold 15 digits
//a long double (10 bytes) can hold 19 digits

//hence all constants match the precision
//of the the type that stores them

const float cnf_epsilon(0.0001f);			//"small" value (considered 0)

//pi constants
//precalced often used pi values

//theroy (radians)
//
//  half_pi radians = 90 degrees
//  pi radians      = 180 degrees
//  two_pi radians  = 360 degeres

const float			 pi(3.141592f);			//(* See ref 1)
const float		half_pi(1.570796f);			//(* See ref 1)
const float		 two_pi(6.283185f);			//(* See ref 1)
const float _pi_over_180(0.0174532f);
const float _180_over_pi(57.295779f);

//---------------------------- 
//----------------------------



//------ ivector3 constants -----
//-------------------------------

extern const ivector3 xaxis(1.0,0.0,0.0);
extern const ivector3 yaxis(0.0,1.0,0.0);
extern const ivector3 zaxis(0.0,0.0,1.0);
extern const ivector3 zero (0.0,0.0,0.0);

//a few colors
extern const ivector3 red  (1.0,0.0,0.0);  
extern const ivector3 green(0.0,1.0,0.0);
extern const ivector3 blue (0.0,0.0,1.0);
extern const ivector3 black(0.0,0.0,0.0);

extern const ivector3 cyan   (0.0,1.0,1.0);
extern const ivector3 magenta(1.0,0.0,1.0);
extern const ivector3 yellow (1.0,1.0,0.0);

//----------------------------
//----------------------------



//-- lookup table constants --
//----------------------------

const float_lookuptable	sintable(5000,0,two_pi,&sinf);
const float_lookuptable	costable(5000,0,two_pi,&cosf);
const float_lookuptable	tantable(5000,0,two_pi,&tanf);

//----------------------------
//----------------------------
 

/*-----------------------------------------------------------------------
		
	Inline Math library : Non member function implementation
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

//-------- headers etc -------
//---------------------------- 

#ifndef _IMATH_CORE_IMP_HEADER_354766003
#define _IMATH_CORE_IMP_HEADER_354766003

//--- non member functions ---
//---------------------------- 
// 
// some may seem silly stuff like
//
// inline float im_sin(float angle) { return sinf(angle) }
// inline float im_sqrt(float val)  { return sqrtf(val)  }
//
// this are interfaces so that these functions can be replaced
// with more efficient versions. (trig functions are replaced with
// lookup tables)
//

//theroy (degrees & radians)
//
//to avoid confusion it is best _NOT_ to convert angle measures.
//Sick with radians untill conversion to degrees is ABSOLUTELY nessesary
//(try not to store degrees, throw them away if appropriate)

//converts radians to degrees
inline float im_r2d(float angle)
{
	return 	(angle * _180_over_pi);
}

//converts degrees to radians
inline float im_d2r(float angle)
{
	return	(angle * _pi_over_180);
}

//random number generation
//seed with the generator (custom seed
inline void im_random_seed(unsigned int seed)
{
	srand( seed );
}
//seed with time
inline void im_random_seed_time()
{
	im_random_seed( GetTickCount() );
}

//random number retrival functions
//returns a float between 0 & 1

inline float im_random_num()
{
	return ((float)rand()/(RAND_MAX+1));
}

//all based on im_random_num
inline float im_random_num(float range)   { return im_random_num() * range;}		//returns a float between 0 & range
inline float im_random_clamp()			  { return -1 + im_random_num() * 2;  }		//returns a float in the range -1.0f -> + 1.0f
inline float im_random_clamp(float range) { return im_random_clamp() * range; }		//returns a float between -range & +range

//sqrt
inline float im_sqrt(float in)
{
	return sqrtf(in);
}

inline float im_sin(float in)
{
	//return sinf(in);
	return sintable.lookup(in);
}

inline float im_cos(float in)
{
	//return cosf(in);
	return costable.lookup(in);
}

inline float im_acos(float in)
{
	return acosf(in);
}

inline float im_tan(float in)
{
	//return tanf(in);
	return tantable.lookup(in);
}

//---------------------------- 
//----------------------------

#endif 

/*-----------------------------------------------------------------------
		
	Inline Math library : Real to Integer conversion by Sree
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

#ifndef _IMATH_REAL2INT_IMP_HEADER_156858806791
#define _IMATH_REAL2INT_IMP_HEADER_156858806791

// Float/Double -> Int
//
// By Sree
//
//  (from web article)
//	Sree's Real2Int
//
//	There are quite a few ways to fix this. One of the best general conversion utilities I have is one
//	that Sree Kotay wrote. It's only safe to use when your FPU is in double-precision mode, so you have
//	to be careful about it if you're switching precision like above. It'll basically return '0' all the
//	time if you're in single precision.
//
//	--- notes by matt ---
//  Somthing is wierd here
//  My lut class attempts to switch to single precision mode
//  but this function must still be working, or (i think) the
//  core im_ functions would all fall down.
//
//	(I've recently discovered that this discussion bears a lot of resemblance to Chris Hecker's article
//	on the topic, so be sure to read that if you can't figure out what I'm saying.)
//	This function is basically an ANSI C compliant conversion -- it always chops, always does the right thing,
//	for positive or negative inputs.
//
//

//-- performace notes by matt --
//
// the blurb states that (int)somefloat casts
// are very slow.  My lookup table class
// didnt seem to show that(the cast was no slower that this function),
// but i have left this function in, and the class uses it.
// (i guess it may be noticable for doubles (which i am avoiding)
// flipcode threads have also mentioned (somewhere) float->int 
// conversions as being slow.
//

typedef double lreal;
typedef float  real;
typedef unsigned long uint32;
typedef long int32;

const lreal _double2fixmagic = 68719476736.0*1.5;     //2^36 * 1.5,  (52-_shiftamt=36) uses limited precisicion to floor
const int32 _shiftamt        = 16;                    //16.16 fixed point representation,

#if BigEndian_
	#define iexp_				0
	#define iman_				1
#else
	#define iexp_				1
	#define iman_				0
#endif //BigEndian_

//--------------------
// im_real2int(double)
//--------------------
inline int32 im_real2int(lreal val)
{
#if DEFAULT_CONVERSION
	return val;
#else
	val		= val + _double2fixmagic;
	return ((int32*)&val)[iman_] >> _shiftamt; 
#endif
}

//-------------------
// im_real2int(float)
//-------------------
inline int32 im_real2int(real val)
{
#if DEFAULT_CONVERSION
	return val;
#else
	return im_real2int ((lreal)val);
#endif
}


#endif 

/*-----------------------------------------------------------------------
		
	Inline Math library : Bounding box implementation (header)
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

/*------------------------------------------------------------------------------------
		
	ibox3 class

	------------
	Description
	------------
		
	Bounding box class based for imath library

	This class performs intersection testing, not "collision detection".
	
	(collisions and intersections are similar, and you need intersections
	 to implement collisions)
	
	-------
	General
	-------

		Everything is public in this class (like the other imath components)
		This class has some non inlined functions

	Orientation
	-----------

		To orient the bbox simply fiddle with
		the quaternion

		eg:		iquaternion4 somerotation;

				//create somerotation with a suitable rotation

				object.GetQuaternion().multiply(somerotation);



---------------------------------------------*/

#ifndef _IBOX3_HEADER_
#define _IBOX3_HEADER_

//#define _IBOX3_TESTWITH_TEAPOT

#ifdef _IBOX3_TESTWITH_TEAPOT
	#include <gl/glut.h>
#else
	#include <gl/gl.h>
	#include <gl/glu.h>
#endif


class ibox3
{
	public:
		
		//assigning
		
		const ibox3 &operator=(const ibox3 &rhs)
		{
			center = rhs.center;
			extents = rhs.extents;	
			pos = rhs.pos;
			quat = rhs.quat;
			bas = rhs.bas;

			return *this;
		}
		
		//------- construction etc ---
		//----------------------------
		
		//constructor
		ibox3() 
		{
			center = zero;
			extents.set(0.01f,0.01f,0.01f);
		}

		//default constructor				
		ibox3(ivector3 &cextents,ivector3 &rotationcenter) 
		{
			center = rotationcenter;
			extents = cextents;
		}

		//copy constructor
		ibox3(const ibox3 ©)
		{
			center = copy.center;
			extents = copy.extents;
			pos = copy.pos;
			quat = copy.quat;
			bas = copy.bas;
		}

		//destructor
		~ibox3()
		{

		}

		//----------------------------
		//----------------------------

		//------ public attributes ---
		//----------------------------
	
		//you can only change pos/angles directly
		//position, or center of box

		ivector3		pos;		//position of box
		imatrix16		rm;			//rotation matrix (made from quaternion as nessesary)
		iquaternion4	quat;		//

		//dont fuck with the basis outside this class (copy)
		inline ibasis9 getBasis()
		{
			return bas;
		}
		//ditto for extents(copy)
		inline ivector3 getExtents()
		{
			return extents;
		}
		//ditto for center(copy)
		inline ivector3 getCenter()
		{
			return center;
		}

		//----------------------------
		//----------------------------

		const bool intersect(const ibox3 &b) const;	
		
		//===========================================
		//inline bool intersectSphere(const ibox3 &b) const
		//
		//	true for bounding sphere collision detection
		//	between this and b.  optimised to no sqrts.
		//
		//
		//	when inlined this give as big speedup
		//  when we need to call it alot of times
		//  (its an quick exit from collision)
		//============================================*/
		//const bool intersectSphere(const ibox3 &b) const;

		inline bool intersectSphere(const ibox3 &b) const
		{

			//float aradsqr =					killed tempories
			//float bradsqr = 
			//float radsqr = aradsqr + bradsqr;

			//float arad = extents.length();
			//float brad = b.extents.length();
			//float rad = arad + brad;

			ivector3 areal = pos;
			ivector3 breal = b.pos;

			ivector3 temp;

			//because bas can be anything
			//there doesnt seem to be a
			//good way to simplfy here

			temp = bas.x;temp.multiply(center.x);areal.add(temp);
			temp = bas.y;temp.multiply(center.y);areal.add(temp);
			temp = bas.z;temp.multiply(center.z);areal.add(temp);

			temp = b.bas.x;temp.multiply(b.center.x);breal.add(temp);
			temp = b.bas.y;temp.multiply(b.center.y);breal.add(temp);
			temp = b.bas.z;temp.multiply(b.center.z);breal.add(temp);

			temp = breal - areal;

			//as with most of my maths, i dont know why we must *2
			//here. its because the two lengthsqr's on the one side are
			//too small
			if ( temp.lengthsqr() < (extents.lengthsqr() + b.extents.lengthsqr() ) *2 ) 
				return true;
			else
				return false;
			

		}

		const void gldraw();			//draws the oobb (bounding box) (shows what will/should be intersected with)
		const void gldrawSphere();		//draws the bounding sphere

		const void update()
		{
			rm.makefromQuaternion(quat);
			bas = rm.getBasis();
		}

		protected:
		
			ibasis9			bas;		//ortho basis
			ivector3		extents;	//the extents of the box are constant
			ivector3		center;		//center offset
};

#endif 

/*-----------------------------------------------------------------------
		
	Inline Math library : Bounding box implementation (cpp)
	
	* See imath_box.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

#include "imath.h"

/*===========================================
const void ibox3::gldrawSphere()

//very dirty way of display the box
============================================*/

const void ibox3::gldrawSphere()
{
	GLUquadricObj* qobj;
	qobj = gluNewQuadric();

	gluQuadricDrawStyle(qobj,GLU_SILHOUETTE);

	ivector3 areal = pos;
	ivector3 temp;

	temp = bas.x;temp.multiply(center.x);areal.add(temp);
	temp = bas.y;temp.multiply(center.y);areal.add(temp);
	temp = bas.z;temp.multiply(center.z);areal.add(temp);


	glPushMatrix();
		
		glTranslatef(areal.x,areal.y,areal.z);
		gluSphere(qobj, extents.length() ,9,9);

	glPopMatrix();

	gluDeleteQuadric(qobj);
}

/*===========================================
const void ibox3::gldraw()

//very dirty way of display the box
============================================*/
const void ibox3::gldraw()
{
	//pos.gldraw();
	
	//draw the box

	ivector3 h,l;

	h.add(extents);
	l.subtract(extents);

	//could do my own rotation around it
	//from here.

	glPushMatrix();
		
		glTranslatef(pos.x,pos.y,pos.z);

		update();

		//push and add the quaternion matrix
		glPushMatrix();
		
			glMultMatrixf(rm.m);

				glTranslatef(center.x,center.y,center.z);
				
//draw test object here
#ifdef _IBOX3_TESTWITH_TEAPOT
					glEnable(GL_LIGHTING);	
					glRotatef(90,1.0,0.0,0.0);
//draw test object here
					glDisable(GL_LIGHTING);	
#endif	
					glDisable(GL_TEXTURE_2D);	

					//draw the bow
					glBegin(GL_LINE_LOOP);							//top loop
 						glVertex3f(h.v[0],h.v[1],h.v[2]);	
						glVertex3f(l.v[0],h.v[1],h.v[2]);
						glVertex3f(l.v[0],l.v[1],h.v[2]);
						glVertex3f(h.v[0],l.v[1],h.v[2]);
					glEnd();
					glBegin(GL_LINE_LOOP);							//bottom loop
						glVertex3f(h.v[0],h.v[1],l.v[2]);
						glVertex3f(l.v[0],h.v[1],l.v[2]);
						glVertex3f(l.v[0],l.v[1],l.v[2]);
						glVertex3f(h.v[0],l.v[1],l.v[2]);
					glEnd();
					glBegin(GL_LINES);									//the lines inbetween the loops
						glVertex3f(h.v[0],h.v[1],h.v[2]);glVertex3f(h.v[0],h.v[1],l.v[2]);
						glVertex3f(l.v[0],h.v[1],h.v[2]);glVertex3f(l.v[0],h.v[1],l.v[2]);
						glVertex3f(l.v[0],l.v[1],h.v[2]);glVertex3f(l.v[0],l.v[1],l.v[2]);
						glVertex3f(h.v[0],l.v[1],h.v[2]);glVertex3f(h.v[0],l.v[1],l.v[2]);
					glEnd();
					
					glBegin(GL_LINES);
						glColor3f(1.0f,0.0,0.0);glVertex3f(-0.5f,0.0,0.0);glVertex3f(1.2f,0.0,0.0);
						glColor3f(0.0,1.0f,0.0);glVertex3f(0.0,-0.5f,0.0);glVertex3f(0.0,1.2f,0.0);
						glColor3f(0.0,0.0,1.0f);glVertex3f(0.0,0.0,-0.5f);glVertex3f(0.0,0.0,1.2f);
					glEnd();

		glPopMatrix();

		//pop off the quat matrix
		//and draw the basis
	
		bas.gldraw();

		glBegin(GL_LINES);
			glColor3f(1.0f,0.0,0.0);glVertex3f(-0.5f,0.0,0.0);glVertex3f(1.2f,0.0,0.0);
			glColor3f(0.0,1.0f,0.0);glVertex3f(0.0,-0.5f,0.0);glVertex3f(0.0,1.2f,0.0);
			glColor3f(0.0,0.0,1.0f);glVertex3f(0.0,0.0,-0.5f);glVertex3f(0.0,0.0,1.2f);
		glEnd();
	glPopMatrix();	

}

/*===========================================
const bool ibox3::intersect(const ibox3 &b)

	true for bounding box collision detection
	between this and b
	
    (seperating axis test based on
	 gamasutra article)
============================================*/

const bool ibox3::intersect(const ibox3 &b) const
{
	
	//<mat2002>
	
	//mat2002 added offset for center
	//
	//we must pass the real center to the intersect test, adjusted for the pivot/center
	//to do this, we need to transform the position by the center with respect to the basis
	//(this just means moving it from its position to its center along its basis)

	ivector3 areal = pos;
	ivector3 breal = b.pos;

	ivector3 temp;

	temp = bas.x;temp.multiply(center.x);areal.add(temp);
	temp = bas.y;temp.multiply(center.y);areal.add(temp);
	temp = bas.z;temp.multiply(center.z);areal.add(temp);

	temp = b.bas.x;temp.multiply(b.center.x);breal.add(temp);
	temp = b.bas.y;temp.multiply(b.center.y);breal.add(temp);
	temp = b.bas.z;temp.multiply(b.center.z);breal.add(temp);

	//works without centers
	//ivector3 areal = (pos);
	//ivector3 breal = (b.pos);

	ivector3 v =  breal - areal;								//translation, in parent frame
	ivector3 T( bas.x.dot(v),bas.y.dot(v),bas.z.dot(v) );		//translation in A's frame

	float R[3][3];												//rotation matrix
	float ra,rb,t;

	int i,k;

	//calculate rotation matrix R[][]
	R[0][0] = bas.x.dot(b.bas.x);												
	R[0][1] = bas.x.dot(b.bas.y);
	R[0][2] = bas.x.dot(b.bas.z);

	R[1][0] = bas.y.dot(b.bas.x);												
	R[1][1] = bas.y.dot(b.bas.y);
	R[1][2] = bas.y.dot(b.bas.z);

	R[2][0] = bas.z.dot(b.bas.x);												
	R[2][1] = bas.z.dot(b.bas.y);
	R[2][2] = bas.z.dot(b.bas.z);



//	for (i=0;i<3;i++)
//		for (k=0;k<3;k++)
//			printf("R[%i][%i] %f\n",i,k,R[i][k]);


	//a's basis vectors
	for (i=0;i<3;i++)
	{
		ra = extents[i];		
		//ra = extents[i];     //mat2002

		rb = b.extents[0] * fabs(R[i][0]) + b.extents[1] * fabs(R[i][1]) + b.extents[2] * fabs(R[i][2]);
		t = fabs(T[i]);
		if (t > ra + rb)
		{
			#ifdef VERBOSE_SEPERATION
				printf("a basis %i seperates a and b\n",i);
			#endif

			return false;
		}
	}

	//safe to here

	//b's basis vectors 
	for (k=0;k<3;k++)
	{
		ra = extents[0] * fabs(R[0][k]) + extents[1] * fabs(R[1][k]) + extents[2] * fabs(R[2][k]);
		rb = b.extents[k];
		t = fabs(T[0] * R[0][k] + T[1] * R[1][k] + T[2] * R[2][k]);
		if (t > ra + rb) 
		{
			#ifdef VERBOSE_SEPERATION
				printf("b basis %i seperates a and b\n",k);
			#endif

			return false;
		}
	}
	
	//9 cross products
	
	//all the a0s
	
	ra = extents[1] * fabs(R[2][0]) + extents[2] * fabs(R[1][0]);
	rb = b.extents[1] * fabs(R[0][2]) + b.extents[2] * fabs(R[0][1]);
	t = fabs(T[2] * R[1][0] - T[1] * R[2][0]);
	if (t > ra + rb) 
	{
		#ifdef VERBOSE_SEPERATION
			printf("a0*b0 seperates a and b\n",k);
		#endif
		
		return false;
	}
	ra = extents[1] * fabs(R[2][1]) + extents[2] * fabs(R[1][1]);
	rb = b.extents[0] * fabs(R[0][2]) + b.extents[2] * fabs(R[0][0]);
	t = fabs(T[2] * R[1][1] - T[1] * R[2][1]);
	if (t > ra + rb)
	{

		#ifdef VERBOSE_SEPERATION
			printf("a0*b1 seperates a and b\n",k);
		#endif

		return false;
	}

	ra = extents[1] * fabs(R[2][2]) + extents[2] * fabs(R[1][2]);
	rb = b.extents[0] * fabs(R[0][1]) + b.extents[1] * fabs(R[0][0]);
	t = fabs(T[2] * R[1][2] - T[1] * R[2][2]);
	if (t > ra + rb) 
	{

		#ifdef VERBOSE_SEPERATION
			printf("a0*b2 seperates a and b\n",k);
		#endif

		return false;
	}


	//all the a1s

	ra = extents[0] * fabs(R[2][0]) + extents[2] * fabs(R[0][0]);
	rb = b.extents[1] * fabs(R[1][2]) + b.extents[2] * fabs(R[1][1]);
	t = fabs(T[0] * R[2][0] - T[2] * R[0][0]);
	if (t > ra + rb) 
	{

		#ifdef VERBOSE_SEPERATION
			printf("a1*b0 seperates a and b\n",k);
		#endif

		return false;
	}

	ra = extents[0] * fabs(R[2][1]) + extents[2] * fabs(R[0][1]);
	rb = b.extents[0] * fabs(R[1][2]) + b.extents[2] * fabs(R[1][0]);
	t = fabs(T[0] * R[2][1] - T[2] * R[0][1]);
	if (t > ra + rb) 
	{
		#ifdef VERBOSE_SEPERATION
			printf("a1*b1 seperates a and b\n",k);
		#endif
			
		return false;
	}

	ra = extents[0] * fabs(R[2][2]) + extents[2] * fabs(R[0][2]);
	rb = b.extents[0] * fabs(R[1][1]) + b.extents[1] * fabs(R[1][0]);
	t = fabs(T[0] * R[2][2] - T[2] * R[0][2]);
	if (t > ra + rb) 
	{
		#ifdef VERBOSE_SEPERATION
			printf("a1*b2 seperates a and b\n",k);
		#endif
			
		return false;
	}


	//all the a2s


	ra = extents[0] * fabs(R[1][0]) + extents[1] * fabs(R[0][0]);
	rb = b.extents[1] * fabs(R[2][2]) + b.extents[2] * fabs(R[2][1]);
	t = fabs(T[1] * R[0][0] - T[0] * R[1][0]);
	if (t > ra + rb) 
	{
		#ifdef VERBOSE_SEPERATION
			printf("a2*b0 seperates a and b\n",k);	
		#endif
			
		return false;
	}

	ra = extents[0] * fabs(R[1][1]) + extents[1] * fabs(R[0][1]);
	rb = b.extents[0] * fabs(R[2][2]) + b.extents[2] * fabs(R[2][0]);
	t = fabs(T[1] * R[0][1] - T[0] * R[1][1]);
	if (t > ra + rb)
	{
		#ifdef VERBOSE_SEPERATION
			printf("a2*b1 seperates a and b\n",k);	
		#endif
			
		return false;
	}

	ra = extents[0] * fabs(R[1][2]) + extents[1] * fabs(R[0][2]);
	rb = b.extents[0] * fabs(R[2][1]) + b.extents[1] * fabs(R[2][0]);
	t = fabs(T[1] * R[0][2] - T[0] * R[1][2]);
	if (t > ra + rb) 
	{
		#ifdef VERBOSE_SEPERATION
			printf("a2*b2 seperates a and b\n",k);
		#endif
			
		return false;
	}

	//disturbed laughter
	return true;
} 

/*-----------------------------------------------------------------------
		
	Inline Math library : Ray and Plane implementation
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/

#ifndef _IMATH_TRACE_IMP_HEADER_529081235
#define _IMATH_TRACE_IMP_HEADER_529081235

/*---------------------------------------------
	
	Raytracing/linetracing and collision detection

	iray6
	iplane6

	-------
	classes
	-------
		
		iray6
		iplane6

	------
	design
	------

		* a ray is defined a point(position) and a normalised vector(direction)
		* a plane is defined as abcd
	
---------------------------------------------*/

//-------------- iplane4 ------------
//	plane class
//------------------------------------

class iplane4
{
	public:
												
		ivector3 norm;
		float d;
		
		//-------------- operators ---
		//---------------------------- 
	    //return types according to Effective C++ 2e (see ref 8)

		//assignment
		inline const iplane4 &operator=(const iplane4 &right)
		{
			if (&right != this) {				//check for self assign
				norm = right.norm;				//copy right to left					
				d = right.d;		
			}
			return *this;						//enables x = y = z						
		}

		//---------------------------- 
		//---------------------------- 

		//------- construction etc ---
		//----------------------------
		
		//inlined ctor			
		inline iplane4()						   
		{
			//no need here
			//vectors construct themselves adequetly
		}					
		
		//inlined ctor	
		inline iplane4(const float in_d,const ivector3 &in_norm)
		{
			norm = in_norm;
			d = in_d;
		}					
		
		//inlined dtor(unnessesary)
		inline ~iplane4()
		{}	
		
		//---------------------------- 
		//---------------------------- 

		//--------- functions --------
		//----------------------------
		//sets components
													
		//dumps the componets (debug aid)
		inline void cdump(const char *caption)  const 
		{
			norm.cdump("iplane4 norm = ");
			printf("iplane4 d = %f\n",d);
		}	
		
		//draw's the plane (debug aid)
		//plane extend to infinity,
		//so the client decideds how far that is (length)
		//(not optimised for drawing lots of planes, which may be needed later)

		//creates a quad from the plane's normal
		//much like a particle is generated7
		inline void gldraw(float length,ivector3 rgbcolor)			
		{
			//some good practice at vector math here
			
			//we need to create an up vector and a left vector
			//from the normal of the plane

			//both are perpendicular to the norm, so crossproducts
			//are an easy way to do this (allthought there may
			//be a faster solution)

			//create any vector perp to the normal
			//(the operand xaxis could be any vector)
			ivector3 left = norm.cross(xaxis);
			left.normalise();
			
			//the up vector needs to be perp to both the norm and the left vector
			ivector3 up =  norm.cross(left);
			up.normalise();
			
			//size it according to length param
			left.multiply(length);
			up.multiply(length);

			//pos is a point on the plane
			ivector3 pos = (norm * d);
			ivector3 temp = pos;

			glDisable(GL_LIGHTING);

			glBegin(GL_LINES);

				glColor3fv(red.v);

				pos.glplot();temp = pos + norm;
				temp.glplot();		
	
				glColor3fv(cyan.v);
			
				pos.glplot();temp = pos + left;
				temp.glplot();

				glColor3fv(yellow.v);

				pos.glplot();temp = pos + up;
				temp.glplot();

			glEnd();

			glColor3fv(rgbcolor.v);

			glBegin(GL_QUADS);
				glVertex3f(pos.x + left.x + up.x,pos.y + left.y + up.y,pos.z + left.z + up.z);
				glVertex3f(pos.x - left.x + up.x,pos.y - left.y + up.y,pos.z - left.z + up.z);
				glVertex3f(pos.x - left.x - up.x,pos.y - left.y - up.y,pos.z - left.z - up.z);
				glVertex3f(pos.x + left.x - up.x,pos.y + left.y - up.y,pos.z + left.z - up.z);
			glEnd();
		
		}		
 
		//---------------------------- 
		//----------------------------

};

//-------------- iray ------------
//
//	ray class
//------------------------------------

class iray6
{
	public:
												
		ivector3 pos;		//treated as a position
		ivector3 dir;		//treated as a direction (magnitude represents light/strength?)
		

		//-------------- operators ---
		//---------------------------- 
	    //return types according to Effective C++ 2e (see ref 8)

		//assignment
		//the comilier should be able to generate this (its only member copy)
		inline const iray6 &operator=(const iray6 &right)
		{
			if (&right != this) {				//check for self assign
				dir = right.dir;				//copy right to left					
				pos = right.pos;		
			}
			return *this;						//enables x = y = z						
		}

		//---------------------------- 
		//---------------------------- 

		//------- construction etc ---
		//----------------------------
		
		//inlined ctor			
		inline iray6()						   
		{
			pos.set(0,0,0);			//we'll start all rays at 0,0,0
			dir.set(0,0,1.0);		//pointing up in the air
		}					
		
		//inlined ctor(initalises to pos,dir)	
		inline iray6(const ivector3 &in_pos,const ivector3 &in_dir)
		{
			pos = in_pos;
			dir = in_dir;
		}	
		
		//inlined cpy ctor
		//the comilier generates this (its only member copy)
		
		//inlined dtor(unnessesary)
		inline ~iray6()
		{}	
		
		//---------------------------- 
		//---------------------------- 

		//--------- functions --------
		//----------------------------

		inline float intersect(class iplane4 &with)
		{
			
			float cosAlpha = dir.dot( with.norm );		//dot product of ray's direction vector (normalised) and plane's normal vector (normalised)
		
			ivector3 diff = (with.norm * with.d) - pos;				//vector from plane origin to ray origin
			float deltaD = diff.dot(with.norm);				//dot product of above and plane's normal
			
			//deltaD is the length of the shortest line to the point 

			return deltaD/cosAlpha;
		}


		//dumps the componets (debug aid)
		inline void cdump(const char *caption)  const 
		{
			pos.cdump("ray pos = ");
			dir.cdump("ray dir = ");
		}	
		
		//draw's the ray (debug aid)
		//rays extend to infinity,
		//so the client decideds how far that is (length)

		//(not optimised for drawing lots of rays, which may be needed later)		
		inline void gldraw(float length,const ivector3 &startcol,const ivector3 &endcol) const			
		{
			ivector3 endpos = dir;
				
			endpos.multiply(length);
			endpos.add(pos);

			glDisable(GL_LIGHTING);		
			//glDisable(GL_BLEND);
			glDisable(GL_TEXTURE_2D);
			glBegin(GL_LINES);
				glColor3fv(startcol.v);
				pos.glplot();
				glColor3fv(endcol.v);
				endpos.glplot();
			glEnd();
		}	
		
		inline void gldraw(float length) const			
		{
			ivector3 endpos = dir;
			endpos.multiply(length);
			endpos.add(pos);
			glDisable(GL_LIGHTING);
			glBegin(GL_LINES);
				pos.glplot();
				endpos.glplot();
			glEnd();
		}	
 
		//---------------------------- 
		//----------------------------

};



#endif 

/*-----------------------------------------------------------------------
		
	Inline Math library : Vector Matrix and Quaternion implementation
	
	* See imath.h for more information

	By Matt Barnett
		(see References in imath.h)
	
------------------------------------------------------------------------*/


#ifndef _IMATH_V_IMP_HEADER_4778097634v
#define _IMATH_V_IMP_HEADER_4778097634v

//-------------- ivector3 ------------
//
//  simple 3d vector with 3 components
//	(there is no 2d version, use this)
//------------------------------------

//enum for auto constructing some vector types
//if you want to contruct a random vector
//or  and random normalised vector etc.
enum ivector3_constructor_enum { IVECTOR3_RANDOM , IVECTOR3_RANDOM_NORM };

class ivector3
{
	public:
			
		//three floats are the data										
		union
		{
			float v[3];
			struct { float x, y, z; };
		};

		//-------------- operators ---
		//---------------------------- 
	    //return types according to Effective C++ 2e (see ref 8)
	
		inline const ivector3 operator-(const ivector3 &right) const 
		{
			return ivector3(x - right.x,y - right.y,z - right.z);
		}
		inline const ivector3 operator+(const ivector3 &right) const
		{
			return ivector3(x + right.x,y + right.y,z + right.z);
		}	
		
		inline const ivector3 operator*(const ivector3 &right) const 
		{
			return ivector3(x * right.x,y * right.y,z * right.z);
		}
		inline const ivector3 operator/(const ivector3 &right) const
		{
			return ivector3(x / right.x,y / right.y,z / right.z);
		}

		inline const ivector3 operator*(const float scalar) const 
		{
			return ivector3(x * scalar,y * scalar,z * scalar);
		}
		inline const ivector3 operator/(const float scalar) const 
		{
			return ivector3(x / scalar,y / scalar,z / scalar);
		}

		//void operators
		inline void operator-=(const ivector3 &right) 
		{
			subtract(right);
		}
		inline void operator+=(const ivector3 &right)
		{
			add(right);
		}

		//assignment
		inline ivector3 &operator=(const ivector3 &right)
		{
			if (&right != this) {				//check for self assign
				x = right.x;					//copy right to left					
				y = right.y;	
				z = right.z;	
			}
			return *this;						//enables x = y = z						
		}

		inline float operator[](int subscript) const
		{ 
			assert( subscript < 3 && subscript >= 0 );
			return v[subscript];
		}
		//---------------------------- 
		//---------------------------- 

		//------- construction etc ---
		//----------------------------
		
		//inlined ctor(initalises to 0,0,0)				
		inline ivector3()						   
		{
			set(0,0,0); 
		}					
		
		//inlined ctor(initalises to x,y,z)	
		inline ivector3(const float setx,const float sety,const float setz)
		{
			set(setx,sety,setz); 
		}	
		
		//ctor for a random normalised vector
		explicit inline ivector3(ivector3_constructor_enum inv_ce)
		{
			if (inv_ce == IVECTOR3_RANDOM_NORM)
			{
				randomize(1.0);
				normalise();
			}
			else if (inv_ce == IVECTOR3_RANDOM)
				randomize(1.0);
		}						

		//inlined dtor(unnessesary)
		inline ~ivector3()
		{}			
		//---------------------------- 
		//---------------------------- 

		//--------- void functions ---
		//----------------------------
		//sets vector components
		
		inline void set(const float setx,const float sety,const float setz)
		{
			x = setx;
			y = sety;
			z = setz;
		}	
		
		//sets z == 0
		inline void make2d()						  
		{ 
			z = 0;
		}															

		//dumps the vectors componets (debug aid)
		inline void cdump(const char *caption)  const 
		{
			cout << caption << x << " " << y << " " << z << endl;
		}	
		
		//draw O -> vector (debug aid)
		inline void gldraw() const			
		{
			glDisable(GL_LIGHTING);
			glBegin(GL_LINES);
			glVertex3f(0.0,0.0,0.0);
			glVertex3fv(v);
			glEnd();
		}														
		//plot the vector (debug aid)
		inline void glplot() const			
		{
			glVertex3fv(v);
		}	

		//multiplies this by scalar
		inline void multiply(const float scalar)			  
		{
			x *= scalar;
			y *= scalar;
			z *= scalar;
		}	
		
		//multiplies this by multiply
		inline void multiply(const ivector3 &multiply)   
		{ 
			x *= multiply.x;
			y *= multiply.y;
			z *= multiply.z; 
		}

		//divides this by scalar
    	inline void divide(const float scalar)			  
		{ 
			x /= scalar; 
			y /= scalar;
			z /= scalar;
		}
		
		//divides this by divide
    	inline void divide(const ivector3 ÷)			  
		{ 
			x /= divide.x;
			y /= divide.y;
			z /= divide.z; 
		}					
	
		//adds add to this		
		inline void add(const ivector3 &add)		
		{ 
			x += add.x;
			y += add.y;
			z += add.z; 
		}			
		
		//subtracts subtract from this
		inline void subtract(const ivector3 &subtract)	 	
		{
			x -= subtract.x;
			y -= subtract.y; 
			z -= subtract.z; 
		}

		//multiplies this by matrix multiply
		//inline void multiply(imatrix16 &multiply)  
		//{ 
		//
		//}	
		
		//normalise
		inline void normalise()			 
		{ 
			const float n = length();
			if (!n)
			{
				#ifdef _IVECTOR3_VERBOSE_NORMALISE_FAIL_
					printf("ivector3::MATH_ERROR ( normalise failed(), length = 0 )\n");
					cdump("vector ");
				#endif

				return;
			}

			x /= n;
			y /= n;
			z /= n; 
		}	
			
		//randomises each component
		inline void randomize(float range)
		{
			x = im_random_clamp(range);
			y = im_random_clamp(range);
			z = im_random_clamp(range);
		}
		
		//makes "this" the normal of the three vectors
		inline void makeTriangleNormal(const ivector3 &a,const ivector3 &b,const ivector3 &c)
		{
			ivector3 ab,ac;

			ab = b - a;
			ac = c - a;

			ivector3 cross = ab.cross(ac);
			cross.normalise();

			*this = cross;
			
		}
		//---------------------------- 
		//---------------------------- 
		
		//-------------- functions ---
		//----------------------------
		//crossproduct
		inline const ivector3 &cross(const ivector3 &in) const
		{
			ivector3 out;
			out.x = y*in.z - z*in.y; 
			out.y = z*in.x - x*in.z;
			out.z = x*in.y - y*in.x;
			return out;
		}
		
		//length/norm
		inline float length() const
		{
			float a = (x*x)+(y*y)+(z*z);

			if (a <= 0)
			{
				#ifdef _IVECTOR3_VERBOSE_LENGTH_FAIL_
					printf("ivector3::MATH_ERROR ( length failed, x*x+y*y+z*z <= 0 )\n");
					cdump("vector ");
					printf("x*x+y*y+z*z = %f\n",a);
				#endif

				return 0;
			}

			return im_sqrt(a); 
		}
		
		//length/norm
		inline float lengthsqr() const
		{
			float a = (x*x)+(y*y)+(z*z);

			if (a <= 0)
			{
				#ifdef _IVECTOR3_VERBOSE_LENGTH_FAIL_
					printf("ivector3::MATH_ERROR ( length failed, x*x+y*y+z*z <= 0 )\n");
					cdump("vector ");
					printf("x*x+y*y+z*z = %f\n",a);
				#endif

				return 0;
			}

			return a; 
		}
		
		//dotproduct
		inline float dot(const ivector3 &in) const
		{
			return x*in.x + y*in.y + z*in.z;
		}			
		
		//angle between this and in
		//
		//inline float angle(const ivector3 &in) const				
		//{	
		//	float tdotin,theta;
		//
		//	tdotin = dot(in);								//dot product
		//	
		//	float li;
		//	li = in.length();
		//
		//	theta = tdotin / (length() * li);				//lengths
		//	theta = (float)im_acos(theta);					//inverse im_cos
		//
		//	return theta;									//return theta in radians
		//}

		inline float cangle(const ivector3 &in) const
		{
			return ( x*in.x + y*in.y + z*in.z ) / ( (length() )*(in.length() ) );
		}

		//matrix based angle axis rotation (radians)
		//matrix isnt stored
		//usefull to rotate a point/vector about an axis
		inline void rotate(const float angle,const ivector3 &axis)
		{	
			//componets of rotation axis
			//float a = axis.x;		
			//float b = axis.y;
			//float c = axis.z;

			//get the matrix
			float f1cos,fasin,fbsin,fcsin;
			
			//fcos = im_cos(angle);		//im_cos of angle
			//fsin = im_sin(angle);		//im_sin of angle
			
			f1cos = (1 - im_cos(angle) );			//1 - im_cos of angle
			
			fasin = (axis.x * im_sin(angle) );			//a * im_sin of angle
			fbsin = (axis.y * im_sin(angle) );			//b * im_sin of angle
			fcsin = (axis.z * im_sin(angle) );			//c * im_sin of angle

			//m[0] = a*a*f1cos + fcos;
			//m[1] = a*b*f1cos - fcsin;
			//m[2] = a*c*f1cos + fbsin;

			//m[3] = a*b*f1cos + fcsin;
			//m[4] = b*b*f1cos + fcos;
			//m[5] = b*c*f1cos - fasin;

			//m[6] = a*c*f1cos - fbsin;
			//m[7] = b*c*f1cos + fasin;
			//m[8] = c*c*f1cos + fcos;

			ivector3 out;

			out.x = x * (axis.x * axis.x * f1cos + im_cos(angle) )  + y * (axis.x * axis.y * f1cos - fcsin) + z * (axis.x * axis.z * f1cos + fbsin);
			out.y = x * (axis.x * axis.y * f1cos + fcsin) + y * (axis.y * axis.y * f1cos + im_cos(angle) ) + z * (axis.y * axis.z * f1cos - fasin);
			out.z = x * (axis.x * axis.z * f1cos - fbsin) + y * (axis.y * axis.z * f1cos + fasin) + z * (axis.z * axis.z * f1cos + im_cos(angle) );

			x = out.x;
			y = out.y;
			z = out.z;

		}
		//---------------------------- 
		//---------------------------- 
};

//------ ivector3 constants -----
//-------------------------------

const ivector3 xaxis(1.0,0.0,0.0);
const ivector3 yaxis(0.0,1.0,0.0);
const ivector3 zaxis(0.0,0.0,1.0);
const ivector3 zero (0.0,0.0,0.0);

//a few colors
const ivector3 red  (1.0,0.0,0.0);  
const ivector3 green(0.0,1.0,0.0);
const ivector3 blue (0.0,0.0,1.0);
const ivector3 black(0.0,0.0,0.0);

const ivector3 cyan   (0.0,1.0,1.0);
const ivector3 magenta(1.0,0.0,1.0);
const ivector3 yellow (1.0,1.0,0.0);

//----------------------------
//----------------------------

//-------------- ibasis9 ------
//
// just for storing and drawing
// orthonominal(??) bases or axes
//
// there is no clever shit here making
// sure that the vectors remain perpendicular
// this is just for storage of 3 vectors.

class ibasis9
{
	public:	
		
		//three vectors are the data
		ivector3 x,y,z;

		ibasis9::ibasis9(const ivector3 &inx,const ivector3 &iny,const ivector3 &inz)
		{
			x = inx;
			y = iny;
			z = inz;
		}
		
		ibasis9::ibasis9()
		{
			x = xaxis;
			y = yaxis;
			z = zaxis;
		}

		//just normalises each vector to save some typing
		inline void normalise()
		{
			x.normalise();
			y.normalise();
			z.normalise();
		}

		//draws the basis (debug aid)
		inline void gldraw() const
		{
			glDisable(GL_LIGHTING);
			glBegin(GL_LINES);
			glColor3fv(cyan.v);   glVertex3fv(zero.v);glVertex3fv(x.v);
			glColor3fv(magenta.v);glVertex3fv(zero.v);glVertex3fv(y.v);
			glColor3fv(yellow.v); glVertex3fv(zero.v);glVertex3fv(z.v);
			glEnd();
		}
};

//---------------------------- 
//---------------------------- 

#endif