SAWilsonBalding does not keep root node number convention

[rbouckaert]

This may be a minor issue, but by convention, trees are assumed to have root node as highest numbered node, but the SAWilsonBalding operator can propose trees with root nodes of lower number since it calls Tree.setRootOnly() instead of Tree.setRoot() and the former does not automatically reorder nodes. The comment in the Tree says exactly that, but it gives no explanation why, so I am not sure why Tree.setRoot() is not called. Please let me know if you know why.

If an SAWilsonBalding proposal is followed by a NodeReight proposal in a StarBeast2 analysis, the latter does call Tree.setRoot() and changes node order so the root is highest numbered.

Perhaps other operators (SAExchange?) behave the same.

I bumped into this issue trying to find out why some caching issue on metadata of the species tree failed and caused a bit of trouble to identify the root, which I expected to be the highest numbered node.

[alexeid]

Wouldn’t it be most robust to either: (1) Make it impossible to call methods on Tree that would leave the root as not the highest number, or (2) Not assume anywhere that root has the highest node number. Otherwise, as it is, Tree objects can definitely have a root that does not have the highest node number, so we shouldn’t assume that is impossible anywhere.

…

On 20/09/2019, at 1:29 PM, Remco Bouckaert ***@***.***> wrote: This may be a minor issue, but by convention, trees are assumed to have root node as highest numbered node, but the SAWilsonBalding operator can propose trees with root nodes of lower number since it calls Tree.setRootOnly() instead of Tree.setRoot() and the former does not automatically reorder nodes. The comment in the Tree says exactly that, but it gives no explanation why, so I am not sure why Tree.setRoot() is not called. Please let me know if you know why. If an SAWilsonBalding proposal is followed by a NodeReight proposal in a StarBeast2 analysis, the latter does call Tree.setRoot() and changes node order so the root is highest numbered. Perhaps other operators (SAExchange?) behave the same. I bumped into this issue trying to find out why some caching issue on metadata of the species tree failed and caused a bit of trouble to identify the root, which I expected to be the highest numbered node. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#16?email_source=notifications&email_token=AAG5MSPXAMCT6UXAGPVUQRDQKQRN3A5CNFSM4IYRT7V2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HMR2YNQ>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAG5MSJVCFHNN6HYEULGNM3QKQRN3ANCNFSM4IYRT7VQ>.

[rbouckaert]

The reason it is convenient to have root as highest number is so that meta data attached to branches can be stored in parameters with indices 0...2N-2 for N taxa, like for branch rates in clock models. If the root cannot be relied upon to have number 2N-1, it will require more effort to map metadata to branches. Also, it means the relaxed clocks with sampled ancestors on a gene tree node may run into problems in the current implementation.

At the moment, Operators manipulate individual nodes, not the tree itself. Since nodes can have children added or removed (I think any reasonable API for nodes should allow this), that allows the tree to be left in an invalid state during the proposal() method. It will be hard to prevent operator implementation not to use such API.

What is not clear to me is why SAWilsonBalding calls Tree.setRootOnly() instead of Tree.setRoot(), which might solve this issue.

[alexeid]

The first problem predates BEAST2 :) I remember writing a constant-cost solution in BEAST1 that went something like: getTraitForNode(double[] values, Node node) { int ind = node.getNr(); int rootInd = node.getTree().getRoot().getNr(); if (ind < rootInd) return values[ind] else return values[ind-1]; } So the root could have any index you want and you could just skip it… Of course that means that which rate/trait matches which node skips around when the root moves, but generally speaking you have to do something to keep the traits/rates in sync when there are tree moves anyway if you want the correct greens ratio.

…

On 20/09/2019, at 2:27 PM, Remco Bouckaert ***@***.***> wrote: The reason it is convenient to have root as highest number is so that meta data attached to branches can be stored in parameters with indices 0...2N-2 for N taxa, like for branch rates in clock models. If the root cannot be relied upon to have number 2N-1, it will require more effort to map metadata to branches. Also, it means the relaxed clocks with sampled ancestors on a gene tree node may run into problems in the current implementation. At the moment, Operators manipulate individual nodes, not the tree itself. Since nodes can have children added or removed (I think any reasonable API for nodes should allow this), that allows the tree to be left in an invalid state during the proposal() method. It will be hard to prevent operator implementation not to use such API. What is not clear to me is why SAWilsonBalding calls Tree.setRootOnly() instead of Tree.setRoot(), which might solve this issue. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#16?email_source=notifications&email_token=AAG5MSOLLVCNTP5MRYB4IPTQKQYKNA5CNFSM4IYRT7V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7FLGIQ#issuecomment-533377826>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAG5MSLWYML32YHPLDKC4NTQKQYKNANCNFSM4IYRT7VQ>.

[gavryushkina]

What is not clear to me is why SAWilsonBalding calls Tree.setRootOnly() instead of Tree.setRoot(), which might solve this issue.

I don't remember why it is so but it was definitely for a reason. So changing to Tree.setRoot() may lead to some problems. Sorry I can't help currently with figuring this out.

[alexeid]

So what unit tests should we run to make sure that changing to Tree.setRoot() doesn’t break anything?

…

On 23/09/2019, at 9:29 AM, Alexandra Gavryushkina ***@***.***> wrote: What is not clear to me is why SAWilsonBalding calls Tree.setRootOnly() instead of Tree.setRoot(), which might solve this issue. I don't remember why it is so but it was definitely for a reason. So changing to Tree.setRoot() may lead to some problems. Sorry I can't help currently with figuring this out. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#16?email_source=notifications&email_token=AAG5MSK26UFAKRIXQIUTZSDQK7PT3A5CNFSM4IYRT7V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7JPVBQ#issuecomment-533920390>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAG5MSNAY5HXPETJNIYIZDLQK7PT3ANCNFSM4IYRT7VQ>.

[rbouckaert]

Removing the SAWilsonBalding operator from the XML solved my particular issue, so other SA operators appear to be no problem.

Tree.setRoot() calls getNodeCount(), which might have returned a different number of nodes in a past implementation where SA-trees were not binary trees with zero branch lengths for ancestral nodes. That might explain why Tree.setRootOnly() is called instead of Tree.setRoot(). Does that ring a bell?

Since SA-trees are now implemented as binary trees, this should not be a problem any more.

[alexeid]

This was also causing problems with Fabio’s recent testing of his new quantitative trait likelihoods on FBD trees...

…

On 25/09/2019, at 1:12 PM, Remco Bouckaert ***@***.***> wrote: Removing the SAWilsonBalding operator from the XML solved my particular issue, so other SA operators appear to be no problem. Tree.setRoot() calls getNodeCount(), which might have returned a different number of nodes in a past implementation where SA-trees were not binary trees with zero branch lengths for ancestral nodes. That might explain why Tree.setRootOnly() is called instead of Tree.setRoot(). Does that ring a bell? Since SA-trees are now implemented as binary trees, this should not be a problem any more. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#16?email_source=notifications&email_token=AAG5MSLQ4W7ZO2U5UA3IBKTQLK3JXA5CNFSM4IYRT7V2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD7QIMEA#issuecomment-534808080>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAG5MSMZBQHDZ3EVV3KHS23QLK3JXANCNFSM4IYRT7VQ>.

[Anaphory]

I ran into this issue these days with the Optimized Relaxed Clock. It took me a while to figure out why I had certain issues.

I have documented my attempts at fixing or at least mitigating the issue in CompEvol/beast2#984, but I got the impression that the core of the issue on the SA side is this order of operations:

sampled-ancestors/src/beast/evolution/operators/SAWilsonBalding.java

Lines 156 to 164 in ff14640

	if (jP != null) {
	jP.removeChild(j); // remove <jP, j>
	jP.addChild(iP); // add <jP, iP>
	jP.makeDirty(Tree.IS_FILTHY);
	} else {
	iP.setParent(null); // completely remove <PiP, iP>
	tree.setRootOnly(iP);
	}
	iP.addChild(j);

It could maybe be solved by pulling the

         iP.addChild(j);

into the two branches of the if, to be able to put it before the tree.setRootOnly(iP); (or then instead tree.setRoot(iP)) in the else block.