function bss_out = bsspdfest2d(x, bnds, n, m, bss)
%BSSPDFEST2D - Estimate a 2-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. Values of X outside the estimation
%       bounds are ignored.
%
%       The number of columns NDIM must be 2.
%
%    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 2.
%
%       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. Reflection based
%       corrections are also made for corners where multiple boundaries are
%       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 2.
%
%       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-08
% Author  : Kevin R. Gehringer
%
%     Changed the method for dealing with bounded domains from using the
%     histogram value at the boundary to using reflection of the data
%     adjacent to the boundary to its other side and into corners and
%     through edges.
%
%     Replaced calls to remat with calls to bsxfun to improve performance.
%
%     Changed the way the coefficient matrix was updated from a point by
%     point summation to a chunked summation based on a sorted list of the
%     unique 2-dimensional indices into the coefficient matrix. This
%     provided a significant speed improvement.
%
% Modified: 2016-10-22
% Author : Kevin R. Gehringer
%
%    Added a check for active bounds on the estimation interval. This is
%    detected by checking to see if any of the bss.m-1 coefficients at
%    either boundary are greater than zero. If so, a histogram is estimate
%    is generated for the bin adjacent to the boundary and used to specify
%    the PDF value at the boundary. The B-spline series coefficients are
%    then modified to produce this value at the boundary.
%
%    Data points that are not in the estimation interval are now ignored.
%    Previously they were evaluated, but always returned zero values, that
%    would stack up on the boundaries. This caused several unnecessary
%    complications when computing boundary values, leading to the change in
%    behavior.
%
% 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.
%   4) Changed the matrix orientation from m by n_p to n_p by m, where m is
%      the B-spline order and n_p is the number of input points for
%      partition interval p.
%   5) Explicitly create indexes and B-spline evaluation matrices for each
%      dimension.
%
% Created: 2014-11-25
% Author : Kevin R. Gehringer
%

%
% Check to be sure that the data are two dimensional.
%

if (size(x,2) ~= 2)
    error('The data must be 2-dimensional with the points as rows.');
end

% 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);

xpn    = (reshape(bss.partition.x(bsxfun(@plus,xidx(:,1),mp1to0),1),npts_x,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),npts_x,bss.m);

xpn    = (reshape(bss.partition.x(bsxfun(@plus,xidx(:,2),mp1to0),2),npts_x,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),npts_x,bss.m);

[u,lidx,uidx,sidx] = bssuniqueidx(xidx);
for iidx = 1:length(lidx)
    idx1 = u(iidx,1)+mp1to0;
    idx2 = u(iidx,2)+mp1to0;
    bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
end

bss.npts = bss.npts+npts_x;

%
% See if we need to do a boundary correction. If so, we estimate the
% boundary value using reflection across the active boundaries and corners.
%

if ( bss.m > 1 )
    lb_active = xidx < 2*bss.m-1;
    ub_active = xidx > repmat(bss.partition.npart'-bss.m,npts_x,1);%change this to ues bsxfun
    if ( any(any([lb_active ub_active],2)) )
        lb = bss.partition.x(bss.m,:);
        ub = [bss.partition.x(bss.partition.npart(1),1) bss.partition.x(bss.partition.npart(2),2)];
        if (any(lb_active(:,1)))
            xtmp = [2*lb(1)-x(lb_active(:,1),1) x(lb_active(:,1),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bsxfun(@plus,bss.m+floor(xtmp(:,2)),mp1to0);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);

            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(ub_active(:,1)))
            xtmp = [2*ub(1)-x(ub_active(:,1),1) x(ub_active(:,1),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bss.partition.npart(1)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bsxfun(@plus,bss.m+floor(xtmp(:,2)),mp1to0);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);

            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(lb_active(:,2)))
            xtmp = [x(lb_active(:,2),1) 2*lb(2)-x(lb_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bsxfun(@plus,bss.m+floor(xtmp(:,1)),mp1to0);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(ub_active(:,2)))
            xtmp = [x(ub_active(:,2),1) 2*ub(2)-x(ub_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bsxfun(@plus,bss.m+floor(xtmp(:,1)),mp1to0);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.partition.npart(2)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        
        if (any(lb_active(:,1) & lb_active(:,2)))
            xtmp = [2*lb(1)-x(lb_active(:,1) & lb_active(:,2),1) 2*lb(2)-x(lb_active(:,1) & lb_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);

            idxp1  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(lb_active(:,1) & ub_active(:,2)))
            xtmp = [2*lb(1)-x(lb_active(:,1) & ub_active(:,2),1) 2*ub(2)-x(lb_active(:,1) & ub_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.partition.npart(2)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(ub_active(:,1) & lb_active(:,2)))
            xtmp = [2*ub(1)-x(ub_active(:,1) & lb_active(:,2),1) 2*lb(2)-x(ub_active(:,1) & lb_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bss.partition.npart(1)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.m + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
        if (any(ub_active(:,1) & ub_active(:,2)))
            xtmp = [2*ub(1)-x(ub_active(:,1) & ub_active(:,2),1) 2*ub(2)-x(ub_active(:,1) & ub_active(:,2),2)];
            n = size(xtmp,1);
            xtmp = bsxfun(@minus,xtmp,bss.bnds(:,1)')*diag(1./bss.partition.h);
            
            idxp1  = bss.partition.npart(1)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp1,1),n,bss.m) - bss.bnds(1,1))/bss.partition.h(1);
            xdiff  = bsxfun(@minus,xtmp(:,1),xpn);
            fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            idxp2  = bss.partition.npart(2)-1 + repmat(mp1to0,n,1);
            xpn    = (reshape(bss.partition.x(idxp2,2),n,bss.m) - bss.bnds(2,1))/bss.partition.h(2);
            xdiff  = bsxfun(@minus,xtmp(:,2),xpn);
            fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);
            
            [u,lidx,uidx,sidx] = bssuniqueidx([idxp1(:,bss.m) idxp2(:,bss.m)]);
            for iidx = 1:length(lidx)
                idx1 = u(iidx,1)+mp1to0;
                idx2 = u(iidx,2)+mp1to0;
                bss.c(idx1,idx2) = bss.c(idx1,idx2) + (fbtmp1(sidx(lidx(iidx):uidx(iidx)),:)')*fbtmp2(sidx(lidx(iidx):uidx(iidx)),:);
            end
        end
    end
end

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

bss_out = bss;
