[epclust.git] / old_C_code / stage2_UNFINISHED / src / 05_cluster2stepWER.r

## File : 05_cluster2stepWER.r
## Description : 

rm(list = ls())

setwd("~/ownCloud/projects/2014_EDF-Orsay-Lyon2/codes/")
  
library(Rwave)       # CWT
library(cluster)     # pam
#library(flexclust)   # kcca
source("aux.r")               # auxiliary clustering functions 

#TODO: [plus tard] alternative à sowa (package disparu) : cwt..
source("sowas-superseded.r")  # auxiliary CWT functions

## 1. Read auxiliar data files ####

identifiants <- read.table("identifs.txt")[ ,1]
dates0       <- read.table("datesall.txt")[, 1]
dates        <- as.character(dates0[grep("2009", dates0)])
rm(dates0)

n     <- length(identifiants)
p <- delta <- length(dates)

synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200RC.txt")))
#synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200-random.txt")))

nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing
synchros09[nas, 1] <- rowMeans(synchros09[nas, 2:4]) #valeurs après 1er janvier

#moyenne pondérée pour compléter deux demi-heures manquantes
imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518])

conso <- synchros09[-201, ];  # series must be on rows 
n     <- nrow(conso)
delta <- ncol(conso)

rm(synchros09, nas)

## 2. Compute WER distance matrix ####

## _.a CWT -- Filtering the lowest freqs (>6m) ####
nvoice   <- 4
# # noctave4 = 2^13 = 8192 half hours ~ 180 days
noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2,
                           tw = 0, noctave = 13)
# # 4 here represent 2^5 = 32 half-hours ~ 1 day
scalevector4  <- 2^(4:(noctave4 * nvoice) / nvoice) * 2
lscvect4      <- length(scalevector4)
lscvect <- lscvect4  # i should clean my code: werFam demands a lscvect


#(Benjamin)
#à partir de là, "conso" == courbes synchrones
n     <- nrow(conso)
delta <- ncol(conso)


#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]

#TODO: une fonction qui fait lignes 59 à 91

#cube:
# Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
#                 scalevector = scalevector4,
#                 lt = delta, smooth = FALSE, 
#                 nvoice = nvoice)      # observations node with CWT
# 
# #matrix:
# ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
# #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
# 
# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
# for(i in 1:n) 
#  Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) 
# 
# #rm(conso, Xcwt4); gc()
# 
# ## _.b WER^2 distances  ########
# Xwer_dist    <- matrix(0.0, n, n)
# for(i in 1:(n - 1)){
#  mat1   <- vect2mat(Xcwt2[i,])
#  for(j in (i + 1):n){
#     mat2 <- vect2mat(Xcwt2[j,])
#     num     <- Mod(mat1 * Conj(mat2))
#     WX      <- Mod(mat1 * Conj(mat1))
#     WY      <- Mod(mat2 * Conj(mat2))
#     smsmnum <- smCWT(num, scalevector = scalevector4)
#     smsmWX  <- smCWT(WX,  scalevector = scalevector4)
#     smsmWY  <- smCWT(WY,  scalevector = scalevector4)
#     wer2    <- sum(colSums(smsmnum)^2)  /
#       sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
#     Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
#     Xwer_dist[j, i] <- Xwer_dist[i, j]
#   }
# }
# diag(Xwer_dist) <- numeric(n)
# 
# save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
# save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")


#lignes 59 à 91 "dépliées" :
Xcwt4   <- toCWT(conso, noctave = noctave4, dt = 1,
                 scalevector = scalevector4,
                 lt = delta, smooth = FALSE, 
                 nvoice = nvoice)      # observations node with CWT
 
 #matrix:
 ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
 Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
 
 #NOTE: delta et lscvect pourraient etre gardés à part (communs)
 for(i in 1:n) 
    Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) 
 
 #rm(conso, Xcwt4); gc()
 
 ## _.b WER^2 distances  ########
 Xwer_dist    <- matrix(0.0, n, n)
 for(i in 1:(n - 1)){
  mat1   <- vect2mat(Xcwt2[i,])

 #NOTE: vect2mat = as.matrix ?! (dans aux.R)
  vect2mat <- function(vect){
                 vect <- as.vector(vect)
                 matrix(vect[-(1:2)], delta, lscvect)
               }
 
 for(j in (i + 1):n){
     mat2 <- vect2mat(Xcwt2[j,])
     num     <- Mod(mat1 * Conj(mat2))
     WX      <- Mod(mat1 * Conj(mat1))
     WY      <- Mod(mat2 * Conj(mat2))
     smsmnum <- smCWT(num, scalevector = scalevector4)
     smsmWX  <- smCWT(WX,  scalevector = scalevector4)
     smsmWY  <- smCWT(WY,  scalevector = scalevector4)
     wer2    <- sum(colSums(smsmnum)^2)  /
       sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
     Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
     Xwer_dist[j, i] <- Xwer_dist[i, j]
   }
 }
 diag(Xwer_dist) <- numeric(n)

#fonction smCWT (dans aux.R)
  smCWT <- function(CWT, sw=  0,  tw=  0, swabs= 0,
                       nvoice= 12, noctave= 2, s0= 2, w0= 2*pi, 
					   lt= 24, dt= 0.5, scalevector )
		 {
#         noctave  <- adjust.noctave(lt, dt, s0, tw, noctave)
#         scalevector  <- 2^(0:(noctave * nvoice) / nvoice) * s0
         wsp     <- Mod(CWT)  
         smwsp   <- smooth.matrix(wsp, swabs)
         smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
         smsmwsp
       }

 #dans sowas.R
smooth.matrix <- function(wt,swabs){
  
  if (swabs != 0)
    smwt <- t(filter(t(wt),rep(1,2*swabs+1)/(2*swabs+1)))
  else
    smwt <- wt
  
  smwt
  
}
smooth.time <- function(wt,tw,dt,scalevector){
  
  smwt <- wt
  
  if (tw != 0){
    for (i in 1:length(scalevector)){
      
      twi <- as.integer(scalevector[i]*tw/dt)
      smwt[,i] <- filter(wt[,i],rep(1,2*twi+1)/(2*twi+1))
      
    }
  } 
  smwt
}

#et filter() est dans stats::

#cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c


load("../res/2009_synchros200WER.Rdata")
#load("../res/2009_synchros200-randomWER.Rdata")

## 3. Cluster using WER distance matrix ####
 
#hc    <- hclust(as.dist(Xwer_dist), method = "ward.D")
#plot(hc)
# 
# #clust <- cutree(hc, 2)
# 
for(K in 2:30){ 
  #K <- 3
  #pamfit <- pam(tdata[-201, ci$selectv], k = K)
  pamfit <- pam(as.dist(Xwer_dist), k = K, diss = TRUE)
  
  #table(pamfit$clustering)

  SC <- matrix(0, ncol = p, nrow = K)

  clustfactor <- pamfit$clustering
#  for(k in 1:K){ 
#    clustk <- which(clustfactor == k)
#  if(length(clustk) > 0) {
#    if(length(clustk) > 1) {
#      SCk <- colSums(synchros09[which(clustfactor == k), ])
#      } else {
#        SCk <- synchros09[which(clustfactor == k), ]
#      }
#      SC[k, ] <- SC[k, ] + SCk
#      rm(SCk)
#  }
#}

#write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt"))
#write.table(clustfactor, file = "~/tmp/clustfactor3.txt")
#write.table(clustfactor, file = paste0("~/tmp/clustfactorWER", K, ".txt"))
write.table(clustfactor, file = paste0("~/tmp/clustfactor-randomWER", K, ".txt"))
}
# 
# # Plots
# layout(1)
# matplot(t(SC)[48*10 + 1:(48*30), ],  type = 'l', ylab = '',col = 1:3, lty = 1)
# matplot(t(SC)[48*100 + 1:(48*30), ], type = 'l', ylab = '', col = 1:3, lty = 1)
# 
# 
#
Commit	Line	Data
ad642dc6 BA	1	## File : 05_cluster2stepWER.r
	2	## Description :
	3
	4	rm(list = ls())
	5
	6	setwd("~/ownCloud/projects/2014_EDF-Orsay-Lyon2/codes/")
	7
	8	library(Rwave) # CWT
	9	library(cluster) # pam
	10	#library(flexclust) # kcca
	11	source("aux.r") # auxiliary clustering functions
572d139a BA	12
572d139a BA	13	#TODO: [plus tard] alternative à sowa (package disparu) : cwt..
ad642dc6 BA	14	source("sowas-superseded.r") # auxiliary CWT functions
	15
	16	## 1. Read auxiliar data files ####
	17
	18	identifiants <- read.table("identifs.txt")[ ,1]
	19	dates0 <- read.table("datesall.txt")[, 1]
	20	dates <- as.character(dates0[grep("2009", dates0)])
	21	rm(dates0)
	22
	23	n <- length(identifiants)
	24	p <- delta <- length(dates)
	25
	26	synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200RC.txt")))
	27	#synchros09 <- t(as.matrix(read.table("~/tmp/2009_synchros200-random.txt")))
	28
	29	nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing
572d139a	30	synchros09[nas, 1] <- rowMeans(synchros09[nas, 2:4]) #valeurs après 1er janvier
ad642dc6	31
572d139a	32	#moyenne pondérée pour compléter deux demi-heures manquantes
ad642dc6 BA	33	imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
	34	synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518])
	35
	36	conso <- synchros09[-201, ]; # series must be on rows
	37	n <- nrow(conso)
	38	delta <- ncol(conso)
	39
	40	rm(synchros09, nas)
	41
	42	## 2. Compute WER distance matrix ####
	43
	44	## _.a CWT -- Filtering the lowest freqs (>6m) ####
47395e4f	45	nvoice <- 4
ad642dc6	46	# # noctave4 = 2^13 = 8192 half hours ~ 180 days
47395e4f	47	noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2,
dc1aa85a	48	tw = 0, noctave = 13)
ad642dc6	49	# # 4 here represent 2^5 = 32 half-hours ~ 1 day
47395e4f BA	50	scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2
	51	lscvect4 <- length(scalevector4)
	52	lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect
572d139a BA	53
572d139a BA	54
dc1aa85a BA	55	#(Benjamin)
	56	#à partir de là, "conso" == courbes synchrones
	57	n <- nrow(conso)
	58	delta <- ncol(conso)
	59
	60
572d139a BA	61	#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]
	62
	63	#TODO: une fonction qui fait lignes 59 à 91
	64
	65	#cube:
ad642dc6 BA	66	# Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
	67	# scalevector = scalevector4,
	68	# lt = delta, smooth = FALSE,
	69	# nvoice = nvoice) # observations node with CWT
	70	#
572d139a BA	71	# #matrix:
572d139a BA	72	# ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
ad642dc6 BA	73	# #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
ad642dc6 BA	74	#
572d139a	75	# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
ad642dc6 BA	76	# for(i in 1:n)
	77	# Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
	78	#
	79	# #rm(conso, Xcwt4); gc()
	80	#
	81	# ## _.b WER^2 distances ########
	82	# Xwer_dist <- matrix(0.0, n, n)
	83	# for(i in 1:(n - 1)){
	84	# mat1 <- vect2mat(Xcwt2[i,])
	85	# for(j in (i + 1):n){
	86	# mat2 <- vect2mat(Xcwt2[j,])
	87	# num <- Mod(mat1 * Conj(mat2))
	88	# WX <- Mod(mat1 * Conj(mat1))
	89	# WY <- Mod(mat2 * Conj(mat2))
	90	# smsmnum <- smCWT(num, scalevector = scalevector4)
	91	# smsmWX <- smCWT(WX, scalevector = scalevector4)
	92	# smsmWY <- smCWT(WY, scalevector = scalevector4)
	93	# wer2 <- sum(colSums(smsmnum)^2) /
	94	# sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
	95	# Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
	96	# Xwer_dist[j, i] <- Xwer_dist[i, j]
	97	# }
	98	# }
	99	# diag(Xwer_dist) <- numeric(n)
	100	#
	101	# save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
	102	# save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")
	103
47395e4f BA	104
	105
	106	#lignes 59 à 91 "dépliées" :
	107	Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
	108	scalevector = scalevector4,
	109	lt = delta, smooth = FALSE,
	110	nvoice = nvoice) # observations node with CWT
	111
	112	#matrix:
	113	############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
	114	Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
	115
	116	#NOTE: delta et lscvect pourraient etre gardés à part (communs)
	117	for(i in 1:n)
	118	Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
	119
	120	#rm(conso, Xcwt4); gc()
	121
	122	## _.b WER^2 distances ########
	123	Xwer_dist <- matrix(0.0, n, n)
	124	for(i in 1:(n - 1)){
	125	mat1 <- vect2mat(Xcwt2[i,])
	126
	127	#NOTE: vect2mat = as.matrix ?! (dans aux.R)
	128	vect2mat <- function(vect){
	129	vect <- as.vector(vect)
	130	matrix(vect[-(1:2)], delta, lscvect)
	131	}
	132
	133	for(j in (i + 1):n){
	134	mat2 <- vect2mat(Xcwt2[j,])
	135	num <- Mod(mat1 * Conj(mat2))
	136	WX <- Mod(mat1 * Conj(mat1))
	137	WY <- Mod(mat2 * Conj(mat2))
	138	smsmnum <- smCWT(num, scalevector = scalevector4)
	139	smsmWX <- smCWT(WX, scalevector = scalevector4)
	140	smsmWY <- smCWT(WY, scalevector = scalevector4)
	141	wer2 <- sum(colSums(smsmnum)^2) /
	142	sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
	143	Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
	144	Xwer_dist[j, i] <- Xwer_dist[i, j]
	145	}
	146	}
	147	diag(Xwer_dist) <- numeric(n)
	148
	149	#fonction smCWT (dans aux.R)
	150	smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0,
	151	nvoice= 12, noctave= 2, s0= 2, w0= 2*pi,
	152	lt= 24, dt= 0.5, scalevector )
	153	{
	154	# noctave <- adjust.noctave(lt, dt, s0, tw, noctave)
	155	# scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
	156	wsp <- Mod(CWT)
	157	smwsp <- smooth.matrix(wsp, swabs)
	158	smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
	159	smsmwsp
	160	}
	161
	162	#dans sowas.R
	163	smooth.matrix <- function(wt,swabs){
	164
	165	if (swabs != 0)
	166	smwt <- t(filter(t(wt),rep(1,2swabs+1)/(2swabs+1)))
	167	else
168	smwt <- wt
169
170	smwt
171
172	}
173	smooth.time <- function(wt,tw,dt,scalevector){
174
175	smwt <- wt
176
177	if (tw != 0){
178	for (i in 1:length(scalevector)){
179
180	twi <- as.integer(scalevector[i]*tw/dt)
181	smwt[,i] <- filter(wt[,i],rep(1,2twi+1)/(2twi+1))
182
183	}
184	}
185	smwt
186	}
187
188	#et filter() est dans stats::
189
190	#cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c
191
192
193
ad642dc6 BA	194	load("../res/2009_synchros200WER.Rdata")
	195	#load("../res/2009_synchros200-randomWER.Rdata")
	196
	197	## 3. Cluster using WER distance matrix ####
	198
	199	#hc <- hclust(as.dist(Xwer_dist), method = "ward.D")
	200	#plot(hc)
	201	#
	202	# #clust <- cutree(hc, 2)
	203	#
	204	for(K in 2:30){
	205	#K <- 3
	206	#pamfit <- pam(tdata[-201, ci$selectv], k = K)
	207	pamfit <- pam(as.dist(Xwer_dist), k = K, diss = TRUE)
	208
	209	#table(pamfit$clustering)
	210
	211	SC <- matrix(0, ncol = p, nrow = K)
	212
	213	clustfactor <- pamfit$clustering
	214	# for(k in 1:K){
	215	# clustk <- which(clustfactor == k)
	216	# if(length(clustk) > 0) {
	217	# if(length(clustk) > 1) {
	218	# SCk <- colSums(synchros09[which(clustfactor == k), ])
	219	# } else {
	220	# SCk <- synchros09[which(clustfactor == k), ]
	221	# }
	222	# SC[k, ] <- SC[k, ] + SCk
	223	# rm(SCk)
	224	# }
	225	#}
	226
	227	#write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt"))
	228	#write.table(clustfactor, file = "~/tmp/clustfactor3.txt")
	229	#write.table(clustfactor, file = paste0("~/tmp/clustfactorWER", K, ".txt"))
	230	write.table(clustfactor, file = paste0("~/tmp/clustfactor-randomWER", K, ".txt"))
	231	}
	232	#
	233	# # Plots
	234	# layout(1)
	235	# matplot(t(SC)[4810 + 1:(4830), ], type = 'l', ylab = '',col = 1:3, lty = 1)
	236	# matplot(t(SC)[48100 + 1:(4830), ], type = 'l', ylab = '', col = 1:3, lty = 1)
	237	#
	238	#
	239	#