Feat:/simulate after local mrca by hossam26644 · Pull Request #2396 · tskit-dev/msprime

hossam26644 · 2025-08-05T09:31:25Z

Here you go: @jeromekelleher
Still missing some tests, and I want to change the flag name of the python algorithm to match the msprime implementation.
As a reminder, this implementation consumes more memory than the when stopping at the youngest root. It could be the node table getting much larger, and could be a memory leak.

codecov · 2025-08-05T09:40:58Z

Codecov Report

❌ Patch coverage is 91.17647% with 3 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.22%. Comparing base (279f480) to head (19dafa7).
⚠️ Report is 4 commits behind head on main.

Files with missing lines	Patch %	Lines
lib/msprime.c	90.00%	1 Missing and 2 partials ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2396   +/-   ##
=======================================
  Coverage   91.21%   91.22%           
=======================================
  Files          20       20           
  Lines       11831    11847   +16     
  Branches     2296     2300    +4     
=======================================
+ Hits        10792    10807   +15     
  Misses        569      569           
- Partials      470      471    +1

Flag	Coverage Δ
C	`91.22% <91.17%> (+<0.01%)`	⬆️
python	`98.63% <100.00%> (+<0.01%)`	⬆️

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

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

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

jeromekelleher · 2025-08-05T10:10:06Z

Implementation looks good, it seems unlikely that it's a memory leak to me. Let's look at the size of the output tables I think. Probably this will be a function of your value of stop_time.

You need to run the C code through clang-format as a minimum to get CI passing.

jeromekelleher · 2025-08-05T10:11:23Z

Note specific requirements for clang-format: https://tskit.dev/tskit/docs/stable/development.html#sec-development-c-requirements

hossam26644 · 2025-08-08T14:40:36Z

It is ready for review now.

Please check:

the tests I added and where I put them.
if we need extra tests to for the implementation (I can not think of another tests for the implementation itself)

jeromekelleher

Looks good!

Some minor comment above.

Also needed:

Some C unit tests. We should be able to add a test or two to test_ancestry.c to cover the functionality. Nothing elaborate, just making sure the code is covered by valgrind.
We also need to cover the other models. This will require updating msp_merge_ancestors in a similar way. Best to start with algorithms.py on this.

jeromekelleher · 2025-08-11T13:36:38Z

        for (avl_node = self->non_empty_populations.head; avl_node != NULL;
             avl_node = avl_node->next) {
-            pop_id = (tsk_id_t)(intptr_t) avl_node->item;
+            pop_id = (tsk_id_t) (intptr_t) avl_node->item;


Seems to be a bunch of formatting noise here - which version of clang-format are you using?

version 18.. I will remove this noise

jeromekelleher · 2025-08-27T15:36:29Z

I'm going to mute this thread @hossam26644 - please ping me when you'd like a review

hossam26644 · 2025-08-28T14:08:37Z

@jeromekelleher check now please, and please check the implementation of merge_ancestors. It works and passes the tests that I put.. but I am not sure if I am missing something

jeromekelleher

Generally looking good. Needs more models though, in particular the pedigree model. You need to construct a pedigree in which full coalescence will happen and verify that you still go through the rest of the nodes (you'll need store_unary as well I guess, but that's true anyway?)

jeromekelleher · 2025-08-28T14:55:59Z


+    if stop_at_local_mrca is None:
+        stop_at_local_mrca = True
+    if not stop_at_local_mrca:  # when set to False


Comment is redundant

jeromekelleher · 2025-09-02T09:02:05Z

+        CU_ASSERT_EQUAL(ret, 0);
+
+        ret = msp_run(&msp, t, ULONG_MAX);
+        CU_ASSERT_EQUAL(ret, 2); // exit on max time


Use the symbolic constant, no need for comment then

hossam26644 · 2025-09-07T13:22:00Z

@jeromekelleher
Ready for another review.

I disabled stop_at_local_mrca with fixed pedigree simulations since it make not sense (if I understand correctly)

One more thing:

I added a condition to stop simulations when there are no events to happen and stop_at_local_mrca is False, i.e., no coalescence, GC, recombination, or migration events; it used to give an error before, since simulations should stop one step before that. This error was raised for smc models only, while for Hudson's model there are always recombination then coalescence events (unless the seq length is 1 base, then it behaves like the smc models). These recombination and coalescence events happen even when the recombination rate is zero (for Hudson).

Now, if we can put a condition when simulations get to this infinite cycle of recombination and coalescence, we won't have to specify an end time; simulations can run until more time means nothing (no unary nodes no recombination no coalescence). I don't understand fully why we see this when recombination rate == zero, so maybe I am missing important details; but if think it will be useful, I can remove the requirement for an end time for simulations above the root.

jeromekelleher

The infinite waiting time problem is in general hard (there's some references to it in the documentation if you look for it). It would be great if we could keep simulating until there are no events for the SMC, but we'll have to be careful not to introduce any unexpected behaviours otherwise.

Re the pedigree model - I don't see why it's not relevant? We can construct a pedigree in which coalescence happens very quickly, but we still want to track how the ancestral material passes through the ancestors. What happens at the moment when partial or full coalescence happens within the pedigree?

hossam26644 · 2025-09-11T14:53:11Z

@jeromekelleher we can go for another round of review:

Simulate above root is implemented with no end_time set. Please check the implementation especially for dtwf (I had to do some gymnastics to make it work)
Test cases for with and without end_time are implemented
Test case for the pedigree model built

jeromekelleher

I'm not following some of the logic, but I think I just need to check it out an play with it myself to get my head around it. I think we should just clean this much up, squash the commits and merge. We can register a few follow-ups issues to close before release.

jeromekelleher · 2025-09-15T12:27:39Z

                        mig_dest = k
            min_time = min(t_re, t_ca, t_gcin, t_gc_left, t_mig)
-            assert min_time != infinity
+            assert (min_time != infinity) or (


Comparing with the False singleton is an antipattern. Use not self.stop_at_local_mrca

jeromekelleher · 2025-09-15T12:29:13Z

            else:
                self.t += min_time
-                if min_time == t_re:
+                if (min_time == infinity) and self.stop_at_local_mrca is False:


Suggested change

if (min_time == infinity) and self.stop_at_local_mrca is False:

if (not self.stop_at_local_mrca) and (min_time == infinity):

Put the constant first so it gets short-circuited. (Not that it matters here, but just as a general principle)

jeromekelleher · 2025-09-15T12:38:58Z

    """


+class PotentialInfiniteSimulationhWarning(UserWarning):


Suggested change

class PotentialInfiniteSimulationhWarning(UserWarning):

class PotentialInfiniteSimulationWarning(UserWarning):

jeromekelleher · 2025-09-15T12:47:03Z

 out:
    msp_safe_free(parents);
    msp_safe_free(segment_mem);
+    if (coalescence_occurred == false && recombination_occurred == false && ret == 0) {


These out blocks are for error condition and cleanup handling only - move this block up to before the goto. Also you don't need to test for ret == 0 then.

Also, change to if ((!coalescence_occurred)) && (!recombination_occurred))

jeromekelleher · 2025-09-15T12:47:35Z

                goto out;
            }
            ret = msp_apply_demographic_events(self, self->next_demographic_event->time);
+            migration_occurred


Same point as before about comments on the same line

jeromekelleher · 2025-09-15T12:48:19Z

        self->time = cur_time;
        ret = msp_dtwf_generation(self);
-        if (ret != 0) {
+        if (ret != 0 && ret != MSP_DTWF_DID_NOTHING) {


Better to change the condition to if (ret < 0) .

jeromekelleher · 2025-09-15T12:48:58Z

            self->next_sampling_event++;
        }
+
+        if (ret == MSP_DTWF_DID_NOTHING && migration_occurred == false


Same style points as above

jeromekelleher · 2025-09-15T12:49:24Z

    msp_safe_free(node_trees);
    msp_safe_free(n);
    msp_safe_free(mig_tmp);
+    if (ret == MSP_DTWF_DID_NOTHING) {


Ditto on putting non cleanup logic in out block

hossam26644 · 2025-12-02T10:11:44Z

up @jeromekelleher

jeromekelleher · 2025-12-02T15:48:57Z

@hossam26644 just FYI - I started looking into this in detail and it seemed easier to me to fix the stopping condition than to explain why we needed the max_time thing. I've fixed it for Hudson and am filling out the rest of the details now.

jeromekelleher · 2025-12-02T17:10:42Z

Hmm - looks like there's a basic flaw in the logic for the multi-merger case, which wasn't really being tested. I'll need to spend more time on this to get to the bottom of it.

hossam26644 · 2025-12-03T09:27:30Z

I can invest some time into it. Actually, I would like to. But I am prioritising the paper draft now. If this is not urgent, maybe it can sit for a while until after the new year.

jeromekelleher · 2025-12-03T09:37:04Z

Leave it with me - I'm going to see how far I can get with it today. There will be some follow-up work to do all right though, don't worry!

jeromekelleher · 2025-12-03T15:26:48Z

OK, I think this is basically done. Can you review please @hossam26644? I'm going to log some follow-up issues.

jeromekelleher · 2025-12-09T09:39:25Z

Something happened to the branch here @hossam26644 - it needs to be rebased

hossam26644

I think it is neat

hossam26644 · 2025-12-09T11:06:23Z

        self.verify_segments()
-        if self.model not in ["fixed_pedigree", "dtwf"]:
-            # The fixed_pedigree model doesn't maintain a bunch of stuff.
-            # It would probably be simpler if it did.


why wasdtwf excluded from verify overlaps? why is it included in the new version?

Good question - it was incorrectly omitted before, and it means that we didn't pick up the basic logic error in the multi-merger case (subtracting min_overlap, instead of len(h)). We also weren't testing this code path at all.

hossam26644 · 2025-12-09T11:08:51Z

 ```{code-cell}
 ts = msprime.sim_ancestry(2, sequence_length=10, recombination_rate=0.1, random_seed=42)
-SVG(ts.draw_svg(y_axis=True, y_ticks=[0, 1, 2, 3, 4, 5]))
+SVG(ts.draw_svg(y_axis=True, y_ticks=[0, 1, 2, 3, 4]))


should this be in a separate commit or is fine here?

Yes, good point, this should be a separate commit.

hossam26644 · 2025-12-09T11:12:18Z

-    with flags should be more intuitive now that we are relying on the 
+    These flags have been deprecated in favour of using
+    :class:`NodeType<NodeType>`. This is not a breaking change. The
+    constant associated with each flag remains the same. However, working


there is a bunch of formatting changes in this file and nothing more, do you want to include them in the PR?

Ah yes, good point. It's just whitespace trimming which my editor does automatically and should probably be enforced, so I think it's fine to include.

jeromekelleher

OK, thanks for the review. I'll push an update and merge.

jeromekelleher · 2025-12-09T12:00:16Z

-    with flags should be more intuitive now that we are relying on the 
+    These flags have been deprecated in favour of using
+    :class:`NodeType<NodeType>`. This is not a breaking change. The
+    constant associated with each flag remains the same. However, working


Ah yes, good point. It's just whitespace trimming which my editor does automatically and should probably be enforced, so I think it's fine to include.

jeromekelleher · 2025-12-09T12:03:22Z

        self.verify_segments()
-        if self.model not in ["fixed_pedigree", "dtwf"]:
-            # The fixed_pedigree model doesn't maintain a bunch of stuff.
-            # It would probably be simpler if it did.


Good question - it was incorrectly omitted before, and it means that we didn't pick up the basic logic error in the multi-merger case (subtracting min_overlap, instead of len(h)). We also weren't testing this code path at all.

jeromekelleher · 2025-12-09T12:04:40Z

 ```{code-cell}
 ts = msprime.sim_ancestry(2, sequence_length=10, recombination_rate=0.1, random_seed=42)
-SVG(ts.draw_svg(y_axis=True, y_ticks=[0, 1, 2, 3, 4, 5]))
+SVG(ts.draw_svg(y_axis=True, y_ticks=[0, 1, 2, 3, 4]))


Yes, good point, this should be a separate commit.

jeromekelleher · 2025-12-09T12:11:32Z

Going to let #2432 go through first, then resolve merge conflicts

Closes tskit-dev#2157 Closes tskit-dev#1898

hossam26644 marked this pull request as ready for review August 8, 2025 14:38

jeromekelleher reviewed Aug 11, 2025

View reviewed changes

hossam26644 commented Aug 28, 2025

View reviewed changes

Comment thread lib/tests/test_ancestry.c Outdated

Comment thread lib/tests/test_ancestry.c Outdated

jeromekelleher reviewed Aug 28, 2025

View reviewed changes

hossam26644 mentioned this pull request Aug 29, 2025

Simulating above the root breaks smc models #2405

Closed

jeromekelleher reviewed Sep 2, 2025

View reviewed changes

jeromekelleher mentioned this pull request Sep 2, 2025

First pass on stop_at_local_mrca #2392

Closed

hossam26644 force-pushed the feat-simulate-after-local-mrca branch from 26f42eb to 640ff6e Compare September 2, 2025 14:27

hossam26644 force-pushed the feat-simulate-after-local-mrca branch from 84de697 to f9af0ec Compare September 7, 2025 13:22

jeromekelleher reviewed Sep 8, 2025

View reviewed changes

Comment thread tests/test_ancestry.py Outdated

jeromekelleher reviewed Sep 15, 2025

View reviewed changes

hossam26644 force-pushed the feat-simulate-after-local-mrca branch 3 times, most recently from 2a250e9 to ec37d9f Compare October 4, 2025 06:37

jeromekelleher force-pushed the feat-simulate-after-local-mrca branch 2 times, most recently from be1b404 to 2ef10be Compare December 3, 2025 15:26

jeromekelleher mentioned this pull request Dec 3, 2025

Rerun ASV benchmarks pre 1.4 release #2430

Closed

hossam26644 force-pushed the feat-simulate-after-local-mrca branch from c4f251f to f1f7635 Compare December 9, 2025 09:42

hossam26644 commented Dec 9, 2025

View reviewed changes

jeromekelleher reviewed Dec 9, 2025

View reviewed changes

jeromekelleher force-pushed the feat-simulate-after-local-mrca branch from f1f7635 to 88fc6b9 Compare December 9, 2025 12:09

jeromekelleher approved these changes Dec 9, 2025

View reviewed changes

jeromekelleher enabled auto-merge December 9, 2025 12:10

jeromekelleher disabled auto-merge December 9, 2025 12:11

jeromekelleher added 4 commits December 9, 2025 12:53

simulations above root implementation

fa8b7e2

Change stopping condition for stop_at_local_mrca=False

96df64f

Add tests for fixed pedigree and clarify semantics; document

7bfefb2

Closes tskit-dev#2157 Closes tskit-dev#1898

Fix trailing tickmark in docs

19dafa7

jeromekelleher force-pushed the feat-simulate-after-local-mrca branch from 88fc6b9 to 19dafa7 Compare December 9, 2025 12:55

jeromekelleher enabled auto-merge December 9, 2025 12:56

jeromekelleher added this pull request to the merge queue Dec 9, 2025

Merged via the queue into tskit-dev:main with commit 1b2f7fd Dec 9, 2025
13 checks passed

jeromekelleher mentioned this pull request Dec 15, 2025

Add option to simulate "above" local roots #2157

Closed

	if (min_time == infinity) and self.stop_at_local_mrca is False:
	if (not self.stop_at_local_mrca) and (min_time == infinity):

	class PotentialInfiniteSimulationhWarning(UserWarning):
	class PotentialInfiniteSimulationWarning(UserWarning):

Conversation

hossam26644 commented Aug 5, 2025

Uh oh!

codecov Bot commented Aug 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

jeromekelleher commented Aug 5, 2025

Uh oh!

jeromekelleher commented Aug 5, 2025

Uh oh!

hossam26644 commented Aug 8, 2025

Uh oh!

jeromekelleher left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

jeromekelleher commented Aug 27, 2025

Uh oh!

Uh oh!

Uh oh!

hossam26644 commented Aug 28, 2025

Uh oh!

jeromekelleher left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

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

Uh oh!

Uh oh!

jeromekelleher left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

hossam26644 commented Sep 11, 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!

hossam26644 commented Dec 2, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jeromekelleher commented Dec 2, 2025

Uh oh!

jeromekelleher commented Dec 2, 2025

Uh oh!

hossam26644 commented Dec 3, 2025

codecov Bot commented Aug 5, 2025 •

edited

Loading

hossam26644 commented Sep 7, 2025 •

edited

Loading

hossam26644 commented Dec 2, 2025 •

edited

Loading