/****************************************************************************************************************
						S E L F - O R G A N I Z I N G   M A P   C L U S T U R I N G

	Author: Pramod Lakshmi Narasimha
			Image Processing and Neural Networks Lab,
			Electrical Engineering Department,
			University of Texas at Arlington.

	This Program requires the Cluster Center vectors file that was generated in Seqential Leader.
*****************************************************************************************************************/

#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>

int IX = 3;
int IY = 4009;
int IZ = 234;


void main()
{
	FILE *df, *wts;
	char Infile[14];
	int Ndim, N, Nc, Nv = 0;
	int i, j, k, p, t, It, Nit, Classid;
	float *x, T1, T2, a1, a2;
	double *meanip, *stdevip, **m, E, Nt, z, dmin;

	double slete(double, double);
	double rand1(int *, int *, int *);
	double d(float *, double *, int);

	puts("\nEnter the name of the input file : ");
	gets(Infile);
	puts("\nEnter the number of elements per vector (Ndim): ");
	scanf("%d", &Ndim);
	puts("\nEnter the number of elements to be clustered (N): ");
	scanf("%d", &N);
	puts("\nEnter the number of clusters (Nc): ");
	scanf("%d", &Nc);

	// Memory Allocation
	x = (float *) malloc(sizeof(float) * N);
	meanip = (double *) malloc(sizeof(double) * N);
	stdevip = (double *) malloc(sizeof(double) * N);
	m = (double **) malloc(sizeof(double *) * Nc);
	for(i = 0; i < Nc; i++)
		m[i] = (double *) malloc(sizeof(double) * N);

	// Initialization of the vectors
	for(i = 0; i < N; i++)
	{
		meanip[i] = 0.0;
		stdevip[i] = 0.0;
	}

	df = fopen(Infile, "r");
	if(df == NULL)
	{
		perror(Infile);
		exit(0);
	}

	while(!feof(df))
	{
		Nv++;
		for(i = 0; i < N; i++)
		{
			fscanf(df, "%f", &x[i]);
			meanip[i] += x[i];
			stdevip[i] += x[i] * x[i];
		}
		for(i = 0; i < Ndim-N; i++)
			fscanf(df, "%f", &x[i]);
	}
	fclose(df);

	for(i = 0; i < N; i++)
	{
		meanip[i] /= Nv;
		stdevip[i] = stdevip[i]/Nv - meanip[i]*meanip[i];
		stdevip[i] = sqrt(stdevip[i]);
	}

	df = fopen("ClusterCenter.txt", "r");

	if(df == NULL)
	{
		perror("ClusterCenter.txt");
		exit(1);
	}

	for(i = 0; i < Nc && !feof(df); i++)
	{
		for(j = 0; j < N; j++)
			fscanf(df, "%lf", &m[i][j]);
		fscanf(df, "%d", &Classid);
	}

	Nc = i;

	fclose(df);

	puts("\nEnter the number of clustering Iterations : ");
	scanf("%d", &Nit);

	T1 = (float)Nv*Nit/3;
	T2 = (float)Nv*Nit/10;
	a1 = (float)Nc/Nv;
	a2 = (float)Nc/10;

	wts = fopen("selfresult.txt", "w");

	if(wts == NULL)
	{
		perror("selfresult.txt");
		exit(1);
	}

	for(It = 1; It <= Nit; It++)
	{
		E = 0;
		df = fopen(Infile, "r");
		for(p = 1; p <= Nv; p++)
		{
			for(i = 0; i < N; i++)
				fscanf(df, "%f", &x[i]);

			t = p + (It-1)*Nv;

			Nt = a2 * exp(-t/T2);
			z = a1 * exp(-t/T1);

			dmin = d(x, m[0], N);
			k = 0;

			for(j = 0; j < Nc; j++)
			{
				if(d(x, m[j], N) < dmin)
				{
					k = j;
					dmin = d(x, m[j], N);
				}
			}

			E += dmin;

			for(j = 0; j < Nc; j++)
			{
				if(abs(k-j) < Nt)
				{
					for(i = 0; i < N; i++)
						m[j][i] += z*(x[i] - m[j][i]);
				}
			}

			for(i = 0; i < Ndim-N; i++)
				fscanf(df, "%f", &x[i]);
		}
		fclose(df);

		printf("Error in %d iteration = %f\n", It, E/Nv);
		fprintf(wts, "Error in %d iteration = %f\n", It, E/Nv);
	}

	fprintf(wts, "\n\nMean Vectors : \n");

	for(i = 0; i < Nc; i++)
	{
		for(j = 0; j < N; j++)
			fprintf(wts, "%f\t", m[i][j]);
		fprintf(wts, "\n\n");
	}
	fclose(wts);
}




double d(float *x, double *m, int N)
{
	int i;
	double d = 0;
	for(i = 0; i < N; i ++)
	{
		d += (x[i] - m[i]) * (x[i] - m[i]);
	}
	return (d);
}




/******************************************************************************

	Subroutine : Random no. Generator
		-the random nos. generated are uniformly distributed
		 between 0 and 1.
		-IX, IY, IZ are the SEEDS

*******************************************************************************/


double rand1(int *ix, int *iy, int *iz)
{
	int ixx, iyy, izz;
	double itemp;
	double temp;


	ixx=(*ix)/177;
	*ix=171*(*ix%177)-2*ixx;

	if(*ix < 0)
		*ix+=30269;

	iyy=(*iy)/176;
	*iy=176*(*iy%176)-2*iyy;

	if(*iy < 0)
		*iy+=30307;

	izz=(*iz)/178;
	*iz=170*(*iz%178)-2*izz;

	if(*iz < 0)
		*iz+=30323;

	
	temp=(double)(*ix)/30629.0+(double)(*iy)/30307.0+(double)(*iz)/30323.0;
	itemp=floor(temp);
	return (temp-itemp);

}




/****************************** Gaussian Random No. generator ******************************/


double slete(double std, double xmean)
{
	double rand1(int *, int *, int *);
	double PI = 3.1415926;
	return (xmean+std*cos(2*PI*rand1(&IX, &IY, &IZ))*sqrt(-2.0*log(rand1(&IX, &IY, &IZ))));
}
