function retRes = glnk(inEY,inLink)
%function retRes = glnk(inEY,inLink)
% 
% Calculate g(E[Y]) for given E[Y]. Note that
%                g(E[Y]) = mu + X*beta.
%
% Input:
%   inEY: E[Y];
%   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);
%
% Output:
%   retRes: g(inEY).
%
% Dabao Zhang    May 21, 2013
%
% Revised by:
%       Dabao Zhang, January 19, 2015 -- Multinomial models.
%
% Copyright (c) 2013 by Dabao Zhang.

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

switch inLink
    case 'logit'    % g(p) = log(p/(1-p))
        % eps<=E[Y]<=1-eps
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = log(inEY./(1-inEY));

    case 'bclogit'    % g(p) = {log(pj/p0)}:nxk
        inEY0 = 1-sum(inEY,2);
        % eps<=E[Yj], j=0,1,...,k
        inEY0(inEY0<eps) = eps;
        inEY(inEY<eps) = eps;
        retRes = log(bsxfun(@rdivide,inEY,inEY0));

    case 'pologit'    % g(p) = {log((p0+...+p(j-1))/(pj+...+pk)):j=1,...,k}
        inEY0 = 1-sum(inEY,2);
        tmpCS = cumsum([inEY0 inEY(:,1:(end-1))],2);
        % eps<=tmpCS<=1-eps
        tmpCS(tmpCS>1-eps) = 1-eps;
        tmpCS(tmpCS<eps) = eps;
        retRes = log(tmpCS./(1-tmpCS));

    case 'crlogit'    % g(p) = {log(p(j-1)/(pj+...+pk)):j=1,...,k}
        inEY =[1-sum(inEY,2) inEY];
        tmpCS = cumsum(inEY(:,end:-1:2),2);
        tmpCS(tmpCS<eps) = eps;
        tmpEY = inEY(:,1:(end-1));
        tmpEY(tmpEY<eps) = eps;
        retRes = log(tmpEY./tmpCS);
    
    case 'probit'   % g(p) = \Phi^{-1}(p)
        % eps<=E[Y]<=1-eps
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = norminv(inEY,0,1);

    case 'comploglog'   % g(p) = log(-log(1-p))
        inEY(inEY>1-eps) = 1-eps;
        inEY(inEY<eps) = eps;
        retRes = log(-log(1-inEY));

    case 'identity'     % g(m) = m
        retRes = inEY;

    case 'log'  % g(m) = log(m)
        inEY(inEY<eps) = eps;
        retRes = log(inEY);

    case 'loglog'   %g(m) = log(log(m))
        inEY(inEY<1+eps) = 1+eps;
        retRes = log(log(inEY));
    
    otherwise
        if( isnumeric(inLink) )
            % g(m) = m^k with k=inLink
            % realmin < m^k < realmax
            minb = realmin^(1/inLink);
            maxb = realmax^(1/inLink);
            if( inLink<0 )
                maxb = max(maxb,realmin);
                minb = min(minb,realmax);
                inEY(inEY<maxb) = maxb;
                inEY(inEY>minb) = minb;
            elseif( inLink==0 )
                error('    gPOCRE: The link is not supported...');
            elseif( inLink==1 )
                error('    gPOCRE: Please use identity link...');
            else
                minb = max(minb,realmin);
                maxb = min(maxb,realmax);
                inEY(inEY<minb) = minb;
                inEY(inEY>maxb) = maxb;
            end
            retRes = inEY.^inLink;
        else
            error('  gPOCRE: The link is not supported...');
        end
end
