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

Example 81 with NodeRef

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

the class AbstractTreeLikelihood method updateNodeAndChildren.

/**
     * Set update flag for a node and its direct children
     */
protected void updateNodeAndChildren(NodeRef node) {
    if (COUNT_TOTAL_OPERATIONS)
        totalRateUpdateSingleCount++;
    updateNode[node.getNumber()] = true;
    for (int i = 0; i < treeModel.getChildCount(node); i++) {
        if (COUNT_TOTAL_OPERATIONS)
            totalRateUpdateSingleCount++;
        NodeRef child = treeModel.getChild(node, i);
        updateNode[child.getNumber()] = true;
    }
    likelihoodKnown = false;
}

Example 82 with NodeRef

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

the class AncestralStateBeagleTreeLikelihood method traverseSample.

public void traverseSample(TreeModel tree, NodeRef node, int[] parentState, int[] rateCategory) {
    int nodeNum = node.getNumber();
    NodeRef parent = tree.getParent(node);
    // This function assumes that all partial likelihoods have already been calculated
    // If the node is internal, then sample its state given the state of its parent (pre-order traversal).
    double[] conditionalProbabilities = new double[stateCount];
    int[] state = new int[patternCount];
    if (!tree.isExternal(node)) {
        if (parent == null) {
            // This is the root node
            getPartials(nodeNum, partials);
            boolean sampleCategory = categoryCount > 1;
            double[] posteriorWeightedCategory = null;
            double[] priorWeightedCategory = null;
            if (sampleCategory) {
                rateCategory = new int[patternCount];
                posteriorWeightedCategory = new double[categoryCount];
                priorWeightedCategory = siteRateModel.getCategoryProportions();
            }
            for (int j = 0; j < patternCount; j++) {
                // Sample across-site-rate-variation, if it exists
                if (sampleCategory) {
                    for (int r = 0; r < categoryCount; r++) {
                        posteriorWeightedCategory[r] = 0;
                        for (int k = 0; k < stateCount; k++) {
                            posteriorWeightedCategory[r] += partials[r * stateCount * patternCount + j * stateCount + k];
                        }
                        posteriorWeightedCategory[r] *= priorWeightedCategory[r];
                    }
                    rateCategory[j] = drawChoice(posteriorWeightedCategory);
                }
                // Sample root character state
                int partialsIndex = (rateCategory == null ? 0 : rateCategory[j]) * stateCount * patternCount;
                System.arraycopy(partials, partialsIndex + j * stateCount, conditionalProbabilities, 0, stateCount);
                // TODO May have more than one set of frequencies
                double[] frequencies = substitutionModelDelegate.getRootStateFrequencies();
                for (int i = 0; i < stateCount; i++) {
                    conditionalProbabilities[i] *= frequencies[i];
                }
                try {
                    state[j] = drawChoice(conditionalProbabilities);
                } catch (Error e) {
                    System.err.println(e.toString());
                    System.err.println("Please report error to Marc");
                    state[j] = 0;
                }
                reconstructedStates[nodeNum][j] = state[j];
                if (!returnMarginalLogLikelihood) {
                    jointLogLikelihood += Math.log(frequencies[state[j]]);
                }
            }
            if (sampleCategory) {
                if (this.rateCategory == null) {
                    this.rateCategory = new int[patternCount];
                }
                System.arraycopy(rateCategory, 0, this.rateCategory, 0, patternCount);
            }
        } else {
            // This is an internal node, but not the root
            double[] partialLikelihood = new double[stateCount * patternCount * categoryCount];
            getPartials(nodeNum, partialLikelihood);
            // Sibon says that this actually works now
            //                if (categoryCount > 1)
            //                    throw new RuntimeException("Reconstruction not implemented for multiple categories yet.");
            getMatrix(nodeNum, probabilities);
            for (int j = 0; j < patternCount; j++) {
                int parentIndex = parentState[j] * stateCount;
                int childIndex = j * stateCount;
                int category = rateCategory == null ? 0 : rateCategory[j];
                int matrixIndex = category * stateCount * stateCount;
                int partialIndex = category * stateCount * patternCount;
                for (int i = 0; i < stateCount; i++) conditionalProbabilities[i] = partialLikelihood[partialIndex + childIndex + i] * probabilities[matrixIndex + parentIndex + i];
                state[j] = drawChoice(conditionalProbabilities);
                reconstructedStates[nodeNum][j] = state[j];
                if (!returnMarginalLogLikelihood) {
                    double contrib = probabilities[parentIndex + state[j]];
                    jointLogLikelihood += Math.log(contrib);
                }
            }
            hookCalculation(tree, parent, node, parentState, state, probabilities, rateCategory);
        }
        // Traverse down the two child nodes
        NodeRef child1 = tree.getChild(node, 0);
        traverseSample(tree, child1, state, rateCategory);
        NodeRef child2 = tree.getChild(node, 1);
        traverseSample(tree, child2, state, rateCategory);
    } else {
        if (useAmbiguities()) {
            getMatrix(nodeNum, probabilities);
            double[] partials = tipPartials[nodeNum];
            for (int j = 0; j < patternCount; j++) {
                final int parentIndex = parentState[j] * stateCount;
                int category = rateCategory == null ? 0 : rateCategory[j];
                int matrixIndex = category * stateCount * stateCount;
                System.arraycopy(probabilities, parentIndex + matrixIndex, conditionalProbabilities, 0, stateCount);
                for (int k = 0; k < stateCount; ++k) {
                    conditionalProbabilities[k] *= partials[j * stateCount + k];
                }
                reconstructedStates[nodeNum][j] = drawChoice(conditionalProbabilities);
                if (!returnMarginalLogLikelihood) {
                    double contrib = probabilities[parentIndex + reconstructedStates[nodeNum][j]];
                    jointLogLikelihood += Math.log(contrib);
                }
            }
        } else {
            getTipStates(nodeNum, reconstructedStates[nodeNum]);
            for (int j = 0; j < patternCount; j++) {
                final int thisState = reconstructedStates[nodeNum][j];
                if (dataType.isAmbiguousState(thisState)) {
                    final int parentIndex = parentState[j] * stateCount;
                    int category = rateCategory == null ? 0 : rateCategory[j];
                    int matrixIndex = category * stateCount * stateCount;
                    getMatrix(nodeNum, probabilities);
                    System.arraycopy(probabilities, parentIndex + matrixIndex, conditionalProbabilities, 0, stateCount);
                    if (useAmbiguities && !dataType.isUnknownState(thisState)) {
                        // Not completely unknown
                        boolean[] stateSet = dataType.getStateSet(thisState);
                        for (int k = 0; k < stateCount; k++) {
                            if (!stateSet[k]) {
                                conditionalProbabilities[k] = 0.0;
                            }
                        }
                    }
                    reconstructedStates[nodeNum][j] = drawChoice(conditionalProbabilities);
                }
                if (!returnMarginalLogLikelihood) {
                    final int parentIndex = parentState[j] * stateCount;
                    getMatrix(nodeNum, probabilities);
                    if (!returnMarginalLogLikelihood) {
                        double contrib = probabilities[parentIndex + reconstructedStates[nodeNum][j]];
                        jointLogLikelihood += Math.log(contrib);
                    }
                }
            }
        }
        hookCalculation(tree, parent, node, parentState, reconstructedStates[nodeNum], null, rateCategory);
    }
}

Example 83 with NodeRef

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

the class BeagleOperationReport method calculateLogLikelihood.

protected double calculateLogLikelihood() {
    if (matrixUpdateIndices == null) {
        matrixUpdateIndices = new int[eigenCount][nodeCount];
        branchLengths = new double[eigenCount][nodeCount];
        branchUpdateCount = new int[eigenCount];
    //            scaleBufferIndices = new int[internalNodeCount];
    //            storedScaleBufferIndices = new int[internalNodeCount];
    }
    if (operations == null) {
        operations = new int[numRestrictedPartials + 1][internalNodeCount * Beagle.OPERATION_TUPLE_SIZE];
        operationCount = new int[numRestrictedPartials + 1];
    }
    recomputeScaleFactors = false;
    for (int i = 0; i < eigenCount; i++) {
        branchUpdateCount[i] = 0;
    }
    operationListCount = 0;
    if (hasRestrictedPartials) {
        for (int i = 0; i <= numRestrictedPartials; i++) {
            operationCount[i] = 0;
        }
    } else {
        operationCount[0] = 0;
    }
    System.out.println(alignmentString.toString());
    final NodeRef root = treeModel.getRoot();
    // Do not flip buffers
    traverse(treeModel, root, null, false);
    for (int i = 0; i < eigenCount; i++) {
        if (branchUpdateCount[i] > 0) {
            if (DEBUG_BEAGLE_OPERATIONS) {
                StringBuilder sb = new StringBuilder();
                sb.append("eval = ").append(new Vector(substitutionModel.getEigenDecomposition().getEigenValues())).append("\n");
                sb.append("evec = ").append(new Vector(substitutionModel.getEigenDecomposition().getEigenVectors())).append("\n");
                sb.append("ivec = ").append(new Vector(substitutionModel.getEigenDecomposition().getInverseEigenVectors())).append("\n");
                sb.append("Branch count: ").append(branchUpdateCount[i]);
                sb.append("\nNode indices:\n");
                if (SINGLE_LINE) {
                    sb.append("int n[] = {");
                }
                for (int k = 0; k < branchUpdateCount[i]; ++k) {
                    if (SINGLE_LINE) {
                        sb.append(" ").append(matrixUpdateIndices[i][k]);
                        if (k < (branchUpdateCount[i] - 1)) {
                            sb.append(",");
                        }
                    } else {
                        sb.append(matrixUpdateIndices[i][k]).append("\n");
                    }
                }
                if (SINGLE_LINE) {
                    sb.append(" };\n");
                }
                sb.append("\nBranch lengths:\n");
                if (SINGLE_LINE) {
                    sb.append("double b[] = {");
                }
                for (int k = 0; k < branchUpdateCount[i]; ++k) {
                    if (SINGLE_LINE) {
                        sb.append(" ").append(branchLengths[i][k]);
                        if (k < (branchUpdateCount[i] - 1)) {
                            sb.append(",");
                        }
                    } else {
                        sb.append(branchLengths[i][k]).append("\n");
                    }
                }
                if (SINGLE_LINE) {
                    sb.append(" };\n");
                }
                System.out.println(sb.toString());
            }
        }
    }
    if (DEBUG_BEAGLE_OPERATIONS) {
        StringBuilder sb = new StringBuilder();
        sb.append("Operation count: ").append(operationCount[0]);
        sb.append("\nOperations:\n");
        if (SINGLE_LINE) {
            sb.append("BeagleOperation o[] = {");
        }
        for (int k = 0; k < operationCount[0] * Beagle.OPERATION_TUPLE_SIZE; ++k) {
            if (SINGLE_LINE) {
                sb.append(" ").append(operations[0][k]);
                if (k < (operationCount[0] * Beagle.OPERATION_TUPLE_SIZE - 1)) {
                    sb.append(",");
                }
            } else {
                sb.append(operations[0][k]).append("\n");
            }
        }
        if (SINGLE_LINE) {
            sb.append(" };\n");
        }
        sb.append("Use scale factors: ").append(useScaleFactors).append("\n");
        System.out.println(sb.toString());
    }
    int rootIndex = partialBufferHelper.getOffsetIndex(root.getNumber());
    System.out.println("Root node: " + rootIndex);
    return 0.0;
}

Example 84 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 85 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;
}