use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.
the class ARGTree method toGraphString.
public String toGraphString() {
StringBuffer sb = new StringBuffer();
for (Node node : nodes) {
sb.append(node.number);
if (node.leftParent != null)
sb.append(" " + node.leftParent.number);
else
sb.append(" 0");
if (node.rightParent != null)
sb.append(" " + node.rightParent.number);
else
sb.append(" 0");
if (node.leftChild != null)
sb.append(" " + node.leftChild.number);
else
sb.append(" 0");
if (node.rightChild != null)
sb.append(" " + node.rightChild.number);
else
sb.append(" 0");
if (node.taxon != null)
sb.append(" " + node.taxon.toString());
sb.append("\n");
}
sb.append("Root = " + ((Node) getRoot()).number + "\n");
return new String(sb);
}
use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.
the class ObsoleteARGNewEventOperator method AddOperation.
private double AddOperation() throws ARGOperatorFailedException {
double logq = 0;
// Draw attachment point for new bifurcation
// System.err.println("Starting add operation.");
ArrayList<NodeRef> potentialBifurcationChildren = new ArrayList<NodeRef>();
ArrayList<NodeRef> potentialReassortmentChildren = new ArrayList<NodeRef>();
double treeHeight = arg.getNodeHeight(arg.getRoot());
double newBifurcationHeight = treeHeight * MathUtils.nextDouble();
double newReassortmentHeight = treeHeight * MathUtils.nextDouble();
if (newReassortmentHeight > newBifurcationHeight) {
double temp = newReassortmentHeight;
newReassortmentHeight = newBifurcationHeight;
newBifurcationHeight = temp;
}
// logq += 2.0 * Math.log(treeHeight); // Uniform before sorting
int totalPotentialBifurcationChildren = findPotentialAttachmentPoints(newBifurcationHeight, potentialBifurcationChildren);
int totalPotentialReassortmentChildren = findPotentialAttachmentPoints(newReassortmentHeight, potentialReassortmentChildren);
if (totalPotentialBifurcationChildren == 0 || totalPotentialReassortmentChildren == 0) {
throw new RuntimeException("Unable to find attachment points.");
}
Node recNode = (Node) potentialReassortmentChildren.get(MathUtils.nextInt(totalPotentialReassortmentChildren));
// logq += Math.log(totalPotentialReassortmentChildren);
Node recParentL = recNode.leftParent;
Node recParentR = recNode.rightParent;
Node recParent = recParentL;
if (recParentL != recParentR) {
if (arg.getNodeHeight(recParentL) < newReassortmentHeight)
recParent = recParentR;
else if (arg.getNodeHeight(recParentR) > newReassortmentHeight && MathUtils.nextDouble() > 0.5)
recParent = recParentR;
}
Node sisNode = (Node) potentialBifurcationChildren.get(MathUtils.nextInt(totalPotentialBifurcationChildren));
// logq += Math.log(totalPotentialBifurcationChildren);
Node sisParentL = sisNode.leftParent;
Node sisParentR = sisNode.rightParent;
Node sisParent = sisParentL;
if (sisParentL != sisParentR) {
if (arg.getNodeHeight(sisParentL) < newBifurcationHeight)
sisParent = sisParentR;
else if (arg.getNodeHeight(sisParentR) > newBifurcationHeight && MathUtils.nextDouble() > 0.5)
sisParent = sisParentR;
}
double newBifurcationRateCategory = 1.0;
double newReassortmentRateCategory = 1.0;
Node newBifurcation = arg.new Node();
newBifurcation.heightParameter = new Parameter.Default(newBifurcationHeight);
newBifurcation.rateParameter = new Parameter.Default(newBifurcationRateCategory);
newBifurcation.setupHeightBounds();
Node newReassortment = arg.new Node();
newReassortment.bifurcation = false;
newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
newReassortment.rateParameter = new Parameter.Default(newReassortmentRateCategory);
newReassortment.setupHeightBounds();
arg.beginTreeEdit();
if (sisParent.bifurcation)
arg.singleRemoveChild(sisParent, sisNode);
else
arg.doubleRemoveChild(sisParent, sisNode);
if (sisNode != recNode) {
if (recParent.bifurcation)
arg.singleRemoveChild(recParent, recNode);
else
arg.doubleRemoveChild(recParent, recNode);
}
if (sisParent.bifurcation)
arg.singleAddChild(sisParent, newBifurcation);
else
arg.doubleAddChild(sisParent, newBifurcation);
if (sisNode != recNode)
arg.singleAddChild(newBifurcation, sisNode);
arg.doubleAddChild(newReassortment, recNode);
Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
// logq +=
drawRandomPartitioning(partitioning);
if (sisNode != recNode) {
arg.addChildAsRecombinant(newBifurcation, recParent, newReassortment, partitioning);
} else {
arg.addChildAsRecombinant(newBifurcation, newBifurcation, newReassortment, partitioning);
}
arg.expandARGWithRecombinant(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
// System.err.println("point 2");
arg.pushTreeSizeChangedEvent();
arg.endTreeEdit();
// try {
// arg.checkTreeIsValid();
// } catch (MutableTree.InvalidTreeException ite) {
// throw new RuntimeException(ite.toString() + "\n" + arg.toString()
// + "\n" + Tree.Utils.uniqueNewick(arg, arg.getRoot()));
// }
// double d1 = findPotentialNodesToRemove(null);
// System.err.printf("d1 = %5.4f\n",d1);
// logq -= Math.log(findPotentialNodesToRemove(null));
// if (!(recNode.leftParent.isBifurcation() && recNode.rightParent.isRoot()) &&
// !(recNode.rightParent.isBifurcation() && recNode.leftParent.isRoot()))
// logq -= Math.log(2.0);
// System.err.println("End add ARG operation.");
// int nodes = arg.getInternalNodeCount() - 1;
//
/* int i = findPotentialNodesToRemove(null);
if (i==0) {
System.err.println("why can't i remove this one?");
System.err.println("graph:"+arg.toGraphString());
System.exit(1);
return Double.NEGATIVE_INFINITY ;
}*/
logq = 0;
// System.err.println("After add:\n"+arg.toGraphString());
return logq;
}
use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.
the class ObsoleteARGNewEventOperator method findPotentialNodesToRemove.
private int findPotentialNodesToRemove(ArrayList<NodeRef> list) {
int count = 0;
int n = arg.getNodeCount();
// int max = arg.getReassortmentNodeCount();
Node root = (Node) arg.getRoot();
for (int i = 0; i < n; i++) {
Node node = (Node) arg.getNode(i);
if (node.isReassortment() && ((node.leftParent != root && node.leftParent.isBifurcation()) || (node.rightParent != root && node.rightParent.isBifurcation()))) // && !(
// (node.leftParent == root && node.rightParent.isReassortment()) ||
// (node.rightParent== root && node.leftParent.isReassortment()))
// && !(node.leftParent == root && node.rightParent == root)
// && (node.leftParent.isBifurcation() || node.rightParent.isBifurcation())
{
if (list != null)
list.add(node);
count++;
}
}
// System.err.printf("max = %d, available = %d\n",max,count);
return count;
}
use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.
the class ObsoleteARGNewEventOperator method sanityCheck.
// if( !recParent1.bifurcation || recParent1.isRoot() ) { // One orientation is valid
public void sanityCheck() {
int len = arg.getNodeCount();
for (int i = 0; i < len; i++) {
Node node = (Node) arg.getNode(i);
if (node.bifurcation) {
boolean equalChild = (node.leftChild == node.rightChild);
if ((equalChild && node.leftChild != null)) {
if (!node.leftChild.bifurcation && ((node.leftChild).leftParent == node))
;
else {
System.err.println("Node " + (i + 1) + " is insane.");
System.err.println(arg.toGraphString());
System.exit(-1);
}
}
} else {
if ((node.leftChild != node.rightChild)) {
System.err.println("Node " + (i + 1) + " is insane.");
System.err.println(arg.toGraphString());
System.exit(-1);
}
}
if (!node.isRoot()) {
double d;
d = node.getHeight();
}
}
}
use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.
the class ARGSwapOperator method doNarrowSwap.
private void doNarrowSwap(NarrowSwap swap) {
arg.beginTreeEdit();
String before = arg.toARGSummary();
if (swap.c == swap.pb) {
Node c = (Node) swap.c;
Node p = (Node) swap.p;
Node gp = (Node) swap.gp;
if (c.leftParent == p) {
c.leftParent = gp;
c.rightParent = p;
} else {
c.leftParent = p;
c.rightParent = gp;
}
} else if (arg.getChild(swap.p, 0) == arg.getChild(swap.p, 1)) {
Node p = (Node) swap.p;
Node c = (Node) swap.c;
if (MathUtils.nextBoolean())
p.leftChild = c.leftParent = null;
else
p.rightChild = c.rightParent = null;
arg.removeChild(swap.gp, swap.pb);
arg.singleAddChild(swap.gp, swap.c);
arg.singleAddChild(swap.p, swap.pb);
} else {
arg.removeChild(swap.gp, swap.pb);
arg.removeChild(swap.p, swap.c);
arg.singleAddChild(swap.gp, swap.c);
arg.singleAddChild(swap.p, swap.pb);
}
assert nodeCheck() : swap + " " + before + " " + arg.toARGSummary();
arg.pushTreeChangedEvent(swap.gp);
arg.pushTreeChangedEvent(swap.p);
arg.endTreeEdit();
try {
arg.checkTreeIsValid();
} catch (MutableTree.InvalidTreeException ite) {
System.out.println(swap);
System.out.println(before);
System.err.println(ite.getMessage());
System.exit(-1);
} catch (NullPointerException e) {
System.out.println(swap);
System.out.println(before);
System.err.println(e.getMessage());
System.exit(-1);
}
}
Aggregations