######################################################################
# Quantile regression with generalized lasso penalty and linear constraints
# Peng Zeng @ Auburn University 
# 08-17-2021
######################################################################

require(CVXR); 

######################################################################
# min sum(rho_{tau}(y - x * b)) + lambda * |D * b|_1
# subject to C * b >= d, E * b = f
# solve it using Linear Programming by CVXR
#
# Input:
#    x      --- n-by-p matrix, predictors 
#    y      --- n-by-1 vector, response
#    lambda --- scalar, nonnegative, tuning parameter in the penalty
#    tau    --- scalar, value in (0, 1), median corresponds to 0.5 
#    Cmat   --- q-by-p matrix, inequality constraints  
#    dvec   --- q-by-1 vector, inequality constraints 
#    Emat   --- s-by-p matrix, equality constraints  
#    fvec   --- s-by-1 vector, equality constraints 
#    Dmat   --- m-by-p matrix, default: identity matrix 
#    ...    --- more parameters for solve() in CVXR
# 
# Output: 
#    b      --- p-by-1 vector, estimated coefficients 
#    Qvalue --- scalar, value of the objective function 
#    cvx    --- list, results of LP problem by CVXR 
######################################################################

quantreg.genlasso.fit.cvxr = function(
    x, y, lambda = 0.0, tau = 0.5, 
    Cmat = NULL, dvec = NULL, Emat = NULL, fvec = NULL, 
    Dmat = diag(rep(1.0, NCOL(x))), ...)
{
    stopifnot(NROW(x) == length(y)); 
    stopifnot(lambda >= 0.0, tau > 0.0, tau < 1.0); 
    stopifnot(NCOL(Dmat) == NCOL(x)); 

    n = NROW(x); p = NCOL(x); m = NROW(Dmat); 

    ######################################################################
    # specify model
    ######################################################################

    bvar = Variable(p);     
    obj = sum(scalene(y - x %*% bvar, tau, 1 - tau)) + lambda * norm1(Dmat %*% bvar); 

    ######################################################################
    # define contraints
    ######################################################################

    constr = list(); 
    if(!is.null(Cmat))
    {
        if(is.null(dvec)) dvec = rep(0.0, NROW(Cmat));
        stopifnot(NROW(Cmat) == length(dvec), NCOL(Cmat) == NCOL(x));
        
        constr = c(constr, list(Cmat %*% bvar >= dvec)); 
    }

    if(!is.null(Emat))
    {
        if(is.null(fvec)) fvec = rep(0.0, NROW(Emat));
        stopifnot(NROW(Emat) == length(fvec), NCOL(Emat) == NCOL(x));

        constr = c(constr, list(Emat %*% bvar == fvec)); 
    }

    prob = Problem(Minimize(obj), constr); 
    result = solve(prob, ...); 

    b = drop( result$getValue(bvar) ); 

    list(x = x, y = y, lambda = lambda, tau = tau,
         Cmat = Cmat, dvec = dvec, Emat = Emat, fvec = fvec, Dmat = Dmat,
         b = b, Qvalue = result$value, cvx = result);
}

######################################################################
# THE END
######################################################################