function r = preferential_median(x, u, n, n_min, n_max)
% PREFERENTIAL_MEDIAN  IF&PA — interval fusion with preference aggregation.
%   r = preferential_median(x, u) returns the preferential median of the
%   intervals [x(k) - u(k), x(k) + u(k)]. u may be a scalar.
%   r = preferential_median(x, u, n) sets the number of discrete values of the
%   range of actual values (RAV); pass [] for the default.
%
%   Fields of r: pm, u_pm, u_grid, u_sampling, n, h, winners (1-based),
%   converged, grid, coverage.
%
%   MATLAB R2016b+ and GNU Octave 6+, no toolboxes.
%
%   References:
%     Muravyov S.V., Khudonogova L.I., Emelyanova E.Yu. Interval data fusion with
%     preference aggregation. Measurement, 2018, 116, 621-630. doi:10.1016/j.measurement.2017.08.045
%     Muravyov S.V., Khudonogova L.I., Shabramov V.S., Ignatyev V.D., Andreyev D.I.
%     Kemeny and Borda rules in constructing central tendency estimators based on the
%     preferential median. Proc. AISP'25, 2025. doi:10.1109/AISP68263.2025.11396113
%     (the Borda rule gives the same preferential median as the Kemeny rule).
%     Muravyov S.V., Khudonogova L.I., Shabramov V.S., Ignatyev V.D. Comparative studies
%     of preferential and weighted medians accuracy in robust estimation of central
%     tendency. Int. J. Data Sci. Anal., 2026, 22(1), 262. doi:10.1007/s41060-026-01222-6
%     (standard uncertainty of the PM due to the RAV partition: 0.41 h).
%
%   Version of the method: IF&PA-L — Borda rule, no self-refinement.
%   Uncertainty: u_pm returned here is the combined standard uncertainty
%     sqrt(u_grid^2 + u_sampling^2), written u(x_pm) on the portal's Method page;
%     u_grid = 0.41 h is what the papers call u_pm; u_sampling = 1.2533 * 1.4826 * MAD / sqrt(m)
%     is an addition of this implementation, not part of the published method.
%
%   Example:
%     r = preferential_median([4.375 4.9], [0.175 0.2], 3);
%     r.pm   % 4.65

if nargin < 3, n = []; end
if nargin < 4 || isempty(n_min), n_min = 11; end
if nargin < 5 || isempty(n_max), n_max = 201; end

x = double(x(:)');
u = double(u(:)');
if isscalar(u), u = repmat(u, size(x)); end
if numel(x) ~= numel(u), error('x and u have different lengths'); end
if isempty(x), error('empty set of intervals'); end
if any(u < 0), error('uncertainties must be non-negative'); end
if ~all(isfinite([x u])), error('NaN or infinity in data'); end

lo = x - u;
hi = x + u;
a1 = min(lo);
an = max(hi);
centres = (lo + hi) / 2;

if an - a1 <= 0
    r = struct('pm', a1, 'u_pm', 0, 'u_grid', 0, 'u_sampling', 0, 'n', 1, ...
               'h', 0, 'winners', 1, 'converged', true, 'grid', a1, ...
               'coverage', numel(lo));
    return
end

if isempty(n)
    width_min = min(hi - lo);            % the narrowest interval, 2 min u_k
    if width_min > 0
        n = ceil((an - a1) / width_min);
    else
        n = n_min;
    end
    n = max(n, n_min);
end
n = max(floor(n), 2);

while true
    h = (an - a1) / (n - 1);
    grid = a1 + h * (0:n-1);
    % tolerance in the units of the data: a share of the grid step and a few ulps of the values
    epsv = max(1e-9 * h, 8 * eps * max(abs(a1), abs(an)));
    % coverage c_i; Borda score z_i = const + n c_i, so winners = max coverage
    % c_i = #{lo_k - eps <= a_i} - #{hi_k + eps < a_i}, by a stable sort of bounds together with the grid:
    % O((n + m) log(n + m)) memory O(n + m) instead of an m-by-n matrix
    coverage = (count_before(lo(:) - epsv, grid(:), true) - count_before(hi(:) + epsv, grid(:), false)).';
    idx = find(coverage == max(coverage));
    adjacent = numel(idx) <= 1 || all(diff(idx) == 1);
    if adjacent || n >= n_max
        if ~adjacent
            % widest block of consecutive winners; tie -> median closest to centres
            starts = [1, find(diff(idx) ~= 1) + 1];
            stops = [starts(2:end) - 1, numel(idx)];
            centre = median(centres);
            best = []; best_d = inf;
            for j = 1:numel(starts)
                b = idx(starts(j):stops(j));
                d = abs(median(grid(b)) - centre);
                if numel(b) > numel(best) || (numel(b) == numel(best) && d < best_d)
                    best = b; best_d = d;
                end
            end
            idx = best;
        end
        u_grid = 0.41 * h;
        m = numel(centres);
        if m < 2
            u_samp = 0;
        else
            u_samp = 1.2533 * 1.4826 * median(abs(centres - median(centres))) / sqrt(m);
        end
        r = struct('pm', median(grid(idx)), 'u_pm', hypot(u_grid, u_samp), ...
                   'u_grid', u_grid, 'u_sampling', u_samp, 'n', n, 'h', h, ...
                   'winners', idx, 'converged', adjacent, 'grid', grid, ...
                   'coverage', coverage);
        return
    end
    n = n + 1;
end
end


function c = count_before(b, g, inclusive)
% For each g(i) of an ascending grid: how many b are <= g(i) (inclusive) or < g(i) (strict).
% sort is stable, so on ties the order of concatenation decides which side counts.
if inclusive
    v = [b; g]; isb = [true(numel(b), 1); false(numel(g), 1)];
else
    v = [g; b]; isb = [false(numel(g), 1); true(numel(b), 1)];
end
[~, o] = sort(v);
cnt = cumsum(isb(o));
c = cnt(~isb(o));
end