'update'

[epclust.git] / epclust / R / stage2.R

#Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1

#(Benjamin)
#à partir de là, "conso" == courbes synchrones
n     <- nrow(conso)
delta <- ncol(conso)


#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]

#TODO: une fonction qui fait lignes 59 à 91

#cube:
# Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
#                 scalevector = scalevector4,
#                 lt = delta, smooth = FALSE, 
#                 nvoice = nvoice)      # observations node with CWT
# 
# #matrix:
# ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
# #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
# 
# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
# for(i in 1:n) 
#  Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) 
# 
# #rm(conso, Xcwt4); gc()
# 
# ## _.b WER^2 distances  ########
# Xwer_dist    <- matrix(0.0, n, n)
# for(i in 1:(n - 1)){
#  mat1   <- vect2mat(Xcwt2[i,])
#  for(j in (i + 1):n){
#     mat2 <- vect2mat(Xcwt2[j,])
#     num     <- Mod(mat1 * Conj(mat2))
#     WX      <- Mod(mat1 * Conj(mat1))
#     WY      <- Mod(mat2 * Conj(mat2))
#     smsmnum <- smCWT(num, scalevector = scalevector4)
#     smsmWX  <- smCWT(WX,  scalevector = scalevector4)
#     smsmWY  <- smCWT(WY,  scalevector = scalevector4)
#     wer2    <- sum(colSums(smsmnum)^2)  /
#       sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
#     Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
#     Xwer_dist[j, i] <- Xwer_dist[i, j]
#   }
# }
# diag(Xwer_dist) <- numeric(n)
# 
# save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
# save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")



#lignes 59 à 91 "dépliées" :
Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
                 scalevector = scalevector4,
                 lt = delta, smooth = FALSE, 
                 nvoice = nvoice)      # observations node with CWT
 
 #matrix:
 ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
 Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
 
 #NOTE: delta et lscvect pourraient etre gardés à part (communs)
 for(i in 1:n) 
    Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) 
 
 #rm(conso, Xcwt4); gc()
 
 ## _.b WER^2 distances  ########
 Xwer_dist    <- matrix(0.0, n, n)
 for(i in 1:(n - 1)){
  mat1   <- vect2mat(Xcwt2[i,])

 #NOTE: vect2mat = as.matrix ?! (dans aux.R)
  vect2mat <- function(vect){
                 vect <- as.vector(vect)
                 matrix(vect[-(1:2)], delta, lscvect)
               }
 
 for(j in (i + 1):n){
     mat2 <- vect2mat(Xcwt2[j,])
     num     <- Mod(mat1 * Conj(mat2))
     WX      <- Mod(mat1 * Conj(mat1))
     WY      <- Mod(mat2 * Conj(mat2))
     smsmnum <- smCWT(num, scalevector = scalevector4)
     smsmWX  <- smCWT(WX,  scalevector = scalevector4)
     smsmWY  <- smCWT(WY,  scalevector = scalevector4)
     wer2    <- sum(colSums(smsmnum)^2)  /
       sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
     Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
     Xwer_dist[j, i] <- Xwer_dist[i, j]
   }
 }
 diag(Xwer_dist) <- numeric(n)

#fonction smCWT (dans aux.R)
  smCWT <- function(CWT, sw=  0,  tw=  0, swabs= 0,
                       nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, 
					   lt= 24, dt= 0.5, scalevector )
		 {
#         noctave  <- adjust.noctave(lt, dt, s0, tw, noctave)
#         scalevector  <- 2^(0:(noctave * nvoice) / nvoice) * s0
         wsp     <- Mod(CWT)  
         smwsp   <- smooth.matrix(wsp, swabs)
         smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
         smsmwsp
       }

 #dans sowas.R
smooth.matrix <- function(wt,swabs){
  
  if (swabs != 0)
    smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1)))
  else
    smwt <- wt
  
  smwt
  
}
smooth.time <- function(wt,tw,dt,scalevector){
  
  smwt <- wt
  
  if (tw != 0){
    for (i in 1:length(scalevector)){
      
      twi <- as.integer(scalevector[i]*tw/dt)
      smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1))
      
    }
  } 
  smwt
}

#et filter() est dans stats::

#cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c