# ========================================
# Simulation Script - 
# Reproduce the Toy Example of Section 2.5
# Works in conjunction with:
# StepGraph R Package
# 
# ========================================

# Clear memory
rm(list=ls())

# Required Libraries


library(corpcor)
library(Rlab)
library(MASS)
library(scatterplot3d)
library(huge)
library(glmnet)
library(clime)
library(sna)
library(Matrix)
library(Tlasso)
library(CVglasso)
library(igraph)
library(StepGraph)

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


# =============================
# Required Auxiliary functions 
# =============================

# Adjacency matrix


Adj_mat = function(Omega){

p = dim(Omega)[1]

Adj_mat = matrix(0,p,p)

 for (i in 1:p){
    for (j in 1:p){
      if (abs(Omega[i,j])<=1e-5){Adj_mat[i,j]=0}
      else (Adj_mat[i,j]=1)
    } 
 }

return(Adj_mat)}


# Adjacency matrix in each step

Adj_mat_step = function(Omega,step){

p = dim(Omega)[1]

Adj_mat = matrix(0,p,p)

 for (i in 1:p){
    for (j in 1:p){
      if (abs(Omega[i,j])==step){Adj_mat[i,j]=1}
      else (Adj_mat[i,j]=0)
    } 
 }

return(Adj_mat)}







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

# -----------------------------------
# Simulate Gaussian Data NN(2)Model
# -----------------------------------


set.seed(1234567)

# Setting the parameters


p = 20 # Dimension

n = 1000 # Sample Size


Omega<- NeighborOmega(p, sd = 1, knn=2, norm.type = 2)

save(Omega, file = "Omega_NN_p_20.RData") 

Sigma<- solve(Omega)




# 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_NN_p20_n1000_", ".RData",sep = ""))


###########################################################
###########################################################
###########################################################

load(file = paste("sample_NN_p20_n1000_", ".RData",sep = ""))


# -----------------------------------
# Estimation of the set of Edges (graph)
# -----------------------------------


res = cv.StepGraph(X.data[[1]],fold=5,alpha_f_min=0.10,alpha_f_max=0.75,n_alpha=10,nei.max=10)
fit.GStepwise = StepGraph(scale(X.data[[1]]),alpha_f=res[[1]][[1]],alpha_b=res[[2]][[1]],nei.max=10)

#Optimal thresholds
print('Optimal alpha forward'); round( res[[1]][[1]],2)

print('Optimal alpha backward'); round( res[[2]][[1]],2)



# -----------------------------------
# True graph and StepGraph step by step 
# -----------------------------------

load(file='Omega_NN_p_20.RData')

# True graph

g.NN.true = graph.adjacency(Adj_mat(Omega), mode = "undirected", weighted = TRUE , diag=F)

pdf("Model_NN_p20_n1000_Graph.pdf",width=8,height=6)

plot(g.NN.true,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

######################################################################
######################################################################

load(file = paste("GStepwise_NN_p20_n1000_", 1,".RData",sep = ""))

g.NN.Gstepwise = graph.adjacency(Adj_mat(fit.GStepwise[[6]]), mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p30_n1000_Graph.pdf",width=8,height=6)

plot(g.NN.Gstepwise,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()


# StepGraph  step by step 

Omega = fit.GStepwise[[4]]

E = Omega[upper.tri(Omega, diag = FALSE)]

Steps = sort(E[E>0])

load(file = paste("GStepwise_NN_p20_n1000_", 1,".RData",sep = ""))

# Step 1

Mat.1 = Adj_mat_step(fit.GStepwise[[4]],step=1)

g.NN.Gstepwise.1 = graph.adjacency(Mat.1, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_1.pdf",width=8,height=6)

plot(g.NN.Gstepwise.1,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 2

Mat.2 = Adj_mat_step(fit.GStepwise[[4]],step=2) + Mat.1

g.NN.Gstepwise.2 = graph.adjacency(Mat.2, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_2.pdf",width=8,height=6)

plot(g.NN.Gstepwise.2,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 3

Mat.3 = Adj_mat_step(fit.GStepwise[[4]],step=3) + Mat.2

g.NN.Gstepwise.3 = graph.adjacency(Mat.3, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_3.pdf",width=8,height=6)

plot(g.NN.Gstepwise.3,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 4

Mat.4 = Adj_mat_step(fit.GStepwise[[4]],step=4) + Mat.3

g.NN.Gstepwise.4 = graph.adjacency(Mat.4, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_4.pdf",width=8,height=6)

plot(g.NN.Gstepwise.4,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 5

Mat.5 = Adj_mat_step(fit.GStepwise[[4]],step=5) + Mat.4

g.NN.Gstepwise.5 = graph.adjacency(Mat.5, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_5.pdf",width=8,height=6)

plot(g.NN.Gstepwise.5,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 6

Mat.6 = Adj_mat_step(fit.GStepwise[[4]],step=6) + Mat.5

g.NN.Gstepwise.6 = graph.adjacency(Mat.6, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_6.pdf",width=8,height=6)

plot(g.NN.Gstepwise.6,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 7

Mat.7 = Adj_mat_step(fit.GStepwise[[4]],step=7) + Mat.6

g.NN.Gstepwise.7 = graph.adjacency(Mat.7, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_7.pdf",width=8,height=6)

plot(g.NN.Gstepwise.7,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 8

Mat.8 = Adj_mat_step(fit.GStepwise[[4]],step=8) + Mat.7

g.NN.Gstepwise.8 = graph.adjacency(Mat.8, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_8.pdf",width=8,height=6)

plot(g.NN.Gstepwise.8,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 9

Mat.9 = Adj_mat_step(fit.GStepwise[[4]],step=9) + Mat.8

g.NN.Gstepwise.9 = graph.adjacency(Mat.9, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_9.pdf",width=8,height=6)

plot(g.NN.Gstepwise.9,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 10

Mat.10 = Adj_mat_step(fit.GStepwise[[4]],step=10) + Mat.9

g.NN.Gstepwise.10 = graph.adjacency(Mat.10, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_10.pdf",width=8,height=6)

plot(g.NN.Gstepwise.10,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 11

Mat.11 = Adj_mat_step(fit.GStepwise[[4]],step=11) + Mat.10

g.NN.Gstepwise.11 = graph.adjacency(Mat.11, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_11.pdf",width=8,height=6)

plot(g.NN.Gstepwise.11,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 12

Mat.12 = Adj_mat_step(fit.GStepwise[[4]],step=12) + Mat.11

g.NN.Gstepwise.12 = graph.adjacency(Mat.12, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_12.pdf",width=8,height=6)

plot(g.NN.Gstepwise.12,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 13

Mat.13 = Adj_mat_step(fit.GStepwise[[4]],step=13) + Mat.12

g.NN.Gstepwise.13 = graph.adjacency(Mat.13, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_13.pdf",width=8,height=6)

plot(g.NN.Gstepwise.13,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 14

Mat.14 = Adj_mat_step(fit.GStepwise[[4]],step=14) + Mat.13

g.NN.Gstepwise.14 = graph.adjacency(Mat.14, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_14.pdf",width=8,height=6)

plot(g.NN.Gstepwise.14,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 15

Mat.15 = Adj_mat_step(fit.GStepwise[[4]],step=15) + Mat.14

g.NN.Gstepwise.15 = graph.adjacency(Mat.15, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_15.pdf",width=8,height=6)

plot(g.NN.Gstepwise.15,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()

# Step 16

Mat.16 = Adj_mat_step(fit.GStepwise[[4]],step=16) + Mat.15

g.NN.Gstepwise.16 = graph.adjacency(Mat.16, mode = "undirected", weighted = TRUE , diag=F)

pdf("GStepwise_NN_p20_n1000_Graph_Step_16.pdf",width=8,height=6)

plot(g.NN.Gstepwise.16,edge.color='black',vertex.color="white",layout=layout_in_circle)

dev.off()








