Examples with NodeRef dr.evolution.tree.NodeRef used on opensource projects

Example 1 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class CheckPointTreeModifier method incorporateAdditionalTaxa.

/**
     * Add the remaining taxa, which can be identified through the TreeDataLikelihood XML elements.
     */
public ArrayList<NodeRef> incorporateAdditionalTaxa(CheckPointUpdaterApp.UpdateChoice choice, BranchRates rateModel) {
    System.out.println("Tree before adding taxa:\n" + treeModel.toString() + "\n");
    ArrayList<NodeRef> newTaxaNodes = new ArrayList<NodeRef>();
    for (String str : newTaxaNames) {
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId().equals(str)) {
                newTaxaNodes.add(treeModel.getExternalNode(i));
                //always take into account Taxon dates vs. dates set through a TreeModel
                System.out.println(treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " with height " + treeModel.getNodeHeight(treeModel.getExternalNode(i)) + " or " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getHeight());
            }
        }
    }
    System.out.println("newTaxaNodes length = " + newTaxaNodes.size());
    ArrayList<Taxon> currentTaxa = new ArrayList<Taxon>();
    for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
        boolean taxonFound = false;
        for (String str : newTaxaNames) {
            if (str.equals((treeModel.getNodeTaxon(treeModel.getExternalNode(i))).getId())) {
                taxonFound = true;
            }
        }
        if (!taxonFound) {
            System.out.println("Adding " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " to list of current taxa");
            currentTaxa.add(treeModel.getNodeTaxon(treeModel.getExternalNode(i)));
        }
    }
    System.out.println("Current taxa count = " + currentTaxa.size());
    //iterate over both current taxa and to be added taxa
    boolean originTaxon = true;
    for (Taxon taxon : currentTaxa) {
        if (taxon.getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("Current taxon " + taxon.getId() + " has node height 0.0");
        }
    }
    for (NodeRef newTaxon : newTaxaNodes) {
        if (treeModel.getNodeTaxon(newTaxon).getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("New taxon " + treeModel.getNodeTaxon(newTaxon).getId() + " has node height 0.0");
        }
    }
    //check the Tree(Data)Likelihoods in the connected set of likelihoods
    //focus on TreeDataLikelihood, which has getTree() to get the tree for each likelihood
    //also get the DataLikelihoodDelegate from TreeDataLikelihood
    ArrayList<TreeDataLikelihood> likelihoods = new ArrayList<TreeDataLikelihood>();
    ArrayList<Tree> trees = new ArrayList<Tree>();
    ArrayList<DataLikelihoodDelegate> delegates = new ArrayList<DataLikelihoodDelegate>();
    for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
        if (likelihood instanceof TreeDataLikelihood) {
            likelihoods.add((TreeDataLikelihood) likelihood);
            trees.add(((TreeDataLikelihood) likelihood).getTree());
            delegates.add(((TreeDataLikelihood) likelihood).getDataLikelihoodDelegate());
        }
    }
    //suggested to go through TreeDataLikelihoodParser and give it an extra option to create a HashMap
    //keyed by the tree; am currently not overly fond of this approach
    ArrayList<PatternList> patternLists = new ArrayList<PatternList>();
    for (DataLikelihoodDelegate del : delegates) {
        if (del instanceof BeagleDataLikelihoodDelegate) {
            patternLists.add(((BeagleDataLikelihoodDelegate) del).getPatternList());
        } else if (del instanceof MultiPartitionDataLikelihoodDelegate) {
            MultiPartitionDataLikelihoodDelegate mpdld = (MultiPartitionDataLikelihoodDelegate) del;
            List<PatternList> list = mpdld.getPatternLists();
            for (PatternList pList : list) {
                patternLists.add(pList);
            }
        }
    }
    if (patternLists.size() == 0) {
        throw new RuntimeException("No patterns detected. Please make sure the XML file is BEAST 1.9 compatible.");
    }
    //aggregate all patterns to create distance matrix
    //TODO What about different trees for different partitions?
    Patterns patterns = new Patterns(patternLists.get(0));
    if (patternLists.size() > 1) {
        for (int i = 1; i < patternLists.size(); i++) {
            patterns.addPatterns(patternLists.get(i));
        }
    }
    //set the patterns for the distance matrix computations
    choice.setPatterns(patterns);
    //add new taxa one at a time
    System.out.println("Adding " + newTaxaNodes.size() + " taxa ...");
    for (NodeRef newTaxon : newTaxaNodes) {
        treeModel.setNodeHeight(newTaxon, treeModel.getNodeTaxon(newTaxon).getHeight());
        System.out.println("\nadding Taxon: " + newTaxon + " (height = " + treeModel.getNodeHeight(newTaxon) + ")");
        //check if this taxon has a more recent sampling date than all other nodes in the current TreeModel
        double offset = checkCurrentTreeNodes(newTaxon, treeModel.getRoot());
        System.out.println("Sampling date offset when adding " + newTaxon + " = " + offset);
        //AND set its current node height to 0.0 IF no originTaxon has been found
        if (offset < 0.0) {
            if (!originTaxon) {
                System.out.println("Updating all node heights with offset " + Math.abs(offset));
                updateAllTreeNodes(Math.abs(offset), treeModel.getRoot());
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        } else if (offset == 0.0) {
            if (!originTaxon) {
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        }
        //get the closest Taxon to the Taxon that needs to be added
        //take into account which taxa can currently be chosen
        Taxon closest = choice.getClosestTaxon(treeModel.getNodeTaxon(newTaxon), currentTaxa);
        System.out.println("\nclosest Taxon: " + closest + " with original height: " + closest.getHeight());
        //get the distance between these two taxa
        double distance = choice.getDistance(treeModel.getNodeTaxon(newTaxon), closest);
        System.out.println("at distance: " + distance);
        //TODO what if distance == 0.0 ??? how to choose closest taxon then (in absence of geo info)?
        //find the NodeRef for the closest Taxon (do not rely on node numbering)
        NodeRef closestRef = null;
        //careful with node numbering and subtract number of new taxa
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)) == closest) {
                closestRef = treeModel.getExternalNode(i);
            }
        }
        System.out.println(closestRef + " with height " + treeModel.getNodeHeight(closestRef));
        //System.out.println("trying to set node height: " + closestRef + " from " + treeModel.getNodeHeight(closestRef) + " to " + closest.getHeight());
        //treeModel.setNodeHeight(closestRef, closest.getHeight());
        double timeForDistance = distance / rateModel.getBranchRate(treeModel, closestRef);
        System.out.println("timeForDistance = " + timeForDistance);
        //get parent node of branch that will be split
        NodeRef parent = treeModel.getParent(closestRef);
        //determine height of new node
        double insertHeight;
        if (treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon)) {
            insertHeight = treeModel.getNodeHeight(closestRef) + timeForDistance / 2.0;
            System.out.println("treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon): " + insertHeight);
            if (insertHeight >= treeModel.getNodeHeight(parent)) {
                insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
            }
        } else {
            double remainder = (timeForDistance - Math.abs(treeModel.getNodeHeight(closestRef) - treeModel.getNodeHeight(newTaxon))) / 2.0;
            if (remainder > 0) {
                insertHeight = Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)) + remainder;
                System.out.println("remainder > 0: " + insertHeight);
                if (insertHeight >= treeModel.getNodeHeight(parent)) {
                    insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
                }
            } else {
                insertHeight = EPSILON * (treeModel.getNodeHeight(parent) - Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)));
                insertHeight += Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon));
                System.out.println("remainder <= 0: " + insertHeight);
            }
        }
        System.out.println("insert at height: " + insertHeight);
        //pass on all the necessary variables to a method that adds the new taxon to the tree
        addTaxonAlongBranch(newTaxon, parent, closestRef, insertHeight);
        //option to print tree after each taxon addition
        System.out.println("\nTree after adding taxon " + newTaxon + ":\n" + treeModel.toString());
        //add newly added Taxon to list of current taxa
        currentTaxa.add(treeModel.getNodeTaxon(newTaxon));
    }
    return newTaxaNodes;
}

Example 2 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ARGAddRemoveEventOperator method AddOperation.

private double AddOperation() {
    double logHastings = 0;
    double treeHeight = arg.getNodeHeight(arg.getRoot());
    double newBifurcationHeight = Double.POSITIVE_INFINITY;
    double newReassortmentHeight = Double.POSITIVE_INFINITY;
    double theta = probBelowRoot / treeHeight;
    while (newBifurcationHeight > treeHeight && newReassortmentHeight > treeHeight) {
        newBifurcationHeight = MathUtils.nextExponential(theta);
        newReassortmentHeight = MathUtils.nextExponential(theta);
    }
    logHastings += theta * (newBifurcationHeight + newReassortmentHeight) - LOG_TWO - 2.0 * Math.log(theta) + Math.log(1 - Math.exp(-2.0 * treeHeight * theta));
    if (newBifurcationHeight < newReassortmentHeight) {
        double temp = newBifurcationHeight;
        newBifurcationHeight = newReassortmentHeight;
        newReassortmentHeight = temp;
    }
    //2. Find the possible re-assortment and bifurcation points.
    ArrayList<NodeRef> potentialBifurcationChildren = new ArrayList<NodeRef>();
    ArrayList<NodeRef> potentialReassortmentChildren = new ArrayList<NodeRef>();
    int totalPotentialBifurcationChildren = findPotentialAttachmentPoints(newBifurcationHeight, potentialBifurcationChildren);
    int totalPotentialReassortmentChildren = findPotentialAttachmentPoints(newReassortmentHeight, potentialReassortmentChildren);
    assert totalPotentialBifurcationChildren > 0;
    assert totalPotentialReassortmentChildren > 0;
    logHastings += Math.log((double) potentialBifurcationChildren.size() * potentialReassortmentChildren.size());
    //3.  Choose your re-assortment location.
    Node reassortChild = (Node) potentialReassortmentChildren.get(MathUtils.nextInt(totalPotentialReassortmentChildren));
    Node reassortLeftParent = reassortChild.leftParent;
    Node reassortRightParent = reassortChild.rightParent;
    Node reassortSplitParent = reassortChild.leftParent;
    boolean splitReassortLeftParent = true;
    if (!reassortChild.bifurcation) {
        boolean[] tester = { arg.getNodeHeight(reassortLeftParent) > newReassortmentHeight, arg.getNodeHeight(reassortRightParent) > newReassortmentHeight };
        if (tester[0] && tester[1]) {
            if (MathUtils.nextBoolean()) {
                splitReassortLeftParent = false;
                reassortSplitParent = reassortRightParent;
            }
        //        		 logHastings += LOG_TWO;
        } else if (tester[0]) {
        //nothing to do, setup above
        } else if (tester[1]) {
            splitReassortLeftParent = false;
            reassortSplitParent = reassortRightParent;
        } else {
            assert false;
        }
    }
    //        if (recParentL != recParentR) {
    //            boolean[] tester = {arg.getNodeHeight(recParentL) > newReassortmentHeight,
    //                    arg.getNodeHeight(recParentR) > newReassortmentHeight};
    //
    //            if (tester[0] && tester[1]) {
    //                if (MathUtils.nextBoolean()) {
    //                    recParent = recParentR;
    //                }
    //
    //                logHastings += LOG_TWO;
    //            } else if (tester[0]) {
    //                recParent = recParentL;
    //            } else {
    //                recParent = recParentR;
    //            }
    //        }
    //4. Choose your bifurcation location.
    Node bifurcationChild = (Node) potentialBifurcationChildren.get(MathUtils.nextInt(potentialBifurcationChildren.size()));
    Node bifurcationLeftParent = bifurcationChild.leftParent;
    Node bifurcationRightParent = bifurcationChild.rightParent;
    Node bifurcationSplitParent = bifurcationLeftParent;
    boolean splitBifurcationLeftParent = true;
    if (!bifurcationChild.bifurcation) {
        boolean[] tester = { arg.getNodeHeight(bifurcationLeftParent) > newBifurcationHeight, arg.getNodeHeight(bifurcationRightParent) > newBifurcationHeight };
        if (tester[0] && tester[1]) {
            if (MathUtils.nextBoolean()) {
                splitBifurcationLeftParent = false;
                bifurcationSplitParent = bifurcationRightParent;
            }
        //    		  logHastings += LOG_TWO;
        } else if (tester[0]) {
        //nothing to do
        } else if (tester[1]) {
            splitBifurcationLeftParent = false;
            bifurcationSplitParent = bifurcationRightParent;
        } else {
            assert false;
        }
    }
    boolean attachNewReassortNewBifurcationThroughLeft = MathUtils.nextBoolean();
    //During the delete step, this guy gets cancelled out!
    logHastings += LOG_TWO;
    //        if (sisParentL != sisParentR) {
    //            boolean[] tester = {arg.getNodeHeight(sisParentL) > newBifurcationHeight,
    //                    arg.getNodeHeight(sisParentR) > newBifurcationHeight};
    //
    //            if (tester[0] && tester[1]) {
    //                if (MathUtils.nextBoolean()) {
    //                    sisParent = sisParentR;
    //                }
    //                logHastings += LOG_TWO;
    //            } else if (tester[0]) {
    //                sisParent = sisParentL;
    //            } else {
    //                sisParent = sisParentR;
    //            }
    //        }
    //5. Make the new nodes.
    //Note: The height stuff is taken care of below.
    Node newReassortment = arg.new Node();
    newReassortment.bifurcation = false;
    double[] reassortmentValues = { 1.0 };
    if (relaxed) {
        reassortmentValues = ratePrior.generateValues();
    }
    //        logHastings += ratePrior.getAddHastingsRatio(reassortmentValues);
    newReassortment.rateParameter = new Parameter.Default(reassortmentValues);
    newReassortment.rateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0, reassortmentValues.length));
    newReassortment.number = arg.getNodeCount() + 1;
    Node newBifurcation = arg.new Node();
    double[] bifurcationValues = { 1.0 };
    if (relaxed) {
        bifurcationValues = ratePrior.generateValues();
        logHastings += ratePrior.getAddHastingsRatio(bifurcationValues);
    }
    newBifurcation.rateParameter = new Parameter.Default(bifurcationValues);
    newBifurcation.rateParameter.addBounds(new Parameter.DefaultBounds(Double.POSITIVE_INFINITY, 0, bifurcationValues.length));
    newBifurcation.number = arg.getNodeCount();
    //6. Begin editing the tree.
    arg.beginTreeEdit();
    adjustRandomPartitioning();
    if (newBifurcationHeight < treeHeight) {
        if (bifurcationChild == reassortChild) {
            if (bifurcationChild.bifurcation) {
                bifurcationChild.leftParent = bifurcationChild.rightParent = newReassortment;
                newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                if (bifurcationSplitParent.bifurcation) {
                    if (bifurcationSplitParent.leftChild == bifurcationChild) {
                        bifurcationSplitParent.leftChild = newBifurcation;
                    } else {
                        bifurcationSplitParent.rightChild = newBifurcation;
                    }
                } else {
                    bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                }
                logHastings -= LOG_TWO;
            } else {
                if (splitBifurcationLeftParent && splitReassortLeftParent) {
                    bifurcationChild.leftParent = newReassortment;
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                    newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                } else if (splitBifurcationLeftParent) {
                    //bifurcation on left, reassortment on right
                    bifurcationChild.leftParent = newBifurcation;
                    bifurcationChild.rightParent = newReassortment;
                    newBifurcation.leftChild = bifurcationChild;
                    newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    if (attachNewReassortNewBifurcationThroughLeft) {
                        newReassortment.leftParent = newBifurcation;
                        newReassortment.rightParent = reassortSplitParent;
                    } else {
                        newReassortment.rightParent = newBifurcation;
                        newReassortment.leftParent = reassortSplitParent;
                    }
                    if (reassortSplitParent.bifurcation) {
                        if (reassortSplitParent.leftChild == reassortChild) {
                            reassortSplitParent.leftChild = newReassortment;
                        } else {
                            reassortSplitParent.rightChild = newReassortment;
                        }
                    } else {
                        reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
                    }
                } else if (splitReassortLeftParent) {
                    //bifurcation on right, reassortment on left
                    bifurcationChild.rightParent = newBifurcation;
                    bifurcationChild.leftParent = newReassortment;
                    newBifurcation.leftChild = bifurcationChild;
                    newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    if (attachNewReassortNewBifurcationThroughLeft) {
                        newReassortment.leftParent = newBifurcation;
                        newReassortment.rightParent = reassortSplitParent;
                    } else {
                        newReassortment.rightParent = newBifurcation;
                        newReassortment.leftParent = reassortSplitParent;
                    }
                    if (reassortSplitParent.bifurcation) {
                        if (reassortSplitParent.leftChild == reassortChild) {
                            reassortSplitParent.leftChild = newReassortment;
                        } else {
                            reassortSplitParent.rightChild = newReassortment;
                        }
                    } else {
                        reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
                    }
                } else {
                    bifurcationChild.rightParent = newReassortment;
                    newReassortment.leftChild = newReassortment.rightChild = bifurcationChild;
                    newReassortment.leftParent = newReassortment.rightParent = newBifurcation;
                    newBifurcation.leftChild = newBifurcation.rightChild = newReassortment;
                    newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
                    if (bifurcationSplitParent.bifurcation) {
                        if (bifurcationSplitParent.leftChild == bifurcationChild) {
                            bifurcationSplitParent.leftChild = newBifurcation;
                        } else {
                            bifurcationSplitParent.rightChild = newBifurcation;
                        }
                    } else {
                        bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
                    }
                    logHastings -= LOG_TWO;
                }
            }
        } else {
            newReassortment.leftChild = newReassortment.rightChild = reassortChild;
            newBifurcation.leftParent = newBifurcation.rightParent = bifurcationSplitParent;
            newBifurcation.leftChild = newReassortment;
            newBifurcation.rightChild = bifurcationChild;
            if (attachNewReassortNewBifurcationThroughLeft) {
                newReassortment.leftParent = newBifurcation;
                newReassortment.rightParent = reassortSplitParent;
            } else {
                newReassortment.rightParent = newBifurcation;
                newReassortment.leftParent = reassortSplitParent;
            }
            if (reassortChild.bifurcation) {
                reassortChild.leftParent = reassortChild.rightParent = newReassortment;
            } else if (splitReassortLeftParent) {
                reassortChild.leftParent = newReassortment;
            } else {
                reassortChild.rightParent = newReassortment;
            }
            if (reassortSplitParent.bifurcation) {
                if (reassortSplitParent.leftChild == reassortChild) {
                    reassortSplitParent.leftChild = newReassortment;
                } else {
                    reassortSplitParent.rightChild = newReassortment;
                }
            } else {
                reassortSplitParent.leftChild = reassortSplitParent.rightChild = newReassortment;
            }
            if (bifurcationChild.bifurcation) {
                bifurcationChild.leftParent = bifurcationChild.rightParent = newBifurcation;
            } else if (splitBifurcationLeftParent) {
                bifurcationChild.leftParent = newBifurcation;
            } else {
                bifurcationChild.rightParent = newBifurcation;
            }
            if (bifurcationSplitParent.bifurcation) {
                if (bifurcationSplitParent.leftChild == bifurcationChild) {
                    bifurcationSplitParent.leftChild = newBifurcation;
                } else {
                    bifurcationSplitParent.rightChild = newBifurcation;
                }
            } else {
                bifurcationSplitParent.leftChild = bifurcationSplitParent.rightChild = newBifurcation;
            }
        }
        Parameter partition = new Parameter.Default(arg.getNumberOfPartitions());
        drawRandomPartitioning(partition);
        newReassortment.partitioning = partition;
        newBifurcation.heightParameter = new Parameter.Default(newBifurcationHeight);
        newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
        newBifurcation.setupHeightBounds();
        newReassortment.setupHeightBounds();
        arg.expandARG(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        //                     nodeRates);
        assert nodeCheck() : arg.toARGSummary();
    } else {
        assert newReassortmentHeight < treeHeight;
        //New bifurcation takes the place of the old root.
        //Much easier to program.
        newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
        newReassortment.setupHeightBounds();
        bifurcationChild = newBifurcation;
        if (arg.isRoot(reassortSplitParent))
            reassortSplitParent = newBifurcation;
        Node root = (Node) arg.getRoot();
        Node rootLeftChild = root.leftChild;
        Node rootRightChild = root.rightChild;
        arg.singleRemoveChild(root, rootLeftChild);
        arg.singleRemoveChild(root, rootRightChild);
        arg.singleAddChild(newBifurcation, rootLeftChild);
        arg.singleAddChild(newBifurcation, rootRightChild);
        if (reassortSplitParent.isBifurcation())
            arg.singleRemoveChild(reassortSplitParent, reassortChild);
        else
            arg.doubleRemoveChild(reassortSplitParent, reassortChild);
        arg.doubleAddChild(newReassortment, reassortChild);
        arg.singleAddChild(root, newBifurcation);
        Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
        drawRandomPartitioning(partitioning);
        arg.addChildAsRecombinant(root, reassortSplitParent, newReassortment, partitioning);
        if (attachNewReassortNewBifurcationThroughLeft) {
            newReassortment.leftParent = root;
            newReassortment.rightParent = reassortSplitParent;
        } else {
            newReassortment.leftParent = reassortSplitParent;
            newReassortment.rightParent = root;
        }
        newBifurcation.heightParameter = new Parameter.Default(root.getHeight());
        newBifurcation.setupHeightBounds();
        root.heightParameter.setParameterValue(0, newBifurcationHeight);
        arg.expandARG(newBifurcation, newReassortment, internalNodeParameters, internalAndRootNodeParameters, nodeRates);
        assert nodeCheck();
    }
    //6a. This is when we do not create a new root.
    //        if (newBifurcationHeight < treeHeight) {
    //            newBifurcation.heightParameter = new Parameter.Default(newBifurcationHeight);
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            newBifurcation.setupHeightBounds();
    //            newReassortment.setupHeightBounds();
    //
    //            if (sisParent.bifurcation)
    //                arg.singleRemoveChild(sisParent, bifurcationChild);
    //            else
    //                arg.doubleRemoveChild(sisParent, bifurcationChild);
    //            if (bifurcationChild != reassortChild) {
    //                if (recParent.bifurcation)
    //                    arg.singleRemoveChild(recParent, reassortChild);
    //                else
    //                    arg.doubleRemoveChild(recParent, reassortChild);
    //            }
    //            if (sisParent.bifurcation)
    //                arg.singleAddChild(sisParent, newBifurcation);
    //            else
    //                arg.doubleAddChild(sisParent, newBifurcation);
    //            if (bifurcationChild != reassortChild)
    //                arg.singleAddChild(newBifurcation, bifurcationChild);
    //            arg.doubleAddChild(newReassortment, reassortChild);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //            if (bifurcationChild != reassortChild) {
    //                arg.addChildAsRecombinant(newBifurcation, recParent,
    //                        newReassortment, partitioning);
    //            } else {
    //                arg.addChildAsRecombinant(newBifurcation, newBifurcation,
    //                        newReassortment, partitioning);
    //            }
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters,
    //                    internalAndRootNodeParameters,
    //                    nodeRates);
    //            assert nodeCheck();
    //
    //            //6b. But here we do.
    //        } else if (newReassortmentHeight < treeHeight) {
    //
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            newReassortment.setupHeightBounds();
    //
    //            bifurcationChild = newBifurcation;
    //            if (arg.isRoot(recParent))
    //                recParent = newBifurcation;
    //
    //
    //            Node root = (Node) arg.getRoot();
    //            Node rootLeftChild = root.leftChild;
    //            Node rootRightChild = root.rightChild;
    //
    //            arg.singleRemoveChild(root, rootLeftChild);
    //            arg.singleRemoveChild(root, rootRightChild);
    //            arg.singleAddChild(newBifurcation, rootLeftChild);
    //            arg.singleAddChild(newBifurcation, rootRightChild);
    //
    //            if (recParent.isBifurcation())
    //                arg.singleRemoveChild(recParent, reassortChild);
    //            else
    //                arg.doubleRemoveChild(recParent, reassortChild);
    //
    //            arg.doubleAddChild(newReassortment, reassortChild);
    //            arg.singleAddChild(root, newBifurcation);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //
    //            arg.addChildAsRecombinant(root, recParent, newReassortment, partitioning);
    //
    //            newBifurcation.heightParameter = new Parameter.Default(root.getHeight());
    //
    //            newBifurcation.setupHeightBounds();
    //            root.heightParameter.setParameterValue(0, newBifurcationHeight);
    //
    //
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters, internalAndRootNodeParameters,
    //                    nodeRates);
    //
    //            assert nodeCheck();
    //
    //        } else {
    //
    //            Node root = (Node) arg.getRoot();
    //            Node rootLeftChild = root.leftChild;
    //            Node rootRightChild = root.rightChild;
    //
    //            arg.singleRemoveChild(root, rootLeftChild);
    //            arg.singleRemoveChild(root, rootRightChild);
    //            arg.singleAddChild(newBifurcation, rootLeftChild);
    //            arg.singleAddChild(newBifurcation, rootRightChild);
    //
    //            arg.doubleAddChild(newReassortment, newBifurcation);
    //            arg.doubleAddChild(root, newReassortment);
    //
    //            Parameter partitioning = new Parameter.Default(arg.getNumberOfPartitions());
    //            drawRandomPartitioning(partitioning);
    //
    //            newReassortment.partitioning = partitioning;
    //
    //            newBifurcation.heightParameter = new Parameter.Default(arg.getNodeHeight(root));
    //            newReassortment.heightParameter = new Parameter.Default(newReassortmentHeight);
    //            root.heightParameter.setParameterValueQuietly(0, newBifurcationHeight);
    //
    //            newBifurcation.setupHeightBounds();
    //            newReassortment.setupHeightBounds();
    //
    //            arg.expandARGWithRecombinant(newBifurcation, newReassortment,
    //                    internalNodeParameters, internalAndRootNodeParameters,
    //                    nodeRates);
    //
    //            assert nodeCheck();
    //
    //        }
    arg.pushTreeSizeIncreasedEvent();
    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();
    logHastings -= Math.log((double) findPotentialNodesToRemove(null));
    assert nodeCheck();
    assert !Double.isNaN(logHastings) && !Double.isInfinite(logHastings);
    if (newReassortment.leftParent.bifurcation && newReassortment.rightParent.bifurcation && newReassortment.leftParent != newReassortment.rightParent) {
        logHastings -= LOG_TWO;
    }
    //You're done, return the hastings ratio!
    //		System.out.println(logHastings);
    logHastings += getPartitionAddHastingsRatio(newReassortment.partitioning.getParameterValues());
    return logHastings;
}

Example 3 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ARGTraceAnalysis method analyzeARG.

/**
	 * Actually analyzes a particular tree using the trace given the burnin
	 */
public final Tree analyzeARG(String target) {
    int n = getTreeCount();
    FlexibleTree meanTree = null;
    for (int i = 0; i < n; i++) {
        Tree tree = getARG(i);
        if (TreeUtils.uniqueNewick(tree, tree.getRoot()).equals(target)) {
            meanTree = new FlexibleTree(tree);
            break;
        }
    }
    if (meanTree == null)
        throw new RuntimeException("No target tree in trace");
    int m = meanTree.getInternalNodeCount();
    for (int j = 0; j < m; j++) {
        double[] heights = new double[n];
        NodeRef node1 = meanTree.getInternalNode(j);
        Set<String> leafSet = TreeUtils.getDescendantLeaves(meanTree, node1);
        for (int i = 0; i < n; i++) {
            Tree tree = getARG(i);
            NodeRef node2 = TreeUtils.getCommonAncestorNode(tree, leafSet);
            heights[i] = tree.getNodeHeight(node2);
        }
        meanTree.setNodeHeight(node1, dr.stats.DiscreteStatistics.mean(heights));
        meanTree.setNodeAttribute(node1, "upper", new Double(dr.stats.DiscreteStatistics.quantile(0.975, heights)));
        meanTree.setNodeAttribute(node1, "lower", new Double(dr.stats.DiscreteStatistics.quantile(0.025, heights)));
    }
    return meanTree;
}

Example 4 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ARGCoalescentLikelihood method calculateIntervals.

public void calculateIntervals() {
    intervals.clear();
    intervals.ensureCapacity(arg.getNodeCount());
    NodeRef x;
    for (int i = 0; i < arg.getInternalNodeCount(); i++) {
        x = arg.getInternalNode(i);
        if (arg.isReassortment(x)) {
            intervals.add(new CoalescentInterval(arg.getNodeHeight(x), RECOMBINATION));
        } else {
            intervals.add(new CoalescentInterval(arg.getNodeHeight(x), COALESCENT));
        }
    }
    for (int i = 0; i < arg.getExternalNodeCount(); i++) {
        x = arg.getExternalNode(i);
        if (arg.getNodeHeight(x) > 0.0) {
            intervals.add(new CoalescentInterval(arg.getNodeHeight(x), NEW_SAMPLE));
        }
    }
    dr.util.HeapSort.sort(intervals);
    double a = 0, b = 0;
    for (int i = 0; i < intervals.size(); i++) {
        b = intervals.get(i).length;
        intervals.get(i).length = intervals.get(i).length - a;
        a = b;
    }
    intervalsKnown = true;
}

Example 5 with NodeRef

use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.

the class ClusterLabelsVirusesStatistic method determine_membership_v2.

//traverse down the tree, top down, do calculation
int[] determine_membership_v2(TreeModel treeModel) {
    NodeRef root = treeModel.getRoot();
    int numClusters = 1;
    LinkedList<NodeRef> list = new LinkedList<NodeRef>();
    list.addFirst(root);
    int[] membership = new int[treeModel.getNodeCount()];
    for (int i = 0; i < treeModel.getNodeCount(); i++) {
        membership[i] = -1;
    }
    //root always given the first cluster
    membership[root.getNumber()] = 0;
    while (!list.isEmpty()) {
        //do things with the current object
        NodeRef curElement = list.pop();
        //String content = "node #" + curElement.getNumber() +", taxon=" + treeModel.getNodeTaxon(curElement) + " and parent is = " ;
        String content = "node #" + curElement.getNumber() + ", taxon= ";
        if (treeModel.getNodeTaxon(curElement) == null) {
            content += "internal node\t";
        } else {
            content += treeModel.getNodeTaxon(curElement).getId() + "\t";
        //content += treeModel.getTaxonIndex(treeModel.getNodeTaxon(curElement)) + "\t";
        }
        if (treeModel.getParent(curElement) == null) {
        //content += "no parent";
        } else {
        //content += "parent node#=" + treeModel.getParent(curElement).getNumber();
        }
        //cluster assignment:
        if (!treeModel.isRoot(curElement)) {
            if ((int) indicators.getParameterValue(curElement.getNumber()) == 1) {
                numClusters++;
                membership[curElement.getNumber()] = numClusters - 1;
            } else {
                //inherit from parent's cluster assignment
                membership[curElement.getNumber()] = membership[treeModel.getParent(curElement).getNumber()];
            }
        }
        //is not Root
        content += " cluster = " + membership[curElement.getNumber()];
        for (int childNum = 0; childNum < treeModel.getChildCount(curElement); childNum++) {
            list.addFirst(treeModel.getChild(curElement, childNum));
        }
    }
    return (membership);
}