# File containing functions used to estimate cepstral coefficients
source("Estim_cep_coef.R")
# File containing functions used to estimate similarity measure
source("Estim_similarity.R")

# Number of Simulations
NSIM <- 1000

# Order of the model
p<-c(1:5)

# Degree of Smoothing
nspans <- c(0,3,5,7)

# Number of time series
G<- 31

# Time series length
N<-12

# Cluster Sizes (They should sum up to G-1)
cluster.sizes <- c(15,15)

# Stop if cluster sizes do not sum up to G-1
stopifnot(sum(cluster.sizes)==G-1)

# Vector of clustering methods to be used
# Available methods are
# "hclust.average", "hclust.single","hclust.complete","hclust.ward","diana","pam","kmeans"
clust.method<- c("hclust.average", "hclust.single","hclust.complete","hclust.ward","diana","pam","kmeans")
nclust.methods <- length(clust.method)

# Use the following to get the average of the simulations
# apply(Similarities.array,c(1:3,5),mean)
# apply(Kmeans.Similarities,1:2,mean)

# Function that generates time series sample for each simulation
get.sample <- function(G,N)
{
    Y<-matrix(NA, nrow=G, ncol=N)

    for(i in 1:((G-1)/2))
    {
      Y[i,] <- rnorm(N)
    }
    
    for(i in 1:((G-1)/2+1))
    {
      omega <- runif(1,0.97,1.03)
      Y[i+((G-1)/2),] <- cos(omega*(2*N*(i-1)+ c(1:N)))+sin(omega*(2*N*(i-1)+ c(1:N)))+rnorm(N,0,2)
    }
                              
    return(Y)
}


#############################################
#      Do not modify after this point!      #
#############################################



do_kmeans <-FALSE

if(sum(clust.method=="kmeans")>0)
{
  # Array used to store the similarity measure for kmeans
  Kmeans.Similarities <- array(NA,dim=c(length(p),length(nspans),NSIM))
  dimnames(Kmeans.Similarities) <- list(p,nspans,c(1:NSIM))
  clust.method <- clust.method[order(clust.method=="kmeans")]
  clust.method <- clust.method[-nclust.methods]
  nclust.methods <- nclust.methods-1
  do_kmeans <-TRUE
}

if(nclust.methods>0)
{
  # Array used to store the similarity measure for all clustering methods except kmeans
  Similarities.array <- array(NA,dim=c(length(p),length(nspans),5,NSIM,nclust.methods))
  dimnames(Similarities.array) <- list(p,nspans,c("Euclidean_Distance_Similarities","Absolute_Value_Distance_Similarities",
                                "Kullback_Leibler_Distance_Similarities","Log_Normalized_Periodogram_Distance_Similarities",
                                "Log_Normalized_Periodogram_Absolute_Value_Distance_Similarities"),c(1:NSIM),clust.method)
}

for(i in 1:NSIM)
{
  Y <- get.sample(G,N)
  if(do_kmeans)
  {
    kmeans_SIMILARITIES <- kmeans_spectral.sims(Y,p,nspans,cluster.sizes)
    Kmeans.Similarities[,,i]<-kmeans_SIMILARITIES[[1]]
  }
  if(nclust.methods>0)
  {
    for(k in 1:nclust.methods)
    {
      SIMILARITIES <- spectral.similarities(Y,p,nspans,cluster.sizes,clust.method[k])
      Similarities.array[,,1,i,k] <- SIMILARITIES[[1]]
      Similarities.array[,,2,i,k] <- SIMILARITIES[[2]]
      Similarities.array[,,3,i,k] <- SIMILARITIES[[3]]
      Similarities.array[,,4,i,k] <- SIMILARITIES[[4]]
      Similarities.array[,,5,i,k] <- SIMILARITIES[[5]]
    }
  }
}

# Use the following to get the average of the simulations
#apply(Similarities.array,c(1:3,5),mean)
#apply(Kmeans.Similarities,1:2,mean)