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

Example 16 with MutableTreeModel

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

the class MultivariateTraitUtils method computeTreeVariance.

public static double[][] computeTreeVariance(FullyConjugateMultivariateTraitLikelihood trait, boolean conditionOnRoot) {
    MutableTreeModel treeModel = trait.getTreeModel();
    final int tipCount = treeModel.getExternalNodeCount();
    double[][] variance = new double[tipCount][tipCount];
    for (int i = 0; i < tipCount; i++) {
        // Fill in diagonal
        double marginalTime = trait.getRescaledLengthToRoot(treeModel.getExternalNode(i));
        variance[i][i] = marginalTime;
        for (int j = i + 1; j < tipCount; j++) {
            NodeRef mrca = findMRCA(trait, i, j);
            if (DEBUG) {
                System.err.println(trait.getTreeModel().getRoot().getNumber());
                System.err.print("Taxa pair: " + i + " : " + j + " (" + mrca.getNumber() + ") = ");
            }
            double length = trait.getRescaledLengthToRoot(mrca);
            if (DEBUG) {
                System.err.println(length);
            }
            variance[i][j] = length;
        }
    }
    // Make symmetric
    for (int i = 0; i < tipCount; i++) {
        for (int j = i + 1; j < tipCount; j++) {
            variance[j][i] = variance[i][j];
        }
    }
    if (DEBUG) {
        System.err.println("");
        System.err.println("New tree conditional variance:\n" + new Matrix(variance));
    }
    if (!conditionOnRoot) {
        double priorSampleSize = trait.getPriorSampleSize();
        for (int i = 0; i < tipCount; ++i) {
            for (int j = 0; j < tipCount; ++j) {
                variance[i][j] += 1.0 / priorSampleSize;
            }
        }
        if (DEBUG) {
            System.err.println("");
            System.err.println("New tree unconditional variance:\n" + new Matrix(variance));
        }
    }
    return variance;
}

Example 17 with MutableTreeModel

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

the class DiffusionRateCovarianceStatistic method getStatisticValue.

/**
 * @return whatever Philippe wants
 */
public double getStatisticValue(int dim) {
    String traitName = traitLikelihoods.get(0).getTraitName();
    for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
        MutableTreeModel tree = traitLikelihood.getTreeModel();
        int counter = 0;
        int index = 0;
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef child = tree.getNode(i);
            NodeRef parent = tree.getParent(child);
            if (parent != null & !tree.isRoot(parent)) {
                double[] childTrait = traitLikelihood.getTraitForNode(tree, child, traitName);
                double[] parentTrait = traitLikelihood.getTraitForNode(tree, parent, traitName);
                double childTime = tree.getBranchLength(child);
                double parentTime = tree.getBranchLength(parent);
                NodeRef grandParent = tree.getParent(parent);
                double[] grandParentTrait = traitLikelihood.getTraitForNode(tree, grandParent, traitName);
                if (useGreatCircleDistances && (childTrait.length == 2)) {
                    // Great Circle distance
                    SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(childTrait[0], childTrait[1]);
                    SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(parentTrait[0], parentTrait[1]);
                    double childDistance = coord1.distance(coord2);
                    SphericalPolarCoordinates coord3 = new SphericalPolarCoordinates(grandParentTrait[0], grandParentTrait[1]);
                    double parentDistance = coord2.distance(coord3);
                    if (!diffusionCoefficient) {
                        childRate[index] = childDistance / childTime;
                        parentRate[index] = parentDistance / parentTime;
                    } else {
                        childRate[index] = Math.pow(childDistance, 2) / (4 * childTime);
                        parentRate[index] = Math.pow(parentDistance, 2) / (4 * parentTime);
                    }
                } else {
                    double childDistance = getNativeDistance(childTrait, parentTrait);
                    double parentDistance = getNativeDistance(parentTrait, grandParentTrait);
                    if (!diffusionCoefficient) {
                        childRate[index] = childDistance / childTime;
                        parentRate[index] = parentDistance / parentTime;
                    } else {
                        childRate[index] = Math.pow(childDistance, 2) / (4 * childTime);
                        parentRate[index] = Math.pow(parentDistance, 2) / (4 * parentTime);
                    }
                }
                index += 1;
            }
        }
    }
    return DiscreteStatistics.covariance(childRate, parentRate);
}