Shortest subtree that contains half the nodes of a population

[stsmall]

Hi All,
I was reading Hejase et al. 2020, and they use a statistic, TMRCAH, which was defined as the shortest subtree that contains half the nodes of a population. What would be an efficient way to calculate the TMRCAH statistic? Would this be a case where I would need to trace all the paths from the leaves to the root? Maybe combining that with simplify?
thanks,
@stsmall

[petrelharp]

Trees already keep track of the number of samples below each node (and efficiently update this), so you could do something like:

ts = msprime.sim_ancestry(2, recombination_rate=0.1, sequence_length=40)
n0 = ts.num_samples / 2
for t in ts.trees():
    tmrcah = min([t.time(n) for n in t.nodes() if t.num_samples(n) > n0])
    print(f"TMRCAH on {t.interval} is {tmrcah}")

or maybe more efficiently

for t in ts.trees():
    tmrcah = np.inf
    for n in t.nodes(order='timeasc'):
        if t.num_samples(n) > n0:
          tmrcah = t.time(n)
          break
    print(f"TMRCAH on {t.interval} is {tmrcah}")

[stsmall]

cool, I didnt know about the timeasc order. I see it in the docs now. Thanks!