wkmean.c

#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>
#include <math.h>

/* a wavelet-based, k-means 2D clustering algorithm */

#include "evaluate.h"

#ifdef DEBUG_CLUSTER
// for debugging and performance analysis
static void dump_clusters(evaluate * e)
{
    int i;
    
    printf("X clusters\n");
    for (i = 0; i < e->current_num_clusters; i++) {
	printf("%i: %lf +/- %lf (%lu)\n", i, e->clusters[0][i].mean,
	       sqrt(e->clusters[0][i].variance),
	       e->clusters[0][i].count);
    }
    
    printf("F clusters\n");
    for (i = 0; i < e->current_num_clusters; i++) {
	printf("%i: %lf +/- %lf (%lu)\n", i, e->clusters[1][i].mean,
	       sqrt(e->clusters[0][i].variance),
	       e->clusters[1][i].count);
    }
}
#else
#define dump_clusters(e)
#endif

// max number of refinement passes to take
#define MAX_REFINE_PASSES 10

void evaluate_cluster_kmean_wavelet(evaluate * e)
{
    uint32_t i;
    uint32_t j;
    uint32_t k;
    uint32_t l;
    uint32_t m;
    uint32_t n;
    
    int a = 0;
    int ai;
    int b;
    
    // decompose dataset
    // array is iteratively replaced by tuples of the average and deviance 
    // from avg for two given, likely disjoint values in the array.
    while (e->current_num_datapoints >> (a + 1) >= e->current_num_clusters) 
    {
	// proper array positioning
	l = 1 << a;
	j = l << 1;
	k = e->current_num_datapoints - l;
	a++;
	
	for (i = 0; i < k; i += j) {
	    // calculate midpoints and average deviance
	    float xa = 0.5 * (e->datapoints[i].x + e->datapoints[i+l].x);
	    float xd = 0.5 * (e->datapoints[i].x - e->datapoints[i+l].x);
	    float fa = 0.5 * (e->datapoints[i].f + e->datapoints[i+l].f);
	    float fd = 0.5 * (e->datapoints[i].f - e->datapoints[i+l].f);
	    // replace existing values with avg and dev
	    e->datapoints[i].x = xa;
	    e->datapoints[i+l].x = xd;
	    e->datapoints[i].f = fa;
	    e->datapoints[i+l].f = fd;
	}
    }
    
    ai = a;
    
    // reset cluster values
    for (m = 0; m < e->current_num_datapoints >> a; m++) 
    {
	e->datapoints[m << a].cluster_x =
	    e->datapoints[m << a].cluster_f = 0xff;
    }
    
    // compose and refine clusters
    for (;;) {
	int refine = 1;
	int passes;
	// refine current (or define initial) clusters
	for (passes = 0; refine && passes < MAX_REFINE_PASSES; passes++)
	{
	    for (b = 0; b < e->current_num_clusters; b++) {
		// x-axis
		// determine current mean and variance for this cluster
		if (e->clusters[0][b].count) {
		    e->clusters[0][b].mean = 
			e->clusters[0][b].sum / (float) e->clusters[0][b].count;
		    if (e->clusters[0][b].count > 1)
			e->clusters[0][b].variance =
			    (e->clusters[0][b].ssum -
			     e->clusters[0][b].sum * e->clusters[0][b].sum /
			     (double) e->clusters[0][b].count) /
			    ((double) e->clusters[0][b].count - 1.0);
		    else
			e->clusters[0][b].variance = 0.0;
		} else {
		    // no datapoints for this cluster, so assign a random mean.
		    n = random() % (e->current_num_datapoints >> ai);
		    e->clusters[0][b].mean = e->datapoints[n << ai].x;
		    e->clusters[0][b].variance = 0.0;
		}
		// f-axis
		// ditto: find mean and variance
		if (e->clusters[1][b].count) {
		    e->clusters[1][b].mean = 
			e->clusters[1][b].sum / (float) e->clusters[1][b].count;
		    if (e->clusters[1][b].count > 1)
			e->clusters[1][b].variance =
			    (e->clusters[1][b].ssum -
			     e->clusters[1][b].sum * e->clusters[1][b].sum /
			     (double) e->clusters[1][b].count) /
			    ((double) e->clusters[1][b].count - 1.0);
		    else
			e->clusters[1][b].variance = 0.0;
		} else {
		    // or just make something up for now
		    n = random() % (e->current_num_datapoints >> ai);
		    e->clusters[1][b].mean = e->datapoints[n << ai].f;
		    e->clusters[1][b].variance = 0.0;
		}
		e->clusters[0][b].sum = e->clusters[1][b].sum = 0.0;
		e->clusters[0][b].ssum = e->clusters[1][b].ssum = 0.0;
		e->clusters[0][b].count = e->clusters[1][b].count = 0;
	    }
	    
	    refine = 0;

	    // for each data point...
	    for (m = 0; m < e->current_num_datapoints >> a; m++) 
	    {
		float x = e->datapoints[m << a].x;
		float f = e->datapoints[m << a].f;
		float t;
		
		int c_x = 0;
		int c_f = 0;

		// find deviance from the initial cluster's mean
		float d_x = fabs(e->clusters[0][0].mean - x);
		float d_f = fabs(e->clusters[1][0].mean - f);

		// find which cluster it is closest to...
		for (b = 1; b < e->current_num_clusters; b++) {
		    t = fabs(e->clusters[0][b].mean - x);
		    if (t < d_x) { d_x = t; c_x = b; }
		    t = fabs(e->clusters[1][b].mean - f);
		    if (t < d_f) { d_f = t; c_f = b; }
		}
		
		// alter the closest cluster to include this data point
		
		// on the x-axis
		e->clusters[0][c_x].sum += x;
		e->clusters[0][c_x].ssum += x * x;
		e->clusters[0][c_x].count++;
		if (e->datapoints[m << a].cluster_x != c_x) refine = 1;
		e->datapoints[m << a].cluster_x = c_x;
		
		// and on the f-axis
		e->clusters[1][c_f].sum += f;
		e->clusters[1][c_f].ssum += f * f;
		e->clusters[1][c_f].count++;
		if (e->datapoints[m << a].cluster_f != c_f) refine = 1;
		e->datapoints[m << a].cluster_f = c_f;
	    }
	}
	
	// we can't recompose the dataset to a higher level than we
	// started with initially...
	if (a == 0 || a <= ai >> 2) break;
	
	// array positioning again
	a--;
	l = 1 << a;
	j = l << 1;
	k = e->current_num_datapoints - l;
	
	// and we recompose the datapoints
	for (i = 0; i < k; i += j) {
	    // the x
	    float x1 = e->datapoints[i].x + e->datapoints[i+l].x;
	    float x2 = e->datapoints[i].x - e->datapoints[i+l].x;
	    // the f
	    float f1 = e->datapoints[i].f + e->datapoints[i+l].f;
	    float f2 = e->datapoints[i].f - e->datapoints[i+l].f;
	    
	    // and make it so:
	    e->datapoints[i].x = x1;
	    e->datapoints[i+l].x = x2;
	    e->datapoints[i+l].cluster_x = e->datapoints[i].cluster_x;
	    e->datapoints[i].f = f1;
	    e->datapoints[i+l].f = f2;
	    e->datapoints[i+l].cluster_f = e->datapoints[i].cluster_f;
	}
    }

    // normalize
    for (i = 0; i < e->current_num_clusters; i++) {
	e->clusters[0][i].count <<= a;
	e->clusters[1][i].count <<= a;
    }
    
    // reset the return value
    memset(e->histogram, 0, sizeof(uint32_t) * 
	   e->current_num_clusters * e->current_num_clusters);

    // and assign weights to the resulting clusters
    for (i = 0; i < e->current_num_datapoints >> a; i++) {
	e->histogram[e->datapoints[i << a].cluster_x * e->current_num_clusters +
		     e->datapoints[i << a].cluster_f] += 1 << a;
    }
}