| 1 | ## File : 05_cluster2step.r |
| 2 | ## Description : |
| 3 | |
| 4 | rm(list = ls()) |
| 5 | |
| 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') |
| 9 | library(cluster) |
| 10 | |
| 11 | ## 1. Read auxiliar data files #### |
| 12 | |
| 13 | identifiants <- read.table("identifs.txt")[ ,1] |
| 14 | dates0 <- read.table("datesall.txt")[, 1] |
| 15 | dates <- as.character(dates0[grep("2009", dates0)]) |
| 16 | rm(dates0) |
| 17 | |
| 18 | n <- length(identifiants) |
| 19 | p <- length(dates) |
| 20 | |
| 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]) |
| 25 | |
| 26 | #synchros09 <- t(synchros09) # series on rows, use if necessary |
| 27 | |
| 28 | ## 2. Compute contributions #### |
| 29 | auxDWT <- t(apply(synchros09, 1, toDWT)) |
| 30 | matcontrib <- t(apply(auxDWT, 1, contrib, rel = TRUE, logit = TRUE)) |
| 31 | rm(auxDWT) |
| 32 | |
| 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] |
| 36 | |
| 37 | hc <- hclust(dist(tdata[, ci$selectv]), method = "ward.D") |
| 38 | plot(hc) |
| 39 | |
| 40 | #clust <- cutree(hc, 2) |
| 41 | |
| 42 | for(K in 2:30){ |
| 43 | #K <- 3 |
| 44 | pamfit <- pam(tdata[-201, ci$selectv], k = K, stand = FALSE) |
| 45 | |
| 46 | #table(pamfit$clustering) |
| 47 | |
| 48 | SC <- matrix(0, ncol = p, nrow = K) |
| 49 | |
| 50 | clustfactor <- pamfit$clustering |
| 51 | # for(k in 1:K){ |
| 52 | # clustk <- which(clustfactor == k) |
| 53 | # if(length(clustk) > 0) { |
| 54 | # if(length(clustk) > 1) { |
| 55 | # SCk <- colSums(synchros09[which(clustfactor == k), ]) |
| 56 | # } else { |
| 57 | # SCk <- synchros09[which(clustfactor == k), ] |
| 58 | # } |
| 59 | # SC[k, ] <- SC[k, ] + SCk |
| 60 | # rm(SCk) |
| 61 | # } |
| 62 | #} |
| 63 | |
| 64 | write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt")) |
| 65 | #write.table(clustfactor, file = "~/tmp/clustfactor3.txt") |
| 66 | } |
| 67 | |
| 68 | # Plots |
| 69 | layout(1) |
| 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) |
| 72 | |
| 73 | |
| 74 | |