scalablekmeans.cpp

Scalable k-means software and test datasets This package (a Unix tar file, gzipped) contains the sou

//		printf("Tried %d %d\n", mi, mj);
		allcluster[mi]->addCluster(allcluster[mj]);
		allcluster[mi]->computeStdDev();
		if(allcluster[mi]->isTight(beta))
		{
			printf("Joined %d %d %d\n", mi, mj,
					allcluster[mi]->getNumPoints());
			for(i = 0; i < totalclusters; i++)
			{
				if(distances[mi*totalclusters + i] >= 0)
					distances[mi*totalclusters + i] = -1;
				if(distances[i*totalclusters + mi] >= 0)
					distances[i*totalclusters + mi] = -1;
				distances[i*totalclusters + mj] = -2;
				distances[mj*totalclusters + i] = -2;
//				printf("distances[%d,%d]=-2\n", i, mj);
//				printf("distances[%d,%d]=-1\n", mi, i);
			}
			distances[mi*totalclusters + mj] = -2;
//			printf("distances[%d,%d]=-2\n", mi, mj);
//			printf("A\n");
			if(mj >= oldnumclusters)
			{	
/*
Point in "retained set" belongs to a _new_ tight secondary cluster
Remove point if that cluster is joined with cluster in compressed
set.
Remove point if that cluster can be added to compressed set.
*/

//				delete allcluster[mj];
				((MainCluster*)allcluster[mj])->id = 3;
				if(mi >= oldnumclusters)
				{	for(p = buffer.getRetainedSet(); p; p = p->getNext())
						if(p->getAssociation() == allcluster[mj])
							p->associate(allcluster[mi]);
				}

				allcluster[mj] = 0;
//				printf("B\n");
			}
			else
			{	buffer.removeFromCompressedSet(allcluster[mj]);
				delete allcluster[mj];
				allcluster[mj] = 0;
//				printf("C\n");
			}
		}
		else
		{	allcluster[mi]->subCluster(allcluster[mj]);
			allcluster[mi]->computeStdDev();
			distances[mi*totalclusters + mj] = -2;
		}
	}

/*
0 Tight, but doesn't fit in buffer.
   Remove cluster
1 Put in compressed set
   Remove points
2 not tight from beginning
   Remove cluster
3 Merged with cluster in compressed set, or points relabeled
   Remove cluster, remove points
*/

	fflush(stdout);

	for(i = 0; i < numclusters; i++)
	{
		if(!newcluster[i]->id)
		{	if(newcluster[i]->getMemUsed() <
				newcluster[i]->getNumPoints()*pointmem)
			{	newcluster[i]->id = 1;
//				printf("Marked cluster.\n");
			}
			else
			{	printf("Secondary cluster doesn't fit.\n");
			}
		}
	}

	for(p = buffer.getRetainedSet(); p; p = np)
	{
		np = p->getNext();
		j = ((MainCluster*)p->getAssociation())->id;
		if(j == 1 || j == 3)
		{	buffer.removeFromRetainedSet(p);
			delete p;
//			printf("Compressed.\n");
		}
	}

	for(i = 0; i < numclusters; i++)
	{
		if(newcluster[i] && newcluster[i]->id == 1)
		{	buffer.addSecondaryCluster(newcluster[i]);
			newcluster[i] = 0;
			printf("Added cluster; %d.\n", buffer.getCompressedSetSize());
		}
	}

/*	for(i = oldnumclusters; i < totalclusters; i++)
	{
		if(allcluster[i] && ((MainCluster*)allcluster[i])->id != 1)
			delete allcluster[i];
	}
*/
	for(i = 0; i < numclusters; i++)
	{
		if(newcluster[i] && newcluster[i]->id != 1)
			delete newcluster[i];
	}

	delete []distances;
	delete []allcluster;
//	delete joincluster;


	// Try to join new tight clusters with each other and clusters
	// in compressed set.

//	for(i = 0; i < numclusters; i++)
//		if(newcluster[i])
//			delete newcluster[i];

	delete []newcluster;

	fflush(stdout);
}

int ScalableKMeans::updateModel(void)
{
	Subcluster *c, *c2;
	Singleton *p, *mp;
	int dim = buffer.getDimensions(), k;
	float moved, div = ((float)dim)*dim*numClusters;
	float d;
	int l = 0;
	int q = 0;

	do
	{
		// Take each point in retained set, new data points, and compressed
		// sets and assign to closest main cluster from prev. iteration.

		// Update association

		for(p = buffer.getRetainedSet(); p; p = p->getNext())
			p->associate(closestMainCluster(p));

		for(c = buffer.getCompressedSet(); c; c = c->getNext())
			c->associate(closestMainCluster(c));

		// Recompute main cluster means.
		// Use retained set, new data points, compressed sets and discard
		// sets.

		for(c = mainClusters; c; c = c->getNext())
			c->clear();

		for(c = buffer.getDiscardSet(); c; c = c->getNext())
			c->getAssociation()->addCluster(c);

		for(p = buffer.getRetainedSet(); p; p = p->getNext())
			p->getAssociation()->addPoint(p);

		for(c = buffer.getCompressedSet(); c; c = c->getNext())
			c->getAssociation()->addCluster(c);

		moved = 0;
		k = 0;
		for(c = mainClusters; c; c = c->getNext())
		{
			if(!c->getNumPoints())
			{	fprintf(stderr, "Warning: Cluster has no points.\n");
				printf("Warning: Cluster has no points.\n");
			}
			moved += c->updateMean();
//			printf("%d ", k++);
			c->print();
		}
		moved /= div;

//		printf("%g\n", moved);

		numMainIterations++;
		q++;
	} while(moved > convergenceThreshold && ++l < 100);

//	printf("\n");

	numUpdates++;

	c2 = oldClusters;
	moved = 0;
	for(c = mainClusters; c; c = c->getNext())
	{
		if(!c->getNumPoints())
		{
			d = 0;
			for(p = buffer.getRetainedSet(); p; p = p->getNext())
				if(p->clusterDistance > d)
				{	d = p->clusterDistance;
					mp = p;
				}
			c->addPoint(mp);
			c->updateMean();
			mp->associate(c);
			mp->clusterDistance = 0;
		}
		c2->clear();
		c2->addCluster(c);
		moved += c2->updateMean();
		c2 = c2->getNext();
	}

	if(logFile)
		fprintf(logFile, "%ld %g %d\n", buffer.getTotalPoints(), moved/div, q);

//	fprintf(stderr, "Moved: %f\n", moved/div);

	// Return 1 if converged, else 0

#if EARLYSTOPPING == 0
	return 0;
#else
	return moved/div < 1e-8;
#endif
}


Subcluster* ScalableKMeans::closestMainCluster(Subcluster *c)
{
	float d = 1e20, t;
	Subcluster *m, *r;

	for(m = mainClusters; m; m = m->getNext())
		if((t = m->distanceSquared(c)) < d)
		{	r = m;
			d = t;
		}

	return r;
}


Subcluster* ScalableKMeans::closestMainCluster(Singleton *p)
{
	float d = 1e20, t;
	Subcluster *m, *r;

	for(m = mainClusters; m; m = m->getNext())
		if((t = m->distanceSquared(p)) < d)
		{	r = m;
			d = t;
		}

	p->clusterDistance = d;

	return r;
}


ScalableKMeans::ScalableKMeans(long bufferSize, Database *db, int k,
		int dim) :
		buffer(bufferSize, dim), database(db), numClusters(k)
{
	setup();
	beta = defaultBeta;
	convergenceThreshold = defaultConvergenceThreshold;
	secondaryConvergenceThreshold = defaultSecondaryCT;
	discardFraction = defaultDiscardFraction;
	confidenceStdDev = defaultConfidenceStdDev;
	primary1 = 1;
	primary2 = 1;
	secondary = 1;
	naive = 0;
	logFile = 0;
	numNewSecondaryClusters = k*4;
//	numNewSecondaryClusters = k*8;
}


ScalableKMeans::~ScalableKMeans()
{
	Subcluster *p, *np;

	for(p = mainClusters; p; p = np)
	{	np = p->getNext();
		delete p;
	}

	for(p = oldClusters; p; p = np)
	{	np = p->getNext();
		delete p;
	}

	if(logFile)
		fclose(logFile);
}


void ScalableKMeans::exportCompressedSet(void)
{
	Subcluster *c;

	FILE *f = fopen("compressed.m", "w");

	fprintf(f, "c=[");
	for(c = buffer.getCompressedSet(); c; c = c->getNext())
	{
//		fprintf(stderr, "C\n");
		c->computeStdDev();
		c->print(f);
	}
	fprintf(f, "];\n");

	fprintf(f, "ci=[");
	for(c = buffer.getCompressedSet(); c; c = c->getNext())
	{
		fprintf(f, "%d ", c->getNumPoints());
		c->printConfidence(f);
	}
	fprintf(f, "];\n");

	fprintf(f, "for i=1:size(c,1)\n");
	fprintf(f, "sd=sqrt(ci(i,1));\n");
	fprintf(f, "rectangle('Curvature',[1 1],'Position', [c(i,1)-sd*ci(i,2) "\
			"c(i,2)-sd*ci(i,3) 2*sd*ci(i,2) 2*sd*ci(i,3)]);\n");
	fprintf(f, "end;\n");

	fclose(f);
}


void ScalableKMeans::naiveCompression(void)
{
	Singleton *p, *np;

	for(p = buffer.getRetainedSet(); p; p = np)
	{
		np = p->getNext();
		p->getAssociation()->getAssociation()->addPoint(p);
		buffer.removeFromRetainedSet(p);
		delete p;
	}

}


void ScalableKMeans::cluster(void)
{
	Singleton *p;
	Singleton *nq;
	Subcluster *c;

	numMainIterations = numSecondaryIterations = numUpdates = 0;

	initModel();

	database->reset();
	fillBuffer();

	int k = 0;
	do
	{
//		for(Singleton *q = buffer.getRetainedSet(); q; q = q->getNext())
//			q->print(stdout);

#if OUTPUTFILES == 1
		exportCompressedSet();
#endif
		if(updateModel())
			break;
		if(naive)
			naiveCompression();
		else
		{	primaryCompression();
			if(secondary)
				secondaryCompression();
		}

/*		for(Singleton *q = buffer.getRetainedSet(); q; q = nq)
		{	nq = q->getNext();
			q->print(stdout);
			buffer.removeFromRetainedSet(q);
			delete q;
		}
*/
		for(p = buffer.getRetainedSet(); p; p = p->getNext())
			p->retain();

	} while(fillBuffer());
//	} while(fillBuffer() && ++k < 100);

//	fprintf(stderr, "Total points: %d\n", buffer.getTotalPoints());
//	printf("Total points: %d\n", buffer.getTotalPoints());

	k = 0;
	for(c = buffer.getDiscardSet(); c; c = c->getNext())
		k += c->getNumPoints();
	printf("Points in discard set: %d\n", k);
	k = 0;
	for(c = buffer.getCompressedSet(); c; c = c->getNext())
		k += c->getNumPoints();
	printf("Points in compressed set: %d\n", k);
	k = 0;
	for(p = buffer.getRetainedSet(); p; p = p->getNext())
		k++;
	printf("Points in retained set: %d\n", k);
}


double ScalableKMeans::logLikelihood(void)
{
	Subcluster *c;
	float *x;
	Singleton p(buffer.getDimensions());
	double ll, d, ll2;

	for(c = mainClusters; c; c = c->getNext())
	{
		c->computeStdDev();
	}

	database->reset();
	ll = 0;
	while(x = database->getPoint())
	{
		p.set(x);
		c = closestMainCluster(&p);

		d = c->distanceSquared(&p);
		ll2 = ll;
		ll += d;
		if(ll < ll2)
		{	fprintf(stderr, "Overflow\n");
			printf("Overflow\n");
		}
	}
	return ll;
}


double ScalableKMeans::distanceToTrueClusters(SyntheticData *db)
{
	int i, k, l;
	int permutation[numClusters];
	int used[numClusters];
	double md = 1e30, d;

	for(l = 0; l < numClusters; l++)
	{	used[l] = 0;
		permutation[l] = -1;
	}

	l = 0;
	k = 0;
	while(l >= 0)
	{
		if(l == numClusters)
		{
//			for(i = 0; i < numClusters; i++)
//				printf("%d ", permutation[i]);
//			printf("\n");
			d = permutationDistance(permutation, db);
			if(d < md)
				md = d;
			l--;
			continue;
		}
		if(permutation[l] >= 0)
			used[permutation[l]] = 0;
		while(++permutation[l] < numClusters && used[permutation[l]]);
		if(permutation[l] >= numClusters)
			l--;
		else
		{	used[permutation[l]] = 1;
			l++;
			if(l < numClusters)
				permutation[l] = -1;
		}
	}

	return md;
}


double ScalableKMeans::permutationDistance(int *permut, SyntheticData *db)
{
	Subcluster *c;
	double distance = 0;
	Singleton p(buffer.getDimensions());
	int i = 0;

	for(c = mainClusters; c; c = c->getNext())
	{
		p.set(db->getMeans(permut[i++]));
		distance += sqrt(c->distanceSquared(&p));
	}
	return distance;
}


void ScalableKMeans::logUpdates(char *filename)
{
	if(logFile)
		fclose(logFile);
	if(!(logFile = fopen(filename, "a")))
		fprintf(stderr, "Error: Can't open logfile '%s'.\n", filename);
	else
		fprintf(logFile, "*****\n");
}


/* End of file scalablekmeans.cpp */