Search in sources :

Example 56 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class AbstractImportanceDistributionOperator method removeChildren.

/**
     * @param parent
     * @warning assumes strictly bifurcating trees
     */
private void removeChildren(NodeRef parent) {
    // assumes strictly bifurcating trees
    NodeRef child = tree.getChild(parent, 0);
    if (child != null) {
        tree.removeChild(parent, child);
    }
    child = tree.getChild(parent, 1);
    if (child != null) {
        tree.removeChild(parent, child);
    }
}
Also used : NodeRef(dr.evolution.tree.NodeRef)

Example 57 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ColouredExchangeOperator method narrow.

/**
     * WARNING: Assumes strictly bifurcating tree.
     */
public void narrow() {
    NodeRef i = null, iP = null, j = null, jP = null;
    int tries = 0;
    while (tries < MAX_TRIES) {
        i = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        while (tree.getRoot() == i || tree.getParent(i) == tree.getRoot()) {
            i = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        }
        iP = tree.getParent(i);
        jP = tree.getParent(iP);
        j = tree.getChild(jP, 0);
        if (j == iP) {
            j = tree.getChild(jP, 1);
        }
        if ((tree.getNodeHeight(j) < tree.getNodeHeight(iP)) && (tree.getNodeHeight(i) < tree.getNodeHeight(jP))) {
            break;
        }
        tries += 1;
    }
    //Eupdate
    if (tries < MAX_TRIES) {
        eupdate(i, j, iP, jP);
        tree.pushTreeChangedEvent(iP);
        tree.pushTreeChangedEvent(jP);
    } else
        throw new RuntimeException("Couldn't find valid narrow move on this tree!!");
}
Also used : NodeRef(dr.evolution.tree.NodeRef)

Example 58 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ColouredExchangeOperator method wide.

/**
     * WARNING: Assumes strictly bifurcating tree.
     */
public void wide() {
    NodeRef i = null, iP = null, j = null, jP = null;
    int tries = 0;
    while (tries < MAX_TRIES) {
        i = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        while (tree.getRoot() == i) {
            i = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        }
        j = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        while (j == i || j == tree.getRoot()) {
            j = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        }
        iP = tree.getParent(i);
        jP = tree.getParent(j);
        if ((iP != jP) && (i != jP) && (j != iP) && (tree.getNodeHeight(j) < tree.getNodeHeight(iP)) && (tree.getNodeHeight(i) < tree.getNodeHeight(jP))) {
            break;
        }
        tries += 1;
    }
    //Eupdate
    if (tries < MAX_TRIES) {
        eupdate(i, j, iP, jP);
    } else {
        throw new RuntimeException("Couldn't find valid wide move on this tree!");
    }
}
Also used : NodeRef(dr.evolution.tree.NodeRef)

Example 59 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ColouredSubtreeSlideOperator method intersectingEdges.

private int intersectingEdges(Tree tree, NodeRef node, double height, List<NodeRef> directChildren) {
    NodeRef parent = tree.getParent(node);
    if (tree.getNodeHeight(parent) < height)
        return 0;
    if (tree.getNodeHeight(node) < height) {
        if (directChildren != null)
            directChildren.add(node);
        return 1;
    }
    int count = 0;
    for (int i = 0; i < tree.getChildCount(node); i++) {
        count += intersectingEdges(tree, tree.getChild(node, i), height, directChildren);
    }
    return count;
}
Also used : NodeRef(dr.evolution.tree.NodeRef)

Example 60 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ColouredSubtreeSlideOperator method doOperation.

/**
     * Do a probablistic subtree slide move.
     *
     * @return the log-transformed hastings ratio
     */
public double doOperation() {
    double logP = colouringModel.getTreeColouringWithProbability().getLogProbabilityDensity();
    double logq;
    NodeRef i, newParent, newChild;
    // 1. choose a random node avoiding root
    do {
        i = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
    } while (tree.getRoot() == i);
    NodeRef iP = tree.getParent(i);
    NodeRef CiP = getOtherChild(tree, iP, i);
    NodeRef PiP = tree.getParent(iP);
    // 2. choose a delta to move
    double delta = getDelta();
    double oldHeight = tree.getNodeHeight(iP);
    double newHeight = oldHeight + delta;
    // 3. if the move is up
    if (delta > 0) {
        // 3.1 if the topology will change
        if (PiP != null && tree.getNodeHeight(PiP) < newHeight) {
            // find new parent
            newParent = PiP;
            newChild = iP;
            while (tree.getNodeHeight(newParent) < newHeight) {
                newChild = newParent;
                newParent = tree.getParent(newParent);
                if (newParent == null)
                    break;
            }
            tree.beginTreeEdit();
            // 3.1.1 if creating a new root
            if (tree.isRoot(newChild)) {
                tree.removeChild(iP, CiP);
                // remove iP from in between PiP and CiP
                tree.removeChild(PiP, iP);
                // add new edge from iP to newChild (old root)
                tree.addChild(iP, newChild);
                // formerly two edges
                tree.addChild(PiP, CiP);
                tree.setRoot(iP);
            //System.err.println("Creating new root!");
            } else // 3.1.2 no new root
            {
                tree.removeChild(iP, CiP);
                // remove iP from in between PiP and CiP
                tree.removeChild(PiP, iP);
                // split edge: first remove old one
                tree.removeChild(newParent, newChild);
                // split edge: put iP in middle
                tree.addChild(iP, newChild);
                // formerly two edges
                tree.addChild(PiP, CiP);
                // split edge: put iP in middle
                tree.addChild(newParent, iP);
            //System.err.println("No new root!");
            }
            tree.setNodeHeight(iP, newHeight);
            tree.endTreeEdit();
            // 3.1.3 count the hypothetical sources of this destination.
            int possibleSources = intersectingEdges(tree, newChild, oldHeight, null);
            //System.out.println("possible sources = " + possibleSources);
            logq = Math.log(1.0 / (double) possibleSources);
        } else {
            // just change the node height
            tree.setNodeHeight(iP, newHeight);
            logq = 0.0;
        }
    } else // 4 if we are sliding the subtree down.
    {
        // 4.0 is it a valid move?
        if (tree.getNodeHeight(i) > newHeight) {
            return Double.NEGATIVE_INFINITY;
        }
        // 4.1 will the move change the topology
        if (tree.getNodeHeight(CiP) > newHeight) {
            List<NodeRef> newChildren = new ArrayList<NodeRef>();
            int possibleDestinations = intersectingEdges(tree, CiP, newHeight, newChildren);
            // if no valid destinations then return a failure
            if (newChildren.size() == 0) {
                return Double.NEGATIVE_INFINITY;
            }
            // pick a random parent/child destination edge uniformly from options
            int childIndex = MathUtils.nextInt(newChildren.size());
            newChild = newChildren.get(childIndex);
            newParent = tree.getParent(newChild);
            tree.beginTreeEdit();
            // 4.1.1 if iP was root
            if (tree.isRoot(iP)) {
                // new root is CiP
                tree.removeChild(iP, CiP);
                tree.removeChild(newParent, newChild);
                tree.addChild(iP, newChild);
                tree.addChild(newParent, iP);
                tree.setRoot(CiP);
            //System.err.println("DOWN: Creating new root!");
            } else {
                tree.removeChild(iP, CiP);
                tree.removeChild(PiP, iP);
                tree.removeChild(newParent, newChild);
                tree.addChild(iP, newChild);
                tree.addChild(PiP, CiP);
                tree.addChild(newParent, iP);
            //System.err.println("DOWN: no new root!");
            }
            tree.setNodeHeight(iP, newHeight);
            tree.endTreeEdit();
            logq = Math.log((double) possibleDestinations);
        } else {
            tree.setNodeHeight(iP, newHeight);
            logq = 0.0;
        }
    }
    if (swapRates) {
        NodeRef j = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        if (j != i) {
            double tmp = tree.getNodeRate(i);
            tree.setNodeRate(i, tree.getNodeRate(j));
            tree.setNodeRate(j, tmp);
        }
    }
    if (swapTraits) {
        NodeRef j = tree.getNode(MathUtils.nextInt(tree.getNodeCount()));
        if (j != i) {
            double tmp = tree.getNodeTrait(i, "trait");
            tree.setNodeTrait(i, "trait", tree.getNodeTrait(j, "trait"));
            tree.setNodeTrait(j, "trait", tmp);
        }
    }
    if (logq == Double.NEGATIVE_INFINITY)
        throw new RuntimeException("invalid slide");
    colouringModel.resample();
    double logQ = colouringModel.getTreeColouringWithProbability().getLogProbabilityDensity();
    return logq + logP - logQ;
}
Also used : NodeRef(dr.evolution.tree.NodeRef) ArrayList(java.util.ArrayList)

Aggregations

NodeRef (dr.evolution.tree.NodeRef)426 ArrayList (java.util.ArrayList)38 LinkedList (java.util.LinkedList)20 Taxon (dr.evolution.util.Taxon)18 Tree (dr.evolution.tree.Tree)14 TreeModel (dr.evomodel.tree.TreeModel)14 Parameter (dr.inference.model.Parameter)14 Clade (dr.evolution.tree.Clade)13 MutableTree (dr.evolution.tree.MutableTree)9 Node (dr.evomodel.arg.ARGModel.Node)9 MultivariateTraitTree (dr.evolution.tree.MultivariateTraitTree)8 BranchMapModel (dr.evomodel.epidemiology.casetocase.BranchMapModel)8 BitSet (java.util.BitSet)8 FixedBitSet (jebl.util.FixedBitSet)8 FlexibleTree (dr.evolution.tree.FlexibleTree)7 AbstractCase (dr.evomodel.epidemiology.casetocase.AbstractCase)7 Matrix (dr.math.matrixAlgebra.Matrix)6 CompoundParameter (dr.inference.model.CompoundParameter)5 TimeScale (dr.evolution.util.TimeScale)4 MicrosatelliteSamplerTreeModel (dr.evomodel.tree.MicrosatelliteSamplerTreeModel)4