epclust/R/stage2.R

   1 #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1
   2
   3 #(Benjamin)
   4 #à partir de là, "conso" == courbes synchrones
   5 n     <- nrow(conso)
   6 delta <- ncol(conso)
   7
   8
   9 #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]
  10
  11 #TODO: une fonction qui fait lignes 59 à 91
  12
  13 #cube:
  14 # Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
  15 #                 scalevector = scalevector4,
  16 #                 lt = delta, smooth = FALSE,
  17 #                 nvoice = nvoice)      # observations node with CWT
  18 #
  19 # #matrix:
  20 # ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
  21 # #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
  22 #
  23 # #NOTE: delta et lscvect pourraient etre gardés à part (communs)
  24 # for(i in 1:n)
  25 #  Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
  26 #
  27 # #rm(conso, Xcwt4); gc()
  28 #
  29 # ## _.b WER^2 distances  ########
  30 # Xwer_dist    <- matrix(0.0, n, n)
  31 # for(i in 1:(n - 1)){
  32 #  mat1   <- vect2mat(Xcwt2[i,])
  33 #  for(j in (i + 1):n){
  34 #     mat2 <- vect2mat(Xcwt2[j,])
  35 #     num     <- Mod(mat1 * Conj(mat2))
  36 #     WX      <- Mod(mat1 * Conj(mat1))
  37 #     WY      <- Mod(mat2 * Conj(mat2))
  38 #     smsmnum <- smCWT(num, scalevector = scalevector4)
  39 #     smsmWX  <- smCWT(WX,  scalevector = scalevector4)
  40 #     smsmWY  <- smCWT(WY,  scalevector = scalevector4)
  41 #     wer2    <- sum(colSums(smsmnum)^2)  /
  42 #       sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
  43 #     Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
  44 #     Xwer_dist[j, i] <- Xwer_dist[i, j]
  45 #   }
  46 # }
  47 # diag(Xwer_dist) <- numeric(n)
  48 #
  49 # save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
  50 # save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")
  51
  52
  53
  54 #lignes 59 à 91 "dépliées" :
  55 Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
  56                  scalevector = scalevector4,
  57                  lt = delta, smooth = FALSE,
  58                  nvoice = nvoice)      # observations node with CWT
  59
  60  #matrix:
  61  ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
  62  Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
  63
  64  #NOTE: delta et lscvect pourraient etre gardés à part (communs)
  65  for(i in 1:n)
  66     Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
  67
  68  #rm(conso, Xcwt4); gc()
  69
  70  ## _.b WER^2 distances  ########
  71  Xwer_dist    <- matrix(0.0, n, n)
  72  for(i in 1:(n - 1)){
  73   mat1   <- vect2mat(Xcwt2[i,])
  74
  75  #NOTE: vect2mat = as.matrix ?! (dans aux.R)
  76   vect2mat <- function(vect){
  77                  vect <- as.vector(vect)
  78                  matrix(vect[-(1:2)], delta, lscvect)
  79                }
  80
  81  for(j in (i + 1):n){
  82      mat2 <- vect2mat(Xcwt2[j,])
  83      num     <- Mod(mat1 * Conj(mat2))
  84      WX      <- Mod(mat1 * Conj(mat1))
  85      WY      <- Mod(mat2 * Conj(mat2))
  86      smsmnum <- smCWT(num, scalevector = scalevector4)
  87      smsmWX  <- smCWT(WX,  scalevector = scalevector4)
  88      smsmWY  <- smCWT(WY,  scalevector = scalevector4)
  89      wer2    <- sum(colSums(smsmnum)^2)  /
  90        sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
  91      Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
  92      Xwer_dist[j, i] <- Xwer_dist[i, j]
  93    }
  94  }
  95  diag(Xwer_dist) <- numeric(n)
  96
  97 #fonction smCWT (dans aux.R)
  98   smCWT <- function(CWT, sw=  0,  tw=  0, swabs= 0,
  99                        nvoice= 12, noctave= 2, s0= 2, w0= 2*pi,
 100                        lt= 24, dt= 0.5, scalevector )
 101          {
 102 #         noctave  <- adjust.noctave(lt, dt, s0, tw, noctave)
 103 #         scalevector  <- 2^(0:(noctave * nvoice) / nvoice) * s0
 104          wsp     <- Mod(CWT)
 105          smwsp   <- smooth.matrix(wsp, swabs)
 106          smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
 107          smsmwsp
 108        }
 109
 110  #dans sowas.R
 111 smooth.matrix <- function(wt,swabs){
 112
 113   if (swabs != 0)
 114     smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1)))
 115   else
 116     smwt <- wt
 117
 118   smwt
 119
 120 }
 121 smooth.time <- function(wt,tw,dt,scalevector){
 122
 123   smwt <- wt
 124
 125   if (tw != 0){
 126     for (i in 1:length(scalevector)){
 127
 128       twi <- as.integer(scalevector[i]*tw/dt)
 129       smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1))
 130
 131     }
 132   }
 133   smwt
 134 }
 135
 136 #et filter() est dans stats::
 137
 138 #cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c
 139