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

Example 96 with Node

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

the class TypedSubtreeExchange method proposal.

@Override
public double proposal() {
    double logHR = 0.0;
    // Select source and destination nodes:
    Node srcNode, srcNodeParent, destNode, destNodeParent;
    if (isNarrowInput.get()) {
        // Narrow exchange selection:
        do {
            srcNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (srcNode.isRoot() || srcNode.getParent().isRoot());
        srcNodeParent = srcNode.getParent();
        destNode = getOtherChild(srcNodeParent.getParent(), srcNodeParent);
        destNodeParent = destNode.getParent();
    } else {
        // Wide exchange selection:
        do {
            srcNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (srcNode.isRoot());
        srcNodeParent = srcNode.getParent();
        do {
            destNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (destNode == srcNode || destNode.isRoot() || destNode.getParent() == srcNode.getParent());
        destNodeParent = destNode.getParent();
    }
    // Reject if substitution would result in negative branch lengths:
    if (destNode.getHeight() > srcNodeParent.getHeight() || srcNode.getHeight() > destNodeParent.getHeight())
        return Double.NEGATIVE_INFINITY;
    // Record probability of old colours:
    logHR += getBranchTypeProb(srcNode) + getBranchTypeProb(destNode);
    // Make changes to tree topology:
    replace(srcNodeParent, srcNode, destNode);
    replace(destNodeParent, destNode, srcNode);
    // Recolour branches involved:
    try {
        logHR -= retypeBranch(srcNode) + retypeBranch(destNode);
    } catch (NoValidPathException e) {
        return Double.NEGATIVE_INFINITY;
    }
    return logHR;
}

Example 97 with Node

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

the class TypedSubtreeExchangeRandom method proposal.

@Override
public double proposal() {
    double logHR = 0.0;
    // Select source and destination nodes:
    Node srcNode, srcNodeParent, destNode, destNodeParent;
    if (isNarrowInput.get()) {
        // Narrow exchange selection:
        do {
            srcNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (srcNode.isRoot() || srcNode.getParent().isRoot());
        srcNodeParent = srcNode.getParent();
        destNode = getOtherChild(srcNodeParent.getParent(), srcNodeParent);
        destNodeParent = destNode.getParent();
    } else {
        // Wide exchange selection:
        do {
            srcNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (srcNode.isRoot());
        srcNodeParent = srcNode.getParent();
        do {
            destNode = mtTree.getNode(Randomizer.nextInt(mtTree.getNodeCount()));
        } while (destNode == srcNode || destNode.isRoot() || destNode.getParent() == srcNode.getParent());
        destNodeParent = destNode.getParent();
    }
    // lengths:
    if (destNode.getHeight() > srcNodeParent.getHeight() || srcNode.getHeight() > destNodeParent.getHeight())
        return Double.NEGATIVE_INFINITY;
    // Record probability of old colours:
    logHR += getBranchTypeProb(srcNode) + getBranchTypeProb(destNode);
    // Make changes to tree topology:
    replace(srcNodeParent, srcNode, destNode);
    replace(destNodeParent, destNode, srcNode);
    // Recolour branches involved:
    logHR -= retypeBranch(srcNode) + retypeBranch(destNode);
    // Force rejection if colouring inconsistent:
    if ((((MultiTypeNode) srcNode).getFinalType() != ((MultiTypeNode) destNodeParent).getNodeType()) || (((MultiTypeNode) destNode).getFinalType() != ((MultiTypeNode) srcNodeParent).getNodeType()))
        return Double.NEGATIVE_INFINITY;
    return logHR;
}

Example 98 with Node

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

the class TypedWilsonBaldingEasy 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 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);
        ((MultiTypeNode) destNode).clearChanges();
        connectBranchToRoot(srcNode, destNode, newTime);
        mtTree.setRoot(srcNodeP);
        // Abort if colouring inconsistent:
        if (((MultiTypeNode) srcNode).getFinalType() != ((MultiTypeNode) destNode).getFinalType())
            return Double.NEGATIVE_INFINITY;
        // Update colour of root node:
        ((MultiTypeNode) srcNodeP).setNodeType(((MultiTypeNode) srcNode).getFinalType());
        // Final test of tree validity:
        if (!mtTree.isValid())
            return Double.NEGATIVE_INFINITY;
        // height changes:
        return Math.log(alpha * t_destNode) + (1.0 / alpha) * (newTime / t_destNode - 1.0) - Math.log(t_srcNodeG - Math.max(t_srcNode, t_srcNodeS));
    }
    if (srcNodeP.isRoot()) {
        // changes. (This would be an irreversible move.)
        if (((MultiTypeNode) srcNodeS).getChangeCount() > 0 || (((MultiTypeNode) srcNodeS).getNodeType() != ((MultiTypeNode) srcNode).getFinalType()))
            return Double.NEGATIVE_INFINITY;
        // 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);
        // Abort if new colouring is inconsistent:
        if (((MultiTypeNode) srcNodeP).getNodeType() != ((MultiTypeNode) srcNode).getFinalType())
            return Double.NEGATIVE_INFINITY;
        // Final test of tree validity:
        if (!mtTree.isValid())
            return Double.NEGATIVE_INFINITY;
        // height changes:
        return Math.log(t_destNodeP - Math.max(t_srcNode, t_destNode)) - Math.log(alpha * t_srcNodeS) - (1.0 / alpha) * (oldTime / t_srcNodeS - 1.0);
    }
    // NON-ROOT MOVE
    // 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);
    // Reject outright if new colouring inconsistent:
    if (((MultiTypeNode) srcNodeP).getNodeType() != ((MultiTypeNode) srcNode).getFinalType())
        return Double.NEGATIVE_INFINITY;
    // Final test of tree validity:
    if (!mtTree.isValid())
        return Double.NEGATIVE_INFINITY;
    // height changes:
    return Math.log(t_destNodeP - Math.max(t_srcNode, t_destNode)) - Math.log(t_srcNodeG - Math.max(t_srcNode, t_srcNodeS));
}

Example 99 with Node

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

the class BeerliFelsenstein method getPartialEventList.

/**
 * Assemble and return list of events excluding those on the edge between
 * node and its parent.
 *
 * @param excludedNode Tree node indicating edge to exclude.
 * @return event list
 */
private List<Event> getPartialEventList(Node excludedNode) {
    List<Event> eventList = Lists.newArrayList();
    // Collect all events
    for (Node node : mtTree.getNodesAsArray()) {
        if (node == excludedNode)
            continue;
        MultiTypeNode mtNode = (MultiTypeNode) node;
        Event event = new Event();
        event.time = node.getHeight();
        event.node = node;
        event.thisDeme = mtNode.getNodeType();
        if (node.isLeaf())
            event.type = EventType.SAMPLE;
        else {
            if (!node.getChildren().contains(excludedNode))
                event.type = EventType.COALESCENCE;
        }
        eventList.add(event);
        int thisDeme = mtNode.getNodeType();
        int prevDeme;
        for (int i = 0; i < mtNode.getChangeCount(); i++) {
            prevDeme = thisDeme;
            thisDeme = mtNode.getChangeType(i);
            Event changeEvent = new Event();
            changeEvent.type = EventType.MIGRATION;
            changeEvent.time = mtNode.getChangeTime(i);
            changeEvent.node = node;
            changeEvent.thisDeme = thisDeme;
            changeEvent.prevDeme = prevDeme;
            eventList.add(changeEvent);
        }
    }
    // Sort events according to times
    Collections.sort(eventList, new Comparator<Event>() {

        @Override
        public int compare(Event e1, Event e2) {
            if (e1.time < e2.time)
                return -1;
            if (e1.time > e2.time)
                return 1;
            return 0;
        }
    });
    return eventList;
}

Example 100 with Node

use of beast.evolution.tree.Node 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;
}