function [iqr,percentiles] = datafuniqr(x,lowup,dim,omitnan)
% DATAFUNIQR Compute interquartile range for histogram, normalize and outlier functions
% Inputs:
%   x       - data
%   lowup   - lower and upper bounds (default is [25 75])
%   dim     - dimension argument, must be a scalar (default is 1)
%   omitnan - true/false (default is false)

%   Copyright 1984-2020 The MathWorks, Inc.

% if not passed, use default
if nargin < 2
    lowup = [25 75];
end
if nargin < 3
    dim = 1;
end
if nargin < 4
    omitnan = false;
end
sz = size(x);
if dim > 1
    % permute and reshape so we always act on columns
    dim = min(dim, ndims(x)+1);
    sz(end+1:dim) = 1;
    sizeiqr = sz;
    sizeiqr(dim) = 1;
    perm = [dim setdiff(1:max(ndims(x),dim),dim)];
    x = permute(x, perm);
    x = reshape(x, [sz(dim), prod(sizeiqr)]);
else
    sizeiqr = sz;
    sizeiqr(dim) = 1;
end
lowup = double(lowup)/100;

if ~isfloat(x)
    x = double(x);
end

x = sort(x,1);
if omitnan && ~isempty(x)
    numNonNaN = sum(~isnan(x),1);
    numNonNaN(numNonNaN == 0) = 1;
    if all(numNonNaN == numNonNaN(1),'all')
        % We need this indexing in case x is ND and dim = 1
        ind = repmat({':'},ndims(x)-1,1);
        x = x(1:numNonNaN(1),ind{:});
        [iqr, percentiles] = iqrImpl(x, lowup);
    else
        iqr = zeros(sizeiqr);
        percentiles = zeros(2*sizeiqr(1),sizeiqr(2));
        for j = 1:prod(size(x,2:ndims(x)))
            [iqr(j), percentiles(:,j)] = iqrImpl(x(1:numNonNaN(j),j), lowup);
        end
    end
else
    [iqr, percentiles] = iqrImpl(x, lowup);
end
if dim > 1
    iqr = reshape(iqr, sizeiqr);
end
end

function [iqr, percentiles] = iqrImpl(x, lowup)

n = size(x,1);
F = ((1:n)'-.5) / n;

if ~isempty(x) && (n > 1)
    if ismatrix(x)
        % interp1q for 2-D because it supports a sparse second input.
        % max/min allows to deal with the case when given percentiles are outside of the [F(1)..F(n)] range
        percentiles = interp1q(F, x, [max(lowup(1),F(1)); min(lowup(2),F(n))]);
    else
        % Multi-valued interp1 for columns of n-D array.
        percentiles = interp1(F, x, [max(lowup(1),F(1)); min(lowup(2),F(n))]);
    end
    iqr = diff(percentiles,1,1);
else
    xsiz = size(x);
    percentiles = NaN([2 xsiz(2:end)],'like',x);
    iqr = NaN([1 xsiz(2:end)],'like',x);
end

end
