From 95c777cbc43f93d53500ea407066ca2abd6ba32d Mon Sep 17 00:00:00 2001 From: Chris Fields Date: Fri, 28 Aug 2026 23:16:55 -0500 Subject: [PATCH 01/10] dev: script the #152 thread-sweep analysis Replaces the ad-hoc Python that has been pasted at these logs all session, and encodes the one thing that must not be got wrong. Segments are aligned by CLUSTER INDEX, never by wall time. At the same t=, a 96-thread run and a 48-thread run sit at different cluster positions, and the serial fraction collapses over a run -- so fixed wall-clock windows compare different phases. Doing exactly that produced a confidently reported conclusion with the sign of the trend inverted before it was caught. Reports four things: - run_dada per thread count, with replicate spread beside every mean, and an explicit note when an arm has one replicate that differences under ~2% are not distinguishable from run-to-run variation. - The serial block (store+shuffle+bud+p_update), which is invariant in seconds and rising in share -- 16S R1 is 208.5s at 48 threads and 207.3s at 96, so map work walks toward the Amdahl wall rather than away from it. That is the figure to quote when weighing bandwidth work (#28) against serial work (#154). - Marginal seconds per cluster segment, showing where in the run the benefit or penalty lands. - The map's screen/align split, which predicts the knee. On that last point the calibration is printed rather than a threshold, because four arms is not a rule: 50.8% screen still scales at 96, 67.7% turns over near 64, and 83-86% turns over below 48. It predicts WITHIN a pool as well as across -- 16S's two reads differ by 17 points of screen share and have different knees -- so it tracks arithmetic intensity, not the dataset. The thread numbers are specific to the EPYC 7713; the ordering is what travels. Co-Authored-By: Claude Opus 5 --- dev/analyze_thread_sweep.py | 261 ++++++++++++++++++++++++++++++++++++ 1 file changed, 261 insertions(+) create mode 100755 dev/analyze_thread_sweep.py diff --git a/dev/analyze_thread_sweep.py b/dev/analyze_thread_sweep.py new file mode 100755 index 0000000..27b10ca --- /dev/null +++ b/dev/analyze_thread_sweep.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python3 +"""Analyse a dada2-rs thread sweep from `--verbose` logs (issue #152). + +Answers three questions a thread sweep is run to answer, and refuses to answer +them the wrong way: + +1. **Does this pool benefit from more threads?** `run_dada` per thread count, + averaged over replicates, with the best count marked. On soil ITS2 the + answer is *no* -- 48 -> 96 threads makes it slower -- while soil 16S keeps + improving, so this is genuinely per-dataset and not a property of the + hardware. + +2. **Where in the run does the benefit (or penalty) land?** Marginal seconds + per fixed *cluster* segment, from the progress lines (#150). + + Segments are aligned by **cluster index, never by wall time**. This is the + whole reason the script exists: at the same `t=`, a 96-thread run and a + 48-thread run are at different cluster positions and therefore in different + phases of a workload whose serial fraction collapses over the run. Comparing + fixed wall-clock windows across thread counts produced a confidently reported + conclusion with the *sign of the trend inverted* before this was caught. + +3. **Why?** The map's screen/align split. The k-mer screen streams k-mer + vectors (bandwidth-bound); the aligner runs DP (compute-bound). A + screen-dominated pool saturates at low thread counts; an align-dominated one + keeps scaling. Measured on four arms, the ordering is monotonic: + + 16S R1 50.8% screen -> still scaling at 96 threads + 16S R2 67.7% screen -> knee near 64 + ITS2 R2 83.2% screen -> knee below 48 + ITS2 R1 85.8% screen -> knee below 48 + + It predicts *within* a pool as well as across pools -- 16S's two reads differ + by 17 points of screen share and have different knees -- so this is a property + of the workload's arithmetic intensity, not of the dataset's name. The split + is printed by every verbose run, so it is usable up front. The thread numbers + are specific to this machine (EPYC 7713); the ordering is what travels. + +Also reports the serial block, which is invariant to thread count -- its *share* +rises as the map gets faster, so map optimisation walks toward the Amdahl wall +rather than away from it. That is the number to quote when deciding between +bandwidth work and serial work. + +Usage: + dev/analyze_thread_sweep.py tmp/issue-152/full-pooling-novaseq-ITS + dev/analyze_thread_sweep.py --segments 6 --read R1 + +Expected layout (as produced by the sweep job scripts): + /rep/threads/dada/dada..log +""" + +import argparse +import collections +import glob +import os +import re +import statistics +import sys + +RE_PROGRESS = re.compile(r"progress t=(\d+)s cluster (\d+) ") +RE_RUNDADA = re.compile(r"run_dada=([\d.]+)s") +RE_PHASES = re.compile( + r"compare=([\d.]+)s \(map=([\d.]+)s parallel, store=([\d.]+)s serial\)\s+" + r"shuffle=([\d.]+)s\s+bud=([\d.]+)s\s+p_update=([\d.]+)s" +) +RE_MAPEFF = re.compile(r"map parallel efficiency: (\d+)% \(busy=(\d+)s") +RE_SCREEN = re.compile(r"kmer screen\s+([\d.]+)s \(\s*([\d.]+)%\)") +RE_ALIGN = re.compile(r"align total\s+([\d.]+)s \(\s*([\d.]+)%\).*?\(([\d.]+)% passed") + + +class Run: + """One log: one (rep, threads, read).""" + + def __init__(self, path): + self.path = path + self.progress = [] # (wall_s, clusters) + self.run_dada = None + self.phases = None # compare, map, store, shuffle, bud, pupdate + self.map_eff = self.busy = None + self.screen_pct = self.align_pct = self.pass_pct = None + self._parse() + + def _parse(self): + with open(self.path, errors="ignore") as fh: + for line in fh: + m = RE_PROGRESS.search(line) + if m: + self.progress.append((float(m.group(1)), float(m.group(2)))) + continue + m = RE_RUNDADA.search(line) + if m: + self.run_dada = float(m.group(1)) + continue + m = RE_PHASES.search(line) + if m: + self.phases = tuple(float(x) for x in m.groups()) + continue + m = RE_MAPEFF.search(line) + if m: + self.map_eff, self.busy = int(m.group(1)), float(m.group(2)) + continue + m = RE_SCREEN.search(line) + if m: + self.screen_pct = float(m.group(2)) + continue + m = RE_ALIGN.search(line) + if m: + self.align_pct, self.pass_pct = float(m.group(2)), float(m.group(3)) + + @property + def serial(self): + """Serial block: store + shuffle + bud + p_update (everything but the map).""" + if not self.phases: + return None + _, _, store, shuffle, bud, pupd = self.phases + return store + shuffle + bud + pupd + + def time_to(self, clusters): + """Wall seconds to reach `clusters`, linearly interpolated between points.""" + prev = (0.0, 0.0) + for wall, cl in self.progress: + if cl >= clusters: + span = cl - prev[1] + f = (clusters - prev[1]) / span if span else 0.0 + return prev[0] + f * (wall - prev[0]) + prev = (wall, cl) + return None + + +def load(sweep_dir): + runs = collections.defaultdict(list) # (read, threads) -> [Run] + pat = os.path.join(sweep_dir, "rep*", "threads*", "dada", "dada.*.log") + for path in sorted(glob.glob(pat)): + parts = path.split(os.sep) + threads = int([p for p in parts if p.startswith("threads")][0][7:]) + read = os.path.basename(path).split(".")[1] + runs[(read, threads)].append(Run(path)) + if not runs: + sys.exit(f"no logs matched {pat}") + return runs + + +def mean(xs): + xs = [x for x in xs if x is not None] + return statistics.mean(xs) if xs else None + + +def spread(xs): + """Max-min as a percentage of the mean -- the replicate noise floor. + + Printed next to every averaged figure so a difference smaller than the + replicate spread is visibly not a difference. + """ + xs = [x for x in xs if x is not None] + if len(xs) < 2: + return None + m = statistics.mean(xs) + return (max(xs) - min(xs)) / m * 100 if m else None + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("sweep_dir") + ap.add_argument("--segments", type=int, default=4, + help="number of equal cluster segments (default 4)") + ap.add_argument("--read", help="restrict to one read (R1/R2)") + args = ap.parse_args() + + runs = load(args.sweep_dir) + reads = sorted({r for r, _ in runs} if not args.read else {args.read}) + threads = sorted({t for _, t in runs}) + + for read in reads: + present = [t for t in threads if (read, t) in runs] + if not present: + continue + print(f"\n{'=' * 78}\n{os.path.basename(args.sweep_dir.rstrip('/'))} {read}\n{'=' * 78}") + + # --- 1. does it scale? ------------------------------------------- + print("\n-- run_dada by thread count (mean of reps; spread = max-min of reps)") + print(f"{'threads':>8} {'run_dada':>10} {'spread':>8} {'vs best':>9} {'reps':>5}") + vals = {t: mean([r.run_dada for r in runs[(read, t)]]) for t in present} + best = min((v, t) for t, v in vals.items() if v is not None)[1] + for t in present: + sp = spread([r.run_dada for r in runs[(read, t)]]) + d = (vals[t] - vals[best]) / vals[best] * 100 if vals[t] else 0 + mark = " <- best" if t == best else "" + spc = f"{sp:>7.1f}%" if sp is not None else " --" + print(f"{t:>8} {vals[t]:>9.1f}s {spc} {d:>+8.1f}% {len(runs[(read, t)]):>5}{mark}") + if all(spread([r.run_dada for r in runs[(read, t)]]) is None for t in present): + print(" NOTE: single replicate per arm -- differences below ~2% are not" + " distinguishable from run-to-run variation.") + + # --- 2. serial block --------------------------------------------- + print("\n-- serial block (store+shuffle+bud+p_update): invariant in seconds," + " rising in share") + print(f"{'threads':>8} {'serial':>10} {'share':>8} {'map':>10} {'busy':>10} {'map eff':>8}") + for t in present: + rs = runs[(read, t)] + s, rd = mean([r.serial for r in rs]), vals[t] + mp = mean([r.phases[1] for r in rs if r.phases]) + bz, me = mean([r.busy for r in rs]), mean([r.map_eff for r in rs]) + if None in (s, rd, mp): + continue + print(f"{t:>8} {s:>9.1f}s {s / rd * 100:>7.1f}% {mp:>9.1f}s " + f"{bz:>9.0f} {me:>7.0f}%") + print(" Amdahl cap on further map work = run_dada / serial.") + + # --- 3. where does it land? -------------------------------------- + maxcl = min(max(cl for _, cl in r.progress) + for t in present for r in runs[(read, t)] if r.progress) + edges = [round(maxcl * i / args.segments) for i in range(args.segments + 1)] + print(f"\n-- marginal seconds per cluster segment (aligned by CLUSTER, not by t=)") + head = "".join(f"{t:>9}" for t in present) + print(f"{'segment':>15}{head} {present[-1]} vs {present[0]}") + for a, b in zip(edges, edges[1:]): + row, ok = {}, True + for t in present: + xs = [] + for r in runs[(read, t)]: + ta = r.time_to(a) if a else 0.0 + tb = r.time_to(b) + if ta is None or tb is None: + ok = False + else: + xs.append(tb - ta) + row[t] = mean(xs) + ok = ok and row[t] is not None + if ok: + d = (row[present[-1]] - row[present[0]]) / row[present[0]] * 100 + print(f"{a:>7}-{b:<7}" + "".join(f"{row[t]:>9.1f}" for t in present) + + f" {d:>+7.1f}%") + + # --- 4. why? ----------------------------------------------------- + r0 = runs[(read, present[0])][0] + if r0.screen_pct is not None: + print(f"\n-- map composition at {present[0]} threads (predicts the knee)") + print(f" k-mer screen {r0.screen_pct:>5.1f}% of busy (streams k-mer" + " vectors: bandwidth-bound)") + print(f" alignment {r0.align_pct:>5.1f}% of busy (DP kernel:" + " compute-bound)") + print(f" screen pass {r0.pass_pct:>5.2f}%") + # Calibration, not a rule: four arms measured on one machine + # (EPYC 7713, 2 NUMA domains, 8 CCDs). The ORDERING has held on all + # four; the thread numbers are hardware-specific and will not travel. + print(" observed on this machine -- knee vs screen share:") + print(" 16S R1 50.8% screen -> still scaling at 96 threads") + print(" 16S R2 67.7% screen -> knee near 64") + print(" ITS2 R2 83.2% screen -> knee below 48") + print(" ITS2 R1 85.8% screen -> knee below 48") + near = min( + [(50.8, "16S R1"), (67.7, "16S R2"), (83.2, "ITS2 R2"), (85.8, "ITS2 R1")], + key=lambda x: abs(x[0] - r0.screen_pct), + ) + print(f" => this arm at {r0.screen_pct:.1f}% sits nearest {near[1]}" + f" ({near[0]}%)") + + +if __name__ == "__main__": + main() From b05c03239a8b665a6955a2685a9881791f519b9d Mon Sep 17 00:00:00 2001 From: Chris Fields Date: Sat, 29 Aug 2026 00:57:46 -0500 Subject: [PATCH 02/10] dev: correct the thread-sweep calibration and pin its reference count Two corrections, both to errors the 24-thread ITS2 arm exposed. ITS2's knee is AT 48 threads, not below it. From 48/64/96 alone the pool looked monotonically degrading and was reported as "already past the knee below 48". Adding 24 shows 48 is the peak, beating 24 by 12-17% on wall. A sweep that does not bracket the peak cannot locate it, and reading a trend off three points on one side of a maximum is how that happened. The screen/align split is itself THREAD-DEPENDENT: ITS2 R1 reads 80.0% screen at 24 threads and 85.8% at 48, because the bandwidth-bound half degrades faster under contention. The script took the composition from the lowest thread count in the sweep, so adding 24 to ITS2 silently made its number incomparable to 16S's (still measured at 48) -- in the same table that presents them as a calibration. Now pinned via --ref-threads (default 48), with a warning when that count is absent from a sweep. The calibration ordering survives both corrections: 50.8% screen -> best at 96, 67.7% -> 64, 83.2% and 85.8% -> 48. Co-Authored-By: Claude Opus 5 --- dev/analyze_thread_sweep.py | 43 +++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/dev/analyze_thread_sweep.py b/dev/analyze_thread_sweep.py index 27b10ca..09b26b3 100755 --- a/dev/analyze_thread_sweep.py +++ b/dev/analyze_thread_sweep.py @@ -25,10 +25,18 @@ screen-dominated pool saturates at low thread counts; an align-dominated one keeps scaling. Measured on four arms, the ordering is monotonic: - 16S R1 50.8% screen -> still scaling at 96 threads - 16S R2 67.7% screen -> knee near 64 - ITS2 R2 83.2% screen -> knee below 48 - ITS2 R1 85.8% screen -> knee below 48 + 16S R1 50.8% screen -> best at 96 (still scaling there) + 16S R2 67.7% screen -> best at 64 + ITS2 R2 83.2% screen -> best at 48 + ITS2 R1 85.8% screen -> best at 48 + + Two cautions learned by getting them wrong. The screen share is itself + THREAD-DEPENDENT -- ITS2 R1 reads 80.0% at 24 threads and 85.8% at 48, because + the bandwidth-bound half degrades faster under contention -- so the predictor + is only comparable when pinned to one reference count (`--ref-threads`). And + a sweep that does not bracket the peak cannot locate it: from 48/64/96 alone, + ITS2 looked monotonically degrading and was reported as "already past the knee + below 48". Adding 24 showed 48 is the peak, beating 24 by 12-17%. It predicts *within* a pool as well as across pools -- 16S's two reads differ by 17 points of screen share and have different knees -- so this is a property @@ -165,6 +173,13 @@ def main(): ap.add_argument("--segments", type=int, default=4, help="number of equal cluster segments (default 4)") ap.add_argument("--read", help="restrict to one read (R1/R2)") + ap.add_argument("--ref-threads", type=int, default=48, + help="thread count at which to report the map's screen/align " + "split (default 48). The split is itself thread-dependent " + "-- ITS2 R1 reads 80.0%% screen at 24 threads and 85.8%% at " + "48, because the bandwidth-bound half degrades faster under " + "contention -- so the predictor is only comparable across " + "arms when pinned to one reference count.") args = ap.parse_args() runs = load(args.sweep_dir) @@ -233,9 +248,14 @@ def main(): + f" {d:>+7.1f}%") # --- 4. why? ----------------------------------------------------- - r0 = runs[(read, present[0])][0] + ref = args.ref_threads if (read, args.ref_threads) in runs else present[0] + r0 = runs[(read, ref)][0] if r0.screen_pct is not None: - print(f"\n-- map composition at {present[0]} threads (predicts the knee)") + print(f"\n-- map composition at {ref} threads (predicts the knee)") + if ref != args.ref_threads: + print(f" WARNING: {args.ref_threads} threads not in this sweep; " + f"using {ref}. The split is thread-dependent, so this is NOT " + "comparable to the calibration below.") print(f" k-mer screen {r0.screen_pct:>5.1f}% of busy (streams k-mer" " vectors: bandwidth-bound)") print(f" alignment {r0.align_pct:>5.1f}% of busy (DP kernel:" @@ -244,11 +264,12 @@ def main(): # Calibration, not a rule: four arms measured on one machine # (EPYC 7713, 2 NUMA domains, 8 CCDs). The ORDERING has held on all # four; the thread numbers are hardware-specific and will not travel. - print(" observed on this machine -- knee vs screen share:") - print(" 16S R1 50.8% screen -> still scaling at 96 threads") - print(" 16S R2 67.7% screen -> knee near 64") - print(" ITS2 R2 83.2% screen -> knee below 48") - print(" ITS2 R1 85.8% screen -> knee below 48") + print(" observed on this machine -- best thread count vs screen share") + print(" (all measured at 48 threads; 2 reps per arm):") + print(" 16S R1 50.8% screen -> 96 (still scaling at 96)") + print(" 16S R2 67.7% screen -> 64") + print(" ITS2 R2 83.2% screen -> 48") + print(" ITS2 R1 85.8% screen -> 48") near = min( [(50.8, "16S R1"), (67.7, "16S R2"), (83.2, "ITS2 R2"), (85.8, "ITS2 R1")], key=lambda x: abs(x[0] - r0.screen_pct), From 79ad37b1a6e1c109b776f64b62e3a2de32c18e80 Mon Sep 17 00:00:00 2001 From: Chris Fields Date: Sat, 29 Aug 2026 11:17:05 -0500 Subject: [PATCH 03/10] dev: script the one-wide-job vs two-narrow-jobs test (#152) The #152 sweep says every arm buys wall time at ~4x the core-seconds -- 16S R1 goes 1072.5s at 24 threads to 594.3s at 96, 45% parallel efficiency at the top -- which implies "pack jobs, don't scale threads". Nobody has measured two jobs sharing a node, so that is currently an inference and not a result. The two jobs would contend for exactly the memory bandwidth the k-mer screen is already limited by. Three arms. The interesting one is NUMA policy, because it may invert an existing finding: dev/numa_pin.sh documents that binding to one domain is the wrong choice for a single job, costing the parallel map 21% by forcing every thread through one node's controllers. With two concurrent jobs on a 2-domain node that reasoning may reverse -- each job gets private controllers and neither disturbs the other, where interleaving both has them share every controller. Not predictable from the single-job result, so both are measured. Runs the SAME read twice rather than R1+R2. Pairing 16S R1 (713.8s) with R2 (767.8s) at 48 threads lets the shorter job finish first, leaving the longer one's tail uncontended -- which flatters the pair and blurs the contention being measured. Times via epoch timestamps rather than /usr/bin/time, which is not installed on the benchmark node. Drops the split arm automatically on a single-domain machine rather than reporting a meaningless comparison, and reminds the operator that outputs must be byte-identical across arms -- concurrency and NUMA policy change timing only. Co-Authored-By: Claude Opus 5 --- dev/run_concurrency_test.sh | 139 ++++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100755 dev/run_concurrency_test.sh diff --git a/dev/run_concurrency_test.sh b/dev/run_concurrency_test.sh new file mode 100755 index 0000000..d6549b6 --- /dev/null +++ b/dev/run_concurrency_test.sh @@ -0,0 +1,139 @@ +#!/usr/bin/env bash +# run_concurrency_test.sh — is it better to run one wide job or two narrow ones? +# --------------------------------------------------------------------------- +# WHY: the #152 thread sweep found that every arm measured buys wall time at +# roughly 4x the core-seconds — 16S R1 goes 1072.5s at 24 threads to 594.3s at +# 96, i.e. 45% parallel efficiency at the top. That says "pack jobs, don't scale +# threads". But nobody has measured two jobs sharing a node, and they would +# contend for exactly the memory bandwidth the k-mer screen is already limited +# by. The packing recommendation is currently an inference, not a result. +# +# THE NUMA QUESTION, which is the interesting half: dev/numa_pin.sh documents +# that binding to a single domain is the WRONG choice for one job — it forces +# every thread through one node's memory controllers and costs the parallel map +# 21%. With two concurrent jobs on a 2-domain node that reasoning may invert: +# each job gets a private set of controllers and neither disturbs the other. +# Interleaving both jobs, by contrast, has them share every controller. Which +# wins is not predictable from the single-job result, so both are measured. +# +# WHY THE SAME READ TWICE: pairing R1 with R2 (713.8s vs 767.8s at 48 threads +# on 16S) lets the shorter job finish first, leaving the longer one's tail +# running uncontended — which flatters the pair and blurs the contention it is +# meant to measure. Two copies of one read contend end to end. +# +# ARMS +# solo one job, N threads, --interleave=all (the #152 baseline) +# both two jobs, N threads each, both interleaved (naive packing) +# split two jobs, N threads each, each bound to its own NUMA domain +# +# Contention penalty = (arm wall) / (solo wall) - 1, per job. Ideal is 0%: +# two jobs finishing in the time one takes means the node was not the limit. +# +# Usage: +# dev/run_concurrency_test.sh \ +# --bin ./target/release-native/dada2-rs \ +# --error-model err.json \ +# --out-root tmp/issue-152/concurrency \ +# --threads 48 \ +# derep/*.json +set -euo pipefail + +BIN=./target/release-native/dada2-rs +ERRMODEL= +OUT_ROOT=concurrency-test +THREADS=48 +ARMS="solo both split" + +while [[ $# -gt 0 ]]; do + case "$1" in + --bin) BIN=$2; shift 2 ;; + --error-model) ERRMODEL=$2; shift 2 ;; + --out-root) OUT_ROOT=$2; shift 2 ;; + --threads) THREADS=$2; shift 2 ;; + --arms) ARMS=$2; shift 2 ;; + --) shift; break ;; + -*) echo "unknown option: $1" >&2; exit 2 ;; + *) break ;; + esac +done +INPUTS=("$@") + +[[ -n "$ERRMODEL" ]] || { echo "--error-model is required" >&2; exit 2; } +[[ ${#INPUTS[@]} -gt 0 ]] || { echo "no input derep files given" >&2; exit 2; } +command -v numactl >/dev/null || { echo "numactl not found" >&2; exit 2; } + +NODES=$(numactl --hardware | awk '/^available:/ {print $2}') +echo "== node has ${NODES} NUMA domain(s); threads/job = ${THREADS}" +if [[ "$NODES" -lt 2 ]]; then + echo "== only one NUMA domain: the 'split' arm is meaningless here, dropping it" + ARMS=${ARMS//split/} +fi +mkdir -p "$OUT_ROOT" + +# One job. $1 label, $2 numactl args. Writes