#ignore jupyter generated file (HTML vignette, and reports)
*.ipynb.html
+
+#ignore object files
+*.o
+*.so
DBI
License: MIT + file LICENSE
RoxygenNote: 6.0.1
-Collate:
+Collate:
'main.R'
'clustering.R'
'de_serialize.R'
'a_NAMESPACE.R'
+ 'plot.R'
#' @include clustering.R
#' @include main.R
#'
+#' @useDynLib epclust
+#'
#' @importFrom Rwave cwt
#' @importFrom cluster pam
#' @importFrom parallel makeCluster clusterExport parLapply stopCluster
#' @importFrom wavelets dwt wt.filter
-#' @importFrom stats filter spline
-#' @importFrom utils tail
+#' @importFrom stats spline
#' @importFrom methods is
-#' @importFrom bigmemory as.big.matrix is.big.matrix
+#' @importFrom bigmemory big.matrix as.big.matrix is.big.matrix
NULL
indices, getContribs, K1, nb_series_per_chunk, ncores_clust=1, verbose=FALSE, parll=TRUE)
{
if (verbose)
- cat(paste("*** Clustering task on ",length(indices)," lines\n", sep=""))
-
- wrapComputeClusters1 = function(inds) {
- if (parll)
- require("epclust", quietly=TRUE)
- if (verbose)
- cat(paste(" computeClusters1() on ",length(inds)," lines\n", sep=""))
- inds[ computeClusters1(getContribs(inds), K1) ]
- }
+ cat(paste("*** Clustering task 1 on ",length(indices)," lines\n", sep=""))
if (parll)
{
while (length(indices) > K1)
{
indices_workers = .spreadIndices(indices, nb_series_per_chunk)
- if (parll)
- indices = unlist( parallel::parLapply(cl, indices_workers, wrapComputeClusters1) )
- else
- indices = unlist( lapply(indices_workers, wrapComputeClusters1) )
+ indices <-
+ if (parll)
+ {
+ unlist( parallel::parLapply(cl, indices_workers, function(inds) {
+ require("epclust", quietly=TRUE)
+ inds[ computeClusters1(getContribs(inds), K1, verbose) ]
+ }) )
+ }
+ else
+ {
+ unlist( lapply(indices_workers, function(inds)
+ inds[ computeClusters1(getContribs(inds), K1, verbose) ]
+ ) )
+ }
}
if (parll)
parallel::stopCluster(cl)
clusteringTask2 = function(medoids, K2,
getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
{
+ if (verbose)
+ cat(paste("*** Clustering task 2 on ",nrow(medoids)," lines\n", sep=""))
+
if (nrow(medoids) <= K2)
return (medoids)
synchrones = computeSynchrones(medoids,
getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll)
distances = computeWerDists(synchrones, ncores_clust, verbose, parll)
# PAM in package 'cluster' cannot take big.matrix in input: need to cast it
- mat_dists = matrix(nrow=K1, ncol=K1)
- for (i in seq_len(K1))
- mat_dists[i,] = distances[i,]
- medoids[ computeClusters2(mat_dists,K2), ]
+ medoids[ computeClusters2(distances[,],K2,verbose), ]
}
#' @rdname clustering
#' @export
-computeClusters1 = function(contribs, K1)
+computeClusters1 = function(contribs, K1, verbose=FALSE)
+{
+ if (verbose)
+ cat(paste(" computeClusters1() on ",nrow(contribs)," lines\n", sep=""))
cluster::pam(contribs, K1, diss=FALSE)$id.med
+}
#' @rdname clustering
#' @export
-computeClusters2 = function(distances, K2)
+computeClusters2 = function(distances, K2, verbose=FALSE)
+{
+ if (verbose)
+ cat(paste(" computeClusters2() on ",nrow(distances)," lines\n", sep=""))
cluster::pam(distances, K2, diss=TRUE)$id.med
+}
#' computeSynchrones
#'
computeSynchrones = function(medoids, getRefSeries,
nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
{
+ if (verbose)
+ cat(paste("--- Compute synchrones\n", sep=""))
+
computeSynchronesChunk = function(indices)
{
- if (verbose)
- cat(paste("--- Compute synchrones for ",length(indices)," lines\n", sep=""))
ref_series = getRefSeries(indices)
+ nb_series = nrow(ref_series)
#get medoids indices for this chunk of series
- for (i in seq_len(nrow(ref_series)))
+
+ #TODO: debug this (address is OK but values are garbage: why?)
+# mi = .Call("computeMedoidsIndices", medoids@address, ref_series, PACKAGE="epclust")
+
+ #R-equivalent, requiring a matrix (thus potentially breaking "fit-in-memory" hope)
+ mat_meds = medoids[,]
+ mi = rep(NA,nb_series)
+ for (i in 1:nb_series)
+ mi[i] <- which.min( rowSums( sweep(mat_meds, 2, ref_series[i,], '-')^2 ) )
+ rm(mat_meds); gc()
+
+ for (i in seq_len(nb_series))
{
- j = which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
if (parll)
synchronicity::lock(m)
- synchrones[j,] = synchrones[j,] + ref_series[i,]
- counts[j,1] = counts[j,1] + 1
+ synchrones[mi[i],] = synchrones[mi[i],] + ref_series[i,]
+ counts[mi[i],1] = counts[mi[i],1] + 1
if (parll)
synchronicity::unlock(m)
}
}
- K = nrow(medoids)
+ K = nrow(medoids) ; L = ncol(medoids)
# Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in //
# TODO: if size > RAM (not our case), use file-backed big.matrix
- synchrones = bigmemory::big.matrix(nrow=K,ncol=ncol(medoids),type="double",init=0.)
- counts = bigmemory::big.matrix(nrow=K,ncol=1,type="double",init=0)
+ synchrones = bigmemory::big.matrix(nrow=K, ncol=L, type="double", init=0.)
+ counts = bigmemory::big.matrix(nrow=K, ncol=1, type="double", init=0)
# synchronicity is only for Linux & MacOS; on Windows: run sequentially
parll = (requireNamespace("synchronicity",quietly=TRUE)
&& parll && Sys.info()['sysname'] != "Windows")
}
indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk)
+ browser()
ignored <-
if (parll)
- parallel::parLapply(indices_workers, computeSynchronesChunk)
+ parallel::parLapply(cl, indices_workers, computeSynchronesChunk)
else
lapply(indices_workers, computeSynchronesChunk)
#' @export
computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
{
-
-
-
-#TODO: re-organize to call computeWerDist(x,y) [C] (in //?) from two indices + big.matrix
-
+ if (verbose)
+ cat(paste("--- Compute WER dists\n", sep=""))
n <- nrow(synchrones)
delta <- ncol(synchrones)
s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
totnoct = noctave + as.integer(s0log/nvoice) + 1
+ Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
+ fcoefs = rep(1/3, 3) #moving average on 3 values
+
+ # Generate n(n-1)/2 pairs for WER distances computations
+ pairs = list()
+ V = seq_len(n)
+ for (i in 1:n)
+ {
+ V = V[-1]
+ pairs = c(pairs, lapply(V, function(v) c(i,v)))
+ }
+
computeCWT = function(i)
{
- if (verbose)
- cat(paste("+++ Compute Rwave::cwt() on serie ",i,"\n", sep=""))
ts <- scale(ts(synchrones[i,]), center=TRUE, scale=scaled)
totts.cwt = Rwave::cwt(ts, totnoct, nvoice, w0, plot=FALSE)
ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
sqres / max(Mod(sqres))
}
+ computeDistancesIJ = function(pair)
+ {
+ i = pair[1] ; j = pair[2]
+ if (verbose && j==i+1)
+ cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
+ cwt_i = computeCWT(i)
+ cwt_j = computeCWT(j)
+ num <- .Call("filter", Mod(cwt_i * Conj(cwt_j)), PACKAGE="epclust")
+ WX <- .Call("filter", Mod(cwt_i * Conj(cwt_i)), PACKAGE="epclust")
+ WY <- .Call("filter", Mod(cwt_j * Conj(cwt_j)), PACKAGE="epclust")
+ wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
+ Xwer_dist[i,j] <- sqrt(delta * ncol(cwt_i) * (1 - wer2))
+ Xwer_dist[j,i] <- Xwer_dist[i,j]
+ Xwer_dist[i,i] = 0.
+ }
+
if (parll)
{
cl = parallel::makeCluster(ncores_clust)
envir=environment())
}
- # list of CWT from synchrones
- # TODO: fit in RAM, OK? If not, 2 options: serialize, compute individual distances
- Xcwt4 <-
+ ignored <-
if (parll)
- parallel::parLapply(cl, seq_len(n), computeCWT)
+ parallel::parLapply(cl, pairs, computeDistancesIJ)
else
- lapply(seq_len(n), computeCWT)
+ lapply(pairs, computeDistancesIJ)
if (parll)
parallel::stopCluster(cl)
- Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
- fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!)
- if (verbose)
- cat("*** Compute WER distances from CWT\n")
-
- #TODO: computeDistances(i,j), et répartir les n(n-1)/2 couples d'indices
- #là c'est trop déséquilibré
-
- computeDistancesLineI = function(i)
- {
- if (verbose)
- cat(paste(" Line ",i,"\n", sep=""))
- for (j in (i+1):n)
- {
- #TODO: 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
- num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
- WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
- WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
- wer2 <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) )
- if (parll)
- synchronicity::lock(m)
- Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2))
- Xwer_dist[j,i] <- Xwer_dist[i,j]
- if (parll)
- synchronicity::unlock(m)
- }
- Xwer_dist[i,i] = 0.
- }
-
- parll = (requireNamespace("synchronicity",quietly=TRUE)
- && parll && Sys.info()['sysname'] != "Windows")
- if (parll)
- m <- synchronicity::boost.mutex()
-
- ignored <-
- if (parll)
- {
- parallel::mclapply(seq_len(n-1), computeDistancesLineI,
- mc.cores=ncores_clust, mc.allow.recursive=FALSE)
- }
- else
- lapply(seq_len(n-1), computeDistancesLineI)
Xwer_dist[n,n] = 0.
Xwer_dist
}
indices_all[((i-1)*nb_series_per_task+1):upper_bound]
})
if (verbose)
- cat(paste("...Run ",ntasks," x stage 1 in parallel\n",sep=""))
+ {
+ message = paste("...Run ",ntasks," x stage 1", sep="")
+ if (WER=="mix")
+ message = paste(message," + stage 2", sep="")
+ cat(paste(message,"\n", sep=""))
+ }
if (WER=="mix")
{synchrones_file = paste(bin_dir,"synchrones",sep="") ; unlink(synchrones_file)}
if (parll && ntasks>1)
indices_medoids = clusteringTask1(
indices, getContribs, K1, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose, parll)
medoids1 = bigmemory::as.big.matrix( getSeries(indices_medoids) )
- medoids2 = computeClusters2(medoids1, K2,
+ medoids2 = clusteringTask2(medoids1, K2,
getRefSeries, nb_curves, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose, parll)
# Cleanup
+++ /dev/null
-#include <stdlib.h>
-#include <math.h>
-#include <stdbool.h>
-
-#ifndef M_PI
-#define M_PI 3.14159265358979323846
-#endif
-
-// n: number of synchrones, m: length of a synchrone
-float computeWerDist(float* s1, float* s2, int n, int m)
-{
- //TODO: automatic tune of all these parameters ? (for other users)
- int nvoice = 4;
- //noctave 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones))
- int noctave = 13
- // 4 here represent 2^5 = 32 half-hours ~ 1 day
- //NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?)
- //R: scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1)
- int* scalevector = (int*)malloc( (noctave*nvoice-4 + 1) * sizeof(int))
- for (int i=4; i<=noctave*nvoice; i++)
- scalevector[i-4] = pow(2., (float)i/nvoice + 1.);
- //condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
- int s0 = 2;
- double w0 = 2*M_PI;
- bool scaled = false;
- int s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
- int totnoct = noctave + as.integer(s0log/nvoice) + 1
-
-
-
-
-
-///TODO: continue
-
-
-
- computeCWT = function(i)
- {
- if (verbose)
- cat(paste("+++ Compute Rwave::cwt() on serie ",i,"\n", sep=""))
- ts <- scale(ts(synchrones[i,]), center=TRUE, scale=scaled)
- totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
- ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
- #Normalization
- sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
- sqres <- sweep(ts.cwt,2,sqs,'*')
- sqres / max(Mod(sqres))
- }
-
- if (parll)
- {
- cl = parallel::makeCluster(ncores_clust)
- parallel::clusterExport(cl,
- varlist=c("synchrones","totnoct","nvoice","w0","s0log","noctave","s0","verbose"),
- envir=environment())
- }
-
- # (normalized) observations node with CWT
- Xcwt4 <-
- if (parll)
- parallel::parLapply(cl, seq_len(n), computeCWT)
- else
- lapply(seq_len(n), computeCWT)
-
- if (parll)
- parallel::stopCluster(cl)
-
- Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
- fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!)
- if (verbose)
- cat("*** Compute WER distances from CWT\n")
-
- #TODO: computeDistances(i,j), et répartir les n(n-1)/2 couples d'indices
- #là c'est trop déséquilibré
-
- computeDistancesLineI = function(i)
- {
- if (verbose)
- cat(paste(" Line ",i,"\n", sep=""))
- for (j in (i+1):n)
- {
- #TODO: 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
- num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
- WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
- WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
- wer2 <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) )
- if (parll)
- synchronicity::lock(m)
- Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2))
- Xwer_dist[j,i] <- Xwer_dist[i,j]
- if (parll)
- synchronicity::unlock(m)
- }
- Xwer_dist[i,i] = 0.
- }
-
- parll = (requireNamespace("synchronicity",quietly=TRUE)
- && parll && Sys.info()['sysname'] != "Windows")
- if (parll)
- m <- synchronicity::boost.mutex()
-
- ignored <-
- if (parll)
- {
- parallel::mclapply(seq_len(n-1), computeDistancesLineI,
- mc.cores=ncores_clust, mc.allow.recursive=FALSE)
- }
- else
- lapply(seq_len(n-1), computeDistancesLineI)
- Xwer_dist[n,n] = 0.
-
- mat_dists = matrix(nrow=n, ncol=n)
- #TODO: avoid this loop?
- for (i in 1:n)
- mat_dists[i,] = Xwer_dist[i,]
- mat_dists
-
--- /dev/null
+#include <stdlib.h>
+#include <float.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+
+#include <stdio.h>
+
+// (K,L): dim(medoids)
+// mi: medoids indices
+SEXP computeMedoidsIndices(SEXP medoids_, SEXP ref_series_)
+{
+ double *medoids = (double*) R_ExternalPtrAddr(medoids_),
+ *ref_series = REAL(ref_series_);
+ int nb_series = INTEGER(getAttrib(ref_series_, R_DimSymbol))[0],
+ K = INTEGER(getAttrib(medoids_, R_DimSymbol))[0],
+ L = INTEGER(getAttrib(ref_series_, R_DimSymbol))[1],
+ *mi = (int*)malloc(nb_series*sizeof(int));
+
+ //TODO: debug this: medoids have same addresses on both sides, but this side fails
+ printf("MED: %x\n",medoids);
+
+ for (int i=0; i<nb_series ; i++)
+ {
+ // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
+ mi[i] = 0;
+ double best_dist = DBL_MAX;
+ for (int j=0; j<K; j++)
+ {
+ double dist_ij = 0.;
+ for (int k=0; k<L; k++)
+ {
+ double delta = ref_series[k*K+i] - medoids[k*K+j];
+ dist_ij += delta * delta;
+ }
+ if (dist_ij < best_dist)
+ {
+ mi[i] = j+1; //R indices start at 1
+ best_dist = dist_ij;
+ }
+ }
+ }
+
+ SEXP R_mi;
+ PROTECT(R_mi = allocVector(INTSXP, nb_series));
+ for (int i = 0; i < nb_series; i++)
+ INTEGER(R_mi)[i] = mi[i];
+ UNPROTECT(1);
+ free(mi);
+ return R_mi;
+}
--- /dev/null
+#include <stdlib.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+
+#include <stdio.h>
+
+SEXP filter(SEXP cwt_)
+{
+ int L = INTEGER(getAttrib(cwt_, R_DimSymbol))[0],
+ D = INTEGER(getAttrib(cwt_, R_DimSymbol))[1];
+ double *cwt = REAL(cwt_);
+ SEXP R_fcwt;
+ PROTECT(R_fcwt = allocMatrix(REALSXP, L, D));
+ double* fcwt = REAL(R_fcwt);
+
+ //TODO: coding style is terrible... no time for now.
+ for (int col=0; col<D; col++)
+ {
+ double v1 = cwt[0];
+ double ma = v1 + cwt[1] + cwt[2];
+ for (int i=1; i<L-2; i++)
+ {
+ fcwt[i] = ma / 3.;
+ ma = ma - v1 + cwt[i+2];
+ v1 = cwt[i];
+ }
+ fcwt[0] = cwt[0];
+ fcwt[L-2] = ma / 3.;
+ fcwt[L-1] = cwt[L-1];
+ cwt = cwt + L;
+ fcwt = fcwt + L;
+ }
+
+ UNPROTECT(1);
+ return R_fcwt;
+}
I = function(i, base)
(i-1) %% base + 1
-test_that("computeClusters1 behave as expected",
+test_that("computeClusters1&2 behave as expected",
{
require("MASS", quietly=TRUE)
if (!require("clue", quietly=TRUE))
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)
+ indices_medoids1 = computeClusters1(coefs, K, verbose=TRUE)
+ indices_medoids2 = computeClusters2(dist(coefs), K, verbose=TRUE)
# Get coefs assignments (to medoids)
- assignment = sapply(seq_len(n), function(i)
- which.min( rowSums( sweep(coefs[indices_medoids,],2,coefs[i,],'-')^2 ) ) )
+ assignment1 = sapply(seq_len(n), function(i)
+ which.min( rowSums( sweep(coefs[indices_medoids1,],2,coefs[i,],'-')^2 ) ) )
+ assignment2 = sapply(seq_len(n), function(i)
+ which.min( rowSums( sweep(coefs[indices_medoids2,],2,coefs[i,],'-')^2 ) ) )
for (i in 1:K)
- expect_equal(sum(assignment==i), cs, tolerance=5)
+ {
+ expect_equal(sum(assignment1==i), cs, tolerance=5)
+ expect_equal(sum(assignment2==i), cs, tolerance=5)
+ }
- costs_matrix = matrix(nrow=K,ncol=K)
+ costs_matrix1 = matrix(nrow=K,ncol=K)
+ costs_matrix2 = 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 )
+ costs_matrix1[i,j] = abs( mean(assignment1[((i-1)*cs+1):(i*cs)]) - j )
+ costs_matrix2[i,j] = abs( mean(assignment2[((i-1)*cs+1):(i*cs)]) - j )
}
}
- permutation = as.integer( clue::solve_LSAP(costs_matrix) )
+ permutation1 = as.integer( clue::solve_LSAP(costs_matrix1) )
+ permutation2 = as.integer( clue::solve_LSAP(costs_matrix2) )
for (i in 1:K)
{
expect_equal(
- mean(assignment[((i-1)*cs+1):(i*cs)]), permutation[i], tolerance=0.05)
+ mean(assignment1[((i-1)*cs+1):(i*cs)]), permutation1[i], tolerance=0.05)
+ expect_equal(
+ mean(assignment2[((i-1)*cs+1):(i*cs)]), permutation2[i], tolerance=0.05)
}
}
})
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)
+ synchrones = computeSynchrones(bigmemory::as.big.matrix(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)
})
+# NOTE: medoids can be a big.matrix
computeDistortion = function(series, medoids)
{
n = nrow(series) ; L = ncol(series)
distortion = 0.
+ if (bigmemory::is.big.matrix(medoids))
+ medoids = medoids[,]
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",
+test_that("clusteringTask1 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) {
+ getSeries = 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)
+ wf = "haar"
+ ctype = "absolute"
+ getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
+ indices1 = clusteringTask1(1:n, getContribs, K1, 75, verbose=TRUE, parll=FALSE)
+ medoids_K1 = getSeries(indices1)
- expect_equal(dim(medoids_K2), c(K2,L))
+ expect_equal(dim(medoids_K1), c(K1,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)
+ # medoids within initial series
+ distorGood = computeDistortion(series, medoids_K1)
for (i in 1:3)
- expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
+ expect_lte( distorGood, computeDistortion(series,series[sample(1:n, K1),]) )
})
-test_that("clusteringTask1 + computeClusters2 behave as expected",
+test_that("clusteringTask2 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)
- getSeries = function(indices) {
+ getRefSeries = 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)
+ # Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
+ medoids_K1 = bigmemory::as.big.matrix(
+ do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) ) )
+ medoids_K2 = clusteringTask2(medoids_K1, K2, getRefSeries, n, 75, 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)
for (i in 1:3)
expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
})
+
+#NOTE: rather redundant test
+#test_that("clusteringTask1 + clusteringTask2 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)
+# require("bigmemory", quietly=TRUE)
+# indices1 = clusteringTask1(1:n, getContribs, K1, 75, verbose=TRUE, parll=FALSE)
+# medoids_K1 = bigmemory::as.big.matrix( getSeries(indices1) )
+# medoids_K2 = clusteringTask2(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),]) )
+#})
projects / epclust.git / commitdiff