1 ## File: extract-features.r
5 ## a. Load data & libraries ####
12 MOJARRITA <- Sys.info()[4] == "mojarrita"
15 setwd("~/Documents/projects/2014_EDF-Orsay-Lyon2/codes/")
17 setwd("~/2014_EDF-Orsay-Lyon2/codes/")
20 #source('http://eric.univ-lyon2.fr/~jcugliari/codes/functional-clustering.r')
23 matcontrib0 <- read.table(file = "~/tmp/2009_contrib.txt")
24 n <- nrow(matcontrib0)
26 sdcontrib <- apply(matcontrib0, 1, sd)
27 lims <- quantile(sdcontrib, probs = c(.005, .995)) # obtain 1%-extreme data
28 is_normal <- which((sdcontrib > lims[1]) & (sdcontrib < lims[2]))
30 matcontri_ext <- matcontrib0[-is_normal, ]
31 matcontrib <- matcontrib0[is_normal, ] # wipe out aberrant data
33 matcontrib <- t(apply(matcontrib, 1, function(x) x / sum(x)))
34 matcontrib <- t(apply(matcontrib, 1, function(p) log(p / (1 - p)) ))
37 ## b. Transform data & compute CI ####
39 tdata <- ci$tdata; rownames(tdata) <- rownames(matcontrib)
42 ## c. Clustering ##########
47 #setup parallel backend to use 8 processors
49 registerDoParallel(cl)
51 clfitlist <- foreach(icount(iters)) %dopar% {
54 clara(x = tdata[, selvar],
63 #save(clfit, file = 'clfit500.Rdata')
64 # save(clfit, file = 'clfit200RC.Rdata')
65 #save(clfitlist, file = 'clfitlist200.Rdata')
66 #rm(ci, matcontrib0, is_normal, lims, selvar)
70 res <- lapply(clfitlist, function(x) x$clustering)
73 save(data.frame(res), file = 'res/clfitdf200.Rdata')
76 ## d. Analyze results ##########
78 #1. Répartition du nombre d'observation par cluster
79 #plot(sort(table(clfit$clustering), decreasing = TRUE),
80 # type = 'l', ylab = 'Fréquence', xlab = 'Classe')
83 #clust <- res$clustering
84 # centres <- aggregate(conso, clust)
87 #sel_veille <- as.Date(rownames(conso)[sel - 1])
88 #sel_lendemain <- as.Date(rownames(conso)[sel + 1])
90 #res_clust <- data.frame(date = rownames(conso),
91 #veille = weekdays(sel_veille),
92 #lendemain = weekdays(sel_lendemain),
96 # assign(paste0("dates_clust", K),
97 # substr(subset(res_clust, clust == k)$date, 1, 7) )
102 #save(file = paste0(dtitle, "_clust.Rdata"),
103 #res_clust, selvar, K, gap)
106 #dates_clust1 <- substr(subset(dates, clust == 1)$date, 1, 7)