1 ## File : 05_cluster2step.r
6 setwd("~/Documents/projects/2014_EDF-Orsay-Lyon2/codes/")
7 source("~/Documents/projects/2014_EDF-Lyon2/R/01_StBr.r") # aux functs
8 source('http://eric.univ-lyon2.fr/~jcugliari/codes/functional-clustering.r')
11 ## 1. Read auxiliar data files ####
13 identifiants <- read.table("identifs.txt")[ ,1]
14 dates0 <- read.table("datesall.txt")[, 1]
15 dates <- as.character(dates0[grep("2009", dates0)])
18 n <- length(identifiants)
21 synchros09 <- as.matrix(read.table("~/tmp/2009_synchros200.txt"))
22 #synchros09 <- as.matrix(read.table("~/tmp/2009_synchros200RC.txt"))
23 nas <- which(is.na(synchros09)[1, ]) # some 1/1/2009 are missing
24 synchros09[nas, 1] <- colMeans(synchros09[2:4, nas])
26 #synchros09 <- t(synchros09) # series on rows, use if necessary
28 ## 2. Compute contributions ####
29 auxDWT <- t(apply(synchros09, 1, toDWT))
30 matcontrib <- t(apply(auxDWT, 1, contrib, rel = TRUE, logit = TRUE))
33 ## 3. Transform data & compute CI ####
34 ci <- CI(matcontrib[-201, ]) # last row has the total
35 tdata <- ci$tdata; rownames(tdata) <- rownames(matcontrib)[-201]
37 hc <- hclust(dist(tdata[, ci$selectv]), method = "ward.D")
40 #clust <- cutree(hc, 2)
44 pamfit <- pam(tdata[-201, ci$selectv], k = K, stand = FALSE)
46 #table(pamfit$clustering)
48 SC <- matrix(0, ncol = p, nrow = K)
50 clustfactor <- pamfit$clustering
52 # clustk <- which(clustfactor == k)
53 # if(length(clustk) > 0) {
54 # if(length(clustk) > 1) {
55 # SCk <- colSums(synchros09[which(clustfactor == k), ])
57 # SCk <- synchros09[which(clustfactor == k), ]
59 # SC[k, ] <- SC[k, ] + SCk
64 write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt"))
65 #write.table(clustfactor, file = "~/tmp/clustfactor3.txt")
70 matplot(t(SC)[48*10 + 1:(48*30), ], type = 'l', ylab = '',col = 1:3, lty = 1)
71 matplot(t(SC)[48*100 + 1:(48*30), ], type = 'l', ylab = '', col = 1:3, lty = 1)