#common
library(Rcpp)
library(mvtnorm)
library(distances)
library(plyr)
library(scclust)
library(combinat)

kmeansf<-function(t){	
  my_dist <- distances(datamatrix)
  my_clustering <- sc_clustering(my_dist, t)
  my_cluster<-as.integer(my_clustering)
  aggdata<- fastAgg(as.matrix(datamatrix), my_cluster)
  kfit <- kmeans(aggdata, 3, nstart=10)
  ft <-fastJoin2(my_cluster,kfit$cluster)
  assign("cluster_threshold",ft,envir = globalenv())
}

kmeansf_mul<-function(input,t){	
  my_dist <- distances(datamatrix)
  my_clustering_old <- sc_clustering(my_dist, t)
  my_cluster_old<- as.integer(my_clustering_old)
  aggdata_old<- fastAgg(as.matrix(datamatrix), my_cluster_old)
  j=1
  Sys.sleep(2)
  while(j<input){
    my_dist_new <- distances(aggdata_old)
    my_clustering_new <- sc_clustering(my_dist_new, t)
    my_cluster_new<-as.integer(my_clustering_new)
    fa<-fastJoin2(my_cluster_old,my_cluster_new)
    aggdata_old<- fastAgg(aggdata_old, my_cluster_new)
    my_cluster_old<-fa
    j <-j+1
  }
  kfit <- kmeans(aggdata_old, 3, nstart=10)
  ft <-fastJoin2(my_cluster_old,kfit$cluster)
  assign("finalcluster_threshold",ft,envir = globalenv())
}

#original km
kmeansforiginal<-function(){	
  fit <- kmeans(datamatrix, 3, nstart=10)
  #Assign a Value to a Name
  assign("finalclustersoriginal",fit$cluster,envir = globalenv())
}

khacf<-function(t){	
  my_dist <- distances(datamatrix)
  my_clustering <- sc_clustering(my_dist, t)
  my_cluster<-as.integer(my_clustering)
  aggdata<- fastAgg(as.matrix(datamatrix), my_cluster)
  kdist<- dist(aggdata, method = "euclidean")
  kfit <- hclust(kdist, method="ward.D") 
  kclust <- cutree(kfit, k=3)
  ft <-fastJoin2(my_cluster,kclust)
  assign("cluster_threshold",ft,envir = globalenv())
}

khacf_mul<-function(input,t){	
  my_dist <- distances(datamatrix)
  my_clustering_old <- sc_clustering(my_dist, t)
  my_cluster_old<- as.integer(my_clustering_old)
  aggdata_old<- fastAgg(as.matrix(datamatrix), my_cluster_old)
  j=1
  Sys.sleep(15)
  while(j<input){
    my_dist_new <- distances(aggdata_old)
    my_clustering_new <- sc_clustering(my_dist_new, t)
    my_cluster_new<-as.integer(my_clustering_new)
    fa<-fastJoin2(my_cluster_old,my_cluster_new)
    aggdata_old<- fastAgg(aggdata_old, my_cluster_new)
    my_cluster_old<-fa
    j <-j+1
  }
  kdist<- dist(aggdata_old, method = "euclidean")
  kfit <- hclust(kdist, method="ward.D") 
  kclust <- cutree(kfit, k=3)
  ft <-fastJoin2(my_cluster_old,kclust)
  assign("finalcluster_threshold",ft,envir = globalenv())
}

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

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;
            }')

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

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]="A"
    }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]="B"
    }else{
      datamatrix[i,] = rmvnorm(1,mean=c(3,5),sigma=matrix(c(3,0,0,4),ncol=2,byrow=T))
      gid[i]="C"
    }
  }
  
  output_matrix <- matrix( nrow = 3, ncol = 9,dimnames = list(c("memory", "time","accuracy"),c("m0","m1", "m2", "m3","m4","m5","m6", "m7","m8")))
  
  #original
  Rprof ( tf <- "log10e4.log",  memory.profiling = TRUE )
  kmeansforiginal()
  #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)
  result0<-table(gid,finalclustersoriginal)
  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[3,1] <- max(sumaccu)/N
  
  #hybrid
  Rprof ( ktf <- "hyb10e4.log",  memory.profiling = TRUE )
  kmeansf(t)
  #khacf(t)
  Rprof ( NULL ) ; 
  temp<- summaryRprof ( ktf,memory = "both" )$by.total
  output_matrix[1,2]<-as.numeric(temp[1,]$mem.total)
  output_matrix[2,2]<-as.numeric(temp[1,]$total.time)
  result0<-table(gid,cluster_threshold)
  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[3,2] <- max(sumaccu)/N
  
  for(input in 2:8) {
    Rprof ( ktf <- "hyb10e4.log",  memory.profiling = TRUE )
    kmeansf_mul(input,t)
    #khacf_mul(input,t)
    Rprof ( NULL ) ; 
    temp<- summaryRprof ( ktf,memory = "both" )$by.total
    print(temp)
    cat("\n")
    tempnumber<-input+1
    output_matrix[1,tempnumber]<-as.numeric(temp[1,]$mem.total)
    output_matrix[2,tempnumber]<-as.numeric(temp[1,]$total.time)
    resultf<-table(gid,finalcluster_threshold)
    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[3,tempnumber]<-max(sumaccu)/N
  }
  write.table(output_matrix, file = paste("fm2_km10e4_",simtime,".csv",sep=""), sep = ",", row.names = TRUE, col.names = TRUE)