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

Example 1 with MultivariateTraitTree

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

the class DiffusionRateStatistic method getStatisticValue.

public double getStatisticValue(int dim) {
    String traitName = traitLikelihoods.get(0).getTraitName();
    double treelength = 0;
    double treeDistance = 0;
    double maxDistanceFromRoot = 0;
    double maxDistanceOverTimeFromRoot = 0;
    //double[] rates =  null;
    List<Double> rates = new ArrayList<Double>();
    //double[] diffusionCoefficients =  null;
    List<Double> diffusionCoefficients = new ArrayList<Double>();
    double waDiffusionCoefficient = 0;
    double lowerHeight = heightLowers[dim];
    double upperHeight = Double.MAX_VALUE;
    if (heightLowers.length == 1) {
        upperHeight = heightUpper;
    } else {
        if (dim > 0) {
            if (!cumulative) {
                upperHeight = heightLowers[dim - 1];
            }
        }
    }
    for (AbstractMultivariateTraitLikelihood traitLikelihood : traitLikelihoods) {
        MultivariateTraitTree tree = traitLikelihood.getTreeModel();
        BranchRateModel branchRates = traitLikelihood.getBranchRateModel();
        for (int i = 0; i < tree.getNodeCount(); i++) {
            NodeRef node = tree.getNode(i);
            if (node != tree.getRoot()) {
                NodeRef parentNode = tree.getParent(node);
                if ((tree.getNodeHeight(parentNode) > lowerHeight) && (tree.getNodeHeight(node) < upperHeight)) {
                    double[] trait = traitLikelihood.getTraitForNode(tree, node, traitName);
                    double[] parentTrait = traitLikelihood.getTraitForNode(tree, parentNode, traitName);
                    double[] traitUp = parentTrait;
                    double[] traitLow = trait;
                    double timeUp = tree.getNodeHeight(parentNode);
                    double timeLow = tree.getNodeHeight(node);
                    double rate = (branchRates != null ? branchRates.getBranchRate(tree, node) : 1.0);
                    MultivariateDiffusionModel diffModel = traitLikelihood.diffusionModel;
                    double[] precision = diffModel.getPrecisionParameter().getParameterValues();
                    if (tree.getNodeHeight(parentNode) > upperHeight) {
                        timeUp = upperHeight;
                        //TODO: implement TrueNoise??
                        traitUp = imputeValue(trait, parentTrait, upperHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    if (tree.getNodeHeight(node) < lowerHeight) {
                        timeLow = lowerHeight;
                        traitLow = imputeValue(trait, parentTrait, lowerHeight, tree.getNodeHeight(node), tree.getNodeHeight(parentNode), precision, rate, false);
                    }
                    double time = timeUp - timeLow;
                    treelength += time;
                    double[] rootTrait = traitLikelihood.getTraitForNode(tree, tree.getRoot(), traitName);
                    if (useGreatCircleDistances && (trait.length == 2)) {
                        // Great Circle distance
                        SphericalPolarCoordinates coord1 = new SphericalPolarCoordinates(traitLow[0], traitLow[1]);
                        SphericalPolarCoordinates coord2 = new SphericalPolarCoordinates(traitUp[0], traitUp[1]);
                        double distance = coord1.distance(coord2);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        SphericalPolarCoordinates rootCoord = new SphericalPolarCoordinates(rootTrait[0], rootTrait[1]);
                        double tempDistanceFromRoot = rootCoord.distance(coord2);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    } else {
                        double distance = getNativeDistance(traitLow, traitUp);
                        treeDistance += distance;
                        double dc = Math.pow(distance, 2) / (4 * time);
                        diffusionCoefficients.add(dc);
                        waDiffusionCoefficient += dc * time;
                        rates.add(distance / time);
                        double tempDistanceFromRoot = getNativeDistance(traitLow, rootTrait);
                        if (tempDistanceFromRoot > maxDistanceFromRoot) {
                            maxDistanceFromRoot = tempDistanceFromRoot;
                            maxDistanceOverTimeFromRoot = tempDistanceFromRoot / (tree.getNodeHeight(tree.getRoot()) - timeLow);
                            //distance between traitLow and traitUp for maxDistanceFromRoot
                            if (timeUp == upperHeight) {
                                maxDistanceFromRoot = distance;
                                maxDistanceOverTimeFromRoot = distance / time;
                            }
                        }
                    }
                }
            }
        }
    }
    if (summaryStat == summaryStatistic.DIFFUSION_RATE) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(rates));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(rates));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(rates));
            return Math.sqrt(DiscreteStatistics.variance(toArray(rates), mean)) / mean;
        } else {
            return treeDistance / treelength;
        }
    } else if (summaryStat == summaryStatistic.DIFFUSION_COEFFICIENT) {
        if (summaryMode == Mode.AVERAGE) {
            return DiscreteStatistics.mean(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.MEDIAN) {
            return DiscreteStatistics.median(toArray(diffusionCoefficients));
        } else if (summaryMode == Mode.COEFFICIENT_OF_VARIATION) {
            // don't compute mean twice
            final double mean = DiscreteStatistics.mean(toArray(diffusionCoefficients));
            return Math.sqrt(DiscreteStatistics.variance(toArray(diffusionCoefficients), mean)) / mean;
        } else {
            return waDiffusionCoefficient / treelength;
        }
    } else if (summaryStat == summaryStatistic.WAVEFRONT_DISTANCE) {
        return maxDistanceFromRoot;
    } else {
        return maxDistanceOverTimeFromRoot;
    }
}

Example 2 with MultivariateTraitTree

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

the class MultivariateTraitUtils method computeLinCombMatrix.

private static double[][] computeLinCombMatrix(FullyConjugateMultivariateTraitLikelihood trait) {
    MultivariateTraitTree treeModel = trait.getTreeModel();
    final int tipCount = treeModel.getExternalNodeCount();
    final int branchCount = 2 * tipCount - 2;
    double[][] linCombMatrix = new double[tipCount][branchCount];
    double tempScalar;
    NodeRef tempNode;
    int tempNodeIndex;
    for (int k = 0; k < tipCount; k++) {
        tempNode = treeModel.getExternalNode(k);
        //check if treeModel.getExternalNode(k).getNumber() == k
        tempScalar = 1;
        tempNodeIndex = k;
        for (int r = 0; r < branchCount; r++) {
            if (r == tempNodeIndex) {
                linCombMatrix[k][r] = tempScalar;
                // tempScalar = tempScalar * (1 - treeModel.getBranchLength(tempNode) * trait.getTimeScaledSelection(tempNode));
                tempScalar = tempScalar * Math.exp(-trait.getTimeScaledSelection(tempNode));
                tempNode = treeModel.getParent(tempNode);
                tempNodeIndex = tempNode.getNumber();
            } else {
                linCombMatrix[k][r] = 0;
            }
        }
    }
    return linCombMatrix;
}

Example 3 with MultivariateTraitTree

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

the class MultivariateTraitUtils method computeTreeVarianceOU.

public static double[][] computeTreeVarianceOU(FullyConjugateMultivariateTraitLikelihood trait, boolean conditionOnRoot) {
    MultivariateTraitTree treeModel = trait.getTreeModel();
    final int tipCount = treeModel.getExternalNodeCount();
    final int branchCount = 2 * tipCount - 2;
    double[][] variance = new double[tipCount][tipCount];
    double[][] tempMatrix = new double[tipCount][branchCount];
    double[][] diagMatrix = new double[branchCount][branchCount];
    /*
        double tempScalar;
        NodeRef tempNode;
        int tempNodeIndex;

        for(int k = 0; k < tipCount; k++){
           tempNode = treeModel.getExternalNode(k);
            //check if treeModel.getExternalNode(k).getNumber() == k
           tempScalar = 1;
           tempNodeIndex = k;

           for(int r = 0; r < branchCount; r++){
               if(r == tempNodeIndex){
                   tempMatrix[k][r] = tempScalar;
                   tempScalar = tempScalar*(1-treeModel.getBranchLength(tempNode)*trait.getTimeScaledSelection(tempNode));
                   tempNode = treeModel.getParent(tempNode);
                   tempNodeIndex = tempNode.getNumber();
               }else{
                   tempMatrix[k][r]= 0;
               }

           }
        }
        */
    tempMatrix = computeLinCombMatrix(trait);
    for (int i = 0; i < branchCount; i++) {
        //  diagMatrix[i][i] = treeModel.getBranchLength(treeModel.getNode(i));
        // diagMatrix[i][i] = (2*strengthOfSelection.getBranchRate(treeModel, node) ) / (1 - Math.exp(-2*getTimeScaledSelection(node)));
        diagMatrix[i][i] = (1 - Math.exp(-2 * trait.getTimeScaledSelection(treeModel.getNode(i)))) / (2 * trait.strengthOfSelection.getBranchRate(treeModel, treeModel.getNode(i)));
    }
    variance = productMatrices(productMatrices(tempMatrix, diagMatrix), transposeMatrix(tempMatrix));
    return variance;
}

Example 4 with MultivariateTraitTree

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

the class MultivariateTraitUtils method computeRootMultipliers.

private static double[] computeRootMultipliers(FullyConjugateMultivariateTraitLikelihood trait) {
    MultivariateTraitTree myTreeModel = trait.getTreeModel();
    final int tipCount = myTreeModel.getExternalNodeCount();
    double[] multiplierVect = new double[tipCount];
    NodeRef tempNode;
    for (int i = 0; i < tipCount; i++) {
        tempNode = myTreeModel.getExternalNode(i);
        multiplierVect[i] = Math.exp(-trait.getTimeScaledSelection(tempNode));
        tempNode = myTreeModel.getParent(tempNode);
        while (!myTreeModel.isRoot(tempNode)) {
            multiplierVect[i] = multiplierVect[i] * Math.exp(-trait.getTimeScaledSelection(tempNode));
            tempNode = myTreeModel.getParent(tempNode);
        }
    }
    return multiplierVect;
}

Example 5 with MultivariateTraitTree

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

the class MultivariateTraitUtils method getShiftContributionToMean.

private static double[] getShiftContributionToMean(NodeRef node, FullyConjugateMultivariateTraitLikelihood trait) {
    MultivariateTraitTree treeModel = trait.getTreeModel();
    double[] shiftContribution = new double[trait.dimTrait];
    if (!treeModel.isRoot(node)) {
        NodeRef parent = treeModel.getParent(node);
        double[] shiftContributionParent = getShiftContributionToMean(parent, trait);
        for (int i = 0; i < shiftContribution.length; ++i) {
            shiftContribution[i] = trait.getShiftForBranchLength(node)[i] + shiftContributionParent[i];
        }
    }
    return shiftContribution;
}