Repository navigation
Faster, exact forward-backward; in-process seed anchoring; array-based progressive MSA - #9
Merged
Merged
Conversation
…d progressive MSA Forward-backward. The posterior probabilities (getPosteriorProbsWithBanding with diagonalCalculationPosteriorMatchProbs or ...PosteriorProbs) now go through a probability-space version of the same banded, chunked recursion. Each diagonal keeps a log scale that is moved only when its values drift out of [1e-100,1e100], the state machine is flattened into arrays once and cached on the StateMachine, and each state of a diagonal is one vectorised loop, with AVX2 gathers for the emissions. The band, the traceback points and which diagonals get posteriors in which chunk are exactly those of the log-space code, which remains for EM. It matches an exact reference to 1e-7; the old logAdd, a piecewise cubic that drops terms under e^-7.5, was off by up to ~0.01. The chunks are planned from the band before any dp, so buffers are sized once, and a chunk whose forward values would pass 16MB keeps only a pair of diagonals before each block and recomputes the block on the way back: an unanchored 10kb x 10kb pair drops from 3.9GB to 116MB. 5-11x faster per pairwise alignment. Seed anchoring (anchorMethod = PAIRWISE_ANCHOR_SEED): lastz's recipe -- 12of19 spaced seeds allowing a transition, ungapped HOXD70 x-drop extension, the heaviest colinear chain of the HSPs -- in process, recursing into the gaps with a threshold that falls with their area. MUM anchors accept "unique" matches of a dozen bases; at 0.3 substitutions per site per branch half of them were off the true alignment, and the banded dp followed them. In simulation at that divergence the MSA goes from 0.68/0.32 precision/recall with MUMs to 0.79/0.48. The default is unchanged (PAIRWISE_ANCHOR_LEGACY); cactus selects seed. Progressive MSA. getMultipleSequenceAlignmentProgressive and filterMultipleAlignedPairs work over arrays indexed by (sequence, position) instead of hashes and AVL trees, with the weights in one block, and no longer maintain the weight-ordered set only the greedy method reads. Same alignment, same order; 2.3x faster and half the memory on 6 x 100kb. column_hashFn no longer sums sequence and position, which made every bucket hold one of each sequence. Also: mum_destruct is a loop rather than a recursion as deep as the MUM chain; a zero-length HSP cannot enter the chain; threshold 0 reports every cell as the log path does. Tests: the fast path against the log-space one (random bands, machines, ragged ends, and unanchored pairs big enough to be recomputed in blocks) and seed anchor validity. test_filterToRemoveOverlap compared each pair with itself, so the set it expected was always empty; it now skips that. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Merged
4 tasks done
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Makes pecan practical again as cactus bar's base aligner: faster and exact posteriors, better anchors, and a cheaper progressive MSA. Used by ComparativeGenomicsToolkit/cactus#2026, which pins this branch's commit (12ce978).
logAddapproximation was off by up to ~0.01. It is 5–11× faster per pairwise alignment. Chunks whose forward matrix would exceed 16 MB recompute blocks on the way back, so an unanchored 10 kb × 10 kb pair uses 116 MB instead of 3.9 GB.PAIRWISE_ANCHOR_SEED): lastz's recipe, in process. It uses 12of19 spaced seeds that allow a transition, ungapped HOXD70 x-drop extension, and the heaviest colinear chain, recursing into the gaps. MUM anchors accept "unique" matches of about a dozen bases; at 0.3 substitutions per site per branch, half of them are off the true alignment. In simulation at that divergence, MSA precision/recall goes from 0.68/0.32 to 0.79/0.48. The library default is unchanged (PAIRWISE_ANCHOR_LEGACY); cactus selects seed in its config.column_hashFnno longer puts one column of each sequence in every bucket.mum_destructrecursed once per chain element (now a loop); a zero-length HSP could enter the chain; threshold 0 now reports every cell, as the log path does.test_filterToRemoveOverlapcompared each pair with itself, so it always expected an empty set.Results in cactus (bar with
baseAligner="pecan")Testing
cPecanLibTestspasses, including under ASan. It also builds with the SSE/AVX macros undefined (the path ARM takes).Please merge with a merge commit rather than squash/rebase, so the commit cactus pins stays reachable.
🤖 Generated with Claude Code