epclust/R/stage2.R

   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 }