#!/usr/bin/env python3
"""
Every number spoken in the long-form film "every chord you have ever heard is out of tune".

Brief: blender/content/briefs/tuning-temperament-01.md   Content Bank #46

Two jobs:
  1. Compute the temperament arithmetic (commas, fifths, Vallotti, predicted beat rates).
  2. RENDER each interval to audio and MEASURE the beat rate that actually comes out,
     because every figure in job 1 is a prediction until something counts it.

Stdlib only, on purpose — the film tells the audience this is short enough to read in one sitting.

  python3 tools/temperament.py --selftest    # prove the measurement works before trusting it
  python3 tools/temperament.py --report      # predicted vs MEASURED for every spoken figure
  python3 tools/temperament.py --render DIR  # write the demo wavs
"""
from __future__ import annotations

import argparse
import array
import math
import os
import sys
import wave
from fractions import Fraction as F

SR = 48000   # the recording convention for demos is 48 kHz (RECORDING-QUEUE.md)
C4 = 261.6255653005986  # A4 = 440, equal temperament

# ---------------------------------------------------------------- the arithmetic

PYTH_COMMA = F(3, 2) ** 12 / 2 ** 7          # 531441:524288
SYNT_COMMA = F(81, 80)
SCHISMA = PYTH_COMMA / SYNT_COMMA


def cents(ratio) -> float:
    return 1200.0 * math.log2(float(ratio))


def quarter_comma_fifth() -> float:
    """1/4-comma meantone fifth, in cents. Four of these minus two octaves is a pure 5:4."""
    return cents(5 ** 0.25)


def vallotti_scale() -> dict:
    """Vallotti as it is USUALLY implemented: six fifths narrowed by 1/6 PYTHAGOREAN comma.

    ⚠ This is NOT what Vallotti wrote. His own text (Trattato della moderna musica, Padova
    1950, Book II pp.195-197) gives the D-A fifth as 27:40, which is exactly one SYNTONIC
    comma from pure — see vallotti_original_da() below. The Pythagorean-comma version is
    Thomas Young's (1800); Barbour p.163 describes it as Young's and never mentions Vallotti.
    """
    pc, p5 = cents(PYTH_COMMA), cents(F(3, 2))
    chain = ["F", "C", "G", "D", "A", "E", "B", "F#", "C#", "G#", "D#", "A#"]
    tempered = {("F", "C"), ("C", "G"), ("G", "D"), ("D", "A"), ("A", "E"), ("E", "B")}
    pitch = {"F": 0.0}
    for a, b in zip(chain, chain[1:]):
        pitch[b] = pitch[a] + p5 - (pc / 6 if (a, b) in tempered else 0.0)
    c0 = pitch["C"]
    return {n: (p - c0) % 1200 for n, p in pitch.items()}


def vallotti_original_da():
    """Vallotti's OWN D-A fifth, and how far it sits from pure. Returns (cents, comma_ratio)."""
    da = F(40, 27)               # 27:40 as a length ratio is 40:27 as a frequency ratio
    return cents(da), F(3, 2) / da   # the quotient is exactly 81:80


def beat_rate(root_hz: float, third_cents: float) -> float:
    """Beats/sec between the root's 5th harmonic and the third's 4th harmonic."""
    third = root_hz * 2 ** (third_cents / 1200.0)
    return abs(5 * root_hz - 4 * third)


def fifth_beat_rate(root_hz: float, fifth_cents: float) -> float:
    """Beats/sec between the root's 3rd harmonic and the fifth's 2nd harmonic."""
    fifth = root_hz * 2 ** (fifth_cents / 1200.0)
    return abs(3 * root_hz - 2 * fifth)


# ---------------------------------------------------------------- synthesis

def complex_tone(f0: float, n: int, dur: float, rolloff: float = 0.75) -> list:
    out = [0.0] * int(SR * dur)
    for k in range(1, n + 1):
        if f0 * k >= SR / 2:
            break
        a, w = rolloff ** (k - 1), 2 * math.pi * f0 * k / SR
        for i in range(len(out)):
            out[i] += a * math.sin(w * i)
    return out


def interval(root_hz: float, interval_cents: float, dur: float = 3.0, n: int = 8) -> list:
    up = root_hz * 2 ** (interval_cents / 1200.0)
    a, b = complex_tone(root_hz, n, dur), complex_tone(up, n, dur)
    mix = [(x + y) for x, y in zip(a, b)]
    peak = max(abs(v) for v in mix) or 1.0
    return [0.7 * v / peak for v in mix]


def write_wav(path: str, samples: list) -> None:
    data = array.array("h", (int(max(-1.0, min(1.0, s)) * 32767) for s in samples))
    with wave.open(path, "wb") as w:
        w.setnchannels(1)
        w.setsampwidth(2)
        w.setframerate(SR)
        w.writeframes(data.tobytes())


# ---------------------------------------------------------------- measurement

def envelope(x: list, cutoff: float = 60.0, decim: int = 64) -> tuple:
    """Rectify, one-pole lowpass, decimate. Returns (env, effective_rate)."""
    a = math.exp(-2 * math.pi * cutoff / SR)
    y, out = 0.0, []
    for i, s in enumerate(x):
        y = (1 - a) * abs(s) + a * y
        if i % decim == 0:
            out.append(y)
    return out, SR / decim


def dominant_modulation(env: list, rate: float, lo: float = 0.5, hi: float = 60.0,
                        step: float = 0.02) -> float:
    """Strongest modulation frequency in the envelope, by Goertzel scan. Hz."""
    mean = sum(env) / len(env)
    x = [v - mean for v in env]
    n = len(x)
    win = [0.5 - 0.5 * math.cos(2 * math.pi * i / (n - 1)) for i in range(n)]
    x = [v * w for v, w in zip(x, win)]
    best_f, best_m = 0.0, -1.0
    f = lo
    while f <= hi:
        w = 2 * math.pi * f / rate
        c = 2 * math.cos(w)
        s1 = s2 = 0.0
        for v in x:
            s0 = v + c * s1 - s2
            s2, s1 = s1, s0
        m = s1 * s1 + s2 * s2 - c * s1 * s2
        if m > best_m:
            best_m, best_f = m, f
        f += step
    return best_f


def resonator(x: list, centre: float, bw: float) -> list:
    """Two-pole bandpass. Isolates ONE coinciding-harmonic region out of the full mix.

    Why this exists: the first version of measure_beat_rate() enveloped the whole mix and
    read 35 beats/s off a PURE 5:4, which does not beat at all. With 8 partials per tone the
    envelope carries a difference tone from every pair — the loudest being the gcd of the two
    fundamentals (65 Hz for a 5:4). The film does not claim the mix wobbles; VO_04c claims a
    specific pair collides: the root's 5th harmonic against the third's 4th. So measure THAT.
    """
    r = math.exp(-math.pi * bw / SR)
    w = 2 * math.pi * centre / SR
    a1, a2 = 2 * r * math.cos(w), -r * r
    b0, b2 = 1 - r, -(1 - r)
    y1 = y2 = x1 = x2 = 0.0
    out = []
    for s in x:
        y = b0 * s + b2 * x2 + a1 * y1 + a2 * y2
        x2, x1 = x1, s
        y2, y1 = y1, y
        out.append(y)
    return out


def measure_beat(root_hz: float, interval_cents: float, dur: float = 3.0,
                 kind: str = "third") -> tuple:
    """Returns (beats_per_sec, modulation_depth) at the harmonic pair the film names.

    ⚠ DEPTH IS NOT OPTIONAL. When an interval does not beat, the envelope is flat and the
    "dominant modulation frequency" is UNDEFINED, not zero — the scan returns whatever noise
    peak happens to win. An earlier selftest asserted that a pure 5:4 must measure ~0 Hz and
    it read 35 Hz, which looked like a synthesis bug and was actually a bad question.
    Depth answers "is it beating at all"; frequency answers "how fast" and is only meaningful
    once depth says yes. → feedback_a_binary_question_has_no_answer_during_the_gap
    """
    centre = 5 * root_hz if kind == "third" else 3 * root_hz
    x = interval(root_hz, interval_cents, dur)
    band = resonator(x, centre, bw=120.0)
    band = band[int(0.3 * SR):]          # let the resonator settle
    env, rate = envelope(band)
    mean = sum(env) / len(env)
    var = sum((v - mean) ** 2 for v in env) / len(env)
    depth = (var ** 0.5) / mean if mean > 1e-12 else 0.0
    return dominant_modulation(env, rate), depth


def measure_beat_rate(root_hz: float, interval_cents: float, dur: float = 3.0,
                      kind: str = "third") -> float:
    return measure_beat(root_hz, interval_cents, dur, kind)[0]


# ---------------------------------------------------------------- the spoken figures

SPOKEN = [
    # (take, label, interval cents, predicted beats/s, kind)
    ("VO_04c", "Pythagorean major 3rd on C4", 407.8200, None, "third"),
    ("VO_05d", "meantone WOLF fifth on C4", 737.6373, None, "fifth"),
    ("VO_06b", "Vallotti C-E", 392.1800, None, "third"),
    ("VO_06b", "Vallotti F#-A#", 407.8200, None, "third"),
    ("VO_09b", "equal-tempered 5th on C4", 700.0000, None, "fifth"),
    ("VO_09c", "equal-tempered major 3rd on C4", 400.0000, None, "third"),
]


def report() -> int:
    print("PREDICTED vs MEASURED — every beat rate the film speaks\n")
    print(f"{'take':<9}{'interval':<32}{'cents':>10}{'predicted':>11}{'measured':>10}{'err':>8}")
    worst = 0.0
    for take, label, ic, _, kind in SPOKEN:
        pred = beat_rate(C4, ic) if kind == "third" else fifth_beat_rate(C4, ic)
        meas = measure_beat_rate(C4, ic, kind=kind)
        err = abs(meas - pred)
        worst = max(worst, err)
        print(f"{take:<9}{label:<32}{ic:>10.2f}{pred:>11.2f}{meas:>10.2f}{err:>8.2f}")
    print(f"\nworst absolute error: {worst:.3f} Hz")
    et3, et5 = beat_rate(C4, 400.0), fifth_beat_rate(C4, 700.0)
    print(f"\nVO_09c ratio claim: ET 3rd / ET 5th = {et3:.2f} / {et5:.2f} = {et3/et5:.1f}x")
    print("VO_09e octave doubling:", " -> ".join(
        f"{beat_rate(C4 * 2 ** k, 400.0):.1f}" for k in range(3)), "beats/s")
    return 0 if worst < 0.5 else 1


# ---------------------------------------------------------------- selftest

def selftest() -> int:
    fails = []

    # 1. the measurement instrument must recover a KNOWN modulation before it is trusted.
    for want in (3.0, 10.0, 16.0):
        n = int(SR * 3.0)
        sig = [(1.0 + 0.9 * math.sin(2 * math.pi * want * i / SR))
               * math.sin(2 * math.pi * 440 * i / SR) for i in range(n)]
        env, rate = envelope(sig)
        got = dominant_modulation(env, rate)
        if abs(got - want) > 0.2:
            fails.append(f"instrument: planted {want} Hz AM, measured {got:.2f}")
        else:
            print(f"  [ok] instrument recovers a planted {want:.0f} Hz modulation as {got:.2f}")

    # 2. a JUST third must not beat. ask about DEPTH, not frequency (see measure_beat docstring).
    _, pure_depth = measure_beat(C4, cents(F(5, 4)), kind="third")
    _, et_depth = measure_beat(C4, 400.0, kind="third")
    # assert the SEPARATION, not an absolute threshold I guessed. the residual on a pure
    # interval is resonator/rectifier ripple, not beating; what matters is that a beating
    # interval is unmistakably deeper. a constant picked to make the test pass tests nothing.
    if et_depth < 5 * pure_depth:
        fails.append(f"depth does not separate pure ({pure_depth:.3f}) from tempered "
                     f"({et_depth:.3f}) — the instrument cannot tell beating from silence")
    else:
        print(f"  [ok] pure 5:4 is FLAT (depth {pure_depth:.4f}) while ET beats "
              f"(depth {et_depth:.4f}, {et_depth/max(pure_depth,1e-9):.0f}x)")

    # 3. MUTATION — detune the third and demand the measurement follows.
    m = measure_beat_rate(C4, cents(F(5, 4)) + 13.686, kind="third")   # ET's error
    if abs(m - beat_rate(C4, 400.0)) > 0.5:
        fails.append(f"mutation: detuned third measured {m:.2f}, predicted "
                     f"{beat_rate(C4, 400.0):.2f}")
    else:
        print(f"  [ok] mutation caught: detuning to ET gives {m:.2f} beats/s")

    # 4. the arithmetic anchors the film states outright.
    if F(3, 2) / F(40, 27) != SYNT_COMMA:
        fails.append("Vallotti's 27:40 is no longer exactly one syntonic comma from pure")
    else:
        print("  [ok] Vallotti D-A 27:40 is EXACTLY 81:80 from pure (VO_07c)")
    m3 = 4 * quarter_comma_fifth() - 2400.0
    if abs(m3 - cents(F(5, 4))) > 1e-6:
        fails.append(f"meantone identity broke: {m3} vs {cents(F(5, 4))}")
    else:
        print(f"  [ok] 4 meantone 5ths - 2 octaves = {m3:.6f} c = pure 5:4 (VO_05b)")

    if fails:
        print("\nSELFTEST FAILED:")
        for f in fails:
            print("  -", f)
        return 1
    print("\nSELFTEST PASSED")
    return 0


def render(outdir: str) -> int:
    os.makedirs(outdir, exist_ok=True)
    # ⚠ NAMING IS NOT COSMETIC. content/audio/<slug>/ is where JONAS drops his voice takes, and
    # VO_*.wav is reserved for them. Audio examples MUST be DEMO_* or his VO_04c.wav collides with
    # ours. Convention: content/planning/RECORDING-QUEUE.md, "Where you PUT what you record".
    for _, label, ic, _, _ in SPOKEN:
        name = "DEMO_" + label.upper().replace(" ", "_").replace("-", "_").replace("#", "SHARP")
        path = os.path.join(outdir, f"{name}.wav")
        write_wav(path, interval(C4, ic, dur=4.0))
        print("wrote", os.path.basename(path))
    write_wav(os.path.join(outdir, "DEMO_JUST_THIRD_REFERENCE.wav"),
              interval(C4, cents(F(5, 4)), dur=4.0))
    print("wrote DEMO_JUST_THIRD_REFERENCE.wav")
    return 0


def main() -> int:
    ap = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--selftest", action="store_true")
    ap.add_argument("--report", action="store_true")
    ap.add_argument("--render", metavar="DIR")
    a = ap.parse_args()
    if a.selftest:
        return selftest()
    if a.render:
        return render(a.render)
    return report()


if __name__ == "__main__":
    sys.exit(main())