function [ sol_coef, time, best_eps, failure, Objective, Aineq, bineq ] = SQR_Decomposition_func( Xfull, Y, alpha)
%Authors:
%Dharmashankar Subramanian, last revised on 03/19/2014
%Pavithra Harsha, last revised on 03/19/2014

% This code is a method for CVaR regression.
% More specifically, it provides superquantile estimates at alpha for the response
% variable Y as a function of the covariates. Superquantile is well-known
% as the conditional value-at-risk (CVaR).

% This method in particular is the decomposition algorithm proposed in the
% paper by "Pavithra Harsha, Ramesh Natarajan and Dharmashankar Subramanian:
% A data-driven approach for the price-setting newsvendor problem"
% submitted in May 2015 and can be found on the authors' webpages.

% This MATLAB function uses CPLEX as the optimization solver.

% Do you want to scale or normalize the parameters? This is important for
% larger N and small alpha. This improves run times in those instances and
% avoids CPLEX from throwing a scaling related exception.
scaling = true;

%first column of Xfull is all ones
% Store in X all but the first column
X = Xfull(:,2:size(Xfull,2));

% dimension of the problem
% N = number of sample points
% P = number of covariates
N = size(X,1);
P = size(X,2);

N_Y = size(Y,1);
M_Y = size(Y,2);
if ( (N_Y ~= N) || (M_Y ~= 1) )
    error ('Dimensional mismatch in your input Y and X');
end


%beta: column vector
Valpha = ceil(N*alpha);
if(Valpha == N*alpha)
    beta = 1/N*(Valpha:1:N)';
else
    beta = [alpha; 1/N*(Valpha:1:N)'];%i = valpha -1, valpha, valpha +1,...,v-1,v
end

sb = size(beta,1);
A = log(ones(sb-2,1) - beta(1:sb-2))- log(ones(sb-2,1) - beta(2:sb-1));
C = beta(2:sb-1) - beta(1:sb-2);


num_UT = sb-2;
num_W = 1;
num_coef = P;

terminateWHILE = false;
initial_iter = 0;
iterate = 0;
coef_sol = zeros(P,1);
U_sol = zeros(num_UT,1);
T_sol = zeros(num_UT,1);
W_sol = 0;
coef_sol = zeros(P,1);
obj_sol = 100000;
epsilon = 0.000001;
if (scaling)
    sigma_y = std(Y);
    mu_y = mean(Y);
    sigma_x = std(X);
    mu_x = mean(X);
    X = (X - repmat(mean(X), size(X,1),1)).*repmat(1./std(X), size(X,1),1);
    Y = (Y - repmat(mean(Y), size(Y,1),1)).*repmat(1./std(Y), size(Y,1),1);
end
W_set = randperm(N);
initialized = false;
tic

%construct the objective
%     myObj = 1/(1-alpha)* C'*U + 1/(N*(1-alpha))* A'*T + W/(N*(1-alpha)) - 1/N* myUnity'*(Y - X*coef);
f_col = [(1/(1-alpha))*C; (1/(N*(1-alpha)))*A; 1/(N*(1-alpha)); (1/N)*(sum(X))'];

%construct Aineq
sp_T_nonneg = [sparse(num_UT,num_UT) -1*speye(num_UT) sparse(num_UT, num_W + num_coef)];
sp_T_init = [sparse(1:num_UT, 1:num_UT, -1*N*ones(num_UT,1)) -1*speye(num_UT) sparse(num_UT,1) -1*ones(num_UT,1)*sum(X)];
Aeq = [];
beq = [];
% if(Valpha == N*alpha)
%
% end
Aineq = [sp_T_nonneg; sp_T_init];

bineq = sparse(num_UT,1);
bineq = [bineq; -1*ones(num_UT,1)*sum(Y)];


while (~initialized)
    initial_iter = initial_iter +1;
    delta_Aineq = [sparse(1,num_UT) sparse(1,num_UT) -1 -1*X(W_set(initial_iter),:)];
    delta_bineq = -1*Y(W_set(initial_iter));
    
    Aineq = [Aineq;delta_Aineq];
    bineq = [bineq; delta_bineq];
    
    try
        [mysol, fval, exitflag, output, lambda] = cplexlp(f_col',Aineq,bineq,Aeq,beq,[], []);
        if (exitflag ~= 1 || output.cplexstatus ~= 1)
            %warning(strcat('Not Initialized yet, after iter :', num2str(initial_iter)))
        else
            %warning(strcat('INITIALIZED SUCCESSFULLY, after iter :', num2str(initial_iter)))
            %prepare outputs
            
            U_sol = mysol(1:num_UT,1);
            T_sol = mysol(num_UT+1:2*num_UT,1);
            W_sol = mysol(2*num_UT+1,1);
            coef_sol = mysol(2*num_UT+2:2*num_UT+1+P,1);
            obj_sol = fval;
            initialized = true;
            clear mysol;
            clear lambda;
        end
        
        
    catch m
        disp(m.message);
    end
    
    if ((initial_iter == N) && (~initialized))
        error('Should not reach here as all W constraints are added')
    end
end

%initial_iter
obj_prev = obj_sol;
bool_break = false;

opt = cplexoptimset('cplex');

%opt.read.scale = 1;

clear Xfull;

best_eps = 0;
failure = 1;
iterate = 0;

while (~terminateWHILE)
    iterate = iterate +1;
    %disp(strcat('IN DECOMP: ITERATE : ', num2str(iterate)))
    
    Bool_mat = sparse(( (Y-X*coef_sol)*ones(1,num_UT) - ones(N,1)*U_sol') > epsilon );
    delta_Aineq = [sparse(1:num_UT, 1:num_UT, -1*sum(Bool_mat)) -1*speye(num_UT) sparse(num_UT,1) -1*sparse(Bool_mat')*X];
    delta_bineq = -1*sparse(Bool_mat')*Y;
    
    Aineq = [Aineq;delta_Aineq];
    bineq = [bineq; delta_bineq];
    %
    [max_W, index_W] = max(Y-X*coef_sol);
    delta_Aineq = [sparse(1,num_UT) sparse(1,num_UT) -1 -1*X(index_W,:)];
    delta_bineq = -1*Y(index_W);
    
    Aineq = [Aineq;delta_Aineq];
    bineq = [bineq; delta_bineq];
    
    try
        [mysol, fval, exitflag, output, lambda] = cplexlp(f_col',Aineq,bineq,Aeq,beq,[], [], opt);
        if (exitflag ~= 1 || output.cplexstatus ~= 1)
            disp(strcat('EXITFLAG ROGUE : ', num2str(exitflag)));
            disp(strcat('OUTPUT.CPLEXSTATUS ROGUE: ', num2str(output.cplexstatus)));
            disp(strcat('OUTPUT.CPLEXSTATUSSTRING ROGUE: ', output.cplexstatusstring));
            
            U_sol = mysol(1:num_UT,1);
            T_sol = mysol(num_UT+1:2*num_UT,1);
            W_sol = mysol(2*num_UT+1,1);
            coef_sol = mysol(2*num_UT+2:2*num_UT+1+P,1);
            obj_sol = fval;
            if(iterate > 49)
                bool_break= true;
            end
            
            
        else
            %warning(strcat('INITIALIZED SUCCESSFULLY, after iter :', num2str(initial_iter)))
            %prepare outputs
            
            U_sol = mysol(1:num_UT,1);
            T_sol = mysol(num_UT+1:2*num_UT,1);
            W_sol = mysol(2*num_UT+1,1);
            coef_sol = mysol(2*num_UT+2:2*num_UT+1+P,1);
            obj_sol = fval;
            clear lambda;
            clear mysol;
            if(iterate > 49)
                best_eps = abs(obj_prev - obj_sol);
                failure = -9999;
                terminateWHILE = true;
                disp(strcat('TERMINATING.... CONVERGED! FAILURE: ', num2str(failure), ' BEST EPS : ', num2str(best_eps)))
            end
            if ( (abs(obj_prev - obj_sol)<epsilon) )
                best_eps = abs(obj_prev - obj_sol);
                
                %disp(strcat('EXITFLAG : ', num2str(exitflag)))
                %disp(strcat('OUTPUT.CPLEXSTATUS : ', num2str(output.cplexstatus)))
                %disp(strcat('OUTPUT.CPLEXSTATUSSTRING : ', output.cplexstatusstring))
                %disp(strcat('TERMINATING.... CONVERGED! AFTER ITER = ', num2str(iterate)))
                terminateWHILE = true;
                
                disp(strcat('TERMINATING.... CONVERGED! FAILURE: ', num2str(failure), ' BEST EPS : ', num2str(best_eps)))
            else
                obj_prev = obj_sol;
            end
            
        end
        
        
    catch m
        disp(m.message);
    end
    if (bool_break)
        best_eps = abs(obj_prev - obj_sol);
        failure = -9999;
        terminateWHILE = true;
        disp(strcat('TERMINATING.... CONVERGED! FAILURE: ', num2str(failure), ' BEST EPS : ', num2str(best_eps)))
        break
    end
    
end

time = toc;

sol_coef=-9999;
Objective = -9999;
if(failure > 0)
    Objective = obj_sol;
    residuals = Y - X*coef_sol;
    Sort_Residuals = sort(residuals);
    index = ceil(alpha*N);
    lambda = (index/N - alpha)/(1-alpha);
    if (index > 0)
        AVaRempiricalDist = (lambda*Sort_Residuals(index) + (1-lambda)/(N - index)*sum(Sort_Residuals(index+1:N)));
    else
        AVaRempiricalDist = (1-lambda)/(1- index/N)*(1/N)*sum(Sort_Residuals(index+1:N));
    end
    
    a_0 = AVaRempiricalDist;
    a_1 = coef_sol;
    sol_coef= [a_0;a_1];
    if (scaling)
        t_0 = a_0*sigma_y -(mu_x./sigma_x)*a_1*sigma_y+mu_y;
        t_1 = a_1*sigma_y./sigma_x';
        sol_coef= [t_0;t_1];
    end
    
end
end