# IF&PA — interval fusion with preference aggregation.
# Preferential median (PM) of a set of intervals [x_k - u_k, x_k + u_k].
#
# Base R, no packages.
#
# 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:
#   source("ifpa.R")
#   r <- preferential_median(c(4.375, 4.9), c(0.175, 0.2), n = 3)
#   r$pm   # 4.65

preferential_median <- function(x, u, n = NULL, n_min = 11, n_max = 201) {
  x <- as.numeric(x)
  u <- as.numeric(u)
  if (length(u) == 1) u <- rep(u, length(x))
  if (length(x) != length(u)) stop("x and u have different lengths")
  if (length(x) == 0) stop("empty set of intervals")
  if (any(u < 0)) stop("uncertainties must be non-negative")
  if (!all(is.finite(c(x, u)))) stop("NaN or infinity in data")

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

  if (an - a1 <= 0) {
    return(list(pm = a1, u_pm = 0, u_grid = 0, u_sampling = 0, n = 1, h = 0,
                winners = 1, converged = TRUE, grid = a1, coverage = length(lo)))
  }

  if (is.null(n)) {
    width_min <- min(hi - lo)                      # the narrowest interval, 2 min u_k
    n <- if (width_min > 0) ceiling((an - a1) / width_min) else n_min
    n <- max(n, n_min)
  }
  n <- max(as.integer(n), 2L)

  # median() in R averages the two middle values, as required
  sampling <- function(c) {
    m <- length(c)
    if (m < 2) return(0)
    1.2533 * 1.4826 * median(abs(c - median(c))) / sqrt(m)
  }

  repeat {
    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
    eps <- max(1e-9 * h, 8 * .Machine$double.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}: O(n log m) instead of O(n m)
    coverage <- findInterval(grid, sort(lo - eps)) - findInterval(grid, sort(hi + eps), left.open = TRUE)
    idx <- which(coverage == max(coverage))        # 1-based indices
    adjacent <- length(idx) <= 1 || all(diff(idx) == 1)
    if (adjacent || n >= n_max) {
      if (!adjacent) {
        # widest block of consecutive winners; tie -> median closest to centres
        blocks <- split(idx, cumsum(c(1, diff(idx) != 1)))
        centre <- median(centres)
        best <- NULL
        for (b in blocks) {
          d <- abs(median(grid[b]) - centre)
          if (is.null(best) || length(b) > length(best$b) ||
              (length(b) == length(best$b) && d < best$d)) best <- list(b = b, d = d)
        }
        idx <- best$b
      }
      u_grid <- 0.41 * h
      u_samp <- sampling(centres)
      return(list(pm = median(grid[idx]), u_pm = sqrt(u_grid^2 + u_samp^2),
                  u_grid = u_grid, u_sampling = u_samp, n = n, h = h,
                  winners = idx, converged = adjacent,
                  grid = grid, coverage = coverage))
    }
    n <- n + 1L
  }
}

# ---------------------------------------------------------------------------
# Alternatives for comparison: common practice for combining x_k +/- u_k
# ---------------------------------------------------------------------------

# 0.95 quantile of chi-square, Wilson-Hilferty approximation (error < 1 %);
# kept instead of qchisq() so that all implementations give the same result
chi2_quantile_95 <- function(df) {
  a <- 2 / (9 * df)
  df * (1 - a + 1.6448536269514722 * sqrt(a))^3
}

# Weighted median: smallest x with cumulative weight >= half; exactly half -> midpoint
weighted_median <- function(x, w) {
  o <- order(x)
  x <- x[o]
  w <- w[o]
  total <- sum(w)
  cum <- cumsum(w)
  j <- which(cum >= total / 2 | abs(cum - total / 2) <= 1e-12 * total)[1]
  if (abs(cum[j] - total / 2) <= 1e-12 * total && j < length(x)) return((x[j] + x[j + 1]) / 2)
  x[j]
}

# Hodges-Lehmann estimate: median of all pairwise (Walsh) averages, i <= j
hodges_lehmann <- function(x) {
  s <- outer(x, x, "+") / 2
  median(s[upper.tri(s, diag = TRUE)])
}

# Robust mean by ISO 13528 Algorithm A; u = 1.25 s* / sqrt(m)
algorithm_a <- function(x, tol = 1e-10, max_iter = 1000) {
  m <- length(x)
  xs <- median(x)
  ss <- 1.483 * median(abs(x - xs))
  if (ss == 0 || m < 2) return(list(value = xs, s = ss, u = 0))
  for (it in seq_len(max_iter)) {
    d <- 1.5 * ss
    cl <- pmin(pmax(x, xs - d), xs + d)
    nx <- mean(cl)
    ns <- 1.134 * sqrt(sum((cl - nx)^2) / (m - 1))
    done <- abs(nx - xs) <= tol * max(abs(xs), 1e-300) && abs(ns - ss) <= tol * ss
    xs <- nx
    ss <- ns
    if (done || ss == 0) break
  }
  list(value = xs, s = ss, u = 1.25 * ss / sqrt(m))
}

# Random-effects consensus (DerSimonian-Laird), as in NIST Consensus Builder
dersimonian_laird <- function(x, u) {
  m <- length(x)
  w <- 1 / u^2
  sw <- sum(w)
  mean_w <- sum(w * x) / sw
  q <- sum(w * (x - mean_w)^2)
  denom <- sw - sum(w^2) / sw
  tau2 <- if (denom > 0) max(0, (q - (m - 1)) / denom) else 0
  ws <- 1 / (u^2 + tau2)
  list(value = sum(ws * x) / sum(ws), u = 1 / sqrt(sum(ws)), tau = sqrt(tau2))
}

# Largest consistent subset (after Cox, 2007), sequential exclusion by chi-square at 0.05
largest_consistent_subset <- function(x, u) {
  keep <- seq_along(x)
  repeat {
    w <- 1 / u[keep]^2
    sw <- sum(w)
    y <- sum(w * x[keep]) / sw
    chi2 <- sum(w * (x[keep] - y)^2)
    if (length(keep) <= 2 || chi2 <= chi2_quantile_95(length(keep) - 1)) {
      return(list(value = y, u = 1 / sqrt(sw), excluded = length(x) - length(keep)))
    }
    d <- abs(x[keep] - y) / sqrt(pmax(u[keep]^2 - 1 / sw, 1e-300))
    keep <- keep[-which.max(d)]
  }
}

# Mode of the mixture of N(x_k, u_k^2) with equal weights (Ciarlini, Cox, Pavese, Regoliosi, 2004; Duewer, 2004)
mixture_mode <- function(x, u, points = 2001, refine = 100) {
  f <- function(t) sum(exp(-0.5 * ((t - x) / u)^2) / u)
  lo <- min(x - 3 * u)
  hi <- max(x + 3 * u)
  step <- (hi - lo) / (points - 1)
  fv <- vapply(lo + step * (0:(points - 1)), f, numeric(1))
  i <- which.max(fv) - 1                       # 0-based, first maximum
  a <- lo + step * max(i - 1, 0)
  b <- lo + step * min(i + 1, points - 1)
  g <- (sqrt(5) - 1) / 2
  cc <- b - g * (b - a)
  dd <- a + g * (b - a)
  for (it in seq_len(refine)) {
    if (f(cc) >= f(dd)) b <- dd else a <- cc
    cc <- b - g * (b - a)
    dd <- a + g * (b - a)
  }
  (a + b) / 2
}

# The preferential median next to common alternatives on the same data.
# Returns a data.frame (name, value, u); u is NA where no standard formula exists.
# Estimators that use u_k need all u_k > 0.
compare_estimators <- function(x, u, n = NULL) {
  x <- as.numeric(x)
  u <- as.numeric(u)
  if (length(u) == 1) u <- rep(u, length(x))
  m <- length(x)
  r <- preferential_median(x, u, n = n)
  med <- median(x)
  u_med <- if (m < 2) 0 else 1.2533 * 1.4826 * median(abs(x - med)) / sqrt(m)
  aa <- algorithm_a(x)
  out <- data.frame(
    name = c("Preferential median (IF&PA)", "Arithmetic mean", "Median", "Hodges-Lehmann",
             "ISO 13528 Algorithm A"),
    value = c(r$pm, mean(x), med, hodges_lehmann(x), aa$value),
    u = c(r$u_pm, if (m > 1) sd(x) / sqrt(m) else NA, u_med, NA, aa$u),
    stringsAsFactors = FALSE)
  add <- function(df, name, value, u) rbind(df, data.frame(name = name, value = value, u = u))
  if (all(u > 0)) {
    w <- 1 / u^2
    out <- add(out, "Weighted mean (1/u^2)", sum(w * x) / sum(w), 1 / sqrt(sum(w)))
    out <- add(out, "Weighted median (1/u^2)", weighted_median(x, w), NA)
    if (m >= 2) {
      dl <- dersimonian_laird(x, u)
      out <- add(out, "DerSimonian-Laird", dl$value, dl$u)
      lcs <- largest_consistent_subset(x, u)
      out <- add(out, "Largest consistent subset (Cox)", lcs$value, lcs$u)
    }
    out <- add(out, "Mixture model mode", mixture_mode(x, u), NA)
  }
  out
}