How to generate offspring from a given set of individuals?

[petrelharp]

@gregorgorjanc asks: "Say I have a tree sequence with nodes associated with individuals. Then I would like to simulate their possible progeny, so a forward step." How to do this?

I think what we'd like to do here is: given a tree sequence and a list of pairs of individuals, add on to it one progeny per pair produced by randomly sampling recombination breakpoints from a rate map. Perhaps the resulting tree sequence should have the new progeny marked as samples (as well as whatever already was samples) and it should return the list of those new sample IDs?

[petrelharp]

Here's a way to do it (note: should be tested more):

import numpy as np
import tskit, msprime
rng = np.random.default_rng()

ts = msprime.sim_ancestry(4, sequence_length=10)
parents = [(0,1), (0, 1), (2, 2)]
recmap = msprime.RecombinationMap(positions=[0, 5, 10], rates=[1.0, 0.2, 0.0])
child_time = -1

def sample_breakpoints(r):
    L = r.get_total_recombination_rate()
    n = rng.poisson(lam=L)
    u = np.sort(rng.uniform(low=0.0, high=L, size=n))
    return np.array([r.genetic_to_physical(x) for x in u])

t = ts.dump_tables()
for pp in parents:
    for p in pp:
        child = t.nodes.add_row(time=child_time, flags=tskit.NODE_IS_SAMPLE)
        nodes = ts.individual(p).nodes
        assert len(nodes) == 2
        j = rng.choice([0, 1])
        left = 0.0
        for right in sample_breakpoints(recmap):
            t.edges.add_row(parent=nodes[j], child=child, left=left, right=right)
            left = right
            j = (j + 1) % 2
        t.edges.add_row(parent=nodes[j], child=child, left=left, right=t.sequence_length)

t.sort()
new_ts = t.tree_sequence()

I'm using the RecombinationMap instead of the RateMap since it has the genetic_to_physical method, and adding the new offspring all at child_time; perhaps the input could also specify what time they should be at?