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 (2008) Extended Bayesian information criteria for
%   model selection with large model spaces. Biometrika, 95: 759-771.
%
%       EBIC=-2*log(L)+v(s)*log(n)+2*gamma*log(p choose v(s))
%       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 (by default,
%       gamma = max(1-0.5*log(n)/log(p-1),0)).
%
% Input:
%   inRes: output from POCREPath(), fPOCREPath(), or gPOCREPath();
%   inMSC: model selection criterion,
%           0 ~ BIC,
%           (0,1] ~ coefficient used in EBIC,
%           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    August 5, 2016
%
% Revised by:
%       Dabao Zhang, September 22, 2016 -- Allow inMSC=[].
%       Dabao Zhang, December 15, 2017 -- Modify the tuning parameter.
%
% Copyright (c) 2016 by Dabao Zhang.

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

n = inRes(1).n;
p = inRes(1).p;
if( nargin<2 ) inMSC = []; end
if( isempty(inMSC) ) inMSC = max(1-0.499999*log(n)/log(p),0); end

switch inRes(1).id
    case {'pocre','fpocre'}
        retRes = SelectLinearModel(inRes,inMSC);
    case 'gpocre'
        retRes = SelectGLModel(inRes,inMSC);
    otherwise
        error('  SelectModel: The model is not supported...');
end



%-------------------------------------------------------------------------
%
function retRes = SelectLinearModel(inRes,inMSC)
%function retRes = SelectLinearModel(inRes,inMSC)
%
% Select the optimal linear model on the basis of prespecified criterion,
% such as EBIC, BIC, AIC, AICc etc.
%
% Input:
%   inRes: output from POCREPath();
%   inMSC: model selection criterion,
%           0 ~ BIC,
%           (0,1] ~ coefficient used in EBIC (Chen & Chen, 2008),
%           2 ~ AIC (default),
%           3 ~ AICc.
%
% Output:
%   retRes: lambda, nCmp, beta, loglik, bSparse(is beta sparse?),
%           omega, varpi, theta, vartheta
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

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

n = inRes(1).n;
p = inRes(1).p;
k = [inRes(:).effp]';
loglik = [inRes(:).loglik]';

switch( inMSC )
    case 2  % AIC: deviance+2*k
        vmsc = -2*loglik+2*sum(k,2);
        strMSC = 'AIC';
    case 3  % AICc: deviance+2*k*n/(n-k-1)
        vmsc = -2*loglik+2*sum(k*n./(n-k-1),2);
        strMSC = 'AICc';
    otherwise   % (E)BIC: deviance+k*log(n)+2*gamma*log([p choose k])
        vmsc = -2*loglik+sum(k*log(n)+2*inMSC*lognchoosek(p+1,k),2);
        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);



%-------------------------------------------------------------------------
%
function retRes = SelectGLModel(inRes,inMSC)
%function retRes = SelectGLModel(inRes,inMSC)
%
% Select the optimal model on the basis of prespecified criterion, such as
% EBIC, BIC, AIC, AICc etc.
%
% Input:
%   inRes: output from gPOCREPath();
%   inMSC: model selection criterion,
%           0 ~ BIC,
%           (0,1] ~ coefficient used in EBIC (Chen & Chen, 2012),
%           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    July 27, May 16, 2013
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

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

n = inRes(1).n;
p = inRes(1).p;
k = [inRes(:).effp]';
deviance = [inRes(:).deviance]';

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);
        vmsc = deviance+sum(k*log(n)+2*inMSC*lognchoosek(p+1,k),2);
        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);



%--------------------------------------------------------------------------
%
function retRes = lognchoosek(n,k)
%function retRes = lognchoosek(n,k)
%
% Calculate logarithm of the binomial coefficient (n, k) = n!/k!(n-k)!.
% 
% Input:
%   n: non-negative integer;
%   k: non-negative integers, may be a matrix;
%
% Output:
%   retRes: log(n, k) = log(n!/k!(n-k)!).
%
% Dabao Zhang    June 4, 2013
%
% Revised by:
%
%   Copyright 2013 Dabao Zhang.

[nR,nC] = size(k);
vk = reshape(k,nR*nC,1);
idx1 = find(vk>n/2);
vk(idx1) = n-vk(idx1);
[uk,idx1,idx2] = unique(vk);    % vk = uk(idx2);

if( any(uk<0) )   % k may be a matrix???
    error('  lognchoosek: should have 0<=k<=n...'); 
end
for j=1:length(uk)
    tmpN = (n-uk(j)+1):n;
    tmpD = 1:uk(j);
    tmpN = log(tmpN./tmpD);
    uk(j) = sum(tmpN);
end

retRes = reshape(uk(idx2),nR,nC);
