function [retRes,retStat] = gPOCRE(inY,inX,inFamily,inLink,inTP,inOption)
%function [retRes,retStat] = gGOCRE(inY,inX,inFamily,inLink,inTP,inOption)
% 
% Penalized Orthogonal-Component Regression to fit generalized linear
% models (gPOCRE) with high dimensional data. The deviance will be
% reported. The number of non-zero coefficients is limited (under maxvar
% =n/2 by default). Please refer gPOCRE to
%
% 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(.);
%   inTP: tuning parameter;
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%       bctype -- type of bias correction: 0 ~ no bias correction,
%                 1 ~ approximate (default), type 2 -- optimal;
%       maxcmp -- maximum number of components;
%       maxvar -- maximum number of selected variables;
%       W -- prespecified weight vector for the weight matrix;
%       ptype -- type of penalization,
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation (default),
%                '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);
%       covidx -- indices of covariates of no interest from components.
%
% Output:
%   retRes: nCmp, mu, beta, deviance, W, bSparse(is beta sparse), vartheta,
%           omega, varpi, xomega ([X_1*omega_1, ...] = X*varpi), theta,
%           poc(j).vartheta (vartheta of j-th component), poc(j).mu,
%           poc(j).beta.
%   retStat: some statistical/inference results based on the model with all
%            non-zero beta, e.g., pvalue, seqpv, indpv, loglik, effp
%            (effective number of predictors), separated.
%
% Notes on suggested Family & Link combinations,
%   'normal': identity (default), log, -1
%   'binomial': logit (default), probit, comploglog
%   'multinomial':
%   'poisson': log (default), identity, 0.5
%   'gamma': -1 (default), loglog, identity, log
%   'inverse gaussian': -2 (default), -1, identity, log
%   'quasi': logit (default), probit, comploglog, identity, log, loglog,
%            -1,-2,0.5, customized function
%
% Dabao Zhang    May 28, 2013
%
% Revised by:
%       Dabao Zhang, September 16, 2013 -- Add inOption.covidx,
%                               retRes.pvalue, retRes.seqpv, retRes.indpv;
%
%       Dabao Zhang, September 23, 2013 -- Add retStat and calculate effp.
%       Dabao Zhang, Oct. 29, 2013 -- Change the initialization of W.
%       Dabao Zhang, Oct. 31, 2013 -- Add poc(j).beta.
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(5, 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('    gPOCRE: X and Y should have the same number of rows!');
end

if( nargin<6 ) inOption = []; end
if( ~isfield(inOption,'tol') ) inOption.tol = 1e-6; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; end
if( ~isfield(inOption,'bctype') ) inOption.bctype = 1; end
if( ~isfield(inOption,'maxcmp') ) inOption.maxcmp = 50; end
if( ~isfield(inOption,'maxvar') ) inOption.maxvar = n/log(n); end
if( ~isfield(inOption,'W') ) inOption.W = []; end
if( ~isfield(inOption,'ptype') ) inOption.ptype = 'ebtz'; end
if( ~isfield(inOption,'covidx') ) inOption.covidx = []; end

if( isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        otherwise
            error('  gPOCRE: The family is not supported...');
    end
end

% Only implemented for single response!
if( k>1 )
    error('    gPOCRE: Only implemented for single response!');
end

% Initial values
retRes.nCmp = 0;
retRes.mu = glnk(mean(inY),inLink);
retRes.beta = zeros(p,k);
retRes.vartheta = zeros(k,0);
if( isempty(inOption.W) )
    retRes.W = ones(n,1)/n;
else
    retRes.W = inOption.W/sum(inOption.W);
end
retRes.varpi = zeros(p,0);
retRes.xomega = zeros(n,0);
retRes.bSparse = 1;
retRes.poc = [];
retRes.omega = zeros(p,0);
retRes.theta = zeros(p,0);
if( size(retRes.W,2)~=1 ) retRes.W = retRes.W'; end

% Centralize X appropriately according to W:
mX = mean(bsxfun(@times,inX, retRes.W),1);
tmpX = bsxfun(@minus,inX,mX);

% Construct all components
idxCmp = 0;     idxNZBeta = [];
omeganrm = 1;
fprintf('\b  ');
zeta = CalcZeta4BC(tmpX,retRes.W,inOption.bctype);
while( (idxCmp<inOption.maxcmp) && retRes.bSparse && (omeganrm>eps) )
    idxCmp=idxCmp+1;

    % Construct the component
    tmpRes = Build1SOC(inY,tmpX,inLink,retRes,inTP,zeta,inOption);
    omeganrm = norm(tmpRes.omega);
    idxNZBeta = union(idxNZBeta,find(abs(tmpRes.varpi)>eps));
    if( length(idxNZBeta)>inOption.maxvar )
        retRes.bSparse = 0;
    elseif( omeganrm>eps )        % Save the information on constructed component
        retRes.mu = tmpRes.mu;
        retRes.vartheta = tmpRes.vartheta;
        retRes.poc(idxCmp).mu = tmpRes.mu;
        retRes.poc(idxCmp).vartheta = tmpRes.vartheta;
        retRes.omega(:,idxCmp) = tmpRes.omega;
        retRes.varpi(:,idxCmp) = tmpRes.varpi;
        retRes.xomega(:,idxCmp) = tmpX*tmpRes.omega;
        retRes.theta(:,idxCmp) = tmpRes.theta;
        retRes.nCmp = idxCmp;

        % Update tmpX
        tmpX = tmpRes.X;
        fprintf('.');
    end
end
fprintf('\n');

% Calculate beta
retRes.beta = retRes.varpi*retRes.vartheta';
for j=1:retRes.nCmp
    retRes.poc(j).beta = retRes.varpi(:,1:j)*retRes.vartheta(1:j)';
end

% Update mu
retRes.mu = retRes.mu-mX*retRes.beta;
for j=1:retRes.nCmp
    retRes.poc(j).mu = retRes.poc(j).mu-mX*retRes.poc(j).beta;
end

% Calculate p-values and log-likelihood value if requested
if( nargout>1 )
    switch inFamily
        case {'normal','binomial','poisson','gamma','inverse gaussian'}
            % supported only
        case 'quasi'
            error('gPOCRE: quasi is not supported...');
        otherwise
            error('gPOCRE: The family is not supported...');
    end

    tmpX = bsxfun(@minus,inX,mX);

    % Covariates of no interest
    covx = [];
    for j = 1:retRes.nCmp
        if( sum(abs(retRes.varpi(inOption.covidx,j)))>eps )
            covx = [covx ...
              tmpX(:,inOption.covidx)*retRes.varpi(inOption.covidx,j)];
        end
    end
    ncovx = size(covx,2);

    % Calculate the p-values for each component:
    tidx = setdiff(1:p,inOption.covidx);
    tmpX = [covx tmpX(:,tidx)*retRes.varpi(tidx,:)];
    [tb,td,ts] = glmfit(tmpX,inY,inFamily,'link',inLink);
    retStat.pvalue = ts.p((ncovx+2):end);

    retStat.seqpv = retStat.pvalue;
    for j=1:(retRes.nCmp-1)
        [tb,td,ts] = glmfit(tmpX(:,1:(ncovx+j)),inY,inFamily,'link',inLink);
        retStat.seqpv(j) = ts.p(j+ncovx+1);
    end

    retStat.indpv = retStat.seqpv;
    for j=2:retRes.nCmp
        [tb,td,ts] = glmfit(tmpX(:,[1:ncovx j+ncovx]),inY,inFamily,'link',inLink);
        retStat.indpv(j) = ts.p(ncovx+2);
    end

    % Calculate the deviance for SelectModel() & Report 'PerfectSeparation'
    origstat = warning('off','all');
    warning('on','stats:glmfit:PerfectSeparation');
    msgid = '';         lastwarn('');

    tmpX = inX(:,abs(retRes.beta)>eps);
    [tmpR,tmpJ] = rref(tmpX);
    retStat.effp = length(tmpJ)+1;
    tmpX = tmpX(:,tmpJ);
    [tb,retStat.deviance] = glmfit(tmpX,inY,inFamily,'link',inLink);

    [msgstr,msgid] = lastwarn;
    retStat.separated = strcmp(msgid,'stats:glmfit:PerfectSeparation');
    warning(origstat);
end



%-------------------------------------------------------------------------
%
function retRes = glnk(inEY,inLink)
%function retRes = glnk(inEY,inLink)
% 
% Calculate g(E[Y]) for given E[Y]. Note that
%                g(E[Y]) = mu + X*beta.
%
% Input:
%   inEY: E[Y];
%   inLink: the link function g(.), such as 'identity' (default for
%           'normal'), 'log' (default for 'poisson'), 'logit' (default for
%           'binomial'), 'probit','comploglog','reciprocal','loglog'
%           (default for 'gamma'),p (numeric, default for 'inverse
%           gaussian'), user-defined link function (not supported now);
%
% Output:
%   retRes: g(inEY).
%
% Dabao Zhang    May 21, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(2, 2, nargin));

switch inLink
    case 'logit'    % g(p) = log(p/(1-p))
        % eps<=E[Y]<=1-eps
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = log(inEY./(1-inEY));

    case 'probit'   % g(p) = \Phi^{-1}(p)
        % eps<=E[Y]<=1-eps
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = norminv(inEY,0,1);

    case 'comploglog'   % g(p) = log(-log(1-p))
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = log(-log(1-inEY));

    case 'identity'     % g(m) = m
        retRes = inEY;

    case 'log'  % g(m) = log(m)
        inEY(inEY<eps) = eps;
        retRes = log(inEY);

    case 'loglog'   %g(m) = log(log(m))
        inEY(inEY<1+eps) = 1+eps;
        retRes = log(log(inEY));
    
    otherwise
        if( isnumeric(inLink) )
            % g(m) = m^k with k=inLink
            % realmin < m^k < realmax
            minb = realmin^(1/inLink);
            maxb = realmax^(1/inLink);
            if( inLink<0 )
                maxb = max(maxb,realmin);
                minb = min(minb,realmax)
                inEY(inEY<maxb) = maxb;
                inEY(inEY>minb) = minb;
            elseif( inLink==0 )
                error('    gPOCRE: The link is not supported...');
            elseif( inLink==1 )
                error('    gPOCRE: Please use identity link...');
            else
                minb = max(minb,realmin);
                maxb = min(maxb,realmax);
                inEY(inEY<minb) = minb;
                inEY(inEY>maxb) = maxb;
            end
            retRes = inEY.^inLink;
        else
            error('  gPOCRE: The link is not supported...');
        end
end



%-------------------------------------------------------------------------
%
function retZeta = CalcZeta4BC(inX,inW,inType)
%function retZeta = CalcZeta4BC(inX,inW,inType)
% 
% Calculate the diagonal elements of W^{1/2}X(X'WX)^{-1}X'W^{1/2} for bias
% correction.
%
% Input:
%   inType: 0 -- no bias correction, i.e., zeta=0_{nxn};
%           1 -- approximate values, i.e., zeta=(1-1/n)I_{nxn};
%           2 -- exact values
%
% Output:
%   retZeta: a vector.
%
% Dabao Zhang    May 23, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(3, 3, nargin));
[n,p] = size(inX);

retZeta = zeros(n,1);
switch( inType )
    %case 0
        %retZeta = retZeta + 0.25;
    case 1
        retZeta = 1-inW;    % sum(inW)=1.
    case 2
        % X=USV': nxp; X(X'X)^{-1}X'=UU'
        tmpM = bsxfun(@times,[ones(rowX,1) inX],sqrt(inW));
        [tmpU,noU1,noU2] = svd(tmpM,'econ');
        retZeta = diag(tmpU*tmpU');
end
