epclust/tests/testthat/test.clustering.R

   1 context("clustering")
   2
   3 test_that("computeSynchrones behave as expected",
   4 {
   5     # Generate 300 sinusoïdal series of 3 kinds: all series of indices == 0 mod 3 are the same
   6     # (plus noise), all series of indices == 1 mod 3 are the same (plus noise) ...
   7     n = 300
   8     x = seq(0,9.5,0.1)
   9     L = length(x) #96 1/4h
  10     K = 3
  11     s1 = cos(x)
  12     s2 = sin(x)
  13     s3 = c( s1[1:(L%/%2)] , s2[(L%/%2+1):L] )
  14     #sum((s1-s2)^2) == 96
  15     #sum((s1-s3)^2) == 58
  16     #sum((s2-s3)^2) == 38
  17     s = list(s1, s2, s3)
  18     series = matrix(nrow=L, ncol=n)
  19     for (i in seq_len(n))
  20         series[,i] = s[[I(i,K)]] + rnorm(L,sd=0.01)
  21
  22     getRefSeries = function(indices) {
  23         indices = indices[indices <= n]
  24         if (length(indices)>0) as.matrix(series[,indices]) else NULL
  25     }
  26
  27     synchrones = computeSynchrones(bigmemory::as.big.matrix(cbind(s1,s2,s3)), getRefSeries,
  28         n, 100, verbose=TRUE, parll=FALSE)
  29
  30     expect_equal(dim(synchrones), c(L,K))
  31     for (i in 1:K)
  32     {
  33         # Synchrones are (for each medoid) sums of closest curves.
  34         # Here, we expect exactly 100 curves of each kind to be assigned respectively to
  35         # synchrone 1, 2 and 3 => division by 100 should be very close to the ref curve
  36         expect_equal(synchrones[,i]/100, s[[i]], tolerance=0.01)
  37     }
  38 })
  39
  40 test_that("Helper function to spread indices work properly",
  41 {
  42     indices <- 1:400
  43
  44     # bigger nb_per_set than length(indices)
  45     expect_equal(epclust:::.spreadIndices(indices,500), list(indices))
  46
  47     # nb_per_set == length(indices)
  48     expect_equal(epclust:::.spreadIndices(indices,400), list(indices))
  49
  50     # length(indices) %% nb_per_set == 0
  51     expect_equal(epclust:::.spreadIndices(indices,200),
  52         c( list(indices[1:200]), list(indices[201:400]) ))
  53     expect_equal(epclust:::.spreadIndices(indices,100),
  54         c( list(indices[1:100]), list(indices[101:200]),
  55             list(indices[201:300]), list(indices[301:400]) ))
  56
  57     # length(indices) / nb_per_set == 1, length(indices) %% nb_per_set == 100
  58     expect_equal(epclust:::.spreadIndices(indices,300), list(indices))
  59     # length(indices) / nb_per_set == 2, length(indices) %% nb_per_set == 42
  60     repartition <- epclust:::.spreadIndices(indices,179)
  61     expect_equal(length(repartition), 2)
  62     expect_equal(length(repartition[[1]]), 179 + 21)
  63     expect_equal(length(repartition[[1]]), 179 + 21)
  64 })
  65
  66 test_that("clusteringTask1 behave as expected",
  67 {
  68     # Generate 60 reference sinusoïdal series (medoids to be found),
  69     # and sample 900 series around them (add a small noise)
  70     n = 900
  71     x = seq(0,9.5,0.1)
  72     L = length(x) #96 1/4h
  73     K1 = 60
  74     s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
  75     series = matrix(nrow=L, ncol=n)
  76     for (i in seq_len(n))
  77         series[,i] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
  78
  79     getSeries = function(indices) {
  80         indices = indices[indices <= n]
  81         if (length(indices)>0) as.matrix(series[,indices]) else NULL
  82     }
  83
  84     wf = "haar"
  85     ctype = "absolute"
  86     getContribs = function(indices) curvesToContribs(series[,indices],wf,ctype)
  87
  88     require("cluster", quietly=TRUE)
  89     algoClust1 = function(contribs,K) cluster::pam(t(contribs),K,diss=FALSE)$id.med
  90     indices1 = clusteringTask1(1:n, getContribs, K1, algoClust1, 75, verbose=TRUE, parll=FALSE)
  91     medoids_K1 = getSeries(indices1)
  92
  93     expect_equal(dim(medoids_K1), c(L,K1))
  94     # Not easy to evaluate result: at least we expect it to be better than random selection of
  95     # medoids within initial series
  96     distor_good = computeDistortion(series, medoids_K1)
  97     for (i in 1:3)
  98         expect_lte( distor_good, computeDistortion(series,series[,sample(1:n, K1)]) )
  99 })
 100
 101 test_that("clusteringTask2 behave as expected",
 102 {
 103     skip("Unexplained failure")
 104
 105     # Same 60 reference sinusoïdal series than in clusteringTask1 test,
 106     # but this time we consider them as medoids - skipping stage 1
 107     # Here also we sample 900 series around the 60 "medoids"
 108     n = 900
 109     x = seq(0,9.5,0.1)
 110     L = length(x) #96 1/4h
 111     K1 = 60
 112     K2 = 3
 113     #for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
 114     s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
 115     series = matrix(nrow=L, ncol=n)
 116     for (i in seq_len(n))
 117         series[,i] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
 118
 119     getRefSeries = function(indices) {
 120         indices = indices[indices <= n]
 121         if (length(indices)>0) as.matrix(series[,indices]) else NULL
 122     }
 123
 124     # Perfect situation: all medoids "after stage 1" are good.
 125     medoids_K1 = bigmemory::as.big.matrix( sapply( 1:K1, function(i) s[[I(i,K1)]] ) )
 126     algoClust2 = function(dists,K) cluster::pam(dists,K,diss=TRUE)$id.med
 127     medoids_K2 = clusteringTask2(medoids_K1, K2, algoClust2, getRefSeries,
 128         n, 75, 4, 8, "little", verbose=TRUE, parll=FALSE)
 129
 130     expect_equal(dim(medoids_K2), c(L,K2))
 131     # Not easy to evaluate result: at least we expect it to be better than random selection of
 132     # synchrones within 1...K1 (from where distances computations + clustering was run)
 133     synchrones = computeSynchrones(medoids_K1,getRefSeries,n,75,verbose=FALSE,parll=FALSE)
 134     distor_good = computeDistortion(synchrones, medoids_K2)
 135     for (i in 1:3)
 136         expect_lte( distor_good, computeDistortion(synchrones, synchrones[,sample(1:K1,3)]) )
 137 })