diff --git a/.gitignore b/.gitignore index 8c0589e..195e129 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,15 @@ bench/cache/ .coverage .coverage.* htmlcov/ + +# scratch worktrees and probe output +tmp/ + +# Local agent instructions and planning notes +/AGENTS.md +/CLAUDE.md +/ROADMAP.md +/SOURCES.md +/skills/ +/add_bindings.py +/add_native.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ffcbf9..c4bc7dd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,11 +3,57 @@ All notable changes to `seqtree`. Dates are release dates. From 1.0.0 the project follows semantic versioning: breaking changes need a **major** bump. -## [Unreleased] +## [1.1.0] — 2026-10-11 ### Added -- **1.0.2 candidate: bounded exhaustive gapblock reductions.** `topk_batch` keeps +- Native batched edit-count histograms, with an exact bit-vector fast path for bounded unit-cost searches. +- Gap-block length profiles and symmetric per-input positional weights, optional group penalties, + exhaustive balls, linked sum-distance counts, integer linear mass, and streaming group audits. +- TextIndex E-values and per-query mismatch cutoffs using distinct-record counts, with native + reduction of repeated offsets and rejection of truncated results. + +### Fixed + +- Stable small Poisson tails and tiny expectation probabilities. +- Shared native batches respect detected CPU allocations and reject oversubscribed + worker counts; serial budgets execute inline and all pools close before returning. +- Setup combines requested test and benchmark extras and propagates installation errors. +- Local agent instructions, planning notes and scratch scripts are excluded from distributions. + +### API details + +- **E-values for substring hits.** `text_evalues` and `text_threshold_for_evalue` + calibrate a `TextIndex` hit against a background text, the way `evalues` already + does for an `Index`: + + ```python + proteome = seqtree.TextIndex.build(records, "aa", k=5) + decoy = seqtree.TextIndex.build(shuffled, "aa", k=5) + for q, r in zip(peptides, seqtree.text_evalues(proteome, decoy, peptides, max_subs=2)): + if r["p_enrichment"] < 0.01: + ... + ``` + + The arithmetic is `evalue.evalue_result` unchanged; what is specific to text search + is two counting choices. **Counts are per distinct record, not per + offset**: one query typically hits a record at several overlapping offsets, those + events are strongly dependent, and counting them separately gives the + over-dispersed compound-Poisson count `appendix/evalue.tex` collapses away by + counting unique clonotypes rather than reads. `N` and `M` are `num_refs`. **The ball + is `max_subs`, never `TextResult.score`**: that score is a similarity, `0` at an + exact match and more negative for a worse one, the opposite sign to the + non-negative penalty `thetas_from_scores` inverts; `n_subs` is non-negative, + integer and monotone in the scope, so it nests the way the proof requires. + + Two limits, both refusals rather than partial answers. A search whose hits were + capped by `max_hits` raises, because a truncated count is a floor and an E-value + computed from it is simply wrong. And there is no `load_control` for texts: the + background must be composition- and length-matched to the target, that assumption + is what the whole calibration rests on, and nothing in the library can check it, so + it stays the caller's explicit choice. + +- **Bounded exhaustive gapblock reductions.** `topk_batch` keeps stable `(score, ref_id)` top-k rows without a dense score matrix; `count_batch` counts all ties at arbitrary per-query thresholds. Linked two-view variants use the maximum lane penalty and preserve both lane scores. Exact identity exclusion @@ -649,7 +695,7 @@ Numbers that constrain the API, all reproducible from `bench/` and the downstrea ### Docs -- `skills/seqtree/SKILL.md` — public API surface, invariants, and the gotchas that have bitten. +- Public API surface, invariants, and gotchas are documented in `docs/api.rst` and the API guides. - `docs/gapblock.rst` — a worked guide: why one gap block, how to choose its position, why a placement rule is a column frame, and why a fixed score cutoff is not a calibrated one. - README corrected: `seqtrie` is a full-width DP that ignores per-type caps, not a banded one, and diff --git a/CLAUDE.md b/CLAUDE.md deleted file mode 100644 index 4548db1..0000000 --- a/CLAUDE.md +++ /dev/null @@ -1,67 +0,0 @@ -# CLAUDE.md — working in this repo - -`seqtree` is a payload-agnostic fuzzy sequence-search core: C++20 + nanobind, **zero runtime -dependencies**, no domain data. What the library *exposes* is in `skills/seqtree/SKILL.md` — read -that first, it is the API reference and the gotcha list. This file is only how to work here. - -## Where things are - -| | | -|---|---| -| `src/`, `include/seqtree/` | the C++ core; `src/_bindings.cpp` is the whole nanobind surface | -| `python/seqtree/` | pure-Python layers (`evalue`, `gapblock`, `distance`, `pmhc`, `layout`, `seeds`) | -| `tests/cpp/`, `tests/python/` | doctest and pytest | -| `bench/` | benchmarks; `bench/tables/` holds the committed gate baselines | -| `docs/`, `docs/design/` | Sphinx site; `docs/design/*.md` are the design docs with their correction logs | -| `appendix/evalue.tex` | the E-value theory — treat it as the spec, not the code | -| `skills/seqtree/SKILL.md` | API reference for agents | -| `ROADMAP.md` | the contract with downstream `vdjmatch` / `mhcmatch` | -| `docs/roadmap.rst` | seqtree's own internal roadmap | - -## Build and test - -```bash -bash setup.sh --tests --bench -cmake -S . -B build -G Ninja -DSEQTREE_TESTS=ON -DCMAKE_BUILD_TYPE=Release && cmake --build build --parallel -ctest --test-dir build --output-on-failure -pytest tests/python -q -``` - -- **`editable.rebuild = false`** (pyproject). Editing C++ and running pytest tests the *old* - `_core.so` — reinstall (`uv pip install -e .`) after every C++ change or the result is a lie. -- Ninja occasionally misses a same-second edit; `touch src/foo.cpp` if it says "no work to do". -- `tests/cpp/test_codec.cpp` owns `DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN`. New test files must not - define a main, and go in the flat `add_executable` list in `CMakeLists.txt`. -- Gates: `RUN_PERF=1 pytest tests/python/test_perf_regression.py tests/python/test_oracle_tables.py -q -s`. - Regenerate with `bench/tables/make_tables.sh --perf` — **on an idle machine**. A baseline taken - while another benchmark saturates the cores reads as a 1.5x regression on untouched code. -- Heavy benchmarks (proteome downloads, 1M tiers) are gated behind `RUN_BENCHMARK=1`. -- Docs are built with `-W` in CI (`sphinx-build -b html -W --keep-going docs docs/_build/html`); - a warning is a failure. Locally: `env -C docs make html`. -- Windows CI is not optional — the `mmap` paths in `text_index.cpp` are the only platform-specific - code in the tree. - -## Conventions that aren't visible in the code - -- **Version is single-sourced** from `pyproject.toml`; `__version__` reads distribution metadata. - Never add a literal. -- **Gitflow**: feature branch → `dev` → `master`. Commit messages end with the `Co-Authored-By` - trailer. -- **Never publish to PyPI without an explicit release** (`ROADMAP.md` §4). Pushing a GitHub Release - fires `publish.yml` → OIDC upload. Tagging is fine; releasing is the owner's click. -- **seqtree stays generic.** No TCR germline data, no VDJdb, no clustering, no structures — those - live downstream in `vdjmatch` / `mhcmatch` / `tcren`. Don't import domain data into `bench/`. -- **Never fabricate a citation** — verify every DOI with a tool before it goes in a doc or comment. -- `docs/design/*.md` carry a status block listing corrections and *rejected* optimisations with - reasons. When a design decision is superseded, append a correction there rather than editing the - body — the rejections are load-bearing (they stop the same optimisation being re-proposed). - -## Open loops - -- **`v1.0.0` is tagged on `master` (`fc9dde4`) but not released.** CI is green; nothing is on PyPI. - The GitHub Release is the owner's to publish. -- **`TextIndex` at `L < 2k`** genuinely falls back to one block and a full radius-`max_subs` ball — - documented under *Limits*, not fixed. The fix would be spaced/gapped seeds (`docs/design/text_index.md` - §3.7), deferred by an earlier decision. -- Next design items are in `docs/roadmap.rst` (PSSM-graded d_TCR, native local alignment, Flashback - build); downstream wrapper design is in `ROADMAP.md` §2–3 with `(TBD: owner)` markers. diff --git a/ROADMAP.md b/ROADMAP.md deleted file mode 100644 index f7786e9..0000000 --- a/ROADMAP.md +++ /dev/null @@ -1,178 +0,0 @@ -# Downstream wrappers roadmap: `vdjmatch` and `mhcmatch` - -**Status:** living draft. Owner: @mikessh. Extend freely — sections marked _(TBD: owner)_ await -more detail. - -**Purpose.** `seqtree` is the shared substrate — a fast, payload-agnostic fuzzy-search core (C++ + -Python) with a control-calibrated E-value theory (see `appendix/evalue.tex`). It deliberately stops -at the reference-implementation level: the *applied* tools — `vdjmatch` (TCR antigen specificity) and -`mhcmatch` (peptide–MHC) — are separate packages built on top of it. This file is the contract -between them, written so an agent developing either package can pick it up cold: it states what -`seqtree` already provides, where each wrapper plugs in, and which design decisions are still open. - -**How to use this doc (for agents).** Before implementing a feature in `vdjmatch`/`mhcmatch`: -1. Check **§1 Substrate** for the `seqtree` primitive that already covers it — do not reimplement - search, E-values, anchor masking, or k-mer indexing. -2. Read the matching subsection in **§2** (vdjmatch) or **§3** (mhcmatch) for the intended design and - the relevant `appendix/evalue.tex` section. -3. Treat the math in the appendix as the spec; treat formulas here marked _(design)_ as direction, - not final — confirm with the owner before committing to one. - ---- - -## 1. Substrate: what `seqtree` provides - -| Capability | Module / symbol | Notes | -|---|---|---| -| Fuzzy fixed-length search | `seqtree.Index`, engines `seqtm` / `seqtrie` | seqtm: per-type edit caps + Hamming fast path. seqtrie: budget-only full-width DP over the trie (not banded), ignores per-type caps | -| Single-gap-block loop alignment | `seqtree.gapblock` (Python) | anchored CDR3/junction alignment: one contiguous indel + positional gap prior | -| Dense batch gap-block scoring | `seqtree.gapblock.score_matrix` → `ScoreMatrix` | every query × every reference in one GIL-released C++ call (~532 M pairs/s on 16 cores); `numpy.asarray` is zero-copy. The shape a prototype-distance embedding needs, where nothing can be pruned | -| Per-island profile (PWM) | `seqtree.gapblock.IslandProfile` | position weight matrix over a frame-aligned island; `pen(j,a) = round(lam·log(p_max_j/p_j(a)))` is `≥ 0` and `0` on the consensus, so it stays a ball and flows through `thetas_from_scores` unchanged. Beats min-over-members only at a strict cutoff | -| Substitution scoring | `SubstitutionMatrix` (Gram → squared-distance penalty `s_aa+s_bb−2·s_ab`) | payload-agnostic | -| Per-position scoring | `PositionalMatrix` (`penalty(pos,a,b)` = base × per-position weight; weight 0 = free/anchor) | the hook for PSSMs and anchor masking | -| k-mer seed index | `KmerIndex` (C++): unique-k-mer trie + CSR postings + per-peptide allele tag; `seed_and_gather` (GIL-released, parallel) | million-scale candidate generation | -| Control-calibrated E-values | `seqtree.evalues`, `seqtree.load_control` | `Ê = (N/M)·n_C`, Poisson tail, `exclude_exact` | -| Calibrated score cutoffs | `seqtree.threshold_for_evalue`, `thetas_from_scores` | inverts `Ê` to a **per-query** θ. A fixed θ is not calibrated: at `gapblock_score ≤ 60`, 31.7% of *random control* junctions land in a component of ≥5 — structure invented by the threshold, which per-query cutoffs remove entirely | -| Gap-placement priors | `central_prior`, `profile_prior`, `frame_prior`, `positions_prior`, `embed_in_frame` | `prior(i, d, m) ≥ 0` and `= 0` at `d = 0`. Only a **constant-`i`** rule (`frame_prior`) makes a frame transitive — a column index, hence a PWM, which `IslandProfile` now ships | -| Seed significance | `seqtree.seeds`: `core_kmers`, `SeedIndex` | precision, not recall: ~0.5% cross-island coverage | -| Background control | `seqtree.load_control`, `sanitize` | uniform reservoir sample over unique **productive** clonotypes. *Planned:* an out-of-frame control as an empirical `P_gen` sample for the `π̂_gen` fallback — needs the `*.ntvj` tables, since the `.aa` table's `_` has already collapsed a run of positions | -| Anchor / layout model | `seqtree.layout`: `AnchorSpec`, `DEFAULTS`, `mask_anchors`, `kmers`, `presentation_features`, `weight_profile` | parametrized anchors; class-II register trick | -| pMHC homology + reverse | `seqtree.pmhc`: `PMHCStore`, `search_homologs`, `assign_allele`, `find_mimics`, `build_kmer_index` | reference impl; mhcmatch productionizes | -| Presentation-aware E-values | `seqtree.pmhc_evalue.homolog_evalue` | per-allele null | - -**E-value theory** lives in `appendix/evalue.tex`: the empirical-control null (§Setup, §Null), the -Poisson/Chen–Stein bound (§Poisson), multiple testing (§E-value), the closest-hit Gumbel law -(§Gumbel), Karlin–Altschul as the i.i.d. special case (§KA), the pMHC presentation-aware extension -(§Epitopes), and elementary applications — UMI/barcode birthday-collision and -CDR3-nt error clustering (§Related applications). - -**Datasets.** `isalgo/pmhc_data` ships two tiers (see `appendix` Table "Scope of the two pmhc_data -tiers"): **full** (every IEDB-positive epitope–allele assay) and **shortlist** (epitope–allele pairs -with ≥2 supporting publications). `mhcmatch` should expose tier choice; the benchmark uses full. -TCR side: VDJdb (specificity-labelled) + `isalgo/airr_control` (matched control repertoire). - ---- - -## 2. `vdjmatch` — TCR antigen-specificity - -**Goal.** Given a query TCR (CDR3, optionally V/J, optionally paired α/β), score its similarity to -known antigen-specific TCRs and return a control-calibrated E-value per candidate epitope — -generalizing TCRNET-style neighbour counting (appendix §Intro) to a usable annotation tool. - -### 2.1 Train / test / validation splits -- **Split by epitope, not by sequence.** Held-out *epitopes* test generalization to unseen - specificities; held-out *sequences within an epitope* test within-specificity recall. Never let the - same clonotype leak across splits (dedup to unique clonotypes first — appendix Prop. on collapse). -- **Three-way:** train (fit substitution/positional weights, §2.2), validation (tune scope/budget and - E-value thresholds), test (report ROC/PR per epitope, as the MHC-guess benchmark does per allele). -- **Null/control** is independent of the split: `airr_control` supplies the background `P₀`; size it - per appendix §"How large must the control be?". -- _(TBD: owner)_ exact epitope-level CV scheme, minimum cluster size, negative sampling. - -### 2.2 Custom substitution matrices (epitope–paratope interaction-aware) -- **Hook:** `seqtree.PositionalMatrix` (per-position `penalty(pos,a,b)`) and/or a bespoke - `SubstitutionMatrix`. seqtree already wires `PositionalMatrix` into the seqtm Hamming path. -- **Intent:** replace BLOSUM with a CDR3-position-weighted matrix reflecting which residues actually - contact the peptide–MHC (paratope), learned from structural contacts or from specificity data - (residues whose substitution most changes specificity get the most weight; framework positions → - weight 0). Mirrors how `layout.weight_profile` builds anchor/TCR-facing weights for pMHC. -- _(TBD: owner)_ source of contact statistics (structures vs. learned), per-V-gene vs. global matrix. - -### 2.3 Single-chain E-value (have today) -- Per chain (β, or α), this is exactly `seqtree.evalues`: count neighbours of the query among the - epitope-specific set vs. the control, report `Ê` and `p_enrich`. No new theory needed. - -### 2.4 Paired α/β E-value (both chains known) -- _(design)_ Under chain independence in the null, the **joint** null ball-mass factorizes: - `π₀^{αβ} ≈ π₀^α · π₀^β`, so the paired intensity is `λ_αβ = N · π₀^α π₀^β` and the paired enrichment - is a Poisson tail on the count of references matching **both** chains within budget. Equivalent to - Fisher-combining the per-chain `p_enrich` when independence holds; the appendix `b₂` co-occupancy - term measures the dependence to correct for. Implement as a joint `seed_and_gather` keyed on a - paired peptide id. -- This is strictly more specific than either chain alone (joint null mass is tiny) → far smaller - E-values for true pairs. - -### 2.5 Paired estimate for a single chain (rarity of the unknown chain) -- _(design)_ When only one chain (say β) is observed, marginalize the paired E-value over the unknown - α: weight the single-chain evidence by the **rarity** of the partner that would complete a known - pair — a generation-probability / abundance prior on α (`Pgen`, appendix §thymic/selection). A β - matching a rare-α-paired reference is more informative than one matching a common-α pairing. Yields - a "paired-equivalent" E-value from one chain without observing the other. -- _(TBD: owner)_ the exact prior (OLGA/Murugan `Pgen` vs. empirical α-abundance), and calibration. - -### 2.6 API sketch _(non-binding)_ -```python -vm = vdjmatch.Annotator.from_vdjdb(chains=("beta",), matrix=paratope_matrix, control=airr_control) -hits = vm.annotate(cdr3b="CASS...", v="TRBV...", scope=2) # -> [(epitope, Ê, p_enrich), ...] -hits = vm.annotate_paired(cdr3a=..., cdr3b=...) # §2.4 -hits = vm.annotate(cdr3b=..., partner_prior="pgen") # §2.5 -``` - ---- - -## 3. `mhcmatch` — peptide–MHC - -**Goal.** Productionize the `seqtree.pmhc` reference layer: epitope homology / cross-reactivity, -molecular mimicry, allele guessing, and non-binder filtering — with tuned thresholds and the -additions below. - -### 3.1 Reuse from seqtree (have today) -- Homology / mimicry: `search_homologs`, `find_mimics` (anchor-masked TCR-facing k-mers; per-allele - presentation-aware E-values via `pmhc_evalue.homolog_evalue`). Positive control: the Dolton et al. - A\*02:01 cross-reactive trio. -- Allele guessing (reverse problem): `assign_allele` + the vote-fraction ranking / register trick - validated in `bench/bench_mhc_guess.py` (appendix §Epitopes). ROC-AUC 0.90–0.98. -- **Non-binder filter** (appendix §Epitopes, "Non-binders and class-II promiscuity"): _binds no MHC_ → - best-over-panel E-value high / no allele scores above background; _doesn't bind allele a_ → - per-allele `E_a > α`. mhcmatch exposes both thresholds. -- Tier choice (full vs shortlist, §1 Datasets). - -### 3.2 Motif logos -- For each allele (and each class-II core register), render a sequence logo of the anchor / pocket - residues from its presented set — the visual counterpart of `presentation_features`. Use a standard - information-content logo (bits per position). Class II: build the logo over the register-aligned - 9-mer cores (the register trick already picks the core), so anchors P1/P4/P6/P9 are columns. -- _(TBD: owner)_ logo library choice; whether to weight by `n_references` (shortlist confidence). - -### 3.3 MHC pseudosequence → clustering & promiscuity -- **Pseudosequence:** represent each allele by the polymorphic residues that line the peptide-binding - groove (the pseudosequence of peptide-contacting positions, as used for per-allele MHC binding - prediction by Glynn, Ghersi & Singh, *PNAS* 2025, [doi:10.1073/pnas.2405106122](https://doi.org/10.1073/pnas.2405106122)). - This turns "allele" from an opaque label into a short sequence the **same seqtree engine** can compare. -- **Cluster MHCs** by pseudosequence distance (seqtree fuzzy search on the pseudosequence "alphabet" - of groove residues) → groups of functionally similar alleles. This is the principled way to express - **cross-allele similarity**, which `seqtree.pmhc` deliberately does *not* model (appendix §Impl. - limitation: "distinct alleles are distinct nulls"). mhcmatch is where that limitation is lifted. -- **Promiscuity** (esp. class II — appendix §"Class-II promiscuity"): a peptide presented across a - pseudosequence-cluster of alleles is promiscuous; quantify as the spread of its positive alleles - over the MHC clustering, and use the cluster to pool nulls for related alleles when data are thin. -- _(TBD: owner)_ pseudosequence position set per locus (HLA-A/B/C, DR/DQ/DP, mouse H-2), distance - metric, cluster cut. - -### 3.4 API sketch _(non-binding)_ -```python -mm = mhcmatch.Store.from_pmhc(tier="shortlist") -mm.logo("HLA-A*02:01") # §3.2 -clusters = mm.cluster_alleles(class_="mhc1") # §3.3, via pseudosequence -mm.is_binder("SIINFEKL", allele="H-2Kb", alpha=0.05) # §3.1 non-binder filter -mm.promiscuity("PKYVKQNTLKLAT") # §3.3 -``` - ---- - -## 4. Shared conventions -- **seqtree is upstream and stays generic.** New general-purpose primitives (a learned matrix loader, - a logo helper, a pseudosequence comparator) may land in seqtree if reusable; tuned thresholds, - predictors, and domain glue stay in the wrappers. -- **Anchors** are parametrized in `seqtree.layout` (presets per class, overridable) — wrappers pass - `AnchorSpec`, they don't hardcode positions. -- **Citations:** never fabricate. Verify every DOI via a tool before adding it (PubMed/arXiv). -- **Versioning / gitflow:** feature branch → `dev` → `master`; end commit messages with the - `Co-Authored-By` trailer; don't publish to PyPI without an explicit release. - -## 5. Pointers -- E-value theory & all derivations: `appendix/evalue.tex` (compiled `appendix/evalue.pdf`). -- pMHC usage & limitations: `docs/pmhc.rst`. Benchmarks & figures: `docs/benchmarks.rst`, - `bench/bench_mhc_guess.py`. -- seqtree internal roadmap (PSSM-graded d_TCR, native local align, Flashback build, predictor proteome - scans): `docs/roadmap.rst`. diff --git a/SOURCES.md b/SOURCES.md deleted file mode 100644 index a3e291f..0000000 --- a/SOURCES.md +++ /dev/null @@ -1,118 +0,0 @@ -# Data sources - -Every dataset seqtree ships, downloads, or reads in a benchmark. Provenance is stated as -**experimental** (measured from a biological sample), **derived** (computed from experimental data), -or **synthetic** (generated by a model). Never conflate the three. - -## Bundled with the wheel - -### `python/seqtree/data/control_human_trb_aa.txt` - -| | | -|---|---| -| what | 250,000 human TRB junction amino-acid sequences, the null distribution `P₀` | -| provenance | **experimental**, then **derived** by uniform subsampling | -| origin | HuggingFace `antigenomics/vdjdb-cdr3-frequencies`, human TRB pooled clonotype table | -| format | one junction per line, productive only (20 standard amino acids) | -| regenerate | `python bench/gen_control.py --species human --chain trb --size 250000 --seed 0` | - -The upstream table is sorted by clonotype count, descending. The sample is a **seeded reservoir -sample over unique clonotypes**, then shuffled, so any prefix `bundled[:k]` is itself a uniform -sub-sample. Taking the first *k* rows of the upstream table instead yields the abundance head, which -is 25.8× more self-similar and carries 3.1× the ball mass — it violates `ass:indep` in -`appendix/evalue.tex` and corrupts every E-value. See `CHANGELOG.md` [0.3.0] Breaking. - -Out-of-frame clonotypes (VDJtools `_`) and in-frame stops (`*`) are dropped. `_` marks a *collapsed -run* of untranslatable positions, so the residue count is already lost and no amino-acid repair is -possible; and out-of-frame junctions escape thymic selection, making them a sample of `P_gen`, which -`lem:hierarchy` proves is not `P₀`. - -### Substitution matrices — `src/blosum62.inc`, `pam250.inc`, `pam100.inc`, `structural.inc` - -| | | -|---|---| -| provenance | **derived** (published log-odds / interaction-strength matrices) | -| origin | BLOSUM62 (NCBI); EMBOSS EPAM250, EPAM100; Miyazawa–Jernigan interaction strengths, transcribed into `MJ_CONTACT` in `bench/gen_matrices.py` from the `MJ` rows of `tcren-ms/src/tcren/data/MJ_Keskin_potentials.csv` (Miyazawa & Jernigan, J Mol Biol 1996, as tabulated by Keskin et al., Protein Sci 1998) | -| regenerate | `python bench/gen_matrices.py` | - -Stored as non-negative penalties via the Gram transform `pen[a][b] = s[a][a] + s[b][b] − 2·s[a][b]`, -clamped at zero, so `penalty(a, a) == 0` and the score defines a ball. - -**Known corruption in the MJ source.** `MJ_Keskin_potentials.csv` has no `A,N` row: the lower-triangle -slot where `N-A` belongs is mangled to `V,1` (mirrored `1,V`) for *both* the `MJ` and `Keskin` -potentials. `MJ_CONTACT` carries the correct `A-N = 0.15` from that row, so the constant will not -diff cleanly against the CSV at that pair. The `Keskin` rows are unused here. Fix not yet upstreamed -to `tcren-ms`. - -## Downloaded on demand - -### `load_control(name)` → `~/.cache/seqtree/control_*.sqtree` - -| | | -|---|---| -| what | larger control repertoires, e.g. `mouse_trb_aa` (694,241 productive clonotypes) | -| provenance | **experimental**, then **derived** by uniform subsampling | -| origin | HuggingFace `antigenomics/vdjdb-cdr3-frequencies` | -| fetch | `seqtree.load_control("mouse_trb_aa")` (uses `huggingface_hub`, an optional dependency) | - -The cached index is keyed by name only, **not** by content. Delete `~/.cache/seqtree/control_*.sqtree` -after upgrading across a control change. - -## Benchmark inputs — not shipped, passed by path - -These live in sibling repositories. seqtree does not depend on either; the benchmarks that read them -take `--prototypes` / `--vdjdb` arguments and parse the files with the standard library. - -### VDJdb slim dump - -| | | -|---|---| -| what | TCR↔epitope pairs with `v.end` / `j.start` templated-flank annotations | -| provenance | **experimental** (curated from published studies) | -| origin | `../mirpy/airr_benchmark/vdjdb/vdjdb-2025-12-29/vdjdb.slim.txt.gz` | -| upstream | https://github.com/antigenomics/vdjdb-db — release 2025-12-29 | -| covers | human TRA/TRB, mouse TRA/TRB, macaque TRB. Dotted column names (`v.segm`, `antigen.epitope`) | -| used by | `bench/gap_rules_by_locus.py` (criterion A) | - -### Per-locus prototype sets - -| | | -|---|---| -| what | 10,000 clonotypes per (species, locus): `v_call`, `j_call`, `junction_aa` | -| provenance | **experimental** where a real control repertoire exists, otherwise **synthetic** (OLGA) | -| origin | `../mirpy/mir/resources/prototypes/{species}_{locus}.tsv` | -| covers | human TRA TRB TRG TRD IGH IGK IGL; mouse TRA TRB | -| regenerate | `python ../mirpy/mir/resources/prototypes/generate_prototypes.py` (committed; never rebuilt at test time) | -| used by | `bench/gap_rules_by_locus.py` (criteria B and C) | - -Junctions, not IMGT CDR3: they include both anchors (`C…F` or `C…W`). Median lengths run from 12 -(IGK) to 20 (IGH), with IGH reaching 50 — which is why a fixed gap-position set calibrated on -length-14 TRB does not transfer. - -### IMGT gene library - -| | | -|---|---| -| what | germline V/D/J nucleotide sequences, 2,710 (species, allele) entries | -| provenance | **experimental** (IMGT reference directory) | -| origin | `../mirpy/mir/resources/gene_library/imgt_gene_library.txt` | -| caveat | **V genes translate in frame 0; J genes do not** (TRBJ2-7 is frame 1, TRBJ1-1 frame 2). The companion `region_annotations.txt` carries `jcdr3_aa`, which excludes the conserved F/W — it is CDR3 convention, not junction, and it disagrees with the translated library on some alleles (TRBV1\*01, IGHV1-2\*01). | - -Not currently read by any committed benchmark. Recorded because a germline-templated-flank prior was -built against it and **falsified**: it recovers the structurally correct block position 0.4% of the -time on TRA against 42.4% for a central prior. - -### UniProt reference proteomes (`isalgo/pmhc_data`) - -| | | -|---|---| -| what | whole-proteome FASTA: `mouse` (UP000000589, 54,857 records / 23,131,234 residues), `human` (UP000005640, 147,506 records / 69,578,135 residues), plus `ecoli_K12_UP000000625` and `fluA_PR8_UP000009255` as size steps | -| provenance | **experimental** (UniProt reference proteomes, mirrored) | -| origin | HuggingFace dataset `isalgo/pmhc_data`, path `proteome/.fasta.gz` | -| regenerate | `hf_hub_download("isalgo/pmhc_data", "proteome/human.fasta.gz", repo_type="dataset")` | -| used by | `bench/bench_text_index.py` (large tier only, `RUN_BENCHMARK=1`) | - -Carries residues the amino-acid codec does not name — **36 `U`** (selenocysteine) in human and -**33** in mouse — plus 8,417 human and 4,263 mouse `X`. `TextIndex` keeps the unnamed ones as -holes and reports the count, rather than refusing the build or silently treating them as -wildcards; `X`/`B`/`Z`/`*` are real symbols that match only themselves. diff --git a/bench/bench_text_index.py b/bench/bench_text_index.py index 28905ad..6ca9a23 100644 --- a/bench/bench_text_index.py +++ b/bench/bench_text_index.py @@ -22,7 +22,7 @@ caller's answer, it changes it, so recall below 1.000 anywhere means the work has failed. Texts, in size steps. The small ones are generated so the script runs with no downloads; the -real proteomes come from the HuggingFace dataset `isalgo/pmhc_data` (see SOURCES.md) and are +real proteomes come from the HuggingFace dataset `isalgo/pmhc_data` (see docs/data-sources.rst) and are used only when `huggingface_hub` is installed. Usage: @@ -71,7 +71,7 @@ def from_fasta_gz(path: str) -> list[str]: def proteome(name: str) -> list[str] | None: - """`isalgo/pmhc_data` proteomes, if huggingface_hub is available. See SOURCES.md.""" + """`isalgo/pmhc_data` proteomes, if huggingface_hub is available. See docs/data-sources.rst.""" try: from huggingface_hub import hf_hub_download except ImportError: diff --git a/bench/gap_rules_by_locus.py b/bench/gap_rules_by_locus.py index acd1df5..ba4a949 100644 --- a/bench/gap_rules_by_locus.py +++ b/bench/gap_rules_by_locus.py @@ -39,7 +39,7 @@ untemplated-span rule recovers the true block position 0.4% of the time for TRA, against 42.4% for a central prior). -Data (see SOURCES.md): mirpy's per-locus prototype files and a VDJdb slim dump. Neither ships with +Data (see docs/data-sources.rst): mirpy's per-locus prototype files and a VDJdb slim dump. Neither ships with seqtree; pass their paths. This script imports only seqtree and the standard library. Run: diff --git a/docs/api.rst b/docs/api.rst index a98d79b..04e2586 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -3,6 +3,20 @@ API Reference .. currentmodule:: seqtree +Batch CPU budgets +----------------- + +``threads=1`` runs inline. ``threads=0`` selects the detected available CPU budget: +Linux affinity, root cgroup v2 quota and SLURM allocation; macOS physical cores; +Windows process affinity. Legacy nonpositive values also select automatic mode. +Positive counts bound workers; a batch launches at most its query count. A request +above the detected allocation raises ``ValueError``; explicit counts are never clamped. +The caller performs one share of the work: a budget of B uses B compute threads +in total, creating only B-1 additional OS threads when there is enough work. +Pools are joined before the call returns; use one batch rather than an outer pool +of internally threaded searches. Cgroup layouts outside the detected root quota +and Windows CPU-rate job limits require an explicit thread budget. + Grouped by the job you are doing. For a starting point: :class:`Index` searches a *set of sequences*, :class:`TextIndex` searches a *long text*, and :mod:`seqtree.pairwise` aligns *two sequences*. @@ -18,6 +32,16 @@ payload. :members: :undoc-members: +For edit-count summaries, ``Index.edit_histogram_batch(queries, params, threads=0, exclude_exact=False)`` +returns one list of exact counts per query, indexed by total substitutions, +insertions and deletions. Set ``params.max_total_edits`` explicitly to a positive +value; bins run from zero through that cap. Reference duplicates count separately. +``exclude_exact=True`` retains bin zero with count zero. Nonpositive ``threads`` +uses the detected available CPU budget, matching the other batch APIs. This requires ``seqtm`` or ``auto`` +and all-hit mode; the search uses one reusable native hit buffer per worker rather +than materializing all hits in Python. Counts are not cumulative; sum the bins up +to a radius to obtain its ball count. + .. autoclass:: SearchParams :members: :undoc-members: @@ -117,6 +141,14 @@ cutoff that achieves a target false-positive rate. See :doc:`evalue`. .. autofunction:: seqtree.evalue.evalue_result +The same question for substring hits in a text, where counts are per distinct record +rather than per offset and the ball is a mismatch radius rather than a penalty. The +module docstring of :mod:`seqtree.text_evalue` explains both choices. + +.. autofunction:: text_evalues + +.. autofunction:: text_threshold_for_evalue + Pairwise alignment (Needleman-Wunsch / Smith-Waterman) ------------------------------------------------------ diff --git a/docs/data-sources.rst b/docs/data-sources.rst new file mode 100644 index 0000000..86316f3 --- /dev/null +++ b/docs/data-sources.rst @@ -0,0 +1,176 @@ +Data sources +============ + +Every dataset seqtree ships, downloads, or reads in a benchmark. Provenance is stated as +**experimental** (measured from a biological sample), **derived** (computed from experimental data), +or **synthetic** (generated by a model). Never conflate the three. + +Bundled with the wheel +---------------------- + +``python/seqtree/data/control_human_trb_aa.txt`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - 250,000 human TRB junction amino-acid sequences, the null distribution ``P₀`` + * - provenance + - **experimental**, then **derived** by uniform subsampling + * - origin + - HuggingFace ``antigenomics/vdjdb-cdr3-frequencies``, human TRB pooled clonotype table + * - format + - one junction per line, productive only (20 standard amino acids) + * - regenerate + - ``python bench/gen_control.py --species human --chain trb --size 250000 --seed 0`` + + +The upstream table is sorted by clonotype count, descending. The sample is a **seeded reservoir +sample over unique clonotypes**, then shuffled, so any prefix ``bundled[:k]`` is itself a uniform +sub-sample. Taking the first *k* rows of the upstream table instead yields the abundance head, which +is 25.8× more self-similar and carries 3.1× the ball mass — it violates ``ass:indep`` in +``appendix/evalue.tex`` and corrupts every E-value. See ``CHANGELOG.md`` [0.3.0] Breaking. + +Out-of-frame clonotypes (VDJtools ``_``) and in-frame stops (``*``) are dropped. ``_`` marks a *collapsed +run* of untranslatable positions, so the residue count is already lost and no amino-acid repair is +possible; and out-of-frame junctions escape thymic selection, making them a sample of ``P_gen``, which +``lem:hierarchy`` proves is not ``P₀``. + +Substitution matrices — ``src/blosum62.inc``, ``pam250.inc``, ``pam100.inc``, ``structural.inc`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - provenance + - **derived** (published log-odds / interaction-strength matrices) + * - origin + - BLOSUM62 (NCBI); EMBOSS EPAM250, EPAM100; Miyazawa–Jernigan interaction strengths, transcribed into ``MJ_CONTACT`` in ``bench/gen_matrices.py`` from the ``MJ`` rows of ``tcren-ms/src/tcren/data/MJ_Keskin_potentials.csv`` (Miyazawa & Jernigan, J Mol Biol 1996, as tabulated by Keskin et al., Protein Sci 1998) + * - regenerate + - ``python bench/gen_matrices.py`` + + +Stored as non-negative penalties via the Gram transform ``pen[a][b] = s[a][a] + s[b][b] − 2·s[a][b]``, +clamped at zero, so ``penalty(a, a) == 0`` and the score defines a ball. + +**Known corruption in the MJ source.** ``MJ_Keskin_potentials.csv`` has no ``A,N`` row: the lower-triangle +slot where ``N-A`` belongs is mangled to ``V,1`` (mirrored ``1,V``) for *both* the ``MJ`` and ``Keskin`` +potentials. ``MJ_CONTACT`` carries the correct ``A-N = 0.15`` from that row, so the constant will not +diff cleanly against the CSV at that pair. The ``Keskin`` rows are unused here. Fix not yet upstreamed +to ``tcren-ms``. + +Downloaded on demand +-------------------- + +``load_control(name)`` → ``~/.cache/seqtree/control_*.sqtree`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - larger control repertoires, e.g. ``mouse_trb_aa`` (694,241 productive clonotypes) + * - provenance + - **experimental**, then **derived** by uniform subsampling + * - origin + - HuggingFace ``antigenomics/vdjdb-cdr3-frequencies`` + * - fetch + - ``seqtree.load_control("mouse_trb_aa")`` (uses ``huggingface_hub``, an optional dependency) + + +The cached reference index is content-addressed by the bundled asset or download source, +requested size, alphabet and sampling seed. Superseded index files are pruned after rebuilding. + +Benchmark inputs — not shipped, passed by path +---------------------------------------------- + +These live in sibling repositories. seqtree does not depend on either; the benchmarks that read them +take ``--prototypes`` / ``--vdjdb`` arguments and parse the files with the standard library. + +VDJdb slim dump +~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - TCR↔epitope pairs with ``v.end`` / ``j.start`` templated-flank annotations + * - provenance + - **experimental** (curated from published studies) + * - origin + - ``../mirpy/airr_benchmark/vdjdb/vdjdb-2025-12-29/vdjdb.slim.txt.gz`` + * - upstream + - https://github.com/antigenomics/vdjdb-db — release 2025-12-29 + * - covers + - human TRA/TRB, mouse TRA/TRB, macaque TRB. Dotted column names (``v.segm``, ``antigen.epitope``) + * - used by + - ``bench/gap_rules_by_locus.py`` (criterion A) + + +Per-locus prototype sets +~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - 10,000 clonotypes per (species, locus): ``v_call``, ``j_call``, ``junction_aa`` + * - provenance + - **experimental** where a real control repertoire exists, otherwise **synthetic** (OLGA) + * - origin + - ``../mirpy/mir/resources/prototypes/{species}_{locus}.tsv`` + * - covers + - human TRA TRB TRG TRD IGH IGK IGL; mouse TRA TRB + * - regenerate + - ``python ../mirpy/mir/resources/prototypes/generate_prototypes.py`` (committed; never rebuilt at test time) + * - used by + - ``bench/gap_rules_by_locus.py`` (criteria B and C) + + +Junctions, not IMGT CDR3: they include both anchors (``C…F`` or ``C…W``). Median lengths run from 12 +(IGK) to 20 (IGH), with IGH reaching 50 — which is why a fixed gap-position set calibrated on +length-14 TRB does not transfer. + +IMGT gene library +~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - germline V/D/J nucleotide sequences, 2,710 (species, allele) entries + * - provenance + - **experimental** (IMGT reference directory) + * - origin + - ``../mirpy/mir/resources/gene_library/imgt_gene_library.txt`` + * - caveat + - **V genes translate in frame 0; J genes do not** (TRBJ2-7 is frame 1, TRBJ1-1 frame 2). The companion ``region_annotations.txt`` carries ``jcdr3_aa``, which excludes the conserved F/W — it is CDR3 convention, not junction, and it disagrees with the translated library on some alleles (TRBV1\*01, IGHV1-2\*01). + + +Not currently read by any committed benchmark. Recorded because a germline-templated-flank prior was +built against it and **falsified**: it recovers the structurally correct block position 0.4% of the +time on TRA against 42.4% for a central prior. + +UniProt reference proteomes (``isalgo/pmhc_data``) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. list-table:: + :widths: 20 80 + + * - what + - whole-proteome FASTA: ``mouse`` (UP000000589, 54,857 records / 23,131,234 residues), ``human`` (UP000005640, 147,506 records / 69,578,135 residues), plus ``ecoli_K12_UP000000625`` and ``fluA_PR8_UP000009255`` as size steps + * - provenance + - **experimental** (UniProt reference proteomes, mirrored) + * - origin + - HuggingFace dataset ``isalgo/pmhc_data``, path ``proteome/.fasta.gz`` + * - regenerate + - ``hf_hub_download("isalgo/pmhc_data", "proteome/human.fasta.gz", repo_type="dataset")`` + * - used by + - ``bench/bench_text_index.py`` (large tier only, ``RUN_BENCHMARK=1``) + + +Carries residues the amino-acid codec does not name — **36 ``U``** (selenocysteine) in human and +**33** in mouse — plus 8,417 human and 4,263 mouse ``X``. ``TextIndex`` keeps the unnamed ones as +holes and reports the count, rather than refusing the build or silently treating them as +wildcards; ``X``/``B``/``Z``/``*`` are real symbols that match only themselves. diff --git a/docs/design/text_index.md b/docs/design/text_index.md index 297cc91..dc7c271 100644 --- a/docs/design/text_index.md +++ b/docs/design/text_index.md @@ -1,7 +1,7 @@ # `TextIndex` — exact k-mismatch search over a concatenated reference text -**Status: implemented in `seqtree` 1.0.0, 2026-09-06.** Steps 1-6 of §8 shipped; §3.7 gapped -q-grams and §9's `mhcmatch` adapter did not. Six things below were corrected during +**Status: implemented in `seqtree` 1.0.0, 2026-09-06, and consumed by `mhcmatch` 1.15.0 the same +day.** Steps 1-8 of §8 shipped, including §9's `mhcmatch` adapter; §3.7 gapped q-grams did not. Six things below were corrected during implementation and the code, not this document, is authoritative on them: 1. **The amino-acid codec is 24 symbols** (`ARNDCQEGHILKMFPSTWYVBZX*`), not 20. So `A^k` at @@ -64,7 +64,7 @@ with room (0.127 ms/query against <= 0.5) and `L = 8-11, m = 2` lands at **1.24- `docs/text-index.rst` for the measured `k` guidance. Driven by a measured failure in `mhcmatch`, but the primitive is generic and stays generic -(`ROADMAP.md` §4, "seqtree is upstream and stays generic"). The reference implementation to beat is +(see the generic-core boundary in `docs/roadmap.rst`). The reference implementation to beat is **PEPMatch** (§3); this document sets out how, on seven axes, with the numbers for each. --- @@ -554,7 +554,7 @@ unfinished** — target **< 5 min / < 2 GB** on 8 threads. 8. **Only then** the `mhcmatch` adapter (§9), behind a keyword so the trie path stays reachable for one release and the two can be diffed on real data. -Gitflow per `ROADMAP.md` §4: `feature/text-index` → `dev` → `master`, `CHANGELOG.md` entry, 0.8.0. +Gitflow for this repository: `feature/text-index` → `dev` → `master`, `CHANGELOG.md` entry, 0.8.0. --- @@ -621,6 +621,27 @@ What `mhcmatch` then deletes: `Proteome._index`, `_index_key`, `_index_from_disk hand-off — an entire subsystem that exists only because the current index is expensive enough to be worth racing over. +**Shipped in mhcmatch 1.15.0, and it went as written.** `proteome.py` is 741 → 469 lines; everything +in the deletion list above is gone, along with `_Meta`, `MHCMATCH_INDEX_WAIT` and +`mhcmatch bootstrap --index`. Measured on the human proteome (147,506 records, 69,578,135 residues): +**0.7 s build / 0.6 GB peak RSS / 0.35 GB on disk**, against ~65 s / 12.6 GB / ~5.5 GB *per query +length* — and ~82 GB for the fifteen lengths MHC class II admits. Its `resolve_genes` stage went from +a recorded **2 h 18 min to 11.1 s** on 447,776 peptides over 45 distinct lengths, reproducing the +previous output with 0 differing lines of 593,842 rows. + +Two things the design did not anticipate, both worth knowing before writing another adapter: + +1. **`X` is a real symbol in the codec, so the caller must filter what the old index filtered at + build time.** `mhcmatch` dropped any window containing a non-standard residue when it enumerated + windows; `TextIndex` indexes the text as given, so an `X`-bearing window is a legitimate hit one + substitution from a query. Re-applying that filter on read is what made the port hit-for-hit + identical. +2. **`best_only` does not compose with a caller-side filter.** The search picks the nearest + non-empty shell *before* the caller sees it, so a shell the caller then empties yields nothing + rather than falling through to the next one. Real cost: 1 peptide of 445,466 silently lost its + five candidate parent genes. The adapter re-asks those few queries without `best_only`. Worth + considering whether a future `min_hits`/predicate hook belongs in the search itself. + --- ## 10. Out of scope diff --git a/docs/gapblock.rst b/docs/gapblock.rst index 1cb648c..88b0e8b 100644 --- a/docs/gapblock.rst +++ b/docs/gapblock.rst @@ -355,3 +355,80 @@ statistical model or calibration. ``Index.ref_seqs()`` retrieves all stored strings in reference-ID order through one native call, preserving duplicate rows. + +Positional weights in bounded batches +------------------------------------- + +``score_matrix``, ``count_batch`` and ``ball_batch`` accept +``position_weights_by_length``: a mapping from full sequence length to a vector +of nonnegative integer weights of that width. Supply every input length. +The longer sequence supplies the frame for each pair; prefix positions use their +original weights and the shifted suffix uses the corresponding longer-frame +weights. Equal-length scoring is weighted Hamming scoring. This pairwise frame +does not imply a consistent embedding of a whole collection. Gap charges and +position priors remain explicit and must use the same units as the weighted +substitution penalties. Invalid widths, fractional weights and possible int32 +overflow are rejected. Omitting the mapping preserves ordinary scoring. + +``ball_batch`` takes one integer cutoff per query and returns all accepted +reference rows, including ties, in reference order. Exact exclusion and duplicate +rows follow ``count_batch``. It avoids a dense score matrix but still evaluates +every reference; accepted output can itself be large. Positional weights are +currently unavailable in ``topk_batch``, ``paired_topk_batch`` and +``paired_count_batch``; ``paired_sum_count_batch`` accepts lane-specific weights. These operations supply +geometry, without a statistical calibration or a retrieval speed guarantee. + +Additive group distances +------------------------ + +``score_matrix``, ``count_batch`` and ``ball_batch`` also accept optional +``query_group_ids``, ``reference_group_ids`` and ``group_distances``. Supply all +three together: one nonnegative integer ID per sequence on each axis, and a +nonempty rectangular matrix of nonnegative integer penalties. The total pair +score is the ordinary gapblock score plus +``group_distances[query_group_ids[i]][reference_group_ids[j]]``. The table must +use the same units as substitution, gap and position-prior penalties. Group IDs +index rows on the query axis and columns on the reference axis. No symmetry or +zero diagonal is imposed; callers wanting a distance must provide those +properties themselves. + +Dimensions, IDs and possible int32 score overflow are checked before workers +start. The small table and encoded group IDs are shared read-only, with no dense +pair-offset allocation. Bounded reducers skip a pair when its group offset alone +exceeds the requested radius, since all other penalties are nonnegative. Exact +exclusion still punctures full encoded junction identity regardless of group ID +or offset, including when a matching junction has a positive group penalty. +Omitting all three options preserves the existing scorer. Group offsets are not +exposed by ``topk_batch``, ``paired_topk_batch`` or ``paired_count_batch``; +``paired_sum_count_batch`` accepts lane-specific groups. The paired maximum-score +contract is unchanged. These additions supply geometry, not statistical calibration. + +Per-input positional geometry +----------------------------- + +``score_matrix``, ``ball_batch`` and ``count_batch`` accept +``query_position_weights`` and ``reference_position_weights``: nonnegative integer +rows matching each full input sequence. Each aligned substitution uses the larger +of its query and reference weights. This rule is symmetric when the substitution +matrix is symmetric; gap costs and exact identity exclusion remain unchanged. +Supply both axes and do not combine them with ``position_weights_by_length``. + +``paired_sum_count_batch`` accepts separate ``_alpha`` and ``_beta`` versions of +these arguments. Each lane retains its own coordinates; distances sum against the +same linked reference row. These are geometry options, not confidence estimates. + +Streaming component and group audit +----------------------------------- + +``group_audit_batch`` accepts the same geometry as ``ball_batch``, with one +threshold per query. It returns three lists: the component count before group +offsets, the full-score count, and the minimum full score within the component +ball. An empty component ball has minimum ``-1``; duplicate references count +separately and encoded exact exclusion is preserved. + +If the component penalty is c and the nonnegative group offset is g, every +full-score hit c+g <= t is also a component hit c <= t. The reduction visits +references once with shared encoded inputs and per-worker alignment scratch; +output occupies O(Q) storage without sparse hit objects or a dense matrix. +The bounded minimum is not an unbounded nearest result. No statistical model +or cutoff recommendation is introduced. diff --git a/docs/index.rst b/docs/index.rst index b3b1082..c9029b1 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -119,3 +119,4 @@ map ``ref_id`` back to their own payloads (V gene, MHC, read counts) and filter api benchmarks roadmap + data-sources diff --git a/docs/roadmap.rst b/docs/roadmap.rst index 0b4d73f..bed9a6a 100644 --- a/docs/roadmap.rst +++ b/docs/roadmap.rst @@ -13,9 +13,11 @@ E-values and neoantigen mimic discovery (see :doc:`pmhc`). Planned next: 2. **Graded / native pMHC scoring** — ``PositionalMatrix`` already supports per-position PSSMs and graded TCR-hotspot weights; remaining work is a native local/best-window engine mode and predictor-driven proteome-to-presented-peptide pipelines for full mimic scans. -3. **Significance, continued.** Control-set + presentation-aware E-values ship; a *tf-idf redundancy - weighting* — down-weighting hits to common / clustered motifs without an external control set — is - still under consideration. +3. **Significance, continued.** Control-set, presentation-aware and substring (``TextIndex``) + E-values ship. Remaining: a bundled background *text* — there is no ``load_control`` analogue + for :func:`seqtree.text_evalues`, so the caller must supply a composition-matched decoy — and a + *tf-idf redundancy weighting*, down-weighting hits to common / clustered motifs without an + external control set, still under consideration. 4. **Memory** — succinct (LOUDS) trie / packed reference strings if reference counts reach the 10M+ range; memory-mapped frozen index for zero-copy load. 5. **Batch-of-batches & distribution** — streaming query batches and optional process-level fan-out. diff --git a/docs/text-index.rst b/docs/text-index.rst index 99d730b..1b4c30e 100644 --- a/docs/text-index.rst +++ b/docs/text-index.rst @@ -12,7 +12,8 @@ query is short and the reference is long and has no length of its own: length, exact answer. This is peptide-to-proteome mapping. It is not homology search: no scoring matrix in the -predicate, no E-value, no local alignment. +predicate, no local alignment. Significance is a separate layer: +:func:`seqtree.text_evalues` calibrates a hit against a background text. Why not ``Index`` ----------------- diff --git a/include/seqtree/kmer_index.hpp b/include/seqtree/kmer_index.hpp index e277f5f..623ab9d 100644 --- a/include/seqtree/kmer_index.hpp +++ b/include/seqtree/kmer_index.hpp @@ -39,7 +39,7 @@ class KmerIndex { // For each query (its k-mer list), return candidates with >= min_shared shared k-mers, // ranked by (shared_kmers desc, best_score asc). `params` drives the per-k-mer fuzzy // match (max_subs / pos_matrix / engine). allele_filter >= 0 keeps only peptides with - // that allele tag; < 0 means no filter. threads <= 0 => hardware_concurrency. + // that allele tag; < 0 means no filter. threads <= 0 => available CPU budget. std::vector> seed_and_gather( const std::vector>& query_kmers, const SearchParams& params, uint32_t min_shared = 1, int64_t allele_filter = -1, int threads = 0) const; diff --git a/include/seqtree/parallel.hpp b/include/seqtree/parallel.hpp index 6723007..6f9a92e 100644 --- a/include/seqtree/parallel.hpp +++ b/include/seqtree/parallel.hpp @@ -6,15 +6,68 @@ #include #include #include +#include +#include #include +#include #include +#include #include #include +#ifdef __linux__ +#include +#elif defined(__APPLE__) +#include +#elif defined(_WIN32) +#ifndef NOMINMAX +#define NOMINMAX +#endif +#ifndef WIN32_LEAN_AND_MEAN +#define WIN32_LEAN_AND_MEAN +#endif +#include +#endif + namespace seqtree { +// Resolve once per batch; affinity/quota limits can change between calls. +inline unsigned available_cpus() { + unsigned count = std::max(1u, std::thread::hardware_concurrency()); +#ifdef __linux__ + cpu_set_t affinity; + if (sched_getaffinity(0, sizeof(affinity), &affinity) == 0) + count = std::min(count, unsigned(std::max(1, CPU_COUNT(&affinity)))); + std::ifstream limits("/sys/fs/cgroup/cpu.max"); + uint64_t quota, period; + if (limits >> quota >> period && period) + count = unsigned(std::min(count, std::max(1, quota / period))); + for (const char* name : {"SLURM_CPUS_PER_TASK", "SLURM_CPUS_ON_NODE"}) { + if (const char* value = std::getenv(name)) { + char* end = nullptr; + const auto allocated = std::strtoul(value, &end, 10); + if (end != value && !*end && allocated) + count = unsigned(std::min(count, allocated)); + } + } +#elif defined(__APPLE__) + unsigned physical = count; + size_t size = sizeof(physical); + if (sysctlbyname("hw.physicalcpu", &physical, &size, nullptr, 0) == 0 && physical) + count = std::min(count, physical); +#elif defined(_WIN32) + DWORD_PTR process_mask, system_mask; + if (GetProcessAffinityMask(GetCurrentProcess(), &process_mask, &system_mask)) { + unsigned allowed = 0; + for (; process_mask; process_mask &= process_mask - 1) ++allowed; + if (allowed) count = std::min(count, allowed); + } +#endif + return count; +} + // Runs body(i, local) for every i in [0, n), across nt workers pulling contiguous chunks off -// one atomic cursor. `threads <= 0` means hardware concurrency. +// one atomic cursor. `threads <= 0` means the detected available CPU budget. // // `make_local()` is called once per worker and its result handed to every body() that worker // runs -- that is how a Searcher or a scratch buffer gets allocated once per thread rather @@ -31,11 +84,12 @@ namespace seqtree { template void parallel_for(size_t n, int threads, MakeLocal make_local, Body body, size_t chunk_min = 1, size_t chunk_max = 1024) { + const unsigned available = available_cpus(); + unsigned nt = threads > 0 ? unsigned(threads) : available; + if (nt > available || nt >= 1000) + throw std::invalid_argument("worker count exceeds the available CPU budget; request fewer threads"); if (n == 0) return; - - unsigned nt = threads > 0 ? unsigned(threads) - : std::max(1u, std::thread::hardware_concurrency()); - nt = unsigned(std::min(nt, n)); // n >= 1 here, so nt >= 1 + nt = unsigned(std::min(nt, n)); // Do not start idle workers for absent work. // A serial budget runs on the calling thread; spawning one worker only adds overhead. if (nt == 1) { @@ -70,7 +124,9 @@ void parallel_for(size_t n, int threads, MakeLocal make_local, Body body, // A failed thread launch must also join workers already started before unwinding. std::vector pool; pool.reserve(nt); - for (unsigned t = 0; t < nt; ++t) pool.emplace_back(worker); + // The caller is one compute thread: a budget B creates only B-1 OS workers. + for (unsigned t = 1; t < nt; ++t) pool.emplace_back(worker); + worker(); for (auto& th : pool) th.join(); if (err) std::rethrow_exception(err); } diff --git a/include/seqtree/seqtree.hpp b/include/seqtree/seqtree.hpp index 809f887..d8ac3c0 100644 --- a/include/seqtree/seqtree.hpp +++ b/include/seqtree/seqtree.hpp @@ -121,10 +121,15 @@ class Index { Alignment align(std::string_view query, uint32_t ref_id, const SearchParams&) const; // Parallel batch: pure C++ (no Python), one Searcher per worker thread. - // threads <= 0 => hardware_concurrency. Results aligned with input order. + // threads <= 0 => available CPU budget. Results aligned with input order. std::vector> search_batch(const std::vector& queries, const SearchParams&, int threads = 0) const; + // Exact edit-count bins 0..max_total_edits; one reusable hit buffer per worker. + std::vector> edit_histogram_batch( + const std::vector& queries, const SearchParams&, + int threads = 0, bool exclude_exact = false) const; + // Per-query count of seqtm collisions: how often a reference was re-reached via a // different edit path during branch-and-bound (0 for the seqtrie engine / Hamming). std::vector collisions_batch(const std::vector& queries, @@ -195,7 +200,7 @@ Alignment align_pair(std::string_view query, std::string_view ref, const Substit Alphabet, AlignMode, int32_t gap_open, int32_t gap_extend); // Dense N*K similarity matrix, row-major (row i is queries[i] vs every ref). Parallel; -// threads <= 0 => hardware_concurrency. +// threads <= 0 => available CPU budget. std::vector align_score_matrix(const std::vector& queries, const std::vector& refs, const SubstitutionMatrix&, Alphabet, AlignMode, @@ -215,16 +220,55 @@ std::vector align_dist_matrix(const std::vector& queries, // The block of length d = |m-n| is placed at whichever of the min(m,n)+1 positions minimises // substitution cost plus `prior`. `prior` is the gap prior flattened to [m][d][i] with stride // (prior_width+1); pass it empty for no prior. It is read only when d > 0, which is what keeps -// the diagonal zero. `matrix` may be null for unit cost. threads <= 0 => hardware_concurrency. +// the diagonal zero. `matrix` may be null for unit cost. threads <= 0 => available CPU budget. +// Optional integer weights are indexed [longer length][longer position], symmetric in q/r. +// All input frames must be covered; empty tables retain the existing unweighted scorer. +// Alternatively, per-input rows index each sequence's full coordinates. Aligned +// substitution costs use max(query weight, reference weight); gaps are unchanged. +// Both axes are required; these rows cannot combine with length profiles. +// Optional group IDs index a nonnegative rectangular penalty table added per pair. +// Exact exclusion remains encoded junction identity, independent of these groups. std::vector gapblock_matrix(const std::vector& queries, const std::vector& refs, Alphabet, const SubstitutionMatrix* matrix, int32_t gap_open, int32_t gap_extend, const std::vector& prior, - uint32_t prior_width, int threads = 0); + uint32_t prior_width, int threads = 0, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); + +// Exhaustive ball rows in reference order; hit edit fields are unset (zero). +std::vector> gapblock_ball_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); // Exhaustive bounded reductions using the same gapblock cell. Top-k ties use // (penalty, ref_id); paired penalty is max(lane scores). Exact exclusion uses // encoded sequence identity (both lanes for pairs), never zero penalty alone. +// Count the component ball before/after group offsets; minimum full score within it. +// One threshold per query; nearest=-1 denotes an empty component ball. O(Q) output. +std::tuple, std::vector, std::vector> +gapblock_group_audit_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); + std::vector> gapblock_topk_batch(const std::vector& queries, const std::vector& refs, uint32_t k, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector& prior, @@ -237,13 +281,45 @@ std::vector>> gapblock_paired std::vector> gapblock_count_batch(const std::vector& queries, const std::vector& refs, const std::vector>& thresholds, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, - const std::vector& prior, uint32_t width, int threads, bool exclude_exact); + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); +std::pair>, std::vector>> gapblock_count_mass_batch(const std::vector& queries, + const std::vector& refs, const std::vector>& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}); std::vector> gapblock_paired_count_batch( const std::vector& qa, const std::vector& qb, const std::vector& ra, const std::vector& rb, const std::vector>& thresholds, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact); +std::pair>, std::vector>> gapblock_paired_sum_count_batch( + const std::vector& qa, const std::vector& qb, + const std::vector& ra, const std::vector& rb, + const std::vector>& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights_alpha, + const std::vector>& weights_beta, + const std::vector& query_groups_alpha, const std::vector& query_groups_beta, + const std::vector& reference_groups_alpha, const std::vector& reference_groups_beta, + const std::vector>& group_distances_alpha, + const std::vector>& group_distances_beta, + const std::vector>& query_weights_alpha = {}, + const std::vector>& query_weights_beta = {}, + const std::vector>& reference_weights_alpha = {}, + const std::vector>& reference_weights_beta = {}); // --------------------------------------------------------------------------------------- // Plain string edit distances on raw characters, unit costs. Deliberately alphabet-agnostic: @@ -259,7 +335,7 @@ int hamming(std::string_view a, std::string_view b); int levenshtein(std::string_view a, std::string_view b); // Dense N*K distance matrices, row-major (row i is a[i] vs every b[k]). Parallel; -// threads <= 0 => hardware_concurrency. hamming_matrix throws on any length-mismatched pair. +// threads <= 0 => available CPU budget. hamming_matrix throws on any length-mismatched pair. std::vector hamming_matrix(const std::vector& a, const std::vector& b, int threads = 0); std::vector levenshtein_matrix(const std::vector& a, diff --git a/include/seqtree/text_index.hpp b/include/seqtree/text_index.hpp index 9072b84..4a505f9 100644 --- a/include/seqtree/text_index.hpp +++ b/include/seqtree/text_index.hpp @@ -70,6 +70,8 @@ struct TextResult { size_t num_queries() const { return query_begin.empty() ? 0 : query_begin.size() - 1; } size_t num_hits() const { return ref_id.size(); } + // Best substitution distance per distinct record, sorted per query. Reject capped hits. + std::vector> record_distances() const; }; // Exact k-mismatch (Hamming) search over a CONCATENATED reference text. @@ -108,7 +110,7 @@ class TextIndex { // Every query must be at least k long: seeds are only indexed where a whole k-mer fits // inside a record, so a shorter query could match within k-L of a record end and be - // missed. Raises rather than answering incompletely. threads <= 0 => hardware_concurrency. + // missed. Raises rather than answering incompletely. threads <= 0 => available CPU budget. TextResult search_batch(const std::vector& queries, const TextQueryOpts&, int threads = 0) const; diff --git a/pyproject.toml b/pyproject.toml index feb221b..ed02978 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "seqtree" -version = "1.0.2" +version = "1.1.0" description = "Fast fuzzy search over biological sequences (C++ core, Python bindings)" readme = "README.md" requires-python = ">=3.10" @@ -40,9 +40,10 @@ wheel.packages = ["python/seqtree"] cmake.build-type = "Release" # Keep the sdist lean: ship only what builds the package. Anchor with a leading "/". sdist.exclude = [ - "/bench", "/tests", "/docs", "/.github", "/setup.sh", "/build", + "/bench", "/tests", "/docs", "/.github", "/setup.sh", "/build", "/appendix", "/examples", # Local agent config: machine-specific absolute paths and unrelated repo names. - "/.claude", "/CLAUDE.md", + "/.claude", "/CLAUDE.md", "/AGENTS.md", "/ROADMAP.md", "/SOURCES.md", "/skills", + "/tmp", "/add_bindings.py", "/add_native.py", "/appendix/immunogenicity_hooks.md", ] editable.rebuild = false diff --git a/python/seqtree/__init__.py b/python/seqtree/__init__.py index e5e38cb..6067d49 100644 --- a/python/seqtree/__init__.py +++ b/python/seqtree/__init__.py @@ -34,6 +34,7 @@ ) from .control import load_control from .evalue import evalues, thetas_from_scores, threshold_for_evalue +from .text_evalue import text_evalues, text_threshold_for_evalue from . import distance, gapblock, layout, pairwise, pmhc, seeds from .gapblock import ( GapBlockIndex, IslandProfile, ScoreMatrix, central_prior, embed_in_frame, frame_prior, @@ -87,4 +88,6 @@ "evalues", "thetas_from_scores", "threshold_for_evalue", + "text_evalues", + "text_threshold_for_evalue", ] diff --git a/python/seqtree/evalue.py b/python/seqtree/evalue.py index 3a7dfb6..81822b9 100644 --- a/python/seqtree/evalue.py +++ b/python/seqtree/evalue.py @@ -22,17 +22,31 @@ def _poisson_sf(k, lam): - """P(Poisson(lam) >= k) via the complementary CDF (stable for small k).""" + """P(Poisson(lam) >= k), summing the smaller tail from a log-PMF seed.""" if k <= 0: return 1.0 if lam <= 0.0: return 0.0 - term = math.exp(-lam) - cdf = term - for i in range(1, k): - term *= lam / i - cdf += term - return min(1.0, max(0.0, 1.0 - cdf)) + if lam == math.inf: + return 1.0 + # Sum the upper tail directly when small: 1-CDF loses significant digits. + upper = k > lam + i = k if upper else k - 1 + log_pmf = i * math.log(lam) - lam - math.lgamma(i + 1) + # Normalize before summation so a subnormal seed cannot erase a representable tail. + term = total = 1.0 + while term > total * 1e-15: + if upper: + i += 1 + term *= lam / i + else: + if i == 0: + break + term *= i / lam + i -= 1 + total += term + probability = math.exp(log_pmf + math.log(total)) + return min(1.0, max(0.0, probability if upper else 1.0 - probability)) def evalue_result(n_target, n_control, n_ref, m_control): @@ -61,8 +75,7 @@ def evalue_result(n_target, n_control, n_ref, m_control): "n_target": n_target, "n_control": n_control, "E": E, - # exp(-E) underflows to 0.0 well before E = 700, but guard the call anyway. - "p_any": 1.0 - math.exp(-E) if E < 700 else 1.0, + "p_any": -math.expm1(-E), "p_enrichment": _poisson_sf(n_target, E), "rule_of_three": rule3, } diff --git a/python/seqtree/gapblock.py b/python/seqtree/gapblock.py index 913ce64..ab2e7bf 100644 --- a/python/seqtree/gapblock.py +++ b/python/seqtree/gapblock.py @@ -44,11 +44,15 @@ from ._core import gapblock_matrix as _gapblock_matrix from ._core import gapblock_topk_batch as _gapblock_topk_batch from ._core import gapblock_count_batch as _gapblock_count_batch +from ._core import gapblock_count_mass_batch as _gapblock_count_mass_batch +from ._core import gapblock_ball_batch as _gapblock_ball_batch +from ._core import gapblock_group_audit_batch as _gapblock_group_audit_batch from ._core import gapblock_paired_topk_batch as _gapblock_paired_topk_batch from ._core import gapblock_paired_count_batch as _gapblock_paired_count_batch +from ._core import gapblock_paired_sum_count_batch as _gapblock_paired_sum_count_batch __all__ = [ - "gapblock_score", "score_matrix", "topk_batch", "count_batch", "paired_topk_batch", "paired_count_batch", "deletion_variants", "central_prior", "profile_prior", + "gapblock_score", "score_matrix", "ball_batch", "group_audit_batch", "topk_batch", "count_batch", "paired_topk_batch", "paired_count_batch", "paired_sum_count_batch", "deletion_variants", "central_prior", "profile_prior", "frame_prior", "positions_prior", "embed_in_frame", "gap_cost", "GapBlockIndex", "ScoreMatrix", "IslandProfile", ] @@ -522,6 +526,9 @@ def score_matrix( gap_prior: GapPrior | None = None, alphabet: str = "aa", threads: int = 0, + position_weights_by_length=None, + query_group_ids=None, reference_group_ids=None, group_distances=None, + query_position_weights=None, reference_position_weights=None, ) -> ScoreMatrix: """Gap-block penalty of every query against every reference, in C++ with the GIL released. @@ -533,6 +540,18 @@ def score_matrix( Args: queries: Query sequences (the matrix rows). refs: Reference sequences (the matrix columns). + position_weights_by_length: Optional mapping of full sequence lengths to nonnegative + integer weights. Weights index the longer sequence; both axes must be covered. + None or an empty mapping preserves unweighted defaults. + query_position_weights: Optional nonnegative integer rows, one per query and residue. + Supply with reference_position_weights; incompatible with length profiles. + Substitutions use max(query_weight, reference_weight) at aligned coordinates. + Gap costs remain unchanged. + reference_position_weights: Per-reference rows with the same contract. + query_group_ids: Optional integer IDs, one per query; supply all group options together. + reference_group_ids: Optional integer IDs, one per reference. + group_distances: Nonempty rectangular nonnegative integer penalty table. + Added to each pair score in the same units; no symmetry/zero diagonal imposed. matrix: Substitution penalties; ``None`` means unit cost. gap_open: Block-opening cost. Defaults to ``2 * matrix.scale()``. See the module note: leaving this at 1 with a real matrix makes gaps ~14x cheaper than substitutions. @@ -565,7 +584,10 @@ def score_matrix( raise ValueError("gap_open and gap_extend must be >= 0") width = max((len(s) for s in q + r), default=0) cube = _prior_cube(gap_prior, width) if gap_prior is not None else [] - return _gapblock_matrix(q, r, alphabet, matrix, gap_open, gap_extend, cube, width, threads) + return _gapblock_matrix(q, r, alphabet, matrix, gap_open, gap_extend, cube, width, threads, + _position_weights(position_weights_by_length, q, r), + *_group_options(query_group_ids, reference_group_ids, group_distances, q, r), + *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r)) def _batch_integer(value, name, lower=0, upper=(1 << 31) - 1): @@ -590,6 +612,62 @@ def _batch_options(axes, matrix, gap_open, gap_extend, gap_prior, threads): return go, ge, cube, width, threads +def _group_options(query_ids, reference_ids, distances, queries, refs): + """Validate one small group table; native preparation rechecks its safety.""" + if query_ids is None and reference_ids is None and distances is None: + return [], [], [] + if query_ids is None or reference_ids is None or distances is None: + raise ValueError("query_group_ids, reference_group_ids and group_distances are required together") + qg = [_batch_integer(v, "query group ID") for v in query_ids] + rg = [_batch_integer(v, "reference group ID") for v in reference_ids] + if len(qg) != len(queries) or len(rg) != len(refs): + raise ValueError("group IDs must match their query/reference axes") + table = [[_batch_integer(v, "group distance") for v in row] for row in distances] + if not table or not table[0] or any(len(row) != len(table[0]) for row in table): + raise ValueError("group distances must be a nonempty rectangular matrix") + if any(v >= len(table) for v in qg) or any(v >= len(table[0]) for v in rg): + raise ValueError("group ID outside group distances") + return qg, rg, table + + +def _position_weights(weights, *axes): + """Validated integer tables in the longer-junction frame; never round/truncate.""" + if weights is None: + return [] + if not hasattr(weights, "items"): + raise ValueError("position_weights_by_length must map lengths to integer weights") + tables = {} + for length, row in weights.items(): + length = _batch_integer(length, "position frame", upper=65535) + row = list(row) + if len(row) != length: + raise ValueError("positional weight length must match its frame") + tables[length] = [_batch_integer(w, "position weight") for w in row] + if not tables: + return [] + for axis in axes: + for sequence in axis: + if len(sequence) and len(sequence) not in tables: + raise ValueError("missing positional weights for an input frame") + return [tables.get(length, []) for length in range(max(tables) + 1)] + + +def _sequence_weights(query_weights, reference_weights, by_length, queries, refs): + """Integer full-coordinate weights; aligned substitutions use the larger weight.""" + if query_weights is None and reference_weights is None: + return [], [] + if query_weights is None or reference_weights is None: + raise ValueError("both query and reference positional weights are required") + if by_length: + raise ValueError("length and per-input weights cannot be combined") + def validate(rows, axis): + rows = [list(row) for row in rows] + if len(rows) != len(axis) or any(len(row) != len(seq) for row, seq in zip(rows, axis)): + raise ValueError("per-input weights must match their sequence axis and lengths") + return [[_batch_integer(value, "position weight") for value in row] for row in rows] + return validate(query_weights, queries), validate(reference_weights, refs) + + def _threshold_rows(thresholds, n_queries): rows = [list(row) for row in thresholds] if len(rows) != n_queries: @@ -622,19 +700,83 @@ def topk_batch(queries, refs, k=10, matrix=None, gap_open=None, gap_extend=1, def count_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, - gap_prior=None, alphabet="aa", threads=0, exclude_exact=False): + gap_prior=None, alphabet="aa", threads=0, exclude_exact=False, position_weights_by_length=None, + query_group_ids=None, reference_group_ids=None, group_distances=None, linear_mass=False, + query_position_weights=None, reference_position_weights=None): """Count all reference scores <= each query's requested integer thresholds. + Optional positional weights and additive group distances match :func:`score_matrix`. + Exact exclusion remains full encoded identity, irrespective of group distance. ``thresholds`` has one row per query; rows may be empty, unsorted or contain duplicates. Return integer counts in their original order, including ALL score ties. Negative cutoffs count zero. Same scorer/exclusion as :func:`topk_batch`; one native batch with O(Q*S) output, no dense score matrix. Worker scratch is O(max_sequence_length + maximum_threshold_row_length). + With ``linear_mass=True``, return ``(counts, masses)`` from the same traversal; + each mass is ``sum(max(0, threshold - score))`` over nonexcluded references. + Boundary ties still count but contribute zero mass. Scores are nonnegative; + integer mass is uint64, with overflow bounds checked before workers. Negative + and zero thresholds have zero mass. No normalization or statistical model is + applied; default counts and paired APIs remain unchanged. """ + if not isinstance(linear_mass, bool): + raise TypeError("linear_mass must be a boolean") q, r = list(queries), list(refs) thresholds = _threshold_rows(thresholds, len(q)) go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads) - return _gapblock_count_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact) + reduction = _gapblock_count_mass_batch if linear_mass else _gapblock_count_batch + return reduction(q, r, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact, + _position_weights(position_weights_by_length, q, r), + *_group_options(query_group_ids, reference_group_ids, group_distances, q, r), + *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r)) + + +def group_audit_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, + gap_prior=None, alphabet="aa", threads=0, exclude_exact=False, + position_weights_by_length=None, + query_group_ids=None, reference_group_ids=None, group_distances=None, + query_position_weights=None, reference_position_weights=None): + """Stream component-ball count, full-score count and bounded nearest in one batch. + + Geometry/offset options match ball_batch. Component excludes the nonnegative + group offset; full score includes it. For each query threshold, return three + lists: component count, full count and minimum full score among component hits. + Empty component balls have nearest=-1. Duplicate references count separately; + full encoded exact puncturing is identical. O(Q) output; no sparse hit objects. + """ + q, r = list(queries), list(refs) + thresholds = [_batch_integer(v, "threshold", -(1 << 31)) for v in thresholds] + if len(thresholds) != len(q): + raise ValueError("group audit thresholds must have one value per query") + go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads) + return _gapblock_group_audit_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width, + threads, exclude_exact, _position_weights(position_weights_by_length, q, r), + *_group_options(query_group_ids, reference_group_ids, group_distances, q, r), + *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r)) + + +def ball_batch(queries, refs, thresholds, matrix=None, gap_open=None, gap_extend=1, + gap_prior=None, alphabet="aa", threads=0, exclude_exact=False, + position_weights_by_length=None, + query_group_ids=None, reference_group_ids=None, group_distances=None, + query_position_weights=None, reference_position_weights=None): + """Exhaustive bounded ball rows, reference order, unset edit fields; no Q*N matrix. + + One integer cutoff per query. Negative cutoffs yield empty rows. Optional integer + positional weights index the longer full sequence, symmetrically, at every length. + Optional additive group distances match :func:`score_matrix`; exact exclusion + remains full encoded identity, irrespective of group distance. + This reduces output memory; it is not a trie-based candidate filter. + """ + q, r = list(queries), list(refs) + thresholds = [_batch_integer(v, "threshold", -(1 << 31)) for v in thresholds] + if len(thresholds) != len(q): + raise ValueError("ball thresholds must have one value per query") + go, ge, cube, width, threads = _batch_options((q, r), matrix, gap_open, gap_extend, gap_prior, threads) + return _gapblock_ball_batch(q, r, thresholds, alphabet, matrix, go, ge, cube, width, + threads, exclude_exact, _position_weights(position_weights_by_length, q, r), + *_group_options(query_group_ids, reference_group_ids, group_distances, q, r), + *_sequence_weights(query_position_weights, reference_position_weights, position_weights_by_length, q, r)) def paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10, matrix=None, @@ -673,6 +815,43 @@ def paired_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, matrix= return _gapblock_paired_count_batch(qa, qb, ra, rb, thresholds, alphabet, matrix, go, ge, cube, width, threads, exclude_exact) +def paired_sum_count_batch(queries_a, queries_b, refs_a, refs_b, thresholds, matrix=None, + gap_open=None, gap_extend=1, gap_prior=None, alphabet="aa", + threads=0, exclude_exact=False, + position_weights_by_length_alpha=None, position_weights_by_length_beta=None, + query_group_ids_alpha=None, query_group_ids_beta=None, + reference_group_ids_alpha=None, reference_group_ids_beta=None, + group_distances_alpha=None, group_distances_beta=None, + query_position_weights_alpha=None, query_position_weights_beta=None, + reference_position_weights_alpha=None, reference_position_weights_beta=None): + """Return (counts, linear masses) for linked rows under summed lane distances. + + Shared matrix/gap costs/prior and separate full-coordinate lane weights and + group tables use the existing kernel. Sum distance is computed in int64; + thresholds remain int32. Count every score <= cutoff and sum max(0, cutoff + - score); ties count with zero mass. Negative cutoffs return zero. Preserve + duplicate reference observations and input threshold order. Exact exclusion + removes a row only when both encoded full strings match. One GIL-released + native traversal, O(Q*S) outputs and O(max_length + max_threshold_row) worker + scratch; no matrix/edges. The existing paired-max API is unchanged. + """ + qa, qb, ra, rb = map(list, (queries_a, queries_b, refs_a, refs_b)) + if len(qa) != len(qb) or len(ra) != len(rb): + raise ValueError("paired query/reference axes must agree") + thresholds = _threshold_rows(thresholds, len(qa)) + go, ge, cube, width, threads = _batch_options((qa, qb, ra, rb), matrix, gap_open, gap_extend, gap_prior, threads) + ag = _group_options(query_group_ids_alpha, reference_group_ids_alpha, group_distances_alpha, qa, ra) + bg = _group_options(query_group_ids_beta, reference_group_ids_beta, group_distances_beta, qb, rb) + aw = _sequence_weights(query_position_weights_alpha, reference_position_weights_alpha, position_weights_by_length_alpha, qa, ra) + bw = _sequence_weights(query_position_weights_beta, reference_position_weights_beta, position_weights_by_length_beta, qb, rb) + return _gapblock_paired_sum_count_batch(qa, qb, ra, rb, thresholds, alphabet, matrix, go, ge, + cube, width, threads, exclude_exact, + _position_weights(position_weights_by_length_alpha, qa, ra), + _position_weights(position_weights_by_length_beta, qb, rb), + ag[0], bg[0], ag[1], bg[1], ag[2], bg[2], + aw[0], bw[0], aw[1], bw[1]) + + class GapBlockIndex: """Search a reference set under the single-gap-block model, reusing the Hamming engine. diff --git a/python/seqtree/text_evalue.py b/python/seqtree/text_evalue.py new file mode 100644 index 0000000..b97e701 --- /dev/null +++ b/python/seqtree/text_evalue.py @@ -0,0 +1,142 @@ +"""Control-calibrated E-values for substring hits in a text. + +:func:`seqtree.evalues` answers "is this whole-sequence neighbourhood bigger than +background?" for an :class:`~seqtree.Index`. This module asks the same question of a +:class:`~seqtree.TextIndex`, where a short query is matched *inside* long records -- +a peptide against a proteome, a read against a transcript set. The arithmetic is +:func:`seqtree.evalue.evalue_result`; only the ball and the counting unit differ. + + E(q) = (N / M) * n_control(q) + p_any(q) = 1 - exp(-E) + p_enrichment = P(Poisson(E) >= n_target(q)) + +Three things are specific to text search, and each is a decision rather than a detail. + +**The ball is ``max_subs``, never ``score``.** ``TextResult.score`` is a similarity -- +``0`` for an exact match, more negative for a worse one -- the opposite sign to the +non-negative penalty :func:`seqtree.thetas_from_scores` expects. ``n_subs`` is a +non-negative integer, zero on an exact match, and monotone in the scope, so it is the +quantity that satisfies the nesting property the derivation rests on. + +**Counts are per distinct record, not per offset.** One query usually hits the same +record at several offsets -- a repeat, a low-complexity stretch, a tandem domain -- and +those events are strongly dependent. Counting them separately gives an over-dispersed, +compound-Poisson count and an E-value that is too large; collapsing to distinct records +avoids counting dependent offsets repeatedly. The Poisson approximation still +requires independent records and rare background hits. This is the same reason +``appendix/evalue.tex`` counts unique clonotypes rather than reads, and the same reason +it says to count the union of a dependent family rather than correct for its size. +``N`` and ``M`` are therefore ``num_refs``, the number of records in each text. + +**The control has to be a matched text.** The estimator assumes target and control are +drawn from one background: same residue composition, same record-length distribution, +same low-complexity content. A shuffled or reversed copy of the target text, or an +unrelated proteome of similar composition, is the usual choice. There is no +``load_control`` for texts -- that choice is the caller's, because it is the assumption +the whole calibration rests on, and nothing here can check it. +""" +from .evalue import evalue_result, thetas_from_scores + +__all__ = ["text_evalues", "text_threshold_for_evalue"] + + +def _search(index, queries, max_subs, threads): + """One native batch search and record reduction; refuse incomplete counts.""" + return index.search_batch(queries, max_subs=max_subs, threads=threads).record_distances() + + +def _counts(distances, exclude_exact): + """Distinct records per query, optionally dropping exact-match records.""" + return [sum(d > 0 for d in row) if exclude_exact else len(row) for row in distances] + + +def text_evalues(target, control, queries, max_subs=0, threads=0, exclude_exact=False): + """Control-calibrated E-values for each query's substring neighbourhood. + + Args: + target: :class:`~seqtree.TextIndex` to score hits in, e.g. a proteome. + control: background :class:`~seqtree.TextIndex`, matched to ``target`` in + composition and record-length distribution. See the module note -- this + choice is the assumption the calibration rests on. + queries: query strings, each at least ``target.k`` long. + max_subs: substitution radius defining the ball, searched identically in both + texts. + threads: worker threads for the batch searches (0 = all cores). + exclude_exact: drop 0-mismatch hits from both counts. Set this when a query may + occur verbatim in either text. A record with any exact occurrence is + excluded, even if it also contains inexact occurrences. + + Returns: + One dict per query with ``n_target``, ``n_control``, ``E``, ``p_any``, + ``p_enrichment`` and ``rule_of_three``. Counts are distinct records, not offsets. + + Raises: + ValueError: If the control text is empty, or if a search was truncated, which + would make a count a floor rather than a count. + + Example: + >>> import seqtree + >>> proteome = seqtree.TextIndex.build(["MKVLAAGIVGLNLGGK"], "aa", k=4) + >>> shuffled = seqtree.TextIndex.build(["GKLNVGALMGIVKGLA"], "aa", k=4) + >>> seqtree.text_evalues(proteome, shuffled, ["AAGIVGL"])[0]["n_target"] + 1 + """ + queries = list(queries) + if target.alphabet != control.alphabet: + raise ValueError("target and control alphabets must match") + n_ref, m_control = target.num_refs, control.num_refs + if m_control == 0: + raise ValueError("control text is empty") + + target_hits = _counts( + _search(target, queries, max_subs, threads), exclude_exact + ) + control_hits = _counts( + _search(control, queries, max_subs, threads), exclude_exact + ) + return [ + evalue_result(n_t, n_c, n_ref, m_control) + for n_t, n_c in zip(target_hits, control_hits) + ] + + +def text_threshold_for_evalue(target, control, queries, max_subs, e_target, threads=0, + exclude_exact=False): + """Per-query mismatch cutoff achieving ``e_target`` against ``control``. + + A fixed radius is not calibrated: at ``max_subs=2`` a low-complexity query picks up + thousands of chance occurrences while a compositionally unusual one picks up none, + so the radius that buys a given E-value has to be found per query. One control + search at ``max_subs`` supplies every cutoff. + + Args: + target: :class:`~seqtree.TextIndex` the E-value is expressed against. Only + ``num_refs`` is used. + control: background :class:`~seqtree.TextIndex`. + queries: query strings. + max_subs: radius the control is searched at, and the ceiling for the returned + cutoffs. + e_target: desired E-value, e.g. ``0.05``. + threads: worker threads for the batch search (0 = all cores). + exclude_exact: drop 0-mismatch hits, matching :func:`text_evalues`. + + Returns: + One integer cutoff per query in ``[0, max_subs]``, or ``-1`` where no radius + achieves ``e_target``. That happens when ``e_target < 3N/M``: with only ``M`` + control records, even an empty ball certifies no better than the rule-of-three + bound. Enlarge the control text. + + Raises: + ValueError: If ``max_subs`` is negative, or a search was truncated. + """ + if max_subs < 0: + raise ValueError("max_subs must be >= 0") + if target.alphabet != control.alphabet: + raise ValueError("target and control alphabets must match") + if control.num_refs == 0: + raise ValueError("control text is empty") + distances = _search(control, list(queries), max_subs, threads) + return thetas_from_scores( + distances, target.num_refs, control.num_refs, e_target, max_subs, + exclude_exact=exclude_exact, + ) diff --git a/setup.sh b/setup.sh index 4bb3fc6..ed91494 100755 --- a/setup.sh +++ b/setup.sh @@ -4,6 +4,7 @@ # --bench also install benchmark extras (huggingface_hub, psutil) # # Needs uv (https://docs.astral.sh/uv/); `brew install uv` or the standalone installer. +set -eu cd "$(dirname "$0")" || exit 1 if ! command -v uv >/dev/null 2>&1; then @@ -21,11 +22,12 @@ fi extras="" for arg in "$@"; do case "$arg" in - --tests) extras="[test]" ;; - --bench) extras="[bench]" ;; + --tests) extras="${extras}test," ;; + --bench) extras="${extras}bench," ;; esac done +if [ -n "$extras" ]; then extras="[${extras%,}]"; fi uv pip install -e ".$extras" echo "" diff --git a/skills/seqtree/SKILL.md b/skills/seqtree/SKILL.md deleted file mode 100644 index e455836..0000000 --- a/skills/seqtree/SKILL.md +++ /dev/null @@ -1,364 +0,0 @@ ---- -name: seqtree -description: Fast fuzzy search over biological sequences — C++20 arena trie + nanobind, with control-calibrated E-values, single-gap-block alignment, and seed significance. ---- - -# seqtree - -Payload-agnostic fuzzy sequence search. Build an immutable index once; search single queries or -batches in parallel. Everything downstream (`vdjmatch`, `mhcmatch`) maps `ref_id` back to its own -payload — seqtree ships **zero runtime dependencies** and no domain data. - -Repo: `/Users/mikesh/vcs/code/seqtree`. Venv: `.venv` (`bash setup.sh`). Docs build: -`env -C docs make html` (must be warning-free). - -## Engines — pick with care - -| engine | enforces `max_subs`/`max_ins`/`max_dels` | honours `max_penalty` | use when | -|---|---|---|---| -| `seqtm` | **yes** | yes | always, unless you specifically want budget-only | -| `seqtrie` | **no** | yes | the score budget is the entire specification | - -`engine="auto"` resolves to `seqtm`. It does **not** choose per query; `seqtrie` silently ignores -edit caps, so `auto` can never select it. The trap: `SearchParams(matrix=..., max_penalty=12)` -without `engine="seqtrie"` still runs `seqtm`, where `max_subs` is 0 by default -- exact matches -only, and nothing in the result says why. Passing a matrix to `seqtrie` without an explicit -`max_penalty` raises rather than scanning the whole table. - -## Core API - -```python -Index.build(refs: list[str], alphabet: str = "aa") -> Index # 'aa' | 'nt' | 'iupac' -Index.ref_seqs() -> list[str] # one native call, ref_id order -Index.search(query, params) -> list[Hit] # Hit(ref_id, score, n_subs, n_ins, n_dels) -Index.search_top(query, params, k=1) -> list[Hit] -Index.search_batch(queries, params, threads=0) -> list[list[Hit]] # releases the GIL -Index.align(ref_id, query, params) -> Alignment # true Gotoh affine, on demand -Index.save(path) / Index.load(path) -pairwise_batch(a, b, params, alphabet="aa", threads=0) -> list[list[Hit]] # a-major - -SearchParams(max_subs=0, max_ins=0, max_dels=0, max_total_edits=0, max_penalty=0, - matrix="", gap_open=1, gap_extend=1, engine="auto", mode="all") -SearchParams.pos_matrix = PositionalMatrix | None -``` - -**`pairwise_batch` is a bounded trie search, not an aligner.** It returns nothing outside -`max_penalty`. The only affine aligner is `Index.align`. - -## Substitution matrices - -```python -SubstitutionMatrix.blosum62() / .blosum45() / .blosum80() / .pam250() / .pam100() - / .structural() / .unit(n) -SubstitutionMatrix.from_similarity(grid: list[list[int]]) # Gram: pen = s_aa + s_bb - 2*s_ab, clamped >= 0 -SubstitutionMatrix.penalty(a, b) -> int # pen(a, a) == 0 always -SubstitutionMatrix.similarity(a, b) -> int # raw signed log-odds -SubstitutionMatrix.scale() -> int # median mismatch; BLOSUM62 == 14 -PositionalMatrix.from_weights(base, weights) / .from_tables(size, width, data, masked=[]) -``` - -`scale()` exists because gap costs must live on the matrix's scale. `gap_open=1` against a -BLOSUM62 mismatch of 14 makes gaps ~14× cheaper than substitutions and every alignment -degenerates to gaps. **Use `gap_open = 2 * matrix.scale()`.** - -Because `pen(a, a) == 0`, a positional weight scales *mismatch* cost only. It carries **zero -match evidence** — evidence lives in the control-counted E-value, never in a weight. - -## E-values and calibrated cutoffs — `seqtree.evalue` - -```python -evalues(target, control, queries, params, threads=0, exclude_exact=False) -> list[dict] - # keys: n_target, n_control, E, p_any, p_enrichment, rule_of_three - # E = (N/M) * n_control ; rule of three (3N/M) when the control ball is empty -evalue_result(n_target, n_control, n_ref, m_control) -> dict # that record, for counts you - # already have. What `evalues` and `pmhc_evalue.homolog_evalue` both call. - -threshold_for_evalue(target, control, queries, params, e_target, threads=0, - exclude_exact=False) -> list[int] # per query; -1 == unreachable -thetas_from_scores(control_scores, n_target, m_control, e_target, theta_max, - *, exclude_exact=False) -> list[int] # the pure-data core -load_control(name="human_trb_aa", size=None) -> Index # 250k bundled; larger via HuggingFace - # Cold-cache safe under a multi-process fan-out (xdist / Snakemake / Nextflow): Index.save - # writes a temp file and renames it into place, so a reader never sees a partial index. - # A corrupt/stale cache is rebuilt, not raised. filelock (optional) additionally dedupes - # the build work; correctness does not depend on it. -``` - -**A fixed score cutoff is not a calibrated cutoff.** The control is dense near germline and -sparse among rare junctions. Measured on human TRB: at `gapblock_score <= 60`, **31.7%** of random -control junctions land in a connected component of ≥5 — structure the threshold invents. Per-query -cutoffs cut that to 0.000 while raising the real signal to 2.334 edges/node. Always invert an -E-value; never pick θ by hand. `-1` means `e_target < 3N/M` — enlarge the control rather than -lowering the bar. - -`threshold_for_evalue` needs `params.max_penalty > 0` as the ceiling to search the control at. -One control scan supplies every query's cutoff. - -## Pairwise alignment — `seqtree.pairwise` (replaces BioPython) - -Everything else here MINIMISES a non-negative penalty. This module MAXIMISES a raw log-odds -similarity, like BLAST/BioPython, because that is what pairwise alignment means. - -```python -from seqtree.pairwise import score, align, score_matrix, dist_matrix -score(q, r, matrix, mode="global", gap_open=11, gap_extend=1) -> int # NW; "local" = SW -align(...) -> Alignment # + aligned strings/ops; .score is a SIMILARITY, not a penalty -score_matrix(queries, refs, ...) -> ScoreMatrix # dense n x K, GIL released -dist_matrix(queries, refs, ...) -> ScoreMatrix # d = s(a,a)+s(b,b)-2s(a,b), 0 on diagonal -``` - -- Gap run of length L costs `gap_open + (L-1)*gap_extend`. **`gap_open == gap_extend` = linear - gaps** (no separate mode). Global charges end gaps (true NW, not semi-global). Local floors at 0. -- **Drop-in for `Bio.Align.PairwiseAligner`**: 0 disagreements over 3 matrices x 10 gap/mode - settings x 60 shapes, and on real germline V genes. **65-87x faster.** BioPython is a TEST-only - dep; seqtree still has zero runtime deps. -- `SubstitutionMatrix` now stores BOTH views: `penalty()` (>=0, Gram, for search/E-values) and - `similarity()` (signed log-odds, for these aligners). The Gram transform is lossy -- it zeroes - the diagonal -- so a similarity cannot be recovered from a penalty. - -## Plain distances and Hamming balls — `seqtree.distance` - -```python -hamming(a, b) -> int # equal length only; RAISES otherwise. Case-SENSITIVE. -levenshtein(a, b) -> int -hamming_matrix(a, b, threads=0) / levenshtein_matrix(...) -> ScoreMatrix # GIL released - -neighbourhood(seq, r=1, alphabet=None, include_self=True, shell=False) -> list[str] -neighbourhood_union(seqs, r=1, ...) -> list[str] # dedup DURING the walk -union_size(seqs, r=1, ...) -> int # size the job first -``` - -Enumeration is **substitution only, fixed length** — no indels, because Hamming is undefined -across lengths and `P_gen` is length-conditioned. `alphabet=None` is the 20 standard residues, -which is `amino_acids()` **minus** `B`/`Z`/`X`/`*` — passing `amino_acids()` itself gives 23·L -neighbours, not 19·L. `shell=True` is the sphere at exactly `r` from the *nearest* centre; -`include_self=False` drops the whole `r = 0` shell, i.e. every centre. - -**Take the union, never the sum of the balls.** 200 length-14 junctions drawn within distance `d` -of a common centre, `r = 1`: at `d = 1` the per-sequence balls double-count **41.7%** -(53,400 → 31,122), at `d = 2` 4.5%, at `d ≥ 3` nothing. 300 junctions of length 14 at `r = 1` is -80,100 distinct sequences in **23 ms** on one M3 core (`union_size` 11 ms), so this is pure -Python and does not need the C++ core. `r = 2` over the same 300 is 9.9 M sequences, 6.8 s, -~1.8 GB — there the dedup is what keeps it to one copy. - -## Text search — `seqtree.TextIndex` - -```python -ix = TextIndex.build(records, alphabet="aa", k=4, group_ids=None) # WHOLE records, not windows -res = ix.search_batch(queries, max_subs=2, max_indels=0, exclude_exact=False, best_only=False, - group_by=False, max_hits=0, matrix=None, threads=0) # GIL released - -len(res); res.num_hits; res[i] # TextHit list for query i, built on demand -res[i][0].mismatches # [(pos, query_aa, text_aa), ...] -- the PAIR -res[i][0].n_ins / .n_dels / .length # indels, and the width of the match in the text -res.groups(i) # [(group_id, min_subs, n_hits), ...] with group_by -res.truncated # per query, 1 if max_hits capped it -res.arrays(); res.to_numpy() # zero-copy views over the flat CSR -ix.save(path); TextIndex.load(path, mmap=True) -``` - -**Use this, not `Index`, when the reference is a text.** `Index` is a trie over reference -*strings*, so a text question means enumerating every length-`L` window as its own string — and -again for every distinct query length. The human proteome has 68,389,335 nine-mer windows; -a query set spanning 45 lengths costs 45 multi-gigabyte builds, which is what ran > 2 h without -finishing and filled 225 GB of cache in `mhcmatch`. Here **`k` belongs to the index**: one -build answers every length and every `max_subs`. - -Exact, not heuristic, by **one search scheme** over the single seed table. Split the query into -`b = min(max_subs + 1, L // k)` disjoint blocks of width `L // b`, and probe block `j`'s leading -k-mer at radius `c_j`. It is **lossless exactly when `sum(c_j) >= max_subs - b + 1`** (the -lightest error vector that could hide is `e_j = c_j + 1`, of weight `sum(c_j) + b`), and the -cheapest legal choice spreads `r = max(0, max_subs - b + 1)` budget units evenly, onto the -*last* blocks. `b = max_subs + 1` is plain pigeonhole; `b = 1` is a single ball over `q[0:k]`; -the useful schemes are in between. Completeness is pinned by brute-force **set equality** over -L 6–30 x `max_subs` 0–3 x k ∈ {3,4,5}, and by the answer being identical across `k`. - -**Choosing `k`: the crossover is `L = 2k`.** Below it only one block fits and the query pays a -full radius-`max_subs` ball, which is one to two orders of magnitude dearer. So the default -`k = 4` is right down to length 8, and `k = 5` is worth a second index when the query set starts -at 10 (human proteome, `L = 12, max_subs = 3`: **1.46 ms/query at k=4, 0.32 at k=5**). Building -is ~0.6 s either way, so a corpus spanning both lengths can hold both indexes. - -**`max_indels` allows gaps, and costs nothing in the index.** The two caps are independent: -`max_subs=2, max_indels=1` is two substitutions AND one gap, not three edits. Seeding is -unchanged -- with `b = max_subs + max_indels + 1` blocks the edit count is below the block -count, so pigeonhole still leaves a **zero-error** block, which matches exactly whether or not -the occurrence has gaps. Only verification changes, to a banded DP over the -`2*max_indels + 1` possible starts. Three things to know: a match is **no longer `len(query)` -wide** (read `.length`, which is `len(query) + n_dels - n_ins`); every block must be exact, so a -query must be at least **`(max_subs + max_indels + 1) * k`** long or it raises; and -`.mismatches` is empty and `matrix=` ignored, because the per-column detail would need a -traceback the verifier does not keep. Costs 25-37x an ungapped search at the same length, still -under 0.1 ms/query on 5 M residues. Pinned by set equality against brute force over -k in {3,4,5} x max_subs 0-3 x max_indels 1-2, 6,732 occurrences, zero missing and zero extra. - -- **`matrix=` only scores** hits the Hamming predicate already accepted (similarity summed over - the mismatched positions) — it never changes which hits return. -- **`best_only=True`** walks the distance upward and stops at the first non-empty shell, - returning **all** of it. No second index, unlike PEPMatch's log2(L) databases. -- **`group_ids` are arbitrary** — gene ids, species, clusters. `group_by=True` folds hits onto - them and a second row is a **tie**, which is what a caller refusing ambiguous parents needs. -- **`max_hits` caps after the sort**, so the best hits survive, and always sets `truncated`. -- Order is `(n_subs, ref_id, offset)`, stable across runs and thread counts. - -## Gap-block alignment — `seqtree.gapblock` - -A V(D)J junction's length variation is **one** contiguous indel event, so restrict alignment to -one gap block of length `d = |m - n|`. - -```python -gapblock_score(q, r, matrix=None, gap_open=None, gap_extend=1, gap_prior=None) -> (score, block_pos) -gap_cost(d, gap_open, gap_extend) -> int # 0 when d == 0 -- MANDATORY guard -deletion_variants(q, d) -> list[(block_pos, variant)] -GapBlockIndex(refs, alphabet="aa", d_max=1).search(query, max_penalty, matrix, gap_open, - gap_extend=1, gap_prior=None) - # -> list[(ref_id, score, block_len, block_pos)] - -GapPrior = Callable[[block_start_i, block_length_d, longer_length_m], int] -central_prior(lam) # lam * |block_midpoint - m/2|; lam ~ 1.5 * scale() -profile_prior(lam, w) # lam * sum(w(j, m) for j in block); w in [0,1], callable -frame_prior(lam, c) # lam * |i - c|; the block start does not depend on d -positions_prior(starts) # 0 at each start (negatives from the end), UNREACHABLE elsewhere -embed_in_frame(seq, width, c, gap="-") -> str - -# Dense n x K, for prototype-distance embeddings. C++, GIL released, buffer protocol. -score_matrix(queries, refs, matrix=None, gap_open=None, gap_extend=1, gap_prior=None, - alphabet="aa", threads=0) -> ScoreMatrix # numpy.asarray(sm) is zero-copy - -# Per-island PWM. pen(j,a) = lam*log(p_max_j / p_j(a)): >= 0 and 0 on the consensus, so it is -# still a ball and still feeds thetas_from_scores. c defaults to the entropy-optimal column. -IslandProfile.fit(members, c=None, lam=1000, pseudocount=0.5) -> IslandProfile - .score(seq) -> int # UNREACHABLE if it does not fit the frame - .score_batch(seqs) -> list[int] - .consensus() -> str # score(consensus()) == 0 -``` - -**Invariants a prior must satisfy** (nothing else): `>= 0`, and `== 0` when `d == 0`. Otherwise -`s(q, q) != 0` and the score stops defining a ball, breaking the E-value theory and admissible -trie pruning. Monotonicity in `d` is **not** required and is false for `central_prior`. - -### What the measurements say - -- **One block is enough.** Against a model-independent structural oracle (iterative superposition - + unrestricted affine DP) over 3,049 crystal junction pairs, the true correspondence is a single - contiguous block in **95.2–100%** of cases for every `d = 1..4`. Forcing one block costs no - median CA-RMSD. -- **The restriction is free on related sequences, and protective on unrelated ones.** At - `gap_open = 2*scale`, gap-block equals unrestricted affine on **98.8%** of related pairs; on - unrelated pairs affine undercuts it by a median of 106 penalty units, i.e. affine *invents* an - alignment. Free gapping buys manufactured similarity. -- **A sequence score cannot place the block.** A hard central pin agrees with the flat, - score-only choice on only 10.6% of pairs. The block sits at the **loop apex** (Cys-offset 6 for - both TRA and TRB) and does not drift with `d`. `central_prior` hits it 42.4% (TRA) / 30.1% - (TRB); a germline-untemplated-span rule hits 0.4% / 19.8% and was rejected. -- **A frame is transitive iff the block start is constant in `d`.** That is `frame_prior`. - `central_prior`'s start drifts, so embedding two shorter members into a common frame relates - them by *two* blocks — no consistent column index, no PWM. Use `embed_in_frame(seq, W, c)`: - left-anchor the first `c` residues, right-anchor the rest. -- **Scoring several candidate placements and keeping the best loses.** At a matched FPR on human - TRB, candidate starts `(3, 4, mid)` reach precision 0.156; a hard-pinned centre 0.414. That is - what `positions_prior` implements, and it exists to reproduce other aligners' conventions — - `mir.distances.aligner.JunctionAligner` hardcodes `(3, 4, -4, -3)` for all seven loci — not - because it is the right rule. - -- **A per-island PWM beats min-over-members only at a strict cutoff, and "strict" depends on `N`.** - The E-value's `k = floor(E*·M/N)` is how many control neighbours the cutoff admits, so FPR = `k/M`. - 108 calibrated VDJdb islands >= 10, held-out members, 250k control negatives, paired bootstrap - over islands: - - | regime | FPR | min | PWM | diff [95% CI] | - |---|---|---|---|---| - | loose reference | 1% | **99.5%** | 99.1% | −0.40 [−1.09, +0.14] | - | per-epitope islands (`N` = group size, median 88) | 0.0568% | 88.3% | **89.3%** | +0.93 [−0.80, +2.79] | - | repertoire annotation (`N` ≈ 20k) | 0.0012% | 37.6% | **48.5%** | +10.90 [+7.69, +14.21] | - - No significant difference *while building* islands; a large one when annotating a repertoire with - them (islands >= 50: 9.8% vs 22.6%). Note `E*=0.05` at `N≈20k` gives `k=0` → `thetas_from_scores` - returns `-1`; the rule of three forbids any `E < 3N/M = 0.236`, and that is the third row's cutoff. - Neither representation generalises to same-epitope junctions in a *different* island (3.7% vs 3.5% - at 1% FPR, ~0% at either operating point). - -### Batch scoring — `score_matrix` - -A search prunes; an embedding cannot, because the distance to every prototype *is* the output. -Measured on an M3, 3,000 refs, human TRB: pure-Python `gapblock_score` **0.41 M pairs/s**, -`score_matrix` **51.3 M** (1 thread) and **532.7 M** (16 threads). The prior costs nothing — it -is flattened once into an `[m][d][i]` cube. Reproduce with `bench/bench_score_matrix.py`. - -The result *is* the distance: the Gram transform is applied per residue when the matrix is built, -so there is no `d = s(a,a) + s(b,b) - 2·s(a,b)` step, and non-negativity, symmetry and a zero -diagonal hold by construction. Budget `4 * n * K` bytes and chunk the queries. - -### Bounded exhaustive reductions - -`gapblock.topk_batch(queries, refs, k=10, matrix=None, gap_open=None, gap_extend=1, -gap_prior=None, alphabet="aa", threads=0, exclude_exact=False)` returns nested Hit -rows sorted by `(score, ref_id)`, scoring every reference with the same native -single-gap-block cell. Only score/ref_id are provided; edit fields remain unset. -Results use O(Q*k), worker scratch O(k+maxlen); encoded inputs/prior cube shared. -Exact encoded identity is excluded before selection; zero-cost nonidentity remains. - -`gapblock.count_batch(queries, refs, thresholds, **same_scoring_args)` returns -nested integer counts for scores<=each threshold, including ALL ties. Thresholds -have one row/query, arbitrary order/duplicates/empty rows allowed, output order -preserved, negative cutoffs count zero. No dense Q*N allocation. - -`gapblock.paired_topk_batch(queries_a, queries_b, refs_a, refs_b, k=10, -**same_scoring_args)` ranks linked rows by `(max(score_a,score_b),ref_id)` and -returns nested `(ref_id,score_a,score_b)` tuples. Matching lane axis lengths required. -`paired_count_batch(..., thresholds, **same_scoring_args)` counts all linked rows -at those maximum-score cutoffs. Exclusion removes only both-view identity. -These are generic exhaustive score reductions, not statistical calibration. - -### Performance - -Profiled over the bundled 250k control: 91% of query time is the query-deletion-variant branch, -9% the auxiliary indices (9.8M entries at `d_max=3`). Netting the prior out of each variant's -budget cuts that branch from ~15 sub-searches to 2.5. **Variant dedup and length-bucketing are -not worth building** — measured, not guessed. At budget 40 over 250k refs, `d_max=2` gap-block -search costs 2,562 µs/query, *less* than the plain Hamming ball at the same budget. - -## Seed E-values — `seqtree.seeds` - -```python -core_kmers(seq, k, flank=4) -> set[str] # drops the germline-framed ends -SeedIndex(seqs, k=5, flank=4) / SeedIndex.from_index(index, k, flank) - .count(seed) / .evalue(seed, n_target) / .seed_evalues(query, n_target) - .significant(query, n_target, alpha=1.0) / .union_evalue(...) / .gather(query, seeds=None) -``` - -Seeds buy **precision, not recall**: a shared rare central k-mer is real evidence (≈4× enriched), -but only ~0.5% of co-specific cross-island pairs share one. Never assume a shared k-mer is -significant — a D-gene run like `LAGG` sits in 3,471/250,000 control sequences (`E_seed` = 1,388). -Count it in the control. - -## Layout / pMHC — `seqtree.layout`, `seqtree.pmhc` - -`AnchorSpec`, `mask_anchors`, `kmers`, `presentation_features`, `weight_profile` (a -*mismatch-tolerance* profile, not information weighting), `PMHCStore`, `find_mimics`. - -## Gotchas - -1. `gap_open` defaults to **1**, which is wrong for any real matrix. Pass `2 * matrix.scale()`. -2. `Index.build` rejects non-alphabet characters — filter `B/Z/X/*` out of CDR3s first. -3. The bundled control is a **uniform reservoir sample** of unique productive clonotypes, shuffled, - so `bundled[:k]` is a valid sub-sample. It used to be the abundance head (25.8x more - self-similar, 3.1x the ball mass); if you have a cached `~/.cache/seqtree/control_*.sqtree` - from before that fix, delete it. -4. OLGA generates from `P_gen`, which is **not** the post-selection background `P_0`. Do not use - an OLGA sample as the null for an operational cutoff. -5. Rebuild after pulling: a stale `_core.so` has silently shadowed new Python-visible C++ methods - more than once (`bash setup.sh`, or `uv pip install -e .`). -6. `TextIndex` treats the two sides differently on purpose. A residue outside the alphabet in the - **text** (36 `U` in the human proteome, 33 in mouse) becomes a hole — no hit crosses it — - and is counted in `ix.num_unknown`; the same residue in a **query** raises. `B`/`Z`/`X`/`*` - are real symbols, not wildcards: `X` matches only `X`. And every query must be at least `k` - long, or it raises rather than being answered incompletely — as does one over 65,535 long, - since mismatch positions are 16-bit. `ref_seq` renders a hole as `'X'`, so it reports what - the index holds, not what was handed to `build`. - -## Layering - -seqtree stays generic and dependency-free. TCR germline data, clustering, structures and -benchmarks against VDJdb live downstream in `vdjmatch` / `tcren-ms`. Do not import domain data -into `bench/`. diff --git a/src/_bindings.cpp b/src/_bindings.cpp index 2052c76..cc57712 100644 --- a/src/_bindings.cpp +++ b/src/_bindings.cpp @@ -186,6 +186,17 @@ nb::list py_search_batch(const Index& idx, const std::vector& queri return out; } +std::vector> py_edit_histogram_batch( + const Index& idx, const std::vector& queries, + const PyParams& pp, int threads, bool exclude_exact) { + if (pp.max_total_edits < 1 || pp.max_total_edits > 65535) + throw nb::value_error("edit histogram requires max_total_edits in [1, 65535]"); + auto mat = make_matrix(pp, idx.alphabet()); + SearchParams cp = to_cpp(pp, mat ? &*mat : nullptr); + nb::gil_scoped_release release; + return idx.edit_histogram_batch(queries, cp, threads, exclude_exact); +} + std::vector py_collisions_batch(const Index& idx, const std::vector& queries, const PyParams& pp, int threads) { auto mat = make_matrix(pp, idx.alphabet()); @@ -349,7 +360,8 @@ ScoreMatrix py_gapblock_matrix(const std::vector& queries, const std::vector& refs, const std::string& alphabet, const std::optional& matrix, int32_t gap_open, int32_t gap_extend, const std::vector& prior, - uint32_t prior_width, int threads) { + uint32_t prior_width, int threads, + const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { Alphabet alph = parse_alphabet(alphabet); ScoreMatrix out; out.rows = queries.size(); @@ -357,7 +369,7 @@ ScoreMatrix py_gapblock_matrix(const std::vector& queries, { nb::gil_scoped_release release; out.data = gapblock_matrix(queries, refs, alph, matrix ? &*matrix : nullptr, gap_open, - gap_extend, prior, prior_width, threads); + gap_extend, prior, prior_width, threads, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); } return out; } @@ -664,6 +676,13 @@ NB_MODULE(_core, m) { nb::arg("threads") = 0, "Search many queries in parallel (releases the GIL). ``threads=0`` uses all " "cores. Returns one hit list per query, in input order.") + .def("edit_histogram_batch", &py_edit_histogram_batch, + nb::arg("queries"), nb::arg("params"), nb::arg("threads") = 0, + nb::arg("exclude_exact") = false, + "Exact reference counts by total n_subs+n_ins+n_dels, bins 0..max_total_edits. " + "Requires explicit positive max_total_edits and uncapped seqtm/auto all-hit search. " + "Duplicates count separately; exclude_exact retains bin 0 with count zero. threads<=0 uses the detected available CPU budget. Releases the GIL; " + "one reusable hit buffer per worker, no all-query hit vectors.") .def("align", &py_align, nb::arg("ref_id"), nb::arg("query"), nb::arg("params"), "Compute a global alignment between ``query`` and a reference, on demand.") .def("collisions_batch", &py_collisions_batch, nb::arg("queries"), nb::arg("params"), @@ -781,7 +800,8 @@ NB_MODULE(_core, m) { m.def("gapblock_matrix", &py_gapblock_matrix, nb::arg("queries"), nb::arg("refs"), nb::arg("alphabet") = "aa", nb::arg("matrix") = std::nullopt, nb::arg("gap_open") = 1, nb::arg("gap_extend") = 1, nb::arg("prior") = std::vector{}, - nb::arg("prior_width") = 0, nb::arg("threads") = 0, + nb::arg("prior_width") = 0, nb::arg("threads") = 0, nb::arg("position_weights_by_length") = std::vector>{}, + nb::arg("query_group_ids") = std::vector{}, nb::arg("reference_group_ids") = std::vector{}, nb::arg("group_distances") = std::vector>{}, nb::arg("query_position_weights") = std::vector>{}, nb::arg("reference_position_weights") = std::vector>{}, "Exhaustive single-gap-block penalties for every (query, ref) pair, GIL released. " "`prior` is the gap prior flattened to [m][d][i]; see seqtree.gapblock.score_matrix, " "which builds it for you."); @@ -792,11 +812,29 @@ NB_MODULE(_core, m) { return gapblock_topk_batch(queries, refs, k, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact); }, nb::arg("queries"), nb::arg("refs"), nb::arg("k"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, "Bounded exhaustive gapblock reduction, GIL released."); - m.def("gapblock_count_batch", [](const std::vector& queries, const std::vector& refs, const std::vector>& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact) { + m.def("gapblock_count_batch", [](const std::vector& queries, const std::vector& refs, const std::vector>& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { const auto alph = parse_alphabet(alphabet); nb::gil_scoped_release release; - return gapblock_count_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact); - }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, "Bounded exhaustive gapblock reduction, GIL released."); + return gapblock_count_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Bounded exhaustive gapblock reduction, GIL released."); + + m.def("gapblock_count_mass_batch", [](const std::vector& queries, const std::vector& refs, const std::vector>& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto alph = parse_alphabet(alphabet); + nb::gil_scoped_release release; + return gapblock_count_mass_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Counts and integer linear mass in one bounded traversal, GIL released."); + + m.def("gapblock_group_audit_batch", [](const std::vector& queries, const std::vector& refs, const std::vector& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto alph = parse_alphabet(alphabet); + nb::gil_scoped_release release; + return gapblock_group_audit_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Streaming component/full ball counts and minimum full score within the component ball, GIL released."); + + m.def("gapblock_ball_batch", [](const std::vector& queries, const std::vector& refs, const std::vector& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights, const std::vector& query_groups, const std::vector& reference_groups, const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto alph = parse_alphabet(alphabet); + nb::gil_scoped_release release; + return gapblock_ball_batch(queries, refs, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + }, nb::arg("queries"), nb::arg("refs"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length")=std::vector>{}, nb::arg("query_group_ids")=std::vector{}, nb::arg("reference_group_ids")=std::vector{}, nb::arg("group_distances")=std::vector>{}, nb::arg("query_position_weights")=std::vector>{}, nb::arg("reference_position_weights")=std::vector>{}, "Bounded exhaustive gapblock reduction, GIL released."); m.def("gapblock_paired_topk_batch", [](const std::vector& qa, const std::vector& qb, const std::vector& ra, const std::vector& rb, uint32_t k, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact) { const auto alph = parse_alphabet(alphabet); @@ -810,6 +848,12 @@ NB_MODULE(_core, m) { return gapblock_paired_count_batch(qa, qb, ra, rb, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact); }, nb::arg("qa"), nb::arg("qb"), nb::arg("ra"), nb::arg("rb"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, "Bounded exhaustive gapblock reduction, GIL released."); + m.def("gapblock_paired_sum_count_batch", [](const std::vector& qa, const std::vector& qb, const std::vector& ra, const std::vector& rb, const std::vector>& thresholds, const std::string& alphabet, const std::optional& matrix, int32_t go, int32_t ge, const std::vector& prior, uint32_t width, int threads, bool exclude_exact, const std::vector>& weights_alpha, const std::vector>& weights_beta, const std::vector& query_groups_alpha, const std::vector& query_groups_beta, const std::vector& reference_groups_alpha, const std::vector& reference_groups_beta, const std::vector>& group_distances_alpha, const std::vector>& group_distances_beta, const std::vector>& query_weights_alpha, const std::vector>& query_weights_beta, const std::vector>& reference_weights_alpha, const std::vector>& reference_weights_beta) { + const auto alph = parse_alphabet(alphabet); + nb::gil_scoped_release release; + return gapblock_paired_sum_count_batch(qa, qb, ra, rb, thresholds, alph, matrix ? &*matrix : nullptr, go, ge, prior, width, threads, exclude_exact, weights_alpha, weights_beta, query_groups_alpha, query_groups_beta, reference_groups_alpha, reference_groups_beta, group_distances_alpha, group_distances_beta, query_weights_alpha, query_weights_beta, reference_weights_alpha, reference_weights_beta); + }, nb::arg("qa"), nb::arg("qb"), nb::arg("ra"), nb::arg("rb"), nb::arg("thresholds"), nb::arg("alphabet")="aa", nb::arg("matrix")=std::nullopt, nb::arg("gap_open")=1, nb::arg("gap_extend")=1, nb::arg("prior")=std::vector{}, nb::arg("prior_width")=0, nb::arg("threads")=0, nb::arg("exclude_exact")=false, nb::arg("position_weights_by_length_alpha")=std::vector>{}, nb::arg("position_weights_by_length_beta")=std::vector>{}, nb::arg("query_group_ids_alpha")=std::vector{}, nb::arg("query_group_ids_beta")=std::vector{}, nb::arg("reference_group_ids_alpha")=std::vector{}, nb::arg("reference_group_ids_beta")=std::vector{}, nb::arg("group_distances_alpha")=std::vector>{}, nb::arg("group_distances_beta")=std::vector>{}, nb::arg("query_position_weights_alpha")=std::vector>{}, nb::arg("query_position_weights_beta")=std::vector>{}, nb::arg("reference_position_weights_alpha")=std::vector>{}, nb::arg("reference_position_weights_beta")=std::vector>{}, "Linked sum counts and linear mass in one bounded traversal, GIL released."); + m.def("hamming", &hamming, nb::arg("a"), nb::arg("b"), "Hamming distance: the number of positions at which two EQUAL-length sequences differ. " "Raises ValueError on a length mismatch. Case-sensitive, byte-for-byte."); @@ -992,6 +1036,10 @@ NB_MODULE(_core, m) { .def("__len__", [](const TextResult& r) { return r.num_queries(); }) .def_prop_ro("num_hits", [](const TextResult& r) { return r.num_hits(); }, "Total hits across every query.") + .def("record_distances", &TextResult::record_distances, + nb::call_guard(), + "Sorted best substitution distances per distinct record for each query. " + "Reduces offsets in one native call; rejects truncated results.") .def( "__getitem__", [](const TextResult& r, Py_ssize_t i) { diff --git a/src/count_unit.hpp b/src/count_unit.hpp new file mode 100644 index 0000000..2b01f6e --- /dev/null +++ b/src/count_unit.hpp @@ -0,0 +1,87 @@ +// Internal exact unit histogram acceleration for the bounded unit-cost profile. +#pragma once +#include "seqtree/seqtree.hpp" +#include +#include +#include +#include +#include + +namespace seqtree::count_unit { +struct Pattern { + std::array eq{}; + std::array codes{}; + size_t length; + Pattern(std::string_view query, const Codec& codec):length(query.size()) { + if (!length || length>63) throw std::invalid_argument("pattern length must be 1..63"); + for (size_t i=0;iint(r.size()) || s+ins+del>cap) continue; + if (i==int(q.length) && j==int(r.size())) best=std::min(best,s+ins+del); + auto put=[&](int ii,int ni,int nd,int ns) { + if (ns+ni+nd<=cap) dp[ii][ni][nd]=std::min(dp[ii][ni][nd],uint8_t(ns)); + }; + if (iq.length+2 || q.length>r.size()+2) return cap+1; + const int delta=int(r.size())-int(q.length); + const int d=q.distance(r); + if (d>cap) return cap+1; + // For any path, ins+del<=d and |ins-del|=|delta|. Thus max(ins,del) + // <=floor((d+|delta|)/2); both per-type caps are automatic when <=2. + if (d+std::abs(delta)<=5) return d; + return restricted(q,r,codec,cap); +} +// Bounded installed-API profiles favor trie traversal at total caps1/2. +inline bool eligible(const SearchParams& p) { + return !p.matrix && !p.pos_matrix && p.gap_open==1 && + p.max_total_edits>=3 && p.max_total_edits<=5 && + p.max_substitutions>=p.max_total_edits && p.max_insertions==2 && p.max_deletions==2 && + (p.max_score_penalty<=0 || p.max_score_penalty>=p.max_total_edits) && + p.engine!=Engine::SeqTrie && p.mode==Mode::AllHits && !p.max_hits; +} +inline bool fill_histogram(const Index& idx,std::string_view query,const SearchParams& p, + bool exclude_exact,std::vector& out) { + if (!eligible(p) || query.empty() || query.size()>63) return false; + Pattern pattern(query,idx.codec()); + for(uint32_t id=0;id pen_table(const SubstitutionMatrix* m, uint8_t A) { // contractually zero -- hence the early return, which also keeps s(q, q) == 0. inline int32_t cell(const uint8_t* q, uint32_t m, const uint8_t* r, uint32_t n, const int32_t* pen, uint8_t A, int32_t gap_open, int32_t gap_extend, - const int32_t* prior_row, int32_t* suf) { + const int32_t* prior_row, int32_t* suf, const int32_t* weights = nullptr, + const int32_t* qw = nullptr, const int32_t* rw = nullptr) { const uint32_t L = std::min(m, n); const uint32_t d = (m > n ? m - n : n - m); const bool q_longer = m >= n; @@ -47,7 +48,7 @@ inline int32_t cell(const uint8_t* q, uint32_t m, const uint8_t* r, uint32_t n, for (uint32_t j = L; j-- > 0;) { const uint8_t a = q_longer ? q[j + d] : q[j]; const uint8_t b = q_longer ? r[j] : r[j + d]; - suf[j] = suf[j + 1] + pen[size_t(a) * A + b]; + suf[j] = suf[j + 1] + pen[size_t(a) * A + b] * (qw && rw ? std::max(qw[j + (q_longer ? d : 0)], rw[j + (q_longer ? 0 : d)]) : (weights ? weights[j + d] : 1)); } if (d == 0) return suf[0]; @@ -56,7 +57,7 @@ inline int32_t cell(const uint8_t* q, uint32_t m, const uint8_t* r, uint32_t n, for (uint32_t i = 0; i <= L; ++i) { const int32_t cand = pre + suf[i] + (prior_row ? prior_row[i] : 0); if (cand < best) best = cand; - if (i < L) pre += pen[size_t(q[i]) * A + r[i]]; + if (i < L) pre += pen[size_t(q[i]) * A + r[i]] * (qw && rw ? std::max(qw[i], rw[i]) : (weights ? weights[i] : 1)); } return best + gap_open + int32_t(d - 1) * gap_extend; } @@ -66,6 +67,9 @@ inline int32_t cell(const uint8_t* q, uint32_t m, const uint8_t* r, uint32_t n, struct GapInput { std::vector> q, r; std::vector pen; + std::vector> weights, query_weights, reference_weights; + std::vector query_groups, reference_groups; + std::vector> group_distances; uint32_t longest = 0; uint8_t A = 0; size_t W1 = 0; @@ -73,7 +77,13 @@ struct GapInput { GapInput prepare(const std::vector& queries, const std::vector& refs, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, - const std::vector& prior, uint32_t width, int threads) { + const std::vector& prior, uint32_t width, int threads, + const std::vector>& weights = {}, + const std::vector& query_groups = {}, + const std::vector& reference_groups = {}, + const std::vector>& group_distances = {}, + const std::vector>& query_weights = {}, + const std::vector>& reference_weights = {}) { if (go < 0 || ge < 0 || threads < 0) throw std::invalid_argument("gap costs and threads must be >= 0"); if (refs.size() > std::numeric_limits::max()) @@ -109,23 +119,95 @@ GapInput prepare(const std::vector& queries, const std::vector= weights.size() || weights[length].size() != length)) + throw std::invalid_argument("missing positional weights for an input frame"); + } + } + if (!query_weights.empty() || !reference_weights.empty()) { + if (!weights.empty()) throw std::invalid_argument("length and per-input weights cannot be combined"); + auto validate = [&](const auto& rows, const auto& axis) { + if (rows.size() != axis.size()) + throw std::invalid_argument("per-input weight rows must match their sequence axis"); + for (size_t i = 0; i < axis.size(); ++i) { + if (rows[i].size() != axis[i].size()) + throw std::invalid_argument("per-input weight length must match its sequence"); + for (auto value : rows[i]) { + if (value < 0) throw std::invalid_argument("positional weights must be nonnegative"); + largest_weight = std::max(largest_weight, value); + } + } + }; + largest_weight = 0; + validate(query_weights, in.q); + validate(reference_weights, in.r); + in.query_weights = query_weights; + in.reference_weights = reference_weights; + } + int32_t largest_offset = 0; + if (!query_groups.empty() || !reference_groups.empty() || !group_distances.empty()) { + if (query_groups.size() != queries.size() || reference_groups.size() != refs.size()) + throw std::invalid_argument("group IDs must match their query/reference axes"); + if (group_distances.empty() || group_distances.front().empty()) + throw std::invalid_argument("group distances must be a nonempty rectangular matrix"); + const auto columns = group_distances.front().size(); + for (const auto& row : group_distances) { + if (row.size() != columns) + throw std::invalid_argument("group distances must be rectangular"); + for (const auto value : row) { + if (value < 0) throw std::invalid_argument("group distances must be nonnegative"); + largest_offset = std::max(largest_offset, value); + } + } + for (const auto id : query_groups) + if (id < 0 || size_t(id) >= group_distances.size()) + throw std::invalid_argument("query group ID outside group distances"); + for (const auto id : reference_groups) + if (id < 0 || size_t(id) >= columns) + throw std::invalid_argument("reference group ID outside group distances"); + in.query_groups = query_groups; + in.reference_groups = reference_groups; + in.group_distances = group_distances; + } const auto largest_pen = *std::max_element(in.pen.begin(), in.pen.end()); if (*std::min_element(in.pen.begin(), in.pen.end()) < 0) throw std::invalid_argument("matrix penalties must be nonnegative"); - const uint64_t bound = uint64_t(in.longest) * uint32_t(largest_pen) + - uint32_t(largest_prior) + uint32_t(go) + + const uint64_t weighted_pen = uint64_t(uint32_t(largest_pen)) * uint32_t(largest_weight); + if (in.longest && weighted_pen > uint64_t(std::numeric_limits::max()) / in.longest) + throw std::invalid_argument("gapblock weighted score can overflow int32"); + const uint64_t bound = uint64_t(in.longest) * weighted_pen + + uint32_t(largest_prior) + uint32_t(largest_offset) + uint32_t(go) + uint64_t(in.longest ? in.longest - 1 : 0) * uint32_t(ge); if (bound > std::numeric_limits::max()) throw std::invalid_argument("gapblock score can overflow int32"); return in; } +inline int32_t group_offset(const GapInput& in, size_t i, size_t j) { + return in.group_distances.empty() ? 0 : in.group_distances[in.query_groups[i]][in.reference_groups[j]]; +} + int32_t score_cell(const GapInput& in, size_t i, size_t j, int32_t go, int32_t ge, const std::vector& prior, int32_t* scratch) { const uint32_t m = uint32_t(in.q[i].size()), n = uint32_t(in.r[j].size()); const uint32_t M = std::max(m, n), d = m > n ? m - n : n - m; const int32_t* prow = prior.empty() ? nullptr : prior.data() + (size_t(M) * in.W1 + d) * in.W1; - return cell(in.q[i].data(), m, in.r[j].data(), n, in.pen.data(), in.A, go, ge, prow, scratch); + return group_offset(in, i, j) + cell(in.q[i].data(), m, in.r[j].data(), n, in.pen.data(), in.A, go, ge, prow, scratch, in.weights.empty() || !M ? nullptr : in.weights[M].data(), + in.query_weights.empty() ? nullptr : in.query_weights[i].data(), + in.reference_weights.empty() ? nullptr : in.reference_weights[j].data()); } struct Ranked { @@ -179,12 +261,26 @@ std::vector> ranked(const GapInput& a, const GapInput* b, ui std::vector> counted(const GapInput& a, const GapInput* b, const std::vector>& thresholds, int32_t go, int32_t ge, - const std::vector& prior, int threads, bool exclude_exact) { + const std::vector& prior, int threads, bool exclude_exact, + std::vector>* masses = nullptr, bool sum_pair_distance = false) { if (thresholds.size() != a.q.size()) throw std::invalid_argument("thresholds must have one row per query"); if (b && (b->q.size() != a.q.size() || b->r.size() != a.r.size())) throw std::invalid_argument("paired query/reference axes must agree"); + if (sum_pair_distance && !b) + throw std::invalid_argument("sum distance requires two linked views"); std::vector> out(a.q.size()); + if (masses) { + masses->resize(a.q.size()); + if (std::any_of(a.pen.begin(), a.pen.end(), [](int32_t cost) { return cost < 0; }) || + (b && std::any_of(b->pen.begin(), b->pen.end(), [](int32_t cost) { return cost < 0; }))) + throw std::invalid_argument("linear mass requires nonnegative scores"); + // Every score admitted to a row is <= its largest positive threshold. + // Bound both score sums and threshold*count before workers start. + for (const auto& row : thresholds) for (auto threshold : row) + if (threshold > 0 && a.r.size() > std::numeric_limits::max() / uint64_t(threshold)) + throw std::overflow_error("gapblock linear mass exceeds uint64"); + } const auto longest = b ? std::max(a.longest, b->longest) : a.longest; parallel_for(a.q.size(), threads, [&] { return std::vector(size_t(longest) + 1); }, [&](size_t i, std::vector& scratch) { @@ -193,20 +289,38 @@ std::vector> counted(const GapInput& a, const GapInput* b, for (size_t t = 0; t < thresholds[i].size(); ++t) sorted.emplace_back(thresholds[i][t], t); std::sort(sorted.begin(), sorted.end()); std::vector bins(sorted.size()); + std::vector score_bins(masses ? sorted.size() : 0); out[i].resize(sorted.size()); + if (masses) (*masses)[i].resize(sorted.size()); if (sorted.empty()) return; for (size_t j = 0; j < a.r.size(); ++j) { if (exclude_exact && a.q[i] == a.r[j] && (!b || b->q[i] == b->r[j])) continue; - auto score = score_cell(a, i, j, go, ge, prior, scratch.data()); - if (b) score = std::max(score, score_cell(*b, i, j, go, ge, prior, scratch.data())); + const int64_t offset = int64_t(group_offset(a, i, j)) + + (sum_pair_distance ? int64_t(group_offset(*b, i, j)) : 0); + if (offset > sorted.back().first) continue; + int64_t score = score_cell(a, i, j, go, ge, prior, scratch.data()); + if (sum_pair_distance) { + if (score > sorted.back().first) continue; + score += int64_t(score_cell(*b, i, j, go, ge, prior, scratch.data())); + } else if (b) { + score = std::max(score, score_cell(*b, i, j, go, ge, prior, scratch.data())); + } const auto pos = std::lower_bound(sorted.begin(), sorted.end(), score, - [](const auto& threshold, int32_t value) { return threshold.first < value; }); - if (pos != sorted.end()) ++bins[size_t(pos - sorted.begin())]; + [](const auto& threshold, int64_t value) { return threshold.first < value; }); + if (pos != sorted.end()) { + const auto bin = size_t(pos - sorted.begin()); + ++bins[bin]; + if (masses) score_bins[bin] += uint64_t(score); + } } - uint64_t cumulative = 0; + uint64_t cumulative = 0, score_sum = 0; for (size_t t = 0; t < sorted.size(); ++t) { cumulative += bins[t]; out[i][sorted[t].second] = cumulative; + if (masses) { + score_sum += score_bins[t]; + (*masses)[i][sorted[t].second] = uint64_t(std::max(sorted[t].first, 0)) * cumulative - score_sum; + } } }); return out; @@ -215,60 +329,82 @@ std::vector> counted(const GapInput& a, const GapInput* b, } // namespace std::vector gapblock_matrix(const std::vector& queries, - const std::vector& refs, Alphabet alphabet, - const SubstitutionMatrix* matrix, int32_t gap_open, - int32_t gap_extend, const std::vector& prior, - uint32_t prior_width, int threads) { - if (gap_open < 0 || gap_extend < 0) - throw std::invalid_argument("gap_open and gap_extend must be >= 0"); - - const Codec codec(alphabet); - const uint8_t A = codec.size(); - if (matrix && matrix->size() != A) - throw std::invalid_argument("matrix size does not match the alphabet"); - - const size_t N = queries.size(), K = refs.size(); - std::vector out(N * K); - if (N == 0 || K == 0) return out; - - const std::vector pen = pen_table(matrix, A); - - std::vector> qc(N), rc(K); - uint32_t longest = 0; - for (size_t i = 0; i < N; ++i) { - qc[i] = encode(codec, queries[i], "queries", i); - longest = std::max(longest, uint32_t(qc[i].size())); - } - for (size_t k = 0; k < K; ++k) { - rc[k] = encode(codec, refs[k], "refs", k); - longest = std::max(longest, uint32_t(rc[k].size())); - } + const std::vector& refs, Alphabet alphabet, + const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto in = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, std::max(threads, 0), weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + if (!refs.empty() && queries.size() > std::numeric_limits::max() / refs.size()) + throw std::invalid_argument("gapblock matrix size overflows size_t"); + std::vector out(queries.size() * refs.size()); + parallel_for(queries.size(), threads, [&] { return std::vector(size_t(in.longest) + 1); }, + [&](size_t i, std::vector& scratch) { + for (size_t j = 0; j < refs.size(); ++j) + out[i * refs.size() + j] = score_cell(in, i, j, go, ge, prior, scratch.data()); + }); + return out; +} - const size_t W1 = size_t(prior_width) + 1; - if (!prior.empty()) { - if (prior.size() != W1 * W1 * W1) - throw std::invalid_argument("prior table must have (prior_width + 1)^3 entries"); - if (longest > prior_width) - throw std::invalid_argument("a sequence is longer than the prior table's width"); - } - const int32_t* P = prior.empty() ? nullptr : prior.data(); +std::vector> gapblock_ball_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + if (thresholds.size() != queries.size()) + throw std::invalid_argument("ball thresholds must have one value per query"); + const auto in = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + std::vector> out(queries.size()); + parallel_for(queries.size(), threads, [&] { return std::vector(size_t(in.longest) + 1); }, + [&](size_t i, std::vector& scratch) { + if (thresholds[i] < 0) return; + for (size_t j = 0; j < refs.size(); ++j) { + if (exclude_exact && in.q[i] == in.r[j]) continue; + if (group_offset(in, i, j) > thresholds[i]) continue; + const auto score = score_cell(in, i, j, go, ge, prior, scratch.data()); + if (score <= thresholds[i]) out[i].push_back(Hit{uint32_t(j), score, 0, 0, 0}); + } + }); + return out; +} - // Every symbol was validated above and the kernel is arithmetic, so no worker throws -- - // parallel_for's exception plumbing simply never fires here. `suf` is the per-worker - // suffix-score scratch, allocated once per thread rather than once per row. - parallel_for(N, threads, [&] { return std::vector(size_t(longest) + 1); }, - [&](size_t i, std::vector& suf) { - const uint8_t* q = qc[i].data(); - const uint32_t m = uint32_t(qc[i].size()); - int32_t* row = out.data() + i * K; - for (size_t k = 0; k < K; ++k) { - const uint32_t n = uint32_t(rc[k].size()); - const uint32_t M = std::max(m, n), d = (m > n ? m - n : n - m); - const int32_t* prow = P ? P + (size_t(M) * W1 + d) * W1 : nullptr; - row[k] = cell(q, m, rc[k].data(), n, pen.data(), A, gap_open, gap_extend, - prow, suf.data()); - } - }); +std::tuple, std::vector, std::vector> +gapblock_group_audit_batch(const std::vector& queries, + const std::vector& refs, const std::vector& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, + const std::vector>& query_weights, + const std::vector>& reference_weights) { + if (thresholds.size() != queries.size()) + throw std::invalid_argument("group audit thresholds must have one value per query"); + const auto in = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads, + weights, query_groups, reference_groups, group_distances, + query_weights, reference_weights); + std::tuple, std::vector, std::vector> out{ + std::vector(queries.size()), std::vector(queries.size()), + std::vector(queries.size(), -1)}; + auto& [before, after, nearest] = out; + parallel_for(queries.size(), threads, + [&] { return std::vector(size_t(in.longest) + 1); }, + [&](size_t i, std::vector& scratch) { + if (thresholds[i] < 0) return; + int32_t minimum = std::numeric_limits::max(); + for (size_t j = 0; j < refs.size(); ++j) { + if (exclude_exact && in.q[i] == in.r[j]) continue; + const auto score = score_cell(in, i, j, go, ge, prior, scratch.data()); + if (score - group_offset(in, i, j) > thresholds[i]) continue; + ++before[i]; + after[i] += score <= thresholds[i]; + minimum = std::min(minimum, score); + } + if (before[i]) nearest[i] = minimum; + }); return out; } @@ -303,11 +439,27 @@ std::vector>> gapblock_paired std::vector> gapblock_count_batch(const std::vector& queries, const std::vector& refs, const std::vector>& thresholds, Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, - const std::vector& prior, uint32_t width, int threads, bool exclude_exact) { - const auto a = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads); + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto a = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); return counted(a, nullptr, thresholds, go, ge, prior, threads, exclude_exact); } +std::pair>, std::vector>> gapblock_count_mass_batch(const std::vector& queries, + const std::vector& refs, const std::vector>& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights, + const std::vector& query_groups, const std::vector& reference_groups, + const std::vector>& group_distances, const std::vector>& query_weights, const std::vector>& reference_weights) { + const auto a = prepare(queries, refs, alphabet, matrix, go, ge, prior, width, threads, weights, query_groups, reference_groups, group_distances, query_weights, reference_weights); + std::vector> masses; + auto counts = counted(a, nullptr, thresholds, go, ge, prior, threads, exclude_exact, &masses); + return {std::move(counts), std::move(masses)}; +} + std::vector> gapblock_paired_count_batch( const std::vector& qa, const std::vector& qb, const std::vector& ra, const std::vector& rb, @@ -319,4 +471,29 @@ std::vector> gapblock_paired_count_batch( return counted(a, &b, thresholds, go, ge, prior, threads, exclude_exact); } +std::pair>, std::vector>> gapblock_paired_sum_count_batch( + const std::vector& qa, const std::vector& qb, + const std::vector& ra, const std::vector& rb, + const std::vector>& thresholds, + Alphabet alphabet, const SubstitutionMatrix* matrix, int32_t go, int32_t ge, + const std::vector& prior, uint32_t width, int threads, bool exclude_exact, + const std::vector>& weights_alpha, + const std::vector>& weights_beta, + const std::vector& query_groups_alpha, const std::vector& query_groups_beta, + const std::vector& reference_groups_alpha, const std::vector& reference_groups_beta, + const std::vector>& group_distances_alpha, + const std::vector>& group_distances_beta, + const std::vector>& query_weights_alpha, + const std::vector>& query_weights_beta, + const std::vector>& reference_weights_alpha, + const std::vector>& reference_weights_beta) { + const auto a = prepare(qa, ra, alphabet, matrix, go, ge, prior, width, threads, + weights_alpha, query_groups_alpha, reference_groups_alpha, group_distances_alpha, query_weights_alpha, reference_weights_alpha); + const auto b = prepare(qb, rb, alphabet, matrix, go, ge, prior, width, threads, + weights_beta, query_groups_beta, reference_groups_beta, group_distances_beta, query_weights_beta, reference_weights_beta); + std::vector> masses; + auto counts = counted(a, &b, thresholds, go, ge, prior, threads, exclude_exact, &masses, true); + return {std::move(counts), std::move(masses)}; +} + } // namespace seqtree diff --git a/src/index.cpp b/src/index.cpp index 8e87f9d..0d043f0 100644 --- a/src/index.cpp +++ b/src/index.cpp @@ -2,6 +2,7 @@ #include "seqtree/parallel.hpp" #include "atomic_write.hpp" #include "trie.hpp" +#include "count_unit.hpp" #include #include @@ -187,6 +188,30 @@ std::vector> Index::search_batch(const std::vector return results; } +std::vector> Index::edit_histogram_batch( + const std::vector& queries, const SearchParams& p, + int threads, bool exclude_exact) const { + if (!p.max_total_edits) + throw std::invalid_argument("edit histogram requires explicit positive max_total_edits"); + if (p.engine == Engine::SeqTrie || p.mode != Mode::AllHits || p.max_hits) + throw std::invalid_argument("edit histogram requires uncapped seqtm/auto all-hit search"); + std::vector> out(queries.size(), + std::vector(size_t(p.max_total_edits) + 1)); + struct Local { Searcher s; std::vector hits; }; + parallel_for(queries.size(), threads, [&] { return Local{Searcher(*this), {}}; }, + [&](size_t i, Local& worker) { + if (count_unit::fill_histogram(*this,queries[i],p,exclude_exact,out[i])) return; + worker.s.search_into(queries[i], p, worker.hits); + for (const auto& h : worker.hits) { + const uint32_t edits = uint32_t(h.n_subs) + h.n_ins + h.n_dels; + if (edits > p.max_total_edits) + throw std::logic_error("search hit exceeds max_total_edits"); + if (!exclude_exact || edits) ++out[i][edits]; + } + }); + return out; +} + std::vector Index::collisions_batch(const std::vector& queries, const SearchParams& p, int threads) const { const size_t n = queries.size(); diff --git a/src/text_index.cpp b/src/text_index.cpp index 41221a8..dafe2b3 100644 --- a/src/text_index.cpp +++ b/src/text_index.cpp @@ -42,6 +42,7 @@ #include "seqtree/text_index.hpp" #include "seqtree/parallel.hpp" +#include #include "atomic_write.hpp" @@ -65,6 +66,24 @@ #endif namespace seqtree { + +std::vector> TextResult::record_distances() const { + if (std::any_of(truncated.begin(), truncated.end(), [](uint8_t flag) { return flag != 0; })) + throw std::invalid_argument("truncated text hits cannot provide complete record counts"); + std::vector> out(num_queries()); + std::unordered_map best; + for (size_t i = 0; i < num_queries(); ++i) { + best.clear(); + for (uint32_t h = query_begin[i]; h < query_begin[i + 1]; ++h) { + auto [it, inserted] = best.emplace(ref_id[h], n_subs[h]); + if (!inserted) it->second = std::min(it->second, n_subs[h]); + } + out[i].reserve(best.size()); + for (const auto& item : best) out[i].push_back(item.second); + std::sort(out[i].begin(), out[i].end()); + } + return out; +} namespace { constexpr char kMagic[4] = {'S', 'Q', 'T', 'X'}; diff --git a/tests/cpp/test_parallel.cpp b/tests/cpp/test_parallel.cpp index aecf709..1b3609e 100644 --- a/tests/cpp/test_parallel.cpp +++ b/tests/cpp/test_parallel.cpp @@ -14,9 +14,9 @@ TEST_CASE("serial batches run inline and create one scratch") { }); CHECK(factories == 1); CHECK(seen == 17); - seqtree::parallel_for(0, 8, [&] { ++factories; return 0; }, [](size_t, int&) {}); + seqtree::parallel_for(0, 1, [&] { ++factories; return 0; }, [](size_t, int&) {}); CHECK(factories == 1); - seqtree::parallel_for(1, 8, [] { return 0; }, [&](size_t, int&) { + seqtree::parallel_for(1, 1, [] { return 0; }, [&](size_t, int&) { CHECK(std::this_thread::get_id() == caller); }); } @@ -27,23 +27,24 @@ TEST_CASE("parallel batches visit each item once and join every scratch") { explicit Local(std::atomic& a) : alive(a) { ++alive; } ~Local() { --alive; } }; - std::atomic alive{0}, factories{0}; + std::atomic alive{0}, factories{0}, caller_factories{0}; std::vector> visits(1031); const auto caller = std::this_thread::get_id(); - seqtree::parallel_for(visits.size(), 4, [&] { + seqtree::parallel_for(visits.size(), 2, [&] { ++factories; + if (std::this_thread::get_id() == caller) ++caller_factories; return std::make_unique(alive); }, [&](size_t i, const auto&) { - CHECK(std::this_thread::get_id() != caller); ++visits[i]; }); - CHECK(factories == 4); + CHECK(factories == 2); + CHECK(caller_factories == 1); CHECK(alive == 0); for (const auto& n : visits) CHECK(n == 1); } TEST_CASE("factory and body exceptions propagate after workers join") { - for (int threads : {1, 4}) { + for (int threads : {1, 2}) { CHECK_THROWS_WITH_AS(seqtree::parallel_for(100, threads, []() -> int { throw std::runtime_error("scratch failed"); }, [](size_t, int&) {}), "scratch failed", std::runtime_error); @@ -58,3 +59,11 @@ TEST_CASE("factory and body exceptions propagate after workers join") { CHECK(alive == 0); } } + +TEST_CASE("parallel batches reject an excessive live worker budget") { + const unsigned available = seqtree::available_cpus(); + REQUIRE(available >= 1); + CHECK_THROWS_AS(seqtree::parallel_for(size_t(available) + 1, int(available) + 1, + [] { return 0; }, [](size_t, int&) {}), + std::invalid_argument); +} diff --git a/tests/cpp/test_text_index.cpp b/tests/cpp/test_text_index.cpp index 47a43f6..026a9aa 100644 --- a/tests/cpp/test_text_index.cpp +++ b/tests/cpp/test_text_index.cpp @@ -249,7 +249,7 @@ TEST_CASE("text index: hit order is deterministic across thread counts") { TextQueryOpts opts; opts.max_subs = 2; auto one = ix->search_batch(queries, opts, 1); - for (int t : {2, 4, 8}) { + for (int t : {1, 2, 0}) { auto many = ix->search_batch(queries, opts, t); CHECK(many.query_begin == one.query_begin); CHECK(many.ref_id == one.ref_id); @@ -459,7 +459,7 @@ TEST_CASE("many threads may search one index at once") { std::vector callers; std::atomic mismatches{0}; - for (int t = 0; t < 8; ++t) + for (int t = 0; t < 2; ++t) callers.emplace_back([&] { for (int rep = 0; rep < 4; ++rep) { auto got = ix->search_batch(queries, opts, 1); diff --git a/tests/python/test_batch_threads.py b/tests/python/test_batch_threads.py index aaabf05..77aaf6e 100644 --- a/tests/python/test_batch_threads.py +++ b/tests/python/test_batch_threads.py @@ -1,4 +1,9 @@ import random +import os +import subprocess +import sys + +import pytest import seqtree @@ -20,7 +25,7 @@ def test_batch_equals_serial(): p = seqtree.SearchParams(max_subs=2, engine="seqtm") serial = [idx.search(q, p) for q in queries] - batch = idx.search_batch(queries, p, threads=8) + batch = idx.search_batch(queries, p, threads=2) assert len(serial) == len(batch) for s, b in zip(serial, batch): assert _key(s) == _key(b) @@ -33,7 +38,7 @@ def test_threads_determinism(): p = seqtree.SearchParams(max_subs=1, max_ins=1, max_dels=1, engine="seqtm") one = idx.search_batch(queries, p, threads=1) - many = idx.search_batch(queries, p, threads=16) + many = idx.search_batch(queries, p, threads=2) assert [_key(r) for r in one] == [_key(r) for r in many] @@ -48,5 +53,27 @@ def test_pairwise_symmetry_of_membership(): # cross-check against explicit index-on-b idx_b = seqtree.Index.build(b, alphabet="aa") - expected = idx_b.search_batch(a, p, threads=4) + expected = idx_b.search_batch(a, p, threads=2) assert [_key(r) for r in res] == [_key(r) for r in expected] + + +@pytest.mark.skipif(not hasattr(os, "sched_setaffinity"), reason="Linux affinity contract") +def test_native_batch_respects_restricted_allocation(): + subprocess.run([sys.executable, "-c", ''' +import os +os.sched_setaffinity(0, {min(os.sched_getaffinity(0))}) +import seqtree +index = seqtree.Index.build(["AAAA", "AAAC"]) +params = seqtree.SearchParams(max_subs=1, max_total_edits=1) +queries = ["AAAA"] * 8 +assert index.edit_histogram_batch(queries, params, threads=0) == [[1, 1]] * 8 +one = index.search_batch(queries, params, threads=1) +all_available = index.search_batch(queries, params, threads=0) +assert [[(h.ref_id, h.score) for h in row] for row in one] == [[(h.ref_id, h.score) for h in row] for row in all_available] +try: + index.search_batch(queries, params, threads=2) +except ValueError: + pass +else: + raise AssertionError("explicit oversubscription was accepted") +'''], check=True, timeout=30) diff --git a/tests/python/test_distance.py b/tests/python/test_distance.py index bd618c9..e2f1414 100644 --- a/tests/python/test_distance.py +++ b/tests/python/test_distance.py @@ -153,7 +153,7 @@ def test_matrix_thread_counts_agree(): a = [rand_seq(rng, "ACGT", 0, 15) for _ in range(30)] b = [rand_seq(rng, "ACGT", 0, 15) for _ in range(20)] one = np.asarray(levenshtein_matrix(a, b, threads=1)) - many = np.asarray(levenshtein_matrix(a, b, threads=8)) + many = np.asarray(levenshtein_matrix(a, b, threads=2)) assert (one == many).all() diff --git a/tests/python/test_doc_coverage.py b/tests/python/test_doc_coverage.py index da9694e..a402464 100644 --- a/tests/python/test_doc_coverage.py +++ b/tests/python/test_doc_coverage.py @@ -15,7 +15,7 @@ DOCS = pathlib.Path(__file__).resolve().parents[2] / "docs" SUBMODULES = ["distance", "gapblock", "pairwise", "seeds", "evalue", "layout", "pmhc", "control", - "pmhc_evalue"] + "pmhc_evalue", "text_evalue"] #: Modules with no public functions/classes of their own need no page of their own. NO_PAGE_REQUIRED = {"control", "pmhc_evalue"} diff --git a/tests/python/test_edit_histogram.py b/tests/python/test_edit_histogram.py new file mode 100644 index 0000000..f07cc29 --- /dev/null +++ b/tests/python/test_edit_histogram.py @@ -0,0 +1,31 @@ +"""Native histogram must preserve the existing search's edit counts and multiplicity.""" +import pytest +from seqtree import Index, SearchParams + + +def test_edit_histogram_matches_search(): + index = Index.build(['ACDE', 'ACDE', 'ACDF', 'ACD', 'ACGDE', 'GGGG']) + queries = ['ACDE', 'ACD', 'GGGG', ''] + params = SearchParams(max_subs=2, max_ins=2, max_dels=2, max_total_edits=2) + for exclude in (False, True): + expected = [] + for row in index.search_batch(queries, params, threads=1): + bins = [0, 0, 0] + for hit in row: + total = hit.n_subs + hit.n_ins + hit.n_dels + if not (exclude and total == 0): + bins[total] += 1 + expected.append(bins) + assert index.edit_histogram_batch(queries, params, threads=1, exclude_exact=exclude) == expected + for threads in (2, 0, -1): + assert index.edit_histogram_batch(queries, params, threads=threads, exclude_exact=exclude) == expected + assert Index.build([]).edit_histogram_batch(['ACDE'], params, threads=1) == [[0, 0, 0]] + assert index.edit_histogram_batch([], params, threads=1) == [] + + +@pytest.mark.parametrize('params', [SearchParams(max_subs=2), + SearchParams(max_total_edits=2, engine='seqtrie'), + SearchParams(max_subs=2, max_total_edits=2, mode='top')]) +def test_edit_histogram_rejects_incomplete_counts(params): + with pytest.raises(ValueError): + Index.build(['ACDE']).edit_histogram_batch(['ACDE'], params, threads=1) diff --git a/tests/python/test_engine_guards.py b/tests/python/test_engine_guards.py index 938816b..1eb041c 100644 --- a/tests/python/test_engine_guards.py +++ b/tests/python/test_engine_guards.py @@ -1,6 +1,6 @@ """Regression tests for the engine-routing and scoring bugs fixed on feature/engine-bugfix. -Each test pins one previously-silent wrong answer. See ROADMAP.md and the commit message. +Each test pins one previously-silent wrong answer. See docs/roadmap.rst and the commit message. """ import pytest diff --git a/tests/python/test_evalue.py b/tests/python/test_evalue.py index f160182..6e7389e 100644 --- a/tests/python/test_evalue.py +++ b/tests/python/test_evalue.py @@ -3,7 +3,7 @@ import pytest import seqtree -from seqtree.evalue import _poisson_sf, evalues, thetas_from_scores, threshold_for_evalue +from seqtree.evalue import _poisson_sf, evalue_result, evalues, thetas_from_scores, threshold_for_evalue AA = "ACDEFGHIKLMNPQRSTVWY" @@ -32,10 +32,35 @@ def _hamming_params(theta, qlen): def test_poisson_sf_basics(): assert _poisson_sf(0, 5.0) == 1.0 assert _poisson_sf(1, 0.0) == 0.0 + assert _poisson_sf(4, float("inf")) == 1.0 # P(Poisson(1) >= 1) = 1 - e^-1 assert abs(_poisson_sf(1, 1.0) - (1 - 2.718281828 ** -1)) < 1e-6 +@pytest.mark.parametrize("k,lam,expected", [ + (1, 1e-20, 1e-20), + (4, 1e-5, 4.1666333334722077e-22), + (2000, 1000.0, 3.058192080169065e-170), + (1000, 1000.0, 0.5042052441802155), + (1001, 1000.0, 0.491590632831494), +]) +def test_poisson_tail_avoids_cancellation_and_underflow(k, lam, expected): + # Independent scipy.stats.poisson.sf(k-1, lam) reference values; no runtime dependency. + assert _poisson_sf(k, lam) == pytest.approx(expected, rel=1e-10, abs=0) + + +def test_poisson_tail_survives_subnormal_pmf(): + # Independent high-precision Decimal references, rounded to binary64. + assert _poisson_sf(14086, 10000.0) == 5e-324 + assert _poisson_sf(112398, 100000.0) == 2e-323 + + +def test_public_evalue_keeps_tiny_nonzero_probabilities(): + result = evalue_result(4, 1, 1, 100000) + assert result["p_enrichment"] == pytest.approx(4.1666333334722077e-22, rel=1e-10, abs=0) + assert evalue_result(1, 1, 1, 10**20)["p_any"] == pytest.approx(1e-20, rel=1e-10, abs=0) + + def test_planted_cluster_is_significant_background_is_not(): rng = random.Random(0) control = seqtree.Index.build(_rand(20000, rng), alphabet="aa") diff --git a/tests/python/test_gapblock.py b/tests/python/test_gapblock.py index b7884b1..a7bb6b3 100644 --- a/tests/python/test_gapblock.py +++ b/tests/python/test_gapblock.py @@ -515,7 +515,7 @@ def test_score_matrix_is_thread_invariant(): r = [rand_seq(rng, rng.randint(6, 18)) for _ in range(40)] prior = central_prior(21) ref = [score_matrix(q, r, BLOSUM, 28, 1, prior, threads=1).row(i) for i in range(len(q))] - for threads in (2, 4, 0): + for threads in (1, 2, 0): got = score_matrix(q, r, BLOSUM, 28, 1, prior, threads=threads) assert [got.row(i) for i in range(len(q))] == ref diff --git a/tests/python/test_gapblock_group_audit.py b/tests/python/test_gapblock_group_audit.py new file mode 100644 index 0000000..0d011fe --- /dev/null +++ b/tests/python/test_gapblock_group_audit.py @@ -0,0 +1,21 @@ +"""Streaming group audit preserves exact component/full balls without hit materialization.""" +import numpy as np +import pytest +from seqtree import gapblock + +def test_group_audit_matches_dense_components_and_threads(): + q=['AA','CA','aa',''];r=['AA','CC','CAA','AA',''] + opts=dict(gap_open=2,gap_extend=1,position_weights_by_length={0:[],2:[0,3],3:[2,1,0]}) + group=dict(query_group_ids=[0,1,0,1],reference_group_ids=[0,1,0,1,1],group_distances=[[0,5],[3,0]]) + base=np.asarray(gapblock.score_matrix(q,r,**opts));full=np.asarray(gapblock.score_matrix(q,r,**opts,**group)) + limits=[0,3,10,-1] + for threads in [1,2]: + for excluded in [False,True]: + before,after,nearest=gapblock.group_audit_batch(q,r,limits,threads=threads,exclude_exact=excluded,**opts,**group) + for i,t in enumerate(limits): + valid=np.array([not excluded or q[i].upper()!=s.upper() for s in r]);selected=valid&(base[i]<=t) + assert before[i]==int(selected.sum()) + assert after[i]==int((valid&(full[i]<=t)).sum()) + assert nearest[i]==(int(full[i,selected].min()) if selected.any() else -1) + assert gapblock.group_audit_batch(['AA'],[],[0])==([0],[0],[-1]) + with pytest.raises(ValueError):gapblock.group_audit_batch(q,r,[]) diff --git a/tests/python/test_gapblock_group_offsets.py b/tests/python/test_gapblock_group_offsets.py new file mode 100644 index 0000000..df1a031 --- /dev/null +++ b/tests/python/test_gapblock_group_offsets.py @@ -0,0 +1,91 @@ +"""Generic additive group penalties share the existing gapblock cell and reducers.""" +import numpy as np +import pytest +from seqtree import gapblock + + +def test_group_offsets_matrix_reducers_and_puncture(): + q=['AC','ACC','','AC']; r=['AC','CA','ACC','CA',''] + qg=[0,1,1,0]; rg=[1,0,1,0,0]; distances=[[7,2],[0,11]] + opts=dict(gap_open=3,gap_extend=2,position_weights_by_length={0:[],2:[0,2],3:[3,1,2]}) + base=np.asarray(gapblock.score_matrix(q,r,threads=1,**opts)) + expected=base+np.array([[distances[a][b] for b in rg] for a in qg]) + group=dict(query_group_ids=qg,reference_group_ids=rg,group_distances=distances) + thresholds=[[-1,0,2,5,20,5],[],[0,3,20],[20]] + for t in (1,2): + np.testing.assert_array_equal(np.asarray(gapblock.score_matrix(q,r,threads=t,**opts,**group)),expected) + for excluded in (False,True): + counts=gapblock.count_batch(q,r,thresholds,threads=t,exclude_exact=excluded,**opts,**group) + balls=gapblock.ball_batch(q,r,[5]*len(q),threads=t,exclude_exact=excluded,**opts,**group) + for i in range(len(q)): + eligible=lambda j:not(excluded and q[i]==r[j]) + assert counts[i]==[sum(eligible(j) and expected[i,j]<=c for j in range(len(r))) for c in thresholds[i]] + assert [(h.ref_id,h.score) for h in balls[i]]==[(j,int(expected[i,j])) for j in range(len(r)) if eligible(j) and expected[i,j]<=5] + # A full-junction exact match remains punctured even with a positive group penalty. + assert gapblock.count_batch(['AC'],['AC'],[[100]],exclude_exact=True,query_group_ids=[0],reference_group_ids=[1],group_distances=distances)==[[0]] + # Offset alone rejects the entire group even when its junction distance is zero. + assert gapblock.ball_batch(['AC'],['AC'],[1],query_group_ids=[0],reference_group_ids=[1],group_distances=distances)==[[]] + + +@pytest.mark.parametrize('api',['score_matrix','count_batch','ball_batch']) +@pytest.mark.parametrize('groups',[ + dict(query_group_ids=[0]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[]), + dict(query_group_ids=[],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[-1],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[1],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[0],reference_group_ids=[1],group_distances=[[0]]), + dict(query_group_ids=[0.5],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[False],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[0],[]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[-1]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[1.5]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[2**31-1]]), +]) +def test_group_validation(api,groups): + args=(['AC'],['CA']) + if api=='count_batch':args+=([[10]],) + if api=='ball_batch':args+=([10],) + with pytest.raises((ValueError,TypeError,OverflowError)): + getattr(gapblock,api)(*args,threads=1,**groups) + + +def test_defaults_and_empty_axes(): + q=['AC',''];r=['AC','CA'] + base=np.asarray(gapblock.score_matrix(q,r,threads=1)) + zero=dict(query_group_ids=[0,0],reference_group_ids=[0,0],group_distances=[[0]]) + np.testing.assert_array_equal(base,np.asarray(gapblock.score_matrix(q,r,threads=-1,**zero))) + assert gapblock.count_batch([],r,[],query_group_ids=[],reference_group_ids=[0,0],group_distances=[[5]])==[] + assert gapblock.count_batch(q,[],[[0],[5]],query_group_ids=[0,0],reference_group_ids=[],group_distances=[[5]])==[[0],[0]] + + +@pytest.mark.parametrize('api',['gapblock_matrix','gapblock_count_batch','gapblock_ball_batch']) +@pytest.mark.parametrize('groups',[ + dict(query_group_ids=[],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[-1],reference_group_ids=[0],group_distances=[[0]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[0],[]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[-1]]), + dict(query_group_ids=[0],reference_group_ids=[0],group_distances=[[2**31-1]]), +]) +def test_native_validation_before_workers(api,groups): + from seqtree import _core + args=(['AC'],['CA']) + if api=='gapblock_count_batch':args+=([[10]],) + if api=='gapblock_ball_batch':args+=([10],) + with pytest.raises(ValueError): + getattr(_core,api)(*args,threads=2,**groups) + + +def test_rectangular_group_table(): + opts=dict(query_group_ids=[0],reference_group_ids=[2,0],group_distances=[[7,4,1]]) + observed=np.asarray(gapblock.score_matrix(['AC'],['AC','AC'],threads=1,**opts)) + np.testing.assert_array_equal(observed,[[1,7]]) + + +def test_int32_boundary_group_only_score(): + maximum=2**31-1 + opts=dict(gap_open=0,gap_extend=0,query_group_ids=[0],reference_group_ids=[0],group_distances=[[maximum]],threads=2) + assert gapblock.score_matrix([''],[''],**opts)[0,0]==maximum + assert gapblock.count_batch([''],[''],[[maximum-1,maximum]],**opts)==[[0,1]] + assert [(h.ref_id,h.score) for h in gapblock.ball_batch([''],[''],[maximum],**opts)[0]]==[(0,maximum)] + assert gapblock.count_batch([''],[''],[[maximum]],exclude_exact=True,**opts)==[[0]] diff --git a/tests/python/test_gapblock_linear_mass.py b/tests/python/test_gapblock_linear_mass.py new file mode 100644 index 0000000..856b76e --- /dev/null +++ b/tests/python/test_gapblock_linear_mass.py @@ -0,0 +1,67 @@ +"""Integer linear taper shares exact counts, scoring and full-identity puncture.""" +import numpy as np +import pytest +from seqtree import gapblock, SubstitutionMatrix + + +def scalar(scores, thresholds, queries, refs, exclude_exact): + counts=[];masses=[] + for i,row in enumerate(thresholds): + eligible=[int(scores[i,j]) for j in range(len(refs)) + if not(exclude_exact and queries[i]==refs[j])] + counts.append([sum(d<=t for d in eligible) for t in row]) + masses.append([sum(max(0,t-d) for d in eligible) for t in row]) + return counts,masses + + +@pytest.mark.parametrize('exclude_exact',[False,True]) +def test_linear_mass_scalar_oracle_and_threads(exclude_exact): + q=['AC','ACC','','AC'];r=['AC','CA','ACC','CA',''] + thresholds=[[-1,0,2,5,20,5],[],[0,3,20],[20,0]] + options=dict(gap_open=3,gap_extend=2,position_weights_by_length={0:[],2:[0,2],3:[3,1,2]}, + query_group_ids=[0,1,1,0],reference_group_ids=[1,0,1,0,0],group_distances=[[7,2],[0,11]]) + scores=np.asarray(gapblock.score_matrix(q,r,threads=1,**options)) + expected=scalar(scores,thresholds,q,r,exclude_exact) + for threads in (1,2): + actual=gapblock.count_batch(q,r,thresholds,threads=threads,exclude_exact=exclude_exact, + linear_mass=True,**options) + assert actual==expected + assert actual[0]==gapblock.count_batch(q,r,thresholds,threads=threads, + exclude_exact=exclude_exact,**options) + + +def test_linear_mass_boundaries_and_uint64(): + assert gapblock.count_batch(['AC'],['AC','CA'],[[0,2,1,-1,2]],linear_mass=True)==( + [[1,2,1,0,2]],[[0,2,1,0,2]]) + maximum=2**31-1 + assert gapblock.count_batch(['AC'],['AC']*4,[[maximum]],linear_mass=True)==([[4]],[[4*maximum]]) + assert gapblock.count_batch(['AC'],['AC'],[[maximum]],exclude_exact=True,linear_mass=True)==([[0]],[[0]]) + assert gapblock.count_batch(['AC'],['AC'],[[5]],linear_mass=True, + query_group_ids=[0],reference_group_ids=[0],group_distances=[[6]])==([[0]],[[0]]) + + +def test_linear_mass_empty_axes_and_prior(): + assert gapblock.count_batch([],['AC'],[],linear_mass=True)==([],[]) + assert gapblock.count_batch(['AC'],[],[[-1,0,2]],linear_mass=True)==([[0,0,0]],[[0,0,0]]) + q=['AC','ACC'];r=['ACC','AC'] + opts=dict(gap_prior=gapblock.positions_prior((1,)),matrix=SubstitutionMatrix.blosum62(),gap_open=28,gap_extend=14) + scores=np.asarray(gapblock.score_matrix(q,r,**opts)) + cut=[[0,28,50],[50,28,0]] + assert gapblock.count_batch(q,r,cut,linear_mass=True,**opts)==scalar(scores,cut,q,r,False) + + +@pytest.mark.parametrize('value',[1,None,'yes']) +def test_linear_mass_requires_explicit_boolean(value): + with pytest.raises((ValueError,TypeError)): + gapblock.count_batch(['AC'],['CA'],[[2]],linear_mass=value) + + +@pytest.mark.parametrize('options',[ + {'gap_open':-1}, {'gap_extend':-1}, + {'query_group_ids':[0],'reference_group_ids':[0],'group_distances':[[-1]]}, + {'position_weights_by_length':[[],[],[-1,1]]}, +]) +def test_linear_mass_native_validation(options): + from seqtree import _core + with pytest.raises((ValueError,TypeError,OverflowError)): + _core.gapblock_count_mass_batch(['AC'],['CA'],[[2]],**options) diff --git a/tests/python/test_gapblock_paired_sum.py b/tests/python/test_gapblock_paired_sum.py new file mode 100644 index 0000000..091bbec --- /dev/null +++ b/tests/python/test_gapblock_paired_sum.py @@ -0,0 +1,92 @@ +"""Linked sum reduction agrees with dense lane sums without int32 truncation.""" +import numpy as np +import pytest +from seqtree import gapblock + + +@pytest.mark.parametrize('exclude_exact',[False,True]) +def test_paired_sum_dense_oracle_and_threads(exclude_exact): + qa=['AC','ACC','','AC'];qb=['CA','AC','AC','CA'] + ra=['AC','AC','ACC','','CA'];rb=['CA','AC','AC','AC','CA'] + cut=[[-1,0,3,20,5,20],[],[0,3,20],[20,0]] + common=dict(gap_open=3,gap_extend=2,gap_prior=gapblock.positions_prior((1,))) + alpha=dict(position_weights_by_length={0:[],2:[0,2],3:[3,1,2]}, + query_group_ids=[0,1,1,0],reference_group_ids=[1,0,1,0,0],group_distances=[[7,2],[0,11]]) + beta=dict(position_weights_by_length={2:[1,3]}, + query_group_ids=[1,0,1,1],reference_group_ids=[1,0,0,1,0],group_distances=[[0,6],[4,0]]) + A=np.asarray(gapblock.score_matrix(qa,ra,**common,**alpha)).astype(np.int64) + B=np.asarray(gapblock.score_matrix(qb,rb,**common,**beta)).astype(np.int64) + expected_counts=[];expected_mass=[] + for i,row in enumerate(cut): + scores=[int(A[i,j]+B[i,j]) for j in range(len(ra)) + if not(exclude_exact and qa[i]==ra[j] and qb[i]==rb[j])] + expected_counts.append([sum(d<=t for d in scores) for t in row]) + expected_mass.append([sum(max(0,t-d) for d in scores) for t in row]) + opts={**common,**{key+'_alpha':value for key,value in alpha.items()}, + **{key+'_beta':value for key,value in beta.items()}} + for threads in (1,2): + assert gapblock.paired_sum_count_batch(qa,qb,ra,rb,cut,threads=threads, + exclude_exact=exclude_exact,**opts)==(expected_counts,expected_mass) + + +def test_paired_sum_full_identity_puncture_and_ties(): + # Alpha alone matches every reference; only the first linked full pair is exact. + assert gapblock.paired_sum_count_batch(['AC'],['CA'],['AC','AC'],['CA','AC'],[[0,2,3]], + exclude_exact=True)==([[0,1,1]],[[0,0,1]]) + # Existing max contract remains different, unchanged. + assert gapblock.paired_count_batch(['AC'],['AC'],['CA'],['CA'],[[2]])==[[1]] + assert gapblock.paired_sum_count_batch(['AC'],['AC'],['CA'],['CA'],[[2]])==([[0]],[[0]]) + + +def test_paired_sum_int32_overflow_and_empty_axes(): + maximum=2**31-1 + groups=dict(query_group_ids_alpha=[0],reference_group_ids_alpha=[0],group_distances_alpha=[[maximum]], + query_group_ids_beta=[0],reference_group_ids_beta=[0],group_distances_beta=[[maximum]]) + assert gapblock.paired_sum_count_batch([''],[''],[''],[''],[[-1,0,maximum]], + gap_open=0,gap_extend=0,**groups)==([[0,0,0]],[[0,0,0]]) + groups['group_distances_beta']=[[0]] + assert gapblock.paired_sum_count_batch([''],[''],[''],[''],[[maximum]], + gap_open=0,gap_extend=0,**groups)==([[1]],[[0]]) + assert gapblock.paired_sum_count_batch([],[],['AC'],['CA'],[])==([],[]) + assert gapblock.paired_sum_count_batch(['AC'],['CA'],[],[],[[0,3]])==([[0,0]],[[0,0]]) + + +@pytest.mark.parametrize('args,opts',[ + ((['AC'],[],['AC'],['CA'],[[2]]),{}), + ((['AC'],['CA'],['AC'],[],[[2]]),{}), + ((['AC'],['CA'],['AC'],['CA'],[]),{}), + ((['AC'],['CA'],['AC'],['CA'],[[2]]),{'position_weights_by_length_beta':{2:[-1,1]}}), + ((['AC'],['CA'],['AC'],['CA'],[[2]]),{'query_group_ids_beta':[0]}), +]) +def test_paired_sum_validation(args,opts): + with pytest.raises((ValueError,TypeError)): + gapblock.paired_sum_count_batch(*args,**opts) + + +@pytest.mark.parametrize('args,opts',[ + ((['AC'],[],['AC'],['CA'],[[2]]),{}), + ((['AC'],['CA'],['AC'],['CA'],[[2]]),{'gap_open':-1}), + ((['AC'],['CA'],['AC'],['CA'],[[2]]),{'query_group_ids_beta':[0]}), + ((['AC'],['CA'],['AC'],['CA'],[[2]]),{'position_weights_by_length_beta':[[],[],[-1,1]]}), +]) +def test_paired_sum_direct_native_validation(args,opts): + from seqtree import _core + with pytest.raises((ValueError,TypeError)): + _core.gapblock_paired_sum_count_batch(*args,**opts) + + +def test_paired_sum_mass_uint64_and_puncture_with_zero_flanks(): + maximum=2**31-1 + assert gapblock.paired_sum_count_batch(['AC'],['CA'],['AC']*4,['CA']*4,[[maximum]])==( + [[4]],[[4*maximum]]) + # Identical weighted cores are not identical full strings in both lanes. + assert gapblock.paired_sum_count_batch(['AC'],['CA'],['CC'],['CA'],[[1]], + position_weights_by_length_alpha={2:[0,1]},exclude_exact=True)==([[1]],[[1]]) + + +def test_paired_sum_cell_addition_is_int64(): + # No group bound or first-lane bound rejects: only the actual lane sum exceeds int32. + weight=2**30 + options=dict(gap_open=0,gap_extend=0,position_weights_by_length_alpha={1:[weight]}, + position_weights_by_length_beta={1:[weight]}) + assert gapblock.paired_sum_count_batch(['A'],['A'],['C'],['C'],[[2**31-1]],**options)==([[0]],[[0]]) diff --git a/tests/python/test_gapblock_position_weights.py b/tests/python/test_gapblock_position_weights.py new file mode 100644 index 0000000..c7fe6ee --- /dev/null +++ b/tests/python/test_gapblock_position_weights.py @@ -0,0 +1,92 @@ +"""Synthetic correctness checks for the queued positional single-block API only.""" +import numpy as np +import pytest +from seqtree import Index,SearchParams,SubstitutionMatrix,PositionalMatrix,gapblock +B=SubstitutionMatrix.blosum62() + +def exhaustive(q,r,weights,go=37,ge=11,prior=None): + length=min(len(q),len(r));longer=max(len(q),len(r));d=longer-length + scores=[] + for start in range(length+1): + value=(go+(d-1)*ge if d else 0)+(prior(start,d,longer) if prior and d else 0) + for j in range(length): + qi=j+(d if len(q)>len(r) and j>=start else 0) + ri=j+(d if len(r)>len(q) and j>=start else 0) + value+=weights[longer][max(qi,ri)]*B.penalty(q[qi],r[ri]) + scores.append(value) + return min(scores) + +@pytest.mark.parametrize('prior',[None,lambda i,d,m:3*abs(i-2)]) +def test_positions_symmetry_reducers(prior): + queries=['','A','AC','ACA','CACA','ACAAA'] + refs=queries+['ACA'] + weights={n:[1+j*17 for j in range(n)] for n in range(6)} + options=dict(matrix=B,gap_open=37,gap_extend=11,gap_prior=prior,position_weights_by_length=weights) + expected=np.array([[exhaustive(q,r,weights,prior=prior) for r in refs] for q in queries]) + for threads in (1,2): + observed=np.asarray(gapblock.score_matrix(queries,refs,threads=threads,**options)) + np.testing.assert_array_equal(observed,expected) + np.testing.assert_array_equal(np.asarray(gapblock.score_matrix(refs,queries,threads=threads,**options)).T,expected) + thresholds=[[-1,0,37,100,500,100,10000] for q in queries] + counts=gapblock.count_batch(queries,refs,thresholds,threads=threads,exclude_exact=True,**options) + balls=gapblock.ball_batch(queries,refs,[100]*len(queries),threads=threads,exclude_exact=True,**options) + for i,q in enumerate(queries): + assert counts[i]==[sum(q!=r and expected[i,j]<=cut for j,r in enumerate(refs)) for cut in thresholds[i]] + assert [(h.ref_id,h.score) for h in balls[i]]==[(j,int(expected[i,j])) for j,r in enumerate(refs) if q!=r and expected[i,j]<=100] + assert all(h.n_subs==h.n_ins==h.n_dels==0 for h in balls[i]) + + +def test_equal_length_positional_matrix(): + queries=['CASSLGQAYEQYF','CASALGQAYEQYF'] + refs=['CASSLGQAYEQYF','AASSLGQAYEQYF','CASSAGQAYEQYF','CASSLGQAYEQAF'] + length=len(queries[0]);weights={length:[117,62,56,85,117,156,153,152,98,62,71,74,98]} + p=SearchParams(max_subs=5,max_ins=0,max_dels=0,max_total_edits=5,engine='seqtm');p.pos_matrix=PositionalMatrix.from_weights(B,weights[length]) + hits=Index.build(refs,'aa').search_batch(queries,p,threads=1) + observed=np.asarray(gapblock.score_matrix(queries,refs,matrix=B,position_weights_by_length=weights,threads=1)) + for i,row in enumerate(hits): + for hit in row: + assert observed[i,hit.ref_id]==hit.score + +@pytest.mark.parametrize('weights',[{2:[1]},{2:[1,-1]},{2:[1,1.5]},{2:[True,1]},{2:[1,2**31]},{3:[1,1,1]}, {2:[2**31-1]*2}]) +@pytest.mark.parametrize('api',['score_matrix','count_batch','ball_batch']) +def test_validation(weights,api): + args=(['AA'],['AC']) + if api=='count_batch':args+=([[100]],) + if api=='ball_batch':args+=([100],) + with pytest.raises((ValueError,TypeError,OverflowError)): + getattr(gapblock,api)(*args,matrix=B,position_weights_by_length=weights,threads=1) + + +def test_defaults_zero_weight_identity_duplicates_empty(): + queries=['AA'];refs=['AA','AC','AC'] + plain=np.asarray(gapblock.score_matrix(queries,refs,matrix=B,threads=1)) + np.testing.assert_array_equal(plain,np.asarray(gapblock.score_matrix(queries,refs,matrix=B,position_weights_by_length={},threads=2))) + assert gapblock.count_batch(queries,refs,[[0]],matrix=B,position_weights_by_length={2:[0,0]},exclude_exact=True,threads=1)==[[2]] + assert gapblock.count_batch([],refs,[],matrix=B,threads=1)==[] + assert gapblock.count_batch(queries,[],[[0]],matrix=B,threads=1)==[[0]] + assert gapblock.ball_batch(queries,refs,[-1],matrix=B,threads=1)==[[]] + with pytest.raises(ValueError):gapblock.ball_batch(queries,refs,[],threads=1) + + +def test_matrix_negative_threads_preserves_automatic_concurrency(): + q=['AC','CA'];r=['AA','AC','CA'] + for weights in (None,{2:[3,7]}): + args=dict(matrix=B,position_weights_by_length=weights) + expected=np.asarray(gapblock.score_matrix(q,r,threads=1,**args)) + np.testing.assert_array_equal(expected,np.asarray(gapblock.score_matrix(q,r,threads=-1,**args))) + + +@pytest.mark.parametrize('api', ['score_matrix', 'count_batch', 'ball_batch']) +def test_zero_weights_accept_large_safe_matrix_penalties(api): + from seqtree import amino_acids + width = len(amino_acids()) + matrix = SubstitutionMatrix.from_similarity( + [[0 if i == j else -600_000_000 for j in range(width)] for i in range(width)]) + kwargs = dict(matrix=matrix, gap_open=0, gap_extend=0, + position_weights_by_length={2: [0, 0]}, threads=1) + if api == 'score_matrix': + assert np.asarray(gapblock.score_matrix(['AA'], ['AC'], **kwargs)).tolist() == [[0]] + elif api == 'count_batch': + assert gapblock.count_batch(['AA'], ['AC'], [[0]], **kwargs) == [[1]] + else: + assert [(h.ref_id, h.score) for h in gapblock.ball_batch(['AA'], ['AC'], [0], **kwargs)[0]] == [(0, 0)] diff --git a/tests/python/test_gapblock_ranked.py b/tests/python/test_gapblock_ranked.py index fc993ef..5311998 100644 --- a/tests/python/test_gapblock_ranked.py +++ b/tests/python/test_gapblock_ranked.py @@ -20,13 +20,13 @@ def test_topk_and_counts_equal_exhaustive_matrix(matrix,prior,exclude): for k in (1,5,10,100): expected=[sorted([(j,sm[i,j]) for j,s in enumerate(r) if not exclude or s.upper()!=x.upper()], key=lambda h:(h[1],h[0]))[:k] for i,x in enumerate(q)] - for threads in (1,4,0): + for threads in (1,2,0): assert tuples(topk_batch(q,r,k,**opts,threads=threads,exclude_exact=exclude))==expected thresholds=[[30,0,-1,9,9,200] if i%2 else [] for i in range(len(q))] expected=[[sum(sm[i,j]<=t for j,s in enumerate(r) if not exclude or s.upper()!=x.upper()) for t in thresholds[i]] for i,x in enumerate(q)] assert count_batch(q,r,thresholds,**opts,threads=1,exclude_exact=exclude)==expected - assert count_batch(q,r,thresholds,**opts,threads=4,exclude_exact=exclude)==expected + assert count_batch(q,r,thresholds,**opts,threads=2,exclude_exact=exclude)==expected def test_paired_topk_max_lane_oracle_and_joint_exact_puncture(): @@ -38,7 +38,7 @@ def test_paired_topk_max_lane_oracle_and_joint_exact_puncture(): expected=[sorted([(j,a[i,j],b[i,j]) for j in range(len(ra)) if not exclude or (qa[i].upper(),qb[i].upper())!=(ra[j].upper(),rb[j].upper())], key=lambda h:(max(h[1],h[2]),h[0]))[:3] for i in range(len(qa))] - for threads in (1,4): + for threads in (1,2): assert paired_topk_batch(qa,qb,ra,rb,3,**opts,threads=threads,exclude_exact=exclude)==expected @@ -98,7 +98,7 @@ def test_paired_count_matches_all_matrix_ties_and_threshold_order(prior): expected=[[sum(max(a[i,j],b[i,j])<=t for j in range(len(ra)) if not exclude or (qa[i],qb[i])!=(ra[j],rb[j])) for t in ts] for i,ts in enumerate(thresholds)] - for threads in (1,4): + for threads in (1,2): assert paired_count_batch(qa,qb,ra,rb,thresholds,**opts,threads=threads,exclude_exact=exclude)==expected selected=paired_topk_batch(qa,qb,ra,rb,3,**opts,exclude_exact=exclude) expected_top=[[j for j in sorted(range(len(ra)),key=lambda j:(max(a[i,j],b[i,j]),j)) diff --git a/tests/python/test_gapblock_sequence_weights.py b/tests/python/test_gapblock_sequence_weights.py new file mode 100644 index 0000000..df0574a --- /dev/null +++ b/tests/python/test_gapblock_sequence_weights.py @@ -0,0 +1,54 @@ +"""Per-input positional geometry: exhaustive symmetric oracle and shared reductions.""" +import numpy as np +import pytest +from seqtree import gapblock + + +def oracle(q, r, qw, rw, go=7, ge=3): + d = abs(len(q)-len(r)) + scores = [] + for start in range(min(len(q), len(r))+1): + score = go+(d-1)*ge if d else 0 + for i in range(min(len(q), len(r))): + qi = i+(d if len(q)>len(r) and i>=start else 0) + ri = i+(d if len(r)>len(q) and i>=start else 0) + score += (q[qi]!=r[ri])*max(qw[qi], rw[ri]) + scores.append(score) + return min(scores) + + +def test_sequence_weights_oracle_and_linked_reductions(): + q = ['', 'AC', 'ACA', 'CAAA'] + r = ['', 'CA', 'ACA', 'ACAAA', 'CA'] + qw = [[], [10,100], [100,10,10], [10,100,100,10]] + rw = [[], [100,10], [10,100,10], [10,10,100,100,10], [10,100]] + expected = np.array([[oracle(a,b,aw,bw) for b,bw in zip(r,rw)] for a,aw in zip(q,qw)]) + cuts = [[-1,0,10,100,250,100] for _ in q] + opts = dict(gap_open=7, gap_extend=3, query_position_weights=qw, reference_position_weights=rw) + for threads in (1,2): + np.testing.assert_array_equal(np.asarray(gapblock.score_matrix(q,r,threads=threads,**opts)), expected) + reverse = gapblock.score_matrix(r,q,gap_open=7,gap_extend=3,threads=threads, + query_position_weights=rw,reference_position_weights=qw) + np.testing.assert_array_equal(np.asarray(reverse).T,expected) + counts,masses = gapblock.count_batch(q,r,cuts,threads=threads,linear_mass=True,exclude_exact=True,**opts) + paircounts,pairmasses = gapblock.paired_sum_count_batch(q,q,r,r,cuts,gap_open=7,gap_extend=3, + threads=threads,exclude_exact=True,query_position_weights_alpha=qw, + query_position_weights_beta=qw,reference_position_weights_alpha=rw,reference_position_weights_beta=rw) + for i,a in enumerate(q): + eligible = np.array([a!=b for b in r]) + assert counts[i] == [int(np.sum(eligible & (expected[i]<=c))) for c in cuts[i]] + assert masses[i] == [int(np.sum(np.maximum(0,c-expected[i][eligible]))) for c in cuts[i]] + assert paircounts[i] == [int(np.sum(eligible & (2*expected[i]<=c))) for c in cuts[i]] + assert pairmasses[i] == [int(np.sum(np.maximum(0,c-2*expected[i][eligible]))) for c in cuts[i]] + uniform = dict(query_position_weights=[[100]*len(a) for a in q],reference_position_weights=[[100]*len(b) for b in r]) + old = gapblock.score_matrix(q,r,gap_open=7,gap_extend=3,position_weights_by_length={n:[100]*n for n in range(6)}) + np.testing.assert_array_equal(np.asarray(old),np.asarray(gapblock.score_matrix(q,r,gap_open=7,gap_extend=3,**uniform))) + + +@pytest.mark.parametrize('qw,rw,by_length', [([[1]],[[1,1]],None), ([[1,1]],None,None), + ([[1,-1]],[[1,1]],None), ([[True,1]],[[1,1]],None), ([[1,1.5]],[[1,1]],None), + ([[2**31-1]*2],[[1,1]],None), ([[1,1]],[[1,1]],{2:[1,1]})]) +def test_sequence_weight_validation(qw,rw,by_length): + with pytest.raises((ValueError,TypeError,OverflowError)): + gapblock.count_batch(['AA'],['AC'],[[100]],query_position_weights=qw, + reference_position_weights=rw,position_weights_by_length=by_length) diff --git a/tests/python/test_pairwise.py b/tests/python/test_pairwise.py index 06150be..54bbdcd 100644 --- a/tests/python/test_pairwise.py +++ b/tests/python/test_pairwise.py @@ -282,7 +282,7 @@ def test_score_matrix_is_thread_invariant(): q = [rand_seq(rng, 8, 20) for _ in range(24)] r = [rand_seq(rng, 8, 20) for _ in range(12)] ref = [score_matrix(q, r, m, threads=1).row(i) for i in range(len(q))] - for threads in (2, 4, 0): + for threads in (1, 2, 0): got = score_matrix(q, r, m, threads=threads) assert [got.row(i) for i in range(len(q))] == ref @@ -340,7 +340,7 @@ def test_bad_arguments_raise(): # Two crashes found by an adversarial audit, both of them regressions from the overflow fix # --------------------------------------------------------------------------------------------- -@pytest.mark.parametrize("threads", [1, 0, 4]) +@pytest.mark.parametrize("threads", [1, 0, 2]) def test_an_overflow_in_a_worker_raises_instead_of_killing_the_process(threads): """The batch paths run on std::thread. An exception escaping a thread entry function calls std::terminate -- so the int64 narrowing guard, which raises cleanly on the scalar path, was diff --git a/tests/python/test_setup.py b/tests/python/test_setup.py new file mode 100644 index 0000000..a925d28 --- /dev/null +++ b/tests/python/test_setup.py @@ -0,0 +1,27 @@ +"""Setup combines extras and never reports success after a failed install.""" +import os +from pathlib import Path +import subprocess + +import pytest + + +@pytest.mark.skipif(os.name == "nt", reason="POSIX setup script") +def test_setup_extras_and_install_failure(tmp_path): + uv = tmp_path / "uv" + uv.write_text('#!/bin/sh\nprintf "%s\\n" "$*"\nexit "${TEST_UV_STATUS:-0}"\n') + uv.chmod(0o755) + env = dict(os.environ, PATH=str(tmp_path) + os.pathsep + os.environ["PATH"]) + root = tmp_path / "project" + root.mkdir() + (root / "setup.sh").write_text((Path(__file__).resolve().parents[2] / "setup.sh").read_text()) + activation = root / ".venv" / "bin" / "activate" + activation.parent.mkdir(parents=True) + activation.write_text("") + ok = subprocess.run(["sh", "setup.sh", "--tests", "--bench"], cwd=root, + env=env, capture_output=True, text=True, timeout=10) + assert ok.returncode == 0 and "pip install -e .[test,bench]" in ok.stdout + failed = subprocess.run(["sh", "setup.sh", "--tests"], cwd=root, + env=dict(env, TEST_UV_STATUS="42"), capture_output=True, + text=True, timeout=10) + assert failed.returncode == 42 and "Done." not in failed.stdout diff --git a/tests/python/test_text_evalue.py b/tests/python/test_text_evalue.py new file mode 100644 index 0000000..4ff5ca5 --- /dev/null +++ b/tests/python/test_text_evalue.py @@ -0,0 +1,149 @@ +"""E-values for substring hits. + +The C++ suite owns the completeness of the search itself. These tests own the +calibration layer: that a planted motif is significant while background is not, that +counts collapse dependent offsets onto records, that the ball nests, and that a count +which is really a floor is refused rather than reported. +""" +import random + +import pytest + +import seqtree + +ALPHABET = "ACDEFGHIKLMNPQRSTVWY" +MOTIF = "WYWMHCWY" + + +def _records(n, length, seed): + rng = random.Random(seed) + return ["".join(rng.choice(ALPHABET) for _ in range(length)) for _ in range(n)] + + +@pytest.fixture(scope="module") +def planted(): + """A target where 40 of 400 records carry a motif, and a matched control.""" + target = _records(400, 120, seed=1) + for i in range(0, 40): + target[i] = target[i][:50] + MOTIF + target[i][50 + len(MOTIF):] + control = _records(400, 120, seed=2) + return (seqtree.TextIndex.build(target, "aa", k=4), + seqtree.TextIndex.build(control, "aa", k=4)) + + +def test_planted_motif_is_significant_and_background_is_not(planted): + target, control = planted + background = _records(1, 8, seed=99)[0] + results = seqtree.text_evalues(target, control, [MOTIF, background], max_subs=1) + + planted_result, background_result = results + assert planted_result["n_target"] >= 40 + assert planted_result["p_enrichment"] < 1e-6 + assert background_result["p_enrichment"] > 0.01 + + +def test_counts_are_records_not_offsets(): + """A motif repeated inside one record counts once, not once per occurrence. + + Overlapping occurrences in a repeat are strongly dependent; counting them + separately is what makes the Poisson estimator over-disperse. + """ + repeat = "AAAA" + MOTIF + "AAAA" + MOTIF + "AAAA" + MOTIF + "AAAA" + target = seqtree.TextIndex.build([repeat], "aa", k=4) + control = seqtree.TextIndex.build(_records(50, 60, seed=3), "aa", k=4) + + hits = target.search_batch([MOTIF], max_subs=0) + assert hits.num_hits == 3 + + result = seqtree.text_evalues(target, control, [MOTIF], max_subs=0)[0] + assert result["n_target"] == 1 + + +def test_evalue_is_monotone_in_the_radius(planted): + target, control = planted + previous = -1.0 + for radius in (0, 1, 2): + result = seqtree.text_evalues(target, control, [MOTIF], max_subs=radius)[0] + assert result["E"] >= previous + previous = result["E"] + + +def test_exclude_exact_drops_only_zero_mismatch_records(planted): + target, control = planted + with_exact = seqtree.text_evalues(target, control, [MOTIF], max_subs=1)[0] + without = seqtree.text_evalues( + target, control, [MOTIF], max_subs=1, exclude_exact=True + )[0] + assert without["n_target"] < with_exact["n_target"] + + +def test_truncated_search_is_refused(planted): + """A capped count is a floor; an E-value computed from it would be wrong.""" + target, control = planted + capped = target.search_batch([MOTIF], max_subs=1, max_hits=1) + assert any(memoryview(capped.arrays()["truncated"]).tolist()) + with pytest.raises(ValueError, match="truncated"): + capped.record_distances() + + +def test_empty_control_is_rejected(planted): + target, _ = planted + empty = seqtree.TextIndex.build([], "aa", k=4) + with pytest.raises(ValueError, match="control text is empty"): + seqtree.text_evalues(target, empty, [MOTIF], max_subs=0) + + +def test_threshold_round_trips_to_the_requested_evalue(planted): + target, control = planted + queries = [MOTIF, _records(1, 10, seed=7)[0]] + cutoffs = seqtree.text_threshold_for_evalue( + target, control, queries, max_subs=2, e_target=1.0 + ) + assert len(cutoffs) == len(queries) + for query, cutoff in zip(queries, cutoffs): + if cutoff < 0: + continue + achieved = seqtree.text_evalues(target, control, [query], max_subs=cutoff)[0] + assert achieved["E"] <= 1.0 + 1e-9 + + +def test_threshold_is_unreachable_when_the_control_is_too_small(planted): + """e_target below 3N/M cannot be certified however tight the ball.""" + target, control = planted + cutoffs = seqtree.text_threshold_for_evalue( + target, control, [MOTIF], max_subs=2, e_target=1e-9 + ) + assert cutoffs == [-1] + + +def test_negative_radius_is_rejected(planted): + target, control = planted + with pytest.raises(ValueError, match="max_subs must be >= 0"): + seqtree.text_threshold_for_evalue( + target, control, [MOTIF], max_subs=-1, e_target=1.0 + ) + + +def test_record_distances_and_exact_record_exclusion(): + # Build-time groups do not change the counting unit; a record containing an + # exact AND an inexact occurrence is wholly excluded by exclude_exact. + target = seqtree.TextIndex.build(["AAAAACAA", "AAAC", "AAAC"], k=2, group_ids=[0, 0, 0]) + control = seqtree.TextIndex.build(["CCCC", "CCCC", "CCCC"], k=2) + hits = target.search_batch(["AAAA"], max_subs=1, group_by=True) + assert hits.record_distances() == [[0, 1, 1]] + assert seqtree.text_evalues(target, control, iter(["AAAA"]), max_subs=1)[0]["n_target"] == 3 + assert seqtree.text_evalues(target, control, ["AAAA"], max_subs=1, exclude_exact=True)[0]["n_target"] == 2 + assert target.search_batch([], max_subs=1).record_distances() == [] + + +def test_calibration_requires_matching_alphabets_and_nonempty_control(): + target = seqtree.TextIndex.build(["AAAA"], k=2) + control = seqtree.TextIndex.build(["AAAA"], alphabet="nt", k=2) + for call in (seqtree.text_evalues, seqtree.text_threshold_for_evalue): + opts = dict(max_subs=1) + if call is seqtree.text_threshold_for_evalue: + opts["e_target"] = 1.0 + with pytest.raises(ValueError, match="alphabets"): + call(target, control, ["AAAA"], **opts) + with pytest.raises(ValueError, match="control text is empty"): + seqtree.text_threshold_for_evalue(target, seqtree.TextIndex.build([], k=2), ["AAAA"], 1, 1.0) diff --git a/tests/python/test_text_index.py b/tests/python/test_text_index.py index 5b9ecd5..008a0c5 100644 --- a/tests/python/test_text_index.py +++ b/tests/python/test_text_index.py @@ -94,7 +94,7 @@ def test_hit_order_is_deterministic_across_thread_counts(): ixr = TextIndex.build(refs, alphabet="aa", k=4) queries = [r[:12] for r in refs] base = ixr.search_batch(queries, max_subs=2, threads=1).to_numpy() - for t in (2, 4, 8): + for t in (1, 2, 0): got = ixr.search_batch(queries, max_subs=2, threads=t).to_numpy() for key in ("query_begin", "ref_id", "offset", "n_subs", "mm_pos"): assert (got[key] == base[key]).all(), key diff --git a/tests/python/test_unit_histogram_fastpath.py b/tests/python/test_unit_histogram_fastpath.py new file mode 100644 index 0000000..360fd8d --- /dev/null +++ b/tests/python/test_unit_histogram_fastpath.py @@ -0,0 +1,49 @@ +"""Exact installed-API oracle for the private unit histogram fast path.""" +import random +import pytest +from seqtree import Index, SearchParams, SubstitutionMatrix, PositionalMatrix + + +def _expected(index, queries, params, exclude): + out=[] + for hits in index.search_batch(queries,params,threads=1): + bins=[0]*(params.max_total_edits+1) + for hit in hits: + total=hit.n_subs+hit.n_ins+hit.n_dels + if not exclude or total:bins[total]+=1 + out.append(bins) + return out + + +def test_exhaustive_binary_and_duplicate_case_counts(): + sequences=['']+[''.join('A' if mask&(1<