[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


#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,
ad642dc6 BA	48	# tw = 0, noctave = 13)
ad642dc6 BA	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
	54
	55	#17000 colonnes coeff 1, puis 17000 coeff 2... [non : dans chaque tranche du cube]
	56
	57	#TODO: une fonction qui fait lignes 59 à 91
	58
	59	#cube:
ad642dc6 BA	60	# Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
	61	# scalevector = scalevector4,
	62	# lt = delta, smooth = FALSE,
	63	# nvoice = nvoice) # observations node with CWT
	64	#
572d139a BA	65	# #matrix:
572d139a BA	66	# ############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
ad642dc6 BA	67	# #Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
ad642dc6 BA	68	#
572d139a	69	# #NOTE: delta et lscvect pourraient etre gardés à part (communs)
ad642dc6 BA	70	# for(i in 1:n)
	71	# Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
	72	#
	73	# #rm(conso, Xcwt4); gc()
	74	#
	75	# ## _.b WER^2 distances ########
	76	# Xwer_dist <- matrix(0.0, n, n)
	77	# for(i in 1:(n - 1)){
	78	# mat1 <- vect2mat(Xcwt2[i,])
	79	# for(j in (i + 1):n){
	80	# mat2 <- vect2mat(Xcwt2[j,])
	81	# num <- Mod(mat1 * Conj(mat2))
	82	# WX <- Mod(mat1 * Conj(mat1))
	83	# WY <- Mod(mat2 * Conj(mat2))
	84	# smsmnum <- smCWT(num, scalevector = scalevector4)
	85	# smsmWX <- smCWT(WX, scalevector = scalevector4)
	86	# smsmWY <- smCWT(WY, scalevector = scalevector4)
	87	# wer2 <- sum(colSums(smsmnum)^2) /
	88	# sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
	89	# Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
	90	# Xwer_dist[j, i] <- Xwer_dist[i, j]
	91	# }
	92	# }
	93	# diag(Xwer_dist) <- numeric(n)
	94	#
	95	# save(Xwer_dist, file = "../res/2009_synchros200WER.Rdata")
	96	# save(Xwer_dist, file = "../res/2009_synchros200-randomWER.Rdata")
	97
47395e4f BA	98
	99
	100	#lignes 59 à 91 "dépliées" :
	101	Xcwt4 <- toCWT(conso, noctave = noctave4, dt = 1,
	102	scalevector = scalevector4,
	103	lt = delta, smooth = FALSE,
	104	nvoice = nvoice) # observations node with CWT
	105
	106	#matrix:
	107	############Xcwt2 <- matrix(0.0, nrow= n, ncol= 2 + delta * lscvect)
	108	Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length((c(Xcwt4[,,1]))))
	109
	110	#NOTE: delta et lscvect pourraient etre gardés à part (communs)
	111	for(i in 1:n)
	112	Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) )
	113
	114	#rm(conso, Xcwt4); gc()
	115
	116	## _.b WER^2 distances ########
	117	Xwer_dist <- matrix(0.0, n, n)
	118	for(i in 1:(n - 1)){
	119	mat1 <- vect2mat(Xcwt2[i,])
	120
	121	#NOTE: vect2mat = as.matrix ?! (dans aux.R)
	122	vect2mat <- function(vect){
	123	vect <- as.vector(vect)
	124	matrix(vect[-(1:2)], delta, lscvect)
	125	}
	126
	127	for(j in (i + 1):n){
	128	mat2 <- vect2mat(Xcwt2[j,])
	129	num <- Mod(mat1 * Conj(mat2))
	130	WX <- Mod(mat1 * Conj(mat1))
	131	WY <- Mod(mat2 * Conj(mat2))
	132	smsmnum <- smCWT(num, scalevector = scalevector4)
	133	smsmWX <- smCWT(WX, scalevector = scalevector4)
	134	smsmWY <- smCWT(WY, scalevector = scalevector4)
	135	wer2 <- sum(colSums(smsmnum)^2) /
	136	sum( sum(colSums(smsmWX) * colSums(smsmWY)) )
	137	Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2))
	138	Xwer_dist[j, i] <- Xwer_dist[i, j]
	139	}
	140	}
	141	diag(Xwer_dist) <- numeric(n)
	142
	143	#fonction smCWT (dans aux.R)
	144	smCWT <- function(CWT, sw= 0, tw= 0, swabs= 0,
	145	nvoice= 12, noctave= 2, s0= 2, w0= 2*pi,
	146	lt= 24, dt= 0.5, scalevector )
	147	{
	148	# noctave <- adjust.noctave(lt, dt, s0, tw, noctave)
	149	# scalevector <- 2^(0:(noctave * nvoice) / nvoice) * s0
	150	wsp <- Mod(CWT)
	151	smwsp <- smooth.matrix(wsp, swabs)
	152	smsmwsp <- smooth.time(smwsp, tw, dt, scalevector)
	153	smsmwsp
	154	}
	155
	156	#dans sowas.R
	157	smooth.matrix <- function(wt,swabs){
	158
	159	if (swabs != 0)
	160	smwt <- t(filter(t(wt),rep(1,2swabs+1)/(2swabs+1)))
	161	else
162	smwt <- wt
163
164	smwt
165
166	}
167	smooth.time <- function(wt,tw,dt,scalevector){
168
169	smwt <- wt
170
171	if (tw != 0){
172	for (i in 1:length(scalevector)){
173
174	twi <- as.integer(scalevector[i]*tw/dt)
175	smwt[,i] <- filter(wt[,i],rep(1,2twi+1)/(2twi+1))
176
177	}
178	}
179	smwt
180	}
181
182	#et filter() est dans stats::
183
184	#cf. filters en C dans : https://svn.r-project.org/R/trunk/src/library/stats/src/filter.c
185
186
187
ad642dc6 BA	188	load("../res/2009_synchros200WER.Rdata")
	189	#load("../res/2009_synchros200-randomWER.Rdata")
	190
	191	## 3. Cluster using WER distance matrix ####
	192
	193	#hc <- hclust(as.dist(Xwer_dist), method = "ward.D")
	194	#plot(hc)
	195	#
	196	# #clust <- cutree(hc, 2)
	197	#
	198	for(K in 2:30){
	199	#K <- 3
	200	#pamfit <- pam(tdata[-201, ci$selectv], k = K)
	201	pamfit <- pam(as.dist(Xwer_dist), k = K, diss = TRUE)
	202
	203	#table(pamfit$clustering)
	204
	205	SC <- matrix(0, ncol = p, nrow = K)
	206
	207	clustfactor <- pamfit$clustering
	208	# for(k in 1:K){
	209	# clustk <- which(clustfactor == k)
	210	# if(length(clustk) > 0) {
	211	# if(length(clustk) > 1) {
	212	# SCk <- colSums(synchros09[which(clustfactor == k), ])
	213	# } else {
	214	# SCk <- synchros09[which(clustfactor == k), ]
	215	# }
	216	# SC[k, ] <- SC[k, ] + SCk
	217	# rm(SCk)
	218	# }
	219	#}
	220
	221	#write.table(clustfactor, file = paste0("~/tmp/clustfactorRC", K, ".txt"))
	222	#write.table(clustfactor, file = "~/tmp/clustfactor3.txt")
	223	#write.table(clustfactor, file = paste0("~/tmp/clustfactorWER", K, ".txt"))
	224	write.table(clustfactor, file = paste0("~/tmp/clustfactor-randomWER", K, ".txt"))
	225	}
	226	#
	227	# # Plots
	228	# layout(1)
	229	# matplot(t(SC)[4810 + 1:(4830), ], type = 'l', ylab = '',col = 1:3, lty = 1)
	230	# matplot(t(SC)[48100 + 1:(4830), ], type = 'l', ylab = '', col = 1:3, lty = 1)
	231	#
	232	#
	233	#