function retRes = gPOCRE1POC(inY,inX,inFamily,inLink,inPreRes,inTP,inZeta,inOption)
%function retRes = gPOCRE1POC(inY,inX,inFamily,inLink,inPreRes,inTP,inZeta,inOption)
% 
% Construct one orthogonal component in gPOCRE. Please refer details in
%
% Lin, Y., Zhang, M., and Zhang, D. (2011) Penalized orthogonal-component
% regression for fitting high-dimensional generalized linear models. Purdue
% Statistics Technical Report ???.
%
% Model: g(E[Y]) = mu + X*beta = mu + \sum_j (X*varpi_j)*vartheta_j'
%                              = mu + \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   inPreRes: mu,W,omega,varpi,poc(j).vartheta(varthea of j-th component).
%   inTP: tuning parameter (lambda);
%   inZeta: zeta for bias correction;
%   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: mu, omega, [vartheta_1,...], X, varpi, theta;
%
% Dabao Zhang    May 23, 2013
%
% Revised by:
%       Dabao Zhang, July 28, 2014 -- Use robust working values of Z.
%       Min Ren, March 17, 2015 -- Implement Multinomial proportinal odds &
%                                  continuation ratio models
%       Dabao Zhang, July 27, 2016 -- Weighting for PO & CR models
%       Dabao Zhang, Oct. 3, 2018 -- Add inFamily.
%
% Copyright (c) 2013-2018 by Dabao Zhang.

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

bLatent = 0;    % Latent Variable Model: multiple responses share identical X*beta.
if( sum(strcmp(inLink,{'pologit','crlogit'})) )
    bLatent = 1;
end

% Initial Values:
retRes.mu = inPreRes.mu;
if( bLatent )
    retRes.vartheta = [inPreRes.vartheta zeros(1,1)];
else
    retRes.vartheta = [inPreRes.vartheta zeros(k,1)];
end
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';
    eta = bsxfun(@plus,[inPreRes.xomega inX*retRes.omega]...
                 *retRes.vartheta',retRes.mu);
    
    % Update Z
    if( k==1 )  
        Z = eta + (inY+inZeta/2-(1+inZeta).*ginv(eta,inLink))...
            ./((1+inZeta).*dginv(eta,inLink));

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

    %elseif( bLatent ) % How to correct the bias???
        %Z = arrayfun(@(i)(eta(i,:)+(inY(i,:)-ginv(eta(i,:),inLink))...
        %            /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        
        %Z = cell2mat(Z);
    else
        %Z = arrayfun(@(i)(eta(i,:)+(inY(i,:)+inZeta(i)/2 ...
        %           -(.5*(k+1)*inZeta(i)+1).*ginv(eta(i,:),inLink)) ...
        %           /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        Z = arrayfun(@(i)(eta(i,:)+((inY(i,:)+inZeta(i)/2)./ ...
                   (1+.5*(k+1)*inZeta(i))-ginv(eta(i,:),inLink)) ...
                   /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        %Z = arrayfun(@(i)(eta(i,:)+((inY(i,:)+inZeta(i)/(k+1))./ ...
        %           (1+inZeta(i))-ginv(eta(i,:),inLink)) ...
        %           /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        Z = cell2mat(Z);
        
        % Robust values of Z
        Z = arrayfun(@(j)(robustize(Z(:,j))),1:k,'UniformOutput',false);
        Z = cell2mat(Z);
    end

    % 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
    if( bLatent )
        % Any better way to handle here?
        optw = 1./sum(CalcIRLSResponseVar(eta,inFamily,inLink),1);
        optw = optw/sum(optw);
        
        optZ = Z*optw';
        colZ = 1;
    else
        optZ = Z;
        colZ = k;
    end
    %retRes.omega = CalcSLEV(Z,inX,inTP,inOption);
    retRes.omega = CalcSLEV(optZ,inX,inTP,inOption);
    
    %if( norm(retRes.omega)<inOption.tol )
    if( norm(retRes.omega)<eps )
        omegaDiff = 0;
        %retRes.vartheta(:,idxCmp) = zeros(k,1);
        retRes.vartheta(:,idxCmp) = zeros(colZ,1);
    else
        % Update all vartheta
        for j=1:(idxCmp-1)
            %retRes.vartheta(:,j) = Z'*inPreRes.xomega(:,j)/tmpDen(j);
            retRes.vartheta(:,j) = optZ'*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);
            retRes.vartheta(:,idxCmp) = zeros(colZ,1);
        else
            %retRes.vartheta(:,idxCmp) = Z'*inX*retRes.omega/tmpDen(idxCmp);
            retRes.vartheta(:,idxCmp) = optZ'*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 j=(idxCmp-1):-1:1
        retRes.varpi = retRes.varpi-inPreRes.omega(:,j)...
                       *(inPreRes.theta(:,j)'*retRes.varpi);
    end
else
    retRes.mu = inPreRes.mu;
end
