function [varargout] = gevcdf(x,varargin)
%GEVCDF Generalized extreme value cumulative distribution function, optimised for GPU.
%   P = GECCDF(X,K,SIGMA,MU)
%
%   See also GEVCDF, GPUARRAY.

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

%   Copyright 2021 The MathWorks, Inc.

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

if nargin < 1
    error(message('stats:gevcdf: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)}] = localgevcdf(x,uflag,varargin{:});
end

function p = localgevcdf(x,uflag,k,sigma,mu)
% Do the actual gevcdf calculation
if nargin < 3, k = 0;     end
if nargin < 4, sigma = 1; end
if nargin < 5, mu = 0;    end

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


function p = gevcdf_scalar(z, k, uflag)
% Process one element of the GEVCDF output

% Find the k==0 cases and fill them in.
if abs(k) < eps
    if uflag == true
        p = -expm1(-exp(-z));
    else
        p = exp(-exp(-z));
    end
else
    % Find the k~=0 cases and fill them in.
    t = z*k;

    % 1 + k.*(x-mu)/sigma > 0 is the support, we'll force 0 or 1 outside that.
    if t <= -1
        if uflag == true
            p = cast(k >= 0, "like", z);
        else
            p = cast(k < 0, "like", z);
        end
    else
        if uflag == true
            p = -expm1(-exp(-(1./k).*log1p(t)));
        else
            p = exp(-exp(-(1./k).*log1p(t)));% exp(-(1 + k.*z).^(-1./k));
        end
    end
end
end