# =============================================================================
#  Project   : anomaly-atlas
#  File      : src/anomaly_atlas/stats/multiple_testing.py
#  Purpose   : FDR/Bonferroni corrections (White RC / SPA / DSR arrive in expF)
#  Author    : Simon-Pierre Boucher
#  Contact   : contact@spboucher.ai
#  Data src  : hfmarketdata.io (sole data source)
#  Created   : 2026-08-12
#  Modified  : 2026-08-12
#  Platform  : macOS / Apple Silicon (arm64)
#  License   : All rights reserved (research code)
# =============================================================================
"""Multiple-testing corrections for the scan stage (Level-0 triage).

Scans use Benjamini–Hochberg FDR (the right error rate when a controlled
fraction of false leads into the next stage is acceptable). Level-1
promotion uses SPA/StepM against artifact nulls — implemented with expF.
Validated on synthetic ground truth (§8.1 gate).
"""

from __future__ import annotations

import numpy as np


def benjamini_hochberg(pvals: np.ndarray, alpha: float = 0.05) -> np.ndarray:
    """BH step-up FDR procedure. Returns a boolean rejection mask.

    NaN p-values are never rejected but still COUNT toward m (conservative:
    an unevaluable test is a spent test, not a free one).
    """
    p = np.asarray(pvals, dtype=float)
    m = len(p)
    if m == 0:
        return np.zeros(0, dtype=bool)
    finite = np.where(np.isfinite(p))[0]
    reject = np.zeros(m, dtype=bool)
    if len(finite) == 0:
        return reject
    order = finite[np.argsort(p[finite])]
    thresholds = alpha * (np.arange(1, len(order) + 1) / m)
    passed = np.where(p[order] <= thresholds)[0]
    if len(passed):
        reject[order[: passed.max() + 1]] = True
    return reject


def bonferroni(pvals: np.ndarray, alpha: float = 0.05) -> np.ndarray:
    """Bonferroni FWE mask (reported alongside FDR for reference)."""
    p = np.asarray(pvals, dtype=float)
    return np.isfinite(p) & (p <= alpha / max(len(p), 1))


def bootstrap_pvalue(samples: np.ndarray, null_value: float = 0.0) -> float:
    """Two-sided percentile-bootstrap p-value of a statistic vs a null value.

    p = 2 * min(P(boot <= null), P(boot >= null)), with the +1/(B+1)
    correction so p is never exactly 0. Triage-grade inference for scans —
    Level-1 promotion re-tests with SPA machinery.
    """
    s = np.asarray(samples, dtype=float)
    s = s[np.isfinite(s)]
    b = len(s)
    if b < 50:
        return float("nan")
    lo = (np.sum(s <= null_value) + 1) / (b + 1)
    hi = (np.sum(s >= null_value) + 1) / (b + 1)
    return float(min(1.0, 2.0 * min(lo, hi)))