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

Example 31 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class BeerliFelsenstein method proposal.

@Override
public double proposal() {
    double logHR = 0.0;
    // Select non-root node at random
    Node node = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()) - 1);
    MultiTypeNode mtNode = (MultiTypeNode) node;
    // Keep copies of node and its sister for reverse move prob calculation
    MultiTypeNode mtNodeOld = mtNode.shallowCopy();
    MultiTypeNode mtNodeOldSis = (MultiTypeNode) getOtherChild(node.getParent(), node);
    // Assemble partial event list
    List<Event> eventList = getPartialEventList(node);
    double oldRootHeight = mtTree.getRoot().getHeight();
    // Topology changes to turn tree into partial tree
    Node nodeParent = node.getParent();
    if (!nodeParent.isRoot())
        disconnectBranch(node);
    else {
        Node sister = getOtherChild(nodeParent, node);
        nodeParent.removeChild(sister);
    }
    // Pre-calculate total lineage migration propensities
    double[] migProp = new double[migModel.getNTypes()];
    for (int d = 0; d < migModel.getNTypes(); d++) {
        migProp[d] = 0.0;
        for (int dp = 0; dp < migModel.getNTypes(); dp++) {
            if (d == dp)
                continue;
            migProp[d] += migModel.getBackwardRate(d, dp);
        }
    }
    List<Set<Node>> nodesOfType = Lists.newArrayList();
    for (int i = 0; i < migModel.getNTypes(); i++) nodesOfType.add(new HashSet<Node>());
    mtNode.clearChanges();
    double coalTime = Double.NaN;
    for (int eidx = 0; eidx < eventList.size(); eidx++) {
        Event event = eventList.get(eidx);
        double intervalEndTime;
        if (eidx < eventList.size() - 1)
            intervalEndTime = eventList.get(eidx + 1).time;
        else
            intervalEndTime = oldRootHeight;
        switch(event.type) {
            case COALESCENCE:
                nodesOfType.get(event.thisDeme).removeAll(event.node.getChildren());
                nodesOfType.get(event.thisDeme).add(event.node);
                break;
            case SAMPLE:
                nodesOfType.get(event.thisDeme).add(event.node);
                break;
            case MIGRATION:
                nodesOfType.get(event.prevDeme).remove(event.node);
                nodesOfType.get(event.thisDeme).add(event.node);
                break;
        }
        double t = Math.max(event.time, node.getHeight());
        // Early exit
        if (t >= intervalEndTime)
            continue;
        int deme = mtNode.getNodeType();
        while (true) {
            // Calculate coalescent propensity
            double coalProp = nodesOfType.get(deme).size() / migModelSC.getPopSize(deme);
            // Select event time
            double dt = Randomizer.nextExponential(coalProp + migProp[deme]);
            t += dt;
            if (t > intervalEndTime)
                break;
            // HR waiting time contribution
            logHR += -(coalProp + migProp[deme]) * dt;
            double u = Randomizer.nextDouble() * (coalProp + migProp[deme]);
            if (u < coalProp) {
                // Coalescence
                // Select edge to coalesce with
                Node coalNode = (Node) selectRandomElement(nodesOfType.get(deme));
                // HR event contribution
                logHR += Math.log(1.0 / migModelSC.getPopSize(deme));
                // Implement coalescence
                coalTime = t;
                connectBranch(node, coalNode, coalTime);
                break;
            } else {
                // Migration
                u -= coalProp;
                int toDeme;
                for (toDeme = 0; toDeme < migModel.getNTypes(); toDeme++) {
                    if (toDeme == deme)
                        continue;
                    u -= migModel.getBackwardRate(deme, toDeme);
                    if (u < 0)
                        break;
                }
                // HR event contribution
                logHR += Math.log(migModel.getBackwardRate(deme, toDeme));
                // Implelent migration
                mtNode.addChange(toDeme, t);
                deme = toDeme;
            }
        }
        // Continue to next interval if no coalescence has occurred
        if (!Double.isNaN(coalTime))
            break;
    }
    if (Double.isNaN(coalTime)) {
        // Continue simulation beyond old root of tree
        double t = oldRootHeight;
        int deme = mtNode.getFinalType();
        MultiTypeNode mtNodeSis = (MultiTypeNode) eventList.get(eventList.size() - 1).node;
        int demeSis = mtNodeSis.getFinalType();
        while (true) {
            // Calculate coalescent propensity
            double coalProp;
            if (deme == demeSis)
                coalProp = 1.0 / migModelSC.getPopSize(deme);
            else
                coalProp = 0.0;
            double totalProp = coalProp + migProp[deme] + migProp[demeSis];
            double dt = Randomizer.nextExponential(totalProp);
            // HR waiting time contribution
            logHR += -totalProp * dt;
            t += dt;
            double u = Randomizer.nextDouble() * totalProp;
            if (u < coalProp) {
                // Coalescence
                logHR += Math.log(1.0 / migModelSC.getPopSize(deme));
                coalTime = t;
                nodeParent.addChild(mtNodeSis);
                mtTree.setRoot(nodeParent);
                break;
            } else {
                // Migration
                u -= coalProp;
                if (u < migProp[deme]) {
                    // Migration in main lineage
                    int toDeme;
                    for (toDeme = 0; toDeme < migModel.getNTypes(); toDeme++) {
                        if (toDeme == deme)
                            continue;
                        u -= migModel.getBackwardRate(deme, toDeme);
                        if (u < 0)
                            break;
                    }
                    // HR contribution
                    logHR += Math.log(migModel.getBackwardRate(deme, toDeme));
                    mtNode.addChange(toDeme, t);
                    deme = toDeme;
                } else {
                    // Migration in sister lineage
                    int toDeme;
                    for (toDeme = 0; toDeme < migModel.getNTypes(); toDeme++) {
                        if (toDeme == demeSis)
                            continue;
                        u -= migModel.getBackwardRate(demeSis, toDeme);
                        if (u < 0)
                            break;
                    }
                    // HR contribution
                    logHR += Math.log(migModel.getBackwardRate(demeSis, toDeme));
                    mtNodeSis.addChange(toDeme, t);
                    demeSis = toDeme;
                }
            }
        }
    }
    return logHR;
}

Example 32 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class RandomRetypeOperator method retypeBranch.

/**
 * Retype branch between srcNode and its parent with rate fixed by the
 * tuning parameter mu.
 *
 * @param srcNode
 * @return Probability of branch typing
 */
protected double retypeBranch(Node srcNode) {
    double mu = muInput.get();
    Node srcNodeParent = srcNode.getParent();
    double t_srcNode = srcNode.getHeight();
    double t_srcNodeParent = srcNodeParent.getHeight();
    int srcNodeType = ((MultiTypeNode) srcNode).getNodeType();
    // Clear existing changes in preparation for adding replacements:
    ((MultiTypeNode) srcNode).clearChanges();
    double t = t_srcNode;
    int lastType = srcNodeType;
    while (t < t_srcNodeParent) {
        // Determine time to next migration event:
        t += Randomizer.nextExponential(mu);
        if (t < t_srcNodeParent) {
            // Select new colour:
            int newType = Randomizer.nextInt(migModel.getNTypes() - 1);
            if (newType >= lastType)
                newType += 1;
            ((MultiTypeNode) srcNode).addChange(newType, t);
            lastType = newType;
        }
    }
    // Return log of branch type probability:
    return -mu * (t_srcNodeParent - t_srcNode) + ((MultiTypeNode) srcNode).getChangeCount() * Math.log(mu / (migModel.getNTypes() - 1));
}

Example 33 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class TypedWilsonBaldingRandom method retypeRootBranches.

/**
 * Retype branches between srcNode and the root (srcNode's
 * parent) and between the root and srcNode's sister with a rate fixed
 * by the tuning parameter mu.
 *
 * @param srcNode
 * @param nChangesNode
 * @param nChangesSister
 * @return Probability of branch typing.
 */
private double retypeRootBranches(Node srcNode) {
    Node root = srcNode.getParent();
    Node srcNodeSister = getOtherChild(root, srcNode);
    // Recolour first branch:
    double logP = retypeBranch(srcNode);
    // Adjust colour of root node:
    ((MultiTypeNode) root).setNodeType(((MultiTypeNode) srcNode).getFinalType());
    return logP + retypeBranch(srcNodeSister);
}

Example 34 with MultiTypeNode

use of beast.evolution.tree.MultiTypeNode in project MultiTypeTree by tgvaughan.

the class TypedWilsonBaldingRandom method proposal.

@Override
public double proposal() {
    // Check that operator can be applied to tree:
    if (mtTree.getLeafNodeCount() < 3)
        throw new IllegalStateException("Tree too small for" + " ColouredWilsonBaldingRandom operator.");
    // Select source node:
    Node srcNode;
    do {
        srcNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
    } while (invalidSrcNode(srcNode));
    Node srcNodeP = srcNode.getParent();
    Node srcNodeS = getOtherChild(srcNodeP, srcNode);
    double t_srcNode = srcNode.getHeight();
    double t_srcNodeP = srcNodeP.getHeight();
    double t_srcNodeS = srcNodeS.getHeight();
    // Select destination branch node:
    Node destNode;
    do {
        destNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
    } while (invalidDestNode(srcNode, destNode));
    Node destNodeP = destNode.getParent();
    double t_destNode = destNode.getHeight();
    if (destNode.isRoot()) {
        // FORWARD ROOT MOVE
        // Record probability of current configuration:
        double logHR = getBranchTypeProb(srcNode);
        // Record srcNode grandmother height:
        double t_srcNodeG = srcNodeP.getParent().getHeight();
        // Choose new root height:
        double newTime = t_destNode + Randomizer.nextExponential(1.0 / (alpha * t_destNode));
        // Implement tree changes:
        disconnectBranch(srcNode);
        connectBranchToRoot(srcNode, destNode, newTime);
        mtTree.setRoot(srcNodeP);
        // Recolour root branches, incorporating probability of new branch
        // into HR:
        logHR -= retypeRootBranches(srcNode);
        // Abort if colouring inconsistent:
        if (((MultiTypeNode) srcNodeP).getNodeType() != ((MultiTypeNode) destNode).getFinalType())
            return Double.NEGATIVE_INFINITY;
        // Incorporate HR contribution of tree topology and node
        // height changes:
        logHR += Math.log(alpha * t_destNode) + (1.0 / alpha) * (newTime / t_destNode - 1.0) - Math.log(t_srcNodeG - Math.max(t_srcNode, t_srcNodeS));
        return logHR;
    }
    if (srcNodeP.isRoot()) {
        // BACKWARD ROOT MOVE
        // Record probability of current configuration:
        double logHR = getRootBranchTypeProb(srcNode);
        // Record old srcNode parent height:
        double oldTime = t_srcNodeP;
        // Choose height of new attachement point:
        double min_newTime = Math.max(t_srcNode, t_destNode);
        double t_destNodeP = destNodeP.getHeight();
        double span = t_destNodeP - min_newTime;
        double newTime = min_newTime + span * Randomizer.nextDouble();
        // Implement tree changes:
        disconnectBranchFromRoot(srcNode);
        connectBranch(srcNode, destNode, newTime);
        srcNodeS.setParent(null);
        mtTree.setRoot(srcNodeS);
        // Recolour new branch, incorporating probability of new branch
        // into HR:
        logHR -= retypeBranch(srcNode);
        // Abort if new colouring is inconsistent:
        if (((MultiTypeNode) srcNodeP).getNodeType() != ((MultiTypeNode) srcNode).getFinalType())
            return Double.NEGATIVE_INFINITY;
        // Incorporate HR contribution of tree topology and node
        // height changes:
        logHR += Math.log(t_destNodeP - Math.max(t_srcNode, t_destNode)) - Math.log(alpha * t_srcNodeS) - (1.0 / alpha) * (oldTime / t_srcNodeS - 1.0);
        return logHR;
    }
    // NON-ROOT MOVE
    // Record probability of old configuration:
    double logHR = getBranchTypeProb(srcNode);
    // Record srcNode grandmother height:
    double t_srcNodeG = srcNodeP.getParent().getHeight();
    // Choose height of new attachment point:
    double min_newTime = Math.max(t_destNode, t_srcNode);
    double t_destNodeP = destNodeP.getHeight();
    double span = t_destNodeP - min_newTime;
    double newTime = min_newTime + span * Randomizer.nextDouble();
    // Implement tree changes:
    disconnectBranch(srcNode);
    connectBranch(srcNode, destNode, newTime);
    // Recolour new branch:
    logHR -= retypeBranch(srcNode);
    // Reject outright if new colouring inconsistent:
    if (((MultiTypeNode) srcNodeP).getNodeType() != ((MultiTypeNode) srcNode).getFinalType())
        return Double.NEGATIVE_INFINITY;
    // Incorporate HR contribution of tree topology and node
    // height changes:
    logHR += Math.log(t_destNodeP - Math.max(t_srcNode, t_destNode)) - Math.log(t_srcNodeG - Math.max(t_srcNode, t_srcNodeS));
    return logHR;
}