#hazard and cumulative hazard functions for proportional hazards (gamma = 0)/odds models (gamma = 1)
Lambda = function(u,beta,gamma,z){
  if(gamma == 0) return(exp(u+z%*%beta))
  else return(log(1+gamma*exp(u+z%*%beta))/gamma)
}

lambda = function(u,beta,gamma,z){
  if(gamma == 0) return(exp(u+z%*%beta))
  else return(exp(u+z%*%beta)/(1+gamma*exp(u+z%*%beta)))
}

dlambda = function(u,beta,gamma,z){
  if(gamma == 0) return(exp(u+z%*%beta))
  else return(exp(u+z%*%beta)/(1+gamma*exp(u+z%*%beta))^2)
}

ddlambda = function(u,beta,gamma,z){
  if(gamma == 0) return(exp(u+z%*%beta))
  else return(exp(u+z%*%beta)*(1-(gamma^2)*exp(2*(u+z%*%beta)))/(1+gamma*exp(u+z%*%beta))^4)
}

#inverse CDF
invF = function(i,u,gamma){
  if(gamma == 0) return(log(-log(1-u[,i])))
  else return(log(((1/(1-u[,i])^gamma)-1)/gamma))
} 


#delta, z should be ordered according to the ranks of time
lik = function(i,sampH,delta,z,beta,gamma){
  p = length(beta)
  betanull = numeric(p)
  L = ((lambda(sampH[,i],beta,gamma,z)^delta)*exp(-Lambda(sampH[,i],beta,gamma,z)))/
      ((lambda(sampH[,i],betanull,gamma,z)^delta)*exp(-Lambda(sampH[,i],betanull,gamma,z)))
          
  return(prod(L))
} 

liknull = function(i,sampH,delta,z,beta,gamma){
  L = (lambda(sampH[,i],beta,gamma,z)^delta)*exp(-Lambda(sampH[,i],beta,gamma,z))
          
  return(prod(L))
} 

dlik = function(i,sampH,delta,z,beta,gamma){
  a = lik(i,sampH,delta,z,beta,gamma)
  b = as.numeric(delta*dlambda(sampH[,i],beta,gamma,z)/lambda(sampH[,i],beta,gamma,z)-lambda(sampH[,i],beta,gamma,z))*z
  c = apply(b,2,sum)
  return(a*c)
}  

ndlik = function(i,sampH,delta,z,beta,gamma){
  a = lik(i,sampH,delta,z,beta,gamma)
  b = as.numeric(delta*dlambda(sampH[,i],beta,gamma,z)/lambda(sampH[,i],beta,gamma,z)-lambda(sampH[,i],beta,gamma,z))
  return(a*b)
}  

ddlik = function(i,sampH,delta,z,beta,gamma){
  a = lik(i,sampH,delta,z,beta,gamma)
  b = as.numeric(delta*dlambda(sampH[,i],beta,gamma,z)/lambda(sampH[,i],beta,gamma,z)-lambda(sampH[,i],beta,gamma,z))*z
  c = apply(b,2,sum)
  d = c%*%t(c)
  e = as.numeric(delta*(ddlambda(sampH[,i],beta,gamma,z)*lambda(sampH[,i],beta,gamma,z)-dlambda(sampH[,i],beta,gamma,z)*dlambda(sampH[,i],beta,gamma,z))/(lambda(sampH[,i],beta,gamma,z)*lambda(sampH[,i],beta,gamma,z))-dlambda(sampH[,i],beta,gamma,z))*z
  f = t(e)%*%z
  return(a*(d+f))
}  

nddlik = function(i,sampH,delta,z,beta,gamma){
  a = lik(i,sampH,delta,z,beta,gamma)
  b = as.numeric(delta*dlambda(sampH[,i],beta,gamma,z)/lambda(sampH[,i],beta,gamma,z)-lambda(sampH[,i],beta,gamma,z))
  c = b%*%t(b)
  d = as.numeric(delta*(ddlambda(sampH[,i],beta,gamma,z)*lambda(sampH[,i],beta,gamma,z)-dlambda(sampH[,i],beta,gamma,z)*dlambda(sampH[,i],beta,gamma,z))/(lambda(sampH[,i],beta,gamma,z)*lambda(sampH[,i],beta,gamma,z))-dlambda(sampH[,i],beta,gamma,z))
  e = d*diag(length(delta))
  return(a*(c+e))
}  


#sampling u according to ranks using progressive II censoring
spu = function(i,umat,n.cen,n,L){
  tempu = sort(umat[i,])
  rankumat = numeric(n)
  index = seq(1,n,1)
  
  i1=1
  i2=1+n.cen[1]
  if(n.cen[1] == 0){
    rankumat[i1] = tempu[1]
    index = index[-1]
  }
  if(n.cen[1] > 0){
    k = sample(index,n.cen[1])
    index = index[-k]
    rankumat[i1:(i2-1)] = 1.0e-8 ##should be 0, but numerically not feasible
    rankumat[i2] = tempu[index[1]]
    index = index[-1]
  } 
  j=2
  while(j<L){  
    i1 = i2 + 1
    i2 = i2 + n.cen[j] + 1
    if(n.cen[j] == 0){
      rankumat[i2] = tempu[index[1]]
      index = index[-1]
    }
    if(n.cen[j] > 0){
      rankumat[i1:(i2-1)] = rankumat[(i1-1)]
      k = sample(c(1:length(index)),n.cen[j])
      index = index[-k]
      rankumat[i2] = tempu[index[1]]
      index = index[-1]
    }
    j=j+1
  }
  if(n.cen[L] > 0){
    rankumat[(i2+1):n] = rankumat[i2]
  }
  return(rankumat)
}

maglik = function(M,sampH,delta,z,beta,gamma){  
  s = sapply(1:M,lik,sampH=sampH,delta=delta,z=z,beta=beta,gamma=gamma)
  return(mean(s))
}

magliknull = function(M,sampH,delta,z,beta,gamma){  
  s = sapply(1:M,liknull,sampH=sampH,delta=delta,z=z,beta=beta,gamma=gamma)
  return(mean(s))
}


gradmaglik = function(M,sampH,delta,z,beta,gamma){
  s = sapply(1:M,dlik,sampH=sampH,delta=delta,z=z,beta=beta,gamma=gamma)
  return(apply(s,1,mean))
}

hessmaglik = function(M,sampH,delta,z,beta,gamma){
  p = length(beta)
  #p = length(delta)
  s = sapply(1:M,ddlik,sampH=sampH,delta=delta,z=z,beta=beta,gamma=gamma)
  return(matrix(apply(s,1,mean),byrow=T,nrow=p))
}

wshoot <- function(n,p,x,y,init,weight,lambda,maxiter=40,tol=1e-7)
{
 #init=initial beta; weight=weight for each components

  Q = t(x)%*%x
  B = t(x)%*%y
  i=0
  status = 0

  lams =lambda*weight 
  oldbeta <- init
  tmpbeta <- oldbeta

  while (i<maxiter && status==0){ 
    for (j in 1:p){ 
         s<-ss2(j,tmpbeta,Q,B)
         if (s > lams[j]) 
            tmpbeta[j]<-(lams[j]-s)/(2*Q[j,j])
       else if (s < (-lams[j]) ) 
            tmpbeta[j]<-(-lams[j]-s)/(2*Q[j,j])
         else
        tmpbeta[j]<- 0.0
      }
    dx<-max(abs(tmpbeta-oldbeta))
        oldbeta <- tmpbeta
        if (dx<=tol) 
            status <- 1 
        i <- i+1
  }
  tmpbeta
}

ss2 <- function(j,tmpb,Q,B)
{
 a <- sum(tmpb*Q[,j])-tmpb[j]*Q[j,j]
 s <- 2*(a-B[j])
 return(s)
}


lse <- function(x,y)
{ Q <- t(x)%*%x
  P <- t(x)%*%y
  beta <- solve(Q,P)
  beta
}


ginv = function(X, tol = sqrt(.Machine$double.eps)){
  s = svd(X)
  nz = s$d > tol * s$d[1]
  if(any(nz)) s$v[,nz] %*% (t(s$u[,nz])/s$d[nz])
  else X*0
}

normalize = function(x){
  y = (x-mean(x))/sqrt(sum((x-mean(x))^2))
  return(y)
}

#time: the observed event time 
#delta: censoring indicator
#z: the regression predictors
#gamma: the parameter specifies the PH (gamma = 0) or PO (gamma = 1) model 
#lambda: the tunning parameter
#NN: number of iterations in Newto-Ralphson (default = 10)
#M: number of importance sampling used to compute the marginal likelihood (default = 1000)

alasso_odds <- function(NN = 10, M = 1000, gamma = 1, time, delta, z, lambda)
{
    iter = 100
    tol = 1.0e-10
    n = length(time)
    p = length(z[1,])
    
    #incomplete rank calculation
    n.fail = sum(delta)#n.fail denotes number of failures in n subjects
    L = n.fail+1#L denotes number of intervals between failures
    n.cen = numeric(L)#n.cen records the number of censored observations in each interval
    failrank = sort(rank(time)[delta==1])#rank of ordered failure time
    temprank1 = c(0,failrank)
    temprank2 = c(failrank,n+1) 
    n.cen = temprank2 - temprank1 - 1
  
    #resampling ranks based on progressive II censoring using uniform distribution
    uvect = runif(M*n)
    umat = matrix(uvect,nrow=M,byrow=T)
    sampu = sapply(1:M,spu,umat=umat,n.cen=n.cen,n=n,L=L)
    sampH = sapply(1:M,invF,u=sampu,gamma=gamma)
    
    true.sd = sqrt(apply(z,2,var)*(n-1))
    delta = delta[order(time)]
    ordz = z[order(time),] 
    z = apply(ordz,2,normalize)
    sd = sqrt(apply(z,2,var)*(n-1))
    #computing initial estimates
    ii = 0
    beta = numeric(p)
    while(ii < NN){
      fn=maglik(M,sampH,delta,z,beta,gamma)
      G=-gradmaglik(M,sampH,delta,z,beta,gamma)/fn
      H=-(hessmaglik(M,sampH,delta,z,beta,gamma)/fn - G%*%t(G))
     
      X = chol(H) 
      vecY = forwardsolve(t(X),H%*%beta-G)
      lsbeta = lm(vecY~-1+X)$coef
      beta1 = lsbeta
      dx = max(abs(beta1-beta))
      ii = ii + 1
      istop = ii
      if(dx <= 1.0e-5) ii = NN
      beta = beta1
    }
    inibeta = beta
    #computing adaptive-Lasso solutions
    sd = sd/abs(inibeta)
    ii = 0
    beta = numeric(p)
    while(ii < NN){
      fn=maglik(M,sampH,delta,z,beta,gamma)
      G=-gradmaglik(M,sampH,delta,z,beta,gamma)/fn
      H=-(hessmaglik(M,sampH,delta,z,beta,gamma)/fn - G%*%t(G))
     
      X = chol(H) 
      vecY = forwardsolve(t(X),H%*%beta-G)
      lsbeta = lm(vecY~-1+X)$coef
      beta1 = wshoot(p,p,X,vecY,init=lsbeta,sd,lambda,iter,tol)
      dx = max(abs(beta1-beta))
      ii = ii + 1
      istop = ii
      if(dx <= 1.0e-5) ii = NN
      beta = beta1
    }
    beta = beta/true.sd
    inibeta = inibeta/true.sd
    return(rbind(t(inibeta),t(beta)))
}