Fix package, ok for R CMD check - ongoing debug for main function
[epclust.git] / epclust / tests / testthat / test.clustering.R
1 context("clustering")
2
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
6
7 test_that("computeClusters1 behave as expected",
8 {
9 require("MASS", quietly=TRUE)
10 if (!require("clue", quietly=TRUE))
11 skip("'clue' package not available")
12
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):(i*cs)]) - 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):(i*cs)]), permutation[i], tolerance=0.05)
45 }
46 }
47 })
48
49 test_that("computeSynchrones behave as expected",
50 {
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) {
66 indices = indices[indices < n]
67 if (length(indices)>0) series[indices,] else NULL
68 }
69 synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, 100)
70
71 expect_equal(dim(synchrones), c(K,L))
72 for (i in 1:K)
73 expect_equal(synchrones[i,], s[[i]], tolerance=0.01)
74 })
75
76 computeDistortion = function(series, medoids)
77 {
78 n = nrow(series) ; L = ncol(series)
79 distortion = 0.
80 for (i in seq_len(n))
81 distortion = distortion + min( rowSums( sweep(medoids,2,series[i,],'-')^2 ) / L )
82 distortion / n
83 }
84
85 test_that("computeClusters2 behave as expected",
86 {
87 n = 900
88 x = seq(0,9.5,0.1)
89 L = length(x) #96 1/4h
90 K1 = 60
91 K2 = 3
92 #for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
93 s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
94 series = matrix(nrow=n, ncol=L)
95 for (i in seq_len(n))
96 series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
97 getRefSeries = function(indices) {
98 indices = indices[indices < n]
99 if (length(indices)>0) series[indices,] else NULL
100 }
101 # Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
102 medoids_K1 = do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) )
103 medoids_K2 = computeClusters2(medoids_K1, K2, getRefSeries, 75)
104
105 expect_equal(dim(medoids_K2), c(K2,L))
106 # Not easy to evaluate result: at least we expect it to be better than random selection of
107 # medoids within 1...K1 (among references)
108
109 distorGood = computeDistortion(series, medoids_K2)
110 for (i in 1:3)
111 expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
112 })
113
114 test_that("clusteringTask + computeClusters2 behave as expected",
115 {
116 n = 900
117 x = seq(0,9.5,0.1)
118 L = length(x) #96 1/4h
119 K1 = 60
120 K2 = 3
121 s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
122 series = matrix(nrow=n, ncol=L)
123 for (i in seq_len(n))
124 series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
125 getSeries = function(indices) {
126 indices = indices[indices <= n]
127 if (length(indices)>0) series[indices,] else NULL
128 }
129 wf = "haar"
130 ctype = "absolute"
131 getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
132 medoids_K1 = getSeries( clusteringTask(1:n, getContribs, K1, 75, 4) )
133 medoids_K2 = computeClusters2(medoids_K1, K2, getSeries, 120)
134
135 expect_equal(dim(medoids_K1), c(K1,L))
136 expect_equal(dim(medoids_K2), c(K2,L))
137 # Not easy to evaluate result: at least we expect it to be better than random selection of
138 # medoids within 1...K1 (among references)
139 distorGood = computeDistortion(series, medoids_K2)
140 for (i in 1:3)
141 expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
142 })