From d7d55bc1e74711b0da84578ecdebc43eeb259599 Mon Sep 17 00:00:00 2001 From: Benjamin Auder <benjamin.auder@somewhere> Date: Sat, 18 Feb 2017 13:12:01 +0100 Subject: [PATCH] simplify stage2.R --- TODO | 7 +++ epclust/R/main.R | 1 + epclust/R/stage2.R | 138 +++++++++++---------------------------------- 3 files changed, 41 insertions(+), 105 deletions(-) diff --git a/TODO b/TODO index 3a7c13e..96a8221 100644 --- a/TODO +++ b/TODO @@ -72,3 +72,10 @@ Reference values : --> 1000(nbTasks) tâches avec itérations possibles, puis phase 2 en fin de chaqune des 1000 tâches. On obtient 1000xK* médoïdes --> Phase 2 sur les 1000xK* médoïdes + +#point avec Jairo: +#rentrer dans code C cwt continue Rwave +#passer partie sowas à C +#fct qui pour deux series (ID, medoides) renvoie distance WER (Rwave ou à moi) +#transformee croisee , smoothing lissage 3 composantes , + calcul pour WER +#determiner nvoice noctave (entre octave + petit et + grand) diff --git a/epclust/R/main.R b/epclust/R/main.R index 867843b..e18ea7b 100644 --- a/epclust/R/main.R +++ b/epclust/R/main.R @@ -181,3 +181,4 @@ processChunk = function(indice, K, WER) #TODO: difficulté : retrouver courbe à partir de l'identifiant (DB ok mais le reste ?) #aussi : que passe-t-on aux noeuds ? curvesToCoeffs en // ? #enfin : WER ?! +#TODO: bout de code qui calcule les courbes synchrones après étapes 1+2 à partir des ID médoïdes diff --git a/epclust/R/stage2.R b/epclust/R/stage2.R index 9c15a74..fa55356 100644 --- a/epclust/R/stage2.R +++ b/epclust/R/stage2.R @@ -1,102 +1,39 @@ -#point avec Jairo: -#rentrer dans code C cwt continue Rwave -#passer partie sowas à C -#fct qui pour deux series (ID, medoides) renvoie distance WER (Rwave ou à moi) -#transformee croisee , smoothing lissage 3 composantes , + calcul pour WER -#attention : code fait pour des series temps desynchronisees ! (deltat, dt == 1,2 ...) -#determiner nvoice noctave (entre octave + petit et + grand) - library("Rwave") -#Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 -#TODO: bout de code qui calcule les courbes synchrones après étapes 1+2 à partir des ID médoïdes - -#toCWT: (aux) -##NOTE: renvoie une matrice 3D -toCWT <- function(X, sw=0, tw=0, swabs=0, nvoice=12, noctave=5, s0=2, w0=2*pi, - lt=24, dt=0.5, spectra=FALSE, smooth=TRUE, scaled=FALSE, scalevector) +#precondition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1 +toCWT <- function(X, tw=0, swabs=0, nvoice=12, noctave=5, s0=2, w0=2*pi, + spectra=FALSE, smooth=TRUE, scaled=FALSE, scalevector) { - noctave <- adjust.noctave(lt, dt, s0, tw, noctave) if(missing(scalevector)) scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 + s0log=as.integer((log2( s0*w0/(2*pi) )-1)*nvoice+1.5) + totnoct=noctave+as.integer(s0log/nvoice)+1 res <- lapply(1:nrow(X), function(n) { - tsX <- ts( X[n,] ) - tsCent <- tsX - mean(tsX) - if(scaled) - tsCent <- ts(scale(tsCent)) - tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0) - tsCent.cwt + ts <- scale(ts( X[n,] ), 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) + smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) + ts.cwt*smat }) if( spectra ) res <- lapply(res, function(l) Mod(l)^2 ) if( smooth ) - res <- lapply(res, smCWT, swabs = swabs, tw = tw, dt = dt, scalevector = scalevector) + res <- lapply(res, smCWT, swabs = swabs, tw = tw, scalevector = scalevector) resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]), length(res))) for( l in 1:length(res) ) resArray[ , , l] <- res[[l]] resArray } -#from sowas -adjust.noctave <- function(N,dt,s0,tw,noctave) -{ - if (tw>0) - { - dumno <- as.integer((log(N*dt)-log(2*tw*s0))/log(2)) - if (dumno<noctave) - { - cat("# noctave adjusted to time smoothing window \n") - noctave <- dumno - } - } - noctave -} - -#from sowas -cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi) +#smooth cwt result +smCWT <- function(CWT, tw= 0, swabs= 0, nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, + lt= 24, scalevector ) { - if (class(ts)!="ts") - stop("# This function needs a time series object as input. You may construct this by using the function ts(data,start,deltat). Try '?ts' for help.\n") - - t=time(ts) - dt=t[2]-t[1] - s0unit=s0/dt*w0/(2*pi) - s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) - if (s0log<1) - { - cat(paste("# s0unit = ",s0unit,"\n",sep="")) - cat(paste("# s0log = ",s0log,"\n",sep="")) - cat("# s0 too small for w0! \n") - } - totnoct=noctave+as.integer(s0log/nvoice)+1 - - #cwt from package Rwave - totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0) - ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)] - - #Normalization - sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) - smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) - - ts.cwt*smat -} - -#NOTE: vect2mat = as.matrix ?! (dans aux.R) -vect2mat <- function(vect, delta, lscvect) -{ - vect <- as.vector(vect) - matrix(vect[-(1:2)], delta, lscvect) -} - -#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 <- smooth.time(smwsp, tw, scalevector) smsmwsp } @@ -113,66 +50,57 @@ smooth.matrix <- function(wt,swabs) smwt } -smooth.time <- function(wt,tw,dt,scalevector) +smooth.time <- function(wt,tw,scalevector) { smwt <- wt if (tw != 0) { for (i in 1:length(scalevector)) { - twi <- as.integer(scalevector[i]*tw/dt) + twi <- as.integer(scalevector[i]*tw) 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) { - #(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] - # #NOTE: delta et lscvect pourraient etre gardés à part (communs) + m <- ncol(conso) #TODO: automatic tune of these parameters ? (for other users) nvoice <- 4 - # # noctave4 = 2^13 = 8192 half hours ~ 180 days - noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2, tw = 0, noctave = 13) - # # 4 here represent 2^5 = 32 half-hours ~ 1 day - scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2 + # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(conso)) + noctave = 13 + # 4 here represent 2^5 = 32 half-hours ~ 1 day + scalevector4 <- 2^(4:(noctave * nvoice) / nvoice) * 2 lscvect4 <- length(scalevector4) lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect # observations node with CWT - Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, scalevector = scalevector4, lt = delta, + Xcwt4 <- toCWT(conso, noctave = noctave, scalevector = scalevector4, smooth = FALSE, nvoice = nvoice) #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() + Xcwt2[i,] <- c(m, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) + + rm(conso, Xcwt4) ; gc() - #Benjamin: FIX is this OK ? lscvect = dim(Xcwt4)[2] - ## _.b WER^2 distances ######## Xwer_dist <- matrix(0.0, n, n) for(i in 1:(n - 1)) { -#browser() -##ERROR là sans FIX lscvect :: delta lscvect --> taille ??! - mat1 <- vect2mat(Xcwt2[i,], delta, lscvect) + mat1 <- matrix(as.vector(Xcwt2[i,])[-(1:2)], m, lscvect) for(j in (i + 1):n) { - mat2 <- vect2mat(Xcwt2[j,], delta, lscvect) + mat2 <- matrix(as.vector(Xcwt2[j,])[-(1:2)], m, lscvect) num <- Mod(mat1 * Conj(mat2)) WX <- Mod(mat1 * Conj(mat1)) WY <- Mod(mat2 * Conj(mat2)) @@ -180,8 +108,8 @@ step2 = function(conso) 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)) + sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) + Xwer_dist[i, j] <- sqrt(m * lscvect * (1 - wer2)) Xwer_dist[j, i] <- Xwer_dist[i, j] } } -- 2.44.0