LaundrySorcery

laundryclustery.d (3913B)

---

 import std.stdio;
 import std.range;
 import std.algorithm;
 import std.typecons;
 import std.math;
 import std.traits;
 
 alias Point=uint;
 
 /**
  * Returns mean and variance of the points given in the range
  */
 auto gaussian(T=double, R)(R r) if(isForwardRange!R && isNumeric!(ElementType!R)){
 	auto p=r
 		.fold!((a,b){return tuple!(T,T,size_t)(a[0]+b, a[1]+b^^2,a[2]+1);})(tuple!(T,T,size_t)(0,0,0));
 	return tuple!(T,"mean",T,"variance")(p[0]/p[2],p[1]/p[2]-(p[0]/p[2])^^2);
 }
 unittest{
 	auto res=gaussian([3,4,5]);
 	assert(res.mean==4.0);
 	assert(res.variance==2.0);
 }
 
 //struct Cluster(T) if(hasLength!T && isForwardRange!T && is(ElementType!T == Point)){
 //	T points;
 struct Cluster{
 	Point[] points;
 	double mean;
 	double variance;
 
 	void add(Point p){
 		points~=p;
 	}
 	void reset(){
 		points.length=0;
 	}
 	
 	void calculate(){
 		auto g=gaussian(points);
 		mean=g[0];
 		variance=g[1];
 	}
 	
 	float calculate_delta(){
 		auto m=mean,v=variance;
 		calculate();
 		return (m-mean)^^2+(v-variance)^^2;
 	}
 }
 
 /**
  * Returns a pointer to the cluster that is closest to point p
  */
 auto closest(Cluster[] clusters, Point p)
 in{
 	assert(clusters.length>0);
 }
 out(res){
 	assert(res[0]!=null);
 }
 do{
 	Cluster *res;
 	float mindist=float.max;
 	foreach(ref c; clusters){
 		auto dist=abs(p-c.mean);
 		if(dist<mindist){
 			res=&c;
 			mindist=dist;
 		}
 	}
 	return tuple(res,mindist);
 }
 
 /**
  * Adds the point p to the cluster whose mean is closest to it
  */
 void addToClosest(Cluster[] clusters, Point p){
 	closest(clusters,p)[0].add(p);
 }
 
 
 /**
  * The k-means++ algorithm
  */
 void kmeans(Point[] points, ref Cluster[] clusters,uint maxiter=100){
 	import std.random;
 	
 	//K-means++ initialization
 	clusters[0].mean=points.choice();
 	ulong[] distances=new ulong[points.length];
 	for(uint i=1; i<clusters.length; i++){
 		auto initialized=clusters[0..i];
 		ulong sum=0;
 		points
 			.map!(a=>cast(ulong)closest(initialized,a)[1]^^2)
 			.tee!(a=>sum+=a)
 			.copy(distances);
 		auto rand=uniform(0,sum);
 		foreach(ii,dist;enumerate(distances)){
 			if(rand<dist){
 				clusters[i].mean=points[ii];
 				break;
 			}
 			rand-=dist;
 		}
 	}
 	
 	void reset(){
 		foreach(ref c; clusters){
 			c.reset();
 		}
 	}
 	
 	foreach(iteration; iota(0,maxiter)){
 		reset();
 		points.each!(a=>addToClosest(clusters,a));
 		auto s=clusters.map!((ref a)=>a.calculate_delta()).sum;
 		if(s<1e-6){
 			return;
 		}
 	}
 }
 
 /**
  * This does k-means with increasing cluster sizes until the maximal 
  * std-deviation/mean ratio is below cutoff.
  */
 
 Cluster[] autokmeans(Point[] points, float reltol=0.1, float abstol=30*60, uint maxclusters=10){
 	Cluster[] res;
 	foreach(nc; iota(1,maxclusters+1)){
 		res=new Cluster[nc];
 		kmeans(points, res);
 		auto tol=res
 			.map!((a){auto stdev=sqrt(a.variance); return tuple(stdev/a.mean, stdev);})
 			.map!(a=>tuple(a[0]<reltol,a[1]<abstol))
 			.fold!((a,b)=>a = a && b[0] && b[1])(true);
 		if(tol){
 			return res;
 		}
 	}
 	return res;
 }
 unittest{
 	Cluster c;
 	c.points=[5,2];
 	c.calculate();
 	assert(c.mean==3.5);
 	assert(c.variance==2.25);
 }
 
 import std.file;
 
 int main(string[] args){
 	File f;
 	import std.getopt;
 	import std.format;
 	string plotfile;
 	auto helpWanted=getopt(args,
 		"plotfile|p", &plotfile,
 	);
 	if(args.length!=2){
 		stderr.writeln("Usage: ", args[0], " </path/to/log/file>");
 		return 1;
 	}
 	if(args[1]=="-"){
 		f=stdin;
 	}
 	else if(!exists(args[1])){
 		stderr.writeln(args[1], " does not exist");
 		return 1;
 	}
 	else{
 		f.open(args[1]);
 	}
 	
 	auto points=f
 		.byRecord!(uint, uint)("%s %s")
 		.map!(a=>a[1]-a[0])
 		.filter!(a=>a>10*60)
 		.filter!(a=>a<5*60*60)
 		.array;
 		
 	auto res=points.autokmeans();
 	res.sort!"a.mean<b.mean";
 	res
 		.filter!(a=>a.points.length>0)
 		.each!(a=>writeln(a.points.length, "\t", a.mean, "\t", a.variance));
 
 	if(plotfile){
 		File cf=File(plotfile,"w+");
 		scope(exit) cf.close();
 		foreach(i,c;enumerate(res)){
 			c.points.each!(a=>cf.writeln(i,"\t ",a));
 		}
 	}
 	return 0;
 }