function y = gevpdf(x,k,sigma,mu)
%GEVPDF Generalized extreme value probability density function, optimised for GPU.
%   Y = GEVPDF(X,K,SIGMA,MU)
%
%   See also GEVPDF, GPUARRAY.

%   See toolbox/stats/stats/gevpdf for implementation details.

%   Copyright 2021 The MathWorks, Inc.

if nargin < 1
    error(message('stats:gevpdf:TooFewInputs'));
end
if nargin < 2, k = 0;     end
if nargin < 3, sigma = 1; end
if nargin < 4, mu = 0;    end

err = internal.stats.statsizechk(4,x,k,sigma,mu);
if err > 0
    error(message('stats:gevpdf:InputSizeMismatch'));
end

% Return NaN for out of range parameters.
sigma(sigma <= 0) = NaN;

z = cast((x-mu)./sigma, "like", internal.stats.dominantType(x,k,sigma,mu));

% Process each element independently using arrayfun. This causes the
% function to be compiled into an optimised GPU kernel.
y = arrayfun( @gevpdf_scalar, z, k, sigma );
end

function y = gevpdf_scalar(z, k, sigma)
% Process one element of the GEVPDF output

% avoid Inf - Inf
if z==-Inf
    z = -z;
end
% Find the k==0 cases and fill them in.
if abs(k) < eps
    y = exp(-exp(-z) - z);
else
    % Fill the k~=0 cases
    t = z*k;
    if t <= -1
        y = zeros("like",z);
    else
        y = exp(-exp(-(1./k).*log1p(t)) - (1+1./k).*log1p(t));
    end
end
y = y ./ sigma;
end