| 1 | library("Rwave") |
| 2 | |
| 3 | #Entrée : courbes synchrones, soit après étape 1 itérée, soit après chaqure étape 1 |
| 4 | #TODO: bout de code qui calcule les courbes synchrones après étapes 1+2 à partir des ID médoïdes |
| 5 | |
| 6 | #(Benjamin) |
| 7 | #à partir de là, "conso" == courbes synchrones |
| 8 | n <- nrow(conso) |
| 9 | delta <- ncol(conso) |
| 10 | |
| 11 | #17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube] |
| 12 | # #NOTE: delta et lscvect pourraient etre gardés à part (communs) |
| 13 | |
| 14 | #lignes 59 à 91 "dépliées" : |
| 15 | Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1, |
| 16 | scalevector = scalevector4, |
| 17 | lt = delta, smooth = FALSE, |
| 18 | nvoice = nvoice) # observations node with CWT |
| 19 | |
| 20 | #toCWT: (aux) |
| 21 | ##NOTE: renvoie une matrice 3D |
| 22 | toCWT <- function(X, sw= 0, tw= 0, swabs= 0, |
| 23 | nvoice= 12, noctave= 5, |
| 24 | s0= 2, w0= 2*pi, lt= 24, dt= 0.5, |
| 25 | spectra = FALSE, smooth = TRUE, |
| 26 | scaled = FALSE, |
| 27 | scalevector) |
| 28 | { noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
| 29 | if(missing(scalevector)) |
| 30 | scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 |
| 31 | res <- lapply(1:nrow(X), function(n) |
| 32 | { tsX <- ts( X[n,] ) |
| 33 | tsCent <- tsX - mean(tsX) |
| 34 | if(scaled) tsCent <- ts(scale(tsCent)) |
| 35 | tsCent.cwt <- cwt.ts(tsCent, s0, noctave, nvoice, w0) |
| 36 | tsCent.cwt |
| 37 | } ) |
| 38 | if( spectra ) res <- lapply(res, function(l) Mod(l)^2 ) |
| 39 | if( smooth ) res <- lapply(res, smCWT, swabs = swabs, |
| 40 | tw = tw, dt = dt, |
| 41 | scalevector = scalevector) |
| 42 | resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]), |
| 43 | length(res))) |
| 44 | for( l in 1:length(res) ) resArray[ , , l] <- res[[l]] |
| 45 | resArray |
| 46 | } |
| 47 | |
| 48 | #from sowas |
| 49 | cwt.ts <- function(ts,s0,noctave=5,nvoice=10,w0=2*pi){ |
| 50 | |
| 51 | if (class(ts)!="ts"){ |
| 52 | |
| 53 | cat("# 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") |
| 54 | |
| 55 | } |
| 56 | else{ |
| 57 | |
| 58 | t=time(ts) |
| 59 | dt=t[2]-t[1] |
| 60 | |
| 61 | s0unit=s0/dt*w0/(2*pi) |
| 62 | s0log=as.integer((log2(s0unit)-1)*nvoice+1.5) |
| 63 | |
| 64 | if (s0log<1){ |
| 65 | cat(paste("# s0unit = ",s0unit,"\n",sep="")) |
| 66 | cat(paste("# s0log = ",s0log,"\n",sep="")) |
| 67 | cat("# s0 too small for w0! \n") |
| 68 | } |
| 69 | totnoct=noctave+as.integer(s0log/nvoice)+1 |
| 70 | |
| 71 | #cwt from package Rwave |
| 72 | totts.cwt=cwt(ts,totnoct,nvoice,w0,plot=0) |
| 73 | |
| 74 | ts.cwt=totts.cwt[,s0log:(s0log+noctave*nvoice)] |
| 75 | |
| 76 | #Normalization |
| 77 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) |
| 78 | smat <- matrix(rep(sqs,length(t)),nrow=length(t),byrow=TRUE) |
| 79 | |
| 80 | ts.cwt*smat |
| 81 | |
| 82 | } |
| 83 | |
| 84 | } |
| 85 | |
| 86 | #matrix: |
| 87 | ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect) |
| 88 | Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1])))) |
| 89 | |
| 90 | #NOTE: delta et lscvect pourraient etre gardés à part (communs) |
| 91 | for(i in 1:n) |
| 92 | Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) |
| 93 | |
| 94 | #rm(conso, Xcwt4); gc() |
| 95 | |
| 96 | ## _.b WER^2 distances ######## |
| 97 | Xwer_dist <- matrix(0.0, n, n) |
| 98 | for(i in 1:(n - 1)){ |
| 99 | mat1 <- vect2mat(Xcwt2[i,]) |
| 100 | |
| 101 | #NOTE: vect2mat = as.matrix ?! (dans aux.R) |
| 102 | vect2mat <- function(vect){ |
| 103 | vect <- as.vector(vect) |
| 104 | matrix(vect[-(1:2)], delta, lscvect) |
| 105 | } |
| 106 | |
| 107 | for(j in (i + 1):n){ |
| 108 | mat2 <- vect2mat(Xcwt2[j,]) |
| 109 | num <- Mod(mat1 * Conj(mat2)) |
| 110 | WX <- Mod(mat1 * Conj(mat1)) |
| 111 | WY <- Mod(mat2 * Conj(mat2)) |
| 112 | smsmnum <- smCWT(num, scalevector = scalevector4) |
| 113 | smsmWX <- smCWT(WX, scalevector = scalevector4) |
| 114 | smsmWY <- smCWT(WY, scalevector = scalevector4) |
| 115 | wer2 <- sum(colSums(smsmnum)^2) / |
| 116 | sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) |
| 117 | Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) |
| 118 | Xwer_dist[j, i] <- Xwer_dist[i, j] |
| 119 | } |
| 120 | } |
| 121 | diag(Xwer_dist) <- numeric(n) |
| 122 | |
| 123 | #fonction smCWT (dans aux.R) |
| 124 | smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0, |
| 125 | nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, |
| 126 | lt= 24, dt= 0.5, scalevector ) |
| 127 | { |
| 128 | # noctave <- adjust.noctave(lt, dt, s0, tw, noctave) |
| 129 | # scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0 |
| 130 | wsp <- Mod(CWT) |
| 131 | smwsp <- smooth.matrix(wsp, swabs) |
| 132 | smsmwsp <- smooth.time(smwsp, tw, dt, scalevector) |
| 133 | smsmwsp |
| 134 | } |
| 135 | |
| 136 | #dans sowas.R (...donc on ne lisse pas à ce niveau ?) |
| 137 | smooth.matrix <- function(wt,swabs){ |
| 138 | |
| 139 | if (swabs != 0) |
| 140 | smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1))) |
| 141 | else |
| 142 | smwt <- wt |
| 143 | |
| 144 | smwt |
| 145 | |
| 146 | } |
| 147 | smooth.time <- function(wt,tw,dt,scalevector){ |
| 148 | |
| 149 | smwt <- wt |
| 150 | |
| 151 | if (tw != 0){ |
| 152 | for (i in 1:length(scalevector)){ |
| 153 | |
| 154 | twi <- as.integer(scalevector[i]*tw/dt) |
| 155 | smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1)) |
| 156 | |
| 157 | } |
| 158 | } |
| 159 | smwt |
| 160 | } |
| 161 | |
| 162 | #et filter() est dans stats:: |
| 163 | > filter |
| 164 | function (x, filter, method = c("convolution", "recursive"), |
| 165 | sides = 2L, circular = FALSE, init = NULL) |
| 166 | { |
| 167 | method <- match.arg(method) |
| 168 | x <- as.ts(x) |
| 169 | storage.mode(x) <- "double" |
| 170 | xtsp <- tsp(x) |
| 171 | n <- as.integer(NROW(x)) |
| 172 | if (is.na(n)) |
| 173 | stop("invalid value of nrow(x)", domain = NA) |
| 174 | nser <- NCOL(x) |
| 175 | filter <- as.double(filter) |
| 176 | nfilt <- as.integer(length(filter)) |
| 177 | if (is.na(n)) |
| 178 | stop("invalid value of length(filter)", domain = NA) |
| 179 | if (anyNA(filter)) |
| 180 | stop("missing values in 'filter'") |
| 181 | if (method == "convolution") { |
| 182 | if (nfilt > n) |
| 183 | stop("'filter' is longer than time series") |
| 184 | sides <- as.integer(sides) |
| 185 | if (is.na(sides) || (sides != 1L && sides != 2L)) |
| 186 | stop("argument 'sides' must be 1 or 2") |
| 187 | circular <- as.logical(circular) |
| 188 | if (is.na(circular)) |
| 189 | stop("'circular' must be logical and not NA") |
| 190 | if (is.matrix(x)) { |
| 191 | y <- matrix(NA, n, nser) |
| 192 | for (i in seq_len(nser)) y[, i] <- .Call(C_cfilter, |
| 193 | x[, i], filter, sides, circular) |
| 194 | } |
| 195 | else y <- .Call(C_cfilter, x, filter, sides, circular) |
| 196 | } |
| 197 | else { |
| 198 | if (missing(init)) { |
| 199 | init <- matrix(0, nfilt, nser) |
| 200 | } |
| 201 | else { |
| 202 | ni <- NROW(init) |
| 203 | if (ni != nfilt) |
| 204 | stop("length of 'init' must equal length of 'filter'") |
| 205 | if (NCOL(init) != 1L && NCOL(init) != nser) { |
| 206 | stop(sprintf(ngettext(nser, "'init' must have %d column", |
| 207 | "'init' must have 1 or %d columns", domain = "R-stats"), |
| 208 | nser), domain = NA) |
| 209 | } |
| 210 | if (!is.matrix(init)) |
| 211 | dim(init) <- c(nfilt, nser) |
| 212 | } |
| 213 | ind <- seq_len(nfilt) |
| 214 | if (is.matrix(x)) { |
| 215 | y <- matrix(NA, n, nser) |
| 216 | for (i in seq_len(nser)) y[, i] <- .Call(C_rfilter, |
| 217 | x[, i], filter, c(rev(init[, i]), double(n)))[-ind] |
| 218 | } |
| 219 | else y <- .Call(C_rfilter, x, filter, c(rev(init[, 1L]), |
| 220 | double(n)))[-ind] |
| 221 | } |
| 222 | tsp(y) <- xtsp |
| 223 | class(y) <- if (nser > 1L) |
| 224 | c("mts", "ts") |
| 225 | else "ts" |
| 226 | y |
| 227 | } |
| 228 | <bytecode: 0x1b05db8> |
| 229 | <environment: namespace:stats> |
| 230 | |
| 231 | |
| 232 | #cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c |
| 233 | #ifdef HAVE_CONFIG_H |
| 234 | # include <config.h> |
| 235 | #endif |
| 236 | |
| 237 | #include <R.h> |
| 238 | #include "ts.h" |
| 239 | |
| 240 | #ifndef min |
| 241 | #define min(a, b) ((a < b)?(a):(b)) |
| 242 | #define max(a, b) ((a < b)?(b):(a)) |
| 243 | #endif |
| 244 | |
| 245 | // currently ISNAN includes NAs |
| 246 | #define my_isok(x) (!ISNA(x) & !ISNAN(x)) |
| 247 | |
| 248 | #Pour method=="convolution" dans filter() (fonction R) |
| 249 | SEXP cfilter(SEXP sx, SEXP sfilter, SEXP ssides, SEXP scircular) |
| 250 | { |
| 251 | if (TYPEOF(sx) != REALSXP || TYPEOF(sfilter) != REALSXP) |
| 252 | error("invalid input"); |
| 253 | R_xlen_t nx = XLENGTH(sx), nf = XLENGTH(sfilter); |
| 254 | int sides = asInteger(ssides), circular = asLogical(scircular); |
| 255 | if(sides == NA_INTEGER || circular == NA_LOGICAL) error("invalid input"); |
| 256 | |
| 257 | SEXP ans = allocVector(REALSXP, nx); |
| 258 | |
| 259 | R_xlen_t i, j, nshift; |
| 260 | double z, tmp, *x = REAL(sx), *filter = REAL(sfilter), *out = REAL(ans); |
| 261 | |
| 262 | if(sides == 2) nshift = nf /2; else nshift = 0; |
| 263 | if(!circular) { |
| 264 | for(i = 0; i < nx; i++) { |
| 265 | z = 0; |
| 266 | if(i + nshift - (nf - 1) < 0 || i + nshift >= nx) { |
| 267 | out[i] = NA_REAL; |
| 268 | continue; |
| 269 | } |
| 270 | for(j = max(0, nshift + i - nx); j < min(nf, i + nshift + 1) ; j++) { |
| 271 | tmp = x[i + nshift - j]; |
| 272 | if(my_isok(tmp)) z += filter[j] * tmp; |
| 273 | else { out[i] = NA_REAL; goto bad; } |
| 274 | } |
| 275 | out[i] = z; |
| 276 | bad: |
| 277 | continue; |
| 278 | } |
| 279 | } else { /* circular */ |
| 280 | for(i = 0; i < nx; i++) |
| 281 | { |
| 282 | z = 0; |
| 283 | for(j = 0; j < nf; j++) { |
| 284 | R_xlen_t ii = i + nshift - j; |
| 285 | if(ii < 0) ii += nx; |
| 286 | if(ii >= nx) ii -= nx; |
| 287 | tmp = x[ii]; |
| 288 | if(my_isok(tmp)) z += filter[j] * tmp; |
| 289 | else { out[i] = NA_REAL; goto bad2; } |
| 290 | } |
| 291 | out[i] = z; |
| 292 | bad2: |
| 293 | continue; |
| 294 | } |
| 295 | } |
| 296 | return ans; |
| 297 | } |