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;
%
% Note: For multivariate responses such as in multinomial models, dginv(.)
% only calculate the gradient for single observation at each time.
%
% 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:
%       Dabao Zhang, Dec. 9, 2015 -- Postive definite matrix.
%
% 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));

    case 'bclogit'  % inEta: 1xk
        ey = ginv(inEta,inLink);
        %myeps = .00000001/length(ey);
        %ey(ey<myeps) = myeps;
        %sumey = sum(ey,2);
        %if( sumey>1-myeps )
        %    ey = (ey./sumey).*(1-myeps);
        %end
        zeta = 0.1;
        ey = (ey+zeta/(length(ey)+1))/(1+zeta);
        retRes = diag(ey)-ey'*ey;

    case 'pologit'
        tmp = 1./(2+exp(-inEta)+exp(inEta));
        %retRes = diag(-tmp)+diag(tmp(2:end),1);
        retRes = diag(-tmp)+diag(tmp(2:end),-1);
        
    case 'crlogit'
        k = length(inEta);
        pp = 1./(1+exp(inEta)); % inEta: 1xk
        cpp = cumprod(pp)';  % kx1
        mat1 = tril(cpp*ones(1,k))+diag(cpp(2:end),1);
        
        mat2 = tril(-(1-pp(2:end))'*(1-pp(1:(end-1))))+diag(1-pp(2:(end-1)),1);
        mat2 = [mat2 [zeros(k-2,1); 1-pp(end)]];
        mat2 = [mat2; pp-1];
        
        %retRes = mat1.*mat2;
        retRes = (mat1.*mat2)'; % Should exchange row & column in calculations!!!!!
        
    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
