---
layout: post
title: K-Means
categories: [Hacking]
tags: [Software, _Phase4]
---

<figure class="align-center">
<canvas width="700" height="350" id="km_sample1" class="align-center">
</canvas>
<figcaption id='stats1'>Live demo: Click to Pause/Resume</figcaption>
</figure>

<script>
	var kmSample1 = {
		pts: null,
		m: null,
		canvas: null,
		ctx: null,
		thisSample: null,
		iter: 0,
		paused: false,
		drawStyle: 0,
		//
		sat_random_color: function() {
			// returns some strong color. We consider RG&B to be evenly-tempered, no perceptual hijinks
			var c = [Math.random(), Math.random(), Math.random()];
			var v = Math.max.apply(Math,c);
			c = c.map(function(x) { return Math.min(0.9, x/v);});
			var v = 1.0 - Math.min.apply(Math,c);
			c = c.map(function(x) { return Math.max(0.0, (1.0-x)/v);});
			c = c.map(function(x) { return Math.min(255, Math.floor(255*x));});
			return c;
		},
		desat_color: function(color) { // color in in 0-255 space
			var v = Math.max.apply(Math,color);
			c = color.map(function(x) { return Math.floor(x+0.75*(v-x));});
			return c;
		},
		init_sample1_data: function(nPoints, mMeans) {
			// define a bunch of points, ten randomly place a few mean candidates
			var i, j;
			var inset = 0.02;
			var i2 = 1.0 - inset*2;
			// distribute points in device-normalized space, just in case we get resized
			this.pts = [];
			var xf = Math.PI/4 + Math.random()*8;
			var yf = Math.PI/4 + Math.random()*4;
			var xo = Math.PI * Math.random();
			var yo = Math.PI * Math.random();
			var lumpy = (Math.random() > 0.3);
			var dens = 1.0;
			var noise = 0.01;
			var nr = 1.0-noise;
			while (this.pts.length < nPoints) {
				var x = Math.random();
				var y = Math.random();
				if (lumpy) {
					dens = noise + nr * (0.5  + 0.5*Math.cos(xo + x*xf)) * (0.5  + 0.5*Math.cos(yo + y*yf));
				}
				if (Math.random() <= dens) {
					this.pts.push({
						x: inset+i2*x,
						y: inset+i2*y,
						m: 0			// altered by update_membership()
					});
				}
			}
			// likewise means
			this.m = [];
			for (i=0; i<mMeans; i=i+1) {
				var cv = this.sat_random_color();
				var cd = this.desat_color(cv);
				this.m.push({
					x: Math.random(),
					y: Math.random(),
					cv: cv,
					c: ('rgb('+cv.join(',')+')'),
					cd: ('rgb('+cd.join(',')+')'),
					active: true
				});
			}
		},
		distance: function(p,m) {
			// distance betwee a point and a mean.
			// "True" Euclidena distance would include a square root,
			//     but we only care about relative ranking between points, so it's unneeded here.
			return (p.x-m.x)*(p.x-m.x) + (p.y-m.y)*(p.y-m.y); // no sqrt()
		},
		draw_points: function() {
			// draw points AND means
			var ip, im, ct, g;
			var w = this.canvas.width;
			var h = this.canvas.height;
			var mSize = 2;
			var mRad = 4;
			this.ctx.fillStyle = 'white';
			this.ctx.fillRect(0,0,w,h)
			this.ctx.strokeStyle = 'black';
			this.ctx.beginPath();
			this.ctx.rect(0,0,w,h);
			this.ctx.stroke();
			for (im=0; im<this.m.length; im=im+1) {
				if (!this.m[im].active)
					continue;
				var ctr = {x: this.m[im].x*w, y: this.m[im].y*h};
				ct = 0;
				for (ip=0; ip<this.pts.length; ip=ip+1) {
					if (this.pts[ip].m == im) {
						ct += 1;
						var a = 1; // use this later
						var p = {x: (this.pts[ip].x * w),
							     y: (this.pts[ip].y * h)};
						// this.ctx.strokeStyle = ('rgba('+this.m[im].cv.join(',')+','+a+')');
						g = this.ctx.createRadialGradient(ctr.x, ctr.y, mRad, ctr.x, ctr.y, h/4);
						g.addColorStop(0, this.m[im].cd);
						g.addColorStop(1, this.m[im].c);
						this.ctx.strokeStyle = g;
						this.ctx.beginPath();
						this.ctx.moveTo(p.x-mSize, p.y-mSize);
						this.ctx.lineTo(p.x+mSize, p.y+mSize);
						this.ctx.stroke();
						this.ctx.beginPath();
						this.ctx.moveTo(p.x-mSize, p.y+mSize);
						this.ctx.lineTo(p.x+mSize, p.y-mSize);
						this.ctx.stroke();
						if (this.drawStyle == 1) {
							this.ctx.moveTo(p.x, p.y);
							this.ctx.lineTo(ctr.x, ctr.y);
							this.ctx.stroke();							
						}
					}
				}
				this.ctx.strokeStyle = this.m[im].c;
				this.ctx.fillStyle = ('rgba('+this.m[im].cv.join(',')+',0.25)');
				this.ctx.beginPath();
				this.ctx.ellipse(ctr.x, ctr.y, mRad,mRad, Math.PI / 4, 0, 2 * Math.PI);
				this.ctx.fill();
				this.ctx.stroke();
			}
		},
		update_memberships: function() {
			// update points, to see if any have switched affiliations. Return a count of
			//     how many have changed.
			var ip, im;
			var nChanged = 0;
			for (ip=0; ip<this.pts.length; ip=ip+1) {
				var dBest = 2; // some large value beyond our 1x1 space
				var mBest = null;
				for (im=0; im<this.m.length; im=im+=1) {
					if (!this.m[im].active)
						continue;
					var dm = this.distance(this.pts[ip],this.m[im]);
					if (dm<dBest) {
						dBest = dm;
						mBest = im;
					}
				}
				if (mBest != this.pts[ip].m) {
					nChanged += 1;
					this.pts[ip].m = mBest;
				}
			}
			// console.log(nChanged+' changed');
			return(nChanged);
		},
		update_centroids: function() {
			// update mean locations (ignore inactive ones)
			var ip, im, n, c;
			for (im=0; im<this.m.length; im+=1) {
				if (!this.m[im].active)
					continue;
				n = 0;
				c = {x:0, y:0};
				for (ip=0; ip<this.pts.length; ip+=1) {
					if (this.pts[ip].m == im) {
						n+=1;
						c.x += this.pts[ip].x;
						c.y += this.pts[ip].y;
					}
				}
				if (n==0) { // set is EMPTY - can be discarded
					this.m[im].actve = false;
				} else { // move to centroid of member points
					this.m[im].x = c.x/n;
					this.m[im].y = c.y/n;
				}
			}
		},
		update_all: function() {
			// our complete method -- just loop on this until you don't
			var m = this.update_memberships();
			if (m > 0) {
				this.update_centroids();
			}
			return(m);
		},
		remove_one: function() {
			// randomly remove a mean, until we reach some minimum
			var i, ct;
			for (i=0; i<this.m.length; i+=1) {
				if (this.m[i].active) {
					this.m[i].active = false;
					break;
				}
			}
			for (i=0, ct=0; i<this.m.length; i+=1) {
				if (this.m[i].active) {
					ct += 1;
				}
			}
			return (ct > 2);
		},
		looper: function(timestamp) {
			// called by requestAnimationFrame() forever
			if (this.paused) {
				window.requestAnimationFrame(this.looper.bind(this));
				return;
			}
			var ch = this.update_all();
			this.draw_points();
			if (ch > 0) {
				window.requestAnimationFrame(this.looper.bind(this));
			} else if (this.remove_one()) {
				window.setTimeout(this.looper.bind(this),200);
			} else {
				window.setTimeout(this.startup.bind(this),100);
			}
			this.iter += 1;
		},
		startup: function() {
			// also called whenever we re-start
			this.init_sample1_data(2000,12);
			this.iter = 0;
			window.requestAnimationFrame(this.looper.bind(this));
		},
		toggle_pause: function() {
			// user can click to stop/start the animation
			this.paused = ! this.paused;
		},
		main: function(canvID) {
			this.canvas = document.getElementById(canvID);
			var p = this.canvas.parentElement;
			if (p.offsetWidth < (this.canvas.width-4)) {
				this.canvas.width = p.offsetWidth - 4;
			}
			this.ctx = this.canvas.getContext('2d');
			this.startup();
			this.canvas.onclick = this.toggle_pause.bind(this);
		}
	}
	window.addEventListener('load', function s2() {kmSample1.main("km_sample1"); });
</script>

_This is the first post in a brief series on algorithms (simple methods) that I like and that I think more people should know about: either because of their simplicity, their novelty, their usefulness, or all three. Each comes with a live demo. To begin, let's introduce: **k-Means Clustering** ._
{: .notice--info}

<!--more-->

### The Idea

This sample shows _k-Means_ in action. k-Means used to also be known as _Lloyd's method._ It's so simple that I'm surprised it's not used more often.

The basic idea is: for some pile of info, whether its sales records for varying shoe sizes, the location of weeds in your yard, whatever, find a set of "best middles" &mdash; or "means" &mdash; for any given set of information. In the sample, we have a bunch of random, possibly-lumpy collections of little x's on the canvas. The method tries to find some smaller number of categories (the "means," marked by circles) that could best-represent that info.

The "k" in "k-Means" is the number of categories we want to create for dividing-up the samples.

The k-Means method doesn't directly solve for "best." It just keeps trying and adjusting, starting from completely random guesses as to where to put the means:

* Look at every point, and mark it "owned" by the closest "mean."
* Move each mean to the middle of its "owned" points.
* _Just keep doing those two steps,_ until either nothing changes each time or you run out of patience. If a proposed mean has _no_ points, just discard that mean. You didn't need it.

---

### The Example Code

In this sample, we let the means move around until they're stable, then randomly take one away, which forces all the rest to re-examine their situatons.

You can pause/start the animation by clicking on the canvas.

You can see the [source here on GitHub.](https://raw.githubusercontent.com/joker-b/botzo/master/_posts/2018-12-05-K-Means.md) The entire example is contained in a single Javascript module that's imaginatively named `kmSample1`

You can decide for yourself what these results look like: Popping soap bubbles? Growing algae mats in time-lapse? Voting districts? That dream about insects eating your car? I suggest turning on some [Cinematic Orchestra](https://www.cinematicorchestra.com/) for your headphones and stare at this page for a long while to decide.

In future examples, we'll try using k-Means (and other tricks) in different ways. Some obvious questions might be:

* Are there better or [different ways to decide which point belongs to which mean?]({{ site.baseurl}}{% post_url 2018-12-07-Special-K %}) What do we mean by "closeness"?
* What might this have to do with [photography]({{ site.baseurl}}{% post_url 2018-12-10-K-Pix %})?