old_C_code/stage2_UNFINISHED/src/unused/05_cluster2step.r

   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