[epclust.git] / epclust / tests / testthat / test.clustering.R

context("clustering")

#shorthand: map 1->1, 2->2, 3->3, 4->1, ..., 149->2, 150->3, ... (is base==3)
I = function(i, base)
	(i-1) %% base + 1

test_that("computeClusters1 behave as expected",
{
	require("MASS", quietly=TRUE)
	if (!require("clue", quietly=TRUE))
		skip("'clue' package not available")

	# 3 gaussian clusters, 300 items; and then 7 gaussian clusters, 490 items
	n = 300
	d = 5
	K = 3
	for (ndK in list( c(300,5,3), c(490,10,7) ))
	{
		n = ndK[1] ; d = ndK[2] ; K = ndK[3]
		cs = n/K #cluster size
		Id = diag(d)
		coefs = do.call(rbind,
			lapply(1:K, function(i) MASS::mvrnorm(cs, c(rep(0,(i-1)),5,rep(0,d-i)), Id)))
		indices_medoids = computeClusters1(coefs, K)
		# Get coefs assignments (to medoids)
		assignment = sapply(seq_len(n), function(i)
			which.min( rowSums( sweep(coefs[indices_medoids,],2,coefs[i,],'-')^2 ) ) )
		for (i in 1:K)
			expect_equal(sum(assignment==i), cs, tolerance=5)

		costs_matrix = matrix(nrow=K,ncol=K)
		for (i in 1:K)
		{
			for (j in 1:K)
			{
				# assign i (in result) to j (order 1,2,3)
				costs_matrix[i,j] = abs( mean(assignment[((i-1)*cs+1):(i*cs)]) - j )
			}
		}
		permutation = as.integer( clue::solve_LSAP(costs_matrix) )
		for (i in 1:K)
		{
			expect_equal(
				mean(assignment[((i-1)*cs+1):(i*cs)]), permutation[i], tolerance=0.05)
		}
	}
})

test_that("computeSynchrones behave as expected",
{
	n = 300
	x = seq(0,9.5,0.1)
	L = length(x) #96 1/4h
	K = 3
	s1 = cos(x)
	s2 = sin(x)
	s3 = c( s1[1:(L%/%2)] , s2[(L%/%2+1):L] )
	#sum((s1-s2)^2) == 96
	#sum((s1-s3)^2) == 58
	#sum((s2-s3)^2) == 38
	s = list(s1, s2, s3)
	series = matrix(nrow=n, ncol=L)
	for (i in seq_len(n))
		series[i,] = s[[I(i,K)]] + rnorm(L,sd=0.01)
	getRefSeries = function(indices) {
		indices = indices[indices <= n]
		if (length(indices)>0) series[indices,] else NULL
	}
	synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, n, 100,
		verbose=TRUE, parll=FALSE)

	expect_equal(dim(synchrones), c(K,L))
	for (i in 1:K)
		expect_equal(synchrones[i,], s[[i]], tolerance=0.01)
})

computeDistortion = function(series, medoids)
{
	n = nrow(series) ; L = ncol(series)
	distortion = 0.
	for (i in seq_len(n))
		distortion = distortion + min( rowSums( sweep(medoids,2,series[i,],'-')^2 ) / L )
	distortion / n
}

test_that("computeClusters2 behave as expected",
{
	n = 900
	x = seq(0,9.5,0.1)
	L = length(x) #96 1/4h
	K1 = 60
	K2 = 3
	#for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
	s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
	series = matrix(nrow=n, ncol=L)
	for (i in seq_len(n))
		series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
	getRefSeries = function(indices) {
		indices = indices[indices <= n]
		if (length(indices)>0) series[indices,] else NULL
	}
	# Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
	medoids_K1 = do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) )
	medoids_K2 = computeClusters2(medoids_K1, K2, getRefSeries, n, 75,
		verbose=TRUE, parll=FALSE)

	expect_equal(dim(medoids_K2), c(K2,L))
	# Not easy to evaluate result: at least we expect it to be better than random selection of
	# medoids within 1...K1 (among references)
	distorGood = computeDistortion(series, medoids_K2)
	for (i in 1:3)
		expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
})

test_that("clusteringTask1 + computeClusters2 behave as expected",
{
	n = 900
	x = seq(0,9.5,0.1)
	L = length(x) #96 1/4h
	K1 = 60
	K2 = 3
	s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
	series = matrix(nrow=n, ncol=L)
	for (i in seq_len(n))
		series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
	getSeries = function(indices) {
		indices = indices[indices <= n]
		if (length(indices)>0) series[indices,] else NULL
	}
	wf = "haar"
	ctype = "absolute"
	getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
	medoids_K1 = getSeries( clusteringTask1(1:n, getContribs, K1, 75,
		verbose=TRUE, parll=FALSE) )
	medoids_K2 = computeClusters2(medoids_K1, K2, getSeries, n, 120,
		verbose=TRUE, parll=FALSE)

	expect_equal(dim(medoids_K1), c(K1,L))
	expect_equal(dim(medoids_K2), c(K2,L))
	# Not easy to evaluate result: at least we expect it to be better than random selection of
	# medoids within 1...K1 (among references)
	distorGood = computeDistortion(series, medoids_K2)
	for (i in 1:3)
		expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
})
Commit	Line	Data
56857861 BA	1	context("clustering")
56857861 BA	2
8702eb86 BA	3	#shorthand: map 1->1, 2->2, 3->3, 4->1, ..., 149->2, 150->3, ... (is base==3)
	4	I = function(i, base)
	5	(i-1) %% base + 1
56857861 BA	6
	7	test_that("computeClusters1 behave as expected",
	8	{
8702eb86	9	require("MASS", quietly=TRUE)
4bcfdbee BA	10	if (!require("clue", quietly=TRUE))
4bcfdbee BA	11	skip("'clue' package not available")
56857861	12
8702eb86 BA	13	# 3 gaussian clusters, 300 items; and then 7 gaussian clusters, 490 items
	14	n = 300
	15	d = 5
	16	K = 3
	17	for (ndK in list( c(300,5,3), c(490,10,7) ))
	18	{
	19	n = ndK[1] ; d = ndK[2] ; K = ndK[3]
	20	cs = n/K #cluster size
	21	Id = diag(d)
	22	coefs = do.call(rbind,
	23	lapply(1:K, function(i) MASS::mvrnorm(cs, c(rep(0,(i-1)),5,rep(0,d-i)), Id)))
	24	indices_medoids = computeClusters1(coefs, K)
	25	# Get coefs assignments (to medoids)
	26	assignment = sapply(seq_len(n), function(i)
	27	which.min( rowSums( sweep(coefs[indices_medoids,],2,coefs[i,],'-')^2 ) ) )
	28	for (i in 1:K)
	29	expect_equal(sum(assignment==i), cs, tolerance=5)
	30
	31	costs_matrix = matrix(nrow=K,ncol=K)
	32	for (i in 1:K)
	33	{
	34	for (j in 1:K)
	35	{
	36	# assign i (in result) to j (order 1,2,3)
	37	costs_matrix[i,j] = abs( mean(assignment[((i-1)cs+1):(ics)]) - j )
	38	}
	39	}
	40	permutation = as.integer( clue::solve_LSAP(costs_matrix) )
	41	for (i in 1:K)
	42	{
	43	expect_equal(
	44	mean(assignment[((i-1)cs+1):(ics)]), permutation[i], tolerance=0.05)
	45	}
	46	}
56857861 BA	47	})
	48
	49	test_that("computeSynchrones behave as expected",
	50	{
8702eb86 BA	51	n = 300
	52	x = seq(0,9.5,0.1)
	53	L = length(x) #96 1/4h
	54	K = 3
	55	s1 = cos(x)
	56	s2 = sin(x)
	57	s3 = c( s1[1:(L%/%2)] , s2[(L%/%2+1):L] )
	58	#sum((s1-s2)^2) == 96
	59	#sum((s1-s3)^2) == 58
	60	#sum((s2-s3)^2) == 38
	61	s = list(s1, s2, s3)
	62	series = matrix(nrow=n, ncol=L)
	63	for (i in seq_len(n))
	64	series[i,] = s[[I(i,K)]] + rnorm(L,sd=0.01)
	65	getRefSeries = function(indices) {
492cd9e7	66	indices = indices[indices <= n]
8702eb86 BA	67	if (length(indices)>0) series[indices,] else NULL
8702eb86 BA	68	}
492cd9e7 BA	69	synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, n, 100,
492cd9e7 BA	70	verbose=TRUE, parll=FALSE)
56857861	71
8702eb86 BA	72	expect_equal(dim(synchrones), c(K,L))
	73	for (i in 1:K)
	74	expect_equal(synchrones[i,], s[[i]], tolerance=0.01)
56857861 BA	75	})
56857861 BA	76
8702eb86 BA	77	computeDistortion = function(series, medoids)
	78	{
	79	n = nrow(series) ; L = ncol(series)
	80	distortion = 0.
	81	for (i in seq_len(n))
	82	distortion = distortion + min( rowSums( sweep(medoids,2,series[i,],'-')^2 ) / L )
	83	distortion / n
	84	}
	85
56857861 BA	86	test_that("computeClusters2 behave as expected",
56857861 BA	87	{
8702eb86 BA	88	n = 900
	89	x = seq(0,9.5,0.1)
	90	L = length(x) #96 1/4h
	91	K1 = 60
	92	K2 = 3
	93	#for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
	94	s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
	95	series = matrix(nrow=n, ncol=L)
	96	for (i in seq_len(n))
	97	series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
	98	getRefSeries = function(indices) {
492cd9e7	99	indices = indices[indices <= n]
8702eb86 BA	100	if (length(indices)>0) series[indices,] else NULL
	101	}
	102	# Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
	103	medoids_K1 = do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) )
492cd9e7 BA	104	medoids_K2 = computeClusters2(medoids_K1, K2, getRefSeries, n, 75,
492cd9e7 BA	105	verbose=TRUE, parll=FALSE)
56857861	106
8702eb86 BA	107	expect_equal(dim(medoids_K2), c(K2,L))
	108	# Not easy to evaluate result: at least we expect it to be better than random selection of
	109	# medoids within 1...K1 (among references)
8702eb86 BA	110	distorGood = computeDistortion(series, medoids_K2)
	111	for (i in 1:3)
	112	expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
56857861 BA	113	})
56857861 BA	114
492cd9e7	115	test_that("clusteringTask1 + computeClusters2 behave as expected",
56857861	116	{
8702eb86 BA	117	n = 900
	118	x = seq(0,9.5,0.1)
	119	L = length(x) #96 1/4h
	120	K1 = 60
	121	K2 = 3
	122	s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
	123	series = matrix(nrow=n, ncol=L)
	124	for (i in seq_len(n))
	125	series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
	126	getSeries = function(indices) {
	127	indices = indices[indices <= n]
	128	if (length(indices)>0) series[indices,] else NULL
	129	}
	130	wf = "haar"
4bcfdbee BA	131	ctype = "absolute"
4bcfdbee BA	132	getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
492cd9e7 BA	133	medoids_K1 = getSeries( clusteringTask1(1:n, getContribs, K1, 75,
	134	verbose=TRUE, parll=FALSE) )
	135	medoids_K2 = computeClusters2(medoids_K1, K2, getSeries, n, 120,
	136	verbose=TRUE, parll=FALSE)
56857861	137
8702eb86 BA	138	expect_equal(dim(medoids_K1), c(K1,L))
	139	expect_equal(dim(medoids_K2), c(K2,L))
	140	# Not easy to evaluate result: at least we expect it to be better than random selection of
	141	# medoids within 1...K1 (among references)
	142	distorGood = computeDistortion(series, medoids_K2)
	143	for (i in 1:3)
	144	expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
56857861	145	})