perf: b_p_update is memory-bound on a scattered 160-byte Raw read, not compute-bound (94% of repricings never reach calc_pA)

[cjfields]

Why now

The #152 thread sweep closes the thread-scaling axis: on soil ITS2, 48 → 96
threads makes run_dada slower (194 → 204 s R1, 255 → 263 s R2), with
worker-busy time rising +87% for identical work while the map wall gets worse.
The map is bandwidth-saturated; adding cores subtracts.

What is left is the serial block. At 48 threads on ITS2 R1 that is store 16.9 +
shuffle 41.8 + p_update 20.0 + bud 0.4 = 79.0 s of 193.5 s, 41% of
run_dada — an Amdahl cap of 2.45× even if the map became free.

Every other serial phase has been through two or three rounds of work (#132,
#136, #139, #143, #147). b_p_update has had none.

The measurement

pool / read	p_update	run_dada	share	repricings	ns each
ITS2 R1 (48t)	20.0 s	193.5 s	10.3%	355.9 M	56
ITS2 R2 (48t)	31.1 s	257.2 s	12.1%	565.4 M	55
soil 16S R1 (64t)	29.7 s	654.1 s	4.5%	486.3 M	61

~56 ns per repricing is compute, not memory. For contrast the store scan —
after #148 — runs at 4.9 ns/raw-visit. get_pA is over 11× more expensive per
item than a cache-line-bound scan, which says the cost is the p-value
arithmetic rather than the traversal.

That matters for the lever: unlike the map, this is not competing for the
bandwidth that is already saturated, so it should actually scale with threads.

Structure, and the one hazard

b_p_update (src/pval.rs:45) loops over clusters; for each with update_e
set, it walks the members and, per raw:

1. computes get_pA(...) and writes b.raws[raw_idx].p — the expensive part;
2. selects bud_min / bud_min_prior, a per-cluster reduction whose tie-break
   is explicitly first-in-scan (min p, then max reads, then lowest r) and
   documented as matching b_bud's scan order.

Each raw belongs to exactly one cluster, so the p writes are disjoint.

Do not parallelise the outer loop over clusters. It looks like the obvious
split and it fails exactly where it is needed: early in the run there are very
few clusters and cluster 0 holds nearly every raw, so outer-loop parallelism
gives ~no parallelism in precisely the early, serial-heavy phase #152 identified
as the worst. It would optimise the part of the run that is already fine.

Do not parallelise the bud selection. Its tie-break is order-dependent, and
reordering it changes which raw buds — that is an inference change, not a perf
change.

Proposed shape

Split the two halves: compute get_pA for the round's raws in parallel,
then run the existing selection serially over the results. Selection is
comparisons on values already in hand, so it should cost a few ns/item against
the 56 ns it currently guards. get_pA is independent per raw, so there is no
reduction-order or floating-point-associativity concern.

Rough sizing, to be treated as a hypothesis and not a plan: if get_pA is ~80%
of the 56 ns, then at 48 threads 20.0 s → ~4 s, i.e. ~8% of run_dada on
ITS2. Smaller on 16S (4.5% share), so this is a pool-dependent lever and ITS2 is
the place to measure it.

Gates

• Byte-identical output, on both pools and both reads, two reps. Keeping
  selection serial is what makes that achievable; if a design needs to relax it,
  that design is out of scope here.
• Measure at 48 threads, not 64. Pooled run_dada is not uniform: occupancy climbs 22 -> 42 cores across a run, and the serial fraction is concentrated in the early rounds #152 shows 64 is already past the knee on
  ITS2, and a lever measured above the knee is measured in the wrong regime.
• Re-check the phase share per window (DADA2RS_PROGRESS_SECS), not just as
  a run total — Pooled run_dada is not uniform: occupancy climbs 22 -> 42 cores across a run, and the serial fraction is concentrated in the early rounds #152 exists because run totals hid a 2× swing in occupancy.

Related

The other open direction from the same analysis is reducing what the map
touches rather than how many threads it gets — 2-bit sequence encoding (#28)
and replacing the k-mer screen with minimizers. Those move the knee to the
right and re-enable thread scaling, which is higher leverage than this issue but
considerably more work. This one is small, self-contained, and untouched.

[cjfields]

Measured: the premise is wrong — this is a layout problem, not a parallelism one

Soil ITS2, both reads, 48 threads, numactl --interleave=all, exclusive node,
branch perf/pupdate-attribution-154.

	R1	R2
p_update	20.63 s	32.00 s
repricings	355,895,585	565,379,573
ns / repricing	58.0	56.6
reached calc_pA	20,332,845 (5.7%)	26,275,552 (4.6%)
exited as singleton	335,424,353 (94.2%)	538,850,126 (95.3%)
exited as centre	138,387 (0.0%)	250,751 (0.0%)
exited on zero lambda	0	3,144 (0.0%)

Over 94% of repricings never reach calc_pA. This issue proposed
parallelising it because the regularised incomplete gamma is expensive
arithmetic that would scale with threads and, unlike the map, is not competing
for saturated bandwidth. That reasoning was sound; the premise underneath it was
not.

The cost model, and why it is robust

At 58.0 ns per repricing and ~52 ns per actual calc_pA call
(examples/pval_cost.rs), the arithmetic contributes ~3 ns of the 58 — about
5%. Parallelising it perfectly would save ~1 s of a 200 s run_dada: 0.5%.

That estimate comes from a benchmark that under-models, so here is the
sensitivity — the conclusion does not depend on it:

if calc_pA really costs	its share of p_update
52 ns/call (measured)	5%
130 ns/call (2.5×)	13%
260 ns/call (5×)	25%
>500 ns/call	>50%

An incomplete gamma is nowhere near 500 ns. calc_pA cannot be the
bottleneck at 5% of calls, whatever it costs.

What the 58 ns actually buys

~55 of it is the early-exit path: two comparisons, plus a scattered read of
four fields (reads, prior, comp.lambda, comp.hamming) out of a 160-byte
Raw, plus the write-back of p to the same line. That is a cache miss per
repricing, paid whether or not any arithmetic follows.

It is #147's pattern, where the
same shape turned out to be 83–87% of b_compare's serial store and the fix was
a dense parallel array rather than threads.

Proposed retarget

Hoist the hot fields into dense arrays on B, parallel to raw_cluster and
e_minmax. They pack into ~24 bytes against Raw's 160, so on ITS2 R1 the
working set goes 132 MB → ~20 MB, which fits the 7713's 32 MB per-CCD L3 —
the change is not only fewer bytes per access but potentially a resident hot
array instead of a guaranteed miss.

Two constraints on any implementation:

• Do not reorder Bi::raws. Sorting member indices would make the traversal
  far more cache-friendly, and it changes results: the bud tie-break is
  first-in-scan on position r, so reordering members changes which raw buds.
  That is an inference change, not a perf change.
• Byte-identical output, both pools, both reads.

Expected payoff is bounded by the phase itself: 10–12% of run_dada on ITS2,
4.5–6% on 16S. So this is a mid-single-digit change at best, not another #148 —
worth stating up front so it is not oversold.

Two incidental findings

• The greedy lock pass is a non-issue. It walked 1,337,664 raws against
  355,895,585 repricings — 0.4%. I had flagged it as a possible uncosted
  second traversal; it is not.
• Repricing is concentrated in a handful of huge clusters. Only 13,092 of
  6,246,344 cluster visits were dirty (0.21%) — ~3.7 dirty clusters per round —
  yet 100,735 raws were repriced per round, roughly 27,000 per dirty cluster.
  This is cluster 0 and its neighbours, and it independently confirms that
  parallelising the outer loop over clusters would have achieved nothing:
  early in a run there is almost nothing to spread.

Note on reading exit singleton

The 94% is work-weighted, not a biological singleton fraction: a singleton in
a frequently-dirtied cluster is counted on every round. It answers exactly the
question here — what fraction of this phase's work is early-exit — but should not
be quoted as a property of the pool.

Outstanding

Equivalent 16S run, for completeness. It would have to show full calc_pA above
~50% to change the conclusion, which is a long way from 5%.

[cjfields]

16S confirms it, and adds a second independent tell

Soil 16S, both reads, 48 threads, numactl --interleave=all, exclusive node.

	ITS2 R1	ITS2 R2	16S R1	16S R2
p_update	20.63 s	32.00 s	31.80 s	47.97 s
repricings	355.9 M	565.4 M	486.3 M	828.3 M
full calc_pA	5.7%	4.6%	7.9%	4.8%
exit singleton	94.2%	95.3%	91.9%	95.1%
ns / repricing	58.0	56.6	65.4	57.9

full calc_pA is 4.6–7.9% on four arms across two pools with very different
character. Nothing here is close to the ~50% that would be needed to make
parallelising the arithmetic the right lever.

The bottom row is a second, independent line of evidence. Per-repricing cost
sits in a 56.6–65.4 ns band across pools differing 1.5× in nraw and 2.3× in
repricing count. A per-unit cost that barely moves with the working set is the
same signature that identified the serial store in
#147 — flat 18.4–19.6
ns/raw-visit across a 5× range in pool size — and it means the cost is a fixed
per-access price, i.e. a cache miss, not work that scales with the data.

It also rules the arithmetic out directly rather than by estimate: 16S R1 has the
highest calc fraction (7.9%) and the highest ns/repricing, but the
3.3-point calc difference against ITS2 R2 accounts for only ~1.7 ns of the 8.8 ns
gap between them.

The two incidental findings hold on 16S as well:

• Greedy lock pass: 5.65 M raws walked against 486.3 M repricings — 1.2%.
  Still a non-issue.
• Dirty clusters: 136,775 of 47,059,550 visited (0.29%), ~14 per round,
  against 50,134 raws repriced per round. Repricing stays concentrated in a few
  very large clusters, so outer-loop parallelism over clusters would have found
  almost nothing to spread.

Conclusion unchanged, now on two pools: this is a layout problem. The
retarget stands — hoist the hot fields into dense arrays parallel to
raw_cluster and e_minmax, without reordering Bi::raws.