Feat:/simulate after local mrca #2396

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.23%. Comparing base (1ab3389) to head (2ef10be).

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.23%   91.23%           
=======================================
  Files          20       20           
  Lines       11831    11847   +16     
  Branches     2296     2300    +4     
=======================================
+ Hits        10794    10809   +15     
  Misses        568      568           
- Partials      469      470    +1

Flag	Coverage Δ
C	`91.23% <91.17%> (+<0.01%)`	⬆️
python	`98.68% <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.

lib/msprime.c

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

lib/msprime.c

        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

msprime/ancestry.py

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

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

lib/tests/test_ancestry.c

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?)

docs/ancestry.md

lib/tests/test_ancestry.c

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

msprime/ancestry.py


+    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

msprime/ancestry.py

tests/test_ancestry.py

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

lib/tests/test_ancestry.c

+        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?

tests/test_ancestry.py

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

algorithms.py

                        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

algorithms.py

            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

msprime/ancestry.py

    """


+class PotentialInfiniteSimulationhWarning(UserWarning):


Suggested change

class PotentialInfiniteSimulationhWarning(UserWarning):

class PotentialInfiniteSimulationWarning(UserWarning):

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

lib/msprime.c

 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

lib/msprime.c

                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

lib/msprime.c

        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

lib/msprime.c

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


Same style points as above

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

lib/msprime.c

    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!

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

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.

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

lib/tests/test_ancestry.c Outdated Show resolved Hide resolved

lib/tests/test_ancestry.c Outdated Show resolved Hide resolved

jeromekelleher reviewed Aug 28, 2025

View reviewed changes

hossam26644 mentioned this pull request Aug 29, 2025

Simulating above the root breaks smc models #2405

Open

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

tests/test_ancestry.py Outdated Show resolved Hide resolved

jeromekelleher reviewed Sep 15, 2025

View reviewed changes

hossam26644 force-pushed the feat-simulate-after-local-mrca branch 2 times, most recently from f083429 to 2a250e9 Compare October 1, 2025 11:30

simulations above root implementation

ec37d9f

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

Change stopping condition for stop_at_local_mrca=False

8ba4e04

jeromekelleher force-pushed the feat-simulate-after-local-mrca branch from 239b8bb to be1b404 Compare December 3, 2025 15:17

Add tests for fixed pedigree and clarify semantics; document

2ef10be

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

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

jeromekelleher mentioned this pull request Dec 3, 2025

Rerun perf benchmarks #2430

Open

	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):

Feat:/simulate after local mrca #2396

Are you sure you want to change the base?

Feat:/simulate after local mrca #2396

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!

codecov bot commented Aug 5, 2025 •

edited

Loading

hossam26644 commented Sep 7, 2025 •

edited

Loading

hossam26644 commented Dec 2, 2025 •

edited

Loading