Initial steps for internal pedigree samples

[jeromekelleher]

Work towards #1855 - unfortunately it's trickier than I had hoped. The logic around lineages/individuals etc is really quite convoluted and it's hard to figure out where we need to inject the extra bits of ancestral material.

[jeromekelleher]

Just FYI @ThierryAM and @PandaGab these are my first steps here. Unfortunately it's more complicated than I had hoped because the logic is quite convoluted and I haven't been able to figure out quite where the required 3 lines of code should go. I'll pick away at it as I have a chance and keep you posted.

[codecov]

Codecov Report

Attention: Patch coverage is 76.66667% with 7 lines in your changes missing coverage. Please review.

Project coverage is 90.93%. Comparing base (b3b1be8) to head (be1bcc8).
Report is 3 commits behind head on main.

Files with missing lines	Patch %	Lines
lib/msprime.c	76.66%	4 Missing and 3 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2321      +/-   ##
==========================================
- Coverage   90.98%   90.93%   -0.05%     
==========================================
  Files          20       20              
  Lines       12021    12014       -7     
  Branches     2437     2438       +1     
==========================================
- Hits        10937    10925      -12     
- Misses        599      602       +3     
- Partials      485      487       +2     

Flag	Coverage Δ
C	90.93% <76.66%> (-0.05%)	⬇️
python	98.70% <ø> (ø)

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.

[jeromekelleher]

I think this is probably working now @ThierryAM and @PandaGab. I've done no real testing of the output here, so I think the next thing is to merge this much, and for you guys to hopefully put together some unit tests. It's probably worth while having a look through the PR here to see what bits of the code have been updated.

[gregorgorjanc]

I am monitoring this issue since we were bitten by it in the past.

I briefly looked into the changes using this example:

import msprime
import tskit
import io

ped_txt = """\
# id parent0 parent1 time is_sample
child1 mom1  dad1  0 1
mom1   gmom1 gdad1 1 1
dad1   gmom2 gdad2 1 0
child2 mom2  dad2  1 1
mom2   gmom1 gdad1 2 0
dad2   gmom2 gdad2 2 0
gmom1  .     .     3 0
gdad1  .     .     3 0
gmom2  .     .     3 0
gdad2  .     .     3 0
"""
# Note that mom1 is an internal sample
pedigree = msprime.parse_pedigree(io.StringIO(ped_txt), sequence_length=100)
ped_ts = msprime.sim_ancestry(
    initial_state=pedigree, model="fixed_pedigree", random_seed=41,
    recombination_rate=0.001)

Using the published package on conda I get the error (msprime._msprime.LibraryError: Samples that are internal nodes in the pedigree are not currently supported. Please comment on this GitHub issue if you would like to see this feature implemented: https://github.com/tskit-dev/msprime/issues/1855), while installing from this branch I get it to work!

Mom1 had id 1 and I can see her in the individuals table:

As well as in the nodes table marked as an internal sample with nodes 2 and 3:

Edges table has entries where nodes 2 and 3 are the child and an entry where node 3 is a parent:

Here is a plot of this tree sequence showing that node 2 is a sample (square) node, while 3 is sample (square) unary/ancestral node.

So, it works!

I would be happy to contribute to unit tests to learn how you do testing, but might need some guidance on how to do it - maybe best if I run/check someones tests first!

[jeromekelleher]

Superb, thanks @gregorgorjanc! OK, so since two of us think this probably works I'm going to merge it. I'll follow up on the original issue with details of what I think needs to be tested.