# -------------------------------------------------------------------------

get.penalityMatrix=function(adj,X1, y1){

  feature.num = dim(adj)[2] 

  M.c = diag(0,feature.num)

  M.lasso = diag(0,feature.num)

  M.elastic = diag(1,feature.num)  

  M.network = Non.NormalizedLaplacianMatrix(adj)

  M.AdaptNet = AdaptNet.Non.NormalizedLap(adj,X1, y1)

  return(list(M.c=M.c,M.lasso=M.lasso,M.elastic=M.elastic,M.network=M.network,M.AdaptNet=M.AdaptNet))

}

# -------------------------------------------------------------------------

# Normalized Laplacian Matrix from adjacency matrix

laplacianMatrix = function(adj){

  diag(adj) <- 0

  deg <- apply(abs(adj),1,sum)

  p <- ncol(adj)

  L <- matrix(0,p,p)

  nonzero <- which(deg!=0)

  for (i in nonzero){

    for (j in nonzero){

      L[i,j] <- -adj[i,j]/sqrt(deg[i]*deg[j])

    }

  }

  diag(L) <- 1

  return(L)

}

# Non-Normalized Laplacian Matrix from adjacency matrix

Non.NormalizedLaplacianMatrix = function(adj){

  diag(adj) <- 0

  deg <- apply(adj,1,sum)

  D = diag(deg)

  L = D - adj

  return(L)

}

# ------------------------------------------------------

AdaptNet.penality.matrix = function(adj, X, y){

  library(glmnet)

  glmnet.fit = glmnet(X, y, lambda=0, family='binomial')

  Beta = coef(glmnet.fit)

  p <- ncol(adj)

  coeff.sign = sign(Beta)[2:(p+1)]

  diag(adj) <- 0

  deg <- apply(abs(adj),1,sum)

  L <- matrix(0,p,p)

  nonzero <- which(deg!=0)

  for (i in nonzero){

    for (j in nonzero){

      temp.sign = coeff.sign[i]*coeff.sign[j]

      L[i,j] <- -temp.sign*adj[i,j]/sqrt(deg[i]*deg[j])

    }

  }

  diag(L) <- 1

  return(L_star=L)

}

AdaptNet.Non.NormalizedLap = function(adj,X, y){

  library(glmnet)

  glmnet.fit = glmnet(X, y, lambda=0, family='binomial')

  Beta = coef(glmnet.fit)

  p <- ncol(adj)

  coeff.sign = sign(Beta)[2:(p+1)]

  diag(adj) <- 0

  deg <- apply(abs(adj),1,sum)

  L <- matrix(0,p,p)

  nonzero <- which(deg!=0)

  for (i in nonzero){

    for (j in nonzero){

      temp.sign = coeff.sign[i]*coeff.sign[j]

      L[i,j] <- -temp.sign*adj[i,j]

    }

  }

  diag(L) <- deg

  return(L_star=L)

}

# -------------------------------------------------------------------------

sen.spe = function(pred, truth){

  #----------------------------------------

  # 1: True sample

  # 0: False sample

  # truth = abs(sign(sim.data$w))

  # pred = abs(sign(out1$w[-length(tmp)]))

  #----------------------------------------

  sen = length(which(pred[which(truth==1)]==1))/length(which(truth==1))

  spe = length(which(pred[which(truth==0)]==0))/length(which(truth==0))

  result = c(sen,spe);names(result) = c("sensitivity","specificity")

  return(result)

}

# Prior network regularization

Prior.network = function(adj){

  # adj is a p x p matrix

  # L   is a (p+1) x (p+1) matrix

  L = laplacianMatrix(adj)

  L = rbind(L,rep(0,p))

  L = cbind(L,rep(0,p+1))

  return(L)

}

# -------------------------------------------------------------------------

get_sim_prior_Net = function(n,t,p11,p12){

  A = matrix(0,nrow=n,ncol=n) 

  for(i in 1:n){

    for(j in 1:n){

      if(i>j){

        set.seed(10*i+8*j)

        if(i<t&j<t){

          if(runif(1)<p11) A[i,j] = 1}

        else{

          if(runif(1)<p12) A[i,j] = 1}  

      }

    }

  }

  A = A +t(A)

  diag(A)=0

  return(A)

}

# -------------------------------------------------------------------------

GET.SIM.DATA = function(smaple.num, feature.num, random.seed, snrlam=0.05){

  # smaple.num = 700; feature.num = 100; random.seed = 10; snrlam=0.05

  ii = random.seed

  set.seed(30)

  w  <-c(rnorm(40),rep(0,(feature.num-40)));b = 0

  mu <- rep(0,40)

  Sigma <- matrix(.6, nrow=40, ncol=40) + diag(40)*.4

  set.seed(ii*2)

  X1 <- mvrnorm(n=smaple.num, mu=mu, Sigma=Sigma)

  set.seed(ii*3)

  X2 <- matrix(rnorm(smaple.num*(feature.num-40), mean = 0, sd = 1), nrow = smaple.num, ncol = feature.num-40)

  X = cbind(X1,X2)

  Xw <- -X%*%w 

  pp <- 1/(1+exp(Xw))

  y  <- rep(1,smaple.num)

  y[pp<0.5] <- 0

  p = dim(X)[1];q = dim(X)[2]

  set.seed(ii*3)

  XX = X + snrlam*matrix(rnorm(p*q),ncol=q)

  return(list(X=XX,y=y,w=w))

}

# -------------------------------------------------------------------------

GET.SIM.DATA2 = function(smaple.num, feature.num, random.seed, snrlam=0.05){

  # smaple.num = 700; feature.num = 100; random.seed = 10; snrlam=0.05

  ii = 10

  set.seed(30)

  w  <-c(rnorm(40),rep(0,(feature.num-40)));b = 0

  mu <- rep(0,40)

  Sigma <- matrix(.6, nrow=40, ncol=40) + diag(40)*.4

  set.seed(ii*2)

  X1 <- mvrnorm(n=smaple.num, mu=mu, Sigma=Sigma)

  set.seed(ii*3)

  X2 <- matrix(rnorm(smaple.num*(feature.num-40), mean = 0, sd = 1), nrow = smaple.num, ncol = feature.num-40)

  X = cbind(X1,X2)

  Xw <- -X%*%w 

  pp <- 1/(1+exp(Xw))

  y  <- rep(1,smaple.num)

  y[pp<0.5] <- 0

  p = dim(X)[1];q = dim(X)[2]

  set.seed(random.seed)

  XX = X + snrlam*matrix(rnorm(p*q),ncol=q)

  return(list(X=XX,y=y,w=w))

}