%2014/09/03
%dlin5@tulane.edu
%Integrative analysis of multiple diverse omics datasets by sparse group multitask regression
%Multitask-sparse group lasso algorithm
%
%

function [x,err,Allopt]=MultitasksgLRCD_General(strcInput, lamda)
%
%%
% Function MultitasksgLRCD_General
%      Least Squares Loss with the 
%           sparse group Lasso Regularization
%
%% Problem
%
%  min  1/2 sum_|| A x - y||^2 + z_2 \|x\|_1 + z_1 * sum_j w_j
%  ||x_{G_j}||_sglasso
%
%% Input parameters:
%
%  strcInput-    Structure of input data
%                 
%  lambda -   The regularization parameter (z=[z_1,z_2] >=0)
%
%% Output parameters:
%
%  x-         Solution saved by datatype separately
%  err-       sum of residual
%  Allopt-    option including strucutre information of datasets

%% Verify and initialize the parameters
%%
if (nargin <2)
    error('\n Inputs: A, y, z, and opts (.ind) should be specified!\n');
end

N=length(strcInput);

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

%read and convert data

lamda1=lamda(1);% tuning para for group lasso regularizaiton
lamda2=lamda(2);%tuning para for lasso

%% Normalization
% mean of A for each variable(each column)
%     mu=mean(A,1); A=A-repmat(mu,n,1);
% variance of each column of A
%     nu=(sum(A.^2,1)/n).^(0.5);    A=A./repmat(nu,n,1);
%     ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.

%% Group & Others 
Allopt=struct('AllgroupN',[],'AllgroupK',[],'Allgroup',[],'Allcausalgene',[],'AllcausalgeneN',[]);
alphanew=cell(N,1);
wgt_loss=ones(N,1);
tempAllgroupN=[];
tempAllcausalgene=[];
datatype=zeros(N,1);
for i=1:N
    n(i)=length(strcInput(i).pheno);
    datatype(i)=strcInput(i).datatype;
end
uniqdatatype=unique(datatype);
% TotN=sum(n);
for i=1:N
    grpInd=strcInput(i).opt.groupInd;
    tempAllgroupN=[tempAllgroupN;grpInd];
    tempAllcausalgene=[tempAllcausalgene;strcInput(i).opt.causalgene];
    
    alphanew{i}=strcInput(i).init;
%     wgt_loss(i)=n(i)/TotN;
end
for i=1:length(uniqdatatype)
     loca_datatype=find(datatype==uniqdatatype(i));
     wgt_loss(loca_datatype)=(n(loca_datatype)/sum(n(loca_datatype)));
end

Allgroups=unique(tempAllgroupN);
Allcausalgene=unique(tempAllcausalgene);
AllcausalgeneN=length(Allcausalgene);
AllgroupN=length(Allgroups);
AllgroupK=zeros(AllgroupN,1);
for i=1:AllgroupN
    AllgroupK(i)=length(find(tempAllgroupN==Allgroups(i)));
end
Allopt.AllgroupN=AllgroupN;
Allopt.AllgroupK=AllgroupK;
Allopt.Allgroup=Allgroups;
Allopt.Allcausalgene=Allcausalgene;
Allopt.AllcausalgeneN=AllcausalgeneN;
% AllgroupK=opts.AllgroupK;
%% Starting point initialization
%Outer Loop
iterstep=0;
ErrOuter=1;

while(ErrOuter>5*1e-3)%convergent rule for the whole loop

        DInd0=ones(N,1);%the beginning index of group   
        alphaold=alphanew;        

        for i=1:AllgroupN
            rD=cell(N,1);
            t=[];
            Ind_Di=cell(N,1);
            grpIndD=zeros(N,1);
            isgrpInD=zeros(N,1);
            
            for j=1:N
                tempgrpIndD=find(strcInput(j).opt.geneName==Allgroups(i));
                if(tempgrpIndD>0)
                    isgrpInD(j)=1;
                    len_Di=strcInput(j).opt.groupK(tempgrpIndD);
                    tempInd_D=[DInd0(j):DInd0(j)+len_Di-1];
                    tempInd_D2=find(strcInput(j).opt.groupInd==tempgrpIndD);
                    Ind_Di{j}=tempInd_D;
                    
                    DInd0(j)=DInd0(j)+len_Di;
                    
                    A=strcInput(j).data;
                    y=strcInput(j).pheno;
                    n=length(y);
                    tempInit=cell2mat(alphanew(j));
                    alphaT=tempInit; Ai=A(:,tempInd_D);
                    alphaT(tempInd_D)=0;   rP=y-A*alphaT;%the residual by excluding the i-th group
                    a=Ai'*rP;
                    
                    rD{j}=rP;  
                    %for dataset A
                    for k=1:len_Di     %
                        if(abs(a(k)*wgt_loss(j))<=lamda2)
                            t1(k,1)=0;
                        else
                            t1(k,1)=sign(a(k))*(abs(a(k)*wgt_loss(j))-lamda2);
                        end
                    end
                    t=[t;t1];
                 end
            end
                
            betanew=cell(N,1);
            J=norm(t);
            if(J<=lamda1) % If the coefficients of i-th group in beta are all zero
                for j=1:N                 
                    if(isgrpInD(j))
                      tempalpha=cell2mat(alphanew(j));  
                      tempInd_Di=cell2mat(Ind_Di(j));
                      tempalpha(tempInd_Di)=0;
                      betanew{j}=tempalpha(tempInd_Di);
                      alphanew{j}=tempalpha;
                      
                    end
                end
           else
                %Inner Loop
                ErrInn=1;
                %betanew=cell(N,1);
                taonew=[];
                %get betanew and taonew
                for j=1:N
                    if(isgrpInD(j))
                      tempalpha=cell2mat(alphanew(j));  
                      tempInd_Di=cell2mat(Ind_Di(j));
                      betanew{j}=tempalpha(tempInd_Di);
                      taonew=[taonew;tempalpha(tempInd_Di)];
                    end
                end

                while(ErrInn>1e-3)
                    taoold=taonew;
                    tao2=norm(taonew)+eps;
                    Lamtheta=lamda1/tao2;
                    taonew=[];
                    for j=1:N
                            if(isgrpInD(j))
                                 thetanew=cell2mat(betanew(j));
                                 
                                 A=strcInput(j).data;
                                 Ai=A(:,cell2mat(Ind_Di(j)));
                                 rP=cell2mat(rD(j));
                                 n=size(A,1);
                                 
                                for Inn_i=1:length(thetanew)
                                       % Compute the residual by excluding the Inn_i -th variable 
                                        thetaT=thetanew;    thetaT(Inn_i)=0;
                                        ri=rP-Ai*thetaT;

                                        %Soft-thresholding to solve: S[Ak(j)t*ri,lamda2)]+/(Ak(j)t*Ak(j)+lamda1/||theta||)
                                        AkTri=((Ai(:,Inn_i))'*ri)*wgt_loss(j);        AkTri2=((Ai(:,Inn_i))'*(Ai(:,Inn_i)))*wgt_loss(j);                 

                                        if(abs(AkTri)>lamda2)
                                            SoftD=sign(AkTri)*(abs(AkTri)-lamda2);
                                            thetanew(Inn_i)=SoftD/(AkTri2+Lamtheta);
                                        else
                                            thetanew(Inn_i)=0;
                                        end

                                end
                                %update betanew and taonew
                                betanew{j}=thetanew;
                                taonew=[taonew;thetanew];
                            end
                    end
                    
                    ErrInn=norm(taoold-taonew);
                end
                
                %update alphanew
                 for j=1:N
                    if(isgrpInD(j))
                        tempInd_Di=cell2mat(Ind_Di(j));
                        tempalpha=cell2mat(alphanew(j));
                        tempbeta=cell2mat(betanew(j));
                        tempalpha(tempInd_Di)=tempbeta;
                        alphanew{j}=tempalpha;
                        
                    end
                end
               
            end
                %updata All_x and All_xInd
            
        end
        %cal total error
        ErrOuter=0;
        for j=1:N        
            ErrOuter=ErrOuter+norm(cell2mat(alphaold(j))-cell2mat(alphanew(j)));
        end
        ErrOuter=ErrOuter/N;
        %     ErrOuter=norm(alphanew)+norm(betanew)+norm(damanew)
        iterstep=iterstep+1;
end
 x=alphanew;
 for j=1:N        
    tmperr(j)=-1*sum((strcInput(j).pheno-strcInput(j).data*cell2mat(x(j))).^2)/2;
 end
 err=sum(tmperr);