more efficient TMRCA computation?

[feilchenfeldt]

Hello,

I recently started to use tskit in python to iterate over tree sequences and interrogate the trees.
I wrote a function for getting the distance in the tree between two samples. Now I realised that tree.tmrca() is implemented, which basically gives the same thing.

However, I see that my function for getting the tmrca is consistently about 10 times faster than tree.tmrca() on my trees of 486 samples, while giving identical output.

I was trying to understand the implemented method, but the actuall code is burried in _ll objects which are obscure to me.

I post my code below and wanted to ask whether there is anything that speaks against implementing .tmrca that way, given that it seems to be much faster.

Analogously to TMRCA, one could also retrieve the id of the MRCA analogously (just returning i).

Thanks,
Hannes

def get_tmrca(tree, i, j):
    t0 = tree.get_time(i)
    t1 = tree.get_time(j)
    while i != j:
        ti = tree.get_time(i)
        tj = tree.get_time(j)
        if ti < tj:
            i = tree.get_parent(i)
        else:
            j = tree.get_parent(j)
    tmrca = tree.get_time(i) - (t0 - t1)/2
    return tmrca

[petrelharp]

Hm, interesting. That's a very nice and clear function! The .tmrca( ) method is implemented in C (here's the algorithm), which is why it's obscured. I'm surprised that it's faster for you? I did a quick test on a much bigger tree that shows the C method scaling much better:

t = msprime.simulate(1000000).first()
# almost instant
sum([t.tmrca(0,j) for j in range(1000000)])
# takes quite a long time
sum([get_tmrca(t,0,j) for j in range(1000000)])

This isn't a very careful timing test, of course.

In your case, was the built-in method actually taking a long time?

[jeromekelleher]

Hi @feilchenfeldt - just wondering if you're interested in following up on this? I've been looking into this and I can make the change pretty easily, but if someone else was interested in doing it I'd be happy to help them through the process.

[feilchenfeldt]

Hello @jeromekelleher,

We would like to contribute to tskit in the medium term.
There are some other features/improvements which we could not find and would like to see implemented.
(For example, a quick way of subsetting a tree to a subset of the samples and their histories and get the topolgy for these samples. We found that functionality for a tree sequence, but not for a single tree. Is the discussion forum the right place to ask about such things?)

For the current problem at hand, if you have looked into this I am sure you could implement it much faster and better, so please go ahead. Would that also improve the low efficiency with relate trees, or is this something the relate developpers should address? (Should we contact them?)

Best wishes,
Hannes

I actually think the algorithm above could be slightly more efficient like this:

def get_mrca(tree, i, j):
    t0 = tree.get_time(i)
    t1 = tree.get_time(j)
    ti = t0
    tj = t1
    while i != j:
        if ti < tj:
            i = tree.get_parent(i)
            ti = tree.get_time(i)
        else:
            j = tree.get_parent(j)
            tj = tree.get_time(j)
    return i

[jeromekelleher]

We would like to contribute to tskit in the medium term.

Excellent!

We found that functionality for a tree sequence, but not for a single tree. Is the discussion forum the right place to ask about such things?

The discussion forum is the perfect place for these discussions. Please go ahead and start as many threads as you like to explore the functionality you'd like to see.

Would that also improve the low efficiency with relate trees, or is this something the relate developpers should address? (Should we contact them?)

I'm not sure there much that can be done for RELATE trees, tbh - it's a fundamental property of the algorithm that node identities aren't shared between trees, and that's what tskit relies on efficiency. There might be some improvements that could be made, though, so there's no harm in letting Leo know the issues you're hitting.

For the current problem at hand, if you have looked into this I am sure you could implement it much faster and better, so please go ahead.

OK I'll implement the tmrca and ping you for review.

[feilchenfeldt]

Great, my PhD student Jaysmita @JSaha20 will be also interested in understanding the changes and reviewing.

[benjeffery]

Here's the PR that resulted from this for the benefit of future readers: #1313