b4b491fcc1b97bf0e9b422f883c09b598b2765e1
[ppam-mpi.git] / 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)