// written by André Betz 
// http://www.andrebetz.de

#include <windows.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.14279

class SOMClassifier
{
public:
	~SOMClassifier();
	SOMClassifier(	unsigned long* ulSite,unsigned long ulVecLen,unsigned long ulDim,
					unsigned long ulMax,unsigned long ulMin,double dDist,double dEta,
					double dDistDec, double dEtaDec);
    unsigned long Step(double* pdPat);
	double*			m_pdMap;
protected:
	double			m_dDist;
	double			m_dEta;
	double			m_dDistDec;
	double			m_dEtaDec;
	unsigned long*	m_pulSite;
	unsigned long	m_ulDim;
	unsigned long	m_ulVecLen;
	unsigned long	m_ulNeurons;
	unsigned long	m_ulMin;
	unsigned long	m_ulMax;
	unsigned long	m_ulNghbFunc;
			
	unsigned long Learn(double* pdDataPtr,unsigned long ulWinPos);
	unsigned long MinDistance(double* pdDataPtr,unsigned long* pulWinPos);  
	unsigned long Nghb(unsigned long ulWinPos, unsigned long ulNeuPos, double* pdDist);
	unsigned long L2Norm(double* pdVec1, double* pdVec2, double* pdDistance);
};

SOMClassifier::~SOMClassifier()
{
	delete m_pulSite;
	delete m_pdMap;
}

SOMClassifier::	SOMClassifier(unsigned long* ulSite,unsigned long ulVecLen,unsigned long ulDim,
							  unsigned long ulMax,unsigned long ulMin,double dDist,double dEta,
							  double dDistDec, double dEtaDec)
{
	unsigned long	ulCountNeu;

	m_ulNghbFunc= 3;
	m_ulMin		= ulMin;
	m_ulMax		= ulMax;
	m_dEta		= dEta;
	m_dDist		= dDist;
	m_ulDim 	= ulDim;
	m_ulVecLen	= ulVecLen;
	m_dDistDec	= dDistDec;
	m_dEtaDec	= dEtaDec;
	m_pulSite	= new unsigned long[m_ulDim];

	for(ulCountNeu=0;ulCountNeu<m_ulDim;ulCountNeu++)
	{	
		m_pulSite[ulCountNeu] = ulSite[ulCountNeu];
	
		if(ulCountNeu)
		{
			m_ulNeurons	*=	m_pulSite[ulCountNeu];	
		}
		else
		{
			m_ulNeurons = m_pulSite[ulCountNeu];
		}
	}

	m_pdMap = new double[m_ulNeurons*m_ulVecLen];

	for(ulCountNeu=0;ulCountNeu<(m_ulNeurons*ulVecLen);ulCountNeu++)
	{
		m_pdMap[ulCountNeu] = (double)(ulMin + rand()%ulMax);
	}
}

unsigned long SOMClassifier::Step(double* pdPat)
{
	unsigned long	ulWinPos;
	
    MinDistance(pdPat,&ulWinPos);
	Learn(pdPat,ulWinPos);
	
	return ulWinPos;
}

unsigned long SOMClassifier::Learn(double* pdPat,unsigned long ulWinPos)
{
	unsigned long	ulCountVec;
	unsigned long	ulCountNeu;
	double			dDistance;

	m_dDist						*= m_dDistDec;
	if(m_dEta > 0.02) m_dEta	*= m_dEtaDec;

	for(ulCountNeu=0;ulCountNeu<m_ulNeurons;ulCountNeu++)
	{
		Nghb(ulWinPos,ulCountNeu,&dDistance);
		for(ulCountVec=0;ulCountVec<m_ulVecLen;ulCountVec++)
		{
			m_pdMap[ulCountNeu*m_ulVecLen+ulCountVec] += dDistance * m_dEta *
				(pdPat[ulCountVec] - m_pdMap[ulCountNeu*m_ulVecLen+ulCountVec]);
		}		
	}

	return 0;
}

unsigned long SOMClassifier::Nghb(unsigned long ulWinPos,unsigned long ulNeuPos,double* pdDist)
{
	unsigned long	ulCount;
	unsigned long	ulSiteNeu;
	unsigned long	ulSiteWin;
	unsigned long	ulSize;
	double			dTemp;

	*pdDist		= 0.0f;
	ulSiteWin	= ulWinPos;
	ulSiteNeu	= ulNeuPos;
	ulSize		= m_ulNeurons;

	for(ulCount=m_ulDim;ulCount>0;ulCount--)
	{
		ulSize /= m_pulSite[ulCount-1];

		*pdDist +=	(ulSiteWin / ulSize - ulSiteNeu / ulSize) *
					(ulSiteWin / ulSize - ulSiteNeu / ulSize);

		ulSiteWin %= ulSize;
		ulSiteNeu %= ulSize;
	}

	*pdDist = sqrt(*pdDist);

	if(m_ulNghbFunc==1)
	{
		if(*pdDist<m_dDist)	*pdDist = 1.0f;
		else				*pdDist = 0.0f;
	}
	else if(m_ulNghbFunc==2)
	{
		dTemp = (*pdDist / m_dDist)*(*pdDist / m_dDist);
		*pdDist = (1-dTemp)*exp(-dTemp);
	}
	else if(m_ulNghbFunc==3)
	{
		dTemp = (*pdDist / m_dDist)*(*pdDist / m_dDist);
		*pdDist = exp(-dTemp);
	}
	else if(m_ulNghbFunc==4)
	{
		dTemp = 1 - *pdDist / m_dDist;
		if(*pdDist<m_dDist)	*pdDist = dTemp;
		else				*pdDist = 0.0f;
	}
	else if(m_ulNghbFunc==5)
	{
		dTemp = cos((*pdDist/m_dDist)*(PI/2));
		if(*pdDist<m_dDist)	*pdDist = dTemp;
		else				*pdDist = 0.0f;
	}

	return 0;
}

unsigned long SOMClassifier::MinDistance(double* pdPat,unsigned long* pulNeuPos)
{
	unsigned long	ulCountNeu;
	double			dDistance;	
	double			dResult;

	for(ulCountNeu=0;ulCountNeu<m_ulNeurons;ulCountNeu++)
	{
		L2Norm(pdPat,&m_pdMap[ulCountNeu*m_ulVecLen],&dDistance);
		if((ulCountNeu==0)||(dDistance < dResult))
		{
			dResult		= dDistance;
			*pulNeuPos	= ulCountNeu;
		}
	}

	return 0;
}
			
unsigned long SOMClassifier::L2Norm(double* pdVec1, double* pdVec2, double* pdDistance)
{
	unsigned long	ulCountVec;

	*pdDistance = 0.0f;

	for(ulCountVec=0;ulCountVec<m_ulVecLen;ulCountVec++)
	{
		*pdDistance +=	(pdVec1[ulCountVec] - pdVec2[ulCountVec]) *
						(pdVec1[ulCountVec] - pdVec2[ulCountVec]);
	}

	*pdDistance = sqrt(*pdDistance);

	return 0;
}


unsigned long	Sites[]	= {20,20};
unsigned long	VecLen	= 2;
unsigned long	Train	= 0;
unsigned long	TrShw	= 1;
unsigned long	wx		= 400;
unsigned long	wy		= 400;

SOMClassifier	som(Sites,VecLen,2,wx,0,10.0,0.15,1.0,1.0);

int Show(HWND hwnd)
{
	double			Vec[2];
	char			Number[] = "000000";
	unsigned long	x,y;
	HDC				hdcWindow;
	HBRUSH			brush;
	HRGN			hrgn;
	BOOL			t;

	Vec[0] = (double)(rand() % wx);
	Vec[1] = (double)(rand() % wy);

	som.Step(Vec);

	if(!(Train%TrShw))
	{
		hdcWindow = GetDC(hwnd);

		hrgn = CreateRectRgn(0,0,wx,wy);
		brush = CreateHatchBrush(HS_BDIAGONAL,0x00ffffff);
		t = FillRgn(hdcWindow,hrgn,brush);

		for(y=0;y<Sites[1];y++)
		{
			for(x=0;x<Sites[0];x++)
			{
				if(x<Sites[0]-1)
				{
					MoveToEx(hdcWindow,	(int)som.m_pdMap[(x+Sites[1]*y)		*VecLen],
										(int)som.m_pdMap[(x+Sites[1]*y)		*VecLen+1],NULL);
					LineTo(hdcWindow,	(int)som.m_pdMap[(x+Sites[1]*y+1)	*VecLen],
										(int)som.m_pdMap[(x+Sites[1]*y+1)	*VecLen+1]);
				}

				if(y<Sites[1]-1)
				{
					MoveToEx(hdcWindow,	(int)som.m_pdMap[(x+Sites[1]*y)		*VecLen],
										(int)som.m_pdMap[(x+Sites[1]*y)		*VecLen+1],NULL);
					LineTo(hdcWindow,	(int)som.m_pdMap[(x+Sites[1]*(y+1))	*VecLen],
										(int)som.m_pdMap[(x+Sites[1]*(y+1))	*VecLen+1]);
				}
			}
		}
	
		wsprintf(Number,"%d",Train);
		TextOut(hdcWindow,0,0,Number,strlen(Number));

		ReleaseDC(hwnd,hdcWindow);
	}

	Train++;

	return 0;
}

LRESULT CALLBACK WndProc(HWND hwnd,UINT Message,WPARAM wParam,LPARAM lParam)
{
	switch(Message)
	{
		case WM_CLOSE:
			DestroyWindow(hwnd);
			break;
		case WM_DESTROY:
			PostQuitMessage(0);
			break;
		case WM_TIMER:
			Show(hwnd);
			break;
 
		default:
			return DefWindowProc(hwnd, Message, wParam, lParam);
   }
   return 0;
}

int WINAPI WinMain(HINSTANCE hInstance,HINSTANCE hPrevInstance,LPSTR lpCmdLine,int nCmdShow)
{
	WNDCLASSEX		WndClass;
	HWND			hwnd;
	MSG				Msg;
	char			ClassName[] = "MyWindowClass";
	UINT			idTimer1	= 1;
	UINT			nTimerDelay = 1;

	srand((unsigned int)time(NULL));

	WndClass.cbSize        = sizeof(WNDCLASSEX);
	WndClass.style         = NULL;
	WndClass.lpfnWndProc   = WndProc;
	WndClass.cbClsExtra    = 0;
	WndClass.cbWndExtra    = 0;
	WndClass.hInstance     = hInstance;
	WndClass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
	WndClass.hCursor       = LoadCursor(NULL, IDC_ARROW);
	WndClass.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
	WndClass.lpszMenuName  = NULL;
	WndClass.lpszClassName = ClassName;
	WndClass.hIconSm       = LoadIcon(NULL, IDI_APPLICATION);

	RegisterClassEx(&WndClass);
	hwnd = CreateWindowEx(WS_EX_CLIENTEDGE,ClassName,"2D-SOM  Betz99",WS_OVERLAPPEDWINDOW,
					CW_USEDEFAULT,CW_USEDEFAULT,wx,wy,NULL,NULL,hInstance,NULL);

	ShowWindow(hwnd,nCmdShow);
	UpdateWindow(hwnd);

	SetTimer(hwnd,idTimer1,nTimerDelay,NULL);

	while(GetMessage(&Msg,NULL,0,0))
	{
		TranslateMessage(&Msg);
		DispatchMessage(&Msg);
	}

	KillTimer(hwnd, idTimer1);
	
	return Msg.wParam;
}