%function bss = bsspdfest(x, bnds, n, m, bss)
function bss = bsspdfest(x, varargin)
%BSSPDFEST - Estimate a B-Spline series probability density function.
%
% INPUTS:
%
%    x:
%       The data to be used for the B-Spline density estimation. This
%       must be a column oriented array with each row representing a data
%       point, X(NDATA,NDIM) in size. Values of X outside the estimation
%       bounds are ignored.
%
%    varargin name/value arguments:
%
%    'bounds', bnds:
%       The lower and upper bounds of the estimation interval for each
%       dimension. This is an NDIM by 2 matrix, specifying the minimum and
%       maximum values used to create the B-spline series partition for
%       each dimension. All values must be finite, that is, inf and -inf
%       values are not permitted.
%
%       If not present, the default estimation bounds are generated from
%       the data using bsscreatebnds().
%
%       If any of the bounds are active, a boundary correction is made
%       using reflection through each active boundary. For 2- and
%       3-dimensional data reflection based corrections are also made for
%       corners and edges where multiple boundaries are active.
%
%    'partitionsize', n:
%       The nominal partition size or the number of subintervals to use for
%       the partition in each dimension. The partition size for each
%       dimension does NOT include the extra subintervals that are created
%       based on the order of the B-splines that are used. This is a scalar
%       or a row or column vector with NDIM elements. NDIM must be >=1.
%
%       If the partition size is a scalar and NDIM > 1, then this value
%       will be used as the partition size for all dimensions.
%
%       If not present, the default value obtained from the function
%       bsspartitionsize() and the number of data points is used for each
%       dimension.
%
%    'bsplineorder', m:
%       The order of the bsplines to be used. This must be an integer
%       scalar value greater than zero (M >= 1).
%
%       If not present or empty, the default value is obtained from the
%       function bssdefaultorder().
%
%    'update', bss:
%       An existing B-Spline series probability density estimate
%       that is to have more data points added to it. This allows one
%       to build a B-Spline density estimate a little bit at a time.
%
%       If this argument exists, no other name/value pair arguments are
%       allowed and their presence will result in an error.
%
% OUTPUTS:
%
%    bss_struct:
%
%       A B-Spline series data structure with the following fields.
%
%          m         - The order of the B-spline basis functions. This is a scalar
%                      value that is the same for all dimensions.
%          ndim      - The number of dimensions for the B-spline series.
%          npts      - The number of data points used to compute the coefficients.
%                      each dimension, including the extra subintervals that
%                      are created based on the order of the B-splines that
%                      are used.
%          bnds      - The estimation bounds for each dimension.
%          partition - The partition structure containing fields X, N, NPART, and H.
%                      The fields represent the partition of the estimation interval
%                      for each dimension, with X defining the boundary points for
%                      the partition subintervals, N defining the nominal partition
%                      size, NPART defining the number of boundary points, and H
%                      defining the width of the subintervals for each dimension.
%          c         - The coefficient matrix for the B-spline series
%

%
% Copyright 2014-2023 Biometrics Northwest LLC
%

%
% Modified: 2023-02-13
% Author : Kevin R. Gehringer
%
%   Converted the interface over to varargin to improve clarity and ease
%   of use.
%
%   Added logical flags for the presence of the varargin arguments to
%   guarantee they only appear once.
%
%   Allow a scalar value of the partition size for data having more than
%   one dimension. This uses the specified partition size for each
%   dimension.
%
%   Added a check for compatible dimensions when the partition size is a
%   vector.
%
% Modified: 2022-12-17
% Author : Kevin R. Gehringer
%
%    Fixed cut and past error when assigning a reasonable initial value for
%    the B-spline order when the number of input values was 4 and the value
%    of M was empty.
%
% Modified: 2022-12-04
% Author : Kevin R. Gehringer
%
%    Assigned appropriate values to missing inputs that made more sense
%    than assigning an empty value.
%
%    Removed temporary values that were used only once in if statements,
%    replacing them with function calls directly used in the if statements.
%
%    Removed dependencies on a B-spline series output from
%    bsscheckinputs().
%
% Modified: 2016-09-14
% Author : Kevin R. Gehringer
%
%    Updated documentation to indicate that boundary corrections are made
%    for dimensions <= 3.
%
% Modified: 2016-09-14
% Author : Kevin R. Gehringer
%
%    Consolidated the input checking and b-spline series creation into a
%    single function called BSSCHECKINPUTS(). Any arguments that are not
%    supplied are assigned empty values and passed into the input checker.
%
% Modified: 2016-08-11
% Author : Kevin R. Gehringer
%
%    Fixed a bug when checking the B-spline order. The bsspartitionsizeok()
%    function was called by mistake. It should have been the bssorderok()
%    function.
%
% Modified: 2016-05-01
% Author  : Kevin R. Gehringer
%
%    Added a feature to increase the partition size based on the magnitude
%    of the standard deviation. If the standard deviation is bigger than
%    one, multiply the partition size by half the standard deviation,
%    otherwise keep it the same.
%
%    Note: at some point it may be useful to consider the width of the
%    range of the data in addition to, or instead of, the standard
%    deviation.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

if ( isempty(x) )
    error('The input data was empty.');
end

if ( ~isnumeric(x) )
    error('The input data must be 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

haveBounds          = false;
havePartitionSize   = false;
haveBsplineOrder    = false;
haveUpdatePartition = false;

useDefaultBounds        = true;
useDefaultPartitionSize = true;
useDefaultBsplineOrder  = true;
updatePartition         = false;

if ( nargin > 1 )
    nv = length(varargin);
    for i = 1:2:nv
        switch lower(varargin{i})
            case 'bounds'
                if ( haveBounds )
                    error('Only one instance of the ''Bounds'' name value pair is allowed.');
                end
                bssbndsok(varargin{i+1});
                bnds             = varargin{i+1};
                useDefaultBounds = false;
                haveBounds       = true;
            case 'partitionsize'
                if ( havePartitionSize )
                    error('Only one instance of the ''PartitionSize'' name value pair is allowed.');
                end
                bsspartitionsizeok(varargin{i+1});
                if (~isscalar(varargin{i+1}))
                    if ( length(varargin{i+1}) ~= size(x,2) )
                        error('Dimension mismatch: The number of dimensions in the data and partition size did not agree.')
                    end
                end
                if ( isscalar(varargin{i+1}) )
                    if ( size(x,2) > 1 )
                        n = varargin{i+1} * ones(size(x,2),1);
                    else
                        n = varargin{i+1};
                    end
                else
                    n = varargin{i+1};
                end
                useDefaultPartitionSize = false;
                havePartitionSize       = true;
            case 'bsplineorder'
                if ( haveBsplineOrder )
                    error('Only one instance of the ''BsplineOrder'' name value pair is allowed.');
                end
                bssorderok(varargin{i+1});
                m                      = varargin{i+1};
                useDefaultBsplineOrder = false;
                haveBsplineOrder       = true;
            case 'updatepartition'
                if ( haveUpdatePartition )
                    error('Only one instance of the ''UpdatePartition'' name value pair is allowed.');
                end
                if ( nargin ~= 3 )
                    error('When updating a B-spline kernel partition only three arguments are allowed.');
                end
                if ( isempty(varargin{i+1}) )
                    error('The B-spline kernel structure to be updates was empty.');
                end
                if ( ~isstruct(varargin{i+1}) )
                    error('The B-spline kernel structure must be a structure.');
                end
                if ( ~isscalar(varargin{i+1}) )
                    error('The B-spline kernel structure must be a scalar value.');
                end
                if ( ~isfield(varargin{i+1},'m') )
                    error('The structure was not a B-spline kernel structure.');
                end
                updatePartition         = true;
                bss                     = varargin{i+1};
                useDefaultBounds        = false;
                useDefaultPartitionSize = false;
                useDefaultBsplineOrder  = false;
                haveUpdatePartition     = true;
            otherwise
                error('Invalid name ''%s''.',varargin{i});
        end
    end

end

if ( updatePartition )

    switch( size(x,2) )
        case 1
            bss = bsspdfest1d(x, [], [], [], bss);
        case 2
            bss = bsspdfest2d(x, [], [], [], bss);
        case 3
            bss = bsspdfest3d(x, [], [], [], bss);
        otherwise
            bss = bsspdfestnd(x, [], [], [], bss);
    end

else

    if ( useDefaultBounds )
        bnds = bsscreatebnds(x);
    end

    if ( useDefaultPartitionSize )
        n = bsspartitionsize(x);
    end

    if( useDefaultBsplineOrder )
        m  = bssdefaultorder();
    end

    switch( size(x,2) )
        case 1
            bss = bsspdfest1d(x, bnds, n, m, []);
        case 2
            bss = bsspdfest2d(x, bnds, n, m, []);
        case 3
            bss = bsspdfest3d(x, bnds, n, m, []);
        otherwise
            bss = bsspdfestnd(x, bnds, n, m, []);
    end

end
