Tskit: kc_distance doesn't ignore unary nodes

[hyanwong]

We should presumably have the same KC distance regardless of unary nodes. But the example below shows that it doesn't work:

# use record_full_arg to leave in some unary nodes
ts_unary = msprime.simulate(10, random_seed=1, recombination_rate=10, record_full_arg=True)

tree02 = ts_unary.at(0.2)
tree08 = ts_unary.at(0.8)
kc_unary = tree02.kc_distance(tree08)

ts_nounary = ts_unary.simplify()

tree02 = ts_nounary.at(0.2)
tree08 = ts_nounary.at(0.8)
kc_nounary = tree02.kc_distance(tree08)

assert kc_unary == kc_nounary  # currently fails

ping @daniel-goldstein as he's working on KC distance calcs at the moment.

[daniel-goldstein]

I haven't worked with simplify or unary nodes much, so I could easily be wrong here, but would simplifying a unary node change the time from that node to its immediate parent? Aside from pair-wise comparisons, kc also considers the branch length of each leaf.

Although, by default kc_distance should ignore time and only use topology so this case should work...

[hyanwong]

Yes, this is with topology-only metrics, so should work.

[hyanwong]

When you simplify() like this, you remove the unary nodes. They should probably be being ignored in the kc distance. Indeed, we might want to simplify before running the calculations, within the kc_distance() code. Probably worth chatting to @jeromekelleher about it.

[jeromekelleher]

This is interesting. I'm not immediately convinced that unary nodes should be ignored (but I agree that seems like the most likely thing to do). What are the actual values output in your sims above @hyanwong? We should boil this down to the simplest possible example and see what we're actually doing with the presence of unary nodes and then see whether this is contrary to the definition of KC distance.

It might end up that we need a parameter ignore_unary=True, if it turns out that we're currently doing something consistent in the presence of unary nodes.

[hyanwong]

What are the actual values output in your sims above @hyanwong?

They are massively different, by an order of magnitude:

>>> kc_unary
97.8417088975862
>>> kc_nounary
9.643650760992955

We should boil this down to the simplest possible example and see what we're actually doing with the presence of unary nodes and then see whether this is contrary to the definition of KC distance.

Agreed

[hyanwong]

NB: I noticed this when comparing the old tsinfer (which simplified with keep_unary=False) to the new one (which has keep_unary=True). It made it appear that the new versions were way worse in terms of topology reconstruction.

[jeromekelleher]

Huh, sounds like a bug in KC if they're way different. Good catch.

[hyanwong]

I'm inclined to leave this to @daniel-goldstein to fix, if he's happy with that, as he is working on the code at the moment. Is that OK Daniel? Happy to take it on if not, though.

[daniel-goldstein]

I’d be happy to take it!

[daniel-goldstein]

So the removal of unary nodes does change the topology as far as KC is concerned. The value m[u, v] is the depth of the mrca of u and v, which will differ before/after simplifying if there are unary nodes above the mrca.

An easy "fix" would be not to count unary nodes toward the depth count, but this does feel like a deliberate decision that could be a flag, like @jeromekelleher suggested. I tested @hyanwong's example with ignoring unary nodes in the topology and the distances were then equivalent. This also explains why the distance was so high before simplifying.

The KC paper deals only with binary trees, in which case it might make sense to default to ignoring unary nodes.

[jeromekelleher]

Right, because depth is measured in node hops.

Simplest thing to do would be to raise a ValueError for now if we detect a unary node. Easy enough to simplify the tree sequence if you want to do comparisons. Any objections @hyanwong?

[hyanwong]

Yep. Or simplify the TS under the hood and document this. Either is fine by me.

[jeromekelleher]

Yep. Or simplify the TS under the hood and document this. Either is fine by me.

Don't really want to simplify as this is a method of the Tree. Gets complicated.

[daniel-goldstein]

Ok, I'll just throw an error for unary nodes.

[hyanwong]

Of course: a good point about it being a tree method. Not sure what the semantics should be when we make an equivalent TS method, though.