Repository navigation
Conversation
Pure refactor, no behavior change: extracts _dxy_per_site_from_counts (one call per side, dxy's own long-standing rule) out of the logic dxy() had inline, and switches dxy(), _per_variant_dxy, the windowed scatter-sum engine's dxy output, and the fused kernel's eventual dxy accumulator onto it, instead of reusing _site_components_from_counts/_twopop_site_components for both purposes. All existing tests pass unchanged; the per-site dxy values computed here are bit-identical to before (same gate, same formula, just extracted).
Hudson FST and Da were biased upward under missing data when a population has only one call at a site: _site_components_from_counts counted a site in the between-pop term wherever each population had at least one call, but zeroed a population's own within-pop term wherever it had fewer than two -- so Hb stayed inflated relative to Hw. Both now require both populations to have at least two calls (their own pair), the same site set for Hb and Hw; this also fixes fst_tskit and pbs, which share the same per-pair helper. The fused kernel gets the same gate for its t_fst_num/t_fst_den/t_pi1/ t_pi2 accumulators (t_dxy was already decoupled onto its own looser gate by the preceding commit), and both fused da_sum computations (windowed and chunked) switch from the dxy accumulator to the fst one, since da's between term now needs the stricter site set too. Tests: a hand-verified three-site pin (one site drops out entirely once its one-call population is excluded), a reproduction of the reported MCAR-missingness bias staying within tolerance after the fix, and a flip of test_single_sample_population's assertion -- a single-haplotype population can never have its own pair, so fst_hudson is now correctly always NaN there rather than incidentally in [0, 1]. Also updates the now-stale pipeline_map.toml undefined/notes text for fst_hudson, da, and pbs.
The scalar-level fix and its tests (previous commit) don't exercise the engines that recompute their own per-site components independently of _site_components_from_counts, most importantly the fused CUDA kernel, which has its own copy of the joint_pairs logic in C++. Add a class that runs the one-call-pair reproducer (a small, heavily-missing population) through the host scatter engine and both fused kernels (non-chunked and chunked with a forced small chunk size, so the fst/da-vs-dxy accumulator split is checked across chunk boundaries too), and checks each against the scalar fst_hudson/da/dxy on the same window's data. Verified these catch the regression: reverting the kernel's joint_pairs gate alone (leaving the Python-side fix in place) fails both fused tests while leaving the scatter test green, since the latter routes through the same Python helper as the scalar reference and so cannot independently catch a kernel-only divergence.
pbs() combines three pairwise Hudson FSTs into one branch-length formula, which requires all three to be computed over the same site set -- its own comment already said so. The upstream all_three prefilter required only one call per population (matching _site_components_from_counts's old gate), but the preceding commits tightened that shared helper's gate to require each population's own pair (>=2 calls). Left unmatched, a site with exactly one call in one population still passed the stale prefilter, then got asymmetrically dropped from the two pairwise FSTs touching that population while staying in the third -- exactly the "enters some branches and not others" failure the comment warns against, introduced as a side effect of the earlier fix rather than present before it. Bumped all_three to require >1 call per population, matching the gate _hudson_num_den already enforces per pair internally. Tests: parametrized test mirroring the existing wholly-missing-population regression, but leaving exactly one call in the gapped population instead of zero -- confirmed it fails 2 of 3 parametrizations without this fix.
andrewkern
reviewed
Nov 6, 2026
andrewkern
left a comment
andrewkern
left a comment
Member
There was a problem hiding this comment.
This looks great to me. A couple findings:
divergence.dxy_components(divergence.py:47) still carries its own per-site Dxy with the one-call gate. Nothing calls it internally, but having it use_dxy_per_site_from_countswould keep the Dxy rule in one place.- Optional:
missing_data.rstcould say that FST and da count a site only where each population has at least two calls, which explains why a single-haplotype population now gives NaN.
andrewkern
approved these changes
Nov 6, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #338
Hudson FST, Da, and the tskit-style FST estimator were biased upward under missing data whenever a population had only one call at a site, because the shared per-site helper (_site_components_from_counts) let such a site contribute to the between-population term (Hb) while silently zeroing that population's own within-population term (Hw) instead of excluding the site entirely (Hb and Hw were summed over different, inconsistent site sets).
The fix requires each population to have its own pair (at least two calls) before a site counts toward Hb or Hw at all, applied consistently across the scalar functions, the host scatter engine, and both fused CUDA kernels; this transitively fixes fst_hudson, fst_tskit, da, and pbs since they all share the same helper, while dxy is decoupled onto its own long-standing, correctly looser one-call-per-side gate via a newly extracted _dxy_per_site_from_counts.
pbs()'s own prefilter for keeping all three pairwise branches on a common site set needs tightened to match, so a site with exactly one call in one population could still asymmetrically enter the branch that didn't involve it; that's fixed and tested too.