function [y,xg] = bssgridevalc(xgrid,bss,evalfun)
% BSSGRIDEVALC - Generate the grid points and evaluate for the desired
%    evaluation function and B-spline series using a compact grid and
%    evaluation point subsets to conserve memory.
%
% INPUTS:
%
%    xgrid:
%       A cell array containing grid vectors defining a compact grid used
%       for a gridded interpolant.
%
%    bss:
%       The B-Spline series data structure containing the series to be used
%       for evaluation.
%
% OUTPUTS:
%
%    y:
%       The values for the specified evaluation function at the grid points
%       defined by XGRID.
%
%       BSSGRIDEVALC generates and evaluates the points in the same order
%       as NDGRID generates the points except that it generates the gridded
%       data points in subsets, as rows in a matrix, iterating over the
%       dimensions to conserve memory, rather than producing the full grid
%       all at once as separate variables or as rows in a matrix and
%       evaluating them.
%

%
% Copyright 2023 Biometrics Northwest LLC
%

%
% Modified: 2023-05-13
% Author  : Kevin R. Gehringer
%
%    Corrected the algorithm so that it works for all dimension sizes.
%
% Created: 2023-02-18
% Author : Kevin R. Gehringer
%

if ( length(xgrid) ~= bss.ndim )
    error('Dimension mismatch: Dimensions of the compact grid and B-spline series did not agree.')
end

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

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

    if ( nargout == 2 )
        xg = zeros(npoints,ndims);
    end

    y = zeros(npoints,1);
    step = ndim(1);
    ibeg = 1;
    idxdim = ones(1,ndims);
    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;

        % assign the grid values
        for i = 2:ndims
            xtmp(:,i) = xgrid{i}(idxdim(i));
        end

        % compute the grid function
        switch( ndims )
            case 2
                y(ibeg:iend) = bsseval2d(xtmp,bss,'pdf');
            case 3
                y(ibeg:iend) = bsseval3d(xtmp,bss,'pdf');
            otherwise
                y(ibeg:iend) = bssevalnd(xtmp,bss,'pdf');
        end

        % increment idxdim here
        for i = 2:ndims
            if ( (idxdim(i) == ndim(i)) && (i<ndims) )
                idxdim(i) = 1;
                if (idxdim(i+1) < ndim(i+1) )
                    idxdim(i+1) = idxdim(i+1) + 1;
                    break
                end
            else
                idxdim(i) = idxdim(i) + 1;
                break
            end
        end

        if ( nargout == 2 )
            xg(ibeg:iend,:) = xtmp;
        end

        ibeg = iend + 1;
    end

    if ( strcmpi(evalfun,'cdf') )
        h = 1;
        y = reshape(y,ndim);
        for idim = 1:bss.ndim
            y = cumsum(y,idim);
            h = h*(xgrid{idim}(2) - xgrid{idim}(1));
        end
        y = y(:) * h;
    end

else
    if (iscolumn(xgrid{1}) )
        y = bsseval1d(xgrid{1},bss,evalfun);
        if (nargout == 2 )
            xg = xgrid{1};
        end
    else
        y = bsseval1d(xgrid{1}',bss,evalfun);
        if (nargout == 2 )
            xg = xgrid{1}';
        end
    end
end

end
