# ========================================
# Simulation Script - 
# Stepwise Gaussian Graphical Model
# 
# Works in conjunction with:
# StepGraph R Package
# 
# ========================================

# Clear memory
rm(list=ls())

# Required Libraries
library(Rlab)
library(MASS)
library(huge)
library(glmnet)
library(clime)
library(CVglasso)
library(psych)
library(StepGraph)

# Required Source Files
source("Classification_Performance.R")
source("Numerical_Performance.R")

# =======================
# Simulation Experiment 
# =======================

# -----------------------------------
# Simulate Gaussian Data AR(1) Model
# -----------------------------------

# Covariance matrix of the AR(1) model

# SigmaAR1 is a function that compute the covariance matrix of AR(1) model

SigmaAR1 = function(p,phi){
  Sigma = diag(p)
  
  for (i in 1:p) 
  {
    for (j in 1:p) 
    {
      Sigma[i,j] = phi^(abs(i-j))
      
    }
  }
  
  return(Sigma)
  
}

# (!!!) INITIAL TIME
t0 <- Sys.time()

# Setting the parameters
set.seed(1234567)
phi=0.4  # coef of AR(1)
p = 50 # Dimension or number of variables
n = 100 # Sample Size
R = 50 # Number of replications

# Generate the covariance matrix
Sigma <- SigmaAR1(p,phi)

# Compute the precision matrix by inverting the covariance matrix
Omega = solve(Sigma)  


# Generate Data from a Gaussian distribution 
X.data = list()
X.data = lapply(1:100, function(r) mvrnorm(n, mu=rep(0,p), Sigma))
save(X.data, file = paste("sample_AR1_p50_n100_", ".RData",sep = ""))

# -----------------------------------
# Application of Algorithms to Data
# -----------------------------------

# Graphical Stepwise
for (r in 1:R){
  cat(r)
  res = cv.StepGraph(x=X.data[[r]],fold=5,alpha_f_min=0.10,alpha_f_max=0.75,n_alpha=10,nei.max=10)
  fit.GStepwise = StepGraph(scale(X.data[[r]]),alpha_f=res[[1]][[1]],alpha_b=res[[2]][[1]],nei.max=10)
  save(res, file = paste("alpha_AR1_p50_n100_", r,".RData",sep = ""))
  save(fit.GStepwise, file = paste("GStepwise_AR1_p50_n100_", r,".RData",sep = ""))
}

# Graphical lasso with CVglasso
for (r in 1:R){
  cat(r)
  fit.Glasso = CVglasso(scale(X.data[[r]]), K=5)$Omega    
  save(fit.Glasso, file = paste("Glasso_AR1_p50_n100_", r,".RData",sep = ""))
}

# CLIME
for (r in 1:R){
  lambda = cv.clime(clime(X.data[[r]]),loss="likelihood",fold=5)$lambdaopt
  fit.CLIME = clime(X.data[[r]],lambda)$Omega[[1]]
  save(lambda, file = paste("lambda_p50_n100_", r,".RData",sep = ""))
  save(fit.CLIME, file = paste("CLIME_AR1_p50_n100_", r,".RData",sep = ""))
}

# -----------------------------------------
# Evaluation of Classification Performance
# -----------------------------------------

# Classification Performance Measures for ***StepGraph***
performance_list = list()
for (r in 1:R){
  load(file = paste("GStepwise_AR1_p50_n100_", r,".RData",sep = ""))
  performance_list[[r]] = performance(fit.GStepwise[[6]],Omega)
}
# Sensitivity
Sensitivity_mean = mean(unlist(lapply(1:R, function(r) performance_list[[r]][[2]])),NaN.rm=TRUE)
Sensitivity_sd = sd(unlist(lapply(1:R, function(r) performance_list[[r]][[2]])))/sqrt(R)
# Specificity
Specificity_mean = mean(unlist(lapply(1:R, function(r) performance_list[[r]][[1]])),NaN.rm=TRUE)
Specificity_sd = sd(unlist(lapply(1:R, function(r) performance_list[[r]][[1]])))/sqrt(R)
# MCC
MCC_mean = mean(unlist(lapply(1:R, function(r) performance_list[[r]][[3]])),NaN.rm=TRUE)
MCC_sd = sd(unlist(lapply(1:R, function(r) performance_list[[r]][[3]])))/sqrt(R)

# Classification Performance Measures for ***CVglasso***
performance_glasso = list()
for (r in 1:R){
  load(file = paste("Glasso_AR1_p50_n100_", r,".RData",sep = ""))
  performance_glasso[[r]] = performance(fit.Glasso,Omega)
}
# Sensitivity
Sensitivity_glasso_mean = mean(unlist(lapply(1:R, function(r) performance_glasso[[r]][[2]])), NaN.rm=TRUE)
Sensitivity_glasso_sd = sd(unlist(lapply(1:R, function(r) performance_glasso[[r]][[2]])))/sqrt(R)
# Specificity
Specificity_glasso_mean = mean(unlist(lapply(1:R, function(r) performance_glasso[[r]][[1]])),NaN.rm=TRUE)
Specificity_glasso_sd = sd(unlist(lapply(1:R, function(r) performance_glasso[[r]][[1]])))/sqrt(R)
# MCC
MCC_glasso_mean = mean(unlist(lapply(1:R, function(r) performance_glasso[[r]][[3]])),na.rm=TRUE)
MCC_glasso_sd = sd(unlist(lapply(1:R, function(r) performance_glasso[[r]][[3]])), na.rm=TRUE)/sqrt(R)

# Classification Performance Measures for ***CLIME***
performance_CLIME = list()
for (r in 1:R){
  load(file = paste("CLIME_AR1_p50_n100_", r,".RData",sep = ""))
  performance_CLIME[[r]] = performance(fit.CLIME,Omega)
}
# Sensitivity
Sensitivity_CLIME_mean = mean(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[2]])),NaN.rm=TRUE)
Sensitivity_CLIME_sd = sd(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[2]])))/sqrt(R)
# Specificity
Specificity_CLIME_mean = mean(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[1]])),NaN.rm=TRUE)
Specificity_CLIME_sd = sd(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[1]])))/sqrt(R)
# MCC
MCC_CLIME_mean = mean(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[3]])),NaN.rm=TRUE)
MCC_CLIME_sd = sd(unlist(lapply(1:R, function(r) performance_CLIME[[r]][[3]])))/sqrt(R)


# Summary of Classification PerformanceOuput
P = matrix(0,6,3)
P[1,] = c(Sensitivity_mean,Specificity_mean,MCC_mean)
P[2,] = c(Sensitivity_sd,Specificity_sd,MCC_sd)
P[3,] = c(Sensitivity_glasso_mean,Specificity_glasso_mean,MCC_glasso_mean)
P[4,] = c(Sensitivity_glasso_sd,Specificity_glasso_sd,MCC_glasso_sd)
P[5,] = c(Sensitivity_CLIME_mean,Specificity_CLIME_mean,MCC_CLIME_mean)
P[6,] = c(Sensitivity_CLIME_sd,Specificity_CLIME_sd,MCC_CLIME_sd)
colnames(P) = c("Sensitivity","Specificity","MCC")
P <- round(P,digits=3)
P

# -------------------------------------
# Evaluation of Numerical Performance
# -------------------------------------

# Numerical Performance for ***StepGraph***
num_performance_list = list()
for (r in 1:R){
  load(file = paste("GStepwise_AR1_p50_n100_", r,".RData",sep = ""))
  num_performance_list[[r]] = Numerical.Performance(fit.GStepwise[[6]],Omega)
}
# LRT
LRT_mean = mean(unlist(lapply(1:R, function(r) num_performance_list[[r]][[1]])),NaN.rm=TRUE)
LRT_sd = sd(unlist(lapply(1:R, function(r) num_performance_list[[r]][[1]])))/sqrt(R)
# MSE
MSE_mean = mean(unlist(lapply(1:R, function(r) num_performance_list[[r]][[2]])),NaN.rm=TRUE)
MSE_sd = sd(unlist(lapply(1:R, function(r) num_performance_list[[r]][[2]])))/sqrt(R)
# NLRT
NLRT_mean = mean(unlist(lapply(1:R, function(r) num_performance_list[[r]][[3]])),NaN.rm=TRUE)
NLRT_sd = sd(unlist(lapply(1:R, function(r) num_performance_list[[r]][[3]])))/sqrt(R)

# Numerical performance for ***CVglasso***
num_performance_glasso = list()
for (r in 1:R){
  load(file = paste("Glasso_AR1_p50_n100_", r,".RData",sep = ""))
  num_performance_glasso[[r]] = Numerical.Performance(fit.Glasso,Omega)
}
# LRT
LRT_glasso_mean = mean(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[1]])),NaN.rm=TRUE)
LRT_glasso_sd = sd(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[1]])))/sqrt(R)
# MSE
MSE_glasso_mean = mean(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[2]])),NaN.rm=TRUE)
MSE_glasso_sd = sd(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[2]])))/sqrt(R)
# NLRT
NLRT_glasso_mean = mean(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[3]])),NaN.rm=TRUE)
NLRT_glasso_sd = sd(unlist(lapply(1:R, function(r) num_performance_glasso[[r]][[3]])))/sqrt(R)

# Numerical performance for ***CLIME***
num_performance_CLIME = list()
for (r in 1:R){
  load(file = paste("CLIME_AR1_p50_n100_", r,".RData",sep = ""))
  num_performance_CLIME[[r]] = Numerical.Performance(fit.CLIME,Omega)
}
# LRT
LRT_CLIME_mean = mean(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[1]])),NaN.rm=TRUE)
LRT_CLIME_sd = sd(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[1]])))/sqrt(R)
# MSE
MSE_CLIME_mean = mean(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[2]])),NaN.rm=TRUE)
MSE_CLIME_sd = sd(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[2]])))/sqrt(R)
# NLRT
NLRT_CLIME_mean = mean(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[3]])),NaN.rm=TRUE)
NLRT_CLIME_sd = sd(unlist(lapply(1:R, function(r) num_performance_CLIME[[r]][[3]])))/sqrt(R)

# Summary of Classification PerformanceOuput
H = matrix(0,6,3)
H[1,] = c(LRT_mean,NLRT_mean, MSE_mean)
H[2,] = c(LRT_sd,NLRT_sd, MSE_sd)
H[3,] = c(LRT_glasso_mean,NLRT_glasso_mean, MSE_glasso_mean)
H[4,] = c(LRT_glasso_sd, NLRT_glasso_sd, MSE_glasso_sd)
H[5,] = c(LRT_CLIME_mean,NLRT_CLIME_mean, MSE_CLIME_mean)
H[6,] = c(LRT_CLIME_sd,NLRT_CLIME_sd,MSE_CLIME_sd)
colnames(H) = c("LRT", "NormLRT","Frobenious")
H <- round (H,3)
H

# (!!!) FINAL TIME
t1 <- Sys.time()  # final time
tt <- t1-t0 # total time



