Make undefined Fst NaN, and fix multi-population single-stat column names - #311
Conversation
| **kwargs, | ||
| ) | ||
| if df2 is not None and not df2.empty: | ||
| for col in df2.columns: |
There was a problem hiding this comment.
This assumes df2 has the same rows as df1. With 3 populations, ['pi', 'fst'], and one empty window this raises ValueError: Length of values (4) does not match length of index (6). That call worked before this PR so this introduces a new bug
| haplotype_matrix, window_size, step_size, | ||
| single_part, populations, missing_data, span_normalize, | ||
| chrom=chrom) | ||
| rest = sorted(requested - set(single_part)) |
There was a problem hiding this comment.
sorted here breaks on a callable statistic. ['pi', 'fst', my_fn] with two populations raises TypeError: '<' not supported between instances of 'str' and 'function'. Before this PR that request fell through to the fallback, which accepts callables.
| # remainder excludes every scatter_single member by construction, so | ||
| # it can never re-enter this branch. | ||
| single_part = sorted(requested & scatter_single) | ||
| if single_part and n_pops >= 2: |
There was a problem hiding this comment.
This covers scatter_single. garud_h1 or mean_nsl with two populations still comes back as one column holding the first population's value, with no warning. In #200 we say we want one rule in every engine. I found this by running quickstart's everything-at-once example which returns pi_pop1, pi_pop2 next to a bare garud_h1.
| df = _windowed_thetas_scatter( | ||
| haplotype_matrix, window_size, step_size, | ||
| statistics, [pop], missing_data, span_normalize, chrom=chrom) | ||
| if df is None: |
There was a problem hiding this comment.
_windowed_thetas_scatter doesn't ever seem to return None? It returns an empty frame when a population has no sites under 'exclude'. Then this loop just skips that population, so you get pi_p1 and no pi_p2 with no error, i think? We should probably give the missing population a NaN column or raise an error.
| stat: f"{stat}_{pop}" for stat in statistics if stat in df.columns | ||
| }) | ||
| if merged is None: | ||
| merged = df |
There was a problem hiding this comment.
The merged frame keeps only the first population's n_variants. Under 'exclude' each population is filtered to its own site set, so a danger is that that count could be wrong for the other columns.
Also, this copy-columns loop looks to be repeated a few times in this file. This could be a small helper that checks row counts to cover all cases.
There was a problem hiding this comment.
in #310 I made the site set shared across populations (same site set for dxy, pi1, pi2). But perhaps that's the wrong choice?
There was a problem hiding this comment.
ah okay! maybe i should have reviewed #310 first. using the same site set seems a reasonable choice to me. sorry about that
There was a problem hiding this comment.
not sure order matters (there are likely bugs there too!) but it's probably easiest to keep this PR scoped to just the NaNs for now, and think more carefully about the other half of the issue (in another PR)
| if cp.any(valid_mask): | ||
| return float((cp.sum(num[valid_mask]) / cp.sum(den[valid_mask])).get()) | ||
| return 0.0 | ||
| return float('nan') |
There was a problem hiding this comment.
okay so now fst(method='hudson') disagrees with fst_nei and fst_tskit which still return 0.0 for this case - can we add those two as well to this PR?
| assert np.isclose(df["pi_pop1"].iloc[0], pi1_ref, rtol=RTOL, atol=ATOL) | ||
| assert np.isclose(df["pi_pop2"].iloc[0], pi2_ref, rtol=RTOL, atol=ATOL) | ||
| assert np.isclose(df["fst_wc"].iloc[0], fst_wc_ref, rtol=RTOL, atol=ATOL) | ||
| # The regression this guards against silently dropped pop2 and returned |
There was a problem hiding this comment.
This comment describes the old bug. Say what the check is for instead: pi must differ between the two populations or a collapsed column would go unnoticed.
|
Looking into this: the issues here largely stem from trying to merge the python-loop fallback (skips windows that are empty) with the scatter/fused engines (emit empty windows with nan or 0). That's different numbers of rows with empty windows, so problems ensue. Whether there should actually be this asymmetry in behavior is a good question, and I don't think I'll try to fix it here. I'll see if I can salvage this PR for the original NaN-vs-zero fix, if not will close and try to break up into multiple issues. |
|
There current plan here is to rescope to replacing 0.0 with NaN, including fst_tskit, after #310 is merged; and fix the other half of the issue in another PR. |
fst_hudson, fst_weir_cockerham, fst_tskit, and fst_nei all returned 0.0 for an undefined ratio (no site with data in both populations), while every windowed engine already used NaN for the same case. Standardize on NaN.
2a28efb to
44d84c2
Compare
|
I'm merging this as it's trivial, will deal with the harder part in a followup |
2 participants
Closes #200
fst_hudson/fst_weir_cockerhamreturned 0.0 for an undefined ratio (no site with data in both populations, or no polymorphism at all) while every windowed engine's own Fst math already usednanfor the same case; the per-window Python-loop fallback inherited the scalar 0.0, so it silently disagreed with the fused/scatter engines for the identical window. Standardize onnan, matchingdivergence.gmin's existing convention.Separately,
windowed_analysis()collapsed a single-population statistic (e.g. pi) onto one bare column whenever 2+ populations were named alongside a two-population statistic, silently keeping only the first population's value -- both under 'include' (via the fused engine'spopulation=pop1) and under 'exclude' forfst/fst_hudson/dxy/da(via a scatter dispatch branch that also only read the first population). Route any single-pop statistic requested with 2+ named populations through a new per-population loop that suffixes each column "_", matching the naming already used elsewhere, and recurse for whatever else was requested.