/*************** Clustering, Bezdek's Algorithm *****************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX_N    150  /* maximum number of data points x1, x2,....,xN */
#define MAX_DIM  10 /* maximum dimension */
#define MAX_CLUST  30 /* maximum number of clusters */
#define STOP 0.0001  /* max. improvement limit for stopping */

void main()
{
int i, j, l, m, n, c, switch1, select, flag1, seed;
float x[MAX_N][MAX_DIM], u[MAX_CLUST][MAX_N], v[MAX_CLUST][MAX_DIM];
float d[MAX_CLUST][MAX_N], old_u[MAX_CLUST][MAX_N], min[MAX_N];
float sum1, sum2, max, delta, ratio, m_prime;
char file_name[20], c1;

FILE *fp1;


printf("Input data file name\n");
printf("Data should be stored as matrix\n");
printf("row  -> data point (vector) Xi: column -> dimension of data set\n\n");
printf("File name: ");
scanf("%s", file_name);
printf("Data file name: %s \n", file_name);

fp1=fopen(file_name, "r");
if(fp1 == NULL) {
	printf("\nError - Can't open the data file\n");
	printf("Hit any key to exit\n");
	c1=getchar();
	c1=getchar();
	exit(1);
	}


printf("Number of data points(vectors): ");
scanf("%d", &n);
printf("Vector size(dimension): ");
scanf("%d", &m);
printf("Number of clusters: ");
scanf("%d", &c);
printf("Enter 1 for direct clustring ---- 2 for fuzzy clustring: ");
scanf("%d", &switch1);

if(switch1 == 2) {
	do{
		printf("m: ");
		scanf("%f", &m_prime);
		if(m_prime<=1) printf("m' should be larger than 1, try again\n");
		}while(m_prime<=1);
	}

printf("Input random number seed(any integer): ");
scanf("%d", &seed);

srand(seed);

/********* CHECK FOR VIOLATION OF CARDINAL RULES *********/

flag1=1;

if(switch1 == 1 || switch1 == 2) flag1=0;

if(n<c) flag1=1;

if(flag1 == 1) {
	printf("\nError - Violation of a Cardinal Rule\n");
	printf("Hit any key to exit\n");
	c1=getchar();
	c1=getchar();
	exit(1);
	}

/********** SCANNING THE DATA ********/

for(i=0; i<n; ++i) {
	for(j=0; j<m; ++j) fscanf(fp1, "%f ", &x[i][j]);
	fscanf(fp1, "\n");
	}

fclose(fp1);

printf("Printing Data:\n");
for(i=0; i<n; ++i) {
	for(j=0; j<m; ++j) printf("%f\t", x[i][j]);
	printf("\n");
	}

/********* CREATING INITIAL U MATRIX *****/

for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) u[i][j]=0.0;
	}

for(j=0; j<c; ++j) u[j][j]=1.0;
	for(j=c; j<n; ++j) {
		select=rand()%c;
		u[select][j]=1.0;
		}


printf("Initial U matrix\n");
for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) printf("%8.4f", u[i][j]);
	printf("\n");
	}


/************ GENERATING V MATRIX *************/

do{

for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) old_u[i][j]=u[i][j];
	}
for(i=0; i<c; ++i) {
	for(j=0; j<m; ++j) {
		sum1=0;
		sum2=0;
		for(l=0; l<n; ++l) {
			if(switch1 == 2) u[i][l]=pow((u[i][l]),m_prime);
			sum1=sum1+u[i][l]*x[l][j];
			sum2=sum2+u[i][l];
			}
      if (sum2==0.0) break;
		v[i][j]=sum1/sum2;

		}
	}

/*************** COMPUTING D MATRIX ***********/

for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) {
		sum1=0;
		for(l=0; l<m; ++l) {
			sum1=sum1+pow((x[j][l]-v[i][l]),2);
			}

		d[i][j]=sqrt(sum1);
		}
	}


/********** UPDATING U MATRIX ********************/



switch(switch1) {

case 1:  /* direct clustring */

	for(j=0; j<n; ++j) min[j]=d[0][j];
	for(j=0; j<n; ++j) {

		for(i=1; i<c; ++i) {
			if(d[i][j]<min[j]) min[j]=d[i][j];
			}
		}

	for(j=0; j<n; ++j) {
		for(i=0; i<c; ++i) {
			if(d[i][j] == min[j]) u[i][j]=1;
			else u[i][j]=0;
			}
		}
	break;

case 2:
	for(i=0; i<c; ++i) {
		for(j=0; j<n; ++j) {
			sum1=0.0;
			for(l=0; l<c; ++l) {
				if(d[l][j]==0.0) d[l][j]=1e-5;
				sum1=sum1+pow((d[i][j]/d[l][j]),(2.0/(m_prime-1.0)));
				}
			u[i][j]=1.0/sum1;
			}
		}
	break;

} /* end of switch */

printf("\nUpdated U matrix\n");
for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) printf("%7.4f", u[i][j]);
	printf("\n");
	}

/********** CALCULATE IMPROVEMENT ***************/

max=u[0][0]-old_u[0][0];
if(max<0) max=-1.0*max;

for(i=0; i<c; ++i) {
	for(j=0; j<n; ++j) {
		delta=u[i][j]-old_u[i][j];
		if(delta<0) delta=-1.0*delta;
		if(delta>max) max=delta;
		}
	}

printf("Max Improvement: %f\n", max);
printf("\nHit return to continue\n");

c1=getchar();


}while(max>STOP);

printf("\n\tCluster coordinates:\n");
printf("Cluster no:\n");
for(i=0; i<c; ++i) {
	printf("%d\t\t", i);
	for(j=0; j<m; ++j) printf("%f\t", v[i][j]);
	printf("\n");
	}


printf("\nCLUSTER\t\tDATA POINTS\n");
printf("========\t=====================\n");
for(i=0; i<c; ++i){
	printf("%d\t\t", i+1);
	for(j=0; j<n; ++j) {
		if(u[i][j]>0.5) printf("%3d", j+1);
		}
	printf("\n");
	}

c1=getchar();

} /* end of main loop */