use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.
the class LogCombiner method rescaleTree.
private void rescaleTree(Tree tree, double scale) {
if (tree instanceof MutableTree) {
MutableTree mutableTree = (MutableTree) tree;
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef node = tree.getNode(i);
if (node != tree.getRoot()) {
double length = tree.getBranchLength(node);
mutableTree.setBranchLength(node, length * scale);
}
}
} else {
throw new IllegalArgumentException("Tree not mutable");
}
}
use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.
the class DensityMap method addTree.
public void addTree(Tree tree, double sampleTime, String attributeName1, String attributeName2) {
checkCalibration();
double[][] variance = null;
Object[] obj = (Object[]) tree.getAttribute(MultivariateDiffusionModel.PRECISION_TREE_ATTRIBUTE);
if (obj != null) {
variance = new Matrix(MatrixParameter.parseFromSymmetricDoubleArray(obj).getParameterAsMatrix()).inverse().toComponents();
}
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef node = tree.getNode(i);
if (node != tree.getRoot()) {
NodeRef parent = tree.getParent(node);
double t1 = tree.getNodeHeight(node);
double t2 = tree.getNodeHeight(parent);
if (t1 <= sampleTime && t2 >= sampleTime) {
Double valueX1 = transform((Double) tree.getNodeAttribute(node, attributeName1));
Double valueY1 = transform((Double) tree.getNodeAttribute(node, attributeName2));
Double valueX2 = transform((Double) tree.getNodeAttribute(parent, attributeName1));
Double valueY2 = transform((Double) tree.getNodeAttribute(parent, attributeName2));
if (valueX1 != null && valueY1 != null && valueX2 != null && valueY2 != null) {
addPoint(sampleTime, t1, t2, valueX1, valueY1, valueX2, valueY2, variance);
}
}
}
}
}
use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.
the class DensityMap method calibrate.
public void calibrate(Tree tree, String attributeName1, String attributeName2) {
boolean foundAttribute1 = false;
boolean foundAttribute2 = false;
jointDensity = true;
if (isCalibrated) {
throw new RuntimeException("Already calibrated");
}
// }
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef node = tree.getNode(i);
if (node != tree.getRoot()) {
Double value = (Double) tree.getNodeAttribute(node, attributeName1);
if (value != null) {
value = transform(value);
if (value < minX)
minX = value;
if (value > maxX)
maxX = value;
foundAttribute1 = true;
}
value = (Double) tree.getNodeAttribute(node, attributeName2);
if (value != null) {
value = transform(value);
if (value < minY)
minY = value;
if (value > maxY)
maxY = value;
foundAttribute2 = true;
}
}
}
if (!foundAttribute1) {
throw new RuntimeException("Can't find any attributes, " + attributeName1 + ", in tree " + tree.getId());
}
if (!foundAttribute2) {
throw new RuntimeException("Can't find any attributes, " + attributeName2 + ", in tree " + tree.getId());
}
// System.err.printf("Calibrated: minY = %3.2f, maxY = %3.2f, minX = %3.2f, maxX = %3.2f\n",minY,maxY,minX,maxX);
// System.exit(-1);
}
use of dr.evolution.tree.NodeRef in project beast-mcmc by beast-dev.
the class TaxaOriginTrait method getTraitName.
private String getTraitName() {
FlexibleTree tree = trees[0];
String traitName = null;
for (int i = 0; i < tree.getExternalNodeCount(); i++) {
NodeRef node = tree.getExternalNode(i);
for (String taxaName : taxaNames) {
if (tree.getNodeTaxon(node).getId().equals(taxaName)) {
String attributeValue = (String) tree.getNodeAttribute(node, attributeName);
if (traitName != null && !traitName.equals(attributeValue)) {
throw new RuntimeException("Not all taxa given have the same trait value");
} else {
traitName = attributeValue;
}
}
}
}
return traitName;
}
use of dr.evolution.tree.NodeRef 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;
}
}
Aggregations