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

#define NN 10000					/* max # of the input samples */
#define MAX_X 100
#define DIM 1 
#define PI 3.14159265358979
#define DEBUG
#define PREC 0.001				/* threshold */

typedef int TPoint[DIM] ;		/* Point in DIM space when we go farther */

TPoint x [ NN ];				/* samples */
int N = 0;						/* # of data */
int Nc = 0;						/* # of clusters */

double * mean;					/* arrays of means */
double * stdv;					/* arrays of standart deviations */

double * yp;					/* prob. to be in class p after E step */

double Co;						/* Constant for gaus */

double gaussian (TPoint p, int n)	/* point and # of gaus. */
{
  return ((Co*1/sqrt(stdv[n]))*exp(-1*pow((p[0]-mean[n]),2)/(2*pow(stdv[n],2))));
}

void generate_data (TPoint x [],	/* data */
					int Nd,		/* number of data */
					int mean,	/* mean of the distribution */
					int stdv ) {	/* varience */
  int i;
  for (i=0; i<Nd; i++){
	double x1=1.0*rand()/RAND_MAX,x2=1.0*rand()/RAND_MAX, z=(sqrt(-2*log(x1))*cos(2*PI*x2));
	/* z - random normal with mean 0 and var =1 */
	x[N+i][0] = (int) ceil(mean + z*stdv);		/* move and scale */
  }	
  N += Nd;
  Nc++;
}
	   
int main () {
  int i,k,was;
  
  Co	= pow(1/sqrt(2*PI),DIM);
	
  generate_data(x, 4*NN/10, 200, 10);
  generate_data(x, 3*NN/10, 130, 30);
  generate_data(x, 2*NN/10, 100, 30);
  /*x[0][0]=-1000;
	x[1][0]=-900;*/
  //  Nc++;
  mean  = (double *) calloc ( Nc, sizeof(double));
  stdv  = (double *) calloc ( Nc, sizeof(double));
  yp    = (double *) calloc ( N*Nc, sizeof(double));
  // initialize distributions
  for (i=0; i< Nc; i++) {
	mean[i] = i * MAX_X / (Nc+1);
	stdv[i]  = 10;
  }

  do {
	was = 0;
	// E-step
	for (i=0;i<N;i++) {
	  double sum=0;
	  for (k=0; k<Nc; k++)  sum += gaussian(x[i],k); 
	  for (k=0; k<Nc; k++)  {
		if (sum!=0)
		  yp [i*Nc+k] = gaussian(x[i],k)/sum;	/* P(k|x(i)) */
		else
		  yp [i*Nc+k] = 1/Nc;	/* point is very far from all of
								   them... so let's assume that it is
								   equaly probable from all classes */
	  }
	}
	
	// M-step
	for (k=0; k<Nc; k++) {
	  double sum=0, sump=0;
	  for (i=0; i<N; i++) {
		sum += yp[i*Nc+k]*x[i][0];
		sump+= yp[i*Nc+k];
	  }
	  if(fabs(1-(sum/sump)/mean[k])>PREC) was=1;
	  mean[k] = sum/sump;
	  sum = 0;
	  for (i=0; i<N; i++)
		sum += yp[i*Nc+k]*(x[i][0]-mean[k])*(x[i][0]-mean[k]);
	  if ((1-fabs(sqrt(sum/sump)/stdv[k]))>PREC) was=1;
	  stdv[k] = sqrt(sum/sump);
	  
	  if (stdv[k]<1) {
		printf("---!!!!!Tends to have varience -> 0. Artificially broaden!!!---\n");
		stdv[k]=100;
	  }
		 
	
  
	  printf("New distribution #%d mean=%5.2f stdv=%3.2f \n",
			 k,mean[k],stdv[k]);
	  
	}
	printf("\n");
  }  while (was);
  
  free(mean); free(stdv); free(yp);
  return 0;
}