Finding nodes with multiple parents

[hyanwong]

A quick 3-liner for finding those nodes in a tree sequence with multiple parents (see #2765 (comment))

uniq_child_parent = np.unique(np.column_stack((ts.edges_child, ts.edges_parent)), axis=0)
nd, count = np.unique(uniq_child_parent[:, 0], return_counts=True)
multiple_parents = nd[count > 1]
print(f"Nodes with multiple parents are {multiple_parents}")

Change the >1 to e.g. >2 to find nodes with 3 or more parents, etc. Note that in the record_full_arg option of msprime, the recombination nodes are the nodes immediately above any node with >1 parent. See tskit-dev/msprime#1942 (comment) for an explanation of this, which is only required because it allows calculation of likelihoods under the strict Hudson CwR (with max 1 breakpoint per chromosome)

[jeromekelleher]

This is a slightly more efficient/simpler way to do [EDIT] something slightly different it (as it avoids a sort):

node_is_child_count = np.bincount(ts.edges_child, minlength=ts.num_nodes)
multiple_parents = np.where(node_is_child_count > 1)[0]
print(multiple_parents)

In general np.bincount is very useful for working with these arrays. E.g., counting the number of mutations per node or site is easy:

sites_num_mutations = np.bincount(ts.mutations_site, minlength=ts.num_sites)
nodes_num_mutations = np.bincount(ts.mutations_node, minlength=ts.num_nodes)

[hyanwong]

Erm, won't that count cases where there are two edges with the same parent/child (e.g. non SMC or non-squashed) as having multiple parents?

[hyanwong]

You could do a bincount after the first sort, I guess? That would save one pass.

[jeromekelleher]

True, I guess that is something you have to worry about in standard coalescent sims.

[hyanwong]

I imagine you will certainly hit this in SLiM simulations. EDIT -you also see in msprime simulations that len(uniq_child_parent) < ts.num_edges, so as you say, this is also seen in msprime simulations