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

Example 11 with Node

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

the class ObsoleteARGAddRemoveEventOperator method AddOperation.

private double AddOperation() throws ARGOperatorFailedException {
    double logq = 0;
    //	    System.err.println("Starting add ARG operation. ... "+arg.getReassortmentNodeCount());
    // Draw potential places to add a reassortment node
    ArrayList<NodeRef> potentialNodes = new ArrayList<NodeRef>();
    int totalPotentials = findPotentialReassortmentNodes(potentialNodes);
    if (totalPotentials == 0) {
        System.err.println("Should never get here AA");
        System.exit(-1);
        throw new ARGOperatorFailedException("No more nodes to recombine.");
    }
    Node recNode = (Node) potentialNodes.get(MathUtils.nextInt(totalPotentials));
    Node recParentL = recNode.leftParent;
    Node recParentR = recNode.rightParent;
    Node recParent = recParentL;
    if ((recParentL != recParentR) && MathUtils.nextDouble() > 0.5)
        recParent = recParentR;
    logq += Math.log(totalPotentials);
    // Attachment must occur
    double minHeight = arg.getNodeHeight(recNode);
    // previous to this time
    ArrayList<NodeRef> attachments = new ArrayList<NodeRef>();
    int totalAttachments = findPotentialAttachmentSisters(recNode, minHeight, attachments);
    if (totalAttachments == 0) {
        throw new ARGOperatorFailedException("no more attachment points for this recomb");
    }
    Node sisNode = (Node) attachments.get(MathUtils.nextInt(totalAttachments));
    Node sisParentL = sisNode.leftParent;
    Node sisParentR = sisNode.rightParent;
    Node sisParent = sisParentL;
    if (sisParentL != sisParentR) {
        if (arg.getNodeHeight(sisParentL) <= minHeight)
            sisParent = sisParentR;
        else if (arg.getNodeHeight(sisParentR) > minHeight && MathUtils.nextDouble() > 0.5)
            sisParent = sisParentR;
    }
    logq += Math.log(totalAttachments);
    Node newBifurcation = arg.new Node();
    Node newReassortment = arg.new Node();
    newReassortment.bifurcation = false;
    double sisHeight = sisNode.getHeight();
    if (sisHeight < minHeight) {
        sisHeight = minHeight;
    }
    double spHeight = arg.getNodeHeight(sisParent);
    double totalLength = spHeight - sisHeight;
    double newLength = sisHeight + MathUtils.nextDouble() * totalLength;
    // prior ratio
    logq -= Math.log(totalLength);
    newBifurcation.heightParameter = new Parameter.Default(newLength);
    newBifurcation.setupHeightBounds();
    // Uniform[spHeight, sisHeight]
    logq += Math.log(totalLength);
    double topHeight = newLength;
    double recParentHeight = arg.getNodeHeight(recParent);
    if (topHeight > recParentHeight)
        topHeight = recParentHeight;
    double recHeight = arg.getNodeHeight(recNode);
    totalLength = topHeight - recHeight;
    newLength = recHeight + MathUtils.nextDouble() * totalLength;
    // prior ratio
    logq -= Math.log(totalLength);
    newReassortment.heightParameter = new Parameter.Default(newLength);
    newReassortment.setupHeightBounds();
    // Uniform[bifurcationHeight,recNodeHeight]
    logq += Math.log(totalLength);
    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);
    /*		BitSet bitLeft = new BitSet();
		BitSet bitRight = new BitSet();
		if (MathUtils.nextDouble() < 0.5) {
			bitLeft.set(0);
			bitRight.set(1);
		} else {
			bitLeft.set(1);
			bitRight.set(0);
		}
		logq += Math.log(2.0);*/
    Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    if (sisNode != recNode) {
        arg.addChildAsRecombinant(newBifurcation, recParent, newReassortment, partitioning);
    } else {
        arg.addChildAsRecombinant(newBifurcation, newBifurcation, newReassortment, partitioning);
    }
    /* if (sisNode != recNode) {
                        arg.addChildAsRecombinant(newBifurcation, recParent,
                                newReassortment, bitLeft, bitRight);
                    } else {
                        arg.addChildAsRecombinant(newBifurcation, newBifurcation,
                                newReassortment, bitLeft, bitRight);
                    }*/
    arg.expandARGWithRecombinant(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
    // arg.addNewHeightParameter(newReassortment.heightParameter,internalNodeParameters);
    // arg.expandNodesWithRecombinant(arg.getRoot(), recNR);
    // arg.expandNodesWithSubtree(arg.getRoot(), newSubtree);
    // arg.reconstructTrees();
    //
    //		System.err.println("End ARG\n" + arg.toGraphString()); System.err.println("recNode   = " + recNode.number);
    //		System.err.println("recParent = " + recParent.number);
    //		System.err.println("sisNode   = " + sisNode.number);
    //		System.err.println("sisParent = " + sisParent.number);
    //System.exit(-1);
    //		ARGTree tree = new ARGTree(arg);
    // System.err.println("Reassortment nodes =
    // "+arg.getReassortmentNodeCount());
    // //System.exit(-1);
    // //System.err.println("ARGTree = "+Tree.Utils.uniqueNewick(tree,
    // tree.getRoot()));
    //		tree = new ARGTree(arg, 0);
    //		System.err.println("Part 0  = " + tree.toString());
    //		tree = new ARGTree(arg, 1);
    //		System.err.println("Part 1  = " + tree.toString());
    // //if( recNode == sisNode)
    //System.exit(-1);
    //    	System.err.println("Sanity check in Add Operation");
    //		sanityCheck();
    //		System.err.println("End Add Operator sanity check");
    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()));
    //        }
    //	    System.err.println("Checking add validity.");
    // todo -- check all ARGTree.Roots
    //	    if (!arg.validRoot())
    //	        throw new OperatorFailedException("Roots are invalid");
    // System.err.println("Made it thru once!");
    // System.exit(-1);
    // times++;
    // System.err.println("Adds = "+times);
    // if( times >= 4 )
    // throw new OperatorFailedException("Do many tries!");
    // System.err.println("Add a recomb node.");
    //		int cnt = arg.getReassortmentNodeCount();
    //		if (cnt > 20)
    //			throw new OperatorFailedException("No more than X reassortments");
    logq -= Math.log(findCurrentReassortmentNodes(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;
    // TODO move into prior
    logq += lnGamma(nodes) - lnGamma(nodes - 2);
    logq = 0;
    return logq;
}

Example 12 with Node

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

the class ObsoleteARGAddRemoveEventOperator method findPotentialReassortmentNodes.

/*
    private ArrayList<NodeRef> getPotentialSubtreesToMerge() {
        ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
        int n = arg.getNodeCount();
        for(int i=0; i<n; i++) {
            Node node = (Node)arg.getNode(i);
            if( node.dupSister != null ) {
                potentials.add(node);
            }
        }
        return potentials;
    }

    private int countPotentialSubtreesToMerge() {
        int count = 0;
        //ArrayList<NodeRef> potentials = new ArrayList<NodeRef>();
        int n = arg.getNodeCount();
        for(int i=0; i<n; i++) {
            Node node = (Node)arg.getNode(i);
            if( node.dupSister != null ) {
                //potentials.add(node);
                count++;
            }
        }
        //return potentials;
        return count;
    }
    */
/*
                * Generates a list of all branches in the ARG, there each branch leads up the ARG from the saved nodes
                *
                * Should return a deterministic function of # bifurcations and # reassortments
                *
                */
private int findPotentialReassortmentNodes(ArrayList<NodeRef> list) {
    int count = 0;
    int n = arg.getNodeCount();
    for (int i = 0; i < n; i++) {
        Node node = (Node) arg.getNode(i);
        if (node.leftParent != null) {
            // rules out only the root; may subclass for different types of problems
            list.add(node);
            count++;
            if (node.isReassortment()) {
                list.add(node);
                // reassortment nodes have two parent branches
                count++;
            }
        }
    }
    //		}
    return count;
}

Example 13 with Node

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

the class ObsoleteARGAddRemoveEventOperator 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();
        }
    }
}

Example 14 with Node

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

the class ObsoleteARGAddRemoveEventOperator method findCurrentReassortmentNodes.

private int findCurrentReassortmentNodes(ArrayList<NodeRef> list) {
    int count = 0;
    int n = arg.getNodeCount();
    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.rightParent != root)) {
            if (list != null)
                list.add(node);
            count++;
        }
    }
    return count;
}

Example 15 with Node

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

the class ObsoleteARGAddRemoveEventOperator method RemoveOperation.

/*private double RemoveOperation() throws OperatorFailedException {
             arg.beginTreeEdit();
             if (arg.getReassortmentNodeCount() > 0)
                 arg.removeNullCounter();
              try {
                         arg.endTreeEdit();
                     } catch (MutableTree.InvalidTreeException ite) {
                         throw new RuntimeException(ite.toString() + "\n" + arg.toString()
                                 + "\n" + Tree.Utils.uniqueNewick(arg, arg.getRoot()));
                     }
             arg.pushTreeSizeChangedEvent();
             return 0;
         }*/
/*   private double AddOperation() throws OperatorFailedException {
            arg.beginTreeEdit();
            arg.addNullCounter();
            try {
                        arg.endTreeEdit();
                    } catch (MutableTree.InvalidTreeException ite) {
                        throw new RuntimeException(ite.toString() + "\n" + arg.toString()
                                + "\n" + Tree.Utils.uniqueNewick(arg, arg.getRoot()));
                    }
            arg.pushTreeSizeChangedEvent();
            return 0;
        }
    */
private double RemoveOperation() throws ARGOperatorFailedException {
    double logq = 0;
    //	    System.err.println("Starting remove ARG operation.");
    // 1. Draw reassortment node uniform randomly
    ArrayList<NodeRef> potentialNodes = new ArrayList<NodeRef>();
    int totalPotentials = findCurrentReassortmentNodes(potentialNodes);
    if (totalPotentials == 0)
        throw new ARGOperatorFailedException("No reassortment nodes to remove.");
    Node recNode = (Node) potentialNodes.get(MathUtils.nextInt(totalPotentials));
    logq += Math.log(totalPotentials);
    double reverseReassortmentSpan = 0;
    double reverseBifurcationSpan = 0;
    arg.beginTreeEdit();
    boolean doneSomething = false;
    Node recParent = recNode.leftParent;
    Node recChild = recNode.leftChild;
    if (recNode.leftParent == recNode.rightParent) {
        // Doubly linked.
        Node recGrandParent = recParent.leftParent;
        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
            // to root can not be added or removed.
            arg.setRoot(recChild);
        } 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 && !recParent2.bifurcation) || (!recParent1.bifurcation && recParent2.isRoot()) || (!recParent2.bifurcation && recParent1.isRoot())) {
            arg.endTreeEdit();
            //                }
            throw new ARGOperatorFailedException("Not reversible deletion.");
        }
        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);
        }
        //			System.err.println("recNode    ="+recNode.number);
        //			System.err.println("recParent1 ="+recParent1.number);
        //			System.err.println("recParent2 ="+recParent2.number);
        // And currently recParent must be a bifurcatio
        Node recGrandParent = recParent1.leftParent;
        //			System.err.println("recGrand   ="+recGrandParent.number);
        //Node recChild = recNode.leftChild;
        Node otherChild = recParent1.leftChild;
        if (otherChild == recNode)
            otherChild = recParent1.rightChild;
        //			System.err.println("recChild   ="+recChild.number);
        //			System.err.println("otherChild ="+otherChild.number);
        reverseReassortmentSpan = Math.min(arg.getNodeHeight(recParent1), arg.getNodeHeight(recParent2)) - arg.getNodeHeight(recChild);
        reverseBifurcationSpan = arg.getNodeHeight(recGrandParent) - Math.max(arg.getNodeHeight(recChild), arg.getNodeHeight(otherChild));
        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;
        if ((recChild.getHeight() > recParent2.getHeight()) || (otherChild.getHeight() > recGrandParent.getHeight())) {
            arg.endTreeEdit();
            throw new RuntimeException("How did I get here?");
        }
        recParent = recParent1;
    }
    if (doneSomething) {
        //			System.err.println("Trying to remove "+recNode.number+" and "+recParent.number);
        //
        arg.contractARGWithRecombinant(recParent, recNode, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
    }
    //		System.err.println("End ARG\n"+arg.toGraphString());
    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()));
    //        }
    //	    System.err.println("Checking remove validity.");
    // todo -- check all ARGTree.Roots
    //	    if (!arg.validRoot())
    //	        throw new OperatorFailedException("Roots are invalid");
    //         logq -= Math.log(reverseBifurcationSpan) + Math.log(reverseReassortmentSpan);   // TODO removed? because of prior ratio
    // findPotentialRessortmentNodes()
    logq -= Math.log(arg.getNodeCount() + arg.getReassortmentNodeCount() - 1);
    logq -= Math.log(this.findPotentialAttachmentSisters(recChild, arg.getNodeHeight(recChild), null));
    //        System.err.println(drawRandomPartitioning(null));
    // logq -= drawRandomPartitioning(null);
    //	    System.err.println("End remove ARG operation.");
    int nodes = arg.getInternalNodeCount() - 1;
    // TODO move into prior
    logq += lnGamma(nodes) - lnGamma(nodes + 2);
    logq += 3 * Math.log(2);
    logq = 0;
    // 1 / total potentials * 1 / 2 (if valid) * length1 * length2 * attachmentSisters
    return logq;
}