#  Copyright (c) 2026. Jose M. Requena-Plens
"""Tests for ISO/PAS 20065:2016 tonal audibility (engineering method).

Anchored on the Annex E combustion-engine worked example: the critical band
(Formulae (2)-(5)), the critical-band level LG (Formula (12)), the masking
index av (Formula (13)), the per-tone audibility ΔL = LT − LG − av
(Formula (14)) and the energy-mean mean audibility (Formula (20)).
"""

from __future__ import annotations

import math

import matplotlib

matplotlib.use("Agg")
import numpy as np
import pytest

import reference_data as ref
from phonometry import (
    HANNING_BANDWIDTH_FACTOR,
    NO_TONE_AUDIBILITY,
    ToneAudibilityResult,
    analyze_spectrum,
    assess_tones,
    audibility_from_levels,
    combined_tone_level,
    critical_band_corners,
    critical_band_level,
    critical_bandwidth_engineering,
    energy_sum_level,
    masking_index,
    mean_audibility,
    mean_narrowband_level,
    resolve_tones_separately,
    tone_audibility,
    tone_level,
    two_tone_separation_frequency,
)


# ---------------------------------------------------------------------------
# Critical band about the tone (Formulae (2)-(5))
# ---------------------------------------------------------------------------
def test_critical_bandwidth_closed_form() -> None:
    ft = 137.3
    expected = 25.0 + 75.0 * (1.0 + 1.4 * (ft / 1000.0) ** 2) ** 0.69
    assert critical_bandwidth_engineering(ft) == pytest.approx(expected)


def test_critical_bandwidth_grows_with_frequency() -> None:
    assert critical_bandwidth_engineering(2000.0) > critical_bandwidth_engineering(
        100.0
    )


def test_critical_band_corners_geometry() -> None:
    # f1·f2 = fT² and f2 − f1 = Δfc (geometric placement, Formulae (3)-(5)).
    for ft in (55.0, 137.3, 592.2, 1582.7):
        f1, f2 = critical_band_corners(ft)
        assert math.sqrt(f1 * f2) == pytest.approx(ft, rel=1e-12)
        assert f2 - f1 == pytest.approx(critical_bandwidth_engineering(ft), rel=1e-12)
        assert 0.0 < f1 < ft < f2


# ---------------------------------------------------------------------------
# Masking index (Formula (13)) and critical-band level (Formula (12))
# ---------------------------------------------------------------------------
def test_masking_index_closed_form() -> None:
    f = 592.2
    expected = -2.0 - math.log10(1.0 + (f / 502.0) ** 2.5)
    assert masking_index(f) == pytest.approx(expected)


def test_masking_index_matches_annex_e() -> None:
    assert masking_index(137.3) == pytest.approx(ref.ISO20065_AV_137, abs=0.005)
    assert masking_index(592.2) == pytest.approx(ref.ISO20065_AV_592, abs=0.005)


def test_masking_index_decreases_with_frequency() -> None:
    assert masking_index(2000.0) < masking_index(100.0) < 0.0


def test_critical_band_level_matches_annex_e() -> None:
    value = critical_band_level(49.22, 137.3, ref.ISO20065_LINE_SPACING)
    assert value == pytest.approx(ref.ISO20065_LG_137, abs=0.05)


def test_critical_band_level_bandwidth_ratio() -> None:
    ls, ft, df = 50.0, 300.0, 2.7
    dfc = critical_bandwidth_engineering(ft)
    assert critical_band_level(ls, ft, df) == pytest.approx(
        ls + 10.0 * math.log10(dfc / df)
    )


# ---------------------------------------------------------------------------
# Audibility (Formula (14))
# ---------------------------------------------------------------------------
def test_audibility_from_levels_identity() -> None:
    assert audibility_from_levels(67.96, 64.98, -2.02) == pytest.approx(5.0)


def test_tone_audibility_matches_annex_e() -> None:
    for ft, ls, lt, expected in ref.ISO20065_ANNEX_E_TONES:
        value = tone_audibility(lt, ls, ft, ref.ISO20065_LINE_SPACING)
        assert value == pytest.approx(expected, abs=0.03)


# ---------------------------------------------------------------------------
# Energy-sum level with the window correction (Formulae (6)/(8))
# ---------------------------------------------------------------------------
def test_energy_sum_hanning_correction() -> None:
    # Two equal 80 dB lines: 10 lg(2·10^8) = 83.01 dB, minus 1.76 dB Hanning.
    assert energy_sum_level([80.0, 80.0]) == pytest.approx(
        83.0103 - 10.0 * math.log10(1.5), abs=1e-4
    )


def test_energy_sum_rectangular_is_plain_sum() -> None:
    value = energy_sum_level([80.0, 80.0], effective_bandwidth_factor=1.0)
    assert value == pytest.approx(10.0 * math.log10(2.0e8), abs=1e-6)


def test_hanning_factor_constant() -> None:
    assert HANNING_BANDWIDTH_FACTOR == 1.5


# ---------------------------------------------------------------------------
# From-spectrum front-end: LS (Formula (6)) and LT (Formula (8))
# ---------------------------------------------------------------------------
def test_mean_narrowband_level_matches_annex_e() -> None:
    ls = mean_narrowband_level(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3
    )
    assert ls == pytest.approx(ref.ISO20065_E1_LS, abs=0.02)


def test_tone_level_matches_annex_e() -> None:
    ls = mean_narrowband_level(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3
    )
    lt = tone_level(ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3, ls)
    assert lt == pytest.approx(ref.ISO20065_E1_LT, abs=0.02)


def test_from_spectrum_full_chain_matches_annex_e() -> None:
    # Spectrum -> LS -> LT -> ΔL reproduces the 137.3 Hz tone (ΔL ≈ 4.99 dB).
    ls = mean_narrowband_level(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3
    )
    lt = tone_level(ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3, ls)
    assert tone_audibility(lt, ls, 137.3, 2.7) == pytest.approx(4.99, abs=0.03)


def test_mean_narrowband_level_rectangular_higher() -> None:
    # Dropping the Hanning correction (factor 1.0) raises LS by at least the
    # 1.76 dB correction; the iterative exclusion also then keeps more lines,
    # so the gap is not a constant offset.
    ls_h = mean_narrowband_level(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3
    )
    ls_r = mean_narrowband_level(
        ref.ISO20065_E1_LEVELS,
        ref.ISO20065_E1_FREQUENCIES,
        137.3,
        effective_bandwidth_factor=1.0,
    )
    assert ls_r > ls_h + 10.0 * math.log10(1.5) - 1e-6


def test_tone_level_single_line_when_isolated() -> None:
    # A lone peak far above its neighbours yields just itself (minus Hanning).
    freqs = [100.0, 102.7, 105.4, 108.1, 110.8]
    levels = [40.0, 41.0, 70.0, 41.0, 40.0]
    ls = 45.0
    lt = tone_level(levels, freqs, 105.4, ls)
    assert lt == pytest.approx(70.0 - 10.0 * math.log10(1.5), abs=1e-6)


def test_mean_narrowband_level_empty_critical_band_raises() -> None:
    # Lines far below the critical band of a 1 kHz tone (≈922-1085 Hz).
    with pytest.raises(ValueError, match="critical band"):
        mean_narrowband_level([50.0, 51.0, 52.0], [100.0, 110.0, 120.0], 1000.0)


def test_mean_narrowband_level_few_lines_returns_unpruned_mean() -> None:
    # Fewer than five lines on one side of the tone: the iteration stops on the
    # first pass and returns the pre-prune mean (Annex D fallback).
    freqs = [930.0, 950.0, 970.0, 990.0, 1010.0, 1030.0, 1050.0, 1070.0]
    levels = [50.0, 50.0, 50.0, 50.0, 60.0, 50.0, 50.0, 50.0]
    ls = mean_narrowband_level(levels, freqs, 1000.0)
    # Line under investigation = nearest to 1000 Hz (990 Hz); only three lines
    # sit below it, so no pruning happens and LS is the raw mean of the rest.
    band = [50.0, 50.0, 50.0, 60.0, 50.0, 50.0, 50.0]
    expected = 10.0 * math.log10(
        np.mean([10.0 ** (x / 10.0) for x in band])
    ) - 10.0 * math.log10(1.5)
    assert ls == pytest.approx(expected)


def test_mean_narrowband_level_rejects_length_mismatch() -> None:
    with pytest.raises(ValueError, match="share their length"):
        mean_narrowband_level([50.0, 51.0], [100.0], 100.0)


def test_mean_narrowband_level_rejects_unsorted_frequencies() -> None:
    with pytest.raises(ValueError, match="increasing"):
        mean_narrowband_level([50.0, 51.0, 52.0], [100.0, 99.0, 101.0], 100.0)


def test_tone_level_rejects_non_finite_ls() -> None:
    with pytest.raises(ValueError):
        tone_level([50.0, 60.0, 50.0], [100.0, 102.7, 105.4], 102.7, float("nan"))


# ---------------------------------------------------------------------------
# Full-spectrum detection (Clause 5.3.8) and FG combination (Formula (17))
# ---------------------------------------------------------------------------
def test_analyze_spectrum_detects_annex_e_tones() -> None:
    result = analyze_spectrum(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, ref.ISO20065_LINE_SPACING
    )
    found = sorted(round(float(f), 1) for f in result.tone_frequencies)
    assert found == sorted(ref.ISO20065_E1_TONE_FREQUENCIES)


def test_analyze_spectrum_137hz_tone_audibility() -> None:
    # The 137.3 Hz tone's critical band is fully inside Table E.1, so its
    # detected audibility matches the oracle (ΔL ≈ 4.99 dB).
    result = analyze_spectrum(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, ref.ISO20065_LINE_SPACING
    )
    i = int(np.argmin(np.abs(result.tone_frequencies - 137.3)))
    assert result.audibilities[i] == pytest.approx(4.99, abs=0.03)


def test_detect_tone_at_first_line_not_dropped() -> None:
    # Regression: a peak on the very first line must not be skipped because the
    # left-neighbour test wrapped around to the last line (lev[-1]).
    from phonometry.tone_audibility import _detect_tones

    freqs = np.arange(100.0, 100.0 + 2.7 * 40, 2.7)
    levels = np.full(freqs.size, 50.0)
    levels[0] = 90.0     # strong peak on the first line
    levels[-1] = 95.0    # even stronger last line (would trigger the wrap bug)
    detected = _detect_tones(levels, freqs, 2.7, 1.5)
    assert 0 in [peak for peak, _lo, _hi, _ls in detected]


def test_analyze_spectrum_rejects_inaudible_peak() -> None:
    # A single weak peak (below LS+6) yields no audible tone.
    freqs = list(np.arange(900.0, 1100.0, 2.7))
    levels = [50.0] * len(freqs)
    levels[len(freqs) // 2] = 52.0
    with pytest.raises(ValueError, match="No audible tone"):
        analyze_spectrum(levels, freqs, 2.7)


def test_combined_tone_level_matches_annex_e() -> None:
    value = combined_tone_level(
        ref.ISO20065_E1_LEVELS,
        ref.ISO20065_E1_FREQUENCIES,
        ref.ISO20065_E1_TONE_FREQUENCIES,
        ref.ISO20065_E1_TONE_LS,
    )
    assert value == pytest.approx(ref.ISO20065_E1_LT_FG, abs=0.02)


def test_combined_tone_level_single_tone_equals_tone_level() -> None:
    # With one tone the FG combination reduces to that tone's level.
    fg = combined_tone_level(
        ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, [137.3], [49.22]
    )
    lt = tone_level(ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, 137.3, 49.22)
    assert fg == pytest.approx(lt)


def test_combined_tone_level_no_double_counting() -> None:
    # Two tones whose line runs overlap must not sum a shared line twice, so the
    # combined level cannot exceed the plain energy sum of the union of lines.
    fg = combined_tone_level(
        ref.ISO20065_E1_LEVELS,
        ref.ISO20065_E1_FREQUENCIES,
        [134.6, 137.3],  # adjacent tones sharing lines
        [49.22, 49.22],
    )
    assert np.isfinite(fg)


def test_combined_tone_level_rejects_length_mismatch() -> None:
    with pytest.raises(ValueError, match="share their length"):
        combined_tone_level(
            ref.ISO20065_E1_LEVELS, ref.ISO20065_E1_FREQUENCIES, [137.3], [49.0, 50.0]
        )


# ---------------------------------------------------------------------------
# Separate evaluation of two tones below 1000 Hz (Formulae (18)/(19))
#
# No numeric worked example exercises this branch (the Annex E band groups three
# tones, so the "exactly two tones" rule never fires there). These tests fix the
# closed-form Formula (19) and the decision logic against the DIN 45681 Annex J
# reference program; they are logic/consistency checks, not a numeric oracle.
# ---------------------------------------------------------------------------
def test_separation_frequency_minimum_at_reference() -> None:
    # fD = 21·10^(1.2·|lg(fT/212)|^1.8): the |lg| makes fT = 212 Hz the minimum,
    # where the exponent vanishes and fD = 21 Hz exactly.
    assert two_tone_separation_frequency(212.0) == pytest.approx(21.0)


def test_separation_frequency_matches_reference_program() -> None:
    # Reproduce the DIN 45681 Annex J expression verbatim:
    #   fD = 21 * 10 ^ (1.2 * Abs(Log(fT / 212) / Log(10)) ^ 1.8)
    for ft in (88.0, 137.3, 300.0, 500.0, 999.0):
        expected = 21.0 * 10.0 ** (1.2 * abs(math.log10(ft / 212.0)) ** 1.8)
        assert two_tone_separation_frequency(ft) == pytest.approx(expected)


def test_separation_frequency_grows_either_side_of_reference() -> None:
    # The absolute value gives a symmetric-in-log minimum at 212 Hz.
    base = two_tone_separation_frequency(212.0)
    assert two_tone_separation_frequency(120.0) > base
    assert two_tone_separation_frequency(600.0) > base


def test_separation_frequency_rejects_bad_frequency() -> None:
    with pytest.raises(ValueError, match="tone_frequency"):
        two_tone_separation_frequency(0.0)


def test_resolve_separately_true_when_far_apart() -> None:
    # 200/260 Hz, both < 1000 Hz, |Δf| = 60 Hz ≫ fD(200) ≈ 21 Hz → separate.
    assert resolve_tones_separately(200.0, 260.0, 3.0, 2.0) is True


def test_resolve_separately_false_when_close() -> None:
    # Annex E tones 118.4/137.3 Hz: |Δf| = 18.9 Hz < fD(137.3) ≈ 24.1 Hz →
    # combined. This is consistent with the Annex E FG grouping.
    assert resolve_tones_separately(118.4, 137.3, 4.0, 5.0) is False


def test_resolve_separately_false_at_or_above_1000hz() -> None:
    # The rule only applies when BOTH tones lie below 1000 Hz.
    assert resolve_tones_separately(1200.0, 1400.0, 3.0, 2.0) is False
    assert resolve_tones_separately(900.0, 1100.0, 3.0, 2.0) is False


def test_resolve_separately_uses_more_prominent_tone() -> None:
    # fT is the tone with the larger ΔL; which tone dominates picks the threshold
    # and, at a spacing that straddles the two thresholds, flips the decision.
    f1, f2 = 300.0, 323.1  # |Δf| = 23.1 Hz
    d_lo = two_tone_separation_frequency(f1)  # ≈ 23.02 Hz
    d_hi = two_tone_separation_frequency(f2)  # ≈ 23.91 Hz
    assert d_lo < abs(f2 - f1) < d_hi
    assert resolve_tones_separately(f1, f2, 5.0, 3.0) is True  # fT = 300, thr low
    assert resolve_tones_separately(f1, f2, 3.0, 5.0) is False  # fT = 323.1, thr hi


def test_resolve_separately_rejects_bad_input() -> None:
    with pytest.raises(ValueError):
        resolve_tones_separately(-1.0, 200.0, 1.0, 1.0)
    with pytest.raises(ValueError):
        resolve_tones_separately(200.0, 260.0, math.nan, 1.0)


# ---------------------------------------------------------------------------
# Mean audibility (Formula (20)) and the no-tone convention (Formula (21))
# ---------------------------------------------------------------------------
def test_mean_audibility_matches_annex_e() -> None:
    value = mean_audibility(ref.ISO20065_DECISIVE_AUDIBILITIES)
    assert value == pytest.approx(ref.ISO20065_MEAN_AUDIBILITY, abs=0.05)


def test_mean_audibility_is_energy_mean() -> None:
    deltas = [4.0, 6.0, NO_TONE_AUDIBILITY]
    expected = 10.0 * math.log10(np.mean([10.0 ** (d / 10.0) for d in deltas]))
    assert mean_audibility(deltas) == pytest.approx(expected)


def test_mean_audibility_between_min_and_max() -> None:
    deltas = ref.ISO20065_DECISIVE_AUDIBILITIES
    value = mean_audibility(deltas)
    assert min(deltas) <= value <= max(deltas)


# ---------------------------------------------------------------------------
# assess_tones result
# ---------------------------------------------------------------------------
def _annex_e_result() -> ToneAudibilityResult:
    freqs = [t[0] for t in ref.ISO20065_ANNEX_E_TONES]
    ls = [t[1] for t in ref.ISO20065_ANNEX_E_TONES]
    lt = [t[2] for t in ref.ISO20065_ANNEX_E_TONES]
    return assess_tones(freqs, lt, ls, ref.ISO20065_LINE_SPACING)


def test_assess_tones_reproduces_annex_e_column() -> None:
    result = _annex_e_result()
    expected = np.array([t[3] for t in ref.ISO20065_ANNEX_E_TONES])
    assert np.max(np.abs(result.audibilities - expected)) < 0.03


def test_assess_tones_decisive() -> None:
    result = _annex_e_result()
    # Per-tone decisive of spectrum 1 is the 137.3 Hz tone (ΔL ≈ 4.99 dB).
    assert result.decisive_frequency == pytest.approx(137.3)
    assert result.decisive_audibility == pytest.approx(4.99, abs=0.03)
    assert result.decisive_audibility == float(np.max(result.audibilities))


def test_assess_tones_audible_mask() -> None:
    result = _annex_e_result()
    assert result.audible.all()  # all nine Annex E tones have ΔL > 0
    quiet = assess_tones([200.0], [40.0], [50.0], 2.7)
    assert not quiet.audible.any()


def test_assess_tones_fields_consistent() -> None:
    result = _annex_e_result()
    assert result.lower_corners.shape == result.tone_frequencies.shape
    # LG = LS + 10 lg(Δfc/Δf) elementwise.
    expected_lg = result.mean_narrowband_levels + 10.0 * np.log10(
        result.critical_bandwidths / result.line_spacing
    )
    assert np.allclose(result.critical_band_levels, expected_lg)
    # ΔL = LT − LG − av elementwise.
    assert np.allclose(
        result.audibilities,
        result.tone_levels - result.critical_band_levels - result.masking_indices,
    )


def test_plot_smoke() -> None:
    ax = _annex_e_result().plot()
    assert ax.get_ylabel()
    import matplotlib.pyplot as plt

    plt.close("all")


# ---------------------------------------------------------------------------
# Validation
# ---------------------------------------------------------------------------
@pytest.mark.parametrize("bad", [0.0, -1.0, float("nan"), float("inf")])
def test_critical_bandwidth_rejects_bad_frequency(bad: float) -> None:
    with pytest.raises(ValueError):
        critical_bandwidth_engineering(bad)


def test_critical_band_level_rejects_bad_spacing() -> None:
    with pytest.raises(ValueError):
        critical_band_level(50.0, 137.3, 0.0)


def test_audibility_rejects_non_finite() -> None:
    with pytest.raises(ValueError):
        audibility_from_levels(float("nan"), 60.0, -2.0)


def test_energy_sum_rejects_empty() -> None:
    with pytest.raises(ValueError):
        energy_sum_level([])


def test_mean_audibility_rejects_empty() -> None:
    with pytest.raises(ValueError):
        mean_audibility([])


def test_assess_tones_rejects_length_mismatch() -> None:
    with pytest.raises(ValueError, match="share their length"):
        assess_tones([100.0, 200.0], [60.0], [50.0], 2.7)


def test_assess_tones_rejects_empty() -> None:
    with pytest.raises(ValueError):
        assess_tones([], [], [], 2.7)


def test_assess_tones_rejects_non_positive_frequency() -> None:
    with pytest.raises(ValueError, match="positive"):
        assess_tones([0.0], [60.0], [50.0], 2.7)
