code/script_clustering_by_pam.r

   1 # Original R script (not required by the C program)
   2
   3
   4 ##################################################################
   5 ## File: script_clustering_by_pam.r
   6 ##
   7 ## Description:   Using PAM to clustering conso using pam
   8 ## Last modified: jan 2012 by JC
   9 ##
  10 ##################################################################
  11
  12 ## Function: Converts a matrix of curves to the discret wavelet domain
  13 toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm"){
  14   x2   <- spline(x, n = 2^ceiling( log(length(x), 2) ),
  15                  method = 'natural')$y
  16   Dx2 <- wd(x2, family = family, filter.number = filter.number)$D
  17   return(Dx2)
  18 }
  19
  20 ## Function: Computes the absolute contribution of the wavelet's scale
  21 ##           to the total total of the curve.
  22 contrib <- function(x) {
  23   J   <- log( length(x)+1, 2)
  24   nrj <- numeric(J)
  25   t0  <- 1
  26   t1  <- 0
  27   for( j in 1:J ) {
  28     t1     <- t1 + 2^(J-j)
  29     nrj[j] <- sqrt( sum( x[t0:t1]^2 ) )
  30     t0     <- t1 + 1
  31   }
  32   return(nrj)
  33 }
  34
  35 ## 1. Load libraries & data ####
  36 library(wavethresh)
  37 library(cluster)
  38
  39
  40 # powerload is a matrix that contains on each line one observation
  41 #           of length delta
  42
  43 ## 2. DWT  ##################
  44 delta <- ncol(powerload)
  45 n     <- nrow(powerload)
  46 Xdwt  <- t(apply(powerload, 1, toDWT)) # DWT over the lines of powerload.
  47 Xnrj  <- t(apply(Xdwt, 1, contrib))    # Absolute contribution to the energy.
  48
  49 ## 3. Cluster  ##############
  50 K <- 8 # Number of clusters
  51 Xnrj_dist <- dist(Xnrj)
  52 Xnrj_pam  <- pam(as.dist(Xnrj_dist), K)