function y = bsseval3d(x, bss, evalfun)
%BSSEVAL3D  Evaluate the 3-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 matrix with NDIM = 3.
%
%    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 if the B-Spline series represents a
%       probability density function:
%
%          pdf      : The probability density function, f(x) (default)
%          cdf      : The cumulative distribution function, F(x)
%          survivor : The survivor function, S(x) = 1-F(x)
%
% 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.
%
%    Added a function Do_PDF_CG to perform the evaluation of the PDF based
%    on a compact grid to conserve memory when building the grid for the
%    CDF.
%
% Modified: 2018-09-10
% Author : Kevin R. Gehringer
%
%    Fixed a bug in the extrapolation of the values for the CDF in do_cdf.
%    There was an incorrect test of the wrong index for dims 1 and 3  >=
%    bnds([1 3],2) and dim 2 < bnds(2,2). This caused the extrapolated
%    values to be zero.
%
% 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-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 reshaping for B-splines of order 1.
%   3) 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.
%   4) Removed unnecessary temporary variables.
%   5) Explicitly create indexes and B-spline matrices for each dimension.
%
% 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', 'survivor'.
%
    okfuns = {'pdf';'cdf';'survivor'};
    if ( ~any(strcmpi(evalfun,okfuns)) )
        error('The evaluation function ''%s'' was not recognized.',evalfun);
    end
    
%
%   Call the appropriate function.
%
    switch (lower(evalfun))
        case {'pdf'}
            y = do_pdf(x,bss);
        case {'cdf'}
            y = do_cdf(x,bss);
        case {'survivor'}
            y = do_survivor(x,bss);
        otherwise
            error('The function ''%s'' was not recognized.',evalfun);
    end
    
end

function y = do_bss(x,bss)

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

npts_x  = size(x,1);

y = zeros(npts_x,1);

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

if ( n == 0 )
    return;
end

xtmp = (x(inbounds,:) - repmat(bss.bnds(:,1)',n,1))*diag(1./bss.partition.h);
xidx = bss.m + floor(xtmp);

idxp1  = repmat(xidx(:,1),1,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  = repmat(xtmp(:,1),1,bss.m) - xpn;
fbtmp1 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);

idxp2  = repmat(xidx(:,2),1,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  = repmat(xtmp(:,2),1,bss.m) - xpn;
fbtmp2 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);

idxp3  = repmat(xidx(:,3),1,bss.m)+repmat(mp1to0,n,1);
xpn    = (reshape(bss.partition.x(idxp3,3),n,bss.m) - bss.bnds(3,1))/bss.partition.h(3);
xdiff  = repmat(xtmp(:,3),1,bss.m) - xpn;
fbtmp3 = reshape(bssfbspline(xdiff(:),bss.m),n,bss.m);

mndim = bss.m*ones(1,3);
msquared = bss.m*bss.m;

for idata = 1:n

    tmp3  = repmat(fbtmp3(idata,:),msquared,1);
    tmp12  = (fbtmp1(idata,:)')*fbtmp2(idata,:);
    tmp123 = tmp3(:).*repmat(tmp12(:),bss.m,1);
    ytmp  = bss.c(idxp1(idata,:),idxp2(idata,:),idxp3(idata,:)) .* reshape(tmp123,mndim);
    y(inbounds(idata)) = sum(ytmp(:));

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 = do_pdf_cg(xgrid,bss)

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

if ( ~do_interp )
    ndims = length(xgrid);
    ndim = zeros(1,ndims);
    for i = 1:ndims
        ndim(i) = length(xgrid{i});
    end

    nblocks = prod(ndim(2:end));
    npoints = prod(ndim);

    y = zeros(npoints,1);
    step = ndim(1);
    ibeg = 1;
    idxdim = ones(1,ndims-1);
    if ( iscolumn(xgrid{1}) )
        xtmp = [xgrid{1} zeros(ndim(1),ndims-1)];
    else
        xtmp = [xgrid{1}' zeros(ndim(1),ndims-1)];
    end

    for ib = 1:nblocks
        iend = ibeg + step - 1;

        for i = 2:ndims
            xtmp(:,i) = xgrid{i}(idxdim(i-1));
        end

        if ( any(idxdim == ndim(2:end)) )
            idxeq = find(idxdim == ndim(2:end));
            if (max(idxeq) < ndims-1)
                idxdim(idxeq) = 1;
                idxdim(idxeq+1) = idxdim(idxeq+1) + 1;
            else
                idxdim(1) = idxdim(1) + 1;
            end
        else
            idxdim(1) = idxdim(1) + 1;
        end

        y(ibeg:iend) = do_bss(xtmp,bss);

        ibeg = iend + 1;
    end

else
    y = bss.pdf(xgrid);
    y(y<0) = 0;
end

end

function y = do_cdf(x,bss)

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

    do_interp = isfield(bss,'cdf') && ~isempty(bss.cdf);
    
    if ( ~do_interp )
        
        np = floor(bssdefaultpointsperpartition(bss.ndim)*bss.partition.n)+1;
    
        xcdf1 = linspace(bss.bnds(1,1),bss.bnds(1,2),np(1))';
        xcdf2 = linspace(bss.bnds(2,1),bss.bnds(2,2),np(2))';
        xcdf3 = linspace(bss.bnds(3,1),bss.bnds(3,2),np(3))';
    
        ndim = [ length(xcdf1) length(xcdf2) length(xcdf3)];
        h = [xcdf1(2)-xcdf1(1) xcdf2(2)-xcdf2(1) xcdf3(2)-xcdf3(1)];
    
        %ycdf = do_pdf(bssgridpoints({xcdf1; xcdf2; xcdf3}),bss);
        ycdf = do_pdf_cg({xcdf1; xcdf2; xcdf3},bss);
        ycdf = reshape(ycdf,ndim);
        for idim = 1:bss.ndim
            ycdf = cumsum(ycdf,idim);
        end
        ycdf = prod(h)*ycdf;
%
%       fix any roundoff errors that may have made the CDF greater than one.
%
        ycdf(ycdf>1) = 1;
%
%       do the interpolation for points that are within the bounds
%
        f = griddedInterpolant({xcdf1; xcdf2; xcdf3},ycdf,'linear');
        y(inbounds) = f(x(inbounds,:));
    else
        y(inbounds) = bss.cdf(x(inbounds,:));
    end
%
%   fix values for out of bounds points that may make sense. This basically
%   performs extrapolations of the edge values at the boundary of the
%   B-Spline series estimation domain.
%
%   y((x(:,1) < bss.bnds(1,1)) | (x(:,2) < bss.bnds(2,1))) = 0;
    y((x(:,1) >= bss.bnds(1,2)) & (x(:,2) >= bss.bnds(2,2)) & (x(:,3) >= bss.bnds(3,2))) = 1;
    
    idx1gt23lt = find((x(:,1) >= bss.bnds(1,2)) & (x(:,2) < bss.bnds(2,2)) & (x(:,3) < bss.bnds(3,2)));
    if ( ~isempty(idx1gt23lt) )
%        y(idx1gt23lt) = interp3(xx2,xx1,xx3,ycdf,bss.bnds(1,2)*ones(length(idx1gt23lt),1),x(idx1gt23lt,2),x(idx1gt23lt,3),'linear',0);
        y(idx1gt23lt) = f(bss.bnds(1,2)*ones(length(idx1gt23lt),1),x(idx1gt23lt,2),x(idx1gt23lt,3));
    end
    
    idx2gt13lt = find((x(:,1) < bss.bnds(1,2)) & (x(:,2) >= bss.bnds(2,2)) & (x(:,3) < bss.bnds(3,2)));
    if ( ~isempty(idx2gt13lt) )
%        y(idx2gt13lt) = interp3(xx2,xx1,xx3,ycdf,x(idx2gt13lt,1),bss.bnds(2,2)*ones(length(idx2gt13lt),1),x(idx2gt13lt,3),'linear',0);
        y(idx2gt13lt) = f(x(idx2gt13lt,1),bss.bnds(2,2)*ones(length(idx2gt13lt),1),x(idx2gt13lt,3));
    end

    idx3gt12lt = find((x(:,1) < bss.bnds(1,2)) & (x(:,2) < bss.bnds(2,2)) & (x(:,3) >= bss.bnds(3,2)));
    if ( ~isempty(idx3gt12lt) )
%        y(idx3gt12lt) = interp3(xx2,xx1,xx3,ycdf,x(idx3gt12lt,1),x(idx3gt12lt,2),bss.bnds(3,2)*ones(length(idx3gt12lt),1),'linear',0);
        y(idx3gt12lt) = f(x(idx3gt12lt,1),x(idx3gt12lt,2),bss.bnds(3,2)*ones(length(idx3gt12lt),1));
    end

    idx12gt3lt = find((x(:,1) >= bss.bnds(1,2)) & (x(:,2) >= bss.bnds(2,2)) & (x(:,3) < bss.bnds(3,2)));
    if ( ~isempty(idx12gt3lt) )
%        y(idx12gt3lt) = interp3(xx2,xx1,xx3,ycdf,bss.bnds(1,2)*ones(length(idx12gt3lt),1),bss.bnds(2,2)*ones(length(idx12gt3lt),1),x(idx12gt3lt,3),'linear',0);
        y(idx12gt3lt) = f(bss.bnds(1,2)*ones(length(idx12gt3lt),1),bss.bnds(2,2)*ones(length(idx12gt3lt),1),x(idx12gt3lt,3));
    end

    idx13gt2lt = find((x(:,1) >= bss.bnds(1,2)) & (x(:,2) < bss.bnds(2,2)) & (x(:,3) >= bss.bnds(3,2)));
    if ( ~isempty(idx13gt2lt) )
%        y(idx13gt2lt) = interp3(xx2,xx1,xx3,ycdf,bss.bnds(1,2)*ones(length(idx13gt2lt),1),x(idx13gt2lt,2),bss.bnds(3,2)*ones(length(idx13gt2lt),1),'linear',0);
        y(idx13gt2lt) = f(bss.bnds(1,2)*ones(length(idx13gt2lt),1),x(idx13gt2lt,2),bss.bnds(3,2)*ones(length(idx13gt2lt),1));
    end

    idx23gt1lt = find((x(:,1) < bss.bnds(1,2)) & (x(:,2) >= bss.bnds(2,2)) & (x(:,3) >= bss.bnds(3,2)));
    if ( ~isempty(idx23gt1lt) )
%        y(idx23gt1lt) = interp3(xx2,xx1,xx3,ycdf,x(idx23gt1lt,1),bss.bnds(2,2)*ones(length(idx23gt1lt),1),bss.bnds(3,2)*ones(length(idx23gt1lt),1),'linear',0);
        y(idx23gt1lt) = f(x(idx23gt1lt,1),bss.bnds(2,2)*ones(length(idx23gt1lt),1),bss.bnds(3,2)*ones(length(idx23gt1lt),1));
    end
    
    y(y<0) = 0;
    y(y>1) = 1;
    
end

function y = do_survivor(x,bss)

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