##### define a negative composite log likelihood

qf_k <- function(x, K, p,n)
{
  
  lm0 <- 0
  
  lm3 <- matrix(0,p,n)
  
  lm4 <- matrix(0,p,n)
  
  lm5 <- rep(0,p)
  
  for (j in 1:p)
  {
    lm0 = lm0 + log(K[j,j])
    
    for (k in 1:n)
    {
      for (i in 1:p)
      {
        if(i != j)
        {
          lm3[j,k] = lm3[j,k] + K[i,j]*x[k,i]
        }
        
      }
      lm4[j,k] =  K[j,j]*(x[k,j] + 1/K[j,j]*lm3[j,k])^2
      
    }
    lm5[j] = lm5[j] + sum(lm4[j,])
    
  }
  

  f_k <- -1/2*n*lm0 + 1/2*sum(lm5)
  
  return(f_k)
  
}

ccf_k <- qf_k(x, K, p,n) # -l_c(\theta) negative composite log-likelihood
ccf_k

## grisearch function is used to do grid search for different turning parameters

gridsearch <- function(n,gridpoints,p,x)
  
{
  le <- length(gridpoints)
  
  compbic_p <- Inf
  compbic <- Inf
  turning_b <- rep(0,4)
  tau_t = lambda_1 = lambda_2 = lambda_3 = 0
  
  
  for (i in 1:le)
  {
    tau_t = gridpoints[i]
    
    for (j in 1:le)
    {
      lambda_1 = lambda_2 = lambda_3 = gridpoints[j]
      K <- ccselection(x, p, tau_t, lambda_1, lambda_2, lambda_3)[[1]]
      SK <- MBVB(K,tau_t)
      #print(K)  
      compbic <- 2*qf_k(x,  SK[[2]], p, n) + SK[[3]]*log(n)
      
      if(compbic < compbic_p)
      {
        compbic_p = compbic
        turning_b = c(tau_t,lambda_1,lambda_2,lambda_3)
        color <- SK
        #print(compbic_p)
        #print(turning_b)
      }
    }
    
  }
  
  #compbic_p
  #turning_b

  return(list(color[[1]],color[[2]],turning_b))
}


timeg <- proc.time()
rgridsearch <- gridsearch(n,gridpoints,p,x) 
rgridsearch
timege <- proc.time()-timeg
timege

# star p = 10,20,30, tau_t=0.1, lambda_1 = lambda_2 = lambda_3 = 0.01 which we have good estimate and low computational cost

repeatdata <- 1
timestart <- proc.time()
for (q in 1:repeatdata)
{
  mu <- rep(0, p)
  x <- mvrnorm(n,mu,sigma) # observations
  
  print(q)
  
  gridpoints <- c(0.01, 0.05, 0.1, 0.5, 1)
  #gridpoints <- c(0.01, 0.1)
  
  rgridsearch <- gridsearch(n,gridpoints,p,x) 
  
  #getwd()
  #setwd("C:/Users/Administrator/Desktop/results(Xiaoying)")
  #write.table(rgridsearch[[1]], file="bestmodel_star_n1000p30.txt", append=T, row.names=FALSE, col.names=FALSE, sep=" ",eol = "\n")
  #write.table(rgridsearch[[2]],file="bestestimate_star_n1000p30.txt", append=T, row.names=FALSE, col.names=FALSE, sep=" ",eol = "\n")
  #write.table(color[[1]], file="E:/backup (2017.7.10)/from desktop/xiaoying/simulations/coloredmodel.txt", append=T, row.names=FALSE, col.names=FALSE, sep=" ",eol = "\n")
  #write.table(color[[2]], file="E:/backup (2017.7.10)/from desktop/xiaoying/simulations/coloredestimates.txt", append=T, row.names=FALSE, col.names=FALSE, sep=" ",eol = "\n")
  #write.table(mse, file="E:/backup (2017.7.10)/from desktop/xiaoying/simulations/mse.txt", append=T, row.names=FALSE, col.names=FALSE, sep=" ",eol = "\n")
}

timeend <- proc.time()-timestart
timeend

#rgridsearch

mses <- read.table("bestmodel_star_n1000p4.txt", header=F)
mses

coloredmodels <- read.table("E:/backup (2017.7.10)/from desktop/xiaoying/simulations/coloredmodel.txt", header=F)
coloredmodels
coloredmodels[1476:1500,]

setwd("C:/Users/Administrator/Desktop/results(Xiaoying)/new")
#coloredestimates <- read.table("E:/backup (2017.7.10)/from desktop/xiaoying/simulations/star p30 n500/coloredestimates.txt", header=F)
coloredestimates <- read.table("bestestimate_star_n1000p30.txt", header=F)
dim(coloredestimates)
coloredestimates[1:p,]



# computing normalized mean squared errors
nmse <- rep(0,repeatdata)
for (i in 1:repeatdata)
{
  j = 1+p*(i-1)
  k = p+p*(i-1)
  A <- K_0 - as.matrix(coloredestimates[j:k,])
  nmse[i] = sum(A^2)/sum((K_0)^2) #normalized mean squared errors
}
mean(nmse)
sd(nmse)

# computing Specificity, Sensitivity and MCC
TP = TN = FP = FN <- rep(0,repeatdata)
Spe = Sen <- rep(0,repeatdata)
MCC = Fscore = d0 <- rep(0,repeatdata)
for (i in 1:repeatdata)
{
  j = 1+p*(i-1)
  k = p+p*(i-1)
  B <- as.matrix(coloredestimates[j:k,])
  for (i1 in 1:(p-1))
  {
    for (i2 in (i1+1):p)
    {
      if(K_0[i1,i2]!=0 & B[i1,i2]!=0)
      {
        TP[i] = TP[i] + 1 
      }
      
      if(K_0[i1,i2]==0 & B[i1,i2]==0)
      {
        TN[i] = TN[i] + 1
      }
      
      if(K_0[i1,i2]!=0 & B[i1,i2]==0)
      {
        FP[i] = FP[i] + 1
      }
      
      if(K_0[i1,i2]==0 & B[i1,i2]!=0)
      {
        FN[i] = FN[i] + 1
      }
      
    }
  }
  Spe[i] = TN[i]/(TN[i]+FP[i]) #specificity
  Sen[i] = TP[i]/(TP[i]+FN[i]) #sensitivity
  MCC[i] = (TP[i]*TN[i]-FP[i]*FN[i])/sqrt((TP[i]+FP[i])*(TP[i]+FN[i])*(TN[i]+FP[i])*(TN[i]+FN[i]))
  Fscore[i] = 2*TP[i]/(2*TP[i]+FP[i]+FN[i])
  d0[i] = (TN[i]+TP[i])/(p*(p-1)/2)
}

mean(Spe)
sd(Spe)
mean(Sen)
sd(Sen)
mean(MCC)
sd(MCC)
mean(d0)
sd(d0)
mean(Fscore)
sd(Fscore)

vma = max(diag(model))
vmi = min(diag(model)[diag(model)!=1])
ema = max(model[row(model)!=col(model)])
emi = min(model[row(model)!=col(model)& model!=0])

ACC <- rep(0,repeatdata)
vpro <- rep(0,repeatdata)
epro <- rep(0,repeatdata)

for (q in 1:repeatdata)
{
  j = 1+p*(i-1)
  k = p+p*(i-1)
  B <- as.matrix(coloredestimates[j:k,])
  
w = rep(0,p) # the number of correct (positive estimation) for free vertex

for (j in 1:p)
{
  if(model[j,j] == 1) 
  {
    for(k in 1:p)
    {
      if(j!=k & B[j,j] != B[k,k]) 
      {
        w[j] = w[j] + 1
      }
        
    }
  }
    
}

w

ww <- rep(0,max(model))
wu <- rep(0,max(model))

for(i in vmi:vma)
{
 for (j in 1:p)
 {
  if(model[j,j] == i) 
  {
    wu[i] = wu[i] + 1
    for(k in 1:p)
    {
      if(j!=k & model[k,k] == i & B[k,k] == B[j,j])
      {
        ww[i] = ww[i] + 1
      }
      if(j!=k & model[k,k] != i & B[k,k] != B[j,j])
      {
        ww[i] = ww[i] + 1
      } 
    } 
  }
    
 }
}

ww
wu

ww1 <- rep(0,max(model))
wu1 <- rep(0,max(model))

for (i in emi:ema)
{
  for (j in 1:(p-1))
  {
    for (j1 in (j+1):p)
    {
      if(model[j,j1] == i)
      {
        wu1[i] = wu1[i] + 1
        for(k in 1:(p-1))
        {
          for(k1 in (k+1):p)
          {
            if(!all(c(j,j1)==c(k,k1))& model[k,k1] == i & B[k,k1] == B[j,j1])
            {
              ww1[i] = ww1[i] + 1
            }
            if(!all(c(j,j1)==c(k,k1)) & model[k,k1] != i & B[k,k1] != B[j,j1])
            {
              ww1[i] = ww1[i] + 1
            } 
          }
        }
      }
    }
  }
}

ww1
wu1
vpro[q] = sum(ww/wu/(p-1),na.rm=TRUE)+sum(w/(p-1))
epro[q] = sum(ww1/wu1/(p*(p-1)/2-1),na.rm=TRUE)

ACC[q] = (d0[q]+vpro[q]+epro[q])/(1+length(wu[wu!=0])+length(wu1[wu1!=0])+length(w[w!=0]))

}



mean(nmse)#normalized mean squared error
sd(nmse)

mean(Fscore) # F1 score
sd(Fscore)

mean(d0) #measure zero elements
sd(d0)

mean(ACC) # measure all
sd(ACC)

#read.table("E:/backup (2017.7.10)/from desktop/xiaoying/simulations/turnparameter.txt", header=F)