function [y,bss] = bsseval(x, bss, varargin)
%BSSEVAL  Evaluate the 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.
%
%       This value may be empty if a gridded interpolant is being created.
%       In this situation the output Y will be assigned an empty value.
%       This feature allows the creation of a gridded interpolant without
%       the need to evaluate it.
%
%    bss:
%       The B-Spline series data structure containing the series to be used
%       for evaluation.
%
%    varargin name/value arguments:
%
%    'function', evalfun:
%
%       The probability function that is to be evaluated.
%
%       If the B-Spline kernel estimator represents a probability density
%       function, the following output functions are available 'pdf',
%       'cdf', 'icdf', 'survivor', 'cumhazard'. If no evaluation function
%       is specified the B-Spline kernel is simply evaluated without
%       transformations.
%
%       Allowed evaluation functions.
%
%        If the B-Spline series represents a probability density function:
%
%          pdf      : The probability density function, f(x)
%          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))
%
%       The evalution functions 'icdf' and 'cumhazard' are only available
%       for 1-dimensional distributions. The 'pdf','cdf', and 'survivor'
%       evaluation functions are available for all dimensions.
%
%    'zerotolerance', 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.
%
%       The zero tolerance is only used for the 1-dimensional 'icdf'
%       evaluation function.
%
%    'createinterpolant':
%
%       Create a gridded interpolant for fast function evaluation. Allowed
%       functions are: 'pdf' and 'cdf'.
%
%       Once created the gridded interpolants are used to evaluate the
%       'pdf' and 'cdf' functions even when computing the various derived
%       functions.
%
%    'gridpointsperpartition', ngridperpart:
%
%       Specifies the number of grid points to use for each partition in
%       the B-spline series. This may be a scalar value, applying to all
%       dimensions, or a vector of length NDIM, where the numper of points
%       per partition is specified for each dimension.
%
%       The values do not need to be integers, but they need to be greater
%       than one (ngridperpart > 1). They are used to multiply the
%       partition size to obtain the total number of grid points to be used
%       in each dimension for the gridded interpolant.
%
%       If not specified, the default value for the dimensionality of the
%       data from the function bssdefaultpointsperpartition() will be used.
%
% 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.
%
%       bss:
%
%           A B-spline series data structure that has been augmented with a
%           gridded interpolant for the specified evaluation function. The
%           gridded interpolant may be used for fast approximate evaluation
%           of the function represented by the B-spline series.
%
%           Two gridded interpolant functions are allowed depending on
%           the evaluation function:
%
%               pdf - The probability density function.
%
%               cdf - The cumulative distribution function.
%
%           Using any other evaluation function with this output argument
%           present is an error.
%

%
% Copyright 2014-2023 Biometrics Northwest LLC
%

%
% Modified: 2023-02-14
% Author : Kevin R. Gehringer
%
%   Converted to use varargin for the optional arguments.
%
%   Removed 'bss' as a valid evaluation function. It is only used
%   internally.
%
%   Added logical flags for the presence of the varargin arguments to
%   guarantee they only appear once.
%
%   Removed unnecessary testing for the existence of the evaluation
%   function variable: it is guaranteed to exist.
%
% Modified: 2023-01-26
% Author : Kevin R. Gehringer
%
%   Rearranged the data and B-spline series initial tests so that the data
%   variable was checked before the B-spline series structure to agree with
%   their order in the interface.
%
% 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.
%   6) An empty value for X is allowed if NARGOUT == 2. In this case y will
%      be assigned an empty value for output.
%
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
% The number of output arguments must be 1 or 2, anything else is an error.
%
if ( nargout > 2 )
    error('Too many output arguments');
end

%
% Process the varargin arguments first since later actions and argument
% requirements are dependent on some of the varargin input values.
%

haveFunction               = false;
haveZeroTol                = false;
haveCreateInterpolant      = false;
haveGridPointsPerPartition = false;

useDefaultFunction               = true;
useDefaultZeroTol                = true;
useDefaultCreateInterpolant      = true;
useDefaultGridPointsPerPartition = true;

defaultFunction          = 'pdf';
defaultZeroTol           = 'default';
defaultCreateInterpolant = false;

if ( nargin > 2 )
    nv = length(varargin);
    i = 1;
    while ( i <= nv )
        switch lower(varargin{i})
            case 'function'
                if ( haveFunction )
                    error('Only one instance of the ''Function'' name value pair is allowed.');
                end
                if ( isempty(varargin{i+1}) )
                    error('The evaluation function name was empty.');
                end
                if ( ~(ischar(varargin{i+1}) || isstring(varargin{i+1})) )
                    error('The evaluation function name must be a character row vector or a string.');
                end
                if ( ischar(varargin{i+1}) )
                    if ( ~isrow(varargin{i+1}) )
                        error('The evaluation function name must be a character row vector.');
                    end
                else
                    if ( ~isscalar(varargin{i+1}) )
                        error('The evaluation function name must be a scalar string.');
                    end
                end
                okfuns = {'pdf'; 'cdf'; 'survivor'; 'icdf'; 'cumhazard'};
                if ( ~any(strcmpi(char(varargin{i+1}),okfuns),'all') )
                    error('Invalid evaluation function ''%s''.',varargin{i+1});
                end
                evalfun = lower(char(varargin{i+1}));
                useDefaultFunction = false;
                haveFunction = true;
                i = i + 2;
            case 'zerotolerance'% two allowed values: 'default', 'exact'
                if ( haveZeroTol )
                    error('Only one instance of the ''ZeroTolerance'' name value pair is allowed.');
                end
                if ( isempty(varargin{i+1}) )
                    error('The zero tolerance type was empty.');
                end
                if ( ~(ischar(varargin{i+1}) || isstring(varargin{i+1})) )
                    error('The zero tolerance type must be a character row vector or a string.');
                end
                if ( ischar(varargin{i+1}) )
                    if ( ~isrow(varargin{i+1}) )
                        error('The zero tolerance type must be a character row vector.');
                    end
                else
                    if ( ~isscalar(varargin{i+1}) )
                        error('The zero tolerance type must be a scalar string.');
                    end
                end
                okzt = {'default'; 'exact'};
                if ( ~any(strcmpi(varargin{i+1},okzt),'all') )
                    error('Invalid zero tolerance type ''%s''. Only ''%s'' and ''%s'' are allowed.',varargin{i+1},okzt{1},okzt{2});
                end
                zerotol = lower(char(varargin{i+1}));
                useDefaultZeroTol = false;
                haveZeroTol = true;
                i = i + 2;
            case 'createinterpolant'
                if ( haveCreateInterpolant )
                    error('Only one instance of the ''CreateInterpolant'' name value pair is allowed.');
                end
                useDefaultCreateInterpolant = false;
                createInterpolant = true;
                haveCreateInterpolant = true;
                i = i + 1;
            case 'gridpointsperpartition'
                if ( haveGridPointsPerPartition )
                    error('Only one instance of the ''GridPointsPerPartition'' name value pair is allowed.');
                end
                if ( ~isnumeric(varargin{i+1}) )
                    error('The number of grid points per partition must be numeric.')
                end
                if ( ~isscalar(varargin{i+1}) && ~isvector(varargin{i+1}) )
                    error('The number of grid points per partition must be a scalar or vector value.')
                end
                if ( ~isscalar(varargin{i+1}) )
                    if ( length(varargin{i+1}) ~= bss.ndim )
                        error('Dimension mismatch: The points per partition size did not agree with the B-spline series dimension.')
                    end
                end
                if ( ~all(varargin{i+1} > 1,'all') )
                    error('The points per partition values must all be greater than one (1).')
                end
                useDefaultGridPointsPerPartition = false;
                if ( isrow(varargin{i+1}))
                    gridPointsPerPartition = varargin{i+1}';
                else
                    gridPointsPerPartition = varargin{i+1};
                end
                haveGridPointsPerPartition = true;
                i = i + 2;
            otherwise
                error('Invalid name ''%s''.',varargin{i});
        end
    end

end

if ( haveGridPointsPerPartition && ~haveCreateInterpolant )
    error('The grid points per partition argument is only allowed when creatinga  gridded interpolant.')
end

if ( ~useDefaultFunction )
    only1dfuns = {'icdf'; 'cumhazard'};
    if ( any(strcmpi(evalfun,only1dfuns)) )
        if (bss.ndim > 1)
            error('Evaluation function ''%s'' is only valid for one dimensional distributions.',evalfun);
        end
    end
    if (strcmpi(evalfun,'icdf') )
        if ( useDefaultZeroTol )
            zerotol = defaultZeroTol;
        end
    else
        if ( ~useDefaultZeroTol )
            error('The zero tolerance is only allowed for the ''icdf'' function.');
        end
    end
else
    evalfun = defaultFunction;
end

if ( useDefaultCreateInterpolant )
    createInterpolant = defaultCreateInterpolant;
end

if ( useDefaultGridPointsPerPartition )
    gridPointsPerPartition = bssdefaultpointsperpartition(bss.ndim);
end

%
% We gotta have the data array X.
%
% Check data for validity. The data must:
%   exist
%   be nonempty unless nargout == 2
%   be numeric
%   be an NDATA by NDIM matrix
%
if ( ~exist('x','var') )
   error('The evaluation points array was not supplied.');
end

if ( isempty(x) )
    if ( ~createInterpolant )
        error('The evaluation points may only be empty when creating a gridded interpolant.');
    end
    emptyx = true;
else
    if ( ~isnumeric(x) )
        error('The evaluation points array was not numeric.');
    end
    if ( ~(iscolumn(x) || ismatrix(x)) || isscalar(x) )
        error('The data must be a column vector or a matrix with the points as rows.');
    end
    emptyx = false;
end

%
% We gotta have a B-spline series structure.
%
if ( ~exist('bss','var') )
    error('The B-spline series data structure was not provided.');
end

if ( isempty(bss) )
    error('The B-spline series data structure was empty.');
end

%
% Check B-Spline series for validity.
%
if ( bss.ndim < 1 )
    error('The number of dimensions must be at least 1 and a value of %d was found in the B-spline series structure.',bss.ndim);
end

%
% If X is not empty, get its size and make sure it is compatable with the
% B-spline series.
%
if ( ~emptyx )
    sz_x = size(x);

    if ( length(sz_x) ~= 2 )
        error('The data array must be 2-dimensional with rows representing the data points.');
    end
    if ( sz_x(2) ~= bss.ndim )
        error('Dimension mismatch: data dimension (%d) is not equal to B-spline series (%d).',sz_x(2), bss.ndim);
    end
end

%
% Validate the evaluation function and zero tolerance based on the domain
% dimension. PDF, CDF, and survivor are allowed for all dimensions. If
% zerotol is also defined for these functions it is an error.
%
okfuns = {'pdf'; 'cdf'; 'survivor'};
if ( any(strcmpi(evalfun,okfuns)) )
    if ( exist('zerotol','var') )
        error('A zero tolerance is only valid with the one dimensional function ''icdf''.');
    end
end

isicdf = false;
okfuns = {'icdf'; 'cumhazard'};
if ( any(strcmpi(evalfun,okfuns)) )
    if (bss.ndim > 1)
        error('Evaluation functions ''icdf'' and ''cumhazard'' are only valid for one dimensional distributions.');
    end
    %
    %   We know bss.ndim == 1  and the function is 'icdf' or 'cumhazard' so:
    %
    %   Check to see if the evalsuation 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') )
        isicdf = true;
        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

end

if ( createInterpolant )
    if ( nargout == 2 )
        okfuns = {'pdf';'cdf'};
        if ( ~any(strcmpi(evalfun,okfuns),'all') )
            error('Invalid evaluation function ''%s'' for creating an interpolant.',evalfun);
        end
        bss.(evalfun) = bsscreateinterpolant(bss,evalfun,gridPointsPerPartition);
    else
        error('To create a gridded interpolant requires two outputs: [y, bss_new] = bsseval(...)')
    end
end

%
%   Select the right case and compute the values, unless X is empty.
%
if ( emptyx )
    
    y = [];
    
else
    
    switch( bss.ndim )
       case 1
           if ( isicdf )
               y = bsseval1d(x,bss,evalfun,zerotol);
           else
               y = bsseval1d(x,bss,evalfun);
           end
       case 2
           y = bsseval2d(x,bss,evalfun);
       case 3
           y = bsseval3d(x,bss,evalfun);
        otherwise
           y = bssevalnd(x,bss,evalfun);
    end
end

end