Commit | Line | Data |
---|---|---|
ad642dc6 BA |
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 |