fix accent (TODO: fix gitweb...)
[ppam-mpi.git] / code / script_clustering_by_pam.r
CommitLineData
bf3db2b6
BA
1# Original R script (not required by the C program)
2
3
81923e5c
BA
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
13toDWT <- 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.
22contrib <- 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 ####
36library(wavethresh)
37library(cluster)
38
39
40# powerload is a matrix that contains on each line one observation
41# of length delta
42
43## 2. DWT ##################
44delta <- ncol(powerload)
45n <- nrow(powerload)
46Xdwt <- t(apply(powerload, 1, toDWT)) # DWT over the lines of powerload.
47Xnrj <- t(apply(Xdwt, 1, contrib)) # Absolute contribution to the energy.
48
49## 3. Cluster ##############
50K <- 8 # Number of clusters
51Xnrj_dist <- dist(Xnrj)
52Xnrj_pam <- pam(as.dist(Xnrj_dist), K)