function ptsperpart = bsspointsperpartition(ndim,action,n)
% BSSPOINTSPERPARTITION return the number of evaluation points per
%    partition to generate the interpolation grid for a B-Spline series.
%
% INPUTS:
%
%    ndim:
%       The dimension of the data. Only four dimensions are recognized: 1,
%       2, 3, and 4 or more.
%
%    action:
%       The action to take for dimension size NDIM. Allowed actions are
%       'set' to set the points per partition value for dimension size NDIM
%       to N and 'default' to set the points per partition value for
%       dimension size NDIM its default value.
%
%       This argument is optional. If not present the current value for the
%       number of evaluation points per partition is returned.
%
%       If ACTION == set the value N is assigned dimension size NDIM is
%       returned.
%
%    n:
%       The points per partition value to use for dimension size NDIM. This
%       is required if ACTION == 'set'.
%
% OUTPUTS:
%
%    ptsperpart:
%       The number of points per partition to use to use when evaluating
%       tha B-Spline series to generate a grid for interpolation. This is a
%       scalar value that is used for each dimension in the partition, but
%       the partition in each dimension may have different sizes.
%

%
% Copyright 2014-2023 Biometrics Northwest LLC
%

%
% Modified: 2023-03-06
% Author : Kevin R. Gehringer
%
%    Changed the default values for points per partition, making them
%    smaller for faster evaluation. This had minimal impact on the
%    estimation errors.
%
% Created: 2016-04-24
% Author : Kevin R. Gehringer
%

    persistent npts_default
    persistent npts
    
    if (isempty(npts_default))
%        npts_default = [25 3 1.25 1.125];
%        npts_default = [25 5 2.5 2.125];
%        npts_default = [10 5 2.5 2.125];
        npts_default = [5 3.25 2.5 2.125];
        npts = npts_default;
    end

    if ( exist ('ndim','var') )
        if ( ~isnumeric(ndim) )
             error('The number of dimensions was not numeric.');
        end
        ndim = floor(ndim);
        if (ndim <1 )
            error('The number of dimensions must be at least 1.');
        end
    end
    
    if ( exist ('action','var') )
        switch lower(action)
            case 'set'
                if ( ~exist('n','var') )
                    error('Not enough input arguments.');
                end
                if ( ~isnumeric(n) )
                    error('The number of points was not numeric.');
                end
                nf = floor(n);
                if (nf <= 0)
                    error('The number of points must be at least 1.');
                end
                if ( (nf > 5) && (ndim > 1) )
                    error('The number of evaluation points is large and evaluation would take a very long time. Consider a maximum value of 5.');
                end
                if ( ndim <= 3)
                    npts(ndim) = nf;
                else
                    npts(4) = nf;
                end
            case 'default'
                if ( exist('n','var') )
                    error('Too many input arguments.');
                end
                if ( ndim <= 3)
                    npts(ndim) = npts_default(ndim);
                else
                    npts(4) = npts_default(4);
                end
            otherwise
                error('Unrecognized action ''%s''. Allowed values are ''set'' and ''default''',action);
        end
    end
    
    if ( ndim <= 3)
        ptsperpart = npts(ndim);
    else
        ptsperpart = npts(4);
    end
    
end