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

Example 21 with NodeRef

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

the class BeagleOperationReport method traverse.

/**
     * Traverse the tree calculating partial likelihoods.
     *
     * @param tree           tree
     * @param node           node
     * @param operatorNumber operatorNumber
     * @param flip           flip
     * @return boolean
     */
private boolean traverse(Tree tree, NodeRef node, int[] operatorNumber, boolean flip) {
    boolean update = false;
    int nodeNum = node.getNumber();
    NodeRef parent = tree.getParent(node);
    if (operatorNumber != null) {
        operatorNumber[0] = -1;
    }
    // First update the transition probability matrix(ices) for this branch
    if (parent != null && updateNode[nodeNum]) {
        final double branchRate = branchRateModel.getBranchRate(tree, node);
        // Get the operational time of the branch
        final double branchTime = branchRate * (tree.getNodeHeight(parent) - tree.getNodeHeight(node));
        if (branchTime < 0.0) {
            throw new RuntimeException("Negative branch length: " + branchTime);
        }
        if (flip) {
            // first flip the matrixBufferHelper
            matrixBufferHelper.flipOffset(nodeNum);
        }
        // then set which matrix to update
        //branchSubstitutionModel.getBranchIndex(tree, node);
        final int eigenIndex = 0;
        final int updateCount = branchUpdateCount[eigenIndex];
        matrixUpdateIndices[eigenIndex][updateCount] = matrixBufferHelper.getOffsetIndex(nodeNum);
        branchLengths[eigenIndex][updateCount] = branchTime;
        branchUpdateCount[eigenIndex]++;
        update = true;
    }
    // If the node is internal, update the partial likelihoods.
    if (!tree.isExternal(node)) {
        // Traverse down the two child nodes
        NodeRef child1 = tree.getChild(node, 0);
        final int[] op1 = { -1 };
        final boolean update1 = traverse(tree, child1, op1, flip);
        NodeRef child2 = tree.getChild(node, 1);
        final int[] op2 = { -1 };
        final boolean update2 = traverse(tree, child2, op2, flip);
        // If either child node was updated then update this node too
        if (update1 || update2) {
            int x = operationCount[operationListCount] * Beagle.OPERATION_TUPLE_SIZE;
            if (flip) {
                // first flip the partialBufferHelper
                partialBufferHelper.flipOffset(nodeNum);
            }
            final int[] operations = this.operations[operationListCount];
            operations[x] = partialBufferHelper.getOffsetIndex(nodeNum);
            if (useScaleFactors) {
                // get the index of this scaling buffer
                int n = nodeNum - tipCount;
                if (recomputeScaleFactors) {
                    // flip the indicator: can take either n or (internalNodeCount + 1) - n
                    scaleBufferHelper.flipOffset(n);
                    // store the index
                    scaleBufferIndices[n] = scaleBufferHelper.getOffsetIndex(n);
                    // Write new scaleFactor
                    operations[x + 1] = scaleBufferIndices[n];
                    operations[x + 2] = Beagle.NONE;
                } else {
                    operations[x + 1] = Beagle.NONE;
                    // Read existing scaleFactor
                    operations[x + 2] = scaleBufferIndices[n];
                }
            } else {
                if (useAutoScaling) {
                    scaleBufferIndices[nodeNum - tipCount] = partialBufferHelper.getOffsetIndex(nodeNum);
                }
                // Not using scaleFactors
                operations[x + 1] = Beagle.NONE;
                operations[x + 2] = Beagle.NONE;
            }
            // source node 1
            operations[x + 3] = partialBufferHelper.getOffsetIndex(child1.getNumber());
            // source matrix 1
            operations[x + 4] = matrixBufferHelper.getOffsetIndex(child1.getNumber());
            // source node 2
            operations[x + 5] = partialBufferHelper.getOffsetIndex(child2.getNumber());
            // source matrix 2
            operations[x + 6] = matrixBufferHelper.getOffsetIndex(child2.getNumber());
            operationCount[operationListCount]++;
            update = true;
        }
    }
    return update;
}

Example 22 with NodeRef

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

the class MultivariateTraitDebugUtilities method getTreeVariance.

public static double[][] getTreeVariance(final Tree tree, final double normalization, final double priorSampleSize) {
    final int tipCount = tree.getExternalNodeCount();
    int length = tipCount;
    double[][] variance = new double[length][length];
    for (int i = 0; i < tipCount; i++) {
        // Fill in diagonal
        double marginalTime = getLengthToRoot(tree, tree.getExternalNode(i)) * normalization;
        variance[i][i] = marginalTime;
        for (int j = i + 1; j < tipCount; j++) {
            NodeRef mrca = findMRCA(tree, i, j);
            variance[i][j] = getLengthToRoot(tree, mrca);
        }
    }
    // Make symmetric
    for (int i = 0; i < length; i++) {
        for (int j = i + 1; j < length; j++) {
            variance[j][i] = variance[i][j];
        }
    }
    if (!Double.isInfinite(priorSampleSize)) {
        for (int i = 0; i < variance.length; ++i) {
            for (int j = 0; j < variance[i].length; ++j) {
                variance[i][j] += 1.0 / priorSampleSize;
            }
        }
    }
    return variance;
}

Example 23 with NodeRef

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

the class Partition method simulatePartition.

// END: loadBeagleInstance
public void simulatePartition() {
    try {
        NodeRef root = treeModel.getRoot();
        // gamma category rates
        double[] categoryRates = siteRateModel.getCategoryRates();
        beagle.setCategoryRates(categoryRates);
        // probabilities for gamma category rates
        double[] categoryProbs = siteRateModel.getCategoryProportions();
        //			beagle.setCategoryWeights(0, categoryProbs);
        //			Utils.printArray(categoryRates);
        //			Utils.printArray(categoryProbs);
        int[] category = new int[partitionSiteCount];
        for (int i = 0; i < partitionSiteCount; i++) {
            category[i] = randomChoicePDF(categoryProbs, partitionNumber, "categories");
        }
        if (DEBUG) {
            System.out.println("category for each site:");
            Utils.printArray(category);
        }
        //END: DEBUG
        int[] parentSequence = new int[partitionSiteCount];
        // set ancestral sequence for partition if it exists
        if (hasRootSequence) {
            if (rootSequence.getLength() == partitionSiteCount) {
                parentSequence = sequence2intArray(rootSequence);
            } else if (dataType instanceof Codons && rootSequence.getLength() == 3 * partitionSiteCount) {
                parentSequence = sequence2intArray(rootSequence);
            } else {
                throw new RuntimeException("Ancestral sequence length of " + rootSequence.getLength() + " does not match partition site count of " + partitionSiteCount + ".");
            }
        } else {
            double[] frequencies = freqModel.getFrequencies();
            for (int i = 0; i < partitionSiteCount; i++) {
                parentSequence[i] = randomChoicePDF(frequencies, partitionNumber, "root");
            }
        }
        if (DEBUG) {
            synchronized (this) {
                System.out.println();
                System.out.println("root Sequence:");
                Utils.printArray(parentSequence);
            }
        }
        //END: DEBUG
        substitutionModelDelegate.updateSubstitutionModels(beagle);
        traverse(root, parentSequence, category);
        if (DEBUG) {
            synchronized (this) {
                System.out.println("Simulated alignment:");
                printSequences();
            }
        }
        //END: DEBUG
        beagle.finalize();
    } catch (Exception e) {
        e.printStackTrace();
    } catch (Throwable e) {
        System.err.println("BeagleException: " + e.getMessage());
        System.exit(-1);
    }
}

Example 24 with NodeRef

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

the class TreeUtils method setHeightsFromDates.

/**
     * Sets the tip heights from the tip dates
     */
public static void setHeightsFromDates(MutableTree tree) {
    Date mostRecent = null;
    for (int i = 0; i < tree.getExternalNodeCount(); i++) {
        Taxon taxon = tree.getNodeTaxon(tree.getExternalNode(i));
        Date date = (Date) taxon.getAttribute("date");
        if (date != null) {
            if ((mostRecent == null) || date.after(mostRecent)) {
                mostRecent = date;
            }
        }
    }
    TimeScale timeScale = new TimeScale(mostRecent.getUnits(), true, mostRecent.getAbsoluteTimeValue());
    for (int i = 0; i < tree.getExternalNodeCount(); i++) {
        NodeRef node = tree.getExternalNode(i);
        Taxon taxon = tree.getNodeTaxon(node);
        Date date = (Date) taxon.getAttribute("date");
        if (date != null) {
            double height = timeScale.convertTime(date.getTimeValue(), date);
            tree.setNodeHeight(node, height);
        } else {
            tree.setNodeHeight(node, 0.0);
        }
    }
    adjustInternalHeights(tree, tree.getRoot());
    if (mostRecent != null) {
        tree.setUnits(mostRecent.getUnits());
    }
}

Example 25 with NodeRef

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

the class GeoDiffusionSimulator method traverse.

void traverse(NodeRef node, double[] parentSequence, double[][] latLongs) {
    for (int iChild = 0; iChild < m_tree.getChildCount(node); iChild++) {
        NodeRef child = m_tree.getChild(node, iChild);
        //find the branch length
        final double branchRate = m_branchRateModel.getBranchRate(m_tree, child);
        final double branchLength = branchRate * (m_tree.getNodeHeight(node) - m_tree.getNodeHeight(child));
        if (branchLength < 0.0) {
            throw new RuntimeException("Negative branch length: " + branchLength);
        }
        double childLat = MathUtils.nextGaussian() * Math.sqrt(branchLength) + parentSequence[LATITUDE_INDEX];
        double childLong = MathUtils.nextGaussian() * Math.sqrt(branchLength) + parentSequence[LONGITUDE_INDEX];
        int childNum = child.getNumber();
        latLongs[childNum][LATITUDE_INDEX] = childLat;
        latLongs[childNum][LONGITUDE_INDEX] = childLong;
        traverse(m_tree.getChild(node, iChild), latLongs[childNum], latLongs);
    }
}