"""Shared statistical helpers for benchmark reports and scorecards."""

from __future__ import annotations

import math

import random

import statistics

from typing import Callable, Iterable

TAIL_MIN_SAMPLES = {"p95": 20, "p99": 100}

def validate_interval_settings(samples: int, confidence: float) -> None:

    if type(samples) is not int or samples <= 0:

        raise ValueError("bootstrap samples must be a positive integer")

    if (isinstance(confidence, bool) or not isinstance(confidence, (int, float)) or

            not math.isfinite(confidence) or not 0 < confidence < 1):

        raise ValueError("confidence must be finite and between zero and one")

def percentile_nearest(values: Iterable[float], percentile: float) -> float | None:

    ordered = sorted(float(value) for value in values)

    if not ordered:

        return None

    rank = math.ceil((percentile / 100.0) * len(ordered))

    index = min(max(rank, 1), len(ordered)) - 1

    return ordered[index]

def bootstrap_ci(

    values: list[float],

    statistic: Callable[[list[float]], float] | str = "median",

    *,

    samples: int = 2000,

    confidence: float = 0.95,

    seed: int = 20260609,

) -> tuple[float, float]:

    if not values:

        return float("nan"), float("nan")

    if len(values) == 1:

        return float(values[0]), float(values[0])

    if isinstance(statistic, str):

        stat = named_statistic(statistic)

    else:

        stat = statistic

    rng = random.Random(seed)

    n = len(values)

    boot = []

    for _ in range(samples):

        sample = [values[rng.randrange(n)] for _ in range(n)]

        boot.append(float(stat(sample)))

    boot.sort()

    alpha = 1.0 - confidence

    low_index = max(0, min(len(boot) - 1, math.floor((alpha / 2.0) * len(boot))))

    high_index = max(0, min(len(boot) - 1, math.ceil((1.0 - alpha / 2.0) * len(boot)) - 1))

    return boot[low_index], boot[high_index]

def named_statistic(name: str) -> Callable[[list[float]], float]:

    if name == "median":

        return lambda sample: float(statistics.median(sample))

    if name == "mean":

        return lambda sample: float(statistics.mean(sample))

    if name == "p90":

        return lambda sample: _percentile_or_nan(sample, 90.0)

    if name == "p95":

        return lambda sample: _percentile_or_nan(sample, 95.0)

    if name == "p99":

        return lambda sample: _percentile_or_nan(sample, 99.0)

    raise ValueError(name)

def _percentile_or_nan(sample: list[float], percentile: float) -> float:

    value = percentile_nearest(sample, percentile)

    return float("nan") if value is None else float(value)

def ratio_ci(

    numerator: list[float],

    denominator: list[float],

    *,

    statistic: str = "median",

    samples: int = 2000,

    confidence: float = 0.95,

    seed: int = 20260609,

) -> tuple[float, float]:

    if not numerator or not denominator:

        return float("nan"), float("nan")

    stat = named_statistic(statistic)

    rng = random.Random(seed)

    ratios = []

    for _ in range(samples):

        n_sample = [numerator[rng.randrange(len(numerator))] for _ in range(len(numerator))]

        d_sample = [denominator[rng.randrange(len(denominator))] for _ in range(len(denominator))]

        d_value = stat(d_sample)

        ratios.append(float("inf") if d_value == 0 else stat(n_sample) / d_value)

    ratios.sort()

    alpha = 1.0 - confidence

    low_index = max(0, min(len(ratios) - 1, math.floor((alpha / 2.0) * len(ratios))))

    high_index = max(0, min(len(ratios) - 1, math.ceil((1.0 - alpha / 2.0) * len(ratios)) - 1))

    return ratios[low_index], ratios[high_index]

def benjamini_hochberg(p_values: list[float], alpha: float = 0.05) -> list[bool]:

    indexed = sorted((p, i) for i, p in enumerate(p_values))

    cutoff_rank = 0

    m = len(indexed)

    for rank, (p, _) in enumerate(indexed, start=1):

        if p <= (rank / m) * alpha:

            cutoff_rank = rank

    rejected = [False] * m

    for _, original_index in indexed[:cutoff_rank]:

        rejected[original_index] = True

    return rejected

def required_n_for_relative_delta(noise_floor_pct: float, target_delta_pct: float, current_n: int) -> int | None:

    if target_delta_pct <= 0:

        return None

    if noise_floor_pct <= target_delta_pct:

        return max(1, current_n)

    ratio = noise_floor_pct / target_delta_pct

    return max(current_n + 1, math.ceil(max(1, current_n) * ratio * ratio))

# --- Statistical upgrades (additive; everything above is untouched) ---------

# Consistency constant: for normally distributed data, sigma_hat = 1.4826 * MAD,

# i.e. the MAD scaled by 1.4826 is an unbiased estimator of the standard

# deviation. Its reciprocal is the 0.6745 that appears in the classic

# Iglewicz-Hoaglin modified z-score (0.6745 = 1 / 1.4826).

MAD_CONSISTENCY = 1.4826

def mad_outliers(values: list[float], *, threshold: float = 14.826) -> dict:

    """Flag outliers via MAD-based modified z-scores. Never deletes anything.

    Algebra (so the threshold semantics are explicit):

        med  = median(values)

        MAD  = median(|x_i - med|)

        z_i  = (x_i - med) / (1.4826 * MAD)        # modified z-score

             = 0.6745 * (x_i - med) / MAD          # identical: 0.6745 = 1/1.4826

    Because ``1.4826 * MAD`` estimates one standard deviation on normal data,

    ``z_i`` is a robust z-score in sigma units. The hyperfine-style default

    threshold ``14.826`` therefore means "more than ~14.8 robust standard

    deviations from the median", equivalently a threshold of

    ``14.826 * 1.4826 ~= 21.98`` on the raw ratio ``|x_i - med| / MAD``. That

    is deliberately far more conservative than the conventional Iglewicz-

    Hoaglin cut of 3.5: only blatant escapees (stalls, paging storms) are

    flagged, not ordinary timing jitter.

    Degenerate scale: when ``MAD == 0`` (constant data, or more than half the

    samples tied at the median) the estimator is blind — every off-median

    point would divide by zero and be flagged at infinity. We document the

    opposite policy: MAD == 0 flags NOTHING and reports all modified z-scores

    as 0.0. A scale-collapsed sample is a job for eyeballs, not a formula.

    First-run special case: if index 0 is the SOLE outlier and it is on the

    slow side (``values[0] > median``), that is the classic warmup / cold-

    cache signature, reported via ``first_run_warmup_hint``. It is a hint to

    add warmups or run a cold-cache lane, not a reason to discard the sample.

    Returns ``{"outlier_indices", "modified_z", "first_run_warmup_hint",

    "policy": "flag_never_delete"}``. Flagged rows are metadata for humans

    and gates; no caller may delete them (see docs/runner-discipline.md).

    """

    vals = [float(v) for v in values]

    result = {

        "outlier_indices": [],

        "modified_z": [0.0] * len(vals),

        "first_run_warmup_hint": False,

        "policy": "flag_never_delete",

    }

    if not vals:

        return result

    med = float(statistics.median(vals))

    mad = float(statistics.median(abs(v - med) for v in vals))

    if mad == 0:

        return result  # scale collapsed: flag nothing (documented above)

    z_scores = [(v - med) / (MAD_CONSISTENCY * mad) for v in vals]

    outlier_indices = [i for i, z in enumerate(z_scores) if abs(z) > threshold]

    result["modified_z"] = z_scores

    result["outlier_indices"] = outlier_indices

    result["first_run_warmup_hint"] = outlier_indices == [0] and vals[0] > med

    return result

def adaptive_rep_plan(

    samples: list[float],

    *,

    min_runs: int = 10,

    time_budget_s: float,

    per_run_s_estimate: float,

    stddev_pct_threshold: float = 5.0,

    sigma_escapee: float = 3.0,

    cap_multiplier: int = 10,

) -> dict:

    """Pure planning math for adaptive repetitions: floor + top-up + cap.

    Policy (hyperfine floor, PTS-style variance top-up):

    - Floor: every cell gets at least ``min_runs`` repetitions.

    - Budget: ``time_budget_s / per_run_s_estimate`` is the total run count

      the wall-time budget pays for. The cap on total runs is

      ``min(min_runs * cap_multiplier, budget_runs)`` — long trials shrink

      the cap, short trials let the variance loop breathe.

    - Top-up triggers (checked only once the floor is met), in order:

        1. ``stddev_above_threshold``: sample %stddev (stdev/|mean| * 100)

           is at or above ``stddev_pct_threshold``;

        2. ``sigma_escapee``: any sample sits more than ``sigma_escapee``

           sample standard deviations from the mean.

      Either trigger grants one more floor-sized batch (``min_runs`` extra

      runs), clipped to the cap. Callers loop: measure, re-plan, repeat.

    - ``budget_exhausted``: a trigger (or the floor itself) wants more runs

      but the cap grants none. ``satisfied`` is False — the cell is reported

      with the variance it has, never silently blessed.

    - ``satisfied``: floor met and no trigger fires; ``additional_runs`` is 0.

    Returns ``{"satisfied", "additional_runs", "reason", "cap"}`` where

    ``reason`` is one of ``below_min_runs | stddev_above_threshold |

    sigma_escapee | satisfied | budget_exhausted``.

    """

    if per_run_s_estimate <= 0:

        raise ValueError("per_run_s_estimate must be positive")

    n = len(samples)

    budget_runs = int(time_budget_s // per_run_s_estimate)

    cap = max(0, min(min_runs * cap_multiplier, budget_runs))

    def plan(additional: int, reason: str) -> dict:

        return {

            "satisfied": reason == "satisfied",

            "additional_runs": additional,

            "reason": reason,

            "cap": cap,

        }

    if n < min_runs:

        grant = min(min_runs, cap) - n

        if grant <= 0:

            return plan(0, "budget_exhausted")

        return plan(grant, "below_min_runs")

    vals = [float(s) for s in samples]

    mean = statistics.fmean(vals)

    stdev = statistics.stdev(vals) if n >= 2 else 0.0

    if mean != 0:

        stddev_pct = abs(stdev / mean) * 100.0

    else:

        stddev_pct = 0.0 if stdev == 0 else float("inf")

    trigger = None

    if stddev_pct >= stddev_pct_threshold:

        trigger = "stddev_above_threshold"

    elif stdev > 0 and any(abs(v - mean) > sigma_escapee * stdev for v in vals):

        trigger = "sigma_escapee"

    if trigger is None:

        return plan(0, "satisfied")

    grant = min(n + min_runs, cap) - n

    if grant <= 0:

        return plan(0, "budget_exhausted")

    return plan(grant, trigger)

def benjamini_hochberg_detail(p_values: list[float], *, fdr: float = 0.05) -> dict:

    """Benjamini-Hochberg step-up with full detail, for full-matrix claims.

    Additive companion to the legacy ``benjamini_hochberg`` above (which

    returns a bare ``list[bool]`` and is imported elsewhere, so its signature

    is frozen). New code should prefer this one.

    Standard step-up procedure: sort the m p-values ascending (a stable sort,

    so ties keep input order); find the largest rank k (1-based) such that

    ``p_(k) <= (k / m) * fdr``; reject every hypothesis whose p-value is

    ``<= p_(k)``. Rejecting by value rather than by sorted position makes the

    result stable for ties: p-values equal to the critical one are all

    rejected, regardless of how the sort ordered them.

    Returns ``{"rejected": [bool per input, input order], "critical_index":

    k-1 as a 0-based rank in the sorted order (None if nothing rejects),

    "threshold": the critical p-value p_(k) (None if nothing rejects)}``.

    Empty input rejects nothing.

    """

    m = len(p_values)

    if m == 0:

        return {"rejected": [], "critical_index": None, "threshold": None}

    order = sorted(range(m), key=lambda i: float(p_values[i]))

    critical_index = None

    for rank, original in enumerate(order, start=1):

        if float(p_values[original]) <= (rank / m) * fdr:

            critical_index = rank - 1

    if critical_index is None:

        return {"rejected": [False] * m, "critical_index": None, "threshold": None}

    threshold = float(p_values[order[critical_index]])

    rejected = [float(p) <= threshold for p in p_values]

    return {"rejected": rejected, "critical_index": critical_index, "threshold": threshold}

def two_stage_confirmation(first_pct: float, second_pct: float, *, threshold_pct: float) -> dict:

    """Two-stage confirmation for borderline regression verdicts (pure math).

    Stage one is the normal devloop comparison; stage two is an independent

    re-measurement. A delta is ``confirmed`` only if BOTH stages exceed the

    threshold — a single borderline excursion is never a verdict. The first

    stage is ``borderline`` when it lands within +/-25% of the threshold

    (e.g. threshold 10% -> the 7.5%..12.5% band), which is the band where a

    second stage is worth running. Wiring into devloop is future work; this

    helper only encodes the decision rule.

    """

    first_exceeds = first_pct > threshold_pct

    second_exceeds = second_pct > threshold_pct

    borderline = abs(first_pct - threshold_pct) <= 0.25 * abs(threshold_pct)

    return {

        "confirmed": first_exceeds and second_exceeds,

        "borderline": borderline,

        "first_exceeds": first_exceeds,

        "second_exceeds": second_exceeds,

        "threshold_pct": float(threshold_pct),

    }