Time windows in AFS #3266

tforest · 2025-09-09T21:16:05Z

Description

This is a WIP and replaces #2948 as a more specific PR regarding time-windows in AFS.

In fine, users will be able to use ts.allele_frequency_spectrum(time_windows=tw), with tw being a list of any non-overlapping and strictly increasing time intervals in [0, inf].

codecov · 2025-09-09T21:21:53Z

Codecov Report

❌ Patch coverage is 97.00000% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.80%. Comparing base (8b615e5) to head (33f990d).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
c/tskit/trees.c	94.33%	0 Missing and 3 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3266      +/-   ##
==========================================
+ Coverage   89.75%   89.80%   +0.04%     
==========================================
  Files          29       29              
  Lines       32546    32602      +56     
  Branches     5962     5973      +11     
==========================================
+ Hits        29212    29278      +66     
+ Misses       1890     1884       -6     
+ Partials     1444     1440       -4

Flag	Coverage Δ
c-tests	`86.85% <94.73%> (+0.08%)`	⬆️
lwt-tests	`80.38% <ø> (ø)`
python-c-tests	`88.25% <100.00%> (+0.01%)`	⬆️
python-tests	`98.79% <100.00%> (+<0.01%)`	⬆️
python-tests-no-jit	`32.95% <16.66%> (-0.08%)`	⬇️
python-tests-numpy1	`50.57% <16.66%> (-0.13%)`	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
c/tskit/core.c	`95.54% <100.00%> (+0.01%)`	⬆️
c/tskit/core.h	`100.00% <ø> (ø)`
python/_tskitmodule.c	`88.25% <100.00%> (+0.01%)`	⬆️
python/tskit/trees.py	`98.84% <100.00%> (+0.01%)`	⬆️
c/tskit/trees.c	`91.01% <94.33%> (+0.19%)`	⬆️

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

petrelharp · 2025-09-14T15:52:16Z

Note: Along the way, this also allows one to call

ts.allele_frequency_spectrum(sample_sets=[0,1,2])

to produce the same thing as

ts.allele_frequency_spectrum(sample_sets=[[0,1,2]])

This slightly breaks with the convention for other stats but we'd like to do it anyhow because:

for other stats, where this would produce differently-shaped output, there is a dimension to drop (that corresponding to indexes) that the AFS does not have
it's annoying/confusing to have to do the two sets of brackets, and the error message isn't very informative, so this could be a sticking point for newbies

petrelharp · 2025-09-16T20:10:53Z

There's one line with partial coverage that I've confirmed using tsk_bug_assert( )s is being covered by the python tests:

I've added another C test that might satisfy this? But this is fine.

petrelharp · 2025-09-16T23:54:54Z

The last thing to do is to put in some tests to compare versions of the method (eg decap-based and not) for a few different time windows.

jeromekelleher · 2025-09-17T08:05:58Z

python/tests/test_util.py

 # MIT License
 #
-# Copyright (c) 2018-2024 Tskit Developers
+# Copyright (c) 2018-2025 Tskit Developers


Just to note here that the Copyright notice is meant to only be updated when the file itself is changed. So, the copyright years should (strictly speaking) be just those in which significant changes were made to the file.

Not that it matters in practise and there's no need to back out of these changes if you don't want to, but I just wanted to clarify my understanding of the copyright notices.

Oh yeah that makes sense, thank you for the note on that @jeromekelleher!

petrelharp · 2025-09-17T16:54:58Z

I've run into a legit bug: running tests/test_tree_stats.py::TestBranchAlleleFrequencySpectrum::test_wright_fisher[no_deep_history] we encounter a situation where an element of the AFS is -1. The good news is that this is only in the C version; the 'naive' and 'python' versions both have a 0 in that spot. @tforest, might you have a look at this?

tforest · 2025-09-17T18:02:48Z

I've run into a legit bug: running tests/test_tree_stats.py::TestBranchAlleleFrequencySpectrum::test_wright_fisher[no_deep_history] we encounter a situation where an element of the AFS is -1. The good news is that this is only in the C version; the 'naive' and 'python' versions both have a 0 in that spot. @tforest, might you have a look at this?

I'm having a look at this, thank you @petrelharp!

petrelharp · 2025-09-17T23:31:06Z

c/tests/testlib.c

+0.07┊ ┃ ┃ ┃ ┃ ┃ ┊ ┃ ┃  4  ┃ ┊ ┃ ┃  4  ┃ ┊         ┊
+    ┊ ┃ ┃ ┃ ┃ ┃ ┊ ┃ ┃ ┏┻┓ ┃ ┊ ┃ ┃ ┏┻┓ ┃ ┊         ┊
+0.00┊ 0 1 3 9 2 ┊ 0 1 2 3 9 ┊ 0 1 2 3 9 ┊ 0 1 2 3 ┊
+    0           2           7          10        100


I've added a non-sample leaf to this example, which we needed to get better coverage.

petrelharp · 2025-09-17T23:32:04Z

c/tskit/core.c

+            ret = "Time windows must start at zero and be strictly increasing. "
                  "(TSK_ERR_BAD_TIME_WINDOWS)";
            break;
+        case TSK_ERR_BAD_TIME_WINDOWS_END:


we needed to disambiguate the various notions of "bad" as different stats have different requirements now

petrelharp · 2025-09-17T23:34:34Z

python/tests/test_tree_stats.py

+
 np.random.seed(5)

+# Notes for refactoring:


These are notes about what happens in this file, in particular for a near-future me in refactoring the tests here.

petrelharp · 2025-09-17T23:35:49Z

python/tests/test_tree_stats.py

 ):
    mode = "branch"

+    def generate_params(self, ts):


Here is the main code testing the new time windows; this gets called on all the example tree sequences.

(Note that it also determines which combinations of windows and sample sets get tested.)

petrelharp · 2025-09-17T23:38:14Z

This is (if tests pass) ready to merge, we think? I've done some refactoring of the tests in test_tree_stats.py and have more in mind (before adding time windows to other stats functions), but will keep that for a different PR.

petrelharp · 2025-09-17T23:44:10Z

Ah one more thing - perhaps @tforest can benchmark this to see if we've affected performance for non-time-windowed AFS by a lot?

jeromekelleher

LGTM. I haven't gone through the logic, but the changes all seems solid modulo a few nitpicks.

jeromekelleher · 2025-09-18T10:18:37Z

c/tskit/trees.c

+        goto out;
+    }
+
+    if (windows[0] < 0.0) {


Could combine these two into one windows[0] <= 0.0 condition?

why so we could

or rather just !=

Hm, also I was just remembering that in fact times can be negative, so being < 0 is not necessarily an error

jeromekelleher · 2025-09-18T10:19:39Z

c/tskit/trees.c

 * SOFTWARE.
 */

+#include "tskit/core.h"


tskit includes should go after the std library ones

Anyway, isn't this already included already by trees.h ?

jeromekelleher · 2025-09-18T10:29:54Z

c/tskit/trees.c

+            goto out;
+        }
+        if (stat_site
+            && tsk_memcmp(time_windows, default_time_windows, 2 * sizeof(const double))


Took me a while to parse this - so it's saying that only the default time windows are supported for STAT_SITE, yeah? Might be easier to read and less fragile to say

// Site mode does not support time windows if (stat_site && !(time_windows[0] == 0 && isinf(time_windows[1]))

ah that would be MUCH better, thanks. I had thought there was some reason we couldn't do that, but it seems okay.

jeromekelleher · 2025-09-18T10:32:20Z

python/_tskitmodule.c

+    if (parse_windows(time_windows, &time_windows_array, &num_time_windows) != 0) {
+        goto out;
+    }
+    shape = PyMem_Malloc((num_sample_sets + 1 + 1) * sizeof(*shape));


The + 1 + 1s are confusing me - maybe just use + 2 and clarify why with a comment?

jeromekelleher · 2025-09-18T10:35:11Z

python/tskit/trees.py

                sample_sets = [sample_sets]
-                drop_dimension = True
-
+                if ll_method.__name__ != "allele_frequency_spectrum":


I'm assuming this is temporary?

this is not temporary (so maybe there's a better way to do it) - this is because of #3266 (comment)

Ah, I see we're doing this control flow based on introspecting the method name in a few places. This does smell quite badly I'm afraid

Can we refactor this by adding some parameters to __one_way_sample_set_stat which default to the standard behaviour, but we set to "False" when called from the top-level allele_frequency_spectrum function?

petrelharp · 2025-09-18T13:56:18Z

Thanks for the nitpicks, @jeromekelleher! I've done all of them except the "I assume this is temporary" one.

petrelharp · 2025-09-18T18:42:32Z

Hmph: this error now happens on Windows:

../tskit/trees.c:3756:70: error: conversion from 'double' to 'float' may change value [-Werror=float-conversion]
 3756 |         if (stat_site && !(time_windows[0] == 0 && isinf(time_windows[1]))) {
      |                                                          ~~~~~~~~~~~~^~~

Might you be able to see about how to deal with this, @tforest? (note: if debugging this is difficult, then making a temporary branch and PR with most of the testing code removed could be a good way to go)

jeromekelleher · 2025-09-18T18:53:42Z

Try 0.0 instead of 0

petrelharp · 2025-09-18T19:36:25Z

(whoops - missed the stuff in the comment above that clarified the problem is that the argument to isinf( ) is double but it's wanting a float)

tforest · 2025-09-18T19:37:44Z

Yes! I'm trying something @petrelharp :)

petrelharp · 2025-09-18T19:47:17Z

I think the problem is the argument: for Windows it looks like isinf requires a float, not a double?

tforest · 2025-09-18T19:58:59Z

That was it, yes. 👍

petrelharp · 2025-09-18T20:25:13Z

Nice. I'm squashing & rebasing.

jeromekelleher · 2025-09-19T08:08:04Z

I do think we should figure out a different way of dealing with dimensions etc than by introspecting the function pointer argument --- this is a really fragile and opaque way of doing control flow. Fine for this to be logged as an issue and done in a follow-up, but it it will bite eventually if we don't fix it.

petrelharp · 2025-09-19T16:40:56Z

Gotcha. I moved that AFS-specific pre-processing step to allele_frequency_spectrum, which seems the right way to do it.

and testing version of time windows for general stat test fixture cleanup remove not-run tests small issue with Tajimas D tests

petrelharp mentioned this pull request Sep 10, 2025

time windows in statistics #2948

Closed

petrelharp force-pushed the time_windows branch from 8a77af8 to e32062e Compare September 14, 2025 14:05

petrelharp mentioned this pull request Sep 16, 2025

time windows for general_stat (branch mode) #3271

Draft

petrelharp force-pushed the time_windows branch 2 times, most recently from f2f1f71 to fadb170 Compare September 16, 2025 20:05

petrelharp force-pushed the time_windows branch from fadb170 to 0af7d01 Compare September 16, 2025 20:28

jeromekelleher reviewed Sep 17, 2025

View reviewed changes

petrelharp force-pushed the time_windows branch from dfaa653 to d2e47c0 Compare September 17, 2025 23:30

petrelharp reviewed Sep 17, 2025

View reviewed changes

petrelharp force-pushed the time_windows branch from d2e47c0 to 051481e Compare September 17, 2025 23:33

petrelharp reviewed Sep 17, 2025

View reviewed changes

petrelharp force-pushed the time_windows branch from 051481e to 921c646 Compare September 17, 2025 23:40

jeromekelleher reviewed Sep 18, 2025

View reviewed changes

petrelharp force-pushed the time_windows branch from bb3eb02 to 11c8daf Compare September 18, 2025 20:25

Implementation of time windows for the AFS,

33f990d

and testing version of time windows for general stat test fixture cleanup remove not-run tests small issue with Tajimas D tests

petrelharp force-pushed the time_windows branch from 8a8c32e to 33f990d Compare September 19, 2025 17:05

jeromekelleher approved these changes Sep 20, 2025

View reviewed changes

jeromekelleher added this pull request to the merge queue Sep 20, 2025

Merged via the queue into tskit-dev:main with commit 44a23e2 Sep 20, 2025
18 checks passed

Time windows in AFS #3266

Time windows in AFS #3266

Uh oh!

Conversation

tforest commented Sep 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Uh oh!

codecov bot commented Sep 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

petrelharp commented Sep 14, 2025

Uh oh!

petrelharp commented Sep 16, 2025

Uh oh!

petrelharp commented Sep 16, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

petrelharp commented Sep 17, 2025

Uh oh!

tforest commented Sep 17, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

petrelharp commented Sep 17, 2025

Uh oh!

petrelharp commented Sep 17, 2025

Uh oh!

jeromekelleher left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

petrelharp commented Sep 18, 2025

Uh oh!

petrelharp commented Sep 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jeromekelleher commented Sep 18, 2025

Uh oh!

petrelharp commented Sep 18, 2025

Uh oh!

tforest commented Sep 18, 2025

Uh oh!

petrelharp commented Sep 18, 2025

Uh oh!

tforest commented Sep 9, 2025 •

edited

Loading

codecov bot commented Sep 9, 2025 •

edited

Loading

petrelharp commented Sep 18, 2025 •

edited

Loading