codes_Jairo_stage2/aux.r

   1 ####################################################################
   2 ##
   3 ## File:        aux.r
   4 ##
   5 ## Description: Miscelaneous functions for clustering with kcca
   6 ##
   7 ## Modified:    june 2010
   8 ##
   9 ####################################################################
  10
  11
  12   #######################################################
  13
  14   # Transforms a matrix of data (one observation by row)
  15   #     into an array where position[ , , i] gives
  16   #     the smoothed modulus of the i-th cwt observation
  17
  18   ########################################################
  19
  20
  21   toCWT  <- function(X, sw=  0,  tw=  0, swabs= 0,
  22                        nvoice= 12, noctave= 5,
  23                        s0= 2, w0= 2*pi, lt= 24, dt= 0.5,
  24                        spectra = FALSE, smooth = TRUE,
  25                        scaled  = FALSE,
  26                      scalevector)
  27      { noctave  <- adjust.noctave(lt, dt, s0, tw, noctave)
  28        if(missing(scalevector))
  29           scalevector  <- 2^(0:(noctave * nvoice) / nvoice) * s0
  30        res <- lapply(1:nrow(X), function(n)
  31            { tsX         <- ts( X[n,] )
  32              tsCent      <- tsX - mean(tsX)
  33              if(scaled)  tsCent <- ts(scale(tsCent))
  34              tsCent.cwt  <- cwt.ts(tsCent, s0, noctave, nvoice, w0)
  35              tsCent.cwt
  36            } )
  37        if( spectra ) res <- lapply(res, function(l) Mod(l)^2 )
  38        if( smooth  ) res <- lapply(res, smCWT, swabs = swabs,
  39                                    tw = tw, dt = dt,
  40                                    scalevector = scalevector)
  41        resArray <- array(NA, c(nrow(res[[1]]), ncol(res[[1]]),
  42                                length(res)))
  43        for( l in 1:length(res) ) resArray[ , , l] <- res[[l]]
  44        resArray
  45      }
  46
  47
  48   # ===============================================================
  49
  50   smCWT <- function(CWT, sw=  0,  tw=  0, swabs= 0,
  51                        nvoice= 12, noctave= 2, s0= 2, w0= 2*pi,
  52                        lt= 24, dt= 0.5, scalevector )
  53          {
  54 #         noctave  <- adjust.noctave(lt, dt, s0, tw, noctave)
  55 #         scalevector  <- 2^(0:(noctave * nvoice) / nvoice) * s0
  56          wsp     <- Mod(CWT)
  57          smwsp   <- smooth.matrix(wsp, swabs)
  58          smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
  59          smsmwsp
  60        }
  61
  62
  63   # ===============================================================
  64
  65   toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm")
  66 { x2   <- spline(x, n = 2^ceiling( log(length(x), 2) ),
  67             method = 'natural')$y
  68   Dx2 <- wd(x2, family = family, filter.number = filter.number)$D
  69         Dx2
  70 }
  71
  72   # ===============================================================
  73
  74   contrib <- function(x)
  75       { J   <- log( length(x)+1, 2)
  76         nrj <- numeric(J)
  77         t0  <- 1
  78         t1  <- 0
  79         for( j in 1:J ) {
  80           t1     <- t1 + 2^(J-j)
  81           nrj[j] <- sqrt( sum( x[t0:t1]^2 ) )
  82           t0     <- t1 + 1
  83         }
  84         return(nrj)
  85       }
  86
  87
  88     # ========================================= distance for coh ===
  89
  90   coherence <- function( x, y)
  91        { J <- log(length(x) + 1, 2)
  92          t0 <- 1
  93          sg2_x <- 0
  94          sg2_y <- 0
  95          sg_xy <- 0
  96          for(j in 0:(J - 1))
  97          {  t1 <- t0 + 2^(J - j)/2  - 1
  98             tt <- t0:t1
  99             sg2_x <- sg2_x + mean(x[t0:t1]^2)
 100             sg2_y <- sg2_y + mean(y[t0:t1]^2)
 101             sg_xy <- sg_xy + mean(x[t0:t1] * y[t0:t1])
 102             t0 <- t1 + 1
 103          }
 104         res <- sg_xy^2 / sg2_x / sg2_y
 105         res
 106        }
 107
 108
 109   vect2mat <- function(vect){
 110                  vect <- as.vector(vect)
 111                  matrix(vect[-(1:2)], delta, lscvect)
 112                }
 113
 114
 115   # =========================================  # myimg for graphics
 116   jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF",
 117                                 "cyan", "#7FFF7F", "yellow",
 118                                 "#FF7F00", "red", "#7F0000"))
 119
 120   myimg <- function(MAT, x = 1:nrow(MAT), y = 1:col(MAT), ... )
 121               filled.contour(  x = x, y = y, z = MAT,
 122                    xlab= 'Time', ylab= 'scale',
 123                    color.palette = jet.colors,
 124                    ... )
 125
 126