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

Example 61 with NodeRef

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

the class DiscretizedLocationOperator method randomizeNodes.

public void randomizeNodes() {
    List<Point2D> listLocations = new ArrayList<Point2D>();
    listLocations.addAll(allLocations);
    for (int i = 0; i < treeModel.getInternalNodeCount(); i++) {
        NodeRef node = treeModel.getInternalNode(i);
        double[] trait = treeModel.getMultivariateNodeTrait(node, traitName);
        Point2D newPt = listLocations.get(MathUtils.nextInt(listLocations.size()));
        trait[0] = newPt.getX();
        trait[1] = newPt.getY();
        recursivelySetTrait(node, trait, null);
    //            treeModel.setMultivariateTrait(node, traitName, trait);            
    }
    System.err.println("Done with randomization");
//        System.exit(-1);
}

Example 62 with NodeRef

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

the class LatentLiabilityGibbs method doPostOrderTraversal.

//Fill out partial mean and precision values in post order
public void doPostOrderTraversal(NodeRef node) {
    // TODO This is already computed IntegratedMultivariateTraitLikelihood
    final int thisNumber = node.getNumber();
    if (treeModel.isExternal(node)) {
        // writes trait values and precision values for tips
        double[] traitValue = getNodeTrait(node);
        for (int j = 0; j < dim; j++) {
            postMeans[thisNumber][j] = traitValue[j];
        }
        postP[thisNumber] = 1 / traitModel.getRescaledBranchLengthForPrecision(node);
        return;
    }
    final NodeRef childNode0 = treeModel.getChild(node, 0);
    final NodeRef childNode1 = treeModel.getChild(node, 1);
    doPostOrderTraversal(childNode0);
    doPostOrderTraversal(childNode1);
    if (!treeModel.isRoot(node)) {
        final int childNumber0 = childNode0.getNumber();
        final int childNumber1 = childNode1.getNumber();
        // precision values
        final double precision0 = postP[childNumber0];
        final double precision1 = postP[childNumber1];
        final double thisPrecision = 1 / traitModel.getRescaledBranchLengthForPrecision(node);
        double tp = precision0 + precision1;
        postP[thisNumber] = tp * thisPrecision / (tp + thisPrecision);
        //mean values
        for (int j = 0; j < dim; j++) {
            postMeans[thisNumber][j] = (precision0 * postMeans[childNumber0][j] + precision1 * postMeans[childNumber1][j]) / (precision0 + precision1);
        }
    }
}

Example 63 with NodeRef

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

the class SubtreeJumpOperator method getReverseProbability.

private double getReverseProbability(Tree tree, NodeRef originalNode, NodeRef targetNode, double height, double maxAge, List<NodeRef> intersectingEdges, double alpha) {
    double[] weights = new double[intersectingEdges.size()];
    double sum = 0.0;
    int i = 0;
    int originalIndex = -1;
    for (NodeRef node1 : intersectingEdges) {
        assert (node1 != targetNode);
        double age = tree.getNodeHeight(TreeUtils.getCommonAncestor(tree, targetNode, node1)) - height;
        age = age / maxAge;
        weights[i] = getJumpWeight(age, alpha);
        sum += weights[i];
        if (node1 == originalNode) {
            originalIndex = i;
        }
        i++;
    }
    return weights[originalIndex] /= sum;
}

Example 64 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 65 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;
}