"""Standalone measurement functions for the published chart-codec study.

Snapshot from projects/audio.py, 2026-09-26. The application upload-size guard
is omitted so these functions require no Django settings. Calculations are
unchanged. Keep this file beside the other three downloaded scripts.
Requires ffmpeg, numpy and soundfile. All input audio remains local.
"""
import os
import logging
from contextlib import contextmanager


def _resample_loudness_series(pairs, n=2048, max_gap_s=1.0):
    """Bucket (time, short-term-LUFS) pairs into ``n`` evenly-timed buckets
    (averaging within each), aligned to the waveform's time axis. Buckets with
    no data; silent stretches, the initial 3 s gating ramp (values < -70 are
    dropped); are ``None`` so the front-end breaks the line there rather than
    diving to -120. Returns a list of ``n`` floats/None, or None if no usable
    data.

    ffmpeg emits a short-term value only ~every 0.1 s, so a short track has far
    fewer source points than ``n`` buckets; naive bucketing then leaves ~half
    the buckets empty, interspersed, and the front-end breaks the line at every
    hole → the curve renders as a dotted line. So we linearly interpolate across
    SHORT interior gaps (resolution artefacts, < ``max_gap_s``) and keep only
    long gaps (real silence) as ``None`` for the line to break there."""
    if not pairs:
        return None
    dur = pairs[-1][0]
    if not dur or dur <= 0:
        return None
    sums = [0.0] * n
    counts = [0] * n
    for t, s in pairs:
        if s <= -70:           # gating ramp / silence placeholder; skip
            continue
        bi = int(t / dur * n)
        if bi < 0:
            bi = 0
        elif bi >= n:
            bi = n - 1
        sums[bi] += s
        counts[bi] += 1
    if not any(counts):
        return None
    out = [round(sums[i] / counts[i], 1) if counts[i] else None for i in range(n)]

    # Bridge short interior gaps so the contour is continuous through dense
    # audio; leave long gaps (real silence) broken.
    bucket_s = dur / n
    max_run = max(1, int(round(max_gap_s / bucket_s))) if bucket_s > 0 else n
    i = 0
    while i < n:
        if out[i] is not None:
            i += 1
            continue
        j = i
        while j < n and out[j] is None:
            j += 1
        # out[i..j-1] is a None run bounded by out[i-1] and out[j]; bridge it
        # only when it sits between two real values and is short enough.
        if i > 0 and j < n and (j - i) <= max_run:
            lo, hi = out[i - 1], out[j]
            span = j - (i - 1)
            for k in range(i, j):
                f = (k - (i - 1)) / span
                out[k] = round(lo * (1 - f) + hi * f, 1)
        i = j
    return out


def analyze_loudness(input_path, timeout=180):
    """EBU R128 loudness + short-term max via ffmpeg ``ebur128`` filter.

    Returns ``{'lufs', 'true_peak', 'lra', 'short_term_max'}`` or None.

    We deliberately do NOT chain ``loudnorm`` here; when its dynamic-norm
    output is piped into ebur128 in the same graph, the short-term LUFS
    stream reflects the *normalized* signal, which understates the master's
    real short-term loudness by ~15 dB. ebur128 alone gives us the original
    integrated I, range LRA, true peak (oversampled) AND the per-time S
    stream we scan for short_term_max.
    """
    import subprocess
    import shutil
    import math
    import re
    log = logging.getLogger('soneam.loudness')

    if not shutil.which('ffmpeg'):
        return None
    try:
        result = subprocess.run([
            'ffmpeg', '-hide_banner', '-nostats',
            '-i', input_path,
            '-af', 'ebur128=peak=true',
            '-f', 'null', '-',
        ], capture_output=True, text=True, timeout=timeout)
        err = result.stderr or ''

        # ── Summary block at end of stderr (multi-line, use DOTALL) ────────
        def _grab(pattern):
            m = re.search(pattern, err, re.DOTALL)
            if not m:
                return None
            try:
                v = float(m.group(1))
            except (TypeError, ValueError):
                return None
            return v if math.isfinite(v) else None

        lufs      = _grab(r'Integrated loudness:.*?I:\s*(-?[\d.]+)\s*LUFS')
        lra       = _grab(r'Loudness range:.*?LRA:\s*(-?[\d.]+)\s*LU')
        # ebur128 prints "True peak: Peak: X dBFS"; that *is* the oversampled
        # true peak; we keep it stored under true_peak_dbtp for consistency
        # with the rest of the app even though ffmpeg labels it dBFS.
        true_peak = _grab(r'True peak:.*?Peak:\s*(-?[\d.]+)\s*dBFS')

        # ── Per-time short-term LUFS (the loudness-over-time curve) ─────────
        # ebur128 logs a line per ~0.1s carrying both t: (seconds) and S:
        # (3-second short-term LUFS). Collect (t, S), then bucket into a fixed
        # series aligned to the waveform so the front-end can draw it as a
        # polyline over the wave (like the TP envelope).
        st_pairs = []
        for m in re.finditer(r't:\s*([\d.]+)[^\n]*?S:\s*(-?[\d.]+)', err):
            try:
                t = float(m.group(1)); s = float(m.group(2))
            except (TypeError, ValueError):
                continue
            if math.isfinite(t) and math.isfinite(s):
                st_pairs.append((t, s))

        # Loudest moment (ignoring the gating ramp / silence padding near -120).
        st_vals = [s for (_, s) in st_pairs if s > -70]
        st_max = max(st_vals) if st_vals else None
        st_series = _resample_loudness_series(st_pairs, n=2048)

        return {
            'lufs':           lufs,
            'true_peak':      true_peak,
            'lra':            lra,
            'short_term_max': st_max,
            'st_series':      st_series,
        }
    except subprocess.TimeoutExpired:
        log.warning('ebur128 timed out for %s', os.path.basename(input_path))
        return None
    except Exception as e:
        log.warning('ebur128 failed for %s: %s', os.path.basename(input_path), e)
        return None


@contextmanager
def _analysis_soundfile(path):
    """Read WAVs whose ancillary chunks libsndfile rejects, without touching
    the original or converting its PCM. ffmpeg rewrites only the container.
    The temporary copy is used solely for metrics and removed on every exit.
    """
    import soundfile as sf
    import subprocess
    import tempfile
    try:
        audio = sf.SoundFile(path)
    except RuntimeError:
        with open(path, 'rb') as source:
            header = source.read(12)
        if header[:4] not in (b'RIFF', b'RF64') or header[8:12] != b'WAVE':
            raise
        with tempfile.TemporaryDirectory(prefix='soneam-analysis-') as directory:
            normalized = os.path.join(directory, 'analysis.wav')
            result = subprocess.run([
                'ffmpeg', '-nostdin', '-v', 'error', '-y', '-i', str(path),
                '-map', '0:a:0', '-c:a', 'copy', '-map_metadata', '-1',
                '-rf64', 'auto', normalized,
            ], capture_output=True, timeout=120)
            if result.returncode:
                raise RuntimeError('WAV analysis container rewrite failed')
            with sf.SoundFile(normalized) as audio:
                yield audio
        return
    with audio:
        yield audio


def analyze_dynamic_range(input_path):
    """Approximate TT-style Dynamic Range from peak/RMS of the loudest blocks.

    Splits the file into ~3-second blocks, computes peak and RMS per block,
    sorts by RMS, takes the top 20%, and returns 20·log10(mean(peak)/mean(RMS)) for
    that subset. The result is in dB; typical commercial masters land in
    the 4–10 dB range; dynamic mixes show 12–20 dB.

    Not a strict TT-DR (which has additional smoothing & 3 s vs 20 s rules
    that vary by spec version), but close enough for engineer guidance and
    free since we already read the file for peaks generation.
    """
    import math
    log = logging.getLogger('soneam.dr')
    try:
        import soundfile as sf
        import numpy as np
    except ImportError:
        return None

    try:
        with _analysis_soundfile(input_path) as f:
            sr = f.samplerate
            channels = f.channels
            if sr <= 0 or channels <= 0:
                return None
            block = max(1, int(sr * 3))                # 3 s blocks
            peaks, rmses = [], []
            while True:
                chunk = f.read(block, dtype='float32', always_2d=True)
                if chunk.shape[0] == 0:
                    break
                # Combine channels: max(abs) for peak, sqrt(mean(square)) for RMS.
                p = float(np.max(np.abs(chunk)))
                r = float(np.sqrt(np.mean(chunk ** 2)))
                if p > 0 and r > 0:
                    peaks.append(p); rmses.append(r)
                if chunk.shape[0] < block:
                    break

            if not rmses:
                return None
            # Top 20% by RMS (the loudest stretches drive perceived dynamics).
            n = max(1, len(rmses) // 5)
            order = np.argsort(rmses)[::-1][:n]
            top_p = float(np.mean([peaks[i] for i in order]))
            top_r = float(np.mean([rmses[i] for i in order]))
            if top_r <= 0 or top_p <= 0:
                return None
            return 20.0 * math.log10(top_p / top_r)
    except Exception as e:
        log.warning('DR analysis failed for %s: %s', os.path.basename(input_path), e)
        return None


def generate_track_peaks(file_path, n=2048):
    """Return waveform buckets, sample peak and a 4× FFT true-peak estimate.

    Snapshot of the estimator used for the study. ffmpeg ebur128 is recorded
    independently; neither estimator is treated as evidence of audibility.
    """
    import math
    log = logging.getLogger('soneam.peaks')

    try:
        import soundfile as sf
        import numpy as np
    except ImportError:
        log.warning('soundfile/numpy not installed; cannot generate peaks')
        return None, None, None, None

    try:
        with _analysis_soundfile(file_path) as f:
            channels = f.channels
            total_frames = f.frames
            sr = f.samplerate
            if total_frames == 0 or sr <= 0:
                return None, None, None, None
            bucket_size = max(1, total_frames // n)
            actual_n = total_frames // bucket_size
            n_ch = min(channels, 2)

            peaks = [[] for _ in range(n_ch)]
            tp_bucket = []           # max true peak per bucket (linear, max across L/R)
            global_peak = 0.0
            global_tp = 0.0
            UP = 4  # 4× oversample for ITU-R BS.1770-style true peak.
            # Context samples read on each side of a bucket and upsampled with it,
            # then DISCARDED. A per-bucket FFT zero-pad upsample treats the bucket
            # as periodic; the wrap-around discontinuity turns into Gibbs ringing
            # that inflated the true peak by ~2 dB. Upsampling with margins and
            # keeping only the bucket interior removes that artifact.
            MARGIN = 128

            for i in range(actual_n):
                start = i * bucket_size
                rstart = max(0, start - MARGIN)
                pre = start - rstart                       # left context actually available
                f.seek(rstart)
                ext = f.read(bucket_size + pre + MARGIN, dtype='float32', always_2d=True)
                if ext.shape[0] == 0:
                    break
                blen = min(bucket_size, ext.shape[0] - pre)  # real bucket length within ext
                if blen <= 0:
                    break
                real = ext[pre:pre + blen]

                max_abs = np.max(np.abs(real), axis=0)
                for ch_idx in range(n_ch):
                    peaks[ch_idx].append(float(max_abs[ch_idx]))
                local_peak = float(np.max(max_abs[:n_ch]))
                if local_peak > global_peak:
                    global_peak = local_peak

                # True peak per bucket: upsample the EXTENDED block, read the peak
                # only from the real-bucket interior (margins discarded). Floor at
                # the bucket's sample peak (TP can never be below it).
                local_tp = local_peak
                if ext.shape[0] >= 8:
                    lo, hi = pre * UP, (pre + blen) * UP
                    for ch in range(n_ch):
                        col = ext[:, ch].astype(np.float64, copy=False)
                        spec = np.fft.rfft(col)
                        padded = np.zeros(col.size * UP // 2 + 1, dtype=complex)
                        padded[:spec.shape[0]] = spec
                        up = np.fft.irfft(padded, n=col.size * UP) * UP
                        p = float(np.max(np.abs(up[lo:hi])))
                        if p > local_tp:
                            local_tp = p
                    if local_tp > global_tp:
                        global_tp = local_tp
                tp_bucket.append(round(local_tp, 4))

            global_tp = max(global_tp, global_peak)   # TP ≥ sample peak, always

            result = {
                'l': peaks[0],
                'r': peaks[1] if channels >= 2 else None,
                'n': len(peaks[0]),
            }
            tp_result = {'t': tp_bucket, 'n': len(tp_bucket)} if tp_bucket else None
            sample_peak_dbfs = 20.0 * math.log10(global_peak) if global_peak > 0 else -120.0
            true_peak_dbtp   = 20.0 * math.log10(global_tp)   if global_tp   > 0 else -120.0
            return result, sample_peak_dbfs, true_peak_dbtp, tp_result
    except Exception as e:
        log.warning('peaks generation failed for %s: %s', os.path.basename(file_path), e)
        return None, None, None, None