function [varargout] = gpcdf(x,varargin)
%GPCDF Generalized Pareto cumulative distribution function, optimised for GPU.
%   P = GPCDF(X,K,SIGMA,THETA)
%
%   See also GPCDF, GPUARRAY.

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

%   Copyright 2021 The MathWorks, Inc.

if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

if nargin < 1
    error(message('stats:gpcdf:TooFewInputs'));
end
if nargin>1 && strcmpi(varargin{end},'upper')
    %Compute upper tail and remove 'upper' flag
    uflag=true;
    varargin(end) = [];
elseif nargin>1 && ischar(varargin{end})&& ~strcmpi(varargin{end},'upper')
    error(message('stats:cdf:UpperTailProblem'));
else
    uflag=false;
end
[varargout{1:max(1,nargout)}] = localgpcdf(x,uflag,varargin{:});
end

function p = localgpcdf(x,uflag,k,sigma,theta)
% Do the actual gpcdf calculation
if nargin < 3 || isempty(k), k = 0;     end
if nargin < 4 || isempty(sigma), sigma = 1; end
if nargin < 5 || isempty(theta), theta = 0; end

err = internal.stats.statsizechk(4,x,k,sigma,theta);
if err > 0
    error(message('stats:gpcdf: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.
p = arrayfun(@gpcdf_scalar, x, k, sigma, theta, uflag, outProto);
end


function p = gpcdf_scalar(x, k, sigma, theta, uflag, outProto)
% Process one element of the GPCDF 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 = zeros("like",z);
end

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