feature_extractor_cpu.py

# feature_extractor_cpu.py — HF Spaces CPU Basic variant.
# Identical to feature_extractor.py. shifts=5 is preserved as required.
# Expected processing time on CPU Basic (2 vCPU): ~33–136 min per track
# depending on track length. This is expected — do not reduce shifts.
# No GPU code. No spaces.GPU decorator. Do not import the 'spaces' package.

import numpy as np
import librosa, torch, openl3, warnings
from demucs.pretrained import get_model
from demucs.apply import apply_model
from scipy.stats import entropy
from scipy.signal import butter, sosfilt, resample_poly
import pyloudnorm as pyln

warnings.filterwarnings('ignore'); EPS = 1e-10

# ─────────────────────────────────────────────────────────────────────────────
# ENG_DIMS — ordered names for every element of the eng_vec returned by
# extract_mir_features(). Used by lane_model.py and analyze_submissions.py
# so that dimension names and indices are never hardcoded in two places.
#
# DESIGN RATIONALE — SPLITTING COMPOSITE DIMENSIONS
# ──────────────────────────────────────────────────
# Every dimension here is intended to be atomic: one measurable quantity,
# one conceptual meaning. Dimensions that previously bundled distinct concepts
# have been decomposed. The affected cases are documented below.
#
# 1. vocal_prominence (old) → vocal_level + vocal_presence + vocal_brightness
#    The old single ratio (vocal RMS / total RMS) conflated three independent
#    questions: how loud are the vocals when they sing? how much of the song
#    has singing at all? how upfront/forward do the vocals sit spectrally?
#    A track can have loud but infrequent vocals (high level, low presence),
#    or frequent but spectrally buried vocals (high presence, low brightness).
#    These are now three separate dimensions.
#
# 2. drum_punch (kept, modified in v4) + drum_prominence (new in v3)
#    drum_punch is the crest factor of the drum stem gated by drum_confidence.
#    drum_prominence is the drum stem's RMS share of the total mix.
#
#    Before v4, drum_punch was a raw crest factor: peak(drum_stem) /
#    rms(drum_stem). This produced false positives for songs with no drums:
#    htdemucs leaks guitar strum transients into the drum stem, giving the
#    stem a tiny RMS denominator but real peaks → inflated crest factor.
#    drum_prominence did not protect against this because it measures the
#    stem's share of the total mix, not the stem's internal peak-to-RMS ratio.
#
#    Fix (v4): drum_punch is now multiplied by drum_confidence [0,1], a
#    frequency-domain drum detection score. Supporting sub-dimensions:
#
#    drum_kick_presence  — (v5, unchanged in v6) Drum stem's dominance in the
#                          40–120 Hz kick band, measured as:
#                            RMS(drum_stem, 40–120 Hz) / RMS(full_mix, 40–120 Hz)
#                          When a real kick drum is present, the drum stem
#                          provides most of the mix's sub-bass energy (ratio
#                          0.50–0.90). Guitar bleed: kick_dom ≈ 0.05–0.30.
#
#    drum_cymbal_presence — (v6 REDESIGN) HPSS-based sustained energy in the
#                          4–10 kHz band of the drum stem.  Measures crash and
#                          ride cymbal ring-out as opposed to transient hi-hat
#                          clicks (which go to hihat_presence below).
#
#                          Formula:
#                            sustained_frac = H_energy(4-10kHz) / total_energy(4-10kHz)
#                            band_prom      = RMS(mag, 4-10kHz) / RMS(mag, all)
#                            cymbal_presence = (sustained_frac × band_prom) / CYMBAL_SAT
#                          (clamped to [0, 1])
#
#                          WHY v5 dominance ratio (8-16 kHz) FAILED:
#                          The v5 formula = RMS(drum_stem, 8-16kHz) / RMS(mix, 8-16kHz)
#                          collapses when a bright guitar or synth pad adds energy
#                          at 8-16 kHz alongside real hi-hats. The mix denominator
#                          inflates while the drum numerator is unchanged, pushing the
#                          ratio below the lane mean even with cymbals throughout.
#                          This caused "cymbal presence thinner than what I feature"
#                          rejections for tracks verified to have hi-hats and crashes.
#
#                          The HPSS H component captures ring-out / shimmer — the
#                          characteristic of crashes and rides — regardless of what
#                          other sources contribute to the full mix at those frequencies.
#
#    hihat_presence      — (v6 NEW) Onset density of the HPSS percussive (P)
#                          component in 8–14 kHz, as a fraction of HIHAT_ONSET_SAT.
#
#                          Hi-hats are short transients → appear in P component.
#                          Guitar chord rings at 8-14 kHz are sustained → appear in H.
#                          Soft-masking by P/(mag+EPS) then zeroing non-hihat bins
#                          isolates percussive 8-14 kHz content. Onset detection on
#                          this signal counts hi-hat events per second.
#
#                          Guitar staccato transients at 8-14 kHz are much weaker
#                          than real hi-hat transients and fall below the
#                          onset_detect delta threshold.
#
#    drum_confidence     — (v6 updated) Composite [0,1] drum detection score:
#                          0.55 × clamp(kick_dom / DRUM_KICK_SAT, 0, 1)
#                        + 0.30 × hihat_presence   (already in [0,1])
#                        + 0.15 × cymbal_presence  (already in [0,1])
#
#                          Hi-hat weight raised vs. v5 cymbal weight (0.30 vs 0.40)
#                          because onset density at 8-14 kHz is near-impossible to
#                          trigger without real drum content. A drumless acoustic
#                          track cannot produce rapid regular transients in that band.
#
#    drum_punch (index 2): drum_crest_raw × drum_confidence (unchanged gate mechanism).
#
# 3. bass_tightness (kept) + bass_prominence (new in v3)
#    Identical reasoning to drum_punch/drum_prominence above.
#
# 4. guitar_energy (old) → other_prominence (renamed in v3)
#    With htdemucs_6s (v8), the "other" stem is now cleaner: guitar and piano
#    are extracted into dedicated stems (indices 4 and 5). other_prominence now
#    reflects only synths, strings, pads, and miscellaneous non-vocal content.
#
# 5. egtr_presence + egtr_level (v8 NEW)
#    Two new dimensions operating on the dedicated guitar stem from htdemucs_6s.
#    See below for full rationale.
#
# RENAMES (done in v3):
#   complexity  → harmonic_complexity  (chroma entropy ≠ general complexity)
#   dynamics    → lufs_variation       (std(LUFS) ≠ all kinds of dynamics)
#
# Index history:
#   Indices 0–25: identical to v3.
#   Indices 26–27: added in v4, semantics updated in v5 (dominance ratio).
#   Index 26: unchanged in v6 (kick dominance ratio).
#   Index 27: v6 redesign — was 8-16kHz dominance ratio, now HPSS sustained
#             fraction in 4-10kHz. INCOMPATIBLE with v5 cache entries.
#   Index 28: v6 updated weights (was 0.60×kick+0.40×cymbal).
#   Index 29: v6 NEW — hihat_presence.
#   Indices 30–31: v8 NEW — egtr_presence, egtr_level.
#             htdemucs switched from 4-stem (htdemucs) to 6-stem (htdemucs_6s).
#             total_rms now includes guitar_rms + piano_rms (more accurate
#             stem-sum denominator). All v7 cache entries are incompatible.
#   v9 UPDATE — egtr_presence formula revised (same indices 30–31, new values):
#             _compute_guitar_dims gains a third sub-dimension: tonal_sustain.
#             tonal_sustain = H_mid_energy / (H_mid_energy + P_mid_energy + EPS)
#             normalized to EGTR_SUSTAIN_SAT (0.70). Captures sustained guitar
#             notes (legato lines, held power chords) that are tonal → H
#             component → invisible to picking_density (percussive-only).
#             Constants recalibrated for semi-clean and moderate-tempo guitar:
#               EGTR_PICK_ONSET_SAT  4.0 → 3.0  (was too fast; 3/sec = 8th
#                                                  notes at 90 BPM)
#               EGTR_FLATNESS_SAT    0.25 → 0.18 (was calibrated for heavy
#                                                  distortion; clean/crunch
#                                                  electric sits 0.06–0.14)
#             Weights revised: 0.45×picking + 0.25×flatness + 0.30×sustain
#               (was 0.65×picking + 0.35×flatness — over-weighted fast picking).
#             All v8 cache entries are incompatible (egtr_presence values shift).
# ─────────────────────────────────────────────────────────────────────────────
ENG_DIMS = [
    "lufs",                # 0:  Integrated loudness (LUFS)
    "crest_factor",        # 1:  Peak-to-RMS ratio (global mix)
    "drum_punch",          # 2:  Drum stem crest factor × drum_confidence (gated)
    "drum_level",          # 3:  Drum stem RMS / total mix RMS [v10: renamed from drum_prominence]
    "stereo_width",        # 4:  Mid-side energy ratio
    "brightness",          # 5:  Spectral centroid of mix (Hz)
    "harmonic_complexity", # 6:  Chroma entropy (harmonic diversity only)
    "clip_num_events",     # 7:  Number of separate clipping events
    "clip_total",          # 8:  Total clipped samples
    "clip_ratio",          # 9:  Clipped samples / total samples
    "clip_max_run",        # 10: Longest consecutive clipping run (samples)
    "clip_mean_run",       # 11: Mean clipping run length across all events
    "other_level",         # 12: Other stem (synth/pad/strings) RMS / total mix RMS [v10: renamed from other_prominence]
    "vocal_level",         # 13: Vocal stem RMS / total mix RMS (loudness when present)
    "vocal_presence",      # 14: Fraction of frames with active vocal signal (0–1)
    "vocal_brightness",    # 15: Spectral centroid of vocal stem (Hz)
    "chorus_lift",         # 16: Max short-term LUFS minus median (LU)
    "crispness",           # 17: Zero-crossing rate (high-frequency texture)
    "clarity",             # 18: Spectral flatness (tonal vs. noise-like character)
    "length_s",            # 19: Track length in seconds
    "lufs_variation",      # 20: Std of short-term LUFS
    "peak_lufs",           # 21: Max short-term LUFS
    "rolloff",             # 22: Spectral rolloff (Hz)
    "bass_tightness",      # 23: Bass stem crest factor
    "bass_level",          # 24: Bass stem RMS / total mix RMS [v10: renamed from bass_prominence]
    "phase_corr",          # 25: L/R phase correlation
    # ── Drum sub-dimensions (v4+, semantics per version noted) ───────────────
    "drum_kick_level",     # 26: RMS(drum stem, 40–120 Hz) / RMS(full mix, 40–120 Hz) [v10: renamed from drum_kick_presence]
    "drum_cymbal_presence",# 27: HPSS sustained fraction × band prominence in 4–10 kHz [v6]
    "drum_confidence",     # 28: 0.55×kick_conf + 0.30×hihat_conf + 0.15×cymbal_conf [v6]
    "hihat_presence",      # 29: Onset density of percussive 8–14 kHz drum component [v6 NEW]
    # ── Electric guitar dimensions (v8 NEW — htdemucs_6s guitar stem) ────────
    "egtr_presence",       # 30: Composite electric guitar detection score [0,1] [v8 NEW]
    "egtr_level",          # 31: egtr_presence × (guitar_rms / total_rms) — guitar mix level [v8 NEW]
    # ── Stem presence + level dimensions (v10 NEW) ────────────────────────────
    "drum_presence",       # 32: Fraction of 1-sec chunks with active drum stem signal [v10 NEW]
    "other_presence",      # 33: Fraction of 1-sec chunks with active other stem signal [v10 NEW]
    "bass_presence",       # 34: Fraction of 1-sec chunks with active bass stem signal [v10 NEW]
    "drum_kick_presence",  # 35: Onset density in 40–120 Hz kick band / DRUM_KICK_ONSET_SAT [v10 NEW]
    "drum_cymbal_level",   # 36: RMS(drum stem, 4–10 kHz) / total mix RMS [v10 NEW]
    "hihat_level",         # 37: RMS(drum stem, 8–14 kHz) / total mix RMS [v10 NEW]
]

# ── Vocal activity detection threshold ───────────────────────────────────────
# Frames whose RMS exceeds this value are counted as active vocal frames.
# 0.01 linear ≈ -40 dBFS.
VOCAL_ACTIVITY_THRESHOLD = 0.01

# ── Drum band analysis constants (v6) ────────────────────────────────────────
#
# KICK (dominance ratio, unchanged from v5):
#   DRUM_KICK_SAT = 0.55
#     The kick dominance ratio at which kick_conf saturates at 1.0. Real kick
#     drum: drum stem provides 50–90 % of mix kick-band energy → conf = 1.0.
#     Acoustic guitar bleed: kick_dom ≈ 0.05–0.30 → conf ≤ 0.55.
#
# CYMBAL (HPSS sustained fraction, 4–10 kHz, v6):
#   CYMBAL_BAND_FMIN / FMAX = 4000–10000 Hz
#     Crash and ride cymbal fundamental energy. Crashes ring in 3–8 kHz;
#     rides shimmer in 5–12 kHz. 4-10 kHz captures the core sustained content
#     of both without extending into the hi-hat domain (8-14 kHz).
#   CYMBAL_SAT = 0.20
#     Value of (sustained_frac × band_prom) at which cymbal_presence = 1.0.
#     Estimated: crash-heavy track → sustained_frac ≈ 0.60, band_prom ≈ 0.35
#     → product ≈ 0.21. Kick+snare only → product ≈ 0.05–0.10.
#     Tune up if lane refs with strong crashes score < 0.70; tune down if
#     kick-only tracks score > 0.40.
#
# HI-HAT (HPSS onset density, 8–14 kHz, v6 NEW):
#   HIHAT_BAND_FMIN / FMAX = 8000–14000 Hz
#     Closed and open hi-hat fundamental energy. Concentrated above 8 kHz;
#     stopping at 14 kHz avoids excessive noise-floor sensitivity at the
#     top of the drum stem frequency range.
#   HIHAT_ONSET_SAT = 5.0  (onsets/sec)
#     Onset rate at which hihat_presence = 1.0. 16th notes at 75 bpm = 5/sec;
#     8th notes at 150 bpm = 5/sec. Most steady hi-hat patterns in rock/pop
#     fall in the 3–8 onsets/sec range. Tune down if lane refs with active
#     hi-hats score < 0.60; tune up if drumless tracks score > 0.15.
#
# DRUM CONFIDENCE WEIGHTS (v6):
#   DRUM_KICK_WEIGHT   = 0.55  (was 0.60 in v5)
#   DRUM_HIHAT_WEIGHT  = 0.30  (new; replaces cymbal_weight contribution)
#   DRUM_CYMBAL_WEIGHT = 0.15  (new; crash/ride as secondary gate input)
#   Sum = 1.0 — fully normalized confidence score.
#
#   Hi-hat raised to 0.30: onset density at 8-14 kHz is nearly impossible to
#   produce without a real drum kit. Acoustic guitars, hand percussion, and
#   electronic pads without a dedicated hi-hat pattern will all score near 0.
#   This makes drum_confidence a much stronger gate against false positives.
#
DRUM_KICK_FMIN     = 40.0     # Hz — kick fundamental lower bound
DRUM_KICK_FMAX     = 120.0    # Hz — kick fundamental upper bound
DRUM_KICK_SAT      = 0.55     # kick dominance ratio at kick_conf = 1.0

CYMBAL_BAND_FMIN   = 4000.0   # Hz — crash/ride lower bound
CYMBAL_BAND_FMAX   = 10000.0  # Hz — crash/ride upper bound
CYMBAL_SAT         = 0.20     # (sustained_frac × band_prom) at cymbal_presence = 1.0

HIHAT_BAND_FMIN    = 8000.0   # Hz — hi-hat lower bound
HIHAT_BAND_FMAX    = 14000.0  # Hz — hi-hat upper bound
HIHAT_ONSET_SAT    = 5.0      # onsets/sec at hihat_presence = 1.0

DRUM_KICK_WEIGHT   = 0.55     # weight for kick_conf in drum_confidence
DRUM_HIHAT_WEIGHT  = 0.30     # weight for hihat_presence in drum_confidence
DRUM_CYMBAL_WEIGHT = 0.15     # weight for cymbal_presence in drum_confidence

# ── Electric guitar detection constants (v9) ──────────────────────────────────
#
# EGTR_MID_FMIN / FMAX = 200–5000 Hz
#   Core range of electric guitar fundamentals and dominant harmonics.
#   Low E string = 82 Hz; high E string = 1318 Hz. Upper harmonics and
#   distortion products extend to 5 kHz. Starting at 200 Hz avoids bleed
#   from the bass stem into the guitar stem. Stopping at 5 kHz keeps the
#   measure out of the hi-hat / pick-noise band (handled by drum dims).
#
# EGTR_PICK_ONSET_SAT = 3.0 (onsets/sec)  [v9: was 4.0]
#   Onset rate at which picking_density = 1.0. 8th notes at 90 bpm = 3/sec;
#   16th notes at 45 bpm = 3/sec. 4.0 was calibrated for fast-picking rock
#   styles; moderate-tempo rhythm guitar (power chords, palm-muted chugging)
#   sits at 2–3 onsets/sec and was being under-scored.
#   Tune DOWN further if lane refs with active rhythm guitar score < 0.50.
#   Tune UP if piano-only or pad-heavy tracks score above 0.30.
#
# EGTR_FLATNESS_SAT = 0.18  [v9: was 0.25]
#   Mid-band spectral flatness at which mid_saturation = 1.0.
#   Clean electric guitar (Fender/Strat clean tone): 0.06–0.14.
#   Semi-clean / light crunch: 0.08–0.18.
#   Crunchy/overdriven: 0.12–0.28. Full distortion: 0.20–0.40.
#   Piano: 0.02–0.08. Acoustic guitar (clean): 0.03–0.10.
#   0.25 was only saturable by heavily overdriven guitar — clean and semi-clean
#   electric guitar scored mid_saturation ≤ 0.40 even with guitar throughout.
#   0.18 rewards any electric guitar content, not only distorted playing.
#
# EGTR_SUSTAIN_SAT = 0.70  [v9 NEW]
#   H-component energy fraction at which tonal_sustain = 1.0.
#   tonal_sustain = H_mid_energy / (H_mid_energy + P_mid_energy)
#   where H = HPSS sustained component, P = HPSS percussive component,
#   both measured in the 200–5000 Hz mid-band of the guitar stem.
#
#   WHY THIS DIMENSION:
#   picking_density and mid_saturation both measure the *attack/transient*
#   character of guitar. A guitar playing sustained power chords, legato
#   lines, or moderate-tempo chord progressions puts most of its energy into
#   the H (sustained) component — the P (percussive) component only captures
#   the brief pick attack. The entire sustain period is invisible to
#   picking_density and adds no flatness signal. tonal_sustain explicitly
#   captures this: as long as guitar notes are ringing, H_mid_energy is high
#   relative to P_mid_energy, and tonal_sustain stays elevated regardless of
#   tempo or gain stage.
#
#   Expected ranges (guitar stem of htdemucs_6s):
#     Electric guitar (any gain, any tempo): H_frac ≈ 0.60–0.85 → sustain 0.86–1.0
#     Acoustic guitar:                       H_frac ≈ 0.55–0.75 → sustain 0.79–1.0
#     Drum bleed in guitar stem:             H_frac ≈ 0.15–0.35 → sustain 0.21–0.50
#     Low-level bleed / near-silence:        Both near EPS → ratio noisy but
#                                            egtr_level ≈ 0 (guitar_rms near 0)
#
#   0.70 saturates at an H_frac consistent with guitar that is sustaining
#   cleanly. Tracks with only transient bleed (drums) will score 0.21–0.50.
#   Tune UP if non-guitar tracks score tonal_sustain > 0.60 consistently.
#
# WEIGHTS (v9):
#   EGTR_PICK_WEIGHT    = 0.45  (was 0.65) — rhythmic attack events
#   EGTR_FLAT_WEIGHT    = 0.25  (was 0.35) — harmonic saturation / gain stage
#   EGTR_SUSTAIN_WEIGHT = 0.30  (NEW)      — sustained tonal guitar content
#   Sum = 1.00
#
#   Rationale for reducing EGTR_PICK_WEIGHT from 0.65:
#   The previous 0.65 weight made picking_density the near-exclusive gate.
#   A track with guitar throughout but moderate tempo and clean tone could
#   score egtr_presence ≈ 0.62 even with guitar clearly audible in the mix,
#   because picking_density was ~0.70 and mid_saturation was ~0.35.
#   Redistributing weight to tonal_sustain captures guitar presence more
#   robustly across gain stages and tempos.
#
# TUNING GUIDANCE (post first training run):
#   Check the printed standards table for egtr_presence and egtr_level.
#   Lane refs with active electric guitar: egtr_presence ≥ 0.65 (ideally 0.75–0.95)
#   Non-lane refs without guitar: egtr_presence ≤ 0.25
#   If guitar tracks still score 0.30–0.50: lower EGTR_PICK_ONSET_SAT (3.0 → 2.0)
#     and/or lower EGTR_FLATNESS_SAT (0.18 → 0.12).
#   If piano/pad tracks score above 0.40: raise EGTR_SUSTAIN_SAT (0.70 → 0.80)
#     or reduce EGTR_SUSTAIN_WEIGHT (0.30 → 0.20).
#
EGTR_MID_FMIN        = 200.0   # Hz — guitar mid-band lower bound
EGTR_MID_FMAX        = 5000.0  # Hz — guitar mid-band upper bound
EGTR_PICK_ONSET_SAT  = 3.0     # onsets/sec at picking_density = 1.0  [v9: was 4.0]
EGTR_FLATNESS_SAT    = 0.18    # flatness value at mid_saturation = 1.0  [v9: was 0.25]
EGTR_SUSTAIN_SAT     = 0.70    # H-frac at tonal_sustain = 1.0  [v9 NEW]
EGTR_PICK_WEIGHT     = 0.45    # weight for picking_density in egtr_presence  [v9: was 0.65]
EGTR_FLAT_WEIGHT     = 0.25    # weight for mid_saturation in egtr_presence  [v9: was 0.35]
EGTR_SUSTAIN_WEIGHT  = 0.30    # weight for tonal_sustain in egtr_presence  [v9 NEW]

# ── Stem temporal presence + kick onset constants (v10) ───────────────────────
DRUM_ACTIVITY_THRESHOLD  = 0.01   # RMS threshold for an active drum chunk
BASS_ACTIVITY_THRESHOLD  = 0.01   # RMS threshold for an active bass chunk
OTHER_ACTIVITY_THRESHOLD = 0.01   # RMS threshold for an active other chunk
DRUM_KICK_ONSET_SAT      = 3.0    # kick onsets/sec at drum_kick_presence = 1.0


device = "cuda" if torch.cuda.is_available() else "cpu"
demucs_model = get_model('htdemucs_6s')   # v8: 6-stem model (adds guitar + piano stems)
demucs_model.to(device)

# htdemucs_6s stem order: 0=drums, 1=bass, 2=other, 3=vocals, 4=guitar, 5=piano
_STEM_DRUMS  = 0
_STEM_BASS   = 1
_STEM_OTHER  = 2
_STEM_VOCALS = 3
_STEM_GUITAR = 4   # v8 NEW — dedicated electric/acoustic guitar stem
_STEM_PIANO  = 5   # v8 NEW — dedicated piano stem


def get_stems(path):
    wav, _ = librosa.load(path, sr=demucs_model.samplerate, mono=False)
    if wav.ndim == 1: wav = np.stack([wav, wav])
    ref = torch.from_numpy(wav).to(device)
    ref_std = ref.std() + EPS
    ref = (ref - ref.mean()) / ref_std
    with torch.no_grad():
        sources = apply_model(demucs_model, ref.unsqueeze(0), shifts=5, split=True)[0]
    return (sources * ref_std.cpu().item()).cpu().numpy()


def _measure_clipping(y, threshold=0.999):
    """
    Returns all five clipping dimensions as individual floats.

    threshold=0.999 (~-0.009 dBFS) avoids flagging legitimately loud
    transient peaks that just touch near-full-scale without flat-topping.

    Returns: (num_events, total_samples, ratio, max_run, mean_run)
    """
    clipped = np.abs(y) >= threshold
    total   = int(np.sum(clipped))
    if total == 0:
        return 0.0, 0.0, 0.0, 0.0, 0.0

    padded      = np.concatenate(([False], clipped, [False]))
    diff        = np.diff(padded.astype(int))
    run_starts  = np.where(diff == 1)[0]
    run_ends    = np.where(diff == -1)[0]
    run_lengths = run_ends - run_starts

    return (
        float(len(run_lengths)),
        float(total),
        float(total / len(y)),
        float(np.max(run_lengths)),
        float(np.mean(run_lengths)),
    )


def _band_rms_abs(y: np.ndarray, sr: int, fmin: float, fmax: float) -> float:
    """
    Absolute RMS of y in the [fmin, fmax] Hz band (not normalized).
    Used by _compute_drum_bands for kick dominance ratio.
    Returns 0.0 if signal is silent or band is invalid.
    """
    nyq = sr / 2.0
    lo  = max(fmin, 1.0)
    hi  = min(fmax, nyq - 1.0)
    if lo >= hi:
        return 0.0
    sos      = butter(2, [lo, hi], btype='bandpass', fs=sr, output='sos')
    filtered = sosfilt(sos, y)
    return float(np.sqrt(np.mean(filtered**2)))


# ─────────────────────────────────────────────────────────────────────────────
# ITU-R BS.1770-4 / BS.1775-2  LOUDNESS AND TRUE PEAK IMPLEMENTATION
#
# WHY THE OLD FORMULA WAS WRONG
# ───────────────────────────────
# The original approximation used:
#   lufs ≈ -0.691 + 10 × log10(mean(y²))
# This is RMS-based power in dB, not integrated LUFS. It is missing two
# mandatory stages defined in ITU-R BS.1770:
#
#   1. K-weighting pre-filter — a two-stage IIR filter that:
#        Stage 1: high-shelf boost (+4 dB at high frequencies) modelling
#                 the acoustic effect of the listener's head on a microphone.
#        Stage 2: high-pass filter (RLB weighting) removing DC and sub-20 Hz
#                 content that inflates RMS but carries no perceived loudness.
#        Together these raise typical music by +1 to +5 dB before squaring.
#
#   2. Gating — BS.1770-4 integrates only "active" blocks:
#        Absolute gate: blocks below -70 LUFS excluded (silence/noise floor).
#        Relative gate: second pass excludes blocks more than -10 LU below
#                       the ungated mean (rejects very quiet passages).
#        Without gating, quiet intros/outros and noise floors pull the
#        integrated value down, causing systematic under-reading of 2–5 dB.
#
# PRACTICAL IMPACT ON THIS SYSTEM
# ────────────────────────────────
# Error: 2–4 dB under-reading of integrated LUFS, varying by program material.
#   • The Spotify normalization gain calculation was off by that same margin.
#   • Lane engineering standards for lufs and peak_lufs were calibrated against
#     the wrong values (consistently wrong, so relative comparisons partially
#     survived — but absolute statements about headroom were incorrect).
#   • short-term LUFS (used for peak_lufs, lufs_variation, chorus_lift) had the
#     same error AND lacked K-weighting, so spectral character affected the
#     "loudness" reading in a non-perceptual way.
#   • true_peak_dbfs was never returned by the modal worker, making the
#     tp_headroom constraint in apply_streaming_normalization() dead code.
#
# FIX SUMMARY (v7)
# ────────────────
# _compute_integrated_lufs()      — pyloudnorm Meter.integrated_loudness()
#                                   Full BS.1770-4: K-weighting + gating.
# _compute_short_term_lufs_array() — K-weighted sliding 3s/1s windows, no gate.
#                                   Same K-weighting as integrated. Replaces the
#                                   raw librosa.feature.rms()-based st_lufs blocks.
# _true_peak_dbfs()               — 4× polyphase oversampling (BS.1775-2) to
#                                   detect intersample peaks. Returned as third
#                                   element of extract_mir_features() so the modal
#                                   worker can include it in the response dict and
#                                   apply_streaming_normalization() can enforce the
#                                   -1.0 dBFS true peak ceiling correctly.
#
# NOTE: pyloudnorm must be installed:  pip install pyloudnorm
# ─────────────────────────────────────────────────────────────────────────────

def _kweight_sos(sr: int) -> np.ndarray:
    """
    BS.1770 K-weighting filter as a stacked SOS array for scipy.signal.sosfilt.

    Computed directly from the analog prototype parameters in BS.1770 Annex 1
    using the Audio EQ Cookbook (Bristow-Johnson) bilinear transform formulas.
    Does NOT touch pyloudnorm internals — meter._filters internal structure
    has changed between pyloudnorm versions and is not part of its public API.

    Two cascaded stages:
      Stage 1 — Pre-filter (high-shelf, +4 dB at high frequencies):
        Models the free-field to diffuse-field acoustic correction for the
        effect of the listener's head on a mounted microphone. Raises
        high-frequency energy before mean-square integration so that HF
        content is weighted appropriately by the loudness measure.
        Parameters: f0 = 1500 Hz, dBgain = +4 dB, Q = 1/sqrt(2)

      Stage 2 — RLB weighting (2nd-order Butterworth high-pass, ~38.1 Hz):
        Removes DC and very low-frequency energy that carries no perceptual
        loudness but would otherwise inflate the mean-square sum.
        Parameters: fc = 38.13494637408312 Hz

    Returns a (2, 6) ndarray — two SOS sections — ready for sosfilt().
    """
    # Stage 1: Pre-filter high-shelf (Audio EQ Cookbook high-shelf formula)
    f0     = 1500.0              # Hz — shelf midpoint frequency
    dBgain = 4.0                 # dB — shelf gain
    Q      = 1.0 / np.sqrt(2)   # = 0.7071... (Butterworth Q)

    A      = 10.0 ** (dBgain / 40.0)          # sqrt of linear power gain
    w0     = 2.0 * np.pi * f0 / sr            # digital frequency (radians/sample)
    alpha  = np.sin(w0) / (2.0 * Q)           # bandwidth term

    # High-shelf numerator / denominator (Audio EQ Cookbook, high-shelf)
    b0 =       A * ((A + 1) + (A - 1) * np.cos(w0) + 2 * np.sqrt(A) * alpha)
    b1 = -2 * A * ((A - 1) + (A + 1) * np.cos(w0))
    b2 =       A * ((A + 1) + (A - 1) * np.cos(w0) - 2 * np.sqrt(A) * alpha)
    a0 =           ((A + 1) - (A - 1) * np.cos(w0) + 2 * np.sqrt(A) * alpha)
    a1 =  2 *     ((A - 1) - (A + 1) * np.cos(w0))
    a2 =           ((A + 1) - (A - 1) * np.cos(w0) - 2 * np.sqrt(A) * alpha)

    # Second-order section: [b0/a0, b1/a0, b2/a0, 1.0, a1/a0, a2/a0]
    sos_shelf = np.array([[b0/a0, b1/a0, b2/a0, 1.0, a1/a0, a2/a0]])

    # Stage 2: RLB high-pass (2nd-order Butterworth)
    fc     = 38.13494637408312   # Hz
    sos_hp = butter(2, fc / (sr / 2.0), btype='high', output='sos')

    return np.vstack([sos_shelf, sos_hp])   # shape (2, 6)


def _compute_integrated_lufs(y_stereo: np.ndarray, sr: int) -> float:
    """
    ITU-R BS.1770-4 integrated LUFS (K-weighted, gated).

    y_stereo: (2, n_samples) stereo array at sr Hz.
    Returns integrated LUFS as float; returns -70.0 for silent / fully-gated
    signals (pyloudnorm returns -inf in those cases).
    """
    meter  = pyln.Meter(sr)
    data   = y_stereo.T.astype(np.float64)    # pyloudnorm expects (n_samples, n_ch)
    result = meter.integrated_loudness(data)
    return float(result) if np.isfinite(result) else -70.0


def _compute_short_term_lufs_array(y_stereo: np.ndarray, sr: int,
                                    window_s: float = 3.0,
                                    hop_s:    float = 1.0) -> np.ndarray:
    """
    K-weighted short-term loudness in sliding windows (BS.1770, no gating).

    Applies the same K-weighting filters as integrated LUFS but without the
    absolute/relative gating step, matching the BS.1770 short-term definition.
    Uses the BS.1770 stereo channel-sum formula (G_L = G_R = 1.0):
        L_window = -0.691 + 10 × log10(mean_sq_L + mean_sq_R)

    window_s = 3.0, hop_s = 1.0 matches the original librosa.feature.rms()
    block framing that peak_lufs, lufs_variation, and chorus_lift are based on.

    y_stereo: (2, n_samples) stereo array at sr Hz.
    Returns array of per-window LUFS values (float64).
    """
    # Apply K-weighting (BS.1770 pre-filter + RLB high-pass) to each channel.
    # _kweight_sos() returns a single stacked (2, 6) SOS array computed
    # directly from the BS.1770 analog prototype — no pyloudnorm internals used.
    sos  = _kweight_sos(sr)
    data = y_stereo.astype(np.float64).copy()    # (2, n_samples)
    for ch in range(data.shape[0]):
        data[ch] = sosfilt(sos, data[ch])

    n_samples = data.shape[1]
    win_len   = int(window_s * sr)
    hop_len   = max(int(hop_s * sr), 1)

    results = []
    for start in range(0, max(n_samples - win_len + 1, 1), hop_len):
        block  = data[:, start : start + win_len]            # (2, win_len)
        ms_sum = float(np.sum(np.mean(block ** 2, axis=1)))  # sum over L, R
        results.append(-0.691 + 10.0 * np.log10(ms_sum + EPS))

    return np.array(results, dtype=np.float64) if results else np.array([-70.0])


def _true_peak_dbfs(y_stereo: np.ndarray) -> float:
    """
    ITU-R BS.1775-2 true peak level (dBFS), measured via 4× oversampling.

    Standard sample-domain peak measurement misses intersample peaks — peaks
    that occur between samples and can exceed the sample-domain maximum after
    reconstruction. 4× oversampling with a polyphase anti-aliasing filter
    (scipy.signal.resample_poly) captures these peaks reliably.

    Spotify, Apple Music, YouTube, and Tidal all measure true peak at ≥ 4×
    oversampling when deciding whether and by how much to apply limiting during
    normalization. -1.0 dBFS TP is the standard limiting ceiling for most DSPs.

    y_stereo: (2, n_samples) stereo or (n_samples,) mono array.
    Returns peak in dBFS. Typical brick-wall limited master: -0.5 to -0.1 dBFS.
    Silent signal: returns -120.0.
    """
    oversampled = resample_poly(y_stereo, 4, 1)
    peak_linear = float(np.max(np.abs(oversampled)))
    if peak_linear < EPS:
        return -120.0
    return float(20.0 * np.log10(peak_linear))


# ─────────────────────────────────────────────────────────────────────────────
# CANONICAL STREAMING NORMALIZATION
#
# Single source of truth used by both lane_model.py (training) and
# analyze_submissions.py (inference). Having one implementation guarantees
# that training and scoring apply identical gain logic — any future change to
# normalization behaviour only needs to be made in one place.
#
# Previously duplicated as:
#   _norm_eng()                    in lane_model.py       (training)
#   apply_streaming_normalization() in analyze_submissions.py (inference)
# Those two copies were kept manually in sync. They are now replaced by this
# single function, imported by both callers.
# ─────────────────────────────────────────────────────────────────────────────

def apply_streaming_normalization(
    eng: np.ndarray,
    target_lufs: float,
    true_peak_dbfs: "float | None" = None,
) -> "tuple[np.ndarray, float]":
    """
    Simulate streaming-platform loudness normalization on an eng vector.

    Computes the gain the platform would apply to reach target_lufs, then
    shifts the lufs and peak_lufs dimensions by that gain.

    Only attenuates: gain is capped at 0.0 dB (platforms never boost quiet
    masters — Spotify, Apple Music, etc. leave them quieter than the target).

    true_peak_dbfs (BS.1775-2, 4× oversampled):
        If provided and finite, the gain is also capped so that true peak
        never exceeds -1.0 dBFS post-normalization, matching the hardware
        limiter ceiling used by most streaming DSPs. This cap only activates
        for masters with intersample peaks above -1.0 dBFS. When None or
        non-finite, the cap is skipped (equivalent to tp_headroom = 0.0,
        which is harmless because the 0.0 attenuation cap already applies).

    Raises ValueError if "lufs" or "peak_lufs" are absent from ENG_DIMS.
    Callers must treat this as a fatal configuration error — it means
    ENG_DIMS was modified without updating this function.

    Parameters
    ----------
    eng : np.ndarray
        Raw engineering feature vector (32 elements, pre-normalization).
    target_lufs : float
        Streaming loudness target in LUFS (e.g. -14.0 for Spotify).
    true_peak_dbfs : float | None
        BS.1775-2 true peak level of the master in dBFS, or None if
        unavailable (true peak cap is skipped when None).

    Returns
    -------
    eng_normed : np.ndarray
        Copy of eng with lufs and peak_lufs shifted by gain_db.
    gain_db : float
        Gain applied (always <= 0.0).
    """
    try:
        lufs_idx = ENG_DIMS.index("lufs")
        peak_idx = ENG_DIMS.index("peak_lufs")
    except ValueError as exc:
        raise ValueError(
            f"apply_streaming_normalization: required dimension not found in ENG_DIMS: {exc}"
        ) from exc

    # True peak cap: only active when true_peak_dbfs is provided and finite.
    # When true_peak > -1.0 dBFS (intersample peak above limiter ceiling),
    # tp_headroom is negative, further restricting the gain so the normalized
    # true peak lands at -1.0 dBFS. When true_peak <= -1.0 dBFS, tp_headroom
    # is non-negative and the existing 0.0 cap is the binding constraint.
    tp_headroom = (
        (-1.0 - true_peak_dbfs)
        if (true_peak_dbfs is not None and np.isfinite(true_peak_dbfs))
        else 0.0
    )
    gain_db  = min(target_lufs - eng[lufs_idx], tp_headroom, 0.0)
    out      = eng.copy()
    out[lufs_idx] += gain_db
    out[peak_idx] += gain_db
    return out, gain_db


def _compute_drum_bands(drum_mono: np.ndarray, mix_mono: np.ndarray, sr: int,
                        total_rms: float) -> tuple:
    """
    Returns (kick_dom, cymbal_presence, hihat_presence, cymbal_level, hihat_level).

    kick_dom        — dominance ratio (unchanged from v5):
                        RMS(drum_stem, 40–120 Hz) / RMS(mix, 40–120 Hz)
    cymbal_presence — HPSS sustained fraction × band prominence in 4–10 kHz [0,1]
    hihat_presence  — onset density of HPSS percussive 8–14 kHz component [0,1]
    cymbal_level    — RMS(drum_stem, 4–10 kHz) / total_rms [v10 NEW]
    hihat_level     — RMS(drum_stem, 8–14 kHz) / total_rms [v10 NEW]

    WHY HPSS INSTEAD OF DOMINANCE RATIO FOR CYMBAL / HI-HAT (v6):

    The v5 formula drum_cymbal_presence = RMS(drum_8-16kHz) / RMS(mix_8-16kHz)
    conflated hi-hats and crashes and failed in mixed-source tracks. When a bright
    guitar or synth pad adds energy at 8-16 kHz alongside real drum content, the
    mix denominator inflates while the drum numerator is unchanged — the ratio drops
    below the lane mean even with cymbals audibly present throughout the track.

    HPSS (Harmonic-Percussive Source Separation, librosa.decompose.hpss) applies
    a horizontal median filter to the magnitude spectrogram to extract sustained
    tonal content (H) and a vertical median filter for transient broadband content (P):

      H (harmonic/sustained): crash ring-out, ride shimmer, guitar chord sustain.
      P (percussive/transient): hi-hat clicks, kick/snare attacks.

    The H component is not affected by other instruments' level in the full mix —
    it measures what is sustained within the drum stem's own spectrogram. A crash
    ring-out is always captured in H regardless of guitar/synth in the full mix.

    CYMBAL (crash/ride): measured via sustained fraction × band prominence.
      sustained_frac = H_energy(4-10kHz) / total_energy(4-10kHz)
        → high when crash ring-out dominates 4-10kHz drum content
        → lower for snare-only (snare ring has less sustained relative energy in band)
      band_prom = RMS(mag_4-10kHz) / RMS(mag_all_freqs)
        → scales by how much 4-10kHz energy the drum stem actually carries
        → near-silent drum stem (no drums) → band_prom ≈ 0 → product ≈ 0
      Combined: cymbal_presence = (sustained_frac × band_prom) / CYMBAL_SAT

    HI-HAT: measured via onset density of percussive 8-14 kHz signal.
      1. Apply soft percussive mask P/(mag+EPS) to complex STFT D.
      2. Zero all frequency bins outside 8-14 kHz.
      3. Reconstruct time-domain hi-hat signal via iSTFT.
      4. Run librosa onset detection on reconstructed signal.
      5. onset_rate (onsets/sec) / HIHAT_ONSET_SAT → [0,1].

      Guitar chord rings at 8-14 kHz are sustained → H component → soft_p low
      for those bins → soft-masked out. Guitar staccato transients at 8-14 kHz
      are much weaker than real hi-hat transients → below onset_detect threshold.
    """
    # ── Kick dominance (unchanged from v5) ────────────────────────────────────
    kick_mix  = _band_rms_abs(mix_mono, sr, DRUM_KICK_FMIN, DRUM_KICK_FMAX)
    kick_drum = _band_rms_abs(drum_mono, sr, DRUM_KICK_FMIN, DRUM_KICK_FMAX)
    kick_dom  = float(min(kick_drum / (kick_mix + EPS), 1.0)) if kick_mix >= EPS else 0.0

    # ── STFT + HPSS on drum stem ──────────────────────────────────────────────
    n_fft = 2048
    hop   = 512
    D     = librosa.stft(drum_mono, n_fft=n_fft, hop_length=hop)    # complex (1025, n_frames)
    mag   = np.abs(D)                                                  # magnitude spectrogram
    H, P  = librosa.decompose.hpss(mag, kernel_size=31, power=2.0)   # H=sustained, P=transient
    freqs = librosa.fft_frequencies(sr=sr, n_fft=n_fft)              # shape (1025,)

    # ── Cymbal presence: HPSS sustained energy in 4–10 kHz ───────────────────
    cymbal_mask  = (freqs >= CYMBAL_BAND_FMIN) & (freqs <= CYMBAL_BAND_FMAX)
    H_cym        = H[cymbal_mask, :]    # sustained component in cymbal band
    mag_cym      = mag[cymbal_mask, :]  # total magnitude in cymbal band

    H_cym_energy   = float(np.mean(H_cym ** 2))
    band_energy    = float(np.mean(mag_cym ** 2)) + EPS
    sustained_frac = H_cym_energy / band_energy   # fraction of cymbal-band energy that is sustained

    band_rms       = float(np.sqrt(np.mean(mag_cym ** 2)))
    total_mag_rms  = float(np.sqrt(np.mean(mag ** 2))) + EPS
    band_prom      = band_rms / total_mag_rms     # cymbal band's prominence in overall drum stem

    cymbal_presence = float(min((sustained_frac * band_prom) / CYMBAL_SAT, 1.0))

    # ── Hi-hat presence: onset density of percussive 8–14 kHz component ──────
    hihat_mask  = (freqs >= HIHAT_BAND_FMIN) & (freqs <= HIHAT_BAND_FMAX)
    # Soft-mask: amplify percussive content, suppress sustained content per bin.
    # soft_p[f,t] = P[f,t] / (mag[f,t] + EPS) ∈ [0,1] — how percussive each bin is.
    soft_p      = P / (mag + EPS)
    D_hihat     = D * soft_p            # complex STFT weighted toward percussive content
    D_hihat[~hihat_mask, :] = 0.0      # zero non-hihat frequencies
    # iSTFT reconstructs time-domain signal containing only percussive 8-14 kHz content.
    hihat_sig   = librosa.istft(D_hihat, hop_length=hop, length=len(drum_mono))

    oenv        = librosa.onset.onset_strength(y=hihat_sig, sr=sr, hop_length=hop)
    onsets      = librosa.onset.onset_detect(onset_envelope=oenv, sr=sr, hop_length=hop)
    duration    = len(drum_mono) / sr
    onset_rate  = len(onsets) / max(duration, 1.0)
    hihat_presence = float(min(onset_rate / HIHAT_ONSET_SAT, 1.0))

    # ── Cymbal band level: absolute amplitude share in full mix ───────────────
    cymbal_level = float(_band_rms_abs(drum_mono, sr, CYMBAL_BAND_FMIN, CYMBAL_BAND_FMAX)
                         / (total_rms + EPS))

    # ── Hi-hat band level: absolute amplitude share in full mix ──────────────
    hihat_level = float(_band_rms_abs(drum_mono, sr, HIHAT_BAND_FMIN, HIHAT_BAND_FMAX)
                        / (total_rms + EPS))

    return kick_dom, cymbal_presence, hihat_presence, cymbal_level, hihat_level


def _compute_stem_presence(stem_stereo: np.ndarray, sr: int, threshold: float,
                            chunk_s: float = 1.0) -> float:
    """
    Fraction of 1-second chunks in which the stem's RMS exceeds threshold.

    1.0 = stem active throughout the entire track.
    0.5 = stem active in roughly half the track.
    0.0 = stem silent throughout.

    threshold = 0.01 (linear, ≈ -40 dBFS) matches VOCAL_ACTIVITY_THRESHOLD.
    stem_stereo: (2, n_samples) stereo stem array at sr Hz.
    """
    mono      = (stem_stereo[0] + stem_stereo[1]) / 2.0
    chunk_len = max(int(chunk_s * sr), 1)
    n_chunks  = max(len(mono) // chunk_len, 1)
    active    = sum(
        1 for i in range(n_chunks)
        if np.sqrt(np.mean(mono[i * chunk_len : (i + 1) * chunk_len] ** 2)) > threshold
    )
    return float(active / n_chunks)


def _compute_kick_onset_presence(drum_mono: np.ndarray, sr: int) -> float:
    """
    Onset density of percussive events in the kick band (40–120 Hz), normalised
    to DRUM_KICK_ONSET_SAT. Returns a value in [0, 1].

    Measures *rhythmic density* of kick hits — distinct from drum_kick_level
    which measures low-end amplitude. A bass note scores high level, near-zero
    presence. A kick drum scores high on both.
    """
    if np.sqrt(np.mean(drum_mono ** 2)) < EPS:
        return 0.0
    lo = max(DRUM_KICK_FMIN, 1.0)
    hi = min(DRUM_KICK_FMAX, sr / 2.0 - 1.0)
    if lo >= hi:
        return 0.0
    hop        = 512
    sos        = butter(2, [lo, hi], btype='bandpass', fs=sr, output='sos')
    kick_sig   = sosfilt(sos, drum_mono)
    oenv       = librosa.onset.onset_strength(y=kick_sig, sr=sr, hop_length=hop)
    onsets     = librosa.onset.onset_detect(onset_envelope=oenv, sr=sr, hop_length=hop)
    duration   = len(drum_mono) / sr
    onset_rate = len(onsets) / max(duration, 1.0)
    return float(min(onset_rate / DRUM_KICK_ONSET_SAT, 1.0))


def _compute_guitar_dims(guitar_mono: np.ndarray, sr: int) -> tuple:
    """
    Returns (egtr_presence, picking_density, mid_saturation, tonal_sustain).
    All values in [0, 1].

    egtr_presence   — composite [0,1] electric guitar detection score (v9):
                        EGTR_PICK_WEIGHT    × picking_density
                      + EGTR_FLAT_WEIGHT    × mid_saturation
                      + EGTR_SUSTAIN_WEIGHT × tonal_sustain

    picking_density — onset density of the HPSS percussive component in
                      200–5000 Hz, normalized to EGTR_PICK_ONSET_SAT (3.0).
                      Detects pick/strum attack events per second.
                      Electric guitar: 1.5–5+ onsets/sec. Sustained pads: 0–0.5.
                      Piano also produces picking onsets; mid_saturation and
                      tonal_sustain together gate against false-positive piano.

    mid_saturation  — spectral flatness of the guitar stem in 200–5000 Hz,
                      normalized to EGTR_FLATNESS_SAT (0.18).
                      Clean electric guitar: 0.06–0.14 (mid_sat ≈ 0.33–0.78).
                      Semi-clean / light crunch: 0.08–0.18 (mid_sat ≈ 0.44–1.0).
                      Crunchy/overdriven: 0.12–0.28 (mid_sat ≈ 0.67–1.0).
                      Piano: 0.02–0.08 (mid_sat ≈ 0.11–0.44).
                      Acoustic guitar: 0.03–0.10 (mid_sat ≈ 0.17–0.56).
                      Rewards harmonic saturation from magnetic pickups and
                      overdrive — the acoustic signature of electric amplification.

    tonal_sustain   — fraction of mid-band energy in the HPSS H (sustained)
                      component, normalized to EGTR_SUSTAIN_SAT (0.70).
                      Formula:
                        H_mid_energy = mean(H[mid_mask]²)
                        P_mid_energy = mean(P[mid_mask]²)
                        raw_sustain  = H_mid_energy / (H_mid_energy + P_mid_energy + EPS)
                        tonal_sustain = min(raw_sustain / EGTR_SUSTAIN_SAT, 1.0)

                      WHY THIS DIMENSION (v9):
                      picking_density and mid_saturation both measure attack/
                      transient character. A guitar playing sustained power chords,
                      legato lines, or moderate-tempo chord progressions puts most
                      of its energy into the H (sustained) HPSS component — the
                      P (percussive) component only captures the brief pick attack.
                      The sustained note body is invisible to picking_density and
                      adds no flatness signal relative to its RMS contribution.
                      tonal_sustain captures this: as long as guitar notes are
                      ringing, H_mid_energy stays elevated relative to P_mid_energy,
                      regardless of tempo or gain stage.

                      Expected ranges (guitar stem, htdemucs_6s):
                        Electric guitar (any gain, any tempo): H_frac ≈ 0.60–0.85
                        Acoustic guitar:                       H_frac ≈ 0.55–0.75
                        Drum bleed in guitar stem:             H_frac ≈ 0.15–0.35
                        Near-silence (minimal bleed):          ratio noisy but
                                                               egtr_level ≈ 0

    WHY HPSS ON GUITAR STEM (v8, unchanged):
    Even with htdemucs_6s providing a dedicated guitar stem, there is residual
    bleed from other instruments. HPSS separates the percussive (transient
    attack) component from the harmonic (sustained) component within the guitar
    stem itself. Extracting the P component in 200-5000 Hz and counting onsets
    isolates picking/strumming events even with stem bleed, because instrument
    bleed (synths, pads, vocals) tends to be sustained (→ H component) rather
    than transient (→ P component) in this frequency range.

    ACOUSTIC VS ELECTRIC LIMITATION (unchanged from v8):
    Clean acoustic and clean electric guitar score similarly on picking_density.
    mid_saturation partially distinguishes them (electric pickups add resonance
    that raises flatness vs acoustic). tonal_sustain is similar for both.
    The CLAP embedding (already in the feature vector) carries the residual
    discrimination; the SVM learns from CLAP + these dims together.

    Intermediate values (picking_density, mid_saturation, tonal_sustain) are
    returned for debug/tuning but are NOT added to ENG_DIMS — only egtr_presence
    and egtr_level are scoreable dimensions, preserving the atomicity principle.
    """
    if np.sqrt(np.mean(guitar_mono ** 2)) < EPS:
        # Silent guitar stem — no guitar in the track.
        return 0.0, 0.0, 0.0, 0.0

    n_fft    = 2048
    hop      = 512
    D        = librosa.stft(guitar_mono, n_fft=n_fft, hop_length=hop)
    mag      = np.abs(D)
    H, P     = librosa.decompose.hpss(mag, kernel_size=31, power=2.0)
    freqs    = librosa.fft_frequencies(sr=sr, n_fft=n_fft)
    mid_mask = (freqs >= EGTR_MID_FMIN) & (freqs <= EGTR_MID_FMAX)

    # ── Picking density: onset rate of percussive mid-band content ────────────
    # Soft-mask the complex STFT toward percussive content, isolate mid band.
    soft_p     = P / (mag + EPS)        # per-bin percussiveness score [0,1]
    D_pick     = D * soft_p             # complex STFT weighted toward transients
    D_pick[~mid_mask, :] = 0.0         # zero bins outside guitar mid band
    pick_sig   = librosa.istft(D_pick, hop_length=hop, length=len(guitar_mono))
    oenv       = librosa.onset.onset_strength(y=pick_sig, sr=sr, hop_length=hop)
    onsets     = librosa.onset.onset_detect(onset_envelope=oenv, sr=sr, hop_length=hop)
    duration   = len(guitar_mono) / sr
    onset_rate = len(onsets) / max(duration, 1.0)
    picking_density = float(min(onset_rate / EGTR_PICK_ONSET_SAT, 1.0))

    # ── Mid-band spectral flatness: harmonic saturation from pickups/overdrive ─
    # Flatness = geometric_mean / arithmetic_mean of magnitude spectrum.
    # Ranges [0, 1]: 0 = perfectly tonal (sine-like), 1 = flat noise.
    # Electric guitar pickups and overdrive add intermodulation products that
    # raise flatness compared to clean acoustic guitar or piano.
    mag_mid = mag[mid_mask, :]
    if mag_mid.size == 0 or mag_mid.shape[0] == 0:
        mid_saturation = 0.0
    else:
        geo_mean    = np.exp(np.mean(np.log(mag_mid + EPS), axis=0))  # per-frame geo mean
        arith_mean  = np.mean(mag_mid, axis=0) + EPS                  # per-frame arith mean
        frame_flat  = geo_mean / arith_mean                            # flatness [0,1] per frame
        raw_flatness = float(np.mean(frame_flat))
        mid_saturation = float(min(raw_flatness / EGTR_FLATNESS_SAT, 1.0))

    # ── Tonal sustain: fraction of mid-band energy in sustained H component ───
    # Captures guitar notes that are ringing/sustaining — the body of each note
    # after the pick attack, which is entirely invisible to picking_density.
    # H (harmonic/sustained) component of HPSS in 200–5000 Hz.
    # P (percussive/transient) component in same band (pick attacks only).
    # raw_sustain = H_energy / (H_energy + P_energy) ∈ [0, 1].
    # When guitar notes sustain: H >> P, raw_sustain high → tonal_sustain high.
    # When only drum transients bleed in: P >> H, raw_sustain low.
    H_mid        = H[mid_mask, :]
    P_mid        = P[mid_mask, :]
    H_mid_energy = float(np.mean(H_mid ** 2))
    P_mid_energy = float(np.mean(P_mid ** 2))
    raw_sustain  = H_mid_energy / (H_mid_energy + P_mid_energy + EPS)
    tonal_sustain = float(min(raw_sustain / EGTR_SUSTAIN_SAT, 1.0))

    egtr_presence = float(
        EGTR_PICK_WEIGHT    * picking_density +
        EGTR_FLAT_WEIGHT    * mid_saturation  +
        EGTR_SUSTAIN_WEIGHT * tonal_sustain
    )

    return egtr_presence, picking_density, mid_saturation, tonal_sustain


def extract_mir_features(path):
    """
    Returns (l3_vec, eng_vec, true_peak_dbfs).

    eng_vec is a 38-element array whose indices correspond to ENG_DIMS.
    Indices 0–25 are unchanged from cache version 3.
    Indices 26–29 are drum sub-dimensions (v4–v6, unchanged in v8):
      26  drum_kick_level      — kick dominance ratio (v5, renamed v10)
      27  drum_cymbal_presence — HPSS sustained fraction × band prom in 4–10 kHz (v6)
      28  drum_confidence      — 0.55×kick + 0.30×hihat + 0.15×cymbal (v6)
      29  hihat_presence       — onset density of percussive 8–14 kHz (v6 NEW)
    Indices 30–31 are electric guitar dimensions (v8 NEW):
      30  egtr_presence        — composite guitar detection score [0,1]
      31  egtr_level           — egtr_presence × (guitar_rms / total_rms)
    Indices 32–37 are stem presence + level dimensions (v10 NEW):
      32  drum_presence        — fraction of 1-sec chunks with active drum signal
      33  other_presence       — fraction of 1-sec chunks with active other signal
      34  bass_presence        — fraction of 1-sec chunks with active bass signal
      35  drum_kick_presence   — kick onset density in 40–120 Hz / DRUM_KICK_ONSET_SAT
      36  drum_cymbal_level    — RMS(drum, 4–10 kHz) / total_rms
      37  hihat_level          — RMS(drum, 8–14 kHz) / total_rms

    drum_punch (index 2) is gated by drum_confidence:
      drum_punch = crest_factor(drum_stem) × drum_confidence

    v8 changes vs v7:
      htdemucs_6s: model switched from htdemucs (4-stem) to htdemucs_6s (6-stem).
        New stems: guitar (_STEM_GUITAR=4), piano (_STEM_PIANO=5).
        other_prominence (idx 12) now reflects only synths/pads/strings —
        guitar and piano are no longer blended into the other stem.
      total_rms: now includes guitar_rms + piano_rms (previously only 4 stems).
        This makes all stem-share ratios (drum_prominence, bass_prominence,
        other_prominence, vocal_level) slightly different from v7.
      egtr_presence (idx 30): HPSS-based picking density + mid-band spectral
        flatness on the dedicated guitar stem. See _compute_guitar_dims().
      egtr_level (idx 31): egtr_presence × (guitar_rms / total_rms).
        Scales presence by the guitar stem's share of the total mix RMS.
        A guitar that is present but buried low in the mix scores lower than
        one that is both detectable and prominent.
      eng_vec now has 32 elements. All v7 cache entries are incompatible.

    true_peak_dbfs is the ITU-R BS.1775-2 true peak level of the full mix
    (4× oversampled, in dBFS). Returned so the caller can enforce the -1.0 dBFS
    true peak ceiling in apply_streaming_normalization() correctly.

    Expected outcomes by source type (approximate, v9 constants):
      Distorted electric guitar prominent in mix:
        egtr_presence ≈ 0.75–1.00, egtr_level ≈ 0.20–0.50
      Clean/semi-clean electric guitar throughout:
        egtr_presence ≈ 0.55–0.85, egtr_level ≈ 0.10–0.35
      Acoustic guitar only (no electric):
        egtr_presence ≈ 0.40–0.70, egtr_level ≈ 0.05–0.20
      Piano-led track (no guitar):
        egtr_presence ≈ 0.10–0.35, egtr_level ≈ 0.01–0.10
      Pad/synth-led track (no guitar):
        egtr_presence ≈ 0.05–0.25, egtr_level ≈ 0.00–0.08
    """
    y, sr        = librosa.load(path, sr=44100, mono=True)
    y_stereo, _  = librosa.load(path, sr=sr, mono=False)
    if y_stereo.ndim == 1:                               # guard: mono input
        y_stereo = np.stack([y_stereo, y_stereo])
    emb_l3, _    = openl3.get_audio_embedding(
        y, sr, content_type="music", embedding_size=512, verbose=False
    )
    l3_vec = np.mean(emb_l3, axis=0)
    src    = get_stems(path)   # shape: (6, 2, n_samples) — stereo stems at 44100 Hz

    # ── Stem RMS values ───────────────────────────────────────────────────────
    drum_rms   = np.sqrt(np.mean(src[_STEM_DRUMS] **2))
    bass_rms   = np.sqrt(np.mean(src[_STEM_BASS]  **2))
    other_rms  = np.sqrt(np.mean(src[_STEM_OTHER] **2))
    voc_rms    = np.sqrt(np.mean(src[_STEM_VOCALS]**2))
    guitar_rms = np.sqrt(np.mean(src[_STEM_GUITAR]**2))
    piano_rms  = np.sqrt(np.mean(src[_STEM_PIANO] **2))
    # v8: total_rms includes all 6 stems for accurate mix-share ratios.
    total_rms  = drum_rms + bass_rms + other_rms + voc_rms + guitar_rms + piano_rms + EPS

    # ── Drum crest factor (raw, before gating) ────────────────────────────────
    drum_crest_raw = np.max(np.abs(src[_STEM_DRUMS])) / (drum_rms + EPS)

    # ── Drum sub-band analysis (v6, unchanged in v8) ──────────────────────────
    drum_mono = (src[_STEM_DRUMS][0] + src[_STEM_DRUMS][1]) / 2.0

    drum_kick_lvl, drum_cymbal_pres, hihat_pres, drum_cymbal_lvl, hihat_lvl = \
        _compute_drum_bands(drum_mono, y, sr, total_rms)

    # drum_confidence: weighted sum of three [0,1] drum indicators.
    kick_conf    = min(drum_kick_lvl / DRUM_KICK_SAT, 1.0)
    hihat_conf   = hihat_pres           # already in [0, 1]
    cymbal_conf  = drum_cymbal_pres     # already in [0, 1]
    drum_conf    = float(
        DRUM_KICK_WEIGHT   * kick_conf +
        DRUM_HIHAT_WEIGHT  * hihat_conf +
        DRUM_CYMBAL_WEIGHT * cymbal_conf
    )

    # Gate drum_punch by confidence: suppresses false positives from guitar bleed.
    drum_punch_val = drum_crest_raw * drum_conf

    # ── Bass crest factor ─────────────────────────────────────────────────────
    bass_tightness_val = np.max(np.abs(src[_STEM_BASS])) / (bass_rms + EPS)

    # ── Mix-level ratios (stem RMS / total RMS) ───────────────────────────────
    drum_lvl   = drum_rms  / total_rms
    bass_lvl   = bass_rms  / total_rms
    other_lvl  = other_rms / total_rms
    voc_level  = voc_rms   / total_rms

    # ── Vocal sub-dimensions ──────────────────────────────────────────────────
    voc_mono = (src[_STEM_VOCALS][0] + src[_STEM_VOCALS][1]) / 2.0
    voc_rms_frames = librosa.feature.rms(
        y=voc_mono, frame_length=2048, hop_length=512
    )[0]
    voc_presence   = float(np.mean(voc_rms_frames > VOCAL_ACTIVITY_THRESHOLD))
    voc_brightness = float(np.mean(
        librosa.feature.spectral_centroid(y=voc_mono, sr=sr)
    ))

    # ── Electric guitar sub-dimensions (v9 updated) ──────────────────────────
    guitar_mono = (src[_STEM_GUITAR][0] + src[_STEM_GUITAR][1]) / 2.0
    egtr_pres, _picking, _flatness, _sustain = _compute_guitar_dims(guitar_mono, sr)
    egtr_level_val = float(egtr_pres * (guitar_rms / total_rms))

    # ── Stem temporal presence (v10 NEW) ─────────────────────────────────────
    drum_presence_val  = _compute_stem_presence(
        src[_STEM_DRUMS],  sr, DRUM_ACTIVITY_THRESHOLD
    )
    other_presence_val = _compute_stem_presence(
        src[_STEM_OTHER],  sr, OTHER_ACTIVITY_THRESHOLD
    )
    bass_presence_val  = _compute_stem_presence(
        src[_STEM_BASS],   sr, BASS_ACTIVITY_THRESHOLD
    )

    # ── Kick onset density (v10 NEW) ─────────────────────────────────────────
    drum_kick_pres_val = _compute_kick_onset_presence(drum_mono, sr)

    # ── Mid-side stereo width ─────────────────────────────────────────────────
    mid, side = (y_stereo[0] + y_stereo[1]) / 2, (y_stereo[0] - y_stereo[1]) / 2
    ms_ratio  = (np.sqrt(np.mean(side**2)) /
                 (np.sqrt(np.mean(mid**2)) + np.sqrt(np.mean(side**2)) + EPS))

    # ── Short-term loudness (K-weighted, BS.1770 — no gating) ────────────────
    st_lufs = _compute_short_term_lufs_array(y_stereo, sr, window_s=3.0, hop_s=1.0)

    # ── Clipping ─────────────────────────────────────────────────────────────
    c_events, c_total, c_ratio, c_max, c_mean = _measure_clipping(y)

    # ── Assemble eng_vec — must match ENG_DIMS order exactly ─────────────────
    eng_vec = np.array([
        _compute_integrated_lufs(y_stereo, sr),                  # 0  lufs (ITU-R BS.1770-4 — K-weighted, gated)
        np.max(np.abs(y)) / (np.sqrt(np.mean(y**2)) + EPS),     # 1  crest_factor
        drum_punch_val,                                          # 2  drum_punch (gated)
        drum_lvl,                                                # 3  drum_level (v10: renamed from drum_prominence)
        ms_ratio,                                                # 4  stereo_width
        librosa.feature.spectral_centroid(y=y, sr=sr).mean(),   # 5  brightness
        entropy(np.mean(librosa.feature.chroma_cqt(y=y, sr=sr), axis=1)),  # 6  harmonic_complexity
        c_events,                                                # 7  clip_num_events
        c_total,                                                 # 8  clip_total
        c_ratio,                                                 # 9  clip_ratio
        c_max,                                                   # 10 clip_max_run
        c_mean,                                                  # 11 clip_mean_run
        other_lvl,                                               # 12 other_level (v10: renamed from other_prominence)
        voc_level,                                               # 13 vocal_level
        voc_presence,                                            # 14 vocal_presence
        voc_brightness,                                          # 15 vocal_brightness
        np.max(st_lufs) - np.median(st_lufs),                   # 16 chorus_lift
        np.mean(librosa.feature.zero_crossing_rate(y)),         # 17 crispness
        np.mean(librosa.feature.spectral_flatness(y=y)),        # 18 clarity
        len(y) / sr,                                            # 19 length_s
        np.std(st_lufs),                                        # 20 lufs_variation
        np.max(st_lufs),                                        # 21 peak_lufs
        np.mean(librosa.feature.spectral_rolloff(y=y, sr=sr)), # 22 rolloff
        bass_tightness_val,                                     # 23 bass_tightness
        bass_lvl,                                               # 24 bass_level (v10: renamed from bass_prominence)
        np.nan_to_num(np.corrcoef(y_stereo[0], y_stereo[1])[0, 1]),  # 25 phase_corr
        drum_kick_lvl,                                          # 26 drum_kick_level (v10: renamed from drum_kick_presence)
        drum_cymbal_pres,                                       # 27 drum_cymbal_presence (v6)
        drum_conf,                                              # 28 drum_confidence (v6)
        hihat_pres,                                             # 29 hihat_presence (v6)
        egtr_pres,                                              # 30 egtr_presence (v8 NEW)
        egtr_level_val,                                         # 31 egtr_level (v8 NEW)
        drum_presence_val,                                      # 32 drum_presence (v10 NEW)
        other_presence_val,                                     # 33 other_presence (v10 NEW)
        bass_presence_val,                                      # 34 bass_presence (v10 NEW)
        drum_kick_pres_val,                                     # 35 drum_kick_presence (v10 NEW)
        drum_cymbal_lvl,                                        # 36 drum_cymbal_level (v10 NEW)
        hihat_lvl,                                              # 37 hihat_level (v10 NEW)
    ])

    # ── True peak (BS.1775-2, 4× oversampled) ────────────────────────────────
    # Returned as third element so the modal worker / analyze_submissions.py
    # can enforce the -1.0 dBFS true peak ceiling in normalization simulation.
    true_peak = _true_peak_dbfs(y_stereo)

    return l3_vec, eng_vec, true_peak