function retRes = dginv(inEta,inLink)
%function retRes = dginv(inEta,inLink)
% 
% Calculate the value of \Delta g^{-1}(eta) for a generalized linear model
% with link function g(.),
%       g(E[Y}) = mu + X*beta;
%
% Input:
%   inEta: value of eta = mu + X*beta;
%   inLink: the link function g(.).
%
% Output:
%   retRes: \Delta g^{-1}(eta).
%
% Dabao Zhang    May 8, 2013;
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

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

switch inLink
    case 'logit'
        % -36<eta<36 as we want eps<1/(1+exp(eta))<1-eps
        inEta(inEta>36) = 36;
        inEta(inEta<-36) = -36;
        retRes = 1./(1+exp(inEta));
        retRes = retRes.*(1-retRes);

    case 'probit'
        retRes = normpdf(inEta,0,1);
        
    case 'comploglog'
        % -350<eta<350 as we want realmin<exp(eta*2))<realmax
        inEta(inEta>350) = 350;
        inEta(inEta<-350) = -350;
        retRes = exp(inEta-exp(inEta));

    case 'identity'
        retRes = 1;

    case 'log'
        % -350<eta<350 as we want realmin<exp(eta*2))<realmax
        inEta(inEta>350) = 350;
        inEta(inEta<-350) = -350;
        retRes = exp(inEta);
        
    case 'loglog'
        % realmin<exp(eta*2))<exp(2*eta+2*exp(eta))<realmax
        inEta(inEta>5.8) = 5.8;
        inEta(inEta<-350) = -350;
        retRes = exp(inEta+exp(inEta));

    otherwise
        if( isnumeric(inLink) )
            % eta^(1/k) with k=inLink
            % realmin < (eta^(2/k-2))/(k*k) < realmax
            texp = inLink/(2-2*inLink);
            minb = ((realmin)^texp)*(abs(inLink)^(2*texp));
            maxb = ((realmax)^texp)*(abs(inLink)^(2*texp));
            if( inLink<0 )
                maxb = max(maxb,realmin);
                minb = min(minb,realmax)
                inEta(inEta<maxb) = maxb;
                inEta(inEta>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);
                inEta(inEta<minb) = minb;
                inEta(inEta>maxb) = maxb;
            end
            % Q: should it be always positive?
            %retRes = (inEta.^(1/inLink-1))./abs(inLink); %???
            retRes = (inEta.^(1/inLink-1))./inLink;
        else
            error('  gPOCRE: The link is not supported...');
        end
end
