Repository navigation
Fix zero-dilution in genetic_relatedness - #390
Merged
Merged
Conversation
…relatedness A sample set with no call at a site got frequency 0.0 instead of undefined in _relatedness_columns, so that fake 0.0 was centered against the across-set mean and summed into the relatedness matrix -- diluting any pair involving a set with missing data, and giving a misleading 0.0 rather than NaN to a pair sharing no called site at all. - pbar(a,s) now excludes sets with no call at a site from the mean, instead of zero-filling them in. - An invalid set's column entry is forced to exactly 0 after centering, so it contributes nothing to any product -- the same zero-contribution mechanism already used by ibs() and the distance_stats.py/decomposition.py pairwise fix. - A new per-pair joint-validity accumulator tracks whether any column ever had both sets of a pair valid; the result is NaN wherever it never did. Applied in both the eager and streaming paths. - _relatedness_columns now returns an additional valid_mask array. Tests: a unit test on _relatedness_columns' forced-zero/mean-exclusion behavior, an eager and a streaming regression test reproducing a pair with disjoint missingness, and an updated host-reference helper matching the corrected convention. Also fixes the now-stale `undefined` text for relatedness.genetic_relatedness in docs/dev/pipeline_map.toml.
…ck sizing The NaN-for-undefined-pairs fix added a valid_mask array and a second matmul (J += vm @ vm.T) alongside the existing per-chunk cols array and R matmul, but left n_intermediates at 3 in both the eager chunk-size and streaming block-size estimates. That understates the actual per-chunk/per-block GPU memory footprint, so chunks/blocks were sized larger than the intended memory_fraction budget allows. Bumps both to 4 (cols + valid_mask-as-float + the two matmul outputs), matching the same correction already made for grm's analogous fix.
andrewkern
requested changes
Nov 5, 2026
Member
There was a problem hiding this comment.
a few nitpicks
- Pairs still shrink with call rate: each pair is a sum over sites with one global denominator, so a pair with fewer called sites gets a smaller value. With each haplotype as its own set (sample_sets=None) and haplotypes 0-19 missing half their calls, their self-relatedness is half that of the complete haplotypes (ratio 0.50, vs 0.98 with complete data). This is the same dilution #391 fixes for grm with a per-pair count. Since J already counts the columns both sets call, scaling each pair by n_cols / J would remove it.
- With no column at all (e.g. missing_data='exclude' with no complete site), the result is now all NaN, also with proportion=False. tskit gives 0.0 for a tree sequence with no sites (a sum over no sites), and the streaming path still says "relatedness of nothing is 0" (relatedness.py:190). Either keep 0.0 when there is no column, or update that comment and note the change.
- Complete data is about 1.6x slower (4000 haplotypes x 100k sites: 0.62 s on main, 1.00 s here). The same fast path as in #391 applies: skip vm @ vm.T for a chunk where every set has a call.
- add a changelog entry: values change under missing data, and a pair with no shared called site is now NaN.
Addresses review feedback on this PR: - Dilution: R summed each pair over however many (site, allele) columns it jointly called, with no rescaling, so a pair that called fewer of the analyzed columns shrank with its own call rate (the same dilution already fixed for grm's per-pair SNP count, and for pairwise_distance_matrix's per-pair site count). Each pair's sum is now rescaled by n_cols / J -- its own jointly-called column count out of the total analyzed -- before any proportion/span normalization, both in the eager and streaming paths. - No column at all (an empty region, or 'exclude' with no complete site) is an empty analysis, not an undefined pair: it returns 0.0 (a sum over no sites), unless proportion=True, where the segregating-site denominator is also 0 and the result is NaN -- matching tskit's own division (out /= denominator gives nan, not 0, when both are 0). A pair with no shared column while other data exists is still NaN either way. - Complete-data chunks skip the validity matmul (same fast path as the preceding grm/ibs commits), recovering the speed lost to the NaN-for- undefined-pairs fix: 4000 haplotypes x 100k sites went from 1.00s back to 0.69s. Tests: hand-verified pin values updated for the rescaling (worked through allele-by-allele), a self-relatedness dilution regression reproducing the reported 0.50 ratio (now ~1.0), and zero-column regressions for both proportion settings on both the eager and streaming paths. Also fixes the pipeline_map.toml undefined/notes text for genetic_relatedness and adds a changelog note.
…rable span_normalize=True silently returned the raw, un-normalized sum when the span was unmeasurable (e.g. a degenerate chrom_start/chrom_end), instead of NaN. Every other span-normalized statistic in the library goes through diversity._apply_span_normalize, which already returns NaN in that case; genetic_relatedness had its own hand-rolled copy of the same division that skipped the else branch. Fixed in both the eager and streaming paths.
Collaborator
Author
|
ok, done, incidental fix for zero-span case too |
Member
|
this looks GTG |
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 #337
A sample set with no call at a site got frequency 0.0 instead of undefined in _relatedness_columns, so that fake 0.0 was centered against the across-set mean and summed into the relatedness matrix -- diluting any pair involving a set with missing data, and giving a misleading 0.0 rather than NaN to a pair sharing no called site at all.
Tests: a unit test on _relatedness_columns' forced-zero/mean-exclusion behavior, an eager and a streaming regression test reproducing a pair with disjoint missingness, and an updated host-reference helper matching the corrected convention. Also fixes the now-stale
undefinedtext for relatedness.genetic_relatedness in docs/dev/pipeline_map.toml.