function bss_out = bsspdfestnd(x, bnds, n, m, bss)
%BSSPDFESTND - Estimate a N-D 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.
%
%       The number of columns NDIM must be >= 1.
%
%    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.
%
%       NDIM must be >= 1.
%
%       If empty or not present, the 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. No corrections are
%       made for corners and edges where multiple boundaries may be active.
%
%    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 row or
%       column vector with NDIM elements. NDIM must be >= 1.
%
%       If empty or not present, the default value obtained from the
%       function bsspartitionsize() and the number of data points is used
%       for each dimension.
%
%    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().
%
%    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, then BNDS, N, and M are ignored and the
%       values in BSS are used.
%
% 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.
%                      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-2022 Biometrics Northwest LLC
%

%
% 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.
%
%    Fixed a bug where an input B-spline series to be updateed would be
%    overwritten by a new B-spline series.
%
%    Commented out redundant initial value assignmments and input testing
%    as this is done by bsspdfest(). This code will likely be deleted at
%    some point as this function is not intended to be called directly.
%
% 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().
%
%    Added a call to bsscreatestruct() after successful input validation.
%
% Modified: 2018-08-22
% Author  : Kevin R. Gehringer
%
%     Added reflection of the data adjacent to a boundary to its other side
%     to deal with active boundaries. Only reflection through the boundary
%     is used, no reflection through corners when intersecting boundaries
%     are active.
%
%     Replaced calls to remat with calls to bsxfun to improve performance.
%
% Modified: 2018-08-17
% Author : Kevin R. Gehringer
%
%    Fixed an off by one error when generating linear index values.
%    bss.partition.npart was used and it should have been
%    bss.partition.npart-1.
%
% 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 alues 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-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 use of temporary variables.
%   3) Removed an unnecessary transpose of the partition sizes.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

% switch (nargin)
%     case 1
%         bnds = bsscreatebnds(x);
%         n    = bsspartitionsize(x);
%         m    = bssdefaultorder();
%         bss  = [];
%     case 2
%         if ( isempty(bnds) )
%             bnds = bsscreatebnds(x);
%         end
%         n    = bsspartitionsize(x);
%         m    = bssdefaultorder();
%         bss  = [];
%     case 3
%         if ( isempty(bnds) )
%             bnds = bsscreatebnds(x);
%         end
%         if ( isempty(n) )
%             n    = bsspartitionsize(x);
%         end
%         m    = bssdefaultorder();
%         bss  = [];
%     case 4
%         if ( isempty(bnds) )
%             bnds = bsscreatebnds(x);
%         end
%         if ( isempty(n) )
%             n    = bsspartitionsize(x);
%         end
%         if ( isempty(m) )
%             m    = bssdefaultorder(m);
%         end
%         bss  = [];    case 5
%         % All inputs present and will be passed to the input checker.
%     otherwise
%         error('Incorrect number of input arguments');
% end
% 
% bsspdfcheckinputs(x,bnds,n,m,bss);

if ( isempty(bss) )
    bss  = bsscreatestruct( bnds, n, m );
end

%
% If we have a B-spline Series with data already added, we need to unscale
% the coefficients if the number of points is greater than zero, and then
% proceed as usual. If the number of points is zero, then we assume this is
% the first use of the B-Spline Series structure and proceed.
%
if ( bss.npts > 0 )
    bss.c = bss.c * bss.npts * prod(bss.partition.h);
end

%
% Do the estimation of the PDF.
%

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

inbounds = true(size(x,1),1);
for idim = 1:bss.ndim
    inbounds = inbounds & ((x(:,idim) >= bss.bnds(idim,1)) & (x(:,idim) < bss.bnds(idim,2)));
end

x      = x(inbounds,:);
npts_x = size(x,1);

xtmp = bsxfun(@minus,x,bss.bnds(:,1)')*diag(1./bss.partition.h);
xidx = bss.m + floor(xtmp);

idxp = zeros(bss.m,bss.ndim);
xpn  = zeros(bss.m,bss.ndim);

for idata = 1:npts_x
    
    for idim = 1:bss.ndim
        idxp(:,idim) = xidx(idata,idim) + mp1to0;
        xpn(:,idim)  = (bss.partition.x(idxp(:,idim),idim)-bss.bnds(idim,1))./bss.partition.h(idim);
    end
    
    xdiff = bsxfun(@minus,xtmp(idata,:),xpn);
    fbtmp = reshape(bssfbspline(xdiff(:),bss.m),bss.m,bss.ndim);
    
    tp   = bsstensorproduct(fbtmp);
    tpi  = bsstensorindex(idxp);
    tpli = bsssub2ind(bss.partition.npart-1,tpi);
    
    bss.c(tpli) = bss.c(tpli) + tp;
    
end
bss.npts = bss.npts+npts_x;

%
% Check for active boundaries
%
if ( bss.m > 1 )
    %lower bounds active
    active = xidx < 2*bss.m-1;
    for idim = 1:bss.ndim
        if( any(active(:,idim)) )
            idx          = find(active(:,idim));
            xtmp         = x(idx,:);
            n            = size(xtmp,1);
            xtmp(:,idim) = 2*bss.bnds(idim,1) - xtmp(:,idim);
            xtmp         = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            for idata = 1:n
                
                idxp  = bsxfun(@plus,xidx(idx(idata),:),mp1to0);
                idxp(:,idim) = bss.m + mp1to0;
                xpn  = (bss.partition.x(idxp)-repmat(bss.bnds(:,1)',bss.m,1))*(diag(1./bss.partition.h));
                
                xdiff = repmat(xtmp(idata,:),bss.m,1) - xpn;
                fbtmp = reshape(bssfbspline(xdiff(:),bss.m),bss.m,bss.ndim);
                
                tp   = bsstensorproduct(fbtmp);
                tpi  = bsstensorindex(idxp);
                tpli = bsssub2ind(bss.partition.npart-1,tpi);
                
                bss.c(tpli) = bss.c(tpli) + tp;
                
            end
        end
    end
    %upper bounds active
    active = xidx > repmat(bss.partition.npart'-bss.m,npts_x,1);
    for idim = 1:bss.ndim
        if( any(active(:,idim)) )
            idx          = find(active(:,idim));
            xtmp         = x(idx,:);
            n            = size(xtmp,1);
            xtmp(:,idim) = 2*bss.bnds(idim,2) - xtmp(:,idim);
            xtmp         = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            for idata = 1:n
                
                idxp  = bsxfun(@plus,xidx(idx(idata),:),mp1to0);
                idxp(:,idim) = bss.partition.npart(idim)-1 + mp1to0;
                xpn  = (bss.partition.x(idxp)-repmat(bss.bnds(:,1)',bss.m,1))*(diag(1./bss.partition.h));

                xdiff = repmat(xtmp(idata,:),bss.m,1) - xpn;
                fbtmp = reshape(bssfbspline(xdiff(:),bss.m),bss.m,bss.ndim);
                
                tp   = bsstensorproduct(fbtmp);
                tpi  = bsstensorindex(idxp);
                tpli = bsssub2ind(bss.partition.npart-1,tpi);
                
                bss.c(tpli) = bss.c(tpli) + tp;
                
            end
        end
    end
end

bss.c = bss.c./(prod(bss.partition.h)*bss.npts);

bss_out = bss;
