docs: VAD benchmarks (Silero vs the Parakeet head vs whisper.cpp) with reproducible scripts

README.md

@@ -85,6 +85,8 @@ ternary ones.
   functions, as a stream (`parakeet_capi_vad_stream_*`), and as the cutter for any ASR model:
   `parakeet-cli transcribe --vad --vad-model silero.gguf`. See [`docs/vad.md`](docs/vad.md).
 - VAD-only slices of Ultra and Redux (6 to 10 MB, `redux-vad.gguf` and `ultra-vad-q8_0.gguf` in the same repo) hold just the head and its front end. They run `vad` and the `parakeet_capi_vad_*` calls and cannot transcribe. See [`docs/vad.md`](docs/vad.md).
+- Accuracy, speed and size of both detectors, with the method and the scripts to repeat
+  them: [`docs/vad-benchmarks.md`](docs/vad-benchmarks.md).
 ---
 
 ## Performance

docs/vad.md

@@ -12,6 +12,9 @@ Both give one speech probability per frame. A shared segmenter
 transcription segments. The segmenter takes the frame period as a parameter
 (`SegmenterOpts::frame_sec`); only the option defaults differ per detector.
 
+Accuracy and speed numbers for both detectors, with the method and the scripts to
+repeat them, are in [vad-benchmarks.md](vad-benchmarks.md).
+
 ## Standalone VAD API
 
 One call returns the speech regions of a clip as JSON, for either detector. Load
@@ -242,7 +245,7 @@ onnxruntime 1.27.0 and the ONNX file of version 6.2.3:
 
 On one CPU thread, `tests/silero_vad_probe` measured about 40 us per chunk
 (F32 or F16, 16 kHz or 8 kHz) on the development machine. That is one run on one
-machine, not a benchmark.
+machine, not a benchmark. For benchmarks, see [vad-benchmarks.md](vad-benchmarks.md).
 
 ## Tests
 

scripts/vad_bench/README.md

@@ -0,0 +1,75 @@
+# VAD benchmark scripts
+
+Scripts behind [docs/vad-benchmarks.md](../../docs/vad-benchmarks.md). The page says what
+each number means and which caveats apply. No audio and no model files are committed;
+`results/` holds the small result files of the runs on that page.
+
+All paths are arguments. Run from the repository root. `PARAKEET_LIB` is the path of
+`libparakeet.so` (build with `-DPARAKEET_SHARED=ON`; the default is `build/libparakeet.so`).
+
+Python packages: `numpy soundfile librosa datasets torch onnxruntime silero-vad`
+(the Silero reference is `silero-vad` 6.2.3). Models: the ASR GGUFs of
+`moondream/parakeet-ultra` (Q8_0 and F16), `moondream/parakeet-redux` (packed `keep` and
+dequantized F16), a TDT 0.6B v3 F16 GGUF for the reference of the TED talks, and the
+Silero GGUFs made with `scripts/convert_silero_vad_to_gguf.py`.
+
+## Data
+
+| Script | What |
+| --- | --- |
+| `fetch_data.py OUT [--libri-count N] [--ted-seconds S]` | LibriSpeech test-clean (every 13th utterance) as `libri.npy`, and the first 600 s of one TED-LIUM talk |
+| `fetch_ted_talks.py OUT [N]` | the first N TED-LIUM long-form talks under 1500 s as `ted_<name>.wav` |
+| `make_clips.py DATA` | the 120 synthetic clips (clean, white 10 dB, pink 10 dB) of the Silero parity run |
+| `vad_synth.py` | clip builder, noise, and metrics (frame F1, kappa, boundary error) |
+
+## Parakeet heads against Silero
+
+```
+scripts/vad_bench/compare_synthetic.py DATA GGUF_DIR out/synthetic_raw.json 24
+scripts/vad_bench/analyze_synthetic.py out/synthetic_raw.json
+scripts/vad_bench/ted_reference.py DATA tdt-0.6b-v3-f16.gguf      # words of each talk
+scripts/vad_bench/compare_ted.py DATA GGUF_DIR out/ted_metrics.json
+taskset -c 0-7 scripts/vad_bench/speed_heads.py talk.wav GGUF_DIR   # timing, quiet machine only
+```
+
+`GGUF_DIR` holds `ultra-q8_0.gguf`, `ultra-f16.gguf`, `redux-keep.gguf` and
+`redux-deq-f16.gguf`. `compare_ted.py` prints speeds; ignore them, use `speed_heads.py`.
+
+## Silero in whisper.cpp and in parakeet.cpp
+
+Build whisper.cpp (CPU), then the test program, which uses its public `whisper_vad_*` API:
+
+```
+g++ -O2 -std=c++17 scripts/vad_bench/wvad.cpp -I whisper.cpp/include -I whisper.cpp/ggml/include \
+    -L whisper.cpp/build/bin -lwhisper -Wl,-rpath,$PWD/whisper.cpp/build/bin -o wvad
+python3 scripts/vad_bench/silero_collect.py CORE DATA MODEL_DIR build/examples/cli/parakeet-cli ./wvad
+python3 scripts/vad_bench/silero_analyze.py DATA      # run where probs.npz was written
+python3 scripts/vad_bench/speed_silero.py CORE ROUNDS DATA/ted300.wav MODEL_DIR \
+    build/examples/cli/parakeet-cli ./wvad whisper.cpp/build/bin out/speed.json
+```
+
+`ted300.wav` is the first 300 s of `ted_talk.wav`. `MODEL_DIR` holds `silero_vad.onnx`,
+`silero-vad-f32.gguf`, `silero-vad-f16.gguf`, `ggml-silero-v6.2.3-ggml.bin`,
+`ggml-silero-v6.2.0.bin` and, for the end to end job, `ggml-tiny.en.bin`.
+`silero_analyze.py` expects `probs.npz` and `segs_own_default.json` in the current directory.
+
+## Long talks, batched decode
+
+```
+scripts/vad_bench/longform_b1.sh TALKS_DIR GGUF_DIR BUILDS_DIR out.tsv loadlog.txt
+```
+
+`BUILDS_DIR` has `build-base`, `build-b0`, `build-b1`: builds of the three commits to compare.
+
+## Timing rules
+
+Timing scripts serialize on a lock file (`/tmp/pk-bench.lock`, `flock`), pin cores with
+`taskset`, and record the load average at the start of every run. Run them only on a
+quiet machine, and keep the load log with the result. Do not quote a timing from a
+loaded machine as a benchmark.
+
+## Slice only head (PR 87)
+
+`vad_slice_bench.cpp` is built inside the PR 87 tree (it uses internal headers) and driven
+by `slice_bench.py speed|load|rss`; see its docstring for the environment variables.
+Results are in `results/slice_*`.

scripts/vad_bench/analyze_synthetic.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+"""usage: analyze_synthetic.py synthetic_raw.json  -> tables on stdout"""
+import json,sys,numpy as np
+from heads_lib import mask,counts,prf,kappa,bound_err
+R=json.load(open(sys.argv[1]))
+conds=R["conds"];refs=R["refs"];P=R["preds"]
+dur=lambda i: refs[i][-1][1]+1.0  # unused
+def total_sec(i): 
+    return None
+import soundfile
+# clip length: need exact; recompute from the raw: last pred... use max(ref end, pred end)+2 s margin rounding => use mask length from ref end + 1.5 (trailing silence >=0.5)
+def n_sec(i): return max(refs[i][-1][1]+0.5, max([e for s,e in (P['silero_matched'][i] or [(0,0)])]+[0]))+0.0
+order=["clean","white20","white10","white5","white0","pink20","pink10","pink5","pink0"]
+def agg(sysname,cs):
+    c=np.zeros(4);bs=[];be=[];ms=0;me=0;nb=0
+    for i,cn in enumerate(conds):
+        if cn not in cs or P[sysname][i] is None: continue
+        L=n_sec(i)+3
+        ref=mask(refs[i],L);pr=mask(P[sysname][i],L)
+        c+=counts(pr,ref)
+        e1,m1=bound_err(P[sysname][i],refs[i],0);e2,m2=bound_err(P[sysname][i],refs[i],1)
+        bs+=e1;be+=e2;ms+=m1;me+=m2;nb+=len(refs[i])
+    return c,bs,be,ms,me,nb
+out={}
+systems=["silero_matched","silero_default","ultra_q8","redux_packed"]
+print("FRAME-LEVEL (10 ms grid), thr 0.5")
+print("cond | "+" | ".join(f"{s} P/R/F1" for s in systems))
+rows=[]
+for cn in order:
+    cells=[]
+    for s in systems:
+        c,*_=agg(s,[cn]);p,r,f=prf(c);cells.append(f"{p*100:.1f}/{r*100:.1f}/{f*100:.1f}");out[f"{cn}|{s}"]=[p,r,f]
+    print(cn," | ".join(cells))
+print("\nBOUNDARY ERROR ms (median / p90 / miss%) start;end, pooled")
+for cn in ["clean","white10","white0","pink10","pink0"]:
+    cells=[]
+    for s in systems:
+        c,bs,be,ms,me,nb=agg(s,[cn])
+        f=lambda a:(f"{np.median(a):.0f}/{np.percentile(a,90):.0f}" if a else "-")
+        cells.append(f"S {f(bs)} ({100*ms/nb:.1f}%) E {f(be)} ({100*me/nb:.1f}%)")
+    print(cn," | ".join(cells))
+print("all-conditions boundary:")
+for s in systems:
+    c,bs,be,ms,me,nb=agg(s,order); a=bs+be
+    print(s,f"median {np.median(a):.0f} p90 {np.percentile(a,90):.0f} start-med {np.median(bs):.0f} end-med {np.median(be):.0f} miss {100*(ms+me)/(2*nb):.1f}%")
+print("\nTHRESHOLD SWEEP F1 (clean,white10,white0,pink10 pooled)  P/R/F1")
+sub=["clean","white10","white0","pink10"]
+for base in ["silero_matched","ultra_q8","redux_packed"]:
+    cells=[]
+    for t in ("@0.3","","@0.7"):
+        nm=base+t; c,*_=agg(nm,sub); p,r,f=prf(c); cells.append(f"thr{t or '@0.5'} {p*100:.1f}/{r*100:.1f}/{f*100:.1f}")
+    print(base," | ".join(cells))
+print("\nAGREEMENT silero_matched vs parakeet (all synthetic): F1 of speech mask, kappa, IoU")
+for s in ["ultra_q8","redux_packed"]:
+    c=np.zeros(4)
+    for i in range(len(conds)):
+        L=n_sec(i)+3;c+=counts(mask(P[s][i],L),mask(P["silero_matched"][i],L))
+    p,r,f=prf(c);print(s,f"F1 {f*100:.1f} kappa {kappa(c):.3f} IoU {c[0]/(c[0]+c[1]+c[2])*100:.1f}")
+c=np.zeros(4)
+for i in range(len(conds)):
+    L=n_sec(i)+3;c+=counts(mask(P["ultra_q8"][i],L),mask(P["redux_packed"][i],L))
+print("ultra vs redux",f"F1 {prf(c)[2]*100:.1f} kappa {kappa(c):.3f}")

scripts/vad_bench/compare_synthetic.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+"""Synthetic benchmark: Parakeet VAD head (Ultra Q8_0, Redux packed) vs Silero.
+
+usage: PARAKEET_LIB=build/libparakeet.so compare_synthetic.py DATA_DIR GGUF_DIR OUT.json [CLIPS_PER_COND]
+GGUF_DIR holds ultra-q8_0.gguf and redux-keep.gguf. Needs: numpy soundfile torch silero-vad onnxruntime.
+Writes the raw spans of every system; analyze_synthetic.py turns them into the tables.
+"""
+import json,sys,numpy as np,soundfile as sf
+from heads_lib import *
+SYS={}
+def make_systems():
+    S={"silero_matched":lambda y:silero(y,True),"silero_default":lambda y:silero(y,False)}
+    return S
+DATA,G,OUT=sys.argv[1],sys.argv[2],sys.argv[3]   # data dir, gguf dir, output json
+N_CLIPS=int(sys.argv[4]) if len(sys.argv)>4 else 24   # clips per condition
+libri=list(np.load(DATA+"/libri.npy",allow_pickle=True))
+rng=np.random.default_rng(7)
+order=rng.permutation(len(libri))
+sets=[[libri[i] for i in order[k*5:(k+1)*5]] for k in range(len(libri)//5)][:N_CLIPS]
+conds=[("clean",None)]+[(k,s) for k in("white","pink") for s in(20,10,5,0)]
+# build corpus once
+corpus=[]
+for ci,(kind,snr) in enumerate(conds):
+    for k,us in enumerate(sets):
+        r=np.random.default_rng(1000+k)            # same speech layout per clip across conditions
+        y,spans=build(us,r)
+        y=add_noise(y,spans,kind,snr,np.random.default_rng(5000+ci*100+k))
+        corpus.append((f"{kind}{'' if snr is None else snr}",y,spans))
+print("corpus",len(corpus),"wavs",sum(len(c[1]) for c in corpus)/16000,"s",flush=True)
+models={"ultra_q8":"ultra-q8_0.gguf","redux_packed":"redux-keep.gguf"}
+preds={}   # system -> list of spans per corpus item
+for name,fn in make_systems().items():
+    preds[name]=[fn(c[1]) for c in corpus]
+print("silero done",flush=True)
+for mname,f in models.items():
+    pk=PK(G+"/"+f)
+    preds[mname]=[parakeet(pk,c[1]) for c in corpus]
+    # threshold sweep on a subset of conditions
+    for thr in (0.3,0.7):
+        preds[f"{mname}@{thr}"]=[parakeet(pk,c[1],thr) if c[0] in("clean","white10","white0","pink10") else None for c in corpus]
+    print(mname,"done",flush=True)
+# silero thresholds
+for thr in (0.3,0.7):
+    preds[f"silero_matched@{thr}"]=[None]*len(corpus)
+    for i,c in enumerate(corpus):
+        if c[0] in("clean","white10","white0","pink10"):
+            ts=get_speech_timestamps(torch.from_numpy(c[1]),SIL,sampling_rate=16000,return_seconds=True,threshold=thr,min_speech_duration_ms=100,min_silence_duration_ms=200,speech_pad_ms=0)
+            preds[f"silero_matched@{thr}"][i]=[(t["start"],t["end"]) for t in ts]
+json.dump({"conds":[c[0] for c in corpus],"refs":[c[2] for c in corpus],"preds":preds},open(OUT,"w"))
+print("saved",flush=True)

scripts/vad_bench/compare_ted.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+"""usage: compare_ted.py DATA_DIR GGUF_DIR OUT.json
+DATA_DIR holds ted_*.wav and the matching ted_*_words.json from ted_reference.py.
+The speed it prints is NOT a valid benchmark (other load); use speed_heads.py."""
+import json,glob,sys,time,subprocess,numpy as np,soundfile as sf
+from heads_lib import *
+from pk import PK
+DATA,G,OUT=sys.argv[1],sys.argv[2],sys.argv[3]
+print("uptime before:",subprocess.getoutput("uptime"))
+talks=[f for f in sorted(glob.glob(DATA+"/ted_*.wav")) if "S103" not in f]   # the original run left out one talk with S103 in its name; reason not recorded
+ys={f:sf.read(f,dtype="float32")[0] for f in talks}
+refs={}
+for f in talks:
+    W=json.load(open(f.replace(".wav","_words.json")))
+    spans=[];
+    for s,e in W:
+        if spans and s-spans[-1][1]<0.4: spans[-1][1]=max(spans[-1][1],e)
+        else: spans.append([s,e])
+    refs[f]=spans
+res={};preds={}
+def run(name,fn):
+    cs=np.zeros(4);tt=0;ds=0;
+    for f in talks:
+        y=ys[f];t0=time.perf_counter();p=fn(y);dt=time.perf_counter()-t0
+        preds[(name,f)]=p;tt+=dt;ds+=len(y)/16000
+        L=len(y)/16000+1;cs+=counts(mask(p,L),mask(refs[f],L))
+    P,R,F=prf(cs);res[name]=dict(P=P,R=R,F1=F,rtf_x=ds/tt,sec=tt,audio=ds,speech_frac=None)
+    print(name,"P %.1f R %.1f F1 %.1f  VAD speed %.0fx real time (%.1f s for %.0f s)"%(P*100,R*100,F*100,ds/tt,tt,ds),flush=True)
+run("silero_matched",lambda y:silero(y,True))
+run("silero_default",lambda y:silero(y,False))
+for n,f in(("ultra_q8","ultra-q8_0.gguf"),("redux_packed","redux-keep.gguf")):
+    pk=PK(G+"/"+f); run(n,lambda y:parakeet(pk,y)); 
+    out=pk.vad(ys[talks[1]],{"threshold":0.5})
+    del pk
+# reference speech fraction and agreement
+cs=np.zeros(4)
+for f in talks:
+    L=len(ys[f])/16000+1; r=mask(refs[f],L); cs+=np.array([r.sum(),0,0,len(r)-r.sum()])
+print("ref speech fraction %.1f%%"%(100*cs[0]/(cs[0]+cs[3])))
+for n in ["silero_matched","silero_default","ultra_q8","redux_packed"]:
+    t=sum(mask(preds[(n,f)],len(ys[f])/16000+1).sum() for f in talks);a=sum(len(mask([],len(ys[f])/16000+1)) for f in talks)
+    print(n,"speech fraction %.1f%%"%(100*t/a))
+print("AGREEMENT on TED talks")
+for a,b in(("ultra_q8","silero_matched"),("redux_packed","silero_matched"),("ultra_q8","redux_packed")):
+    c=np.zeros(4)
+    for f in talks:
+        L=len(ys[f])/16000+1;c+=counts(mask(preds[(a,f)],L),mask(preds[(b,f)],L))
+    print(a,"vs",b,"F1 %.1f kappa %.3f"%(prf(c)[2]*100,kappa(c)))
+print("talks",[(f,round(len(ys[f])/16000)) for f in talks])
+json.dump(res,open(OUT,"w"),indent=1)
+print("uptime after:",subprocess.getoutput("uptime"))

scripts/vad_bench/fetch_data.py

@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+"""Fetch the public audio used by the VAD benchmarks (streams from the Hugging Face hub).
+
+usage: fetch_data.py OUT_DIR [--libri-count N] [--ted-seconds S]
+
+Writes OUT_DIR/libri.npy (LibriSpeech test-clean, every 13th utterance, N of them) and
+OUT_DIR/ted_talk.wav (the first TED-LIUM long-form talk longer than 700 s, cut to S seconds,
+default 600). The silero parity run used N=200. The vad-compare run used a subset with 24 x 5
+utterances; the script that cut it was not saved, so that corpus is not byte-identical to the
+one in the published tables (see docs/vad-benchmarks.md).
+"""
+import argparse, io, os
+import numpy as np, soundfile as sf, librosa
+from datasets import load_dataset, Audio
+
+ap = argparse.ArgumentParser()
+ap.add_argument("out")
+ap.add_argument("--libri-count", type=int, default=200)
+ap.add_argument("--ted-seconds", type=int, default=600)
+a = ap.parse_args()
+os.makedirs(a.out, exist_ok=True)
+
+ds = load_dataset("openslr/librispeech_asr", "clean", split="test", streaming=True).cast_column("audio", Audio(decode=False))
+U = []
+for i, ex in enumerate(ds):
+    if i % 13:
+        continue
+    y, sr = sf.read(io.BytesIO(ex["audio"]["bytes"]), dtype="float32")
+    assert sr == 16000
+    U.append(y)
+    if len(U) >= a.libri_count:
+        break
+np.save(f"{a.out}/libri.npy", np.array(U, dtype=object), allow_pickle=True)
+print("libri", len(U))
+
+ds = load_dataset("distil-whisper/tedlium-long-form", split="test", streaming=True).cast_column("audio", Audio(decode=False))
+for ex in ds:
+    y, sr = sf.read(io.BytesIO(ex["audio"]["bytes"]), dtype="float32")
+    if y.ndim > 1:
+        y = y.mean(1)
+    if sr != 16000:
+        y = librosa.resample(y, orig_sr=sr, target_sr=16000)
+    if len(y) / 16000 < 700:
+        continue
+    sf.write(f"{a.out}/ted_talk.wav", y[: a.ted_seconds * 16000], 16000, subtype="PCM_16")
+    print("ted", len(y) / 16000)
+    break

scripts/vad_bench/fetch_ted_talks.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python3
+"""usage: fetch_ted_talks.py OUT_DIR [N_TALKS]  -> OUT_DIR/ted_<name>.wav, the first N talks under 1500 s of
+distil-whisper/tedlium-long-form (test). The published run used three talks; which ones is not recorded."""
+import io, re, sys, numpy as np, soundfile as sf, librosa
+from datasets import load_dataset, Audio
+D=sys.argv[1]; N=int(sys.argv[2]) if len(sys.argv)>2 else 3
+# TED-LIUM long-form: first talks
+ds=load_dataset("distil-whisper/tedlium-long-form",split="test",streaming=True).cast_column("audio",Audio(decode=False))
+n=0
+for i,ex in enumerate(ds):
+    y,sr=sf.read(io.BytesIO(ex["audio"]["bytes"]),dtype="float32")
+    if y.ndim>1:y=y.mean(1)
+    if sr!=16000:y=librosa.resample(y,orig_sr=sr,target_sr=16000)
+    if len(y)/16000>1500: continue
+    name=re.sub(r"[^A-Za-z0-9_-]","",(ex["audio"]["path"] or f"talk{i}").split("/")[-1].rsplit(".",1)[0]) or f"talk{i}"
+    sf.write(f"{D}/ted_{name}.wav",y,16000,subtype="PCM_16")
+    print(name,len(y)/16000,flush=True); n+=1
+    if n>=N: break

scripts/vad_bench/heads_lib.py

@@ -0,0 +1,60 @@
+"""Detector wrappers for the Parakeet-head vs Silero comparison.
+Silero reference: the `silero-vad` pip package (6.2.3), ONNX, one thread. Parakeet: libparakeet via pk.py.
+The synthetic-clip helpers are the same code as vad_synth.py (kept apart so each script runs on its own).
+"""
+import numpy as np, torch
+from pk import PK
+torch.set_num_threads(1)
+from silero_vad import load_silero_vad, get_speech_timestamps
+SIL=load_silero_vad(onnx=True)
+GR=0.01  # mask grid 10 ms
+def mask(spans,n_sec):
+    m=np.zeros(int(round(n_sec/GR)),bool)
+    for s,e in spans:
+        m[int(round(s/GR)):int(round(e/GR))]=True
+    return m
+def silero(y,matched=True):
+    kw=dict(threshold=0.5,min_speech_duration_ms=100,min_silence_duration_ms=200,speech_pad_ms=0) if matched else {}
+    ts=get_speech_timestamps(torch.from_numpy(y),SIL,sampling_rate=16000,return_seconds=True,**kw)
+    return [(t["start"],t["end"]) for t in ts]
+def parakeet(pk,y,thr=0.5):
+    return [(s["start"],s["end"]) for s in pk.vad(y,{"threshold":thr})["segments"]]
+def trim(u):
+    fr=len(u)//160; r=np.sqrt((u[:fr*160].reshape(fr,160)**2).mean(1))
+    idx=np.where(r>r.max()*0.01)[0]
+    return u[idx[0]*160:(idx[-1]+1)*160]
+def pink(n,rng):
+    X=np.fft.rfft(rng.standard_normal(n)); f=np.arange(len(X)); f[0]=1
+    y=np.fft.irfft(X/np.sqrt(f),n); return y/y.std()
+def build(utts,rng):
+    g=lambda a,b: rng.uniform(a,b)
+    parts=[];spans=[];t=g(0.5,1.5)
+    parts.append(rng.standard_normal(int(t*16000))*1e-3)
+    for i,u in enumerate(utts):
+        u=trim(u); spans.append((t,t+len(u)/16000)); parts.append(u); t+=len(u)/16000
+        gap=g(0.3,2.0) if i<len(utts)-1 else g(0.5,1.5)
+        parts.append(rng.standard_normal(int(gap*16000))*1e-3); t+=int(gap*16000)/16000
+    return np.concatenate(parts).astype(np.float32),spans
+def add_noise(y,spans,kind,snr,rng):
+    if kind=="clean": return y
+    sp=np.concatenate([y[int(s*16000):int(e*16000)] for s,e in spans]); P=(sp**2).mean()
+    n=rng.standard_normal(len(y)) if kind=="white" else pink(len(y),rng)
+    n=n*np.sqrt(P/10**(snr/10))
+    return np.clip(y+n,-1,1).astype(np.float32)
+def counts(pred,ref):
+    n=min(len(pred),len(ref)); p=pred[:n];r=ref[:n]
+    return np.array([(p&r).sum(),(p&~r).sum(),(~p&r).sum(),(~p&~r).sum()])
+def prf(c):
+    tp,fp,fn,tn=c; P=tp/max(1,tp+fp);R=tp/max(1,tp+fn);return P,R,2*P*R/max(1e-9,P+R)
+def kappa(c):
+    tp,fp,fn,tn=c;n=c.sum();po=(tp+tn)/n;pe=((tp+fp)*(tp+fn)+(fn+tn)*(fp+tn))/n**2;return (po-pe)/(1-pe)
+def bound_err(pred,ref,idx):
+    out=[];miss=0
+    pb=np.array([p[idx] for p in pred]) if pred else np.array([])
+    for r in ref:
+        if len(pb)==0: miss+=1;continue
+        d=np.abs(pb-r[idx]).min()
+        if d>1.0: miss+=1
+        else: out.append(d*1000)
+    return out,miss
+