From b767e23bbe97322db21e6bd4330b04518509ac3f Mon Sep 17 00:00:00 2001 From: Tim Millar Date: Thu, 11 Jul 2019 08:06:25 +1200 Subject: [PATCH 1/2] fix for floating point bug in heterozygosity_expected --- allel/stats/hw.py | 3 ++- allel/test/test_stats.py | 5 ++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/allel/stats/hw.py b/allel/stats/hw.py index 1ef7e20f..1212584e 100644 --- a/allel/stats/hw.py +++ b/allel/stats/hw.py @@ -93,7 +93,8 @@ def heterozygosity_expected(af, ploidy, fill=np.nan): out = 1 - np.sum(np.power(af, ploidy), axis=1) # fill values where allele frequencies could not be calculated - af_sum = np.sum(af, axis=1) + precision = np.finfo(af.dtype).precision + af_sum = np.round(np.sum(af, axis=1), precision) with ignore_invalid(): out[(af_sum < 1) | np.isnan(af_sum)] = fill diff --git a/allel/test/test_stats.py b/allel/test/test_stats.py index b2ee4b6e..6dad2e79 100644 --- a/allel/test/test_stats.py +++ b/allel/test/test_stats.py @@ -425,7 +425,10 @@ def refimpl(f, ploidy, fill=0): """Limited reference implementation for testing purposes.""" # check allele frequencies sum to 1 - af_sum = np.sum(f, axis=1) + # round to suitable precision for the dtype + af_sum = (np.sum(f, axis=1)) + precision = np.finfo(af_sum.dtype).precision + af_sum = np.round(np.sum(f, axis=1), precision) # assume three alleles p = f[:, 0] From 7de90cc6881695f6f44fae7fbeb4d2020a8713aa Mon Sep 17 00:00:00 2001 From: Tim Millar Date: Thu, 11 Jul 2019 08:33:51 +1200 Subject: [PATCH 2/2] whitespace --- allel/test/test_stats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/allel/test/test_stats.py b/allel/test/test_stats.py index 6dad2e79..2825e28a 100644 --- a/allel/test/test_stats.py +++ b/allel/test/test_stats.py @@ -425,7 +425,7 @@ def refimpl(f, ploidy, fill=0): """Limited reference implementation for testing purposes.""" # check allele frequencies sum to 1 - # round to suitable precision for the dtype + # round to suitable precision for the dtype af_sum = (np.sum(f, axis=1)) precision = np.finfo(af_sum.dtype).precision af_sum = np.round(np.sum(f, axis=1), precision)