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

Example 16 with Node

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

the class ARGAddRemoveEventOperator method findPotentialAttachmentPoints.

private int findPotentialAttachmentPoints(double time, ArrayList<NodeRef> list) {
    int count = 0;
    for (int i = 0, n = arg.getNodeCount(); i < n; i++) {
        NodeRef nr = arg.getNode(i);
        if (!arg.isRoot(nr)) {
            if (arg.getNodeHeight(nr) < time) {
                Node left = (Node) arg.getParent(nr, 0);
                Node right = (Node) arg.getParent(nr, 1);
                if (arg.isBifurcation(nr)) {
                    assert left == right;
                    if (arg.getNodeHeight(left) > time) {
                        if (list != null)
                            list.add(nr);
                        count++;
                    }
                } else {
                    if (arg.getNodeHeight(left) > time) {
                        if (list != null)
                            list.add(nr);
                        count++;
                    }
                    if (arg.getNodeHeight(right) > time) {
                        if (list != null)
                            list.add(nr);
                        count++;
                    }
                }
            }
        } else {
            if (arg.getNodeHeight(nr) < time) {
                if (list != null)
                    list.add(nr);
                count++;
            }
        }
    }
    return count;
}

Example 17 with Node

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

the class ARGAddRemoveEventOperator method RemoveOperation.

private double RemoveOperation() throws NoReassortmentEventException {
    double logHastings = 0;
    // 1. Draw reassortment node uniform randomly
    ArrayList<NodeRef> potentialNodes = new ArrayList<NodeRef>();
    int totalPotentials = findPotentialNodesToRemove(potentialNodes);
    if (totalPotentials == 0)
        throw new NoReassortmentEventException();
    //		logHastings += Math.log((double)arg.getReassortmentNodeCount());
    logHastings += Math.log((double) totalPotentials);
    //		double diff =(double)arg.getReassortmentNodeCount() - totalPotentials;
    //
    //		if(MathUtils.nextDouble() < diff/totalPotentials)
    //			throw new NoReassortmentEventException();
    Node recNode = (Node) potentialNodes.get(MathUtils.nextInt(totalPotentials));
    double[] removePartitioningValues = recNode.partitioning.getParameterValues();
    double beforeReassortmentHeight = recNode.getHeight();
    double beforeBifurcationHeight = 0;
    double beforeTreeHeight = arg.getNodeHeight(arg.getRoot());
    arg.beginTreeEdit();
    boolean doneSomething = false;
    Node recParent = recNode.leftParent;
    Node recChild = recNode.leftChild;
    Node beforeReassortChild = recNode.leftChild;
    Node beforeBifurcationChild = recNode.leftParent;
    if (recNode.leftParent == recNode.rightParent) {
        if (!arg.isRoot(recNode.leftParent)) {
            beforeBifurcationHeight = recParent.getHeight();
            Node recGrandParent = recParent.leftParent;
            arg.doubleRemoveChild(recGrandParent, recParent);
            arg.doubleRemoveChild(recNode, recChild);
            if (recGrandParent.bifurcation)
                arg.singleAddChild(recGrandParent, recChild);
            else
                arg.doubleAddChild(recGrandParent, recChild);
            doneSomething = true;
            beforeBifurcationChild = beforeReassortChild;
        } else {
            assert recChild.bifurcation;
            assert false;
        }
        logHastings += LOG_TWO;
    } else {
        Node recDeleteParent = null;
        Node recKeepParent = null;
        if (recNode.leftParent.bifurcation && recNode.rightParent.bifurcation) {
            if (MathUtils.nextBoolean()) {
                recDeleteParent = recNode.rightParent;
                recKeepParent = recNode.leftParent;
            } else {
                recDeleteParent = recNode.leftParent;
                recKeepParent = recNode.rightParent;
            }
            logHastings += LOG_TWO;
        } else if (recNode.rightParent.bifurcation) {
            recDeleteParent = recNode.rightParent;
            recKeepParent = recNode.leftParent;
        } else {
            recDeleteParent = recNode.leftParent;
            recKeepParent = recNode.rightParent;
        }
        beforeBifurcationChild = recDeleteParent.leftChild;
        if (beforeBifurcationChild == recNode) {
            beforeBifurcationChild = recDeleteParent.rightChild;
        }
        beforeBifurcationHeight = recDeleteParent.getHeight();
        if (arg.isRoot(recDeleteParent)) {
            Node oldRoot = (Node) arg.getOtherChild(recDeleteParent, recNode);
            Node oldRootLeft = oldRoot.leftChild;
            Node oldRootRight = oldRoot.rightChild;
            if (oldRoot == recKeepParent) {
                arg.singleRemoveChild(recDeleteParent, recNode);
                arg.singleRemoveChild(recDeleteParent, oldRoot);
                arg.singleRemoveChild(oldRoot, oldRootLeft);
                arg.singleRemoveChild(oldRoot, oldRootRight);
                arg.singleAddChild(recDeleteParent, oldRootLeft);
                arg.singleAddChild(recDeleteParent, oldRootRight);
                arg.singleRemoveChild(recDeleteParent, recNode);
                arg.doubleRemoveChild(recNode, recChild);
                arg.singleAddChild(recDeleteParent, recChild);
                recDeleteParent.setHeight(oldRoot.getHeight());
                recDeleteParent = oldRoot;
            } else {
                arg.singleRemoveChild(recDeleteParent, recNode);
                arg.singleRemoveChild(recDeleteParent, oldRoot);
                arg.singleRemoveChild(oldRoot, oldRootLeft);
                arg.singleRemoveChild(oldRoot, oldRootRight);
                arg.singleAddChild(recDeleteParent, oldRootLeft);
                arg.singleAddChild(recDeleteParent, oldRootRight);
                if (recKeepParent.bifurcation)
                    arg.singleRemoveChild(recKeepParent, recNode);
                else
                    arg.doubleRemoveChild(recKeepParent, recNode);
                arg.doubleRemoveChild(recNode, recChild);
                if (recKeepParent.bifurcation)
                    arg.singleAddChild(recKeepParent, recChild);
                else
                    arg.doubleAddChild(recKeepParent, recChild);
                recDeleteParent.setHeight(oldRoot.getHeight());
                recDeleteParent = oldRoot;
            }
        } else {
            Node recGrandParent = recDeleteParent.leftParent;
            Node otherChild = recDeleteParent.leftChild;
            if (otherChild == recNode)
                otherChild = recDeleteParent.rightChild;
            if (recGrandParent.bifurcation)
                arg.singleRemoveChild(recGrandParent, recDeleteParent);
            else
                arg.doubleRemoveChild(recGrandParent, recDeleteParent);
            arg.singleRemoveChild(recDeleteParent, otherChild);
            if (recKeepParent.bifurcation)
                arg.singleRemoveChild(recKeepParent, recNode);
            else
                arg.doubleRemoveChild(recKeepParent, recNode);
            arg.doubleRemoveChild(recNode, recChild);
            if (otherChild != recChild) {
                if (recGrandParent.bifurcation)
                    arg.singleAddChild(recGrandParent, otherChild);
                else
                    arg.doubleAddChild(recGrandParent, otherChild);
                if (recKeepParent.bifurcation)
                    arg.singleAddChild(recKeepParent, recChild);
                else
                    arg.doubleAddChild(recKeepParent, recChild);
            } else {
                if (recGrandParent.bifurcation)
                    arg.singleAddChildWithOneParent(recGrandParent, otherChild);
                else
                    arg.doubleAddChildWithOneParent(recGrandParent, otherChild);
                if (recKeepParent.bifurcation)
                    arg.singleAddChildWithOneParent(recKeepParent, recChild);
                else
                    arg.doubleAddChildWithOneParent(recKeepParent, recChild);
            }
        }
        doneSomething = true;
        recParent = recDeleteParent;
    }
    if (relaxed) {
        double[] rateValues = recParent.rateParameter.getParameterValues();
        logHastings -= ratePrior.getAddHastingsRatio(rateValues);
    }
    if (doneSomething) {
        try {
            arg.contractARG(recParent, recNode, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        //                arg.contractARGWithRecombinant(recParent, recNode,
        //                        internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        } catch (Exception e) {
            System.err.println(e.getMessage());
            System.err.println(e);
            System.exit(-1);
        }
    }
    int max = Math.max(recParent.getNumber(), recNode.getNumber());
    int min = Math.min(recParent.getNumber(), recNode.getNumber());
    for (int i = 0, n = arg.getNodeCount(); i < n; i++) {
        Node x = (Node) arg.getNode(i);
        if (x.getNumber() > max) {
            x.number--;
        }
        if (x.getNumber() > min) {
            x.number--;
        }
    }
    adjustRandomPartitioning();
    arg.pushTreeSizeDecreasedEvent();
    arg.endTreeEdit();
    try {
        arg.checkTreeIsValid();
    } catch (MutableTree.InvalidTreeException ite) {
        throw new RuntimeException(ite.toString() + "\n" + arg.toString() + "\n" + TreeUtils.uniqueNewick(arg, arg.getRoot()));
    }
    assert nodeCheck() : arg.toARGSummary();
    //Do the backwards stuff now :(
    double afterTreeHeight = arg.getNodeHeight(arg.getRoot());
    //		This is the ugly mixture proposal
    //		double meanRoot = 4.0 / afterTreeHeight;
    //		double case1 = 0.95;
    //
    //		if(beforeBifurcationHeight < afterTreeHeight){
    //			logHastings -= 2.0*Math.log(afterTreeHeight) - Math.log(2.0*case1);
    //		}else{
    //			double additional = beforeBifurcationHeight - afterTreeHeight;
    //			logHastings -= Math.log(afterTreeHeight) + additional*meanRoot -
    //				Math.log((1-case1)*meanRoot);
    //		}
    double theta = probBelowRoot / afterTreeHeight;
    logHastings -= theta * (beforeBifurcationHeight + beforeReassortmentHeight) - LOG_TWO - 2.0 * Math.log(theta) + Math.log(1 - Math.exp(-2.0 * afterTreeHeight * theta));
    logHastings -= Math.log((double) findPotentialAttachmentPoints(beforeBifurcationHeight, null) * findPotentialAttachmentPoints(beforeReassortmentHeight, null));
    logHastings -= getPartitionAddHastingsRatio(removePartitioningValues);
    logHastings -= LOG_TWO;
    //        }
    assert nodeCheck();
    assert !Double.isNaN(logHastings) && !Double.isInfinite(logHastings);
    return logHastings;
}

Example 18 with Node

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

the class ARGAddRemoveEventOperator method sanityCheck.

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 19 with Node

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

the class ARGAddRemoveEventOperator method findPotentialNodesToRemove.

private int findPotentialNodesToRemove(ArrayList<NodeRef> list) {
    int count = 0;
    int n = arg.getNodeCount();
    for (int i = 0; i < n; i++) {
        Node node = (Node) arg.getNode(i);
        Node lp = node.leftParent;
        Node rp = node.rightParent;
        if (node.isReassortment() && (lp.bifurcation || rp.bifurcation)) {
            if (list != null)
                list.add(node);
            count++;
        }
    }
    return count;
}

Example 20 with Node

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

the class ARGRelaxedClock method getBranchRate.

public double getBranchRate(Tree tree, NodeRef nodeRef) {
    Node treeNode = (Node) nodeRef;
    Node argNode = (Node) treeNode.mirrorNode;
    return globalRateParameter.getParameterValue(0) * argNode.getRate(partition);
}