diff --git a/dev/analyze_thread_sweep.py b/dev/analyze_thread_sweep.py new file mode 100755 index 0000000..09b26b3 --- /dev/null +++ b/dev/analyze_thread_sweep.py @@ -0,0 +1,282 @@ +#!/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 -> 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 + 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)") + 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) + 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? ----------------------------------------------------- + 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 {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:" + " 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 -- 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), + ) + print(f" => this arm at {r0.screen_pct:.1f}% sits nearest {near[1]}" + f" ({near[0]}%)") + + +if __name__ == "__main__": + main() diff --git a/dev/benchmark/bench_pooled.py b/dev/benchmark/bench_pooled.py index 8e48d11..302a610 100755 --- a/dev/benchmark/bench_pooled.py +++ b/dev/benchmark/bench_pooled.py @@ -81,6 +81,7 @@ import json import os import re +import shutil import subprocess import sys import time @@ -223,9 +224,99 @@ def capture_version(bin_, outdir): return ver + +# --------------------------------------------------------------------------- +# NUMA policy (issue #152) +# --------------------------------------------------------------------------- +# +# Memory placement is worth 25-30% of `run_dada` on this hardware, so a +# benchmark that pins one stack and not the other is not a benchmark. The policy +# is therefore applied in `numa_wrap`, which BOTH stacks' launches route through +# (`run_step` and `run_phase_concurrent`) — rather than at each call site, where +# an arm could be missed and the speedup silently inflated by ~25%. +# +# R DADA2 threads the same algorithm through `multithread=`, with the same k-mer +# screen and banded NW, so there is no reason to expect it responds differently. +# Whether it gains as much is the interesting question and is measurable only +# with both arms treated alike. +# +# The EFFECTIVE policy is recorded in summary.csv, not the requested one: a run +# that asked for `bind` on a machine without numactl must not be filed as bound. + +NUMA_POLICY = "none" # what the operator asked for +NUMA_EFFECTIVE = "none" # what actually happened +NUMA_PREFIX = [] + + +def numa_init(policy, threads): + """Resolve the NUMA policy once, warning loudly on any downgrade. + + Never fails the run: a missing `numactl` or a single-domain machine + degrades to unpinned execution, which is correct behaviour but must be + visible — an unpinned run recorded as bound would be worse than no feature. + """ + global NUMA_POLICY, NUMA_EFFECTIVE, NUMA_PREFIX + NUMA_POLICY = policy + NUMA_EFFECTIVE = "none" + NUMA_PREFIX = [] + if policy == "none": + return + + if shutil.which("numactl") is None: + print(f" [numa] WARNING: --numa {policy} requested but numactl is not " + "installed; running unpinned. Recorded as numa=none.", file=sys.stderr) + return + try: + hw = subprocess.run(["numactl", "--hardware"], capture_output=True, + text=True, check=True).stdout + except (subprocess.SubprocessError, OSError) as e: + print(f" [numa] WARNING: numactl --hardware failed ({e}); running " + "unpinned. Recorded as numa=none.", file=sys.stderr) + return + + m = re.search(r"^available:\s+(\d+)", hw, re.M) + nodes = int(m.group(1)) if m else 1 + if nodes < 2: + print(f" [numa] single NUMA domain; --numa {policy} has nothing to do. " + "Recorded as numa=none.", file=sys.stderr) + return + + if policy == "interleave": + NUMA_PREFIX = ["numactl", "--interleave=all"] + NUMA_EFFECTIVE = "interleave" + elif policy == "bind": + # Binding is only the faster policy while the thread count FITS inside a + # domain. At or above a domain's core count it saturates that domain's + # controllers and loses to interleaving — which is how dev/numa_pin.sh + # came to recommend the opposite, from a 64-thread measurement on a + # 64-core domain. Warn rather than silently choose for the operator. + cores = re.findall(r"^node \d+ cpus:(.*)$", hw, re.M) + per_domain = len(cores[0].split()) if cores else 0 + if per_domain and threads > per_domain: + print(f" [numa] WARNING: --threads {threads} exceeds the {per_domain} " + f"CPUs in one domain; binding cannot honour it. Falling back to " + "interleave. Recorded as numa=interleave.", file=sys.stderr) + NUMA_PREFIX = ["numactl", "--interleave=all"] + NUMA_EFFECTIVE = "interleave" + else: + NUMA_PREFIX = ["numactl", "--cpunodebind=0", "--membind=0"] + NUMA_EFFECTIVE = "bind" + if per_domain and threads > per_domain * 0.9: + print(f" [numa] note: {threads} threads nearly fills a " + f"{per_domain}-CPU domain; the binding advantage shrinks as " + "the domain saturates.", file=sys.stderr) + print(f" [numa] policy={NUMA_EFFECTIVE} ({' '.join(NUMA_PREFIX)}) " + "— applied to BOTH stacks", file=sys.stderr) + + +def numa_wrap(cmd): + """Prefix a command with the resolved NUMA policy. Both stacks route here.""" + return [*NUMA_PREFIX, *cmd] if NUMA_PREFIX else cmd + + def run_step(name, cmd, logf, results, append_log=False): """Run cmd as one process; record (name, wall_s, maxrss_kb, rc). Returns rc.""" - cmd = maybe_align_backend(maybe_verbose(cmd)) + cmd = numa_wrap(maybe_align_backend(maybe_verbose(cmd))) print(f" ==> {name}: {' '.join(str(c) for c in cmd)}", flush=True) start = time.time() with open(logf, "ab" if append_log else "wb") as lf: @@ -255,7 +346,7 @@ def run_phase_concurrent(name, jobs, results, max_workers): print(f" ==> {name}: {len(jobs)} samples, up to {max_workers} concurrent", flush=True) def one(cmd, logf): - cmd = maybe_align_backend(maybe_verbose(cmd)) + cmd = numa_wrap(maybe_align_backend(maybe_verbose(cmd))) with open(logf, "wb") as lf: proc = subprocess.Popen([str(c) for c in cmd], stdout=subprocess.DEVNULL, stderr=lf) @@ -884,6 +975,14 @@ def main(): p.add_argument("input", help="directory of raw FASTQ files") p.add_argument("--outdir", default="bench_pooled_out") p.add_argument("--threads", type=int, default=1) + p.add_argument("--numa", choices=["none", "interleave", "bind"], default="none", + help="memory placement, applied to BOTH stacks (default: none, " + "matching every historical run). 'bind' pins CPUs and memory " + "to one NUMA domain and is worth 25-30%% of run_dada when the " + "thread count fits inside a domain; 'interleave' spreads pages " + "and is the more reproducible policy. Degrades to unpinned, " + "with a warning, where numactl is unavailable or the machine " + "has one domain — the EFFECTIVE policy is what gets recorded.") p.add_argument("--nbases", type=float, default=1e8) p.add_argument("--pool", choices=["true", "false", "pseudo"], default="true", help="denoising mode: 'true' = pooled (R pool=TRUE / dada-pooled, " @@ -1026,6 +1125,7 @@ def main(): INFER_KDIST = args.kdist_cutoff LEARN_KDIST = args.learn_kdist_cutoff LOESS_PRESET = args.loess_preset + numa_init(args.numa, args.threads) if args.reestimate_err_between_rounds and args.pool != "pseudo": p_err = "--reestimate-err-between-rounds applies to --pool pseudo only" raise SystemExit(p_err) @@ -1100,13 +1200,23 @@ def main(): learn_c = LEARN_KDIST if LEARN_KDIST is not None else infer_c kdmode = (f", kdist learn={learn_c}/infer={infer_c} (decoupled)" if learn_c != infer_c else f", kdist={infer_c}") + nmode = f", numa={NUMA_EFFECTIVE}" + if NUMA_EFFECTIVE != NUMA_POLICY: + nmode += f" (requested {NUMA_POLICY})" print(f"BENCHMARK SUMMARY — {args.platform}, {mode} denoise, " - f"{args.threads} thread(s){bemode}{kdmode}{rmode}") + f"{args.threads} thread(s){bemode}{kdmode}{rmode}{nmode}") print("=" * 56) print(f" cores = CPU/wall (effective cores; ideal ≈ {args.threads} for an " "in-process step, ≈ min(#samples, threads) for fanned steps)") csv_path = outdir / "summary.csv" with open(csv_path, "w") as cf: + # The effective policy is a header comment rather than a column: it is + # a property of the whole run, and a run gathered under one policy is + # not comparable to one gathered under another (25-30% on this + # hardware). Recording the REQUESTED policy too makes a silent + # downgrade — missing numactl, single domain — visible in the artifact. + cf.write(f"# numa_effective={NUMA_EFFECTIVE} numa_requested={NUMA_POLICY} " + f"threads={args.threads}\n") cf.write("stack,step,wall_s,cpu_s,cores,maxrss_kb\n") rs = print_stack("dada2-rs", rust_results, cf) if rust_results else None rr_split = print_stack("R-split", r_split_results, cf) if r_split_results else None diff --git a/dev/numa_pin.sh b/dev/numa_pin.sh index aaa421b..1f6dbca 100755 --- a/dev/numa_pin.sh +++ b/dev/numa_pin.sh @@ -28,6 +28,29 @@ # also sets memory policy only, never CPU affinity, so unlike binding it cannot # silently oversubscribe a thread count sized for the whole node. # +# THREAD-COUNT CONDITIONAL -- READ THIS BEFORE APPLYING THE ABOVE (#152). +# The 21% figure was measured at 64 threads, where 64 threads bound to a 64-core +# domain saturates that domain's controllers. At 48 threads, which leaves +# headroom, binding is the FASTER policy by a wide margin -- and on both pools +# tested, not just one: +# +# run_dada, 48 threads, interleave -> cpunodebind ITS2 R1 ITS2 R2 16S R1 +# one job, --interleave=all 197.8s 259.7s 708.1s +# one job, --cpunodebind=0 --membind=0 143.9s 185.2s 528.1s +# -27% -29% -25% +# +# Every phase gains, including the single-threaded ones (16S: map -27%, store +# -22%, shuffle -21%, p-update -31%), which is memory *latency* rather than +# bandwidth. Outputs byte-identical across all arms. +# +# So: do not read "interleave beats bind" as general. It holds when the thread +# count fills or exceeds a domain; it reverses when the count fits inside one. +# Benchmark numbers gathered under this helper are therefore CONSERVATIVE in +# absolute terms at sub-domain thread counts -- A/B deltas are unaffected, since +# both arms always share the policy. Whether binding is as *reproducible* as +# interleaving is not yet measured, which is the only reason the default here +# has not changed. See docs/tuning-for-your-data.md. +# # NOTE this is a *measurement* tool. Interleaving is not a speed recommendation: # at 64 threads it matches default placement on mean run_dada (372s vs 372s). # What it buys is predictability. diff --git a/dev/run_concurrency_test.sh b/dev/run_concurrency_test.sh new file mode 100755 index 0000000..97640ff --- /dev/null +++ b/dev/run_concurrency_test.sh @@ -0,0 +1,155 @@ +#!/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) +# sbind one job, N threads, bound to ONE NUMA domain (locality, no sharing) +# 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 sbind 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