Examples with Node dr.evomodel.arg.ARGModel.Node used on opensource projects

Example 6 with Node

use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.

the class ARGSwapOperator method bifurcationSwap.

private double bifurcationSwap(NodeRef x) {
    Node startNode = (Node) x;
    //		Node keepChild = startNode.leftChild;
    Node moveChild = startNode.rightChild;
    if (MathUtils.nextBoolean()) {
        //			keepChild = moveChild;
        moveChild = startNode.leftChild;
    }
    ArrayList<NodeRef> possibleNodes = new ArrayList<NodeRef>(arg.getNodeCount());
    findNodesAtHeight(possibleNodes, startNode.getHeight());
    assert !possibleNodes.contains(startNode);
    assert possibleNodes.size() > 0;
    Node swapNode = (Node) possibleNodes.get(MathUtils.nextInt(possibleNodes.size()));
    Node swapNodeParent = swapNode.leftParent;
    arg.beginTreeEdit();
    String before = arg.toARGSummary();
    if (swapNode.bifurcation) {
        swapNodeParent = swapNode.leftParent;
        arg.singleRemoveChild(startNode, moveChild);
        if (swapNodeParent.bifurcation) {
            arg.singleRemoveChild(swapNodeParent, swapNode);
            arg.singleAddChild(swapNodeParent, moveChild);
        } else {
            arg.doubleRemoveChild(swapNodeParent, swapNode);
            arg.doubleAddChild(swapNodeParent, moveChild);
        }
        arg.singleAddChild(startNode, swapNode);
    } else {
        boolean leftSide = true;
        boolean[] sideOk = { swapNode.leftParent.getHeight() > startNode.getHeight(), swapNode.rightParent.getHeight() > startNode.getHeight() };
        if (sideOk[0] && sideOk[1]) {
            if (MathUtils.nextBoolean()) {
                swapNodeParent = swapNode.rightParent;
                leftSide = false;
            }
        } else if (sideOk[1]) {
            swapNodeParent = swapNode.rightParent;
            leftSide = false;
        }
        //Double linked parents
        if (swapNode.leftParent == swapNode.rightParent) {
            arg.singleRemoveChild(startNode, moveChild);
            if (leftSide) {
                swapNode.leftParent = null;
                swapNodeParent.leftChild = null;
            } else {
                swapNode.rightParent = null;
                swapNodeParent.rightChild = null;
            }
            arg.singleAddChild(startNode, swapNode);
            arg.singleAddChild(swapNodeParent, moveChild);
        } else if (swapNode.leftParent == startNode || swapNode.rightParent == startNode) {
            arg.singleRemoveChild(startNode, moveChild);
            if (swapNodeParent.bifurcation) {
                arg.singleRemoveChild(swapNodeParent, swapNode);
                arg.singleAddChild(swapNodeParent, moveChild);
            } else {
                arg.doubleRemoveChild(swapNodeParent, swapNode);
                arg.doubleAddChild(swapNodeParent, moveChild);
            }
            if (startNode.leftChild == null)
                startNode.leftChild = swapNode;
            else
                startNode.rightChild = swapNode;
            if (swapNode.leftParent == null)
                swapNode.leftParent = startNode;
            else
                swapNode.rightParent = startNode;
        } else {
            arg.singleRemoveChild(startNode, moveChild);
            if (swapNodeParent.bifurcation) {
                arg.singleRemoveChild(swapNodeParent, swapNode);
                arg.singleAddChild(swapNodeParent, moveChild);
            } else {
                arg.doubleRemoveChild(swapNodeParent, swapNode);
                arg.doubleAddChild(swapNodeParent, moveChild);
            }
            arg.singleAddChild(startNode, swapNode);
        }
    }
    // TODO Send only changed nodes
    arg.pushTreeChangedEvent();
    assert nodeCheck();
    arg.endTreeEdit();
    try {
        arg.checkTreeIsValid();
    } catch (MutableTree.InvalidTreeException ite) {
        System.out.println(before);
        System.err.println(ite.getMessage());
        System.exit(-1);
    }
    return 0;
}

Example 7 with Node

use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.

the class ARGSwapOperator method reassortmentSwap.

private double reassortmentSwap(NodeRef x) {
    Node startNode = (Node) x;
    Node startChild = startNode.leftChild;
    ArrayList<NodeRef> possibleNodes = new ArrayList<NodeRef>(arg.getNodeCount());
    findNodesAtHeight(possibleNodes, startNode.getHeight());
    assert !possibleNodes.contains(startNode);
    assert possibleNodes.size() > 0;
    Node swapNode = (Node) possibleNodes.get(MathUtils.nextInt(possibleNodes.size()));
    Node swapParent;
    arg.beginTreeEdit();
    if (swapNode.bifurcation) {
        swapParent = swapNode.leftParent;
        arg.doubleRemoveChild(startNode, startChild);
        if (swapParent.bifurcation)
            arg.singleRemoveChild(swapParent, swapNode);
        else
            arg.doubleRemoveChild(swapParent, swapNode);
        arg.doubleAddChild(startNode, swapNode);
        if (startChild.bifurcation) {
            startChild.leftParent = swapParent;
            startChild.rightParent = swapParent;
        } else {
            if (startChild.leftParent == null) {
                startChild.leftParent = swapParent;
            } else {
                startChild.rightParent = swapParent;
            }
        }
        if (!swapParent.bifurcation) {
            swapParent.leftChild = startChild;
            swapParent.rightChild = startChild;
        } else {
            if (swapParent.leftChild == null) {
                swapParent.leftChild = startChild;
            } else {
                swapParent.rightChild = startChild;
            }
        }
    } else {
        boolean leftSide = true;
        boolean[] sideOk = { swapNode.leftParent.getHeight() > startNode.getHeight(), swapNode.rightParent.getHeight() > startNode.getHeight() };
        swapParent = swapNode.leftParent;
        if (sideOk[0] && sideOk[1]) {
            if (MathUtils.nextBoolean()) {
                leftSide = false;
                swapParent = swapNode.rightParent;
            }
        } else if (sideOk[1]) {
            leftSide = false;
            swapParent = swapNode.rightParent;
        }
        if (swapNode.leftParent == swapNode.rightParent) {
            arg.doubleRemoveChild(startNode, startChild);
            if (leftSide) {
                swapParent.leftChild = swapNode.leftParent = null;
                swapParent.leftChild = startChild;
                swapNode.leftParent = startNode;
            } else {
                swapParent.rightChild = swapNode.rightParent = null;
                swapParent.rightChild = startChild;
                swapNode.rightParent = startNode;
            }
            startNode.leftChild = startNode.rightChild = swapNode;
            if (startChild.bifurcation) {
                startChild.leftParent = startChild.rightParent = swapParent;
            } else {
                if (startChild.leftParent == null)
                    startChild.leftParent = swapParent;
                else
                    startChild.rightParent = swapParent;
            }
        } else {
            arg.doubleRemoveChild(startNode, startChild);
            if (swapParent.bifurcation)
                arg.singleRemoveChild(swapParent, swapNode);
            else
                arg.doubleRemoveChild(swapParent, swapNode);
            startNode.leftChild = startNode.rightChild = swapNode;
            if (leftSide)
                swapNode.leftParent = startNode;
            else
                swapNode.rightParent = startNode;
            if (swapParent.bifurcation) {
                if (swapParent.leftChild == null)
                    swapParent.leftChild = startChild;
                else
                    swapParent.rightChild = startChild;
            } else {
                swapParent.leftChild = swapParent.rightChild = startChild;
            }
            if (startChild.bifurcation) {
                startChild.leftParent = startChild.rightParent = swapParent;
            } else {
                if (startChild.leftParent == null)
                    startChild.leftParent = swapParent;
                else
                    startChild.rightParent = swapParent;
            }
        }
    }
    // TODO Limit tree hit
    arg.pushTreeChangedEvent();
    arg.endTreeEdit();
    try {
        arg.checkTreeIsValid();
    } catch (MutableTree.InvalidTreeException ite) {
        System.err.println(ite.getMessage());
        System.exit(-1);
    }
    return 0;
}

Example 8 with Node

use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.

the class ARGSwapOperator method findAllNarrowSwaps.

public int findAllNarrowSwaps(ArrayList<NarrowSwap> moves) {
    for (int i = 0, n = arg.getInternalNodeCount(); i < n; i++) {
        Node x = (Node) arg.getInternalNode(i);
        if (x.bifurcation && !x.isRoot() && x.leftParent.bifurcation) {
            NarrowSwap a = new NarrowSwap(x.leftChild, x, x.leftParent);
            NarrowSwap b = new NarrowSwap(x.rightChild, x, x.leftParent);
            if (a.isValid())
                moves.add(a);
            if (b.isValid())
                moves.add(b);
        }
    }
    return moves.size();
}

Example 9 with Node

use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.

the class ObsoleteARGAddRemoveEventOperator method findPotentialAttachmentSisters.

/*private void checkAllPartitionLabels() {
                   int len = arg.getExternalNodeCount();
                   System.err.println("Tip Partitions:");
                   for(int i=0; i<len; i++) {
                       Node node = (Node)arg.getExternalNode(i);
                       BitSet partSet = node.partitionSet;
                       System.err.print(Tree.Utils.uniqueNewick(arg, node)+":");
                       for(int j=partSet.nextSetBit(0); j>=0; j=partSet.nextSetBit(j+1)) {
                           System.err.print(" "+j);
                       }
                       System.err.println();
                   }
               }*/
/*	private ArrayList<NodeRef> getPotentialTipsToMerge() {
        ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
        int len = arg.getExternalNodeCount();
        for(int i=0; i<len; i++){
            NodeRef nr = arg.getExternalNode(i);
            BitSet pS = ((Node)nr).partitionSet;
            if( pS.cardinality() == 1)
                potentials.add(nr);
        }
        return potentials;
    }

    private ArrayList<Node> getSameTaxonTips(Node asMe) {
        ArrayList<Node> potentials = new ArrayList<Node>();
        int len = arg.getExternalNodeCount();
        Taxon same = asMe.taxon;
        for(int i=0; i<len; i++){
            Node n = (Node)arg.getExternalNode(i);
            if( (n != asMe) && (n.taxon == same) )
                potentials.add(n);
        }
        return potentials;
    }

    private ArrayList<NodeRef> getPotentialTipsToSplit() {
        ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
        int len = arg.getExternalNodeCount();
        for(int i=0; i<len; i++) {
            NodeRef nr = arg.getExternalNode(i);
            BitSet pS = ((Node)nr).partitionSet;
            if( pS.cardinality() > 1 )  // can be split
                potentials.add(nr);
        }
        return potentials;
    }*/
/*;private int drawOnePartition(BitSet partitionSet) {
             int draw = MathUtils.nextInt(partitionSet.cardinality());
     //		System.err.println("draw = "+draw);
             int result = partitionSet.nextSetBit(0);
             for (int i = 0; i < draw; i++)
                 result = partitionSet.nextSetBit(result + 1);
             return result;
         }*/
//    private BitSet bsTransportor = null;
/*
               private ArrayList<NodeRef> getPotentialSubtreesToSplit() {
                   ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
                   // Start with all possible nodes except the root
                   BitSet bsTmp = null;
                   int len = arg.getNodeCount();
                   for(int i=0; i<len; i++) {
                       NodeRef nr = arg.getNode(i);
                       if( ! arg.isRoot(nr) ) {
                           // All tips of nr must have at least two partitions
                           ArrayList<NodeRef> allTips = arg.getDescendantTipNodes(nr);
                           int n = allTips.size();
                           boolean valid = true;
                           for(int j=0; valid && j<n; j++) {
                               BitSet pS = ((Node)allTips.get(j)).partitionSet;
                               if( j == 0 ) {
                                   bsTmp = (BitSet)pS.clone();
                               }
                               else
                                   bsTmp.and(pS);
                               if( pS.cardinality() == 1 )  // can be split
                                   valid = false;
                           }
                           if( valid && (bsTmp.cardinality()>1) )
                               potentials.add(nr);
                       }
                   }
                   bsTransportor = bsTmp;
                   return potentials;
               }
               */
/*	private int countPotentialSubtreesToSplit() {
//		ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
		int count = 0;
		// Start with all possible nodes except the root
		BitSet bsTmp = null;
		int len = arg.getNodeCount();
		for(int i=0; i<len; i++) {
			NodeRef nr = arg.getNode(i);
			if( ! arg.isRoot(nr) ) {
				// All tips of nr must have at least two partitions
				ArrayList<NodeRef> allTips = arg.getDescendantTipNodes(nr);
				int n = allTips.size();
				boolean valid = true;
				for(int j=0; valid && j<n; j++) {
					BitSet pS = ((Node)allTips.get(j)).partitionSet;
					if( j == 0 ) {
						bsTmp = (BitSet)pS.clone();
					}
					else
						bsTmp.and(pS);
					if( pS.cardinality() == 1 )  // can be split
						valid = false;
				}
				if( valid && (bsTmp.cardinality()>1) )
					//potentials.add(nr);
					count++;
			}
		}
		bsTransportor = bsTmp;
		//return potentials;
		return count;
	}*/
//	private ArrayList<NodeRef> getPotentialReattachments(double min) {
//		ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
//		int len = arg.getNodeCount();
//		for(int i=0; i<len; i++) {
//			NodeRef nr = arg.getNode(i);
//			if( !arg.isRoot(nr) && (arg.getMinParentNodeHeight(nr) > min) ) {
//				// Do not reroot and reattach only above current height
//
//				potentials.add(nr);
//			}
//		}
//		return potentials;
//	}
private int findPotentialAttachmentSisters(NodeRef rec, double min, ArrayList<NodeRef> list) {
    int count = 0;
    int len = arg.getNodeCount();
    for (int i = 0; i < len; i++) {
        Node nr = (Node) arg.getNode(i);
        if (!nr.isRoot()) {
            if (arg.getNodeHeight(nr.leftParent) > min) {
                // can add node somewhere between min and nr.parent
                if (list != null)
                    list.add(nr);
                count++;
            }
            if (nr.isReassortment() && arg.getNodeHeight(nr.rightParent) > min) {
                if (list != null)
                    list.add(nr);
                count++;
            }
        }
    }
    return count;
}

Example 10 with Node

use of dr.evomodel.arg.ARGModel.Node in project beast-mcmc by beast-dev.

the class ObsoleteARGNewEventOperator method RemoveOperation.

private double RemoveOperation() throws ARGOperatorFailedException {
    double logq = 0;
    //	    System.err.println("Starting remove ARG operation.");
    //	    System.err.println("Before remove:\n"+arg.toGraphString());
    // 1. Draw reassortment node uniform randomly
    ArrayList<NodeRef> potentialNodes = new ArrayList<NodeRef>();
    int totalPotentials = findPotentialNodesToRemove(potentialNodes);
    if (totalPotentials == 0)
        throw new ARGOperatorFailedException("No reassortment nodes to remove.");
    //	    System.err.println("potentials exist!");
    Node recNode = (Node) potentialNodes.get(MathUtils.nextInt(totalPotentials));
    //        logq += Math.log(totalPotentials);
    double reverseReassortmentHeight = 0;
    double reverseBifurcationHeight = 0;
    double treeHeight = arg.getNodeHeight(arg.getRoot());
    //	    System.err.println("RecNode to remove = "+recNode.number);
    reverseReassortmentHeight = arg.getNodeHeight(recNode);
    arg.beginTreeEdit();
    boolean doneSomething = false;
    Node recParent = recNode.leftParent;
    Node recChild = recNode.leftChild;
    if (recNode.leftParent == recNode.rightParent) {
        // Doubly linked.
        Node recGrandParent = recParent.leftParent;
        //	        reverseBifurcationHeight = arg.getNodeHeight(recParent);
        //            reverseReassortmentSpan = arg.getNodeHeight(recParent) - arg.getNodeHeight(recChild);
        //            reverseBifurcationSpan = arg.getNodeHeight(recGrandParent) - arg.getNodeHeight(recChild);
        //            if (arg.isRoot(recParent)) { // This case should never happen as double links
        //                arg.setRoot(recChild);    // to root can not be added or removed.
        //	            System.err.println("doubled link root?");
        //	            System.exit(-1);
        //            } else {
        //Node recGrandParent = recParent1.leftParent; // And currently recParent must be a bifurcatio
        //System.err.println("recGrand   ="+recGrandParent.number);
        arg.doubleRemoveChild(recGrandParent, recParent);
        arg.doubleRemoveChild(recNode, recChild);
        if (recGrandParent.bifurcation)
            arg.singleAddChild(recGrandParent, recChild);
        else
            arg.doubleAddChild(recGrandParent, recChild);
        //            }
        doneSomething = true;
    // There are not left/right choices to be made for doubly linked removals.
    // End doubly linked.
    } else {
        // Two different parents.
        Node recParent1 = recNode.leftParent;
        Node recParent2 = recNode.rightParent;
        if (!recParent1.bifurcation || recParent1.isRoot()) {
            // One orientation is valid
            recParent1 = recNode.rightParent;
            recParent2 = recNode.leftParent;
        } else if (recParent2.bifurcation && !recParent2.isRoot()) {
            // Both orientations are valid
            if (MathUtils.nextDouble() < 0.5) {
                // choose equally likely
                recParent1 = recNode.rightParent;
                recParent2 = recNode.leftParent;
            }
            logq += Math.log(2);
        }
        // And currently recParent must be a bifurcatio
        Node recGrandParent = recParent1.leftParent;
        Node otherChild = recParent1.leftChild;
        if (otherChild == recNode)
            otherChild = recParent1.rightChild;
        if (recGrandParent.bifurcation)
            arg.singleRemoveChild(recGrandParent, recParent1);
        else
            arg.doubleRemoveChild(recGrandParent, recParent1);
        arg.singleRemoveChild(recParent1, otherChild);
        if (recParent2.bifurcation)
            arg.singleRemoveChild(recParent2, recNode);
        else
            arg.doubleRemoveChild(recParent2, recNode);
        arg.doubleRemoveChild(recNode, recChild);
        if (otherChild != recChild) {
            if (recGrandParent.bifurcation)
                arg.singleAddChild(recGrandParent, otherChild);
            else
                arg.doubleAddChild(recGrandParent, otherChild);
            if (recParent2.bifurcation)
                arg.singleAddChild(recParent2, recChild);
            else
                arg.doubleAddChild(recParent2, recChild);
        } else {
            if (recGrandParent.bifurcation)
                arg.singleAddChildWithOneParent(recGrandParent, otherChild);
            else
                arg.doubleAddChildWithOneParent(recGrandParent, otherChild);
            if (recParent2.bifurcation)
                arg.singleAddChildWithOneParent(recParent2, recChild);
            else
                arg.doubleAddChildWithOneParent(recParent2, recChild);
        }
        //			System.err.println("Sanity check in Remove Operator");
        //			sanityCheck();
        //			System.err.println("End Remove Operator in sanity check");
        doneSomething = true;
        recParent = recParent1;
    }
    reverseBifurcationHeight = arg.getNodeHeight(recParent);
    if (doneSomething) {
        try {
            arg.contractARGWithRecombinant(recParent, recNode, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        } catch (Exception e) {
            System.err.println("here");
            System.err.println(e);
        }
    }
    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 = findPotentialAttachmentPoints(reverseBifurcationHeight,null);
    //	    double d2 = findPotentialAttachmentPoints(reverseReassortmentHeight,null);
    //	    System.err.printf("d1 = %5.4f, d2 = %5.4f\n",d1,d2);
    //	    logq -= Math.log(findPotentialAttachmentPoints(reverseBifurcationHeight,null));
    //	    logq -= Math.log(findPotentialAttachmentPoints(reverseReassortmentHeight,null));
    //        int nodes = arg.getInternalNodeCount() - 1;
    //
    //        logq += 3*Math.log(2);
    //
    //        logq = 100;
    //	    System.err.println("logq remove = "+logq);
    logq = 0;
    // 1 / total potentials * 1 / 2 (if valid) * length1 * length2 * attachmentSisters
    return logq;
}