--- /dev/null
+# Original R script (not required by the C program)
+
+
+##################################################################
+## File: script_clustering_by_pam.r
+##
+## Description: Using PAM to clustering conso using pam
+## Last modified: jan 2012 by JC
+##
+##################################################################
+
+## Function: Converts a matrix of curves to the discret wavelet domain
+toDWT <- function(x, filter.number = 6, family = "DaubLeAsymm"){
+ x2 <- spline(x, n = 2^ceiling( log(length(x), 2) ),
+ method = 'natural')$y
+ Dx2 <- wd(x2, family = family, filter.number = filter.number)$D
+ return(Dx2)
+}
+
+## Function: Computes the absolute contribution of the wavelet's scale
+## to the total total of the curve.
+contrib <- function(x) {
+ J <- log( length(x)+1, 2)
+ nrj <- numeric(J)
+ t0 <- 1
+ t1 <- 0
+ for( j in 1:J ) {
+ t1 <- t1 + 2^(J-j)
+ nrj[j] <- sqrt( sum( x[t0:t1]^2 ) )
+ t0 <- t1 + 1
+ }
+ return(nrj)
+}
+
+## 1. Load libraries & data ####
+library(wavethresh)
+library(cluster)
+
+
+# powerload is a matrix that contains on each line one observation
+# of length delta
+
+## 2. DWT ##################
+delta <- ncol(powerload)
+n <- nrow(powerload)
+Xdwt <- t(apply(powerload, 1, toDWT)) # DWT over the lines of powerload.
+Xnrj <- t(apply(Xdwt, 1, contrib)) # Absolute contribution to the energy.
+
+## 3. Cluster ##############
+K <- 8 # Number of clusters
+Xnrj_dist <- dist(Xnrj)
+Xnrj_pam <- pam(as.dist(Xnrj_dist), K)