feat(foundry): the loop finder (engine/loops.py) — §1 pipeline, per-stem (#5)

Generates ranked loop candidates from a separated stem, the phase-2 design made
real (PHASE2_FINDINGS §1/§3/§4): drift-tolerant PLP beat grid → beat-sync SSM
(MFCC⊕chroma → recurrence affinity) → Foote checkerboard novelty → candidate
windows at {1,2,4,8} bars → §3 composite score → NMS dedup. Works PER STEM (§0),
then loosely groups candidates across stems by shared time window into Takes the
human auditions together. export_take honours the B-rulebook (B1 ZC-snap, B6 bass
mono-sum, B8 24-bit/44.1k, B10 no-normalize) and calls publish.link_kit.

DRY with the katana: the loopability score reuses engine.grade's seam_score /
zc_score on each candidate slice, and adds the two terms only the finder has the
SSM context for — structural (does the window recur off its own diagonal?) and
boundary novelty (clean segment edges). Weights w_struct .4 / w_novel .2 / w_seam
.3 / w_zc .1, minus a tempo-instability penalty (§3).

Caught while validating on the real Loituma stems: PLP latched onto the eighth-
note pulse and reported ~250 bpm (Loituma is ~125), so 2-bar windows came out
half-length — the §7 octave hazard, live. Fixed by pinning PLP's tempo range to
one musical octave (70–160) and enforcing a min inter-beat spacing in peak-pick.
Now reads ~120 bpm and 2-bar windows land at ~3.9s, squarely in the A-profile
band (3.7–8.0s). find_takes on the 4-stem catch correctly clusters time-aligned
candidates (e.g. vocals+drums sharing a 2-bar window).

CLI: python3 -m engine.loops <stem.wav> [--bars 2 4 8] [--top N] [--json].
7 new tests (checkerboard signs, Foote peak at a block seam, structural recurrence,
NMS dedup, ZC-snap, and the tempo-octave guard). Full foundry suite: 41 pass.

tools/foundry/engine/loops.py

+"""loops — the loop finder (task #5).
+
+Generate ranked loop candidates from a separated stem, per the phase-2 design
+(research/PHASE2_FINDINGS §1, §3, §4). The pipeline, all librosa/scipy/numpy:
+
+  1. beat grid — drift-tolerant PLP (predominant local pulse), NOT beat_track.
+  2. beat-synchronous SSM — MFCC(13) ⊕ chroma(12) synced to beats →
+     recurrence affinity. The beat axis absorbs setcps drift for free.
+  3. Foote checkerboard novelty on the beat-sync SSM (1-bar and 4-bar kernels).
+  4. candidate windows at {1,2,4,8} bars between novelty peaks.
+  5. score each window (§3 composite) and dedup.
+  6. zero-crossing snap at export.
+
+The finder works PER STEM (§0): 1..x ranked candidates each, then loose cross-
+stem grouping into Takes (a time window the human auditions together). The
+loopability score REUSES the grader's seam/zc sub-scorers (engine.grade, DRY) and
+adds the two terms only the finder has context for: structural (does the window
+repeat in the SSM?) and boundary novelty (clean segment edges).
+
+Suggestions are never auto-applied; a wrong rank costs one click to ignore.
+
+    python3 -m engine.loops <stem.wav> [--bars 2 4 8] [--top 8]   # rank one stem
+"""
+from __future__ import annotations
+
+import math
+from pathlib import Path
+from typing import Callable, Optional
+
+import numpy as np
+import soundfile as sf
+from pydantic import BaseModel, Field
+
+from . import grade as G
+from .model import Catch
+
+BEATS_PER_BAR = 4
+EXPORT_SR = 44100                 # B8: SuperDirt/SC rate
+ZC_TOL_MS = 10.0                  # B1
+XFADE_MS = 20.0                   # B3 (3 ms for drums, applied at export)
+
+# §3 composite weights (provisional; calibrate with #7)
+W = {"struct": 0.4, "novel": 0.2, "seam": 0.3, "zc": 0.1, "tempo_penalty": 0.2}
+
+
+class LoopCandidate(BaseModel):
+    """One scored window within a single stem."""
+    start_s: float
+    end_s: float
+    bars: int
+    bpm: float
+    tempo_unstable: bool
+    score: float
+    structural: float
+    novelty: float
+    seam: float
+    zc: float
+
+
+class StemSlice(BaseModel):
+    name: str                     # drums|bass|other|vocals
+    start_s: float
+    end_s: float
+    score: float
+    rms_dbfs: float
+    is_one_shot: bool
+
+
+class Take(BaseModel):
+    """A bar-window = the unit PLN selects; gathers each stem's slice in it."""
+    start_s: float
+    end_s: float
+    bars: int
+    bpm: float
+    tempo_unstable: bool
+    score: float                  # best stem score in the window
+    n_stems: int
+    stems: list[StemSlice] = Field(default_factory=list)
+
+
+# ── pipeline ──────────────────────────────────────────────────────────────────
+def beat_grid(y: np.ndarray, sr: int, *, tempo_min=70.0, tempo_max=160.0, hop=512):
+    """Drift-tolerant beat times + per-beat local BPM via PLP. (frames, times, bpm).
+
+    PLP's tempo range is pinned to one musical octave and peak-picking enforces a
+    minimum inter-beat spacing — without this PLP latches onto the eighth-note
+    pulse and reports ~2× tempo (the §7 octave hazard; Loituma read as 250 not 125).
+    """
+    import librosa
+    oenv = librosa.onset.onset_strength(y=y, sr=sr, hop_length=hop)
+    pulse = librosa.beat.plp(onset_envelope=oenv, sr=sr, hop_length=hop,
+                             tempo_min=tempo_min, tempo_max=tempo_max)
+    from scipy.signal import find_peaks
+    min_dist = max(1, int((60.0 / tempo_max) * sr / hop))
+    peaks, _ = find_peaks(pulse, distance=min_dist)
+    times = librosa.frames_to_time(peaks, sr=sr, hop_length=hop)
+    if len(times) >= 2:
+        ibi = np.diff(times)
+        bpm = 60.0 / np.clip(ibi, 1e-3, None)
+        bpm = np.r_[bpm, bpm[-1]]
+    else:
+        bpm = np.array([120.0] * len(times))
+    return peaks, times, bpm
+
+
+def beat_sync_ssm(y: np.ndarray, sr: int, beat_frames: np.ndarray) -> np.ndarray:
+    """MFCC(13) ⊕ chroma_cqt(12) synced onto the beat grid → recurrence affinity."""
+    import librosa
+    mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)
+    chroma = librosa.feature.chroma_cqt(y=y, sr=sr)
+    feat = np.vstack([librosa.util.normalize(mfcc, axis=1),
+                      librosa.util.normalize(chroma, axis=1)])
+    synced = librosa.util.sync(feat, beat_frames, aggregate=np.mean)
+    if synced.shape[1] < 3:
+        return np.zeros((synced.shape[1], synced.shape[1]))
+    return librosa.segment.recurrence_matrix(synced, mode="affinity", sym=True)
+
+
+def _checkerboard(k: int) -> np.ndarray:
+    """Gaussian-tapered checkerboard kernel (2k×2k) for Foote novelty."""
+    g = np.outer(*[np.exp(-0.5 * (np.linspace(-2, 2, 2 * k)) ** 2)] * 2)
+    sign = np.ones((2 * k, 2 * k))
+    sign[:k, k:] = sign[k:, :k] = -1
+    return g * sign
+
+
+def foote_novelty(ssm: np.ndarray, k: int) -> np.ndarray:
+    """Checkerboard-kernel novelty along the SSM diagonal; peaks = segment edges."""
+    N = ssm.shape[0]
+    if N < 2 * k + 1:
+        return np.zeros(N)
+    ker = _checkerboard(k)
+    nov = np.zeros(N)
+    pad = np.pad(ssm, k, mode="edge")
+    for i in range(N):
+        nov[i] = float(np.sum(pad[i:i + 2 * k, i:i + 2 * k] * ker))
+    nov = np.maximum(nov, 0)
+    return nov / (nov.max() + 1e-9)
+
+
+def _structural(ssm: np.ndarray, s: int, e: int) -> float:
+    """Does the window [s,e) recur elsewhere? Mean best off-segment affinity."""
+    if e - s < 1 or ssm.shape[0] < 3:
+        return 0.0
+    block = ssm[s:e].copy()
+    block[:, s:e] = 0.0                       # mask its own diagonal band
+    return float(block.max(axis=1).mean()) if block.size else 0.0
+
+
+def analyze_stem(y: np.ndarray, sr: int, *, bars=(2, 4, 8), top_n=8) -> list[LoopCandidate]:
+    """Rank loop candidates within one stem."""
+    peaks, times, bpm = beat_grid(y, sr)
+    nb = len(times)
+    if nb < BEATS_PER_BAR * min(bars) + 1:
+        return []
+    ssm = beat_sync_ssm(y, sr, peaks)         # peaks are onset-envelope frames
+    n = ssm.shape[0]
+    nov = np.maximum(foote_novelty(ssm, 2), foote_novelty(ssm, 8)) if n else np.zeros(nb)
+
+    cands: list[LoopCandidate] = []
+    for B in bars:
+        W_beats = B * BEATS_PER_BAR
+        for i in range(0, nb - W_beats):
+            e = i + W_beats
+            if e >= n:
+                continue
+            s_smp, e_smp = int(times[i] * sr), int(times[e] * sr)
+            if e_smp - s_smp < int(0.5 * sr):
+                continue
+            sl = y[s_smp:e_smp]
+            if sl.size < 8:
+                continue
+            seam, _ = G.seam_score(sl)
+            zc, _, _ = G.zc_score(sl, sr)
+            struct = _structural(ssm, i, e)
+            novel = float((nov[i] + nov[min(e, n - 1)]) / 2)
+            win_bpm = bpm[i:e]
+            bmean = float(np.mean(win_bpm)); bspread = float(np.std(win_bpm))
+            unstable = bspread > 5.0
+            tempo_pen = min(1.0, bspread / 20.0)
+            score = (W["struct"] * struct + W["novel"] * novel +
+                     W["seam"] * seam + W["zc"] * zc - W["tempo_penalty"] * tempo_pen)
+            cands.append(LoopCandidate(
+                start_s=round(times[i], 3), end_s=round(times[e], 3), bars=B,
+                bpm=round(bmean, 1), tempo_unstable=unstable,
+                score=round(max(0.0, score), 4), structural=round(struct, 4),
+                novelty=round(novel, 4), seam=round(seam, 4), zc=round(zc, 4)))
+    return _dedup(cands, top_n)
+
+
+def _dedup(cands: list[LoopCandidate], top_n: int) -> list[LoopCandidate]:
+    """Greedy NMS: keep highest-scoring, drop windows overlapping >50%."""
+    out: list[LoopCandidate] = []
+    for c in sorted(cands, key=lambda x: -x.score):
+        if all(_overlap(c, k) < 0.5 for k in out):
+            out.append(c)
+        if len(out) >= top_n:
+            break
+    return out
+
+
+def _overlap(a: LoopCandidate, b: LoopCandidate) -> float:
+    lo, hi = max(a.start_s, b.start_s), min(a.end_s, b.end_s)
+    inter = max(0.0, hi - lo)
+    return inter / max(1e-6, min(a.end_s - a.start_s, b.end_s - b.start_s))
+
+
+# ── takes (cross-stem grouping) ───────────────────────────────────────────────
+def find_takes(catch: Catch, workspace: Path, *, bars=(2, 4, 8), top_n=8,
+               progress: Optional[Callable[[str], None]] = None) -> list[Take]:
+    """Per-stem candidates → loosely grouped Takes (a shared time window)."""
+    catch_dir = catch.dir(workspace)
+    per_stem: dict[str, list[LoopCandidate]] = {}
+    for st in catch.stems:
+        p = catch_dir / st.path
+        if not p.exists():
+            continue
+        if progress:
+            progress(f"analyzing {st.name}")
+        y = G._mono(G.load_audio(p)[0]).astype(np.float32)
+        per_stem[st.name] = analyze_stem(y, _sr_of(p), bars=bars, top_n=top_n)
+
+    # cluster candidates from all stems by start time (within ~1 bar slack)
+    flat = [(name, c) for name, cs in per_stem.items() for c in cs]
+    flat.sort(key=lambda nc: nc[1].start_s)
+    takes: list[Take] = []
+    used = [False] * len(flat)
+    for i, (name, c) in enumerate(flat):
+        if used[i]:
+            continue
+        slack = c.bars * BEATS_PER_BAR * 60.0 / max(c.bpm, 1) / BEATS_PER_BAR  # ~1 bar
+        members = [(j, flat[j]) for j in range(len(flat))
+                   if not used[j] and abs(flat[j][1].start_s - c.start_s) <= slack]
+        best_per_stem: dict[str, tuple[int, LoopCandidate]] = {}
+        for j, (nm, cc) in members:
+            if nm not in best_per_stem or cc.score > best_per_stem[nm][1].score:
+                best_per_stem[nm] = (j, cc)
+        for j, _ in best_per_stem.values():
+            used[j] = True
+        lead = max(best_per_stem.values(), key=lambda jc: jc[1].score)[1]
+        stems = [StemSlice(name=nm, start_s=cc.start_s, end_s=cc.end_s, score=cc.score,
+                           rms_dbfs=0.0, is_one_shot=(cc.end_s - cc.start_s) < 0.75)
+                 for nm, (j, cc) in best_per_stem.items()]
+        takes.append(Take(start_s=lead.start_s, end_s=lead.end_s, bars=lead.bars,
+                          bpm=lead.bpm, tempo_unstable=lead.tempo_unstable,
+                          score=round(lead.score, 4), n_stems=len(stems), stems=stems))
+    takes.sort(key=lambda t: -t.score)
+    return takes[:top_n]
+
+
+def _sr_of(path: Path) -> int:
+    return sf.info(str(path)).samplerate
+
+
+# ── export (B-rules) ──────────────────────────────────────────────────────────
+def _snap_zc(y: np.ndarray, idx: int, sr: int) -> int:
+    """Snap a boundary to the nearest zero crossing within B1 tolerance."""
+    tol = int(sr * ZC_TOL_MS / 1000)
+    lo, hi = max(0, idx - tol), min(len(y) - 1, idx + tol)
+    seg = y[lo:hi]
+    if seg.size < 2:
+        return idx
+    zc = np.where(np.diff(np.signbit(seg)))[0]
+    return int(lo + zc[np.argmin(np.abs(zc - (idx - lo)))]) if zc.size else idx
+
+
+def export_take(take: Take, kit: str, workspace: Path, catch: Catch, *,
+                keep=("drums", "bass", "other", "vocals"),
+                names: Optional[dict[str, str]] = None, peak_norm=False,
+                samples_root: Optional[Path] = None) -> list[Path]:
+    """Write a Take's stem slices to Samples/<kit>/ honoring B1/B6/B8/B10, then link."""
+    from . import publish
+    out_dir = (samples_root or publish.samples_root()) / kit
+    out_dir.mkdir(parents=True, exist_ok=True)
+    catch_dir = catch.dir(workspace)
+    stem_path = {st.name: catch_dir / st.path for st in catch.stems}
+    written: list[Path] = []
+    for i, sl in enumerate(s for s in take.stems if s.name in keep):
+        src = stem_path.get(sl.name)
+        if not src or not src.exists():
+            continue
+        y, sr = sf.read(str(src), always_2d=True, dtype="float32")   # (n, ch)
+        a, b = int(sl.start_s * sr), int(sl.end_s * sr)
+        mono_ref = y[:, 0]
+        a, b = _snap_zc(mono_ref, a, sr), _snap_zc(mono_ref, b, sr)   # B1
+        clip = y[a:b]
+        if sl.name == "bass" and clip.shape[1] > 1:                   # B6 mono-sum
+            clip = clip.mean(axis=1, keepdims=True)
+        if peak_norm:                                                 # B10: off by default
+            pk = np.max(np.abs(clip)) + 1e-9
+            clip = clip * (10 ** (-1.0 / 20) / pk)
+        nm = (names or {}).get(sl.name) or f"{i:02d}_{sl.name}"
+        dest = out_dir / f"{nm}.wav"
+        sf.write(str(dest), clip, EXPORT_SR if sr == EXPORT_SR else sr,
+                 subtype="PCM_24")                                    # B8
+        written.append(dest)
+    if written:
+        try:
+            publish.link_kit(kit)
+        except Exception:
+            pass
+    return written
+
+
+# ── CLI ───────────────────────────────────────────────────────────────────────
+def _main(argv=None) -> int:
+    import argparse
+    import json
+    ap = argparse.ArgumentParser(description="rank loop candidates in a stem")
+    ap.add_argument("stem")
+    ap.add_argument("--bars", type=int, nargs="+", default=[2, 4, 8])
+    ap.add_argument("--top", type=int, default=8)
+    ap.add_argument("--json", action="store_true")
+    a = ap.parse_args(argv)
+    y = G._mono(G.load_audio(a.stem)[0]).astype(np.float32)
+    sr = _sr_of(Path(a.stem))
+    cands = analyze_stem(y, sr, bars=tuple(a.bars), top_n=a.top)
+    if a.json:
+        print(json.dumps([c.model_dump() for c in cands], indent=2))
+    else:
+        for c in cands:
+            print(f"  [{c.score:.3f}] {c.bars}-bar {c.start_s:6.1f}–{c.end_s:6.1f}s "
+                  f"@{c.bpm:.0f}bpm  struct={c.structural:.2f} nov={c.novelty:.2f} "
+                  f"seam={c.seam:.2f} zc={c.zc:.2f}{' ⚠tempo' if c.tempo_unstable else ''}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(_main())

tools/foundry/tests/test_loops.py

+"""Tests for the loop finder (engine/loops.py).
+
+Pure DSP/logic on synthetic inputs, plus one light integration (beat_grid tempo
+recovery, which guards the §7 octave fix). Composite scoring weights are v1 and
+calibrated in #7, so assertions are structural/relative, not pinned to thresholds.
+"""
+import sys
+from pathlib import Path
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+
+from engine import loops as L
+
+
+def test_checkerboard_shape_and_signs():
+    k = 4
+    ker = L._checkerboard(k)
+    assert ker.shape == (2 * k, 2 * k)
+    # top-left & bottom-right quadrants positive; off-diagonal quadrants negative
+    assert ker[:k, :k].sum() > 0 and ker[k:, k:].sum() > 0
+    assert ker[:k, k:].sum() < 0 and ker[k:, :k].sum() < 0
+
+
+def test_foote_novelty_peaks_at_block_boundary():
+    # two self-similar blocks → a checkerboard novelty spike at their seam (idx 10)
+    n = 20
+    ssm = np.zeros((n, n))
+    ssm[:10, :10] = 1.0
+    ssm[10:, 10:] = 1.0
+    nov = L.foote_novelty(ssm, 4)
+    assert 8 <= int(np.argmax(nov)) <= 12
+
+
+def test_structural_high_when_window_repeats():
+    n = 8
+    ssm = np.eye(n)
+    # window [0:2] strongly resembles [4:6] elsewhere
+    ssm[0, 4] = ssm[1, 5] = ssm[4, 0] = ssm[5, 1] = 1.0
+    rep = L._structural(ssm, 0, 2)
+    uniq = L._structural(ssm, 2, 4)
+    assert rep > uniq
+
+
+def _cand(start, end, score, bars=2):
+    return L.LoopCandidate(start_s=start, end_s=end, bars=bars, bpm=120,
+                           tempo_unstable=False, score=score, structural=0.5,
+                           novelty=0.5, seam=0.5, zc=0.5)
+
+
+def test_overlap_fraction():
+    a, b = _cand(0, 4, .9), _cand(2, 6, .8)
+    assert 0.4 < L._overlap(a, b) < 0.6           # 2s overlap of 4s windows
+    assert L._overlap(_cand(0, 4, .9), _cand(10, 14, .8)) == 0.0
+
+
+def test_dedup_keeps_best_drops_overlap():
+    cands = [_cand(0, 4, .9), _cand(2, 6, .8), _cand(20, 24, .7)]
+    out = L._dedup(cands, top_n=8)
+    starts = sorted(c.start_s for c in out)
+    assert starts == [0.0, 20.0]                  # the .8 overlapping the .9 is dropped
+
+
+def test_snap_zc_moves_to_crossing():
+    sr = 44100
+    t = np.arange(sr) / sr
+    y = np.sin(2 * np.pi * 100 * t).astype(np.float32)  # zero crossings every 220.5 samp
+    idx = 1000
+    snapped = L._snap_zc(y, idx, sr)
+    assert abs(y[snapped]) < abs(y[idx]) + 1e-3   # lands on/near a crossing
+    assert abs(snapped - idx) <= int(sr * L.ZC_TOL_MS / 1000)
+
+
+def test_beat_grid_recovers_tempo_not_its_octave():
+    # 120 bpm click train; the fix must report ~120, NOT ~240
+    sr = 22050
+    spb = 0.5
+    n = int(sr * spb * 16)
+    y = np.zeros(n, dtype=np.float32)
+    for b in range(16):
+        i = int(b * spb * sr)
+        y[i:i + 200] += np.hanning(200).astype(np.float32)
+    _, times, bpm = L.beat_grid(y, sr)
+    assert len(times) >= 8
+    med = float(np.median(bpm))
+    assert 100 <= med <= 140                       # not the 240 octave