Skip to content

bathsearch: search every query in one pass over the target - #38

Open
kevinmoran1988 wants to merge 11 commits into
TravisWheelerLab:mainfrom
kevinmoran1988:bathsearch-one-pass
Open

kevinmoran1988 wants to merge 11 commits into
TravisWheelerLab:mainfrom
kevinmoran1988:bathsearch-one-pass

Conversation

@kevinmoran1988

Copy link
Copy Markdown
Contributor

This follows up on #35. There I proposed a second program, bathscan, and you suggested not adding a tool but letting bathsearch pick the loop order itself. This PR goes one step further than a switch: the one-pass driver becomes bathsearch's only driver, so there is no mode to choose and no threshold to tune. That is possible because, with the changes it builds on, it takes the same time as the per-query driver for a single profile and less for two or more.

What it builds on

I'm posting this now so the driver can be reviewed, but it can't be merged yet: #37 and TravisWheelerLab/easel#7 are prerequisites. The driver relies on the ORF block being read-only to the pipeline, and on hits' traces being allocated at their own length. CI here will fail until the easel change is in.

Until those are in, this branch carries their commits, and the one from #36, underneath its own. The commits that belong to this PR are the last three: growing ov5 to the model in p7_pli_Frameshift(), two optional hooks on P7_PIPELINE for building frameshift profiles on demand, and the driver. Once the prerequisites are incorporated to your satisfaction, I'll rebase this on upstream main, which will leave just those three.

The timings below also have #36 and TravisWheelerLab/easel#6 on both sides.

What changes

bathsearch read and translated the whole target again for every query. It now reads and translates the target once, and searches every query against each block.

  • A ring of block slots. The reader fills a small ring of slots with blocks of target DNA, the same blocks as before. A slot's memory is set up the first time it is used, so a small target sets up one.
  • Translate once. A worker translates a slot once, both strands, into ORFs that are shared and read-only.
  • Unit of work. One unit is one query on one translated slot. A free worker takes the oldest slot that still has a query to run, runs p7_Pipeline_BATH() for that pair, and comes back for more. A slot is freed when every query has run on it.
  • One state per query, not per thread. A query's profiles, pipeline and hit list are a lane, used by one worker at a time.
  • Slow queries get more workers on their own. A query starts with one lane, so one worker at a time. When a worker finds nothing to run or translate, because every query still waiting on a slot is busy, it adds a lane to one of them and runs it. So a slow profile picks up extra workers as the search goes, with nothing to tune. In the 200-profile run below, 15 to 18 lanes were added this way.
  • DP memory belongs to the worker. It is lent to the query for each unit. After a unit the worker's full matrices go back to their starting size, so no worker keeps the largest one it has met.
  • Frameshift profiles on demand. With --fs, a lane's frameshift profiles are built through the pipeline hooks the first time a window needs them.
  • Serial mode. With --cpu 0, or in a build without threads, each block is read, translated and searched with every query in turn.

After the pass, each query's hits are thresholded and written as before, including the splice step, which still runs per query.

Results

Per-query driver → one-pass driver, built from the same code apart from the driver commit. 352 Mb genome, 16 threads unless noted, AVX2. Median of 3 runs for the first six rows, of 2 for the 50-profile rows, single runs for 200 profiles.

search time CPU time peak memory
1 profile, 115 positions 0.42 → 0.42 s 4.0 → 4.2 s 66 → 58 MB
1 profile, 647 positions 0.41 → 0.41 s 6.1 → 6.3 s 77 → 109 MB
1 profile, 1,578 positions 0.67 → 0.68 s 10.4 → 10.5 s 125 → 135 MB
1 profile, 115 positions, 1 thread 3.01 → 3.03 s 3.2 → 3.2 s 12 → 24 MB
2 profiles 0.81 → 0.50 s 10.0 → 7.8 s 85 → 74 MB
10 profiles 4.75 → 2.85 s 69 → 43 s 312 → 220 MB
50 profiles 19.0 → 10.2 s 288 → 158 s 417 → 256 MB
50 profiles, --fs 22.8 → 13.4 s 336 → 208 s 904 → 970 MB
200 profiles 87.4 → 49.0 s 1,343 → 771 s 1.41 → 0.99 GB
200 profiles, --fs 111.4 → 69.2 s 1,615 → 1,084 s 3.54 → 3.17 GB

523 small single-profile searches captured from a real run take 7.5–7.6 s at one thread where they took 7.7–8.1 s, and 7.2–7.3 s with --cpu 0 where they took 7.0–7.1 s.

So a single profile costs the same, and two or more take 38–46% less time. Memory is lower from two profiles up, except --fs at 50 profiles, which is about 7% higher. A single profile can use more: 30 MB more for the 647-position one, and 12 MB more at one thread.

What is different for a user

  • --cpu 0 now gives the threaded result. The old serial loop differed slightly from the threaded one. On my test set that was one E-value in the last digit in default mode, and 17 of 7,306 hits with --fs. The new driver gives the old threaded output at every thread count.
  • Several queries against a target that can't be rewound now work. The target is read once, so stdin and gzip input are fine with a multi-query file. i17-stdin.pl expected a failure there; it now checks that the piped target gives the same result as the file.
  • The per-query timing lines report the pass. The # CPU time and Elapsed lines printed for each query now give the time since the pass began, because queries are no longer searched one after another.
  • Memory grows with the number of profiles. Every profile is held for the whole pass, about 1.1 MB each, where the per-query driver's memory did not depend on the number of profiles. In Proposal: bathscan, searching many profiles against a genome in one pass #35 I said I would probably need a safety gate for very large profile sets. I don't think it needs one now: at 200 profiles the peak is below the per-query driver's, with and without --fs, and a search of 5,222 profiles against this genome took 20 minutes and peaked at 5.9 GB.

What I checked

  • The hit tables are identical between the two drivers in every run in the table, and in all 523 small searches.
  • The full output, alignments included, is identical to the old threaded driver's for 23 profiles against 14 Mb in default, --fs and --splice modes, and --cpu 0 matches 16 threads in all three.
  • An MSA query and a sequence query give the same output as before.
  • make check passes (easel 87/87, BATH 161/161), with i17-stdin.pl updated as above.

What I haven't done

  • Run it on macOS or ARM.
  • Built and run a configuration without thread support. That path compiles in a syntax check, and the same code runs under --cpu 0.
  • Compared the 5,222-profile search against the per-query driver. I only ran it with the new one.

This replaces a lot of bathsearch.c, so if you'd rather keep the per-query driver available behind an option, or want this split or arranged differently, I'm happy to do that.

🤖 Generated with Claude Code

@traviswheeler

Copy link
Copy Markdown
Member

The speed/memory gains are noteworthy. This is a substantial change, so I expect it will take a while for Genevieve to work through the evaluation. I won't rush that.

In the meantime, I think it's important to have a gate for very large profile sets. For example, it's not hard to imagine someone wanting to search with geNomad, which has >200,000 family profiles.

@kevinmoran1988

Copy link
Copy Markdown
Contributor Author

Thanks, and no rush at all on the evaluation.

You're right about very large profile sets. I checked how memory grows with the number of profiles, using the same 647-position profile repeated N times so that only the count changes (14 Mb target, 16 threads, peak memory):

profiles per-query driver this driver
50 87 MB 56 MB
200 88 MB 147 MB
800 92 MB 522 MB

The per-query driver stays flat and this one grows by about 0.6 MB per profile (about 1.6 MB with --fs), so something like geNomad would need well over 100 GB. What I wrote in the description, that it didn't need a gate, only holds for the set sizes I tested.

What I'd propose is a limit on the number of profiles searched in one pass. Profiles would be read and searched in batches of at most that many, with the target rewound between batches, and the output would stay in profile order as it is now. Two details:

  • Batch size. I haven't tested this yet, so I don't have a default to suggest. I'd measure time and memory across a range of sizes first and report back, and add an option to change it.
  • Targets that can't be rewound. With more profiles than one batch and a target on stdin or gzipped, I'd fail with a message, as bathsearch does today for any multi-profile search of such a target.

If that sounds reasonable I'll add it to this PR as a separate commit.

To be clear about the scope of this: the risk is only in the new driver. The other open PRs (#36, #37, TravisWheelerLab/easel#6 and TravisWheelerLab/easel#7) are speedups and memory savings in their own right and don't depend on it. The per-query column above already has them in, and it stays flat, so taking them carries no risk of running out of memory with large profile sets.

Separately, I've pushed one more commit here. It fixes a slowdown I missed: my timings were all at 16 threads, and with more threads than that a search with one short profile was slower than with the per-query driver (0.59 s against 0.48 s at 24 threads, 0.70 s against 0.50 s at 32). Every change of state woke all the waiting workers, and they queued for the driver's lock while the reader waited behind them to hand over each block. The commit gives the reader its own condition variable and wakes one worker for new work. That search now takes 0.28 s at 16 threads and 0.32 s at 32. Searches with many profiles are unchanged or 2–3% faster, and the output is identical.

I'm also chasing down a few more possible savings, and will post whatever holds up.

🤖 Generated with Claude Code

@kevinmoran1988

Copy link
Copy Markdown
Contributor Author

I've added one more small commit. With --fs, each worker was building its pair of 5-codon profiles once, sized for the longest query, and keeping it for the whole pass. It now builds the pair for the query it is working on. Peak memory with --fs falls from 974 to 866 MB for 50 profiles and from 3.18 to 2.66 GB for 200, with no change in time or output. This is one of the savings I mentioned looking into.

🤖 Generated with Claude Code

@kevinmoran1988

Copy link
Copy Markdown
Contributor Author

I've added the gate as a commit here.

bathsearch now reads and searches the queries a batch at a time: --qbatch <n> queries per pass over the target, 250 by default. A batch's profiles and hits are freed before the next batch is read, so memory no longer grows with the number of queries. A search with 250 queries or fewer behaves exactly as before.

On the two details I raised:

  • Batch size. I've measured it now. 2,000 profiles against the 352 Mb genome, 16 threads (these times have pipeline: reject ORFs on their SSV maximum #39 applied as well):

    batch passes time peak memory
    all 2,000 1 422 s 3.1 GB
    1,000 2 408 s 2.7 GB
    500 4 412 s 2.0 GB
    250 8 403 s 1.7 GB

    There is no time penalty: reading and translating the target again is small next to searching a few hundred profiles. It only shows with very small batches. With 200 profiles, batches of 25 were 3% slower and batches of 10 were 10% slower. Memory keeps falling as the batch shrinks but flattens out, which is why I settled on 250. With --fs, 500 profiles went from 4.1 GB in one pass to 3.6 GB in batches of 250.

  • Targets that can't be rewound. I said I would fail with a message, but there is no need. A target on stdin or gzipped simply gets every query in its one pass, as it did before this commit. Nothing that worked stops working; that search just isn't bounded.

So a set the size of geNomad should now need a few GB rather than well over 100.

The size of the target is not a concern in the same way. The target is read in blocks, and each block's ORFs go into the same few slots, which are emptied and reused. So that part is a few MB whatever the size of the genome: it was 3 MB in a heap profile of the 352 Mb search. What does grow with the target is the hits, which are held until a batch's pass ends, and a smaller --qbatch reduces that as well. I haven't run a target larger than 352 Mb.

What batching doesn't lower is the floor under those numbers, which is the DP matrices of domain definition, the frameshift ones above all. A long profile against a long stretch of DNA needs a matrix of profile length × DNA length, and the largest I saw was 1.1 GB. That is the same in every driver, old or new. There are savings left there, but they would require shenanigans, so I've left it alone.

The output is identical for every batch size I tried, in default, --fs and --splice modes.

🤖 Generated with Claude Code

kevinmoran1988 and others added 11 commits October 5, 2026 15:53
p7_pli_Frameshift() recorded which DNA window an ORF belongs to in
orfsq->idx, the only thing the pipeline wrote into the ORF block it is
given. Keep it in the per-call workspace instead (pli_tmp->orf_win[],
starting at -1). The ORF block is then read-only to the pipeline, and
an ORF needs no idx field.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Nothing calls it, and with the window index kept in the pipeline's
workspace it was the last code that wrote an ORF's idx.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The pipeline reads only the residues, length and coordinates of an ORF,
so it now takes easel's ESL_ORF_BLOCK in place of a block of ESL_SQ
objects, about 75 bytes per ORF where an ESL_SQ took about 1 kB.
bathsearch creates and empties its ORF blocks with the matching easel
calls, and the seed-window scoring in p7_hmmwindow.c and the pipeline
unit test pass their translated sequence as an ESL_ORF. This needs the
easel commit that adds ESL_ORF_BLOCK.

Output is unchanged. For one profile against a 352 Mb genome with 16
threads, peak memory falls from 133 MB to 66 MB and CPU time by 17%.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The clone took the source's allocation sizes, nalloc and ndomalloc, and
copied that much. The domain definition's working trace keeps the size
of the longest trace it has held, so every hit's copy was that large:
about 2,000 states allocated for an average of 122 used, 34 kB a hit.
The clone now allocates and copies what the trace holds.

Hits are unchanged. This is a prerequisite for the one-pass bathsearch
driver, which holds the hits of every query until the pass ends: with
it, a search of 5,222 profiles against a 352 Mb genome peaks at 5.9 GB
where it peaked at 28.0 GB.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Each step of a band folded all its vectors into one running maximum, so
every max waited for the one before it. Bands of three or more vectors
now keep separate maxima for even and odd vectors and join them when the
band ends, as the NEON kernel already does with six. Bands of one or two
vectors keep a single maximum; splitting there only adds register copies.

Scores are unchanged. For models of 100 or more positions and sequences
of 30 or more residues, SSV calls are 18-48% faster with AVX2 or SSE.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Frameshift domain definition runs om_fs5 through pli->ov5, which was only
ever sized when the pipeline was created. Grow it next to ov3, so a pipeline
created for a smaller model can run a larger one.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Two optional callbacks on P7_PIPELINE. fs_prepare() is called before
om_fs3 is first read in a p7_pli_Frameshift() call, and fs5_get()
supplies om_fs5 and gm_fs5 when a window reaches frameshift domain
definition. A caller can then build those profiles only when a search
gets that far. Both are NULL by default, which leaves the pipeline as it
was; ov5 is grown and om_fs5 configured after the 5-codon pair is known.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
bathsearch read and translated the whole target again for every query.
It now reads and translates it once, and searches every query against
each block.

- The reader fills a small ring of block slots. A worker translates a
  slot once, both strands, into ORFs that are shared and read-only.
- A unit of work is one query on one translated slot. A query's state
  (profiles, pipeline, hit list) is a lane used by one worker at a
  time, so nothing is copied per thread. A worker with nothing to run
  adds a lane to a query the others are waiting on.
- DP memory belongs to the worker and is lent to the query for each
  unit. After a unit the worker's full matrices go back to their
  starting size, so no worker keeps the largest one it has met.
- With --fs, a lane's frameshift profiles are built through the
  pipeline hooks, the first time a window needs them.
- With --cpu 0, or in a build without threads, each block is read,
  translated and searched with every query in turn.

Output is the same as the old threaded driver's at every thread count,
including --cpu 0, where the old serial loop differed slightly from the
threaded one. A target that can't be rewound (stdin, gzip) can now be
searched with several queries, and i17-stdin.pl checks that instead of
expecting a failure.

352 Mb genome, 16 threads, against the per-query driver: one profile
takes the same time; 10 profiles 4.8 s to 2.9 s; 200 profiles 83 s to
47 s, and with --fs 112 s to 67 s. Peak memory for 200 profiles is
1.1 GB against 1.4 GB, and 3.3 GB against 3.9 GB with --fs.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The one-pass driver woke every waiting worker at each change of state.
With more workers than work they all queued for the lock each time, and
the reader waited behind them to hand over a block. One 115-position
profile against a 352 Mb genome took 0.41 s with 16 threads, against
0.32 s with 12, and 0.70 s with 32.

The reader now has its own condition variable, and new work wakes one
worker; a worker that takes work wakes one more, which looks for work in
its turn. The same search takes 0.28 s with 16 threads and 0.32 s with
32. Searches with many profiles are unchanged or 2-3% faster.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
With --fs, each worker built its pair of 5-codon profiles once, for the
longest query, and kept it for the whole pass. It now builds the pair
for the query at hand when the query changes.

Peak memory with --fs against a 352 Mb genome falls from 974 to 866 MB
for 50 profiles and from 3.18 to 2.66 GB for 200. The time and the
output are unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The one-pass driver held every query's profiles and hits for the whole
pass, so its memory grew with the number of queries: about 0.5 MB each,
and about 2 MB with --fs. It now reads and searches the queries --qbatch
at a time (default 250), each batch in one pass over the target, and
frees a batch before it reads the next. A search with no more queries
than one batch is unchanged.

A target that can't be read again (stdin, gzip) gets every query in its
one pass, as before.

For 2,000 profiles against a 352 Mb genome, peak memory falls from 3.1
to 1.7 GB and the time does not change. The output is the same for every
batch size.

Co-Authored-By: Claude Opus 5.5 <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