
#include <stdio.h>
#include <math.h>
#include "clustering.h"

Clustering::Clustering() : nrows(0), ncols(0), n(10) {

}

void Clustering::Allocate(int r, int c) {
	int i;
	nrows = r;
	ncols = c;
	data = new float *[nrows];
	clusternum = new int[nrows];
	for (i=0;i<nrows;i++) {
		data[i] = new float[ncols];
	}
}

void Clustering::DeAllocate() {
	int i;
	if (nrows > 0) {
		for (i=0;i<nrows;i++) {
			delete [] data[i];
		}
		delete [] data;
		delete [] clusternum;
		nrows = 0;
	}
}

float Clustering::FindDist(float *a, float *b) {

	int i;
	float r;

	r = 0.0;
	for (i=0; i<ncols; i++) {
		r = r + (a[i]-b[i])*(a[i]-b[i]);
	}
	r = sqrt(r);

	return r;
}



void Clustering::Cluster() {

	int numleft;
	int lastpoint;
	float nearestdist;
	int nearest;
	int currentcluster;
	int i;
	float dist;
	//FILE *fptr;

	numclusters = 0;
	deficit = n;
	lastpoint = 0;
	clusternum[0] = currentcluster = 0;
	numleft = nrows;

	nearest = 0;
	//fptr = fopen("debug.txt","wt");
	while (numleft > 0) {
		// find nearest neighbor of last point
		nearestdist = -1.0;
		for (i=0;i<nrows;i++) {
			if (clusternum[i] < 0) {
				dist = FindDist(data[lastpoint],data[i]);	
				if ((dist < nearestdist) || (nearestdist < -0.1)) {
					nearestdist = dist;
					nearest = i;
				}
			}
		}
		if (deficit==0) {
			// if nearest neighbor within threshold, then add to cluster
			if (nearestdist < lastndist * 5.0) {
				lastndist = (lastndist * (float)(n-1)/(float)n) +
							(nearestdist/(float)n);
			} else { // else start new cluster
				currentcluster++;
				deficit = n;
			}
		} else {
			lastndist = ((float)(n-deficit)*lastndist + nearestdist)/(float)(n-deficit+1);
			deficit--;
		}
		clusternum[nearest] = currentcluster;
		lastpoint = nearest;
		//fprintf(fptr,"nearest %d cluster %d\n", nearest, currentcluster);
		numleft--;
	}

	printf("num clusters %d\n", currentcluster + 1);
	numclusters = currentcluster + 1;
	//fclose(fptr);

}

// cptr has been allocated
// put in cptr the coords of the centroid of cluster cnum
void Clustering::FindCentroid(int cnum, float *cptr) {
	int i,j,n;
	n = 0;
	if (cnum < numclusters) {
		for (i=0;i<nrows;i++) {
			if (clusternum[i] == cnum) {
				for (j=0;j<ncols;j++) {
					cptr[j] = (data[i][j] + (float)n * cptr[j])/(float)(n+1);
				}
			}
		}
	}
}

