Skip to content

test: A/B bench harness that drops rounds it could not measure honestly - #171

Open
BenjaminDEMAILLE wants to merge 1 commit into
scverse:mainfrom
BenjaminDEMAILLE:bd/bench-ab
Open

test: A/B bench harness that drops rounds it could not measure honestly#171
BenjaminDEMAILLE wants to merge 1 commit into
scverse:mainfrom
BenjaminDEMAILLE:bd/bench-ab

Conversation

@BenjaminDEMAILLE

Copy link
Copy Markdown
Contributor

test/bench_ab.sh: interleaved A/B of two binaries, which drops the rounds it could not measure honestly instead of averaging over them.

Why

Four traps, each of which I fell into at least once while measuring this aligner over one session. Three are now enforced by the script; the fourth is documented in its header because no script can enforce it.

Idle was checked once, at the start. A previous harness sampled CPU idle before the first run and never again. Three separate measurements then drifted below the threshold mid-series and it kept going, leaving the contaminated rounds to be spotted by eye. This samples idle before and after every run, drops the round if any of the four samples falls short, and prints the count:

round1  binA=3.33s binB=3.42s  idle>=84.92%
round3  binA=3.25s binB=3.22s  DROPPED (idle fell to 71.4%)

kept 4 of 6 rounds (2 dropped)

A median over four surviving rounds now cannot be read as a median over six. No silent truncation.

Load average is the wrong signal. It is an exponential average over minutes, so it stays high long after the offending job has gone. An earlier guard refused to measure at 2.24 on a machine whose sixteen cores were all free. The check is on CPU idle percentage, which answers the question actually being asked.

A median can be smaller than the noise. The script prints the spread within each side next to the difference between the medians, and when the difference is smaller it says so:

median  binA=3.305s  binB=3.25s  delta=-1.7%
spread  binA=0.08s   binB=0.20s

the difference between the medians is smaller than the spread within a side.
that is not a result: report it as unmeasured rather than as a speedup.

That message is not hypothetical. Two changes I measured this week looked like wins on medians alone and were inside the spread; both are now described as unmeasured in their PRs (#169, #170).

Both sides must see identical argv. @PG CL: records the command line verbatim, so running ./old against ./new, or writing to differently-named output directories, makes the BAM differ in the header and in its compressed size. That cost me two false "output is not neutral" alarms. Each side is copied to rustar-aligner inside its own directory and run with --outFileNamePrefix ./.

Timing a total hides the part that changed. Only in the header, since it is a judgement the script cannot make: BAM writing is 1-4% of a yeast run, so a total dominated by alignment cannot resolve a change to the writer at any reachable scale. Run one configuration with --outSAMtype None alongside and read the difference. #164 was closed because that difference, once isolated, showed the change was a regression rather than the improvement its description claimed.

Usage

test/bench_ab.sh <binA> <binB> <genomeDir> <reads.fastq> [threads] [rounds]

BENCH_MIN_IDLE=88   percent CPU idle a round must hold to count
BENCH_MODE="None"   --outSAMtype words (quote the two-word forms)
BENCH_ARGS=""       extra flags passed to both sides

It exits non-zero if every round was dropped, rather than printing a median over nothing.

Scope

Test tooling only: no source change, no dependency, nothing user-facing, so no CHANGELOG entry. The numbers in #167 and #170 were produced with the same method; this makes the method reproducible rather than living in my shell history.

Happy to add a pointer to it from CONTRIBUTING.md's verification section if you want one there, but I did not want to edit that file uninvited.

Related: #168 has the allocation measurements taken with this method, including the result that removing 7.6% of the allocations moved nothing.

Every perf measurement on this project needs the same three guards, and
I got each of them wrong at least once in a single session:

An earlier harness sampled CPU idle once, before the first run. Three
separate measurements then drifted below the threshold mid-series and it
kept going, leaving me to spot the contaminated rounds by eye in the
output. This samples idle before and after every run, drops a round if
any of the four samples falls short, and prints how many it dropped, so
a median over four surviving rounds cannot be mistaken for a median over
six.

The check is on CPU idle rather than load average, because load average
is an exponential average over minutes: it refused to measure at 2.24 on
a machine whose cores were all free.

It reports the spread within each side next to the difference between
the medians, and says so in as many words when the difference is
smaller. Two changes I measured looked like wins on medians alone and
were inside the spread.

Both sides run as ./rustar-aligner with --outFileNamePrefix ./ from
inside their own directory, because the @pg CL: line records argv
verbatim: running ./old against ./new is enough to make the output
differ, which cost me two false "output is not neutral" alarms.

The header documents the fourth trap, which no script can enforce:
timing a total hides the part that changed. BAM writing is 1-4% of a
yeast run, so a total dominated by alignment cannot resolve a change to
the writer. Run the None configuration alongside and read the
difference.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Psy-Fer added a commit that referenced this pull request Aug 5, 2026
* fix: --runThreadN 1 ran on every core, not on one

The rayon global pool was configured only when `--runThreadN` was
greater than 1. Skipping the build at 1 does not give one thread: it
leaves rayon's default, which is one worker per logical core. So
`--runThreadN 1` ran the whole machine.

Measured on 200k reads, before: 2.44 s wall at **1300% CPU**. After:
26.23 s wall at 100% CPU. The old figure was not a fast single-threaded
run, it was a sixteen-way run wearing the wrong flag.

This matters beyond the flag reading falsely. A scheduler or a container
given one CPU gets sixteen worker threads; on a shared machine the run
oversubscribes every other job; and with a thread-caching allocator each
of those threads keeps its own heap, which is the very cost the comment
above this code says the pool sizing exists to avoid.

It also means the project's thread-invariance checks were weaker than
they read: the `--runThreadN 1` leg was not a one-thread leg. Verified
now that it is one: records are byte-identical between 1 and 8 threads
on 200k real reads, and byte-identical to the previous binary's output
at `--runThreadN 1`. Only the `@PG` `CL:` line differs between thread
counts, because it records the command line.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the --runThreadN 1 fix

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* refactor: drop Transcript::read_seq, which nothing ever reads

The field was filled with a full copy of the read at every finalised
transcript and assigned again in four more places, and no code in the
crate reads it. Not "reads it rarely": the compiler was asked, and after
deleting the field every one of the 97 resulting errors is a struct
literal or an assignment. There is no read site, in `src` or in tests.

Measured on 200k real reads at 8 threads: about 50 transcripts are
finalised per read, so removing it takes **15.1 million allocations off a
200 million total** (7.6%) and 3.87 GB of copying.

Wall clock does not move, and that is worth recording rather than hiding:
six interleaved rounds at 87-92% CPU idle give medians 20.90 s against
20.65 s with the direction mixed, inside the run-to-run spread. mimalloc
is fast enough that seventy-five small allocations per read do not
surface. The reason to remove it is that it is dead weight, not that it
is slow.

That number also calibrates #168 downward: if removing 7.6% of the
allocations changes nothing measurable, the rest of the allocation
programme is unlikely to be worth a new dependency.

`Transcript::read_seq` is `pub`, so this is an API removal and needs
sign-off. Nothing outside the crate can be relying on its contents being
meaningful, though, since it is only ever written.

Output-neutral: SAM byte-identical on 200k real reads.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the Transcript::read_seq removal

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* test: an A/B harness that refuses to average over a busy machine

Every perf measurement on this project needs the same three guards, and
I got each of them wrong at least once in a single session:

An earlier harness sampled CPU idle once, before the first run. Three
separate measurements then drifted below the threshold mid-series and it
kept going, leaving me to spot the contaminated rounds by eye in the
output. This samples idle before and after every run, drops a round if
any of the four samples falls short, and prints how many it dropped, so
a median over four surviving rounds cannot be mistaken for a median over
six.

The check is on CPU idle rather than load average, because load average
is an exponential average over minutes: it refused to measure at 2.24 on
a machine whose cores were all free.

It reports the spread within each side next to the difference between
the medians, and says so in as many words when the difference is
smaller. Two changes I measured looked like wins on medians alone and
were inside the spread.

Both sides run as ./rustar-aligner with --outFileNamePrefix ./ from
inside their own directory, because the @pg CL: line records argv
verbatim: running ./old against ./new is enough to make the output
differ, which cost me two false "output is not neutral" alarms.

The header documents the fourth trap, which no script can enforce:
timing a total hides the part that changed. BAM writing is 1-4% of a
yeast run, so a total dominated by alignment cannot resolve a change to
the writer. Run the None configuration alongside and read the
difference.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* feat(solo): --soloCellReadStats CB writes CellReads.stats

One row per cell barcode with the fifteen counters STARsolo reports: how the
barcode matched, whether the read mapped to one locus or several, whether it
landed on a feature, where in the gene and on which strand, whether it was
mitochondrial, and whether it reached the matrix. The per-cell UMI and gene
totals come from the raw matrix rather than from the read counters, so they
agree with what the matrix says by construction.

Reads whose barcode never resolved are summed into a single `CBnotInPasslist`
row instead of being dropped. That row is the reason the file is useful: it is
the difference between "these cells look thin" and "most of the input never
reached a cell at all".

The region columns split by strand — an antisense read counts under `exonicAS`
or `intronicAS`, never under `exonic` or `intronic`.

`--genomeChrSetMitochondrial` names the chromosomes behind the `mito` column.
Without it the column is zero throughout, which is honest: no chromosome was
declared mitochondrial.

D24 comes with it. STAR emits these rows by walking a libc++ `unordered_map`,
which at these sizes is the reverse of each barcode's first appearance. That is
reproduced, including across threads: the per-read accumulator merges in read
order, so a threaded run writes the same file as a serial one. It stops being
reproducible past the point where libc++ rehashes, since the order then depends
on the bucket count. The values never differ, only which line they sit on.
Recorded in docs-old/dev/divergences.md.

* docs: record the CellReads.stats row order in DIVERGENCE.md

Section 3.2, in the format CONTRIBUTING.md asks for.

* docs(changelog): keep only this PR's entry

CONTRIBUTING.md requires the description to match the code; the entries for the
other themes split out of #152 belong to their own PRs.

* refactor(solo): drop CellReadStats::merge, which nothing calls

Reads are folded in under a mutex, so there are no per-thread partials to merge;
the function was reachable only from its own test. CONTRIBUTING.md rules out
shipping a function no production path reaches, and the PR description claimed
its test as evidence of thread-safety that the mutex actually provides.

* feat(solo): --runMode soloCellFiltering

Cell-calls an existing raw count matrix without aligning anything, taking the
raw directory and an output prefix as STAR does:
`--runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix`.

Cell calling is a decision about a matrix, not about reads. Re-calling with
different `--soloCellFilter` parameters should not mean re-aligning 400 million
reads, and a matrix produced by another tool should be callable too.

The matrix is streamed into the same temp-body form the align path builds, so
`called_cells` and `emptydrops_called` are the identical code here and there
rather than a second implementation free to drift. Counts are rounded on the way
in: a multimapper matrix carries real values, and the filters work on UMI
totals.

`--runMode` becomes a token list, because that is what STAR's is: the mode
followed by its arguments. The mode itself is now validated rather than falling
back to `alignReads`, so a typo is refused instead of quietly running something
else.

The standalone `emptydrops` binary still exists and still carries its own copy
of the algorithm, which no longer matches this one. Removing it means moving
`test/solo_genefull_compare.py` and `test/solo_genefull_h5_compare.py` to the
new mode first, so it is left alone here rather than broken.

---------

Co-authored-by: Benjamin Demaille <benjamin.demaille@icloud.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants