feat(foundry): correlation engine — recover loop provenance at scale (#12)

engine/correlate.py recovers where each hand-cut loop was cut from: a loop is a
literal slice of a separated stem, so normalized cross-correlation of the waveform
has a razor peak at the true offset (true slice ~1.0, foreign ~0.03). build_provenance.py
runs it over the whole corpus → provenance.json.

Result (honest): 482/1690 loops matched (29%) across 28 kits — only the recent
demucs-derived kits have a source in separated/; the rest predate the workflow, and
we quantify that rather than pretend. Recoveries are exact and name-confirming:
bumbum←MC Fioti "Bum Bum Tam Tam", diams_dj←Diam's "DJ", like_sugar←Chaka Khan
"Like Sugar", praise←Fatboy Slim. Cross-stem alignment 0.34 — the measured answer to
§0's open question (loops share a window across stems about a third of the time).

Two bugs caught and fixed during the build, both load-bearing:
- NAIVE all-pairs full-res NCC took 56 MIN on 1690×120. Fix: coarse-to-fine —
  rank every stem by NCC on a frame-RMS *energy envelope* (phase-robust, ~40× cheaper),
  then run the precise sample-domain NCC only on the top-3 candidates. Minutes, not an hour.
- SILENT-WINDOW BLOWUP: catastrophic cancellation in the sliding-variance cumsum
  drove the normaliser to ~0 in quiet stem regions, so silence scored ∞ and
  masqueraded as a perfect match (an early run reported a bogus 55%/0.69 — ~440 false
  positives). Fix: floor the window norm at 5% of the stem's global energy and clamp
  NCC to [-1,1]. Scores are now all ≤1.0 (max 1.000, median 0.969); the honest 29%
  is what survived. Regression test added.

forkserver gotcha: py3.14's default mp start method re-imports the module per worker
(empty globals), so the preloaded-stems COW trick silently fails — pinned mp fork
context. 8 correlate tests (NCC offset/foreign/gain-invariance/silent-window, envelope
prefilter, locate top-k). Feeds the tierlist cut-critique and the finder recall eval (#7).

armada/tide-table/build_provenance.py

+#!/usr/bin/env python3
+"""build_provenance — recover where each hand-cut loop was cut from (task #12).
+
+Correlate every Samples loop against every separated demucs stem (engine/correlate,
+DRY) to recover its source: track, stem/role, start offset, length. A hand-cut
+loop is a literal slice of a stem, so a clean NCC match is high-confidence
+ground truth (PHASE2_FINDINGS §5·0). Two payoffs beyond the provenance map:
+
+  • coverage, quantified — only recent demucs-derived kits have a source in
+    separated/; we report matched/total, never pretend full.
+  • cross-stem alignment — the actual answer to "how often are a kit's loops cut
+    at the SAME time window across stems?" (the §0 tendency, measured): for each
+    (kit, source-track) we cluster matched offsets and report the share that
+    co-occur within a bar's slack.
+
+Stems are preloaded ONCE (pickle cache) and shared COW across workers; loops are
+short and load per-task.
+
+    python3 build_provenance.py            # → provenance.json
+    python3 build_provenance.py --jobs 8 --rebuild-cache
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import os
+import pickle
+import multiprocessing as mp
+import sys
+from collections import defaultdict
+from concurrent.futures import ProcessPoolExecutor
+from pathlib import Path
+
+# Python 3.14 defaults to forkserver, which re-imports this module in each worker
+# with an EMPTY STEMS global. We rely on fork's copy-on-write to share the
+# preloaded stems set before the pool starts — so pin the fork context.
+_FORK = mp.get_context("fork")
+
+HERE = Path(__file__).resolve().parent
+ROOT = HERE.parent.parent
+FOUNDRY = ROOT / "tools" / "foundry"
+sys.path.insert(0, str(FOUNDRY))
+
+from engine import correlate as C       # noqa: E402
+
+HOME = Path.home()
+SAMPLES = Path(os.environ.get("FOUNDRY_SAMPLES", HOME / "Work/Sound/Samples"))
+# Where whole-track demucs output lives (main store + co-located per-kit).
+STEM_GLOBS = [
+    HOME / "Downloads/separated/htdemucs",
+    SAMPLES,                            # **/separated/htdemucs/*/ caught below
+]
+CACHE = HERE / "_provenance_stems.pkl"
+OUT = HERE / "provenance.json"
+AS_OF = os.environ.get("PROVENANCE_AS_OF", "")
+
+STEMS: dict[str, "object"] = {}        # module global → COW-shared into workers
+
+
+def gather_stem_paths() -> dict[str, Path]:
+    """{ '<track>/<stem>': path } across all known separated/ locations."""
+    out: dict[str, Path] = {}
+    for base in STEM_GLOBS:
+        if not base.exists():
+            continue
+        # base may itself BE an htdemucs/ dir (the main store) or CONTAIN them
+        # nested under each kit (the co-located per-kit separated/ dirs).
+        wavs = base.glob("*/*.wav") if base.name == "htdemucs" else base.rglob("htdemucs/*/*.wav")
+        for wav in wavs:
+            track = wav.parent.name
+            sid = f"{track}/{wav.stem}"
+            out.setdefault(sid, wav)
+    return out
+
+
+def load_stems(rebuild: bool) -> dict[str, "object"]:
+    if CACHE.exists() and not rebuild:
+        with open(CACHE, "rb") as f:
+            return pickle.load(f)
+    paths = gather_stem_paths()
+    print(f"loading + resampling {len(paths)} stems to {C.TARGET_SR}Hz…", file=sys.stderr)
+    stems = {}
+    for i, (sid, p) in enumerate(sorted(paths.items()), 1):
+        try:
+            stems[sid] = C.load_mono(p)
+        except Exception as e:
+            print(f"  skip {sid}: {e}", file=sys.stderr)
+        if i % 20 == 0:
+            print(f"  {i}/{len(paths)}", file=sys.stderr)
+    with open(CACHE, "wb") as f:
+        pickle.dump(stems, f, protocol=pickle.HIGHEST_PROTOCOL)
+    return stems
+
+
+def kit_wavs(kit_dir: Path):
+    return sorted(p for p in kit_dir.iterdir()
+                  if p.is_file() and p.suffix.lower() == ".wav")
+
+
+def _locate_one(args):
+    kit, path = args
+    try:
+        hit = C.locate(path, STEMS)
+    except Exception as e:
+        return {"kit": kit, "loop": Path(path).name, "error": str(e)}
+    row = {"kit": kit, "loop": Path(path).name}
+    if hit:
+        track, _, role = hit["stem"].partition("/")
+        row.update(track=track, role=role, score=hit["score"],
+                   offset_s=hit["offset_s"], length_s=hit["length_s"])
+    return row
+
+
+def cross_stem_alignment(matches: list[dict], bar_slack_s: float = 2.0) -> dict:
+    """Within each (kit, track), what share of matched loops share an offset?
+
+    Answers §0: loops are *often* (not always) cut at the same window across stems.
+    A loop is 'aligned' if another matched loop in the same (kit,track) starts
+    within bar_slack_s. Reports the aligned fraction over all matched loops.
+    """
+    groups = defaultdict(list)
+    for m in matches:
+        if "offset_s" in m:
+            groups[(m["kit"], m["track"])].append(m["offset_s"])
+    aligned = total = 0
+    multi_stem_windows = 0
+    for offs in groups.values():
+        offs = sorted(offs)
+        for i, o in enumerate(offs):
+            total += 1
+            if any(j != i and abs(o - offs[j]) <= bar_slack_s for j in range(len(offs))):
+                aligned += 1
+        # count distinct windows that gather ≥2 stems
+        used = [False] * len(offs)
+        for i in range(len(offs)):
+            if used[i]:
+                continue
+            cluster = [j for j in range(len(offs)) if abs(offs[j] - offs[i]) <= bar_slack_s]
+            if len(cluster) >= 2:
+                multi_stem_windows += 1
+            for j in cluster:
+                used[j] = True
+    return {"aligned_fraction": round(aligned / total, 4) if total else None,
+            "matched_in_groups": total, "multi_stem_windows": multi_stem_windows}
+
+
+def main(argv=None) -> int:
+    global STEMS
+    ap = argparse.ArgumentParser(description="recover loop→stem provenance over the corpus")
+    ap.add_argument("--jobs", type=int, default=max(1, (os.cpu_count() or 4) - 2))
+    ap.add_argument("--rebuild-cache", action="store_true")
+    ap.add_argument("--limit", type=int, default=0)
+    a = ap.parse_args(argv)
+
+    STEMS = C.prepare(load_stems(a.rebuild_cache))   # {id: (full, envelope)} for coarse-to-fine
+    if not STEMS:
+        print("no stems found — nothing to correlate", file=sys.stderr)
+        return 1
+    tracks = sorted({sid.split("/")[0] for sid in STEMS})
+    print(f"{len(STEMS)} stems across {len(tracks)} tracks", file=sys.stderr)
+
+    kits = sorted(d for d in SAMPLES.iterdir() if d.is_dir())
+    if a.limit:
+        kits = kits[: a.limit]
+    jobs = [(k.name, str(w)) for k in kits for w in kit_wavs(k)]
+    print(f"correlating {len(jobs)} loops on {a.jobs} workers…", file=sys.stderr)
+
+    with ProcessPoolExecutor(max_workers=a.jobs, mp_context=_FORK) as ex:
+        rows = list(ex.map(_locate_one, jobs, chunksize=4))
+
+    matches = [r for r in rows if "offset_s" in r]
+    by_kit = defaultdict(lambda: {"loops": 0, "matched": 0})
+    for r in rows:
+        by_kit[r["kit"]]["loops"] += 1
+        if "offset_s" in r:
+            by_kit[r["kit"]]["matched"] += 1
+    kits_with_any = sum(1 for v in by_kit.values() if v["matched"])
+
+    report = {
+        "schema": "foundry provenance v1",
+        "as_of": AS_OF,
+        "engine": "tools/foundry/engine/correlate.py (NCC, sr=%d, min_score=%.2f)"
+                  % (C.TARGET_SR, C.MIN_SCORE),
+        "n_stems": len(STEMS),
+        "n_source_tracks": len(tracks),
+        "n_loops": len(jobs),
+        "n_matched": len(matches),
+        "coverage": round(len(matches) / len(jobs), 4) if jobs else 0,
+        "kits_with_a_source": kits_with_any,
+        "cross_stem_alignment": cross_stem_alignment(matches),
+        "by_kit": {k: v for k, v in sorted(by_kit.items()) if v["matched"]},
+        "matches": sorted(matches, key=lambda m: (m["kit"], m["loop"])),
+    }
+    OUT.write_text(json.dumps(report, indent=2))
+    al = report["cross_stem_alignment"]["aligned_fraction"]
+    print(f"✓ {OUT.name}: {len(matches)}/{len(jobs)} loops matched "
+          f"({report['coverage']*100:.0f}%) across {kits_with_any} kits | "
+          f"cross-stem aligned {al}", file=sys.stderr)
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tools/foundry/engine/correlate.py

+"""correlate — recover a loop's provenance by finding the stem it was cut from.
+
+A hand-cut ParVagues loop is *literally a slice of a separated stem* (PLN keeps
+stems as-cut, no normalize — PHASE2_FINDINGS §0/§5·0). So normalized cross-
+correlation of the loop's waveform against a stem's waveform has a razor-sharp
+peak at the true offset: a real slice scores ~1.0, a foreign clip ~0.03. That
+peak gives us, for free and at scale: which source track, which stem, the start
+offset, and the length — dozens of ground-truth (sample → source) pairs the
+finder (#7) and the tierlist cut-critique (#11) can lean on.
+
+Only the ~recent demucs-derived kits will match (separated/ holds far fewer
+tracks than the kit count); the driver QUANTIFIES coverage rather than pretending
+every loop has a recoverable source.
+
+`ncc_locate` is pure (numpy + scipy, tested on synthetic slices). Runs under
+system python3; no GPU, no network.
+
+    python3 -m engine.correlate <loop.wav> <stem.wav> [more-stems…]   # locate one loop
+"""
+from __future__ import annotations
+
+from pathlib import Path
+from typing import Optional
+
+import numpy as np
+from scipy.signal import fftconvolve
+
+TARGET_SR = 4000          # plenty of waveform shape for a discriminative match; fast
+MIN_SCORE = 0.55          # below this, treat as "no source found" (foreign ≈ 0.03)
+_EPS = 1e-9
+
+
+def load_mono(path: str | Path, sr: int = TARGET_SR) -> np.ndarray:
+    """Load any audio → mono float32 at `sr` (soxr-backed resample via librosa)."""
+    import librosa
+    y, _ = librosa.load(str(path), sr=sr, mono=True)
+    return y.astype(np.float32)
+
+
+def ncc_locate(loop: np.ndarray, stem: np.ndarray) -> tuple[float, int]:
+    """Normalized cross-correlation peak of `loop` within `stem`.
+
+    Returns (score 0..1, offset in samples). score≈1 ⇒ loop is that slice of stem.
+    If the loop is longer than the stem, returns (0, 0).
+    """
+    L = loop.size
+    if L == 0 or stem.size < L:
+        return 0.0, 0
+    loop = loop - loop.mean()
+    ln = float(np.sqrt((loop.astype(np.float64) ** 2).sum())) + _EPS
+    # numerator: sliding dot product (cross-correlation) via FFT
+    num = fftconvolve(stem, loop[::-1], mode="valid")
+    # denominator: sliding L2 norm of the (mean-removed) stem window
+    st = stem.astype(np.float64)
+    c1 = np.cumsum(np.r_[0.0, st])
+    c2 = np.cumsum(np.r_[0.0, st * st])
+    wsum = c1[L:] - c1[:-L]
+    wsq = c2[L:] - c2[:-L]
+    wvar = np.maximum(wsq - wsum * wsum / L, 0.0)      # window L2² (mean-removed)
+    # Floor the window norm relative to the stem's overall energy so near-silent
+    # windows (where catastrophic cancellation drives wvar→0) can't blow the ratio
+    # up — otherwise a quiet gap scores ∞ and masquerades as a perfect match.
+    gstd = float(st.std()) + _EPS
+    floor = np.sqrt(L) * 0.05 * gstd
+    denom = np.maximum(np.sqrt(wvar), floor) * ln
+    ncc = np.clip(num / denom, -1.0, 1.0)
+    i = int(np.argmax(ncc))
+    return float(ncc[i]), i
+
+
+def envelope(y: np.ndarray, sr: int, fps: int = 100) -> np.ndarray:
+    """Frame-RMS energy envelope — a cheap, phase-robust fingerprint for prefiltering.
+
+    A true slice's envelope matches a contiguous span of the source's envelope, so
+    NCC on envelopes ranks candidate stems at ~1/40th the sample-domain cost (and
+    unlike naive decimation, it doesn't depend on the offset being phase-aligned).
+    """
+    hop = max(1, sr // fps)
+    n = y.size // hop
+    if n == 0:
+        return np.sqrt(np.mean(y ** 2, keepdims=True)) if y.size else np.zeros(1)
+    frames = y[: n * hop].reshape(n, hop)
+    return np.sqrt((frames.astype(np.float64) ** 2).mean(axis=1) + _EPS)
+
+
+def prepare(stems: dict[str, np.ndarray], *, sr: int = TARGET_SR) -> dict[str, tuple]:
+    """Precompute (full, envelope) per stem once, so locate() can coarse-rank cheaply."""
+    return {sid: (arr, envelope(arr, sr)) for sid, arr in stems.items()}
+
+
+def locate(loop_path: str | Path, prepared: dict[str, tuple],
+           *, sr: int = TARGET_SR, min_score: float = MIN_SCORE,
+           topk: int = 3) -> Optional[dict]:
+    """Find the stem a loop was cut from. Coarse-to-fine:
+
+    1. rank every stem by NCC on the energy envelope (cheap);
+    2. run the precise sample-domain NCC only on the top-`topk` candidates.
+
+    `prepared` comes from prepare(); returns {stem, score, offset_s, length_s} or None.
+    """
+    loop = load_mono(loop_path, sr)
+    lenv = envelope(loop, sr)
+    coarse = []
+    for sid, (arr, senv) in prepared.items():
+        if senv.size < lenv.size or arr.size < loop.size:
+            continue
+        sc, _ = ncc_locate(lenv, senv)
+        coarse.append((sc, sid))
+    if not coarse:
+        return None
+    coarse.sort(reverse=True)
+    best = None
+    for _, sid in coarse[:topk]:
+        arr = prepared[sid][0]
+        score, off = ncc_locate(loop, arr)
+        if best is None or score > best[1]:
+            best = (sid, score, off)
+    if best is None or best[1] < min_score:
+        return None
+    sid, score, off = best
+    return {"stem": sid, "score": round(score, 4),
+            "offset_s": round(off / sr, 3), "length_s": round(loop.size / sr, 3)}
+
+
+def _main(argv=None) -> int:
+    import argparse
+    import json
+
+    ap = argparse.ArgumentParser(description="locate the stem a loop was cut from")
+    ap.add_argument("loop")
+    ap.add_argument("stems", nargs="+", help="candidate stem wavs")
+    ap.add_argument("--min-score", type=float, default=MIN_SCORE)
+    a = ap.parse_args(argv)
+
+    prepared = prepare({f"{Path(s).parent.name}/{Path(s).stem}": load_mono(s) for s in a.stems})
+    hit = locate(a.loop, prepared, min_score=a.min_score)
+    print(json.dumps(hit, indent=2) if hit else "no source found above threshold")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(_main())

tools/foundry/tests/test_correlate.py

+"""Tests for the provenance correlator (engine/correlate.py).
+
+Synthetic: embed a known slice in a longer 'stem' and recover its offset; confirm
+a foreign clip scores near zero. No audio files, no librosa needed for the math.
+"""
+import sys
+from pathlib import Path
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+
+from engine import correlate as C
+
+
+def _stem(n=200_000, seed=0):
+    return np.random.default_rng(seed).standard_normal(n).astype(np.float32)
+
+
+def test_ncc_recovers_exact_offset():
+    stem = _stem()
+    off, L = 73_210, 16_000
+    loop = stem[off:off + L].copy()
+    score, found = C.ncc_locate(loop, stem)
+    assert found == off
+    assert score > 0.99
+
+
+def test_ncc_foreign_clip_scores_low():
+    stem = _stem(seed=1)
+    foreign = np.random.default_rng(99).standard_normal(16_000).astype(np.float32)
+    score, _ = C.ncc_locate(foreign, stem)
+    assert score < 0.2
+
+
+def test_ncc_loop_longer_than_stem():
+    assert C.ncc_locate(_stem(1000), _stem(500)) == (0.0, 0)
+
+
+def test_ncc_silent_windows_dont_blow_up():
+    # a stem that's mostly silence with one real region; a FOREIGN loop must score
+    # low, not ∞ (silent windows underflow the normaliser → spurious perfect match).
+    rng = np.random.default_rng(0)
+    stem = np.zeros(200_000, dtype=np.float32)
+    stem[80_000:120_000] = rng.standard_normal(40_000).astype(np.float32)
+    foreign = rng.standard_normal(12_000).astype(np.float32)
+    score, _ = C.ncc_locate(foreign, stem)
+    assert score <= 1.0 and score < 0.3
+    # the true slice in the loud region still resolves exactly
+    loop = stem[90_000:102_000].copy()
+    s2, off = C.ncc_locate(loop, stem)
+    assert off == 90_000 and s2 > 0.99
+
+
+def test_ncc_survives_gain_change():
+    # the loop graded as-cut, but a global gain shouldn't break the match (NCC is
+    # scale-invariant) — guards against a future "PLN nudged the level" case.
+    stem = _stem(seed=2)
+    off, L = 40_000, 12_000
+    loop = stem[off:off + L].copy() * 0.5
+    score, found = C.ncc_locate(loop, stem)
+    assert found == off and score > 0.99
+
+
+def test_locate_picks_right_stem(monkeypatch):
+    stems = {"trackA/bass": _stem(seed=3), "trackB/vocals": _stem(seed=4)}
+    off = 50_000
+    true_loop = stems["trackB/vocals"][off:off + 10_000].copy()
+    # bypass file IO: feed the array straight through load_mono
+    monkeypatch.setattr(C, "load_mono", lambda p, sr=C.TARGET_SR: true_loop)
+    hit = C.locate("ignored.wav", C.prepare(stems, sr=4000), sr=4000)
+    assert hit["stem"] == "trackB/vocals"
+    assert hit["score"] > 0.99
+    assert hit["offset_s"] == round(off / 4000, 3)
+
+
+def test_locate_returns_none_below_threshold(monkeypatch):
+    stems = {"trackA/bass": _stem(seed=5)}
+    foreign = np.random.default_rng(7).standard_normal(8_000).astype(np.float32)
+    monkeypatch.setattr(C, "load_mono", lambda p, sr=C.TARGET_SR: foreign)
+    assert C.locate("ignored.wav", C.prepare(stems)) is None
+
+
+def test_envelope_matches_slice_offset():
+    # the coarse prefilter: a slice's envelope should rank its own source highest
+    rng = np.random.default_rng(11)
+    stem = rng.standard_normal(120_000).astype(np.float32)
+    loop = stem[30_000:42_000].copy()
+    se, le = C.envelope(stem, 4000), C.envelope(loop, 4000)
+    score, _ = C.ncc_locate(le, se)
+    assert score > 0.9          # envelope NCC still discriminates the true source