function retRes = POCREPath(inY,inX,inOption)
%function retRes = POCREPath(inY,inX,inOption)
%
% Fitting the Penalized Orthogonal-Component Regression (POCRE) with all
% possible values of the tuning parameter. The loglikelihood function value
% of each model will be reported. 
%
% Input:
%   inY: Y, response values of original data;
%   inX: X;
%   inOption:
%       tol -- tolerance,
%       maxit -- maximum number of iterations,
%       maxcmp -- maximum number of components;
%       maxvar -- maximum number of selected variables;
%       maxtps -- maximum number of tuning parameters;
%       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, varpi, vartheta, omega, theta, ...
%           loglik, bSparse(is beta sparse?), sigmae2.
%
% 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, November 6, 2013 -- p-values for multiple Y (k>1) need
%                                        to be implemented using mvregress.
%       Dabao Zhang, June 16, 2014 -- Output sigmae2, rsq (R^2).
%       Dabao Zhang, October 30, 2014 -- Output xShift.
%       Dabao Zhang, July 26, 2016 -- Output n,p,id to retRes.
%       Dabao Zhang, July 10, 2017 -- Allow preset lambda sets.
%
% Copyright (c) 2013-2017 by Dabao Zhang.

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

if( nargin<3 ) inOption = []; end
if( ~isfield(inOption,'delta') ) inOption.delta = 0.1; end
if( ~isfield(inOption,'lambda0') ) inOption.lambda0 = 1; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; end
if( ~isfield(inOption,'maxtps') ) inOption.maxtps = 500; end

[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('  POCREPath: inX and inY should have the same number of rows!');
end

% Explore all tuning parameter values such that 0<#[beta~=0]<maxvar
idxM = 0;
if( length(inOption.lambda0)>1 )
    inOption.maxtps = length(inOption.lambda0);
else
    lambda = inOption.lambda0;
end
while( idxM<inOption.maxtps )
    idxM = idxM + 1;

    if( length(inOption.lambda0)>1 )
        disp(['    Tuning Parameter: ' num2str(inOption.lambda0(idxM))]);
        [tmpRes,tmpStat] = POCRE(inY,inX,inOption.lambda0(idxM),inOption);
    else
        disp(['    Tuning Parameter: ' num2str(lambda)]);
        [tmpRes,tmpStat] = POCRE(inY,inX,lambda,inOption);
    end
    
    % Save the results
    retRes(idxM).n = tmpRes.n;
    retRes(idxM).p = tmpRes.p;
    retRes(idxM).id = tmpRes.id;
    retRes(idxM).nCmp = tmpRes.nCmp;
    retRes(idxM).mu = tmpRes.mu;
    retRes(idxM).beta = tmpRes.beta;
    retRes(idxM).sigmae2 = tmpRes.sigmae2;
    retRes(idxM).rsq = 1-tmpRes.sigmae2./var(inY);
    retRes(idxM).bSparse = tmpRes.bSparse;
    retRes(idxM).varpi = tmpRes.varpi;
    retRes(idxM).vartheta = tmpRes.vartheta;
    retRes(idxM).xShift = tmpRes.xShift;
    if( length(inOption.lambda0)>1 )
        retRes(idxM).lambda = inOption.lambda0(idxM);
    else
        retRes(idxM).lambda = lambda;
    end
    retRes(idxM).xShift = tmpRes.xShift;
    %retRes(idxM).sigmae2 = tmpStat.sigmae2;
    if( k==1 )
        retRes(idxM).pvalue = tmpStat.pvalue;
        retRes(idxM).seqpv = tmpStat.seqpv;
        retRes(idxM).indpv = tmpStat.indpv;
    end
    retRes(idxM).effp = tmpStat.effp;
    retRes(idxM).loglik = tmpStat.loglik;
    retRes(idxM).nzBeta = sum(sum(abs(retRes(idxM).beta)>eps));

    if( length(inOption.lambda0)<=1 )
        if( lambda>=inOption.lambda0 )
            if( (retRes(idxM).nzBeta==0)&&(retRes(idxM).bSparse) )
                if( retRes(1).bSparse )
                    lambda = inOption.lambda0-inOption.delta;
                else
                    break;
                end
            else
                lambda = lambda+inOption.delta;
            end
        elseif( (lambda<inOption.lambda0)&&(retRes(idxM).bSparse) )
            lambda = lambda-inOption.delta;
        else
            break;
        end
    end
end

if( length(inOption.lambda0)<=1 )
    % Sort retRes according to lambda values:
    lambda = [retRes(1:end).lambda];
    [tmpL, idxL] = sort(lambda);
    retRes = retRes(idxL);
end
