epclust/R/computeWerDists.R

   1 #' computeWerDists
   2 #'
   3 #' Compute the WER distances between the synchrones curves (in columns), which are
   4 #' returned (e.g.) by \code{computeSynchrones()}
   5 #'
   6 #' @param indices Range of series indices to cluster
   7 #' @inheritParams claws
   8 #' @inheritParams computeSynchrones
   9 #'
  10 #' @return A distances matrix of size K x K where K == length(indices)
  11 #'
  12 #' @export
  13 computeWerDists = function(indices, getSeries, nb_series_per_chunk, nvoice, nbytes, endian,
  14     ncores_clust=1, verbose=FALSE, parll=TRUE)
  15 {
  16     n <- length(indices)
  17     L <- length(getSeries(1)) #TODO: not very nice way to get L
  18     noctave = ceiling(log2(L)) #min power of 2 to cover serie range
  19     # Since a CWT contains noctave*nvoice complex series, we deduce the number of CWT to
  20     # retrieve/put in one chunk.
  21     nb_cwt_per_chunk = max(1, floor(nb_series_per_chunk / (nvoice*noctave*2)))
  22
  23     # Initialize result as a square big.matrix of size 'number of medoids'
  24     Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
  25
  26     # Generate n(n-1)/2 pairs for WER distances computations
  27     pairs = list()
  28     V = seq_len(n)
  29     for (i in 1:n)
  30     {
  31         V = V[-1]
  32         pairs = c(pairs, lapply(V, function(v) c(i,v)))
  33     }
  34
  35     cwt_file = ".cwt.bin"
  36     # Compute the getSeries(indices) CWT, and store the results in the binary file
  37     computeSaveCWT = function(indices)
  38     {
  39         if (parll)
  40         {
  41             require("bigmemory", quietly=TRUE)
  42             require("Rwave", quietly=TRUE)
  43             require("epclust", quietly=TRUE)
  44         }
  45
  46         # Obtain CWT as big vectors of real part + imaginary part (concatenate)
  47         ts_cwt <- sapply(indices, function(i) {
  48             ts <- scale(ts(getSeries(i)), center=TRUE, scale=FALSE)
  49             ts_cwt <- Rwave::cwt(ts, noctave, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE)
  50             c( as.double(Re(ts_cwt)),as.double(Im(ts_cwt)) )
  51         })
  52
  53         # Serialization
  54         binarize(ts_cwt, cwt_file, nb_cwt_per_chunk, ",", nbytes, endian)
  55     }
  56
  57     if (parll)
  58     {
  59         cl = parallel::makeCluster(ncores_clust)
  60         Xwer_dist_desc <- bigmemory::describe(Xwer_dist)
  61         parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","L",
  62             "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment())
  63     }
  64
  65     if (verbose)
  66         cat(paste("--- Precompute and serialize synchrones CWT\n", sep=""))
  67
  68     ignored <-
  69         if (parll)
  70             parallel::parLapply(cl, 1:n, computeSaveCWT)
  71         else
  72             lapply(1:n, computeSaveCWT)
  73
  74     # Function to retrieve a synchrone CWT from (binary) file
  75     getCWT = function(index, L)
  76     {
  77         flat_cwt <- getDataInFile(index, cwt_file, nbytes, endian)
  78         cwt_length = length(flat_cwt) / 2
  79         re_part = as.matrix(flat_cwt[1:cwt_length], nrow=L)
  80         im_part = as.matrix(flat_cwt[(cwt_length+1):(2*cwt_length)], nrow=L)
  81         re_part + 1i * im_part
  82     }
  83
  84
  85
  86
  87 #TODO: better repartition here,
  88
  89
  90
  91     # Compute distance between columns i and j in synchrones
  92     computeDistanceIJ = function(pair)
  93     {
  94         if (parll)
  95         {
  96             # parallel workers start with an empty environment
  97             require("bigmemory", quietly=TRUE)
  98             require("epclust", quietly=TRUE)
  99             Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc)
 100         }
 101
 102         i = pair[1] ; j = pair[2]
 103         if (verbose && j==i+1 && !parll)
 104             cat(paste("   Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
 105
 106         # Compute CWT of columns i and j in synchrones
 107         cwt_i <- getCWT(i, L)
 108         cwt_j <- getCWT(j, L)
 109
 110         # Compute the ratio of integrals formula 5.6 for WER^2
 111         # in https://arxiv.org/abs/1101.4744v2 §5.3
 112         num <- filterMA(Mod(cwt_i * Conj(cwt_j)))
 113         WX  <- filterMA(Mod(cwt_i * Conj(cwt_i)))
 114         WY <- filterMA(Mod(cwt_j * Conj(cwt_j)))
 115         wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
 116
 117         Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2))
 118         Xwer_dist[j,i] <- Xwer_dist[i,j]
 119         Xwer_dist[i,i] <- 0.
 120     }
 121
 122     if (verbose)
 123         cat(paste("--- Compute WER distances\n", sep=""))
 124
 125     ignored <-
 126         if (parll)
 127             parallel::parLapply(cl, pairs, computeDistanceIJ)
 128         else
 129             lapply(pairs, computeDistanceIJ)
 130
 131     if (parll)
 132         parallel::stopCluster(cl)
 133
 134     unlink(cwt_file)
 135
 136     Xwer_dist[n,n] = 0.
 137     Xwer_dist[,] #~small matrix K1 x K1
 138 }