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

require(CVXR); 

######################################################################
# min sum(rho_{M}(y - x * b)) + lambda * |D * b|_1
# subject to C * b >= d, E * b = f
# solve it using Quadratic Programming (CVXR)
# 
# Input:
#    x      --- n-by-p matrix, predictors 
#    y      --- n-by-1 vector, response
#    lambda --- scalar, nonnegative, tuning parameter in the penalty
#    M      --- scalar, tuning parameter in Huber loss, default = 1.35 
#    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 gurobi
# 
# Output: 
#    b      --- p-by-1 vector, estimated coefficients 
#    Qvalue --- scalar, value of the objective function 
#    cvx    --- list, results of LP problem by CVXR 
######################################################################

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

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

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

    bvar = Variable(p); 
    obj = 0.5 * sum(huber(y - x %*% bvar, M)) + lambda * p_norm(Dmat %*% bvar, 1); 

    ######################################################################
    # 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, M = M,
         Cmat = Cmat, dvec = dvec, Emat = Emat, fvec = fvec, Dmat = Dmat,
         b = b, Qvalue = result$value, cvx = result);
}

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