X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=epclust%2FR%2Fstage2.R;h=3ccbbad10903b256c005b0476b418bf3bbfdd9a3;hb=db6fc17ddd53fb0c64cf957296dc615ba830db56;hp=f952da2e5dfce72493c757c2d120947c747fe979;hpb=dc1aa85a96bbf815b0d896c22a9b4a539a9e8a9c;p=epclust.git diff --git a/epclust/R/stage2.R b/epclust/R/stage2.R index f952da2..3ccbbad 100644 --- a/epclust/R/stage2.R +++ b/epclust/R/stage2.R @@ -1,139 +1,51 @@ -#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") - - +library("Rwave") -#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 +#Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 +step2 = function(conso) +{ + n <- nrow(conso) + delta <- ncol(conso) + #TODO: automatic tune of all these parameters ? (for other users) + nvoice <- 4 + # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(conso)) + noctave = 13 + # 4 here represent 2^5 = 32 half-hours ~ 1 day + #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?) + scalevector <- 2^(4:(noctave * nvoice) / nvoice) * 2 + #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1 + s0=2 + w0=2*pi + scaled=FALSE + s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 ) + totnoct = noctave + as.integer(s0log/nvoice) + 1 + + # (normalized) observations node with CWT + Xcwt4 <- lapply(seq_len(n), function(i) { + ts <- scale(ts(conso[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,MARGIN=2,sqs,'*') + sqres / max(Mod(sqres)) + }) + + Xwer_dist <- matrix(0., n, n) + fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!) + for (i in 1:(n-1)) + { + for (j in (i+1):n) + { + #TODO: later, compute CWT here (because not enough storage space for 32M series) + # '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)) ) + Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2)) + Xwer_dist[j,i] <- Xwer_dist[i,j] + } + } + diag(Xwer_dist) <- numeric(n) + Xwer_dist } - -#et filter() est dans stats:: - -#cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c -