function y = gppdf(x,k,sigma,theta)
%GPPDF Generalized Pareto probability density function, optimised for GPU.
%   Y = GPPDF(X,K,SIGMA,THETA)
%
%   See also GPPDF, GPUARRAY.

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

%   Copyright 2021 The MathWorks, Inc.

if nargin < 1
    error(message('stats:gpcdf:TooFewInputs'));
end
if nargin < 2 || isempty(k), k = 0;     end
if nargin < 3 || isempty(sigma), sigma = 1; end
if nargin < 4 || isempty(theta), theta = 0; end

err = internal.stats.statsizechk(4,x,k,sigma,theta);
if err > 0
    error(message('stats:gppdf:InputSizeMismatch'));
end
% Work out the output type
outProto = ones("like",internal.stats.dominantType(x,k,sigma,theta));

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


function y = gppdf_scalar(x, k, sigma, theta, outProto)
% Process one element of the GPPDF output

% Return NaN for out of range parameters.
if sigma <= 0
    sigma = nan("like",sigma);
end
% Normalize the input data
z = cast((x-theta)./sigma,"like",outProto);
if z<0
    z = Inf("like",z);
end

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