function retRes = SelectModel(inRes,inMSC)
%function retRes = SelectModel(inRes,inMSC)
%
% Select the optimal model on the basis of prespecified criterion, such as
% EBIC, BIC, AIC, AICc etc.
%
% Considered EBIC:
%   Chen J and Chen Z (2010) Extended BIC for small-n-large-p sparse GLM.
%                            Statistica Sinica. doi: 10.5705/ss.2010.216.
%
%       EBIC=-2*log(L)+v(s)*log(n)+2*v(s)*gamma*log(p)
%       where gamma>max(1-0.5*log(n)/log(p),0), v(s) is #parameters in the
%       model, p is the total #variables in consideration.
%
% Input:
%   inRes: output from gPOCREPath();
%   inMSC: model selection criterion,
%           0 ~ BIC,
%           (0,1] ~ coefficient used in EBIC (Chen & Chen, 2010),
%           2 ~ AIC (default),
%           3 ~ AICc.
%
% Output:
%   retRes: lambda, nCmp, 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),
%           poc(j).mu (mu of j-th component).
%
% Dabao Zhang    May 16, 2013
%
% Revised by:
%       Dabao Zhang, September 23, 2013 -- Calculate effp.
%
% Copyright (c) 2013 by Dabao Zhang.

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

nLambda = length(inRes);
vmsc = zeros(nLambda,1);
k = zeros(nLambda,1);     % # of non-zero beta
deviance = zeros(nLambda,1);
n = length(inRes(end).W);
p = size(inRes(end).varpi,1);
for j=1:nLambda
    %nzb(j) = sum(abs(inRes(j).beta)>eps);
    k(j) = inRes(j).effp;
    deviance(j) = inRes(j).deviance;
end

switch( inMSC )
    case 2  % AIC: deviance+2*k
        vmsc = deviance+2*k;
        strMSC = 'AIC';
    case 3  % AICc: deviance+2*k*n/(n-k-1)
        vmsc = deviance+2*k*n./(n-k-1);
        strMSC = 'AICc';
    otherwise   % (E)BIC: deviance+k*log(n)+2*k*gamma*log(p)
        vmsc = deviance+k*log(n)+2*k*inMSC*log(p);
        strMSC = ['EBIC(gamma=' num2str(inMSC) ')'];
end

idxOpt = find(vmsc==min(vmsc));  idxOpt = idxOpt(1);
disp(['    Tuning Parameter (' strMSC '): ' num2str(inRes(idxOpt).lambda)]);
retRes = inRes(idxOpt);
retRes.vmsc = vmsc(idxOpt);
