// IF&PA — interval fusion with preference aggregation.
// Preferential median (PM) of a set of intervals [x_k - u_k, x_k + u_k].
//
// Header-only, C++17, standard library only.
//
// 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.
//
// Usage:
//   #include "ifpa.hpp"
//   auto r = ifpa::preferential_median({4.375, 4.9}, {0.175, 0.2}, 3);
//   r.pm;   // 4.65
#pragma once

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <limits>
#include <stdexcept>
#include <utility>
#include <vector>

namespace ifpa {

struct Result {
    double pm = 0, u_pm = 0, u_grid = 0, u_sampling = 0, h = 0;
    int n = 0;
    std::vector<int> winners;     // 0-based indices into grid
    bool converged = true;
    std::vector<double> grid;
    std::vector<int> coverage;
};

inline double median(std::vector<double> v) {
    if (v.empty()) throw std::invalid_argument("empty sample");
    std::sort(v.begin(), v.end());
    const std::size_t p = v.size();
    return p % 2 ? v[(p - 1) / 2] : (v[p / 2 - 1] + v[p / 2]) / 2.0;
}

// n <= 0 means "choose automatically"
inline Result preferential_median(const std::vector<double>& x, const std::vector<double>& u,
                                  int n = 0, int n_min = 11, int n_max = 201) {
    if (x.size() != u.size()) throw std::invalid_argument("x and u have different lengths");
    if (x.empty()) throw std::invalid_argument("empty set of intervals");
    const std::size_t m = x.size();
    std::vector<double> lo(m), hi(m), centres(m);
    for (std::size_t k = 0; k < m; ++k) {
        if (!(u[k] >= 0)) throw std::invalid_argument("uncertainties must be non-negative");
        if (!std::isfinite(x[k]) || !std::isfinite(u[k])) throw std::invalid_argument("NaN or infinity in data");
        lo[k] = x[k] - u[k];
        hi[k] = x[k] + u[k];
        centres[k] = (lo[k] + hi[k]) / 2.0;
    }
    const double a1 = *std::min_element(lo.begin(), lo.end());
    const double an = *std::max_element(hi.begin(), hi.end());

    Result r;
    if (an - a1 <= 0) {
        r.pm = a1; r.n = 1; r.winners = {0}; r.grid = {a1}; r.coverage = {int(m)};
        return r;
    }
    if (n <= 0) {
        double width_min = std::numeric_limits<double>::infinity();
        for (std::size_t k = 0; k < m; ++k) width_min = std::min(width_min, hi[k] - lo[k]);
        n = width_min > 0 ? int(std::ceil((an - a1) / width_min)) : n_min;
        n = std::max(n, n_min);
    }
    n = std::max(n, 2);

    const double centre = median(centres);
    for (;;) {
        const double h = (an - a1) / (n - 1);
        std::vector<double> grid(n);
        for (int i = 0; i < n; ++i) grid[i] = a1 + h * i;
        // tolerance in the units of the data: a share of the grid step and a few ulps of the values
        const double eps = std::max(1e-9 * h, 8 * std::numeric_limits<double>::epsilon() * std::max(std::fabs(a1), std::fabs(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}; two pointers: O(n + m log m) instead of O(n m)
        std::vector<double> left(m), right(m);
        for (std::size_t k = 0; k < m; ++k) { left[k] = lo[k] - eps; right[k] = hi[k] + eps; }
        std::sort(left.begin(), left.end());
        std::sort(right.begin(), right.end());
        std::vector<int> coverage(n, 0);
        for (std::size_t i = 0, p = 0, q = 0; i < static_cast<std::size_t>(n); ++i) {
            while (p < m && left[p] <= grid[i]) ++p;
            while (q < m && right[q] < grid[i]) ++q;
            coverage[i] = static_cast<int>(p - q);
        }
        const int best = *std::max_element(coverage.begin(), coverage.end());
        std::vector<int> idx;
        for (int i = 0; i < n; ++i) if (coverage[i] == best) idx.push_back(i);
        bool adjacent = true;
        for (std::size_t j = 1; j < idx.size(); ++j) if (idx[j] - idx[j - 1] != 1) adjacent = false;

        if (adjacent || n >= n_max) {
            if (!adjacent) {
                // widest block of consecutive winners; tie -> median closest to centres
                std::vector<int> bestBlock, cur{idx[0]};
                double bestDist = std::numeric_limits<double>::infinity();
                auto consider = [&](const std::vector<int>& b) {
                    std::vector<double> vals;
                    for (int i : b) vals.push_back(grid[i]);
                    const double d = std::fabs(median(vals) - centre);
                    if (b.size() > bestBlock.size() || (b.size() == bestBlock.size() && d < bestDist)) {
                        bestBlock = b; bestDist = d;
                    }
                };
                for (std::size_t j = 1; j < idx.size(); ++j) {
                    if (idx[j] == cur.back() + 1) cur.push_back(idx[j]);
                    else { consider(cur); cur = {idx[j]}; }
                }
                consider(cur);
                idx = bestBlock;
            }
            std::vector<double> win;
            for (int i : idx) win.push_back(grid[i]);
            r.pm = median(win);
            r.u_grid = 0.41 * h;
            if (m >= 2) {
                std::vector<double> dev(m);
                for (std::size_t k = 0; k < m; ++k) dev[k] = std::fabs(centres[k] - centre);
                r.u_sampling = 1.2533 * 1.4826 * median(dev) / std::sqrt(double(m));
            }
            r.u_pm = std::hypot(r.u_grid, r.u_sampling);
            r.n = n; r.h = h; r.winners = idx; r.converged = adjacent;
            r.grid = grid; r.coverage = coverage;
            return r;
        }
        ++n;
    }
}

// ---------------------------------------------------------------------------
// Alternatives for comparison: common practice for combining x_k +/- u_k
// ---------------------------------------------------------------------------
struct Estimate {
    const char* name;
    double value;
    double u;   // NaN where no standard uncertainty formula exists
};

// 0.95 quantile of chi-square, Wilson-Hilferty approximation (error < 1 %)
inline double chi2_quantile_95(int df) {
    const double a = 2.0 / (9.0 * df);
    return df * std::pow(1.0 - a + 1.6448536269514722 * std::sqrt(a), 3);
}

// Weighted median: smallest x with cumulative weight >= half; exactly half -> midpoint
inline double weighted_median(const std::vector<double>& x, const std::vector<double>& w) {
    std::vector<std::size_t> o(x.size());
    for (std::size_t k = 0; k < o.size(); ++k) o[k] = k;
    std::stable_sort(o.begin(), o.end(), [&](std::size_t a, std::size_t b) { return x[a] < x[b]; });
    double total = 0;
    for (double v : w) total += v;
    double cum = 0;
    for (std::size_t j = 0; j < o.size(); ++j) {
        cum += w[o[j]];
        if (std::fabs(cum - total / 2) <= 1e-12 * total && j + 1 < o.size()) return (x[o[j]] + x[o[j + 1]]) / 2.0;
        if (cum >= total / 2) return x[o[j]];
    }
    return x[o.back()];
}

// Hodges-Lehmann estimate: median of all pairwise (Walsh) averages, i <= j
inline double hodges_lehmann(const std::vector<double>& x) {
    std::vector<double> w;
    for (std::size_t i = 0; i < x.size(); ++i)
        for (std::size_t j = i; j < x.size(); ++j) w.push_back((x[i] + x[j]) / 2.0);
    return median(w);
}

// Robust mean by ISO 13528 Algorithm A; returns {x*, u = 1.25 s*/sqrt(m)}
inline std::pair<double, double> algorithm_a(const std::vector<double>& x, double tol = 1e-10, int max_iter = 1000) {
    const std::size_t m = x.size();
    double xs = median(x);
    std::vector<double> dev(m);
    for (std::size_t k = 0; k < m; ++k) dev[k] = std::fabs(x[k] - xs);
    double ss = 1.483 * median(dev);
    if (ss == 0 || m < 2) return {xs, 0.0};
    for (int it = 0; it < max_iter; ++it) {
        const double d = 1.5 * ss;
        std::vector<double> c(m);
        double s = 0;
        for (std::size_t k = 0; k < m; ++k) { c[k] = std::min(std::max(x[k], xs - d), xs + d); s += c[k]; }
        const double nx = s / double(m);
        double q = 0;
        for (double v : c) q += (v - nx) * (v - nx);
        const double ns = 1.134 * std::sqrt(q / double(m - 1));
        const bool done = std::fabs(nx - xs) <= tol * std::max(std::fabs(xs), 1e-300) && std::fabs(ns - ss) <= tol * ss;
        xs = nx; ss = ns;
        if (done || ss == 0) break;
    }
    return {xs, 1.25 * ss / std::sqrt(double(m))};
}

// Largest consistent subset (after Cox, 2007), sequential exclusion by chi-square at 0.05
inline std::pair<double, double> largest_consistent_subset(const std::vector<double>& x, const std::vector<double>& u) {
    std::vector<std::size_t> keep(x.size());
    for (std::size_t k = 0; k < keep.size(); ++k) keep[k] = k;
    for (;;) {
        double sw = 0, swx = 0;
        for (std::size_t k : keep) { const double w = 1.0 / (u[k] * u[k]); sw += w; swx += w * x[k]; }
        const double y = swx / sw;
        double chi2 = 0;
        for (std::size_t k : keep) chi2 += (x[k] - y) * (x[k] - y) / (u[k] * u[k]);
        if (keep.size() <= 2 || chi2 <= chi2_quantile_95(int(keep.size()) - 1)) return {y, 1.0 / std::sqrt(sw)};
        const double uy2 = 1.0 / sw;
        std::size_t worst = 0;
        double worstD = -1;
        for (std::size_t j = 0; j < keep.size(); ++j) {
            const std::size_t k = keep[j];
            const double dk = std::fabs(x[k] - y) / std::sqrt(std::max(u[k] * u[k] - uy2, 1e-300));
            if (dk > worstD) { worst = j; worstD = dk; }
        }
        keep.erase(keep.begin() + long(worst));
    }
}

// Mode of the mixture of N(x_k, u_k^2) with equal weights (Ciarlini, Cox, Pavese, Regoliosi, 2004; Duewer, 2004)
inline double mixture_mode(const std::vector<double>& x, const std::vector<double>& u, int points = 2001, int refine = 100) {
    auto f = [&](double t) {
        double s = 0;
        for (std::size_t k = 0; k < x.size(); ++k) { const double z = (t - x[k]) / u[k]; s += std::exp(-0.5 * z * z) / u[k]; }
        return s;
    };
    double lo = std::numeric_limits<double>::infinity(), hi = -lo;
    for (std::size_t k = 0; k < x.size(); ++k) { lo = std::min(lo, x[k] - 3 * u[k]); hi = std::max(hi, x[k] + 3 * u[k]); }
    const double step = (hi - lo) / (points - 1);
    int best_i = 0;
    double best_f = -1;
    for (int i = 0; i < points; ++i) { const double fi = f(lo + step * i); if (fi > best_f) { best_i = i; best_f = fi; } }
    double a = lo + step * std::max(best_i - 1, 0), b = lo + step * std::min(best_i + 1, points - 1);
    const double g = (std::sqrt(5.0) - 1.0) / 2.0;
    double c = b - g * (b - a), d = a + g * (b - a);
    for (int it = 0; it < refine; ++it) {
        if (f(c) >= f(d)) b = d; else a = c;
        c = b - g * (b - a); d = a + g * (b - a);
    }
    return (a + b) / 2.0;
}

// The preferential median next to common alternatives on the same data.
// Estimators that use u_k need all u_k > 0.
inline std::vector<Estimate> compare_estimators(const std::vector<double>& x, const std::vector<double>& u,
                                                int n = 0) {
    const double nan = std::numeric_limits<double>::quiet_NaN();
    const std::size_t m = x.size();
    const Result r = preferential_median(x, u, n);
    std::vector<Estimate> out{{"Preferential median (IF&PA)", r.pm, r.u_pm}};

    double mean = 0;
    for (double v : x) mean += v;
    mean /= double(m);
    double ss = 0;
    for (double v : x) ss += (v - mean) * (v - mean);
    out.push_back({"Arithmetic mean", mean, m > 1 ? std::sqrt(ss / double(m - 1)) / std::sqrt(double(m)) : nan});

    const double med = median(x);
    double u_med = 0;
    if (m >= 2) {
        std::vector<double> dev(m);
        for (std::size_t k = 0; k < m; ++k) dev[k] = std::fabs(x[k] - med);
        u_med = 1.2533 * 1.4826 * median(dev) / std::sqrt(double(m));
    }
    out.push_back({"Median", med, u_med});
    out.push_back({"Hodges-Lehmann", hodges_lehmann(x), nan});
    const auto aa = algorithm_a(x);
    out.push_back({"ISO 13528 Algorithm A", aa.first, aa.second});

    bool positive = true;
    for (double v : u) if (!(v > 0)) positive = false;
    if (positive) {
        std::vector<double> w(m);
        double sw = 0, swx = 0, sw2 = 0;
        for (std::size_t k = 0; k < m; ++k) {
            w[k] = 1.0 / (u[k] * u[k]);
            sw += w[k]; swx += w[k] * x[k]; sw2 += w[k] * w[k];
        }
        const double mean_w = swx / sw;
        out.push_back({"Weighted mean (1/u^2)", mean_w, 1.0 / std::sqrt(sw)});
        out.push_back({"Weighted median (1/u^2)", weighted_median(x, w), nan});
        if (m >= 2) {
            // random-effects consensus (DerSimonian-Laird), as in NIST Consensus Builder
            double q = 0;
            for (std::size_t k = 0; k < m; ++k) q += w[k] * (x[k] - mean_w) * (x[k] - mean_w);
            const double denom = sw - sw2 / sw;
            const double tau2 = denom > 0 ? std::max(0.0, (q - double(m - 1)) / denom) : 0.0;
            double sr = 0, srx = 0;
            for (std::size_t k = 0; k < m; ++k) {
                const double wr = 1.0 / (u[k] * u[k] + tau2);
                sr += wr; srx += wr * x[k];
            }
            out.push_back({"DerSimonian-Laird", srx / sr, 1.0 / std::sqrt(sr)});
            const auto lcs = largest_consistent_subset(x, u);
            out.push_back({"Largest consistent subset (Cox)", lcs.first, lcs.second});
        }
        out.push_back({"Mixture model mode", mixture_mode(x, u), nan});
    }
    return out;
}

}  // namespace ifpa