function [Sol_MT] = STM_HFS(Xs, ys, Lambda, opts)

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

%% Implementation of the sequential HFS rule for STM

%% input:

%         Xs:

%            Xs{i} stores the data matrix of the i-th task, each column corresponds to a feature

%            each row corresponds to a data instance

%

%         ys:

%            ys{i} strores the response vector of the i-th task

%

%         Lambda:

%            the parameter values of lambda

%

%         opts:

%            settings for the solver

%% output:

%         Sol:

%              the solution; Sol(:,:,i) stores the the

%              solution for the ith values in Lambda

%

%% For any problem, please contact Weizhong Zhang (zhangweizhongzju@gmail.com)

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

%%

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

p = size(Xs{1}, 2);

T_num = length(Xs); % number of tasks

npar = length(Lambda); % number of parameter values of lambda

ind = opts.ind;

ind_MT = TreeTransform(ind, T_num);

%clear('ind');

opts.init=1;        % starting from a zero point

if opts.tFlag==2

    funVal = opts.funVal;

end

% --------------------- recover tree structure for Multi Task-------------------------- %

eind_MT = find(ind_MT(2,:) == p*T_num);

if ind_MT(1,1) == -1 % find the depth of the tree

    d_MT = length(eind_MT);

    nnl_MT = [p*T_num,eind_MT(1)-1,diff(eind_MT)]; % number of nodes per layer

    nnind_MT = [0,1,eind_MT]; % ind(1,nnind1(i)+1:nnind(i+1)) stores the node from the ith layer

else

    d_MT = length(eind_MT)-1;

    nnl_MT = [eind_MT(1),diff(eind_MT)];

    nnind_MT = [0,eind_MT];

end

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

Sol_MT = zeros(p,T_num,npar);

ind_zf_MT = false(p*T_num,d_MT,npar);

tsolver_MT = zeros(1,npar);%should be changed

tscreen_MT = zeros(1,npar);

% ------------------- compute the effective region of lambda ------------ %

Xsys_MT = Xsys_MT_cal(Xs, ys);

%lambda_max=findLambdaMax(X'*y, p, ind, size(ind,2)); % why?

lambda_max=findLambdaMax(Xsys_MT, p*T_num, ind_MT, size(ind_MT,2));

%lambda_max1=findLambdaMax(Xsys_MT, p*T_num, ind, size(ind,2));

if opts.rFlag == 1

    Lambda = Lambda * lambda_max;

    opts.rFlag = 0;

end

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

% --------- compute the norm of each feature and each submatrix ---------

Xnorm_MT = zeros(size(Xs{1},2)*T_num,1);

idx_row = [0:size(Xs{1},2)-1]*T_num+1;

for idx_t = 1:T_num

    Xnorm_MT(idx_row) = (sqrt(sum(Xs{idx_t}.^2,1)))';

    idx_row = idx_row +1;

end

ng_MT = size(ind_MT,2);

Xgnorm_MT = zeros(1,ng_MT);

gind_MT = zeros(d_MT,p*T_num);

if ind_MT(1,1)==-1

    j = 2;

    l = 2;

else

    l = 1;

    j = 1;

end

k = 1;

for i = j : ng_MT-1

    for idx_t = 1: length(Xs)

        if(idx_t ==1)

            Xgnorm_MT(i) = norm(Xs{idx_t}(:,((ind_MT(1,i)-1)/T_num+1):(ind_MT(2,i)/T_num)));

        else

            Xgnorm_MT(i)= max(Xgnorm_MT(i),norm(Xs{idx_t}(:,((ind_MT(1,i)-1)/T_num+1):(ind_MT(2,i)/T_num))));

        end

    end

    gind_MT(l,ind_MT(1,i):ind_MT(2,i)) = k;

    k = k + 1;

    if ind_MT(2,i)==p*T_num

        k = 1;

        l = l + 1;

    end

end

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

Gind_MT = cell(1,d_MT);

if ind_MT(1,1) == -1

    j = 2;

else

    j = 1;

end

for i = j:d_MT

    Gind_MT{1,i} = sparse(gind_MT(i,:),1:p*T_num,ones(1,p*T_num),nnl_MT(i),p*T_num);

end

% --------------- put Xgnorm and weights in tree structure ---------------

XgnormTree_MT = zeros(p*T_num,d_MT);

weightTree_MT = zeros(p*T_num,d_MT);

for i = 1:d_MT

    if ind_MT(1,1)==-1&&i==1

        XgnormTree_MT(:,i) = Xnorm_MT;

        weightTree_MT(:,i) = ind_MT(3,1);

    else

        G_MT = Gind_MT{1,i};

        XgnormTree_MT(:,i) = G_MT'*(Xgnorm_MT(nnind_MT(i)+1:nnind_MT(i+1)))';

        weightTree_MT(:,i) = G_MT'*(ind_MT(3,nnind_MT(i)+1:nnind_MT(i+1)))';

    end

end

% ----------- solve STM sequentially via HFS ------------------

opts.rFlag = 0; % the input parameters are their true values

s_MT = zeros(p*T_num,1);

c2_MT = zeros(p*T_num,1);

minn_MT = zeros(p*T_num,1);

rLambdav = 1./Lambdav;

lambdap = Lambdav(1);

rlambdap = rLambdav(1);

vnormTree_MT = zeros(p*T_num,d_MT);

tol0 = 1e-12;

for i = 1:npar

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

    lambdac = Lambdav(1,i);

    rlambdac = rLambdav(1,i);

    if lambdac>=lambda_max

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

    else

        starts_screening =tic;

        if lambdap==lambda_max

            theta_MT = [];

            for idx_t = 1: length(ys)

                theta_MT = [theta_MT; ys{idx_t}*rlambdap];

            end

            z_MT = Xsys_MT*rlambdap;

            [u_MT, v_MT] = Hierarchical_Projection( z_MT, ind_MT, nnind_MT, Gind_MT );

            if ind_MT(3,end)==0

                weightd_MT = (ind_MT(3,nnind_MT(d_MT)+1:nnind_MT(d_MT+1)))';

                [~,Xmxind_MT] = min(abs(Gind_MT{1,d_MT}*(v_MT(:,d_MT).*v_MT(:,d_MT))-weightd_MT.*weightd_MT));

                idx_column_MT  = [ind_MT(1,nnind_MT(d_MT)+Xmxind_MT):ind_MT(2,nnind_MT(d_MT)+Xmxind_MT)];

                nv_MT = [];

                idx_column_X = [(idx_column_MT(end)/T_num)-(length(idx_column_MT)/T_num)+1:(idx_column_MT(end)/T_num)];

                for idx_t = 1:length(ys)

                    idx_column = idx_column_MT(idx_t:T_num:end);

                    nv_MT = [nv_MT; Xs{idx_t}(:,idx_column_X)*v_MT(idx_column,d_MT)];

                end

            else

                nv_MT = [];

                for idx_t = 1:length(ys)

                    idx_column = [idx_t:T_num:lenght(v_MT)];

                    nv_MT = [nv_MT; Xs{idx_t}*v_MT(idx_column,end)];

                end

            end

        else

            theta_MT = [];

            y_all = [];

            for idx_t = 1:length(ys)

                theta_MT = [theta_MT;(ys{idx_t} - Xs{idx_t}*Sol_MT(:,idx_t,Lambda_ind(i-1)))*rlambdap];

                y_all = [y_all; ys{idx_t}];

            end

            nv_MT = y_all*rlambdap-theta_MT;

        end

        % ----- estimate the possible region of the dual optimum at lambdac

        nv_MT = nv_MT/norm(nv_MT);

        %rv = y*rlambdac-theta;

        rv_MT = [];

        for idx_t = 1:length(ys)

            rv_MT = [rv_MT;ys{idx_t}*rlambdac];

        end

        rv_MT = rv_MT-theta_MT;

        Prv_MT = rv_MT - (nv_MT'*rv_MT)*nv_MT;

        o_MT = theta_MT + 0.5*Prv_MT;

        r_MT = 0.5*norm(Prv_MT);

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

        c_MT = zeros(p*T_num,1);

        idx_row = (0:size(Xs{1},2)-1)*length(Xs)+1;

        for idx_t = 1:length(ys)

            c_MT(idx_row) = Xs{idx_t}'*o_MT((idx_t-1)*length(ys{1})+1:idx_t*length(ys{1}));

            idx_row = idx_row+1;

        end

        [u_MT, v_MT] = Hierarchical_Projection( c_MT, ind_MT, nnind_MT, Gind_MT );

        v2_MT = v_MT.*v_MT;

        for l = 1:d_MT % compute norm of v for each node and arrange them based on the tree structure

            if l==1&&ind_MT(1,1)==-1

                vnormTree_MT(:,l) = abs(v_MT(:,l));

            else

                G_MT = Gind_MT{1,l};

                vnormTree_MT(:,l)=G_MT'*sqrt(G_MT*v2_MT(:,l));

            end

        end

        csDifwv_MT = cumsum(weightTree_MT-vnormTree_MT);

        T_MT = false(p*T_num,1); % identify non-leaf inactive nodes

        for l = d_MT:-1:2

            Tl_MT = ~T_MT; % find the indices of the remaining features

            % case 1

            Tc_MT = false(p*T_num,1);

            Tc_MT(Tl_MT) = vnormTree_MT(Tl_MT,l)>tol0;

            if nnz(Tc_MT)>0

                s_MT(Tc_MT) = vnormTree_MT(Tc_MT,l)+r_MT*XgnormTree_MT(Tc_MT,l);

            end

            % case 2 & 3

            if nnz(Tc_MT)<nnz(Tl_MT) % if not all remaining nodes in level l fall in case 1

                Tcc_MT = false(p*T_num,1);

                Tcc_MT(Tl_MT) = ~Tc_MT(Tl_MT);

                lind_MT = nnind_MT(l)+1:nnind_MT(l+1);

                G_MT = Gind_MT{1,l};

                Tn_MT = G_MT*Tcc_MT==(ind_MT(2,lind_MT)-ind_MT(1,lind_MT)+1)';

                indl_MT = ind_MT(:,lind_MT);

                indlr_MT = indl_MT(:,Tn_MT);

                for n = 1:nnz(Tn_MT)

                    minn_MT(n)=min(csDifwv_MT(indlr_MT(1,n):indlr_MT(2,n),l-1));

                end

                sdist_MT = (G_MT(Tn_MT,:))'*minn_MT(1:nnz(Tn_MT));

                s_MT(Tcc_MT) = max(0,r_MT*XgnormTree_MT(Tcc_MT,l)-sdist_MT(Tcc_MT));

            end

            ind_zf_MT(Tl_MT,l,Lambda_ind(i))=s_MT(Tl_MT)<weightTree_MT(Tl_MT,l);

            T_MT = T_MT|ind_zf_MT(:,l,Lambda_ind(i));

        end

        Tl_MT = ~T_MT; % identify inactive leaf nodes

        if ind_MT(1,1)==-1

            s_MT(Tl_MT) = abs(c_MT(Tl_MT))+r_MT*Xnorm_MT(Tl_MT);

            ind_zf_MT(Tl_MT,1,Lambda_ind(i))=s_MT(Tl_MT)<ind_MT(3,1);

        else

            G_MT = Gind_MT{1,1};

            lind_MT = nnind_MT(1)+1:nnind_MT(2);

            Tn_MT = G_MT*Tl_MT == (ind_MT(2,lind_MT)-ind_MT(1,lind_MT)+1)';

            c2_MT(Tl_MT) = c_MT(Tl_MT).*c_MT(Tl_MT);

            cnorm_MT = (G_MT(Tn_MT,:))'*sqrt(G_MT(Tn_MT,:)*c2_MT);

            s_MT(Tl_MT) = cnorm_MT(Tl_MT)+r_MT*XgnormTree_MT(Tl_MT,1);

            ind_zf_MT(Tl_MT,1,Lambda_ind(i))=s_MT(Tl_MT)<weightTree_MT(Tl_MT,1);

        end

        T_MT = T_MT|ind_zf_MT(:,1,Lambda_ind(i));

        nT_MT = ~T_MT;

        %Xr = X(:,nT);

        nT_MT_X = nT_MT(1:T_num:end);

        for idx_t = 1:T_num

            Xrs{idx_t} = Xs{idx_t}(:,nT_MT_X);

        end

        if lambdap == lambda_max

            opts.x0 = zeros(nnz(nT_MT_X)*T_num,1);

        else

            x0_Matrix = Sol_MT(nT_MT_X,:,Lambda_ind(i-1));

            x0_temp =zeros(size(x0_Matrix,1)*T_num,1);

            idx_row = [0:size(x0_Matrix,1)-1]*T_num+1;

            for idx_t = 1:T_num

                x0_temp(idx_row) = x0_Matrix(:,idx_t);

                idx_row = idx_row+1;

            end

            opts.x0 = x0_temp;

        end

        % ------------ construct the reduced tree ---------------

        Tind_MT = false(ng_MT,1);

        nnlr_MT = zeros(1,d_MT+1);

        nnlr_MT(end) = 1;

        if ind_MT(1,1)==-1

            j=2;

            nnlr_MT(1)=nnz(nT_MT);

        else

            j=1;

        end

        for l = j:d_MT

            lind_MT = nnind_MT(l)+1:nnind_MT(l+1);

            Tind_MT(lind_MT) = Gind_MT{1,l}*T_MT==(ind_MT(2,lind_MT)-ind_MT(1,lind_MT)+1)';

            nnlr_MT(l)=nnz(~Tind_MT(lind_MT));

        end

        if ind_MT(1,1)==-1

            nnindr_MT=[0,1,cumsum(nnlr_MT(2:end))+1];

        else

            nnindr_MT=[0,cumsum(nnlr_MT)];

        end

        indr_MT = ind_MT(:,~Tind_MT);

        mapinde_MT = cumsum(nT_MT);

        mapinds_MT = nnz(nT_MT)+1-cumsum(nT_MT,'reverse');

        for l=j:d_MT+1

            lind_MT = nnindr_MT(l)+1:nnindr_MT(l+1);

            oind1_MT = indr_MT(1,lind_MT);

            oind2_MT = indr_MT(2,lind_MT);

            indr_MT(1,lind_MT) = mapinds_MT(oind1_MT);

            indr_MT(2,lind_MT) = mapinde_MT(oind2_MT);

        end

        opts.ind_MT = indr_MT;

        % --- solve the STM problem on the reduced data matrix -- %

        if opts.tFlag == 2

            opts.tol = funVal(Lambda_ind(i));

        end

        tscreen_MT(Lambda_ind(i)) = toc(starts_screening);

        starts = tic;

        [x1, ~, ~]= tree_LeastR_MT(Xrs, ys, lambdac, opts);

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

        nT_MT_temp = nT_MT(1:T_num:end);

        idx_row= [1:T_num:length(x1)];

        for idx_t = 1:T_num

            Sol_MT(nT_MT_temp,idx_t,Lambda_ind(i)) = x1(idx_row);

            idx_row = idx_row +1;

        end

    end

    lambdap = lambdac;

    rlambdap = rlambdac;

end

end