function retRes = gPOCREPath(inY,inX,inFamily,inLink,inOption)
%function retRes = gPOCREPath(inY,inX,inFamily,inLink,inOption)
%
% Apply the Penalized Orthogonal-Component Regression to generalized linear
% models (gPOCRE) with all possible values of the tuning parameter. The
% deviance of each model will be reported. Please refer gPOCRE to
%
% Lin, Y., Zhang, M., and Zhang, D. (2011) Penalized orthogonal-component
% regression for high-dimensional generalized linear models. Purdue
% Statistics Technical Report ???.
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inFamily: the distribution family, such as 'normal', 'binomial',
%             'poisson', 'gamma', 'inverse gaussian', 'quasi' (not
%             supported now), 'multibc', 'multipo','multicr';
%   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);
%   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;
%       maxtps -- maximum number of tuning parameters;
%       W -- 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);
%       delta -- pace for choosing value of the tuning parameter;
%       lambda0 -- initial value of the tuning parameter (1 by default);
%       covidx -- indices of covariates of no interest from components.
%
% Output:
%   retRes: lambda, nCmp, mu, beta, deviance, W, bSparse(is beta sparse?),
%           omega, varpi, xomega ([X_1*omega_1, ...] = X*varpi), theta,
%           poc(j).vartheta (varthea of j-th component), pvalue, seqpv,
%           indpv, poc(j).mu (mu of j-th component), separated.
%
% Dabao Zhang    June 4, 2013
%
% Revised by:
%       Dabao Zhang, September 16, 2013 -- Add inOption.covidx,
%                               retRes.pvalue, retRes.seqpv, retRes.indpv;
%
%       Dabao Zhang, September 23, 2013 -- Calculate effp.
%       Dabao Zhang, September 26, 2013 -- Add inOption.maxtps.
%       Dabao Zhang, Oct. 29, 2013 -- Change the initialization of W.
%       Dabao Zhang, August 2, 2014 -- Calculate deviance-based R-squares.
%       Dabao Zhang, July 27, 2016 -- Output n, p, & id='gpocre'.
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(3, 5, nargin));

if( nargin<5 ) inOption = []; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; end
if( ~isfield(inOption,'maxtps') ) inOption.maxtps = 500; end
if( ~isfield(inOption,'delta') ) inOption.delta = 0.1; end
if( ~isfield(inOption,'lambda0') ) inOption.lambda0 = 1; end
if( ~isfield(inOption,'W') ) inOption.W = []; end
if( (nargin<4)||isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'multibc'
            inLink = 'bclogit';
        case 'multipo'
            inLink = 'pologit';
        case 'multicr'
            inLink = 'crlogit';
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        case 'cox'
            inLink = 'cox';
        otherwise
            error('  gPOCREPath: The family is not supported...');
    end
end

[n,p] = size(inX);    % X with n objects p variables;
if( n~=size(inY,1) )
    error('  gPOCREPath: inX and inY should have the same number of rows!');
end

% Special for multinomial models
bMultiN = 0;
if( sum(strcmp(inFamily,{'multibc','multipo','multicr'})) )
    bMultiN = 1;
end

% Calculate Null Deviance
if( bMultiN )
    [tb,dev0] = mnrfit(zeros(n,0),inY);
elseif( strcmp(inFamily,'cox') )
    [tb,dev0] = coxphfit(zeros(n,0),inY(:,1),'censoring',inY(:,2));
    
    % Prepares for fast calculation: inY to [ObsTime CStatus RIdx CIdx]
    % Sort the data
    [inY(:,1),tmpIdx] = sort(inY(:,1));
    inY(:,2) = inY(tmpIdx,2);
    inX = inX(tmpIdx,:);
    
    % Calculate the rank of failture times for fast calculation
    [~,inY(:,3)] = ismember(inY(:,1),flipud(inY(:,1)),'legacy');
    inY(:,3) = length(inY(:,1))+1-inY(:,3);
    [~,inY(:,4)] = ismember(inY(:,1),flipud(inY(:,1)));
    inY(:,4) = length(inY(:,1))+1-inY(:,4);
else
    [tb,dev0] = glmfit(zeros(n,0),inY,inFamily,'link',inLink);
end

% Explore all tuning parameter values such that 0<#[beta~=0]<maxvar
idxM = 0;           lambda = inOption.lambda0;
if( isempty(inOption.W) )
    bInitW = 1;
    inOption.W = ones(n,1)/n;
else
    bInitW = 0;
end
while( idxM<=inOption.maxtps )
    disp(['    Tuning Parameter: ' num2str(lambda(end))]);

    idxM = idxM + 1;
    [tmpRes,tmpS] = gPOCRE(inY,inX,inFamily,inLink,lambda,inOption);
    if( bInitW&&(tmpRes.nCmp>0) ) % Update W if necessary
        fprintf('\b  / ');
        if( strcmp(inFamily,'cox') )
            inOption.W = CalcOptCoxW(inX*tmpRes.beta,inY(:,2),inY(:,3:4));
        else
            %inOption.W = CalcOptW(tmpRes.mu+inX*tmpRes.beta,inFamily,inLink);
            inOption.W = CalcOptW(bsxfun(@plus,inX*tmpRes.beta,tmpRes.mu),inFamily,inLink);
        end
        [tmpRes,tmpS] = gPOCRE(inY,inX,inFamily,inLink,lambda,inOption);
    end
    
    % Save the results
    retRes(idxM).nCmp = tmpRes.nCmp;
    retRes(idxM).mu = tmpRes.mu;
    retRes(idxM).beta = tmpRes.beta;
    retRes(idxM).W = tmpRes.W;
    retRes(idxM).bSparse = tmpRes.bSparse;
    retRes(idxM).n = tmpRes.n;
    retRes(idxM).p = tmpRes.p;
    retRes(idxM).id = tmpRes.id;
    retRes(idxM).varpi = tmpRes.varpi;
    retRes(idxM).poc = tmpRes.poc;
    retRes(idxM).lambda = lambda;
    retRes(idxM).nzBeta = sum(abs(retRes(idxM).beta)>eps);
    retRes(idxM).deviance = tmpS.deviance;
    retRes(idxM).rsq = 1-tmpS.deviance/dev0;
    retRes(idxM).separated = tmpS.separated;
    retRes(idxM).effp = tmpS.effp;
    retRes(idxM).pvalue = tmpS.pvalue;
    retRes(idxM).seqpv = tmpS.seqpv;
    retRes(idxM).indpv = tmpS.indpv;
        
    if( lambda>=inOption.lambda0 )
        if( sum(retRes(idxM).nzBeta==0)&&(retRes(idxM).bSparse) )
            lambda = inOption.lambda0-inOption.delta;
            inOption.W = retRes(1).W;
        else
            lambda = lambda+inOption.delta;
        end
    elseif( (lambda<inOption.lambda0)&&(retRes(idxM).bSparse) )
        lambda = lambda-inOption.delta;
    else
        break;
    end
end

if( idxM>inOption.maxtps )
    warning('gPOCREPath: over the limit to # of tuning parameters!');
end
