add code/_ReadMe. Clean latex/ folder
[ppam-mpi.git] / code / script_clustering_by_pam.r
diff --git a/code/script_clustering_by_pam.r b/code/script_clustering_by_pam.r
new file mode 100644 (file)
index 0000000..b4b491f
--- /dev/null
@@ -0,0 +1,52 @@
+# 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)