Statistical testing with keep_unary

[jeromekelleher]

We should run some statistical tests on the trees output by keep_unary, just to make sure that nothing weird has happened. It's probably sufficient to run some simple downstream stats that shouldn't be affected by these extra nodes.

[jeromekelleher]

Should also do this with some stats on the mutations so that we know sim_mutations is doing the right thing on these topologies (see also #2137)

[sgravel]

My expectation is that the recording of unary nodes should not affect the pseudorandom generation, so we should be able to see that the statistics are exactly equal before and after recording. Is that what you had in mind?
If so we are preparing a few tests of that sort.

[General-Solution]

I was looking at the branch-based diversity statistic between simulations with and without unary nodes in the Hudson model, and found that they are sometimes but not always equal. After investigating more it seems that most tree sequences have the same structure with and without unary nodes but some are very different.

[General-Solution]

For example pairs of tree sequences will turn out like this:

(2 samples, 0.1 replication rate, seed 14)

[General-Solution]

But then some will turn out like this:

(2 samples, 0.1 replication rate, seed 18)

[sgravel]

So I guess the question is whether we expect the recording of unary nodes to change the simulated ts under a given seed. If it is purely a recording process, it should probably not affect the trees. On the other hand, it is possible that the recording of unary nodes can change something (e.g., the ordering of lineages) which ultimately can affect the tree sequence (but not its distribution). It would be nice (but not critical) for the node recording to not affect the simulation outcome.
In this case it seemed interesting behaviour because simulations were identical most of the time, with maybe 10% of simulations showing differences.

[GertjanBisschop]

Thanks @General-Solution for looking into this! I agree @sgravel, ideally, recording nodes should not affect the simulation outcome. This is also what I would have expected the behaviour to be.
Could it be that defragmenting the segment chain (which will be different with or without record_unary) leads to a different structure of the recombination_mass_index @jeromekelleher ?

[sgravel]

That sounds plausible to me. It does happen fairly frequently that adjacent trees have different unary nodes.
, presumably due to a recombination and coalescence within the same lineage. In principle the ancestral segment above the coalescence could be identical and span the two trees, but it is possible that recombination treats the two incoming edges as separate segments, whereas they are a single segment in the non-recording.

[jeromekelleher]

Sorry about the delay in getting back, was at a conference and had an admin mountain to climb...

So I guess the question is whether we expect the recording of unary nodes to change the simulated ts under a given seed.

I don't think we can make any guarantees about equivalence under a given seed. I agree it would be nice if we could, but even the tiniest change in the state can alter the simulation trajectory and past experience has taught me that it's extremely difficult keep the internal states identical when doing slightly different things.

I don't see immediately why the simulations would diverge here, but it would have been surprising to me had the outputs turned out to be identical in all cases (not the other way around!).

[jeromekelleher]

Regarding statistical tests, I think we can probably just redo a few of the existing tests (with subclassing?) with this keep_unary flag set and we should see identical results. For example, any mutation based tests should be unaffected.

[sgravel]

Regarding the seed: Ok with me, it's probably not worth refactoring the code to maintain the seed behaviour. Just wanted to make sure.
Regarding the statistical tests: sounds good. To be clear, when talking about statistical tests, you are referring to adding tests in verification.py, rather than in unit tests, where we are returning plots and eyeballing rather than coming up with formal unit test criteria?

[General-Solution]

Regarding statistical tests, I think we can probably just redo a few of the existing tests (with subclassing?) with this keep_unary flag set and we should see identical results. For example, any mutation based tests should be unaffected.

By mutation based tests, is that in the sense of counting diversity as a site statistic instead of a branch statistic? Because that seems to also be changed by including store-unary.

[sgravel]

@General-Solution
My understanding is that the statistical tests (that cannot be easily factored in a deterministic unit test) are collected in verification.py. You'll find e.g. "test_analytical_pi" that performs a simulation and computes pairwise_diversity and generates some plots that can be inspected to ascertain that things behave as expected. I believe that @jeromekelleher means that we could replace the simulations in these tests with simulations that record the unary node. In this case, the modified tests should give essentially identical results up to fluctuations.
(This is by contrast to tests based on the length of edges or similar metrics could give different results.)

[jeromekelleher]

Regarding the statistical tests: sounds good. To be clear, when talking about statistical tests, you are referring to adding tests in verification.py, rather than in unit tests, where we are returning plots and eyeballing rather than coming up with formal unit test criteria?

Yes, exactly.

[GertjanBisschop]

Solved for StandardCoalescent, SMC, SMCprime and multimergers in #2156.

[GertjanBisschop]

Solved for DTWF and FIXED_PEDIGREES in #2225.