function [x_MT, funVal, ValueL_MT]=tree_LeastR_MT(As, ys, z, opts)
%
%%
% Function tree_LeastR
%      Least Squares Loss with the
%           tree structured group Lasso Regularization in Multi-Task
%           Learning
%
%% Problem
%
%  min  1/2 sum_t || A_t x_t - y_t||^2 + z * sum_j w_j ||x_{G_j}||
%
%  G_j's are nodes with tree structure
%
%  The tree overlapping group information is contained in
%  opts.ind_MT, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind_MT(1,:) contains the starting index
%  opts.ind_MT(2,:) contains the ending index
%  opts.ind_MT(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x( opts.ind_MT(1,j):opts.ind_MT(2,j) ) denotes x_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and n.
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x( opts.G_MT ( opts.ind_MT(1,j):opts.ind_MT(2,j) ) ) denotes x_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind_MT(1:2,:) are within 1 and length(opts.G_MT).
%
% The following example shows how G and ind works:
%
% G={ {1, 2}, {4, 5}, {3, 6}, {7, 8},
%     {1, 2, 3, 6}, {4, 5, 7, 8},
%     {1, 2, 3, 4, 5, 6, 7, 8} }.
%
% ind={ [1, 2, 100]', [3, 4, 100]', [5, 6, 100]', [7, 8, 100]',
%       [9, 12, 100]', [13, 16, 100]', [17, 24, 100]' },
%
% where each node has a weight of 100.
%
%
%% Input parameters:
%
%  As-         A set of Matrix of size m_t x n
%                At can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  ys -        A set of Response vector (of size m_tx1)
%  z -        Tree Structure regularization parameter (z >=0)
%  opts-      optional inputs (default value: opts=[])
%
%% Output parameters:
%
%  xs-         A set of solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2010-2011 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on October 3, 2010.
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%% Related functions:
%
%  sll_opts
%
%%

%% Verify and initialize the parameters
%%
if (nargin <4)
    error('\n Inputs: As, ys, z, and opts.ind_MT should be specified!\n');
end

[m_MT,n_MT]=size(As{1});
T_num = length(ys);

if (length(ys{1}) ~=m_MT)
    error('\n Check the length of ys!\n');
end

if (z<0)
    error('\n z should be nonnegative!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags) 

% restart the program for better efficiency
%  this is a newly added function
if (~isfield(opts,'rStartNum'))
    opts.rStartNum=opts.maxIter;
else
    if (opts.rStartNum<=0)
        opts.rStartNum=opts.maxIter;
    end
end

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%



if (~issparse(As{1})) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & Others

% Initialize ind
if (~isfield(opts,'ind_MT'))
    error('\n In tree_LeastR, the field .ind should be specified');
else
    ind_MT=opts.ind_MT;
    
    if (size(ind_MT,1)~=3)
        error('\n Check opts.ind');
    end
end

GFlag_MT=0;
% if GFlag=1, we shall apply general_altra
if (isfield(opts,'G_MT'))
    GFlag_MT=1;
    
    G_MT=opts.G_MT;
    if (max(G_MT) >n_MT*T_num|| max(G_MT) <1)
        error('\n The input G is incorrect. It should be within %d and %d',1,n);
    end
end

%% Starting point initialization
% compute AT y
if (opts.nFlag==0)
    %ATy =A'*y;
    AsTys_MT = Xsys_MT_cal(As, ys);
end

% process the regularization parameter
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
    %     if (z<0 || z>1)
    %         error('\n opts.rFlag=1, and z should be in [0,1]');
    %     end
    
    %     computedFlag=0;
    %     if (isfield(opts,'lambdaMax'))
    %         if (opts.lambdaMax~=-1)
    %             lambda=z*opts.lambdaMax;
    %             computedFlag=1;
    %         end
    %     end
    %
    %     if (~computedFlag)
    %         if (GFlag==0)
    %             lambda_max=findLambdaMax(ATy, n, ind, size(ind,2));
    %         else
    %             lambda_max=general_findLambdaMax(ATy, n, G, ind, size(ind,2));
    %         end
    %
    %         % As .rFlag=1, we set lambda as a ratio of lambda_max
    %         lambda=z*lambda_max;
    %     end
end


% The following is for computing lambdaMax
% we use opts.lambdaMax=-1 to show that we need the computation.
%
% One can use this for setting up opts.lambdaMax
% if (isfield(opts,'lambdaMax'))
%
%     if (opts.lambdaMax==-1)
%         if (GFlag==0)
%             lambda_max=findLambdaMax(ATy, n, ind, size(ind,2));
%         else
%             lambda_max=general_findLambdaMax(ATy, n, G, ind, size(ind,2));
%         end
%
%         x=lambda_max;
%         funVal=lambda_max;
%         ValueL=lambda_max;
%
%         return;
%     end
% end


% initialize a starting point
if opts.init==2
    x_MT=zeros(n_MT*T_num,1);
else
    if isfield(opts,'x0')
        x_MT=opts.x0;
        if (length(x_MT)~=n_MT*T_num)
            error('\n Check the input .x0');
        end
    else
        x_MT=AsTys_MT;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute A x
if (opts.nFlag==0)
    Ax_MT = [];
    idx_row = [1:T_num:length(x_MT)];
    for idx_t = 1:T_num
        Ax_MT = [Ax_MT;As{idx_t}*x_MT(idx_row)];
        idx_row = idx_row +1;
    end
    %Ax=A* x;
    % elseif (opts.nFlag==1)
    %     invNu=x./nu; mu_invNu=mu * invNu;
    %     Ax=A*invNu -repmat(mu_invNu, m, 1);
    % else
    %     Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
end

if (opts.init==0)
    % ------  This function is not available
    %
    % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail
    %
    
    % Here, we only support starting from zero, due to the complex tree
    % structure
    
    x_MT=zeros(n_MT*T_num,1);
end

%% The main program
% The Armijo Goldstein line search schemes + accelearted gradient descent

bFlag_MT=0; % this flag tests whether the gradient step only changes a little

if (opts.mFlag==0 && opts.lFlag==0)
    L=1;
    % We assume that the maximum eigenvalue of A'A is over 1
    
    % assign xp with x, and Axp with Ax
    xp_MT=x_MT; Axp_MT=Ax_MT; xxp_MT=zeros(n_MT*T_num,1);
    
    alphap=0; alpha=1;
    
    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s_MT=x_MT + beta* xxp_MT;
        
        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x
        
        % compute the gradient (g) at s
        As_MT=Ax_MT + beta* (Ax_MT-Axp_MT);
        
        % compute AT As
        if (opts.nFlag==0)
            
            ATAs_MT = zeros(n_MT*T_num,1);
            idx_row = [0:n_MT-1]*T_num+1;
            for idx_t = 1:T_num
                ATAs_MT(idx_row) = As{idx_t}'*As_MT((idx_t-1)*length(ys{1})+1:idx_t*length(ys{1}));
                idx_row = idx_row+1;
            end
            %ATAs=A'*As;
            %         elseif (opts.nFlag==1)
            %             ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
            %         else
            %             invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end
        
        % obtain the gradient g
        g_MT=ATAs_MT-AsTys_MT;
        
        % copy x and Ax to xp and Axp
        xp_MT=x_MT;    Axp_MT=Ax_MT;
        
        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the L1/Lq-norm regularized projection
            v_MT=s_MT-g_MT/L;
            
            % tree overlapping group Lasso projection
            ind_work_MT(1:2,:)=ind_MT(1:2,:);
            ind_work_MT(3,:)=ind_MT(3,:) * (lambda / L);
            
            if (GFlag_MT==0)
                x_MT=altra(v_MT, n_MT*T_num, ind_work_MT, size(ind_work_MT,2));
            else
                x_MT=general_altra(v_MT, n_MT*T_num, G_MT, ind_work_MT, size(ind_work_MT,2));
            end
            
            v_MT=x_MT-s_MT;  % the difference between the new approximate solution x
            % and the search point s
            
            % compute A x
            if (opts.nFlag==0)
                %Ax=A* x;
                Ax_MT = [];
                idx_row = [1:T_num:n_MT*T_num];
                for idx_t = 1:T_num
                    Ax_MT = [Ax_MT;As{idx_t}*x_MT(idx_row)];
                    idx_row = idx_row +1;
                end
%             elseif (opts.nFlag==1)
%                 invNu=x./nu; mu_invNu=mu * invNu;
%                 Ax=A*invNu -repmat(mu_invNu, m, 1);
%             else
%                 Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end
            
            Av_MT=Ax_MT -As_MT;
            r_sum_MT=v_MT'*v_MT; l_sum_MT=Av_MT'*Av_MT;
            
            if (r_sum_MT <=1e-20)
                bFlag_MT=1; % this shows that, the gradient step makes little improvement
                break;
            end
            
            % the condition is ||Av||_2^2 <= L * ||v||_2^2
            if(l_sum_MT <= r_sum_MT * L)
                break;
            else
                L=max(2*L, l_sum_MT/r_sum_MT);
                %fprintf('\n L=%5.6f',L);
            end
        end
        
        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;
        
        xxp_MT=x_MT-xp_MT;   
        y_all = [];
        for idx_t = 1:T_num
        y_all = [y_all;ys{idx_t}];
        end
        Axy_MT=Ax_MT-y_all;
        
        ValueL_MT(iterStep)=L;
        
        % compute the regularization part
        if (GFlag_MT==0)
            tree_norm_MT=treeNorm(x_MT, n_MT*T_num, ind_MT, size(ind_MT,2));
        else
            tree_norm_MT=general_treeNorm(x_MT, n_MT*T_num, G_MT, ind_MT, size(ind_MT,2));
        end
        
        % function value = loss + regularizatioin
        funVal(iterStep)=Axy_MT'* Axy_MT/2 + lambda * tree_norm_MT;
        
        if (bFlag_MT)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end
        
        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp_MT=sqrt(xxp_MT'*xxp_MT);
                if ( norm_xxp_MT <=opts.tol)
                    break;
                end
            case 4
                norm_xp_MT=sqrt(xp_MT'*xp_MT);    norm_xxp_MT=sqrt(xxp_MT'*xxp_MT);
                if ( norm_xxp_MT <=opts.tol * max(norm_xp_MT,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
        
        % restart the program every opts.rStartNum
        if (~mod(iterStep, opts.rStartNum))
            alphap=0; alpha=1;
            xp_MT=x_MT; Axp_MT=Ax_MT; xxp_MT=zeros(n_MT*T_num,1); L =L/2;
        end
    end
else
    error('\n The function does not support opts.mFlag neq 0 & opts.lFlag neq 0!');
    endfunction [ output_args ] = tree_LeastR_MT( input_args )
    %TREE_LEASTR_MT Summary of this function goes here
    %   Detailed explanation goes here
    
    
end

