#Required packages
# install.package(lpSolve)
# install.package(ggplot2)
# install.package(directlabels)

#Inputs
###  sc=selection criteria, AIC=1, BIC=2
###  gamma=threshold
###  X=design matrix
###  Y=response vector

#Outputs

## model.terms=terms selected based on value of selection criteria (sc)
## estimates=full matrix of coefficient estimates at each value of delta
## plot=profile plot of each variable's estimate vs. delta

##   Dantzig Selector variable selection   ##
dantzig=function(X,Y,delta,gamma,sc, vlength){
  k=dim(X)[2]
  c1=matrix(1,k,1)
  c2=matrix(0,k,1)
  c=rbind(c1,c2)
  A11=t(X)%*%X
  A12=-t(X)%*%X
  A1=cbind(A11,A12)
  A21=-t(X)%*%X
  A22=t(X)%*%X
  A2=cbind(A21,A22)
  A31=2*diag(k)
  A32=-diag(k)
  A3=cbind(A31,A32)
  A=rbind(A1,A2,A3)
  b1=-t(X)%*%Y-delta*matrix(1,k,1)
  b2=t(X)%*%Y-delta*matrix(1,k,1)
  b3=matrix(0,k,1)
  b=rbind(b1,b2,b3)
  dir=matrix(">=",dim(b)[1],1)
  
  vars=lp("min", c, A, dir, b)$solution
  vars=as.vector(vars)
  u=vars[1:k]
  beta=(vars[(k+1):length(vars)]-u)/vlength
  terms=(abs(beta)>gamma)
  answer=list(sum(terms),terms, beta, delta)
  return(answer)
  
  
}



dantzig.selection <- function(X,Y, sc)  {
  require(lpSolve)
  require(ggplot2)
  require(directlabels)
  Yc=Y-mean(Y)
  Xc<-apply(X, 2, function(X) X - mean(X))
  vlength<-sqrt(diag(t(Xc)%*%Xc))
  Xs<-t(t(Xc)*(1/vlength))
  n=length(Y)
  i=0
  d0=round(max(abs((t(Xs)%*%Yc))))
  out=dantzig(Xs,Yc,delta=0,gamma=0, sc=2, vlength)
  gamma<-0.1*max(abs(out[[3]]))
  output=matrix(0,nrow=d0+1,ncol=3)
  est=matrix(0, nrow=d0+1, ncol=ncol(X))

    for (delta in seq(0,d0,by=1)){
    i=i+1
    out=dantzig(Xs,Yc,delta,gamma, sc=2, vlength)
    beta=out[[3]]
    beta=as.vector(beta)
    p=out[[1]]+1
    if (sum(out[[2]])==0) Xmodel=rep(1,dim(X)[1])
    else Xmodel=X[,out[[2]]]
    lin=lm(Y~Xmodel)
    SSE=deviance(lin)
    if (sc==1){
      crit=n*log(SSE/n)+(2*p)+((2*p*(p+1))/(n-p-1)) #AICc
    } else {
      crit<-n*log(SSE/n)+p*log(n) #BIC
    }

    est[i,]=beta
    output[i,1]=delta
    output[i,2]=p
    output[i,3]=crit
  
    
  }
crit.best=min(output[,3])
delta.best=which(output[,3]==min(output[,3]))[1]-1
model.terms=dantzig(Xs,Yc,delta.best,gamma, sc=2, vlength)
delta<-seq(1:nrow(est))
df<-as.data.frame((est))
label=colnames(df)
k<-dim(df)[2]
df<-stack(df)
df$delta<-rep(delta, k)
p<-ggplot(data=df, aes(y=values, x=delta, group=ind, label=ind, color=ind))+geom_line(size=0.75)+
  theme_bw()+ggtitle("Dantzig Coefficient Estimates")+labs(title="Scaled Dantzig Coefficient Estimates", x="\u03B4", y="Estimates")
p<-direct.label(p, list("first.points", cex=0.7))
answer=list(terms=model.terms, estimates=est, plot=p, gamma=gamma)
return(answer)

}
