#options(error = recover)
#options(warn=2)
#library(readr)
#common
library(Rcpp)
library(mvtnorm)
library(distances)
library(plyr)
library(scclust)
library(combinat)
#library(Itis)
library(ggplot2)

cppFunction('NumericMatrix fastAgg(NumericMatrix orgMeans, IntegerVector cats) {
            
            //orgMeans is the original dataset
            //cats is the categories, numbered 0 to n-1
            
            //Store dimensons of the for loop
            long catLeng = cats.length();
            int numCols = orgMeans.ncol();
            
            //initialize number of categories, number of observations in each category, aggregated means
            long numCats = max(cats);
            IntegerVector catSize(numCats+1);
            NumericMatrix aggMeans(numCats+1,numCols);
            
            for(int j = 0; j < numCols; j++){
            for(long i = 0; i < catLeng; i++ ){
            //Different instructions if j = 0 and if j notequal zero
            if(j == 0){
            
            //Increase the category total
            catSize[cats[i]]++;
            
            //Update the means
            aggMeans(cats[i],j) = (double)(catSize[cats[i]]-1)/catSize[cats[i]]*aggMeans(cats[i],j) + (double)1/catSize[cats[i]]*orgMeans(i,j);
            
            }
            else{
            aggMeans(cats[i],j) = (double)aggMeans(cats[i],j) + (double)1/catSize[cats[i]]*orgMeans(i,j);
            }
            }
            }
            return aggMeans;
            
            }')

cppFunction('IntegerVector fastJoin2(IntegerVector mer1, IntegerVector mer2) {
            
            //mer1 is original data matrix
            //mer2 in final cluster
            //This function want to perform inner join
            
            //Store dimensons of the for loop
            
            long numN = mer1.length();
            IntegerVector mer3(numN);         
            
            for(long j = 0; j < numN; j++){
            mer3(j)=mer2( mer1(j) );
            }
            
            return mer3;
            }')

#original
hacforiginal<-function(){	
  matrix_dist<-dist(datamatrix,method="euclidean")
  fit <- hclust(matrix_dist, method="ward.D")
  dendogramgroups <- cutree(fit, k=3)
  return(dendogramgroups)
  #Assign a Value to a Name
  #assign("finalclustersoriginal",dendogramgroups,envir = globalenv())
}

Hac_Thres <- function(dat, t, n){
  my_dist <- distances(dat)
  my_clustering_new <- sc_clustering(my_dist, t)
  aggdata_old <- fastAgg(as.matrix(dat), my_clustering_new)
  my_clustering_old <- as.integer(my_clustering_new)
  clu_num <- max(my_clustering_old)
  clusterlabel <- my_clustering_old
  while(clu_num > t){
    dat <- aggdata_old
    my_dist <- distances(dat)
    my_clustering_new <- sc_clustering(my_dist, t)
    aggdata_old <- fastAgg(as.matrix(dat), my_clustering_new)
    label <- fastJoin2(my_clustering_old, my_clustering_new)
    clusterlabel <- cbind(clusterlabel, label)
    my_clustering_old <- label
    clu_num <- max(label)
  }
  return(clusterlabel)
}

#' Cluster sum of square (between, within, total)
ss <- function(x) sum(scale(x, scale = FALSE)^2)
#' Wss
#' Find the cluster that has smallest Within-cluster sum of square
#' @param dat is datamatrix
#' @param k is number of clusters for data
#' @param cluster is the clustering result
Wss <- function(dat, k, cluster){	
  within_ss <- vector(mode="numeric", length = k)
  for (i in 0:k-1){
    within_ss[i+1] <- ss(dat[which(cluster == i), ])
  }
  min_wss <- min(within_ss)
  min_lab <- which.min(within_ss) - 1
  min_center <- apply(dat[which(cluster == min_lab), ], 2, mean)
  #return(list(within_ss, min_wss, min_lab, min_center))
  return(append(min_center, min_lab))
}

#' Km_wss
#' Within-cluster sum of square
#' @param dat is datamatrix
#' @param k is number of clusters
#' @param cluster is the clustering result
Km_wss <- function(dat, k, cluster){	
  km_within <- vector(mode="numeric", length = k)
  for (i in 1:k){
    km_within[i] <- ss(dat[which(cluster == i), ])
  }
  km_wss <- sum(km_within)
  return(km_wss)
}

Assign_fct <- function(label_old, minimum_old, label_new){
  mer <- rep(NA, length(label_old))
  len <- length(label_new)
  mer[c(which(label_old == minimum_old[3]))] <- label_new[len]
  mer[c(-which(label_old == minimum_old[3]))] <- label_new[-len]
  return(mer)
}

##################Function for method II (new)
#' Function for Method II
#' @param dat is datamatrix
#' @param t is threshold value for TC
#' @param n is size for dataset
Fct_Hatc <- function(dat, t, n){
  finalcluster <- rep(NA, n)
  label_result <- list()
  min_result <- list()
  i=1
  my_dist <- distances(dat)
  label_result[[i]] <- sc_clustering(my_dist, t)
  clu_num <- max(label_result[[i]])+1
  finalcluster <- label_result[[i]]
  newdata <- dat
  while(nrow(newdata)>t){
    min_result[[i]] <- Wss(dat = newdata, k = clu_num, cluster = label_result[[i]])
    discard <- c(which(label_result[[i]] == min_result[[i]][3]))
    newdata <- rbind(newdata[-discard, ], min_result[[i]][-3])
    if(nrow(newdata) <= 2) break
    #dist_mat <- distance_matrix(my_dist)
    #dist_mat <- as.matrix(my_dist)
    #dist_mat_edit <- dist_mat[-discard, ]
    #dist_mat_edit <- dist_mat_edit[, -discard]
    #dist_mat_add <- apply(newdata[-nrow(newdata),], 1, function(x){ sqrt(sum((x-min_result[[i]][-3])^2))})
    #dist_mat_addrow <- rbind(cbind(dist_mat_edit, dist_mat_add), c(dist_mat_add,0))
    #my_dist <- as.dist(dist_mat_addrow)  
    my_dist <- distances(newdata)
    label_result[[i+1]] <- sc_clustering(my_dist, t)
    clu_num <- max(label_result[[i+1]])+1
    middle <- Assign_fct(label_result[[i]], min_result[[i]], label_result[[i+1]])
    if(i >1){
      for(j in i:2){
        mer4 <- Assign_fct(label_result[[j-1]], min_result[[j-1]], middle)
        middle <- mer4}
    }
    finalcluster <- rbind(finalcluster, middle)
    i <- i +1
  }
  return(finalcluster)
}

# Find combination of 3!
z <- permn(c(1:3))
len_accu <- length(z)
sims = 2:10
N<-10^3
t <-2
d <- 2
k <- 3

datamatrix = matrix(NA, nrow = N, ncol = 2)
colnames(datamatrix)<-c("x1","x2")
gid = rep(NA,N)
#for(simtime in sims){
ran =runif(N)

for(i in 1:N){
  if(ran[i]<.5){
    datamatrix[i,] = rmvnorm(1,mean=c(1,2),sigma=matrix(c(1,0,0,0.5),ncol=2,byrow=T))
    gid[i]=0
  }else if(ran[i]<.8){
    datamatrix[i,] = rmvnorm(1,mean=c(7,8),sigma=matrix(c(2,0,0,1),ncol=2,byrow=T))
    gid[i]=1
  }else{
    datamatrix[i,] = rmvnorm(1,mean=c(3,5),sigma=matrix(c(3,0,0,4),ncol=2,byrow=T))
    gid[i]=2
  }
}

output_matrix <- matrix( nrow = 7, ncol = 3, 
                         dimnames = list(c("memory", "time", "within", "total",
                                           "between", "bss/tss", "accuracy"), 
                                         c("HAC", "HTC","HATC")))

#original
Rprof ( tf <- "log10e4.log",  memory.profiling = TRUE )
hac_cluster <- hacforiginal()
Rprof ( NULL ) ; 
temp0<- summaryRprof ( tf ,memory="both")$by.total
output_matrix[1,1]<-as.numeric(temp0[1,]$mem.total)
output_matrix[2,1]<-as.numeric(temp0[1,]$total.time)
output_matrix[3,1] <- Km_wss(dat = datamatrix, k = k, cluster = hac_cluster)
output_matrix[4,1] <- ss(datamatrix)
output_matrix[5,1] <- output_matrix[4,1] - output_matrix[3,1]
output_matrix[6,1] <- output_matrix[5,1]/output_matrix[4,1]
result0<-table(gid,hac_cluster)
sumaccu <- vector(length = len_accu)
for(i in 1:len_accu){
  sumaccu[i] <- result0[1, z[[i]][1]] + result0[2, z[[i]][2]] + result0[3, z[[i]][3]]
}
output_matrix[7,1] <- max(sumaccu)/N

#HTC
input = 1
Rprof ( ktf <- "hyb10e4.log",  memory.profiling = TRUE )
finalcluster_threshold <- Hac_Thres(dat = datamatrix, t = t, n = N)
Rprof ( NULL ) ; 
temp<- summaryRprof ( ktf,memory = "both" )$by.total
tempnumber<-input + 1
output_matrix[1,tempnumber]<-as.numeric(temp[1,]$mem.total)
output_matrix[2,tempnumber]<-as.numeric(temp[1,]$total.time)
cluster_col <- ncol(finalcluster_threshold)
if(max(finalcluster_threshold[, cluster_col])+1 < k){
  result0 <- table(gid,finalcluster_threshold[, cluster_col-1])
}else{
  result0<-table(gid,finalcluster_threshold[, cluster_col])
}
print(result0)
if(ncol(result0) == k){
  sumaccu <- vector(length = len_accu)
  for(i in 1:len_accu){
    sumaccu[i] <- result0[1, z[[i]][1]] + result0[2, z[[i]][2]] + result0[3, z[[i]][3]]
  }
  output_matrix[7,tempnumber] <- max(sumaccu)/N
  output_matrix[3, tempnumber] <- Km_wss(dat = datamatrix, k = k, cluster = finalcluster_threshold[,7])
  output_matrix[4, tempnumber] <- ss(datamatrix)
  output_matrix[5, tempnumber] <- output_matrix[4, tempnumber] - output_matrix[3, tempnumber]
  output_matrix[6, tempnumber] <- output_matrix[5, tempnumber]/output_matrix[4, tempnumber]
}
sumaccu <- vector(length = length(z))
for(i in 1:length(z)){
  sumaccu[i] <- result0[1, z[[i]][1]] + result0[2, z[[i]][2]] + result0[3, z[[i]][3]]
}
output_matrix[7,tempnumber] <- max(sumaccu)/N

#HATC
input = 2
Rprof ( ktf <- "hyb10e4.log",  memory.profiling = TRUE )
finalcluster_hatc <- Fct_Hatc(dat = datamatrix, t = t, n = N)
Rprof ( NULL ) ; 
temp<- summaryRprof ( ktf,memory = "both" )$by.total
tempnumber<-input + 1
output_matrix[1,tempnumber]<-as.numeric(temp[1,]$mem.total)
output_matrix[2,tempnumber]<-as.numeric(temp[1,]$total.time)
count_clu <- rep(NA,nrow(finalcluster_hatc))
for(p in nrow(finalcluster_hatc):1){
  count_clu[i] <- max(finalcluster_hatc[i,])
  if(count_clu[i] == k){
    result0 <- table(gid,finalcluster_hatc[, i])
    break
  }
}
print(result0)
sumaccu <- vector(length = len_accu)
for(i in 1:len_accu){
  sumaccu[i] <- result0[1, z[[i]][1]] + result0[2, z[[i]][2]] + result0[3, z[[i]][3]]
}
output_matrix[7,tempnumber] <- max(sumaccu)/N
output_matrix[3, tempnumber] <- Km_wss(dat = datamatrix, k = k, cluster = finalcluster_hatc[,i])
output_matrix[4, tempnumber] <- ss(datamatrix)
output_matrix[5, tempnumber] <- output_matrix[4, tempnumber] - output_matrix[3, tempnumber]
output_matrix[6, tempnumber] <- output_matrix[5, tempnumber]/output_matrix[4, tempnumber]

simtime = 2
write.table(output_matrix, file = paste("hatc_10e3_",simtime,".csv",sep=""), 
            sep = ",", row.names = TRUE, col.names = TRUE)
#}
output_matrix