CompEvol / beast2

Bayesian Evolutionary Analysis by Sampling Trees
www.beast2.org
GNU Lesser General Public License v2.1
240 stars 84 forks source link

Add case to Tree.adjustTreeNodeHeights method to accommodate zero-length branches #1161

Open jjmccollum opened 3 months ago

jjmccollum commented 3 months ago

Methods like Node.isDirectAncestor and Node.isFake assume a tree in which branches of length 0 can occur, and this is an important case for me to include in a BEAST package I'm writing. In theory, we can allow for zero-length branches by setting the adjustTreeNodeHeights attribute of the tree element in our input to false. But in practice, I find that if I do this, many of the branches in my starting tree have negative lengths, which I still don't want. That said, I'd like to have a version of the Tree.adjustTreeNodeHeights method that cleans up negative branch lengths, but leaves branch lengths of 0 intact.

Currently, this method is defined as follows:

/**
     * Ensure no negative branch lengths exist in tree.  This can occur if
     * leaf heights given as a trait are incompatible with the existing tree.
     */
    final static double EPSILON = 0.0000001;

    protected void adjustTreeNodeHeights(final Node node) {
        if (!node.isLeaf() && adjustTreeNodeHeightsInput.get()) {
            for (final Node child : node.getChildren()) {
                adjustTreeNodeHeights(child);
            }
            for (final Node child : node.getChildren()) {
                final double minHeight = child.getHeight() + EPSILON;
                if (node.getHeight() < minHeight) {
                    node.setHeight(minHeight);
                }
            }
        }
    }

Notably, the description of the method already suggests what I have in mind: "no negative branch lengths". The problem is that the EPSILON constant is strictly positive, which makes branch lengths of 0 impossible.

Would the following modification be acceptable?

/**
     * Ensure no negative branch lengths exist in tree.  This can occur if
     * leaf heights given as a trait are incompatible with the existing tree.
     * If adjustTreeNodeHeightsInput is set to true, then all non-positive branch lengths
     * are coerced to be small positive lengths.
     */
    final static double EPSILON = 0.0000001;

    protected void adjustTreeNodeHeights(final Node node) {
        if (!node.isLeaf()) {
            for (final Node child : node.getChildren()) {
                adjustTreeNodeHeights(child);
            }
            for (final Node child : node.getChildren()) {
                if (adjustTreeNodeHeightsInput.get()) {
                    final double minHeight = child.getHeight() + EPSILON;
                    if (node.getHeight() < minHeight) {
                        node.setHeight(minHeight);
                    }
                } else {
                    final double minHeight = child.getHeight();
                    if (node.getHeight() < minHeight) {
                        node.setHeight(minHeight);
                    }
                }
            }
        }
    }

I don't see any use case for negative branch lengths, but if it's better for backwards compatibility to allow them, then a more precise solution would be to add support for an allowZeroBranchLengths attribute and modify the method as follows:

    /**
     * Ensure no negative branch lengths exist in tree.  This can occur if
     * leaf heights given as a trait are incompatible with the existing tree.
     * If allowZeroBranchLengths is false, then all non-positive branch lengths
     * are coerced to be small positive lengths.
     */
    final static double EPSILON = 0.0000001;

    protected void adjustTreeNodeHeights(final Node node) {
        if (!node.isLeaf() && adjustTreeNodeHeightsInput.get()) {
            for (final Node child : node.getChildren()) {
                adjustTreeNodeHeights(child);
            }
            for (final Node child : node.getChildren()) {
            if (allowZeroBranchLengthsInput.get()) {
                final double minHeight = child.getHeight();
                if (node.getHeight() < minHeight) {
                    node.setHeight(minHeight);
                }
            } else {
                 final double minHeight = child.getHeight() + EPSILON;
                if (node.getHeight() < minHeight) {
                    node.setHeight(minHeight);
                }
            }
        }
    }