R/gurobimodelinit.R

In qhengncsu/BNRPxMCMC: Bayesian l1 Trend Filtering via Proximal MCMC

#' gurobi model initializer for proximal bayesian trendfiltering
#' @param x predictor (must be strictly increasing!)
#' @param w weight of x
#' @param restrictions a collection of shape restrictions,
#' including 'increasing','decreasing','convex','concave'
#' @param k order of piecewise polynomial; 1 linear, 2 quadratic, 3 cubic;k=0 and k>=4 not recommended
gurobimodelinit <- function(x,w,restrictions,k){
  n <- length(x)
  D <- getD(k,n,x)
  D1 <- getD(0,n,x)
  D2 <- getD(1,n,x)
  model <- list()
  model$Q <- Matrix::Diagonal(2*n-k-1,c(rep(0.5*w),rep(0,n-k-1)))
  model$A <- rbind(cbind(D,-Matrix::Diagonal(n-k-1)),
                   cbind(-D,-Matrix::Diagonal(n-k-1)))
  if('increasing' %in% restrictions){
    model$A <- rbind(model$A,cbind(-D1,matrix(0,n-1,n-k-1)))
  }
  if('decreasing' %in% restrictions){
    model$A <- rbind(model$A,cbind(D1,matrix(0,n-1,n-k-1)))
  }
  if('convex' %in% restrictions){
    model$A <- rbind(model$A,cbind(-D2,matrix(0,n-2,n-k-1)))
  }
  if('concave' %in% restrictions){
    model$A <- rbind(model$A,cbind(D2,matrix(0,n-2,n-k-1)))
  }
  model$sense <- rep('<',nrow(model$A))
  model$rhs <- rep(0,nrow(model$A))
  model$lb <- rep(-Inf,2*n-k-1)
  return(model)
}