function y = bsseval1d(x, bss, evalfun, zerotol)
%BSSEVAL1D  Evaluate the 1-dimensional B-Spline series at x.
%
% INPUTS:
%
%    x:
%       The data to be used for evaluating the B-Spline series. This is an
%       NPOINTS by NDIM matrix with the rows representing the points.
%
%       X must be a column vector with NDIM = 1.
%
%    bss:
%       The B-Spline series data structure containing the series to be used
%       for evaluation.
%
%    evalfun:
%       The function that is to be evaluated. If the B-Spline series
%       represents a probability density function, the following output
%       functions are available 'pdf' (the default), 'cdf', 'icdf',
%       'survivor', 'cumhazard'. If no evaluation function is specified the
%       B-Spline series is simply evaluated without transformations.
%
%       Allowed evaluation functions for the B-Spline series probability
%       density function estimate:
%
%          pdf      : The probability density function, f(x) (default)
%          cdf      : The cumulative distribution function, F(x)
%          icdf     : The inverse cumulative distribution function,
%                     F^(-1)(p)
%          survivor : The survivor function, S(x) = 1-F(x)
%          cumhazard: The cumulative hazard function,
%                     -log(1-F(x)) = -log(S(x))
%
%    zerotol:
%
%       A zero tolerance for finding flat spots in the CDF when computing
%       the ICDF. It is also used to determine a minimal support interval
%       for the CDF for interpolation of the ICDF. If a value is not
%       supplied a default value of sqrt(eps) is used. A value of zero is
%       not recommended. A nearly exact zero tolerance may be specified by
%       using the string value 'exact', indicating that a zero tolerance of
%       2048*eps will be used. The string value 'default' will use the
%       default value of sqrt(eps). Large values of the zero tolerance may
%       give inaccurate or incorrect results.
%
% OUTPUTS:
%
%    y:
%       The values for the B-Spline series evaluated at x.
%
%       If the B-Spline series represents a probability density function
%       and one of the allowed functions is specified, then the output is
%       an approximation of the desired function based on the probability
%       density function represented by the B-Spline series.
%

%
% Copyright 2014-2018 Biometrics Northwest LLC
%

%
% Modified: 2023-03-07
% Author : Kevin R. Gehringer
%
%    Commented out the input checks. These have already been performed by
%    BSSEVAL and are not necessary here.
%
%    Removed 'bss' as a valid function name.
%
% Modified: 2018-08-29
% Author : Kevin R. Gehringer
%
%    Removed unnecessary comments and commented out test code.
%
%    Removed unnecessary intermediate variables that were used only once.
%
% Modified: 2016-10-10
% Author : Kevin R. Gehringer
%
%    If the number of inbounds points was zero an error occurred, but the
%    output values should all have had a value of zero. This was fixed. If
%    all points are out of bounds zero values are returned.
%
% Modified: 2016-09-14
% Author : Kevin R. Gehringer
%
%    Fixed a problem in do_bss that caused a size incompatibility when
%    there was only a single point within the bounds. In this instance the
%    evaluation of 
%
%       bss.partition.x(idxp)
%
%    returned a column vector (since bss.partition.x is a column vector)
%    not a row vector (the orientation of idxp).
%
% Modified: 2016-03-12
% Author : Kevin R. Gehringer
%
%   Several enhancements were made to improve the performance:
%
%   1) Replaced outer products with vectors of ones used to replicate data
%      with calls to REPMAT. Since REPMAT is now a built-in function it is
%      significantly faster than using the outer product to do replication.
%   2) Removed unnecessary preallocation of the output values.
%   3) Removed unnecessary reshaping for B-splines of order 1.
%   4) Changed the matrix orientation from m by n to n by m, where m is the
%      B-spline order and n is the total number of evaluation points
%      specified by the input X.
%   5) Removed unnecessary temporary variables.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
%   Validate the evaluation function and zero tolerance based on the domain
%   dimension. Valid values are: 'bss', 'pdf', 'cdf', 'icdf', 'survivor',
%   and 'cumhazard'.
%
    okfuns = {'pdf';'cdf';'icdf';'survivor';'cumhazard'};
    if ( ~any(strcmpi(evalfun,okfuns)) )
        error('The evaluation function ''%s'' was not recognized.',evalfun);
    end
%
%   Check to see if the evauation function is 'icdf' and a zero
%   tolerance was specified. If  the function was 'icdf' and no zero
%   tolerance was specified, assign 'default'. If the function was not
%   'icdf' and a zero tolerance was specified, signal an error.
%
    if (strcmpi(evalfun,'icdf') )
        if ( ~exist('zerotol','var') )
            zerotol = 'default';
        end
    else
        if ( exist('zerotol','var') )
            error('The zero tolerance is only allowed for the ''icdf'' function.');
        end
    end

%
%   Call the appropriate function.
%
    switch (lower(evalfun))
        case {'pdf'}
            y = do_pdf(x,bss);
        case {'cdf'}
            y = do_cdf(x,bss);
        case {'icdf'}
            y = do_icdf(x,bss,zerotol);
        case {'survivor'}
            y = do_survivor(x,bss);
        case {'cumhazard'}
            y = do_cumhazard(x,bss);
        otherwise
            error('The function ''%s'' was not recognized.',evalfun);
    end
    
end

function y = do_bss(x,bss)

mp1to0 = (-bss.m+1:0);

y = zeros(length(x),1);

inbounds = find(x >= bss.bnds(1) & x < bss.bnds(2));
nin  = length(inbounds);

if ( nin == 0 )
    return;
end

xtmp = (x(inbounds) - bss.bnds(1))./bss.partition.h;

idxp = repmat(bss.m + floor(xtmp),1,bss.m)+repmat(mp1to0,nin,1);
if (nin > 1)
    xxp  = repmat(xtmp,1,bss.m) - (bss.partition.x(idxp) - bss.bnds(1))/bss.partition.h;
else
    xxp  = repmat(xtmp,1,bss.m) - (bss.partition.x(idxp)' - bss.bnds(1))/bss.partition.h;
end
tmpb = bssfbspline(xxp(:),bss.m).* bss.c(idxp(:));

if ( bss.m == 1)
    y(inbounds) = tmpb;
else
    y(inbounds) = sum(reshape(tmpb,nin,bss.m),2);
end

end

function y = do_pdf(x,bss)

    do_interp = isfield(bss,'pdf') && ~isempty(bss.pdf);
    
    if ( ~do_interp )
        y = do_bss(x,bss);
    else
        y = bss.pdf(x);
        y(y<0) = 0;
    end

end


function [y,xcdf,ycdf] = do_cdf(x,bss)

% x can be empty if nargout == 3 to skip generating y values when
% evaluating the icdf function. Signalling an error will be misleading, so
% we just return an empty y.
emptyx = isempty(x);
if ( emptyx )
    y = [];
else
    y = zeros(length(x),1);
    y(x >= bss.bnds(2)) = 1;
    inbounds = ((x >= bss.bnds(1)) & (x < bss.bnds(2)));
end

do_interp = isfield(bss,'cdf') && ~isempty(bss.cdf);

if (~do_interp)

    np =  floor(bssdefaultpointsperpartition(bss.ndim)*bss.partition.n)+1;
    xcdf = linspace(bss.bnds(1),bss.bnds(2),np)';
    %ycdf = cumtrapz(xcdf,do_bss(xcdf,bss));
    ycdf = (xcdf(2)-xcdf(1))*cumtrapz(do_bss(xcdf,bss));

    % fix any roundoff errors that may have made the CDF greater than one.
    if ( ~emptyx )
        ycdf(ycdf>1) = 1;
        f = griddedInterpolant({xcdf},ycdf,'linear');
        y(inbounds) = f(x(inbounds));
    end

else
    if ( ~emptyx )
        y(inbounds) = bss.cdf(x(inbounds));
    end
    if (nargout == 3)
        xcdf = bss.cdf.GridVectors{1};
        ycdf = bss.cdf.Values;
    end
end

if ( ~emptyx )
    y(y<0) = 0;
    y(y>1) = 1;
end

end

function y = do_survivor(x,bss)

    y = 1 - do_cdf(x,bss);
    
end

function y = do_cumhazard(x,bss)

    y = zeros(length(x),1);
    s = do_survivor(x,bss);
    sgt0 = s>0;
    y(~sgt0) = nan;
    y(sgt0) = -log(s(sgt0));

end

function x = do_icdf(p,bss,zerotol)
%
%   Define a default zero tolerance and an exact zero tolerance.
%
    defaulttol = default_zero_tol();
    exacttol   = exact_zero_tol();

    if (exist('zerotol','var'))
        if ( ischar(zerotol) )
            switch lower(zerotol)
                case {'exact'}
                    doexact = true;
                    zerotol = exacttol;
                case {'default'}
                    doexact = false;
                    zerotol = defaulttol;
                otherwise
                    error('Unrecognized zero tolerance string ''%s''.',zerotol);
            end
        elseif (isnumeric(zerotol))
            if (zerotol == 0 )
                doexact = true;
                zerotol = exacttol;
            elseif(zerotol < exacttol)
                warning('The zero tolerance %e was less than the exact tolerance, using the exact tolerance %e.',zerotol,exacttol);
                doexact = true;
                zerotol = exacttol;
            else
                doexact = false;
            end
        else
            error('The zero tolerance was not a recognized string or a numeric value.')
        end
    else
        doexact = true;
        zerotol = exacttol;
    end
    
    if ( ~doexact && (zerotol < 0 ) )
        warning('The zero tolerance was negative, using the default value of %e.',defaulttol);
        doexact = false;
        zerotol = defaulttol;
    end
    
    if (~doexact && (zerotol > defaulttol) )
        warning('The zero tolerance was larger than %e, and may give incorrect results for extreme values in [0,1].',defaulttol);
    end
%
%   p must have values in the interval [0,1]
%
    badp = (p<0 | p>1);
%
%   Compute the cdf
%
    [~,xcdf,ycdf] = do_cdf([],bss);
    ncdf = length(xcdf);
%
%   find the beginning and ending index values for 0 and 1 values of the CDF. 
%
    idx0 = find(ycdf<zerotol,1,'last');
    idx1 = find(abs(1-ycdf)<zerotol,1,'first');

    if ( isempty(idx0) || isempty(idx1) )
        if ( doexact )
            warning('A minimal support interval could not be determined using the exact tolerance %e, trying %e.',exacttol,defaulttol);
            zerotol = defaulttol;
        elseif (zerotol < defaulttol) 
            warning('A minimal support interval could not be determined using a tolerance of %e, trying %e.',zerotol,defaulttol);
            zerotol = defaulttol;
        else
            cdfztol = 2*max([min(ycdf) min(abs(1-ycdf))]);
            warning('A minimal support interval could not be determined using a tolerance of %e, trying %e.',zerotol,cdfztol);
            zerotol = cdfztol;
        end
        idx0 = find(ycdf<zerotol,1,'last');
        idx1 = find(abs(1-ycdf)<zerotol,1,'first');
    end
    
    if ( isempty(idx0) || isempty(idx1) )
        error('A minimal support interval could not be determined with a zero tolerance of %e.',zerotol);
    end
%
%   Fixup the leading and trailing points of the ICDF support interval to
%   remove similar small values that break the interpolation.
%
    if (all(diff(ycdf)<zerotol))
        error('The zero tolerance %e was too big. A minimal support interval could not be determined.',zerotol);
    end
    
    while ( (idx0 < ncdf) && (ycdf(idx0+1) < zerotol+exacttol) )
        idx0 = idx0+1;
    end
    
    while ( (idx1>1) && (abs(1-ycdf(idx1)) < zerotol+exacttol) )
        idx1 = idx1-1;
    end
    ycdf = ycdf(idx0:idx1);
    xcdf = xcdf(idx0:idx1);

%
%   Find the flat spots and remove them
%
    idxflat = find(abs(diff([1;ycdf]))<zerotol);
    
    if ( ~isempty(idxflat) )
        idxnot1 = find([0;diff(idxflat)]~=1);
        nnot1 = length(idxnot1);
        for inot1 = 1:nnot1
            if ( inot1 < nnot1)
                xcdf(idxflat(idxnot1(inot1):idxnot1(inot1+1))-1) = nan;
                ycdf(idxflat(idxnot1(inot1):idxnot1(inot1+1))-1) = nan;
            else
                xcdf(idxflat(idxnot1(inot1):end)) = nan;
                ycdf(idxflat(idxnot1(inot1):end)) = nan;
            end
        end
        xcdf = xcdf(~isnan(xcdf));
        ycdf = ycdf(~isnan(ycdf));
    end
%
%   Finally interpolate. We allow extrapolation here because the range for
%   p has already been checked and verified to have values in the interval
%   [0,1]. We just need to be sure that we can get values for the endpoints
%   of 0 and 1.
%
    x = nan(size(p));
    f = griddedInterpolant({ycdf},xcdf,'linear');
    x(~badp) = f(p(~badp));
   
end

function ztol = default_zero_tol()
    ztol = sqrt(eps);
end

function ztol = exact_zero_tol()
    ztol = 2048*eps;
end
