function retRes = Build1SOC4Cox(inY,inDelta,inX,inPreRes,inTP,inOption)
%function retRes = Build1SOC4Cox(inY,inDelta,inX,inPreRes,inTP,inOption)
% 
% Build one orthogonal component in gPOCRE4Cox. 
%
% Model: pl(eta), with eta = X*beta = \sum_j (X*varpi_j)*vartheta_j'
%                              = \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: Y response values of original data;
%   inDelta: observation status, 0 if censored, 1 if not censored;
%   inX: X;
%   inPreRes: mu,W,omega,varpi,poc(j).vartheta(varthea of j-th component).
%   inTP: tuning parameter (lambda);
%   inOption: 
%       tol -- tolerance,
%       maxit -- maximum number of iterations
%       maxvar: maximum number of selected variables;
%       ptype -- type of penalization:
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation,
%                'l1' ~ L1,
%                'scad' ~ SCAD (gamma=3.7 by default),
%                'mcp' ~ MCP (gamma=3.7 by default);
%       gamma -- parameter for SCAD or MCP (=3.7 by default);
%
% Output:
%   retRes: omega, [vartheta_1,...], X, varpi, theta;
%
% Dabao Zhang    January 6, 2015
%
% Revised by:
%
% Copyright (c) 2015 by Dabao Zhang.

error(nargchk(6, 6, nargin));
[n,p] = size(inX);    % X with n objects p variables;
k = size(inY,2);    % Y with n objects k variables;
if( n~=size(inY,1) )
    error('gPOCRE4Cox: X and Y should have the same number of rows!');
end
idxCmp = inPreRes.nCmp+1; % Index of current component

% Initial Values:
%retRes.mu = inPreRes.mu;
retRes.vartheta = [inPreRes.vartheta zeros(k,1)];
retRes.omega = zeros(p,1);
retRes.varpi = zeros(p,1);

tmpDen = ones(idxCmp,1);
for j=1:(idxCmp-1)
    tmpDen(j) = inPreRes.xomega(:,j)'...
                *bsxfun(@times,inPreRes.xomega(:,j),inPreRes.W);
end

omegaPre = retRes.omega;    omegaDiff = 1;  idIter = 0;
while( (omegaDiff>inOption.tol) && (idIter<inOption.maxit) )
    idIter = idIter+1;
    omegaPre = retRes.omega;

    % Update eta
    eta = retRes.mu+[inPreRes.xomega inX*retRes.omega]*retRes.vartheta';
    
    % Update Z
    Z = eta + (inY+inZeta/2-(1+inZeta).*ginv(eta,inLink))...
        ./((1+inZeta).*dginv(eta,inLink));

    % Robust values of Z
    Z = robustize(Z);

    % Weighted Z
    Z = bsxfun(@times,Z,inPreRes.W);

    % Update mu
    retRes.mu = sum(Z,1);   % sum(inPreRew.W) = 1;
    %disp([ 'mu = ' num2str(retRes.mu)]);
    
    % Update omega
    retRes.omega = CalcSLEV(inX,Z,inTP,inOption);
    
    %if( norm(retRes.omega)<inOption.tol )
    if( norm(retRes.omega)<eps )
        omegaDiff = 0;
        retRes.vartheta(:,idxCmp) = zeros(k,1);
    else
        % Update all vartheta
        for j=1:(idxCmp-1)
            %retRes.vartheta(:,j) = Z'*inX*inPreRes.omega(:,j)/tmpDen(j);
            retRes.vartheta(:,j) = Z'*inPreRes.xomega(:,j)/tmpDen(j);
        end
        tmpDen(idxCmp) = retRes.omega'*inX'...
                         *bsxfun(@times,inX,inPreRes.W)*retRes.omega;
        %if( tmpDen(idxCmp) < inOption.tol )
        if( tmpDen(idxCmp) < eps )
            retRes.vartheta(:,idxCmp) = zeros(k,1);
        else
            retRes.vartheta(:,idxCmp) = Z'*inX*retRes.omega/tmpDen(idxCmp);
        end
    
        omegaDiff = norm(omegaPre-retRes.omega);
    end
end

%if( norm(retRes.omega)>=inOption.tol )
if( norm(retRes.omega)>eps )
    % Update X
    retRes.theta = inX'*bsxfun(@times,inX,inPreRes.W)*retRes.omega...
                   /tmpDen(idxCmp);     % px1
    retRes.X = inX-inX*retRes.omega*retRes.theta';
    
    % Calculate varpi
    retRes.varpi = retRes.omega;
    for k=(idxCmp-1):-1:1
        retRes.varpi = retRes.varpi-inPreRes.omega(:,k)...
                       *(inPreRes.theta(:,k)'*retRes.varpi);
    end
else
    retRes.mu = inPreRes.mu;
end



%-------------------------------------------------------------------------
%
function retLEV = CalcSLEV(inX,inZ,inTP,inOption)
%function retLEV = CalcSLEV(inX,inZ,inTP,inOption)
% 
% Calculate the leading sparse eigenvecor of X'ZZ'X, that is, thresholding
% X'Z when Z is a column vector.
%
% Input:
%   inX: X;
%   inZ: Z;
%   inTP: tuning parameter (lambda);
%   inOption: 
%       tol -- tolerance,
%       maxit -- maximum number of iterations
%       ptype -- type of penalization:
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation,
%                'l1' ~ L1,
%                'scad' ~ SCAD (gamma=3.7 by default),
%                'mcp' ~ MCP (gamma=3.7 by default);
%       gamma -- parameter for SCAD or MCP (=3.7 by default);
%
% Output:
%   retLEV: the leading sparse eigenvector of X'ZZ'X.
%
% Dabao Zhang    May 28, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(4, 4, nargin));

retLEV = inX'*inZ;
switch( inOption.ptype )
    case 'ebt' % EBT w/o z-transformationg
        tmpSD = median(abs(retLEV))*inTP/norminv(0.75,0,1);
        resEB = GEBCauchy(retLEV,tmpSD,[1,1,1]);
        retLEV = resEB.mu;

    case 'ebtz' % EBT w/ z-transformationg
        retLEV = retLEV/norm(retLEV);
        retLEV = .5*log((1+retLEV)./(1-retLEV));
        tmpSD = inTP/sqrt(length(retLEV)-3);
        resEB = GEBCauchy(retLEV,tmpSD,[1,1,1]);
        retLEV = 1-2./(exp(2*resEB.mu)+1);
    
    case 'l1' % L1
        tmpSD = median(abs(retLEV))/norminv(0.75,0,1);
        tmpSLEV = max(abs(retLEV)-inTP*tmpSD*sqrt(2*log(size(inX,2))),0);
        retLEV = tmpSLEV.*sign(retLEV);
    
    case 'scad' % SCAD
        if( ~isfield(inOption,'gamma') ) inOption.gamma = 3.7; end
        tmpSD = median(abs(retLEV))/norminv(0.75,0,1);
        tmpT = inTP*tmpSD*sqrt(2*log(size(inX,2))); % Threshold of L1
        
        % |z|<=2*lambda
        tIdx = find(abs(retLEV)<=2*tmpT);
        retLEV(tIdx) = max(abs(retLEV(tIdx))-tmpT,0).*sign(retLEV(tIdx));
    
        % 2*lambda<|z|<=gamma*lambda
        tIdx = setdiff(find(abs(retLEV)<=inOption.gamma*tmpT),tIdx);
        tmpSLEV = max(abs(retLEV(tIdx))-tmpT*inOption.gamma/(inOption.gamma-1),0);
        retLEV(tIdx) = tmpSLEV.*sign(retLEV(tIdx))/(1-1/(inOption.gamma-1));

    case 'mcp' % MCP
        if( ~isfield(inOption,'gamma') ) inOption.gamma = 3.7; end
        tmpSD = median(abs(retLEV))/norminv(0.75,0,1);
        tmpT = inTP*tmpSD*sqrt(2*log(size(inX,2))); % Threshold of L1
        tIdx = find(abs(retLEV)<=inOption.gamma*tmpT);
        tmpSLEV = max(abs(retLEV(tIdx))-tmpT,0);
        retLEV(tIdx) = tmpSLEV.*sign(retLEV(tIdx))/(1-1/inOption.gamma);
    
    otherwise
        error('gPOCRE: The penalization method is not supported!');
end
if( norm(retLEV)>eps )
    retLEV = retLEV/norm(retLEV);
else
    retLEV = zeros(size(inX,2),1);
end



%-------------------------------------------------------------------------
%
function retX = robustize(inX)
%function retX = robustize(inX)
% 
% Reset the values of any outliers in X, and make subsequent statistical
% analysis robust.
%
% Input:
%   inX: X;
%
% Output:
%   retX: robust X.
%
% Dabao Zhang    November 26, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

rmean = trimmean(inX,50);
rstd = median(abs(inX-rmean))/normcdf(0.75,0,1);
cutoff = norminv(1-0.01/size(inX,1),0,1)*rstd;
retX = inX;
retX(retX<rmean-cutoff) = rmean-cutoff;
retX(retX>rmean+cutoff) = rmean+cutoff;
