epclust/R/stage2.R

   1 library("Rwave")
   2
   3 #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1
   4 step2 = function(conso)
   5 {
   6     n <- nrow(conso)
   7     delta <- ncol(conso)
   8     #TODO: automatic tune of all these parameters ? (for other users)
   9     nvoice   <- 4
  10     # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(conso))
  11     noctave = 13
  12     # 4 here represent 2^5 = 32 half-hours ~ 1 day
  13     #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?)
  14     scalevector  <- 2^(4:(noctave * nvoice) / nvoice) * 2
  15     #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
  16     s0=2
  17     w0=2*pi
  18     scaled=FALSE
  19     s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
  20     totnoct = noctave + as.integer(s0log/nvoice) + 1
  21
  22     # (normalized) observations node with CWT
  23     Xcwt4 <- lapply(seq_len(n), function(i) {
  24         ts <- scale(ts(conso[i,]), center=TRUE, scale=scaled)
  25         totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
  26         ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
  27         #Normalization
  28         sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
  29         sqres <- sweep(ts.cwt,MARGIN=2,sqs,'*')
  30         sqres / max(Mod(sqres))
  31     })
  32
  33     Xwer_dist <- matrix(0., n, n)
  34     fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!)
  35     for (i in 1:(n-1))
  36     {
  37         for (j in (i+1):n)
  38         {
  39             #TODO: later, compute CWT here (because not enough storage space for 32M series)
  40             #      'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
  41             num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
  42             WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
  43             WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
  44             wer2    <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) )
  45             Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2))
  46             Xwer_dist[j,i] <- Xwer_dist[i,j]
  47         }
  48     }
  49     diag(Xwer_dist) <- numeric(n)
  50     Xwer_dist
  51 }