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);
}
}
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!!");
}
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!");
}
}
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;
}
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;
}
Aggregations