# ========================================
# Breast Cancer Data Measure Script - 
# Stepwise Gaussian Graphical Model
# 
# Works in conjunction with:
# StepGraph R Package
# 
# ========================================

# Clear memory
rm(list=ls())

# Required Libraries
library(Rlab)
library(MASS)

# ----------------------
# LDA Score Function
# ----------------------

LDA_score = function(data_set,omega){
  
  group_true = data_set[[5]]
  
  x = data_set[[2]]
  
  X = data_set[[1]]
  
  p = dim(x)[2]
  
  n = dim(x)[1]
  
  n_RD = data_set[[6]]
  
  n_pCR = data_set[[7]]
  
  S = apply(X, 2, sd)
  
  mean_RD = apply(X[1:83,], 2, mean)
  
  mean_pCR = apply(X[84:112,], 2, mean)
  
  x_RD = scale(x,center=mean_RD,S)
  x_pCR = scale(x,center=mean_pCR,S)
  
  LDA_RD = -0.5*mahalanobis(x_RD, center=F, cov = omega,inverted = T)^2+log(n_RD/n)
  LDA_pCR = -0.5*mahalanobis(x_pCR, center=F, cov = omega,inverted = T)^2+log(n_pCR/n)
  
  group = ifelse(LDA_pCR>=LDA_RD,1,0)
  
  TP = sum(group[group_true==1]==1)
  
  FN = sum(group[group_true==1]==0)
  
  TN = sum(group[group_true==0]==0)
  
  FP = sum(group[group_true==0]==1)
  
  sensitivity = TP/sum(group_true==1)
  
  specificity = TN/sum(group_true==0)
  
  MCC = ((TP*TN)-(FP*FN))/((TP+FP)*(TP+FN)*(TN+FP)*(TN+FN))^0.5
  
  edges = (sum(abs(omega - diag(diag(omega)))>1e-5))/2
  
  return(list(sensitivity,specificity,MCC,edges))}

# ----------------------
# StepGraph Measures
# ----------------------

R = 100

LDA_list = list()
for (r in 1:R){
  load(file = paste("data_set_", r,".RData",sep = ""))
  load(file = paste("list_opt_", r,".RData",sep = ""))
  LDA_list[[r]] = LDA_score(data_set,list)
}

Sensitivity_mean_list = mean(unlist(lapply(1:R, function(r) LDA_list[[r]][[1]])))
Sensitivity_sd_list = sd(unlist(lapply(1:R, function(r) LDA_list[[r]][[1]])))/sqrt(R)

Specificity_mean_list =  mean(unlist(lapply(1:R, function(r) LDA_list[[r]][[2]])))
Specificity_sd_list =  sd(unlist(lapply(1:R, function(r) LDA_list[[r]][[2]])))/sqrt(R)

MCC_mean_list =   mean(unlist(lapply(1:R, function(r) LDA_list[[r]][[3]])),na.rm = TRUE)
MCC_sd_list =   sd(unlist(lapply(1:R, function(r) LDA_list[[r]][[3]])),na.rm = TRUE)/sqrt(R)

Edges_mean_list =   mean(unlist(lapply(1:R, function(r) LDA_list[[r]][[4]])))
Edges_sd_list =   sd(unlist(lapply(1:R, function(r) LDA_list[[r]][[4]])))/sqrt(R)

# ----------------------
# CVglasso Measures
# ----------------------

LDA_list_glasso = list()

for (r in 1:R){
  load(file = paste("data_set_", r,".RData",sep = ""))
  load(file = paste("glasso_opt_", r,".RData",sep = ""))
  LDA_list_glasso[[r]] = LDA_score(data_set,list_glasso)
}

Sensitivity_mean_list_glasso = mean(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[1]])))
Sensitivity_sd_list_glasso = sd(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[1]])))/sqrt(R)

Specificity_mean_list_glasso =  mean(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[2]])))
Specificity_sd_list_glasso =  sd(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[2]])))/sqrt(R)

MCC_mean_list_glasso =   mean(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[3]])),na.rm = TRUE)
MCC_sd_list_glasso =   sd(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[3]])),na.rm = TRUE)/sqrt(R)

Edges_mean_list_glasso =   mean(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[4]])))
Edges_sd_list_glasso =   sd(unlist(lapply(1:R, function(r) LDA_list_glasso[[r]][[4]])))/sqrt(R)


# ----------------------
# CLIME Measures
# ----------------------

LDA_list_clime = list()

for (r in 1:R){
  load(file = paste("data_set_", r,".RData",sep = ""))
  load(file = paste("clime_opt_", r,".RData",sep = ""))
  LDA_list_clime[[r]] = LDA_score(data_set,list_clime)
}

Sensitivity_mean_list_clime = mean(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[1]])))
Sensitivity_sd_list_clime = sd(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[1]])))/sqrt(R)

Specificity_mean_list_clime =  mean(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[2]])))
Specificity_sd_list_clime =  sd(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[2]])))/sqrt(R)

MCC_mean_list_clime =   mean(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[3]])),na.rm = TRUE)
MCC_sd_list_clime =   sd(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[3]])),na.rm = TRUE)/sqrt(R)

Edges_mean_list_clime =   mean(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[4]])))
Edges_sd_list_clime =   sd(unlist(lapply(1:R, function(r) LDA_list_clime[[r]][[4]])))/sqrt(R)

# ----------------------
# Output Printing
# ----------------------

# StepGraph

round(Sensitivity_mean_list, digits = 3)
round(Sensitivity_sd_list, digits = 3)

round(Specificity_mean_list, digits = 3)
round(Specificity_sd_list, digits = 3)

round(MCC_mean_list, digits = 3)
round(MCC_sd_list, digits = 3)

round(Edges_mean_list, digits = 3)
round(Edges_sd_list, digits = 3)

# CLIME

round(Sensitivity_mean_list_clime, digits = 3)
round(Sensitivity_sd_list_clime, digits = 3)

round(Specificity_mean_list_clime, digits = 3)
round(Specificity_sd_list_clime, digits = 3)

round(MCC_mean_list_clime, digits = 3)
round(MCC_sd_list_clime, digits = 3)

round(Edges_mean_list_clime, digits = 3)
round(Edges_sd_list_clime, digits = 3)

# CVglasso

round(Sensitivity_mean_list_glasso, digits = 3)
round(Sensitivity_sd_list_glasso, digits = 3)

round(Specificity_mean_list_glasso, digits = 3)
round(Specificity_sd_list_glasso, digits = 3)

round(MCC_mean_list_glasso, digits = 3)
round(MCC_sd_list_glasso, digits = 3)

round(Edges_mean_list_glasso, digits = 3)
round(Edges_sd_list_glasso, digits = 3)



