hegeltest 0.27.0

#!/usr/bin/env python3
"""Generate scipy-derived oracle tables for the vendored statistics module.

Run from the repo root in an environment that has scipy installed:

    python scripts/generate_statistics_oracle.py

Emits `tests/embedded/native/statistics_oracle_data.rs` containing arrays of
`(df, t, expected)` and `(df, p, expected)` rows, sourced from
`scipy.special.stdtr` and `scipy.special.stdtrit`. The Rust test compares the
vendored `stdtr` / `stdtrit` against these values.
"""

from __future__ import annotations

import math
import random
from pathlib import Path

import scipy.special

REPO_ROOT = Path(__file__).resolve().parent.parent
OUT_PATH = REPO_ROOT / "tests" / "embedded" / "native" / "statistics_oracle_data.rs"


def stdtr_grid() -> list[tuple[int, float, float]]:
    """Mirror of Hypothesis's `test_stdtr_explicit` plus a random sample.

    Hypothesis's explicit grid is `df in {1, 2, 3, 4, 10, 100}` times
    `t in {-100, -1, -0.5, 0, 0.5, 1, 100}`. We add a randomized fan-out so
    the table also exercises the parity-of-df / large-|t| code paths.
    """
    rng = random.Random(0)
    rows: list[tuple[int, float, float]] = []
    # Explicit grid (mirrors test_stdtr_explicit).
    explicit_df = [1, 2, 3, 4, 10, 100]
    explicit_t = [-100.0, -1.0, -0.5, 0.0, 0.5, 1.0, 100.0]
    for df in explicit_df:
        for t in explicit_t:
            rows.append((df, t, float(scipy.special.stdtr(df, t))))
    # Randomized fan-out across all the branches of stdtr (odd df with df>1,
    # even df, very large |t|, tiny |t|).
    for _ in range(200):
        df = rng.randint(1, 100)
        # mix of small / large magnitudes
        exp = rng.uniform(-3, 6)
        t = math.copysign(10**exp, rng.random() - 0.5)
        rows.append((df, t, float(scipy.special.stdtr(df, t))))
    return rows


def stdtrit_strict_grid() -> list[tuple[int, float, float]]:
    """`df ∈ {1, 2}` rows — closed-form, expected to match scipy to ~last ULP."""
    rng = random.Random(1)
    rows: list[tuple[int, float, float]] = []
    for df in (1, 2):
        # cover the symmetric p=0.5 hinge and both tails
        for p in (0.0001, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.9999):
            rows.append((df, p, float(scipy.special.stdtrit(df, p))))
        # randomized fan-out into the tails
        for _ in range(50):
            # log-uniform p in (1e-12, 1 - 1e-12) to exercise extreme tails
            log_p = rng.uniform(-12, -1e-6)
            p = 10**log_p
            if rng.random() < 0.5:
                p = 1 - p
            rows.append((df, p, float(scipy.special.stdtrit(df, p))))
    return rows


def stdtrit_lax_grid() -> list[tuple[int, float, float]]:
    """`df >= 3` rows. Tolerance is lax (Newton vs Halley)."""
    rng = random.Random(2)
    rows: list[tuple[int, float, float]] = []
    explicit_df = [3, 4, 5, 10, 100]
    for df in explicit_df:
        for p in (0.001, 0.01, 0.1, 0.5, 0.9, 0.99, 0.999):
            rows.append((df, p, float(scipy.special.stdtrit(df, p))))
    for _ in range(100):
        df = rng.randint(3, 100)
        # keep p well away from 0/1 because Newton's accuracy degrades there
        p = rng.uniform(1e-6, 1 - 1e-6)
        rows.append((df, p, float(scipy.special.stdtrit(df, p))))
    return rows


def fmt_float(x: float) -> str:
    # Emit a literal Rust f64 in shortest-round-trip form so a re-parse via
    # f64::from_str recovers the exact same value.
    if math.isnan(x):
        return "f64::NAN"
    if math.isinf(x):
        return "f64::INFINITY" if x > 0 else "f64::NEG_INFINITY"
    return repr(x)


def main() -> None:
    stdtr_rows = stdtr_grid()
    stdtrit_strict_rows = stdtrit_strict_grid()
    stdtrit_lax_rows = stdtrit_lax_grid()

    out_lines: list[str] = []
    out_lines.append(
        "// Auto-generated by scripts/generate_statistics_oracle.py.\n"
        "// Do not edit by hand; regenerate against scipy if the oracle\n"
        "// needs refreshing. Each row is (df, x, expected_scipy_result).\n"
        f"// Source: scipy {scipy.__version__}.\n"
    )
    out_lines.append(
        f"pub const STDTR_CASES: &[(i32, f64, f64)] = &[\n"
    )
    for df, t, expected in stdtr_rows:
        out_lines.append(f"    ({df}, {fmt_float(t)}, {fmt_float(expected)}),\n")
    out_lines.append("];\n\n")

    out_lines.append(
        f"pub const STDTRIT_STRICT_CASES: &[(i32, f64, f64)] = &[\n"
    )
    for df, p, expected in stdtrit_strict_rows:
        out_lines.append(f"    ({df}, {fmt_float(p)}, {fmt_float(expected)}),\n")
    out_lines.append("];\n\n")

    out_lines.append(
        f"pub const STDTRIT_LAX_CASES: &[(i32, f64, f64)] = &[\n"
    )
    for df, p, expected in stdtrit_lax_rows:
        out_lines.append(f"    ({df}, {fmt_float(p)}, {fmt_float(expected)}),\n")
    out_lines.append("];\n")

    OUT_PATH.write_text("".join(out_lines))
    print(f"wrote {OUT_PATH} ({len(stdtr_rows)} stdtr, "
          f"{len(stdtrit_strict_rows)} stdtrit-strict, "
          f"{len(stdtrit_lax_rows)} stdtrit-lax rows)")


if __name__ == "__main__":
    main()