update.diagMatrix = function(X, w, gate = 0){

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

  # X   is a n x (p+1) matrix

  # w   is a (p+1) x 1 vector

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

  n = nrow(X)

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

  if(gate == 1){

    D = diag(n)*0.25

    return(D)

  }

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

  pp = 1/(1 + exp(-X%*%w))

  D = diag(n)

  for(i in 1:n){

    D[i,i] = max(0.00001,pp[i]*(1-pp[i]))

    #D[i,i] = pp[i]*(1-pp[i])

  }

  return(D)

}

update.zVector= function(X, y, w, D){

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

  # X   is a n x (p+1) matrix

  # y   is a n x 1 vector

  # w   is a (p+1) x 1 vector

  # D   is a n x n matrix

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

  n = nrow(X)

  z = matrix(0,nrow  = n)

  Xw = X%*%w

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

  for(i in 1:n){

    z[i] = Xw[i] + (y[i]-pp[i])/D[i,i]

  }

  return(z)

}

# Soft-thresholding function

Soft.thresholding = function(x,lamb){

  if (lamb==0) {return(x)}

  value1=abs(x)-lamb; value1[round(value1,7)<=0]=0

  value2=sign(x)*value1

  return(value2)

}

# Objective function

objective = function(X,y,w){

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

  # X   is a n x (p+1) matrix

  # y   is a n x 1 vector

  # w   is a (p+1) x 1 vector

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

  value = t(y)%*%X%*%w - sum(log(1+exp(X%*%w)))

  value = value/length(y)

}

logreg.obj = function(X,y,w,M,lambda0, alpha){

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

  # X   is a n x (p+1) matrix

  # y   is a n x 1 vector

  # w   is a (p+1) x 1 vector

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

  s = X%*%w

  v = y * s - log(1+exp(s))

  L = mean(v)

  w1 = w[1:(length(w)-1)]

  R1 = lambda0 *alpha* sum(abs(w1))

  R2 = lambda0 *(1-alpha)*t(w)%*%M%*%w/2 

  R1 + R2 - L

}

abs.logreg.obj = function(X,y,w,M,lambda0, alpha){

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

  # X   is a n x (p+1) matrix

  # y   is a n x 1 vector

  # w   is a (p+1) x 1 vector

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

  s = X%*%w

  v = y * s - log(1+exp(s))

  L = mean(v)

  w1 = w[1:(length(w)-1)]

  R1 = lambda0 *alpha* sum(abs(w1))

  R2 = lambda0 *(1-alpha)*abs(t(w))%*%M%*%abs(w)/2 

  R1 + R2 - L

}

gelnet.logreg.obj2 = function(X,y,w,b, M,l1, l2){

  s <- X %*% w + b

  v <- y * s - log(1+exp(s))

  LL <- mean(v)

  R1 = l1*sum(abs(w))

  R2 = l2*t(w)%*%M%*%w/2

  R1 + R2 - LL

}

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

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

# mian function

# Network-regularized Sparse Logistic Regression (NSLR) 

SGNLR = function(X, y, M, lambda0, alpha, niter=100, gate = 0){

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

  # X   is a n x p matrix

  # y   is a n x 1 vector

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

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

  if(!is.matrix(y)){y = matrix(y,ncol = 1)}

  if(!is.matrix(X)){X = matrix(y,ncol = ncol(X))}

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

  # Adjusting penalty parameters

  lambda = n*alpha*lambda0

  eta = n*(1-alpha)*lambda0

  X = cbind(X,rep(1,n))

  # X   is a n x (p+1) matrix

  if(dim(M)[1]!= p+1){M = cbind(M,rep(0,p)); M = rbind(M,rep(0,p+1))} 

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

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

  # Initialization

  w = matrix(0, nrow = p+1)

  err=0.0001

  obj=log.likelihhod = NULL

  for(i in 1:niter){

    w0 = w

    D = update.diagMatrix(X, w, gate)

    z = update.zVector(X, y, w, D)

    XDX = t(X)%*%D%*%X + eta*M

    t   = t(X)%*%D%*%z

    for(j in 1:p){

      w[j] = 0

      w[j] = (Soft.thresholding(t[j]-t(w)%*%XDX[j,], lambda))/XDX[j,j]

    }

    w[p+1] = 0

    w[p+1] = (t[p+1]-t(w)%*%XDX[p+1,])/XDX[p+1,p+1]

    # Calculate J (for testing convergence)

    log.likelihhod = c(log.likelihhod,objective(X,y,w))

    obj = c(obj,logreg.obj(X,y,w,M,lambda0, alpha))

    if(i>3&&abs(obj[i]-obj[i-1])/abs(obj[i-1])<0.00001){break}

  }

  names(w) = c(paste("w", 1:(length(w)-1), sep = ""), "Intercept")

  return(list(w=w, log.likelihhod = log.likelihhod,obj=obj))

}

# Absolute Network-regularized Logistic Regression (AbsNet.LR)

abs.SNLR = function(X, y, M, lambda0, alpha, niter=100, gate = 0){

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

  # X   is a n x p matrix

  # y   is a n x 1 vector

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

  # example:

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

  if(!is.matrix(y)){y = matrix(y,ncol = 1)}

  if(!is.matrix(X)){X = matrix(y,ncol = ncol(X))}

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

  # Adjusting penalty parameters

  lambda = n*alpha*lambda0

  eta = n*(1-alpha)*lambda0

  X = cbind(X,rep(1,n))

  # X   is a n x (p+1) matrix

  if(dim(M)[1]!= p+1){M = cbind(M,rep(0,p)); M = rbind(M,rep(0,p+1))}

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

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

  # Initialization

  w = matrix(0, nrow = p+1)

  err=0.0001

  obj = NULL

  for(i in 1:niter){

    w0 = w

    D = update.diagMatrix(X, w, gate)

    z = update.zVector(X, y, w, D)

    B = t(X)%*%D%*%X

    t   = t(X)%*%D%*%z

    for(j in 1:p){

      w[j] = 0

      soft.value = lambda + 2*eta*abs(t(w))%*%M[j,]

      w[j] = (Soft.thresholding(t[j]-t(w)%*%B[j,], soft.value))/(B[j,j]+2*eta*M[j,j])

    }

    XDX = B + eta*M

    w[p+1] = 0

    w[p+1] = (t[p+1]-t(w)%*%XDX[p+1,])/XDX[p+1,p+1]

    # Calculate J (for testing convergence)

    obj = c(obj,abs.logreg.obj(X,y,w,M,lambda0, alpha))

    if(i>3&&abs(obj[i]-obj[i-1])/abs(obj[i-1])<0.00001){break}

  }

  names(w) = c(paste("w", 1:(length(w)-1), sep = ""), "Intercept")

  return(list(w=w, obj = obj))

}

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

predict.SGNLR = function(Train.dat,True.label,Est.w){

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

  # X   is a n x p matrix

  # w   is a (p+1) x 1 vector

  # w = (w1,w2,...,b)

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

  library(pROC)

  X = Train.dat; y0 = True.label; w = Est.w

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

  if(!is.matrix(w)){w = matrix(c(w),ncol = 1)}

  X = cbind(X,rep(1,n)) # X is a n x (p+1) matrix

  pp = 1/(1+exp(-X%*%w))

  y = rep(1,n)

  y[pp<0.5] = 0

  accuracy.rate = length(which(y==y0))/n 

  # AUC = accuracy.rate

  AUC = auc(y0, c(pp))[1]

  # AUC = roc(as.factor(y0), c(pp))$auc[1]

  return(list(pp=pp, y=y, acc = accuracy.rate, AUC = AUC))

}