#!/usr/bin/Rscript
a <- c(2, 10, 1, 5, 8, 4, 5, 8, 7, 5, 6, 4, 2, 2, 4, 9)
dim(a) <- c(2,8)
a <- t(a)
plot(a, col=1:8, pch=19)
text(a[,1]+0.15, a[,2]+0.15, 1:8)
print(dist(a))
print(dist(a) ^ 2)

medoids <- c(1,4,7)
step <- 0
dirty <- TRUE
while(dirty)
{
	d <- as.matrix(dist(a),"manhattan")[medoids,]
	mapping <- apply(d, 2, which.min)
	dirty <- FALSE
	cat("step", step, ":\n")
	cat("cluster mappings:\n")
	print(mapping)
	plot(a, col=mapping, pch=1)
	text(a[,1]+0.15, a[,2]+0.15, 1:8)
	cat("medoids:\n")
	print(medoids)
	points(a[medoids,], pch=19, col=1:3)
	Sys.sleep(1)
	for(m in medoids)
	{
		points_in_cluster <- which(medoids[mapping]==m)
		cur_d <- sum(as.matrix(dist(a),"manhattan")[points_in_cluster,m])
		for(o in points_in_cluster)
		{
			new_d <- sum(as.matrix(dist(a),"manhattan")[points_in_cluster,o])
			if(new_d < cur_d)
			{
				medoids[which(medoids==m)] <- o
				dirty <- TRUE
			}
		}
	}
	step <- step + 1
}