function retRes = EvalPerf4GLM(inFamily,inLink,inBeta,inRes,tsYX)
%function retRes = EvalPerf4GLM(inFamily,inLink,inBeta,inRes,tsYX)
%
% Performance evaluation of the fitting results on generalized linear
% models.
%
% Model: g(E[Y]) = mu + X*beta
%
% Input:
%   inFamily: family of the distribution;
%   inLink: link function of the generalized linear model;
%   inBeta: true value of beta;
%   inRes: lambda, nCmp, mu, beta, poc(j).mu (mu of j-th component);
%   tsYX: test data.
%
% Output:
%   mr: misclassification rate;
%   press: sum of squared prediction error;
%   ntp: number of true positives (in estimated beta);
%   ntn: number of true negatives (in estimated beta);
%   nfp: number of false positives (in estimated beta);
%   nfn: number of false negatives (in estimated beta);
%   fdr: false discovery rate (=nfp/(ntp+nfp));
%
% Dabao Zhang    May 22, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

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

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

if( nargin<5 )
    retRes.mr = [];
    retRes.press = [];
else
    %tmpEta = repmat(inRes.mu,size(tsYX.Y,1),1)+tsYX.X*inRes.beta;
    tmpEta = bsxfun(@plus,tsYX.X*inRes.beta,inRes.mu);

    %tmpEY = ginv(tmpEta,inLink);
    if( bMultiN )
        tmpEY = arrayfun(@(i)(ginv(tmpEta(i,:),inLink)),...
                         [1:size(tmpEta,1)]','UniformOutput',false);
        tmpEY = cell2mat(tmpEY);
        tmpEY = [1-sum(tmpEY,2) tmpEY];
    elseif( strcmp(inFamily,'cox') )
        % Not implemented here!
        tmpEY = 0;
    else
        tmpEY = ginv(tmpEta,inLink);
    end
    switch inFamily
        case 'binomial' 
            % miss classification
            yPred = round(tmpEY);
            retRes.mr = sum(abs(tsYX.Y-yPred))/size(tsYX.Y,1);
        case {'multibc','multipo','multicr'}
            yPred = bsxfun(@eq,tmpEY,max(tmpEY,[],2));
            retRes.mr = sum(sum(abs(tsYX.Y-yPred)))/(2*size(tsYX.Y,1));
            
        case {'normal','poisson','gamma','inverse gaussian','quasi','cox'}
            retRes.mr = [];
            
        otherwise
            error('    gPOCRE: The family is not supported...');
    end
    
    % PRESS
    retRes.press = sum(sum((tsYX.Y-tmpEY).^2))/size(tsYX.X,1);
end

% # of true positives (in estimated beta)
retRes.ntp = sum(sum((abs(inBeta)>eps).*(abs(inRes.beta)>eps)));

% # of true negative (in estimated beta)
retRes.ntn = sum(sum((abs(inBeta)<=eps).*(abs(inRes.beta)<=eps)));

% # of false positives (in estimated beta)
retRes.nfp = sum(sum((abs(inBeta)<=eps).*(abs(inRes.beta)>eps)));

% # of false negative (in estimated beta)
retRes.nfn = sum(sum((abs(inBeta)>eps).*(abs(inRes.beta)<=eps)));

% FDR
retRes.fdr = retRes.nfp/(retRes.nfp+retRes.ntp);
if( isnan(retRes.fdr) )
    retRes.fdr = 0;
end
