#include "mex.h"
#include <math.h>

/* 
 * Kmeans.c
 * Kmeans clusters given data into defined number of clusters.
 * Authors Anu Niemi, Annemari Auvinen
 */

/*
 * Initializes the middle points by taking data randomly fom
 * data array. Same data can be selected only once.
 * k	Number of klusters
 * n	Number of rows in original data matrice 
 * m 	Number of columns in original data matrice
 * w	Array consisting of the  middle points
 * data	Array of the data
 */
void initialize(int k, int n, int m, double w[], double data[])
{
  int randNumber, found, j, i, l;
  int *used;
  
  used = mxCalloc(k, sizeof(int)); //Includes used random values
  for(i=0; i<k; i++) {
    used[i] = -1;
  }
  
  for (j=0; j<k; j++) {
    found = 0;
    do {
      randNumber = rand() % n;
      found = 1;
      for (i=0; i<=j; i++) {
        if (used[i] == randNumber) 
          found = 0;
      }
      if (found == 1) {
        used[j] = randNumber;
        for ( l=0; l<m; l++) {
          w[j+l*k] = data[randNumber+l*n];
        }
      }
    } while(found == 0);
  }
  mxFree(used);
  
}


/**
 * Calculates distances of the data items and middle points and defines that data item belongs
 * to the cluster which middle point's distance from it is shortest.  
 * w	Array of the middle points
 * data	Array of the data
 * k	Number of clusters
 * n	Number of rows in original data matrice
 * m	Number of columns in original data matrice
 * s	Array consists information which data items belongs to which clusters	
 */
void calculateDistances(double w[], double data[], int k, int n, int m, int s[])
{
  int i, j, l, p, index;
  double minDist, sum;
  double calcDist = 0.0;
  for(i=0; i<n; i++) {
    for(j=0; j<k; j++) { 
      s[i*k+j] = -1;
    }
  }

  for (j=0;j<n; j++ ) {
    for (p=0;p<k; p++) {
       calcDist = 0.0;
       sum = 0.0;
       for (l=0;l<m; l++) {
        sum = sum + pow((data[j+l*n] - w[p+l*k]), 2);
      }
      calcDist = sqrt(sum);
      if (p == 0) {
        minDist = calcDist;
        index = p;
      }
      else if (abs(calcDist) < abs(minDist)) {
        minDist = calcDist;
        index = p;
      }
    }
    for(i=0;i<n; i++) {
      if (s[index+i*k] == -1) {
        s[index+i*k] = j;
        break;
     }
    }
  }
}

/*
 * Calculates middle points of the clusters
 * data	Array of data
 * s	Array, which defines the data items belonging to each cluster
 * k	Number of clusters
 * nw	Array on calculated middle points
 * n	Number of rows in original data matrice
 * m	Number of columns in original data matrice 
 */
void calculateMidPoints(double data[], int s[], int k, double nw[], int n, int m) 
{
  int i, j, l, index;
  double sum;

  for (i=0; i<k; i++) {
    for (j=0; j<m; j++) {
      sum = 0.0;
      for ( l=0; l<n; l++) {
        index = s[i+l*k]; 
        if ( index == -1 && l != 0 ) {                                    
          nw[i+j*k] = sum / l;
          break;
        }
        sum = sum + data[index+j*n];
      }
    }
  }
}

/*
 * Checks if the new distance of all old and new middle points are minor than epsilon. 
 * Returns 1 it is, otherwise 0.
 * w	Array of old middle points
 * nw	Array of new middle points
 * k	Number of clusters 
 * m	Number of columns in the original data matrice
 */
int endCriteria(double w[], double wn[], double epsilon, int k, int m) 
{
  int i, j;  
  double distance;
  double sum;

  for(i=0; i<k; i++)  {
    distance = 0.0;
    sum = 0.0;
    for(j=0; j<m; j++) {  
      sum = sum + pow((w[i+j*k] - wn[i+j*k]),2);
    }
    distance = sqrt(sum);
    if(distance > epsilon) {    
      for(i=0; i<k*m;i++) {
        w[i] = wn[i];   
      }
      return 0;
    }
  }
  for(i=0; i<k*m;i++) {
    w[i] = wn[i];   
  }
  return 1;
}
 
/**
 * Clusters data into defined number of clusters.
 * k	Number of clusters
 * n	Number of rows in the original data matrice
 * m	Number of columns in the original data matrice
 * data	Array of the data
 * w	Array of the middle points
 * s	Array, which defines the data items belonging to each cluster
 */
void kmeans(int k, int n, int m, double data[], double w[], int s[])   
{
  double *nw, epsilon;
  epsilon = 0.0001;
  nw = mxCalloc(k*n,sizeof(double));
  
  initialize(k,n,m, w, data);
  do{
    calculateDistances(w, data, k, n, m,s);  
    calculateMidPoints(data, s, k, nw, n, m);
  } while(endCriteria(w, nw, epsilon, k, m) == 0);
}


/* Gateway routine of the MEX-function. 
 * 
 */
void mexFunction(
  int nlhs, mxArray *plhs[],
  int nrhs, const mxArray *prhs[] )
{
  int k, n, m, j;
  double *x, kt, *w, *s;

  /* Check for proper number of arguments */
  if (nrhs != 2)
    mexErrMsgTxt("Two inputs in the form (k, data) is required.");
  else if (!mxIsNumeric(prhs[0]) || !mxIsNumeric(prhs[1]))
    mexErrMsgTxt("One of the inputs is not numeric.");
  else if (nlhs != 2)
    mexErrMsgTxt("Two outputs in the form [middlepoints, clusters] is required.");

  /* Create a pointer to a copy of the vector data. */
  x = mxGetPr(prhs[1]);

  kt=mxGetScalar(prhs[0]);
  k = (int)kt;
	

  /* Get the dimensions of the data vector */
  m = mxGetN(prhs[1]);
  n = mxGetM(prhs[1]);
  if(k<1 || k>n)
    mexErrMsgTxt("Number of clusters should be between 0 and the rows in the data matrice");

  /* Create the output w and s  */
  plhs[0] = mxCreateDoubleMatrix(k,m,mxREAL);
  w = mxGetPr(plhs[0]);
  plhs[1] = mxCreateNumericMatrix(k,n, mxINT32_CLASS, mxREAL);
  s = mxGetPr(plhs[1]);
    
  kmeans(k, n, m, x, w, s);
 
}