Examples with Node beast.evolution.tree.Node used on opensource projects

Example 56 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class CFOperator method collapseConversions.

/**
 * Take conversions which connect to edge above srcNode at times greater than
 * destTime and attach them instead to the lineage above destNode.
 *
 * Assumes topology has not yet been altered.
 *
 * @param srcNode   source node for move
 * @param destNode  dest node for move
 * @param destTime  new time of attachment of edge above srcNode to edge
 *                  above destNode
 * @return log probability of the collapsed attachments.
 */
protected double collapseConversions(Node srcNode, Node destNode, double destTime) {
    double logP = 0.0;
    boolean reverseRootMove = srcNode.getParent().isRoot();
    Node srcNodeP = srcNode.getParent();
    Node srcNodeS = getSibling(srcNode);
    double maxChildHeight = getMaxRootChildHeight();
    double volatileHeight = Math.max(maxChildHeight, destTime);
    // Collapse non-root conversions
    Node node = destNode;
    while (!node.isRoot() && node.getHeight() < srcNodeP.getHeight()) {
        double lowerBound = Math.max(destTime, node.getHeight());
        double upperBound = Math.min(node.getParent().getHeight(), srcNodeP.getHeight());
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getHeight1() > lowerBound && conv.getHeight1() < upperBound) {
                    if (conv.getNode1() == srcNode)
                        conv.setNode1(node);
                    if (conv.getNode1() == node && (!reverseRootMove || conv.getHeight1() < volatileHeight))
                        logP += Math.log(0.5);
                }
                if (conv.getHeight2() > lowerBound && conv.getHeight2() < upperBound) {
                    if (conv.getNode2() == srcNode)
                        conv.setNode2(node);
                    if (conv.getNode2() == node && (!reverseRootMove || conv.getNode1() != node || conv.getHeight1() < volatileHeight))
                        logP += Math.log(0.5);
                }
            }
        }
        node = node.getParent();
    }
    if (reverseRootMove) {
        double L = 2.0 * (srcNode.getParent().getHeight() - volatileHeight);
        double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
        List<Conversion> toRemove = new ArrayList<>();
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getHeight1() > volatileHeight)
                    toRemove.add(conv);
            }
        }
        logP += -Nexp + toRemove.size() * Math.log(Nexp);
        for (Conversion conv : toRemove) {
            logP += Math.log(1.0 / L) + getAffectedRegionProb(conv) + getEdgeCoalescenceProb(conv);
            acg.deleteConversion(conv);
        }
    }
    // Apply topology modifications.
    disconnectEdge(srcNode);
    connectEdge(srcNode, destNode, destTime);
    if (reverseRootMove && destTime < maxChildHeight)
        acg.setRoot(srcNodeS);
    return logP;
}

Example 57 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class CFSubtreeSlide method proposal.

@Override
public double proposal() {
    double logHGF = 0.0;
    double logHalf = Math.log(0.5);
    // Choose non-root node:
    Node srcNode = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    Node srcNodeP = srcNode.getParent();
    Node srcNodeS = getSibling(srcNode);
    double fMin = Math.min(scaleFactorInput.get(), 1.0 / scaleFactorInput.get());
    double f = fMin + Randomizer.nextDouble() * (1.0 / fMin - fMin);
    logHGF += -Math.log(f);
    double newHeight = srcNodeP.getHeight() * f;
    // Reject invalid moves:
    if (newHeight < srcNode.getHeight())
        return Double.NEGATIVE_INFINITY;
    if (f < 1.0) {
        // Search downwards for new attachment point:
        Node newSister = srcNodeS;
        while (newHeight < newSister.getHeight()) {
            if (newSister.isLeaf())
                return Double.NEGATIVE_INFINITY;
            newSister = Randomizer.nextBoolean() ? newSister.getLeft() : newSister.getRight();
            logHGF -= logHalf;
        }
        // Update topology:
        logHGF += collapseConversions(srcNode, newSister, newHeight);
    } else {
        // Search upwards for new attachment point:
        Node newSister = srcNodeS;
        while (getParentSkip(newSister, srcNodeP) != null && newHeight > getParentSkip(newSister, srcNodeP).getHeight()) {
            newSister = getParentSkip(newSister, srcNodeP);
            logHGF += logHalf;
        }
        // Update topology
        logHGF -= expandConversions(srcNode, newSister, newHeight);
    }
    assert !acg.isInvalid() : "CFSTS proposed invalid state.";
    return logHGF;
}

Example 58 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class CFUniform method proposal.

@Override
public double proposal() {
    double logHGF = 0.0;
    double logHalf = Math.log(0.5);
    // Select internal non-root node at random.
    Node node = acg.getNode(acg.getLeafNodeCount() + Randomizer.nextInt(acg.getInternalNodeCount()));
    Node leftChild = node.getLeft();
    Node rightChild = node.getRight();
    double oldHeight = node.getHeight();
    double maxChildHeight = Math.max(leftChild.getHeight(), rightChild.getHeight());
    // Choose new height:
    double newHeight;
    if (node.isRoot()) {
        double fMin = Math.min(scaleFactorInput.get(), 1.0 / scaleFactorInput.get());
        double fMax = 1.0 / fMin;
        double f = fMin + (fMax - fMin) * Randomizer.nextDouble();
        newHeight = node.getHeight() * f;
        logHGF += Math.log(1.0 / f);
        if (newHeight < maxChildHeight)
            return Double.NEGATIVE_INFINITY;
    } else {
        Node parent = node.getParent();
        newHeight = maxChildHeight + Randomizer.nextDouble() * (parent.getHeight() - maxChildHeight);
    }
    if (newHeight > oldHeight) {
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if (conv.getNode1() == node && conv.getHeight1() < newHeight) {
                    conv.setNode1(Randomizer.nextBoolean() ? leftChild : rightChild);
                    logHGF -= logHalf;
                }
                if (conv.getNode2() == node && conv.getHeight2() < newHeight) {
                    conv.setNode2(Randomizer.nextBoolean() ? leftChild : rightChild);
                    logHGF -= logHalf;
                }
            }
        }
        node.setHeight(newHeight);
        if (node.isRoot()) {
            // Draw a number of conversions
            double L = 2.0 * (newHeight - oldHeight);
            double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * deltaInput.get().getValue());
            int N = (int) Randomizer.nextPoisson(Nexp);
            // N! cancels
            logHGF -= -Nexp + N * Math.log(Nexp);
            for (int i = 0; i < N; i++) {
                Conversion conv = new Conversion();
                double u = L * Randomizer.nextDouble();
                if (u < 0.5 * L) {
                    conv.setNode1(leftChild);
                    conv.setHeight1(oldHeight + u);
                } else {
                    conv.setNode1(rightChild);
                    conv.setHeight1(oldHeight + u - 0.5 * L);
                }
                logHGF -= Math.log(1.0 / L) + coalesceEdge(conv) + drawAffectedRegion(conv);
                acg.addConversion(conv);
            }
        }
    } else {
        List<Conversion> toRemove = new ArrayList<>();
        for (Locus locus : acg.getLoci()) {
            for (Conversion conv : acg.getConversions(locus)) {
                if ((conv.getNode1() == leftChild || conv.getNode1() == rightChild) && conv.getHeight1() > newHeight) {
                    if (node.isRoot()) {
                        toRemove.add(conv);
                        continue;
                    } else {
                        conv.setNode1(node);
                        logHGF += logHalf;
                    }
                }
                if ((conv.getNode2() == leftChild || conv.getNode2() == rightChild) && conv.getHeight2() > newHeight) {
                    conv.setNode2(node);
                    logHGF += logHalf;
                }
            }
        }
        if (node.isRoot()) {
            double L = 2.0 * (oldHeight - newHeight);
            double Nexp = L * rhoInput.get().getValue() * (acg.getTotalSequenceLength() + acg.getLoci().size() * deltaInput.get().getValue());
            // N! cancels
            logHGF += -Nexp + toRemove.size() * Math.log(Nexp);
            for (Conversion conv : toRemove) {
                logHGF += Math.log(1.0 / L) + getEdgeCoalescenceProb(conv) + getAffectedRegionProb(conv);
                acg.deleteConversion(conv);
            }
        }
        node.setHeight(newHeight);
    }
    if (acg.isInvalid())
        throw new IllegalStateException("Something is broken: CFUniform proposed illegal state!");
    return logHGF;
}

Example 59 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class CFWilsonBalding method proposal.

@Override
public double proposal() {
    // Determine whether we can apply this operator:
    if (acg.getLeafNodeCount() < 3)
        return Double.NEGATIVE_INFINITY;
    // Select non-root node:
    Node srcNode;
    do {
        srcNode = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    } while (invalidSrcNode(srcNode));
    Node srcNodeP = srcNode.getParent();
    Node srcNodeS = getSibling(srcNode);
    double t_srcNode = srcNode.getHeight();
    double t_srcNodeP = srcNodeP.getHeight();
    double t_srcNodeS = srcNodeS.getHeight();
    // Select destination branch node:
    Node destNode;
    do {
        destNode = acg.getNode(Randomizer.nextInt(acg.getNodeCount()));
    } while (invalidDestNode(srcNode, destNode));
    Node destNodeP = destNode.getParent();
    double t_destNode = destNode.getHeight();
    if (destNode.isRoot()) {
        if (!includeRootInput.get())
            return Double.NEGATIVE_INFINITY;
        double logHGF = 0.0;
        double t_srcNodeG = srcNodeP.getParent().getHeight();
        logHGF += Math.log(1.0 / (t_srcNodeG - Math.max(t_srcNode, t_srcNodeS)));
        double newTime = t_destNode + Randomizer.nextExponential(1.0 / (alpha * t_destNode));
        logHGF -= Math.log(1.0 / (alpha * t_destNode)) - (1.0 / alpha) * (newTime / t_destNode - 1.0);
        // Randomly reconnect some of the conversions ancestral
        // to srcNode to the new edge above srcNode.
        logHGF -= expandConversions(srcNode, destNode, newTime);
        assert !acg.isInvalid() : "CFWB proposed invalid state.";
        return logHGF;
    }
    if (srcNodeP.isRoot()) {
        if (!includeRootInput.get())
            return Double.NEGATIVE_INFINITY;
        double logHGF = 0.0;
        logHGF += Math.log(1.0 / (alpha * t_srcNodeS)) - (1.0 / alpha) * (t_srcNodeP / t_srcNodeS - 1.0);
        double min_newTime = Math.max(t_srcNode, t_destNode);
        double t_destNodeP = destNodeP.getHeight();
        double newTime = min_newTime + (t_destNodeP - min_newTime) * Randomizer.nextDouble();
        logHGF -= Math.log(1.0 / (t_destNodeP - min_newTime));
        // Reconnect conversions on edge above srcNode older than
        // newTime to edges ancestral to destNode.
        logHGF += collapseConversions(srcNode, destNode, newTime);
        assert !acg.isInvalid() : "CFWB proposed invalid state.";
        return logHGF;
    }
    // Non-root move
    double logHGF = 0.0;
    double t_srcNodeG = srcNodeP.getParent().getHeight();
    logHGF += Math.log(1.0 / (t_srcNodeG - Math.max(t_srcNode, t_srcNodeS)));
    double min_newTime = Math.max(t_destNode, t_srcNode);
    double t_destNodeP = destNodeP.getHeight();
    double newTime = min_newTime + (t_destNodeP - min_newTime) * Randomizer.nextDouble();
    logHGF -= Math.log(1.0 / (t_destNodeP - min_newTime));
    if (newTime < srcNodeP.getHeight())
        logHGF += collapseConversions(srcNode, destNode, newTime);
    else
        logHGF -= expandConversions(srcNode, destNode, newTime);
    assert !acg.isInvalid() : "CFWB proposed invalid state.";
    return logHGF;
}

Example 60 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ConvertedEdgeFlip method proposal.

@Override
public double proposal() {
    if (acg.getTotalConvCount() == 0)
        return Double.NEGATIVE_INFINITY;
    Conversion recomb = chooseConversion();
    Node node1 = recomb.getNode1();
    Node node2 = recomb.getNode2();
    double height1 = recomb.getHeight1();
    double height2 = recomb.getHeight2();
    if (node1 == node2 || node2.isRoot())
        return Double.NEGATIVE_INFINITY;
    if (height1 < node2.getHeight() || height1 > node2.getParent().getHeight() || height2 < node1.getHeight() || height2 > node1.getParent().getHeight())
        return Double.NEGATIVE_INFINITY;
    recomb.setNode1(node2);
    recomb.setNode2(node1);
    return 0.0;
}