function [ Sol, ind_zf, tsol ] = DPC_MTFL_ZWZ( Xs, ys, Lambda, opts )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Changed from the function DPC_MTFL, the difference is that the input Xs and ys are cells instead of matrix and vector

%% Implementation of the sequential DPC rule proposed in

%  Safe Screening for Multi-Task Learning with Multiple Data Matrices

%  Jie Wang, Peter Wonka, and Jieping Ye. 

%

%% input: 

%         X: 

%            stacked data matrix, each column corresponds to a feature 

%            each row corresponds to a data instance

%

%         y: 

%            the response vector

%       

%         Lambda: 

%            the parameters sequence

%

%         opts: 

%            settings for the solver

%% output:

%         Sol: 

%              the solution; the ith column corresponds to the solution

%              with the ith parameter in Lambda

%

%         ind_zf: 

%              the index of the discarded features; the ith column

%              refers to the solution of the ith value in Lambda

%

%         t_sol: 

%              the ith entry of ts_sol records the running time of the

%              solver after screening

%% For any problem, please contact Jie Wang (jiewangustc@gmail.com)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X = [];y = [];

for idx_t = 1:length(ys)

X = [X;Xs{idx_t}];

Xs{idx_t} = [];

y = [y;ys{idx_t}];

ys{idx_t}=[];

end

[N,d] = size(X);

npar = length(Lambda); % number of parameters

% ---------------------------- pass parameters -------------------------- %

opts.init = 0; % set .init for warm start

task_ind = opts.ind;

stask = diff(task_ind); % size of each task

T = length(stask); % number of tasks

% --------------------- initialize the output --------------------------- %

Sol = zeros(d,T,npar);

ind_zf = false(d,npar);

%rej_ratio = zeros(1,npar);

% ------- construct sparse matrix to vectorize the computation ---------- %

ft_ind = zeros(1,N);

for i = 1:T

    ft_ind(1,task_ind(i)+1:task_ind(i+1))=i; % fg_ind(j) is the index of the group 

                                     % which contains the jth feature

end

tS = sparse(ft_ind,1:N,ones(1,N),T,N,N); % ith row refers to the ith task, if the ith row

                   % refers to the i_s to i_e features, then

                   % tS(i,i_s:i_e) = 1.

% ------ compute the norm of each feature of each task ||x_l^t||_2 ---------- %  

Xtnorm = sqrt(tS*(X.*X));

H = -2*Xtnorm;

[xtnmx,xtnmxind] = max(Xtnorm);

% --------------------------- compute lambda_max ------------------------ %

Xtty = tS*(X.*repmat(y,1,d)); % Xtty is of T*d and Xtty(i,j) = <x_j^i,y_i>

[lambda_max2, indmx] = max(sum(Xtty.*Xtty,1));

lambda_max = sqrt(lambda_max2);

opts.lambda_max = lambda_max;

if opts.rFlag == 1

    opts.rFlag = 0; % the parameter value passing to the solver is its 

                       % absolute value rather than a ratio

    Lambda = Lambda*lambda_max;

end

% ----------------- sort the parameters in descend order ---------------- %

[Lambdav,Lambda_ind] = sort(Lambda,'descend');

rLambdav = 1./Lambdav;

% ----------------- solve MTFL sequentially with DPC ------------ %

lambdap = lambda_max;

tsol = zeros(1,npar);

for i = 1:npar

    fprintf('in DPC step: %d\n',i);

    lambdac = Lambdav(1,i);

    if lambdac>=lambda_max

        %Sol(:,Lambda_ind(i)) = 0; % because Sol is initalized to be 0

        ind_zf(:,Lambda_ind(i)) = true;     

        %rej_ratio(Lambda_ind(i)) = 1;

    else

        if lambdap==lambda_max

            theta = y/lambda_max;

            indmx = indmx(1);

            nv = (tS'*Xtty(:,indmx)).*X(:,indmx);

        else

            theta = (y - sum(X.*(Sol(:,:,Lambda_ind(i-1))*tS)',2))*rlambdap;

            nv = y*rlambdap-theta;

        end

        rlambdac = rLambdav(1,i);

        nv = nv/norm(nv);

        rv = y*rlambdac-theta;

        Prv = rv - (nv'*rv)*nv;

        o = theta + 0.5*Prv;

        r = 0.5*norm(Prv);

        % ------- screening by DPC, remove the ith feature if T(i)=1 ----- %

        Xtto = tS*(X.*repmat(o,1,d));

        q = -2*(Xtto.*Xtnorm);

        qpmx = -newton_qp(H,q,r,2*xtnmx,xtnmxind);

        Tmt = 1 > (sum(Xtto.*Xtto,1))' + qpmx + 1e-8;      

        ind_zf(:,Lambda_ind(i)) = Tmt;

        nTmt = ~Tmt;

        Xr = X(:,nTmt);

        if lambdap == lambda_max

            opts.x0 = zeros(size(Xr,2),T);

        else

            opts.x0 = Sol(nTmt,:,Lambda_ind(i-1));

        end

        % --- solve the group Lasso problem on the reduced data matrix -- %

        starts = tic;

        [x1, ~, ~]= mtLeastR(Xr, y, lambdac, opts);   

        tsol(Lambda_ind(i)) = toc(starts);

        Sol(nTmt,:,Lambda_ind(i)) = x1;

        lambdap = lambdac;

        rlambdap = rlambdac;

    end

end

end