use of dr.evolution.tree.Clade in project beast-mcmc by beast-dev.
the class AbstractCladeImportanceDistribution method getNonComplementaryClades.
/**
* creates a list with all clades but just the non-complementary ones.
* ((A,B),(C,D)) just {A,B,C,D} and {A,B} are inserted. {A,B} is
* complementary to {C,D}
*
* @param tree - the tree from which the clades are extracted
* @param node - the starting node. All clades below starting at this branch
* are added
*/
protected Clade getNonComplementaryClades(Tree tree, NodeRef node, List<Clade> parentClades, List<Clade> childClade) {
// create a new bit set for this clade
BitSet bits = new BitSet();
Clade c = null;
// check if the node is external
if (tree.isExternal(node)) {
// if so, the only taxon in the clade is I
int index = node.getNumber();
bits.set(index);
c = new Clade(bits, tree.getNodeHeight(node));
} else {
// otherwise, call all children and add its taxon together to one
// clade
NodeRef childNode = tree.getChild(node, 0);
// add just my first child to the list
// the second child is complementary to the first
Clade leftChild = getNonComplementaryClades(tree, childNode, parentClades, childClade);
bits.or(leftChild.getBits());
childNode = tree.getChild(node, 1);
// add just my first child to the list
// the second child is complementary to the first
Clade rightChild = getNonComplementaryClades(tree, childNode, parentClades, childClade);
bits.or(rightChild.getBits());
c = new Clade(bits, tree.getNodeHeight(node));
if (leftChild.getSize() >= 2) {
parentClades.add(c);
childClade.add(leftChild);
} else if (rightChild.getSize() >= 2) {
parentClades.add(c);
childClade.add(rightChild);
}
}
return c;
}
use of dr.evolution.tree.Clade in project beast-mcmc by beast-dev.
the class WeightedMultiplicativeBinary method calculateTreeProbabilityLog.
/**
* Calculates the probability of a given tree.
*
* @param tree - the tree to be analyzed
* @return estimated posterior probability in log
*/
private double calculateTreeProbabilityLog(Tree tree, HashMap<String, Integer> taxonMap) {
double prob = 0.0;
// calculate the number of possible splits
final double splits = Math.pow(2, tree.getExternalNodeCount() - 1) - 1;
List<Clade> clades = new ArrayList<Clade>();
List<Clade> parentClades = new ArrayList<Clade>();
// get clades contained in the tree
getClades(tree, tree.getRoot(), parentClades, clades, taxonMap);
// tree probability
for (Clade c : clades) {
// set the occurrences to epsilon
double occurrences = EPSILON;
if (cladeProbabilities.containsKey(c.getBits())) {
// if we observed this clade in the trace, add the occurrences
// to epsilon
double cladesInTreeSpace = getTrees(c.getBits().cardinality()) * getTrees(TAXA_COUNT - c.getBits().cardinality() + 1);
occurrences += (cladeProbabilities.get(c.getBits()).getSampleCount() / cladesInTreeSpace);
}
// multiply the conditional clade probability to the tree
// probability
prob += Math.log(occurrences / (samples + (splits * EPSILON)));
}
return prob;
}
use of dr.evolution.tree.Clade in project beast-mcmc by beast-dev.
the class WeightedMultiplicativeBinary method splitClade.
/*
* (non-Javadoc)
*
* @see
* dr.evolution.tree.ImportanceDistribution#splitClade(dr.evolution.tree
* .Clade, dr.evolution.tree.Clade[])
*/
public double splitClade(Clade parent, Clade[] children) {
// the number of all possible clades is 2^n with n the number of tips
// reduced by 2 because we wont consider the clades with all or no tips
// contained
// note: this time we consider each clade of a split separately with its
// own probability because every clade has a different chance for
// itself.
// #splits = 2^(n) - 1
final double splits = Math.pow(2, parent.getSize()) - 1;
double prob = 0;
double sum = 0.0;
List<Clade> childClades = getPossibleChildren(parent);
for (Clade child : childClades) {
sum += child.getSampleCount();
}
sum += EPSILON * splits;
double randomNumber = Math.random() * sum;
for (Clade child : childClades) {
randomNumber -= (child.getSampleCount() + EPSILON);
if (randomNumber < 0) {
children[0] = child;
double chance = (child.getSampleCount() + EPSILON) / samples;
// the other clade which would have resulted into the same split
BitSet secondChild = (BitSet) children[0].getBits().clone();
secondChild.xor(parent.getBits());
if (secondChild.cardinality() > 1) {
Clade counterClade = cladeProbabilities.get(secondChild);
if (counterClade != null) {
chance += (counterClade.getSampleCount() + EPSILON) / samples;
} else {
chance += EPSILON / samples;
}
prob = chance / 2.0;
} else {
prob = chance;
}
break;
}
}
// we take a clade which we haven't seen so far
if (randomNumber >= 0) {
// System.out.println("Random Clade");
BitSet newChild;
do {
do {
newChild = (BitSet) parent.getBits().clone();
int index = -1;
do {
index = newChild.nextSetBit(index + 1);
if (index > -1 && MathUtils.nextBoolean()) {
newChild.clear(index);
}
} while (index > -1);
} while (newChild.cardinality() == 0 || newChild.cardinality() == parent.getSize());
} while (cladeProbabilities.containsKey(newChild));
Clade randomClade = new Clade(newChild, 0.5);
children[0] = randomClade;
BitSet secondChild = (BitSet) children[0].getBits().clone();
secondChild.xor(parent.getBits());
if (cladeProbabilities.containsKey(secondChild)) {
children[1] = cladeProbabilities.get(secondChild);
} else {
children[1] = new Clade(secondChild, 0.5);
}
if (children[0].getSize() > 1 && children[1].getSize() > 1) {
prob = (children[0].getSampleCount() + children[1].getSampleCount() + (2.0 * EPSILON)) / (samples * 2.0);
} else {
if (children[0].getSize() > 1) {
prob = (children[0].getSampleCount() + EPSILON) / samples;
} else {
prob = (children[1].getSampleCount() + EPSILON) / samples;
}
}
} else {
BitSet secondChild = (BitSet) children[0].getBits().clone();
secondChild.xor(parent.getBits());
children[1] = cladeProbabilities.get(secondChild);
// children[1] = childClades.get(secondChild);
if (children[1] == null) {
children[1] = new Clade(secondChild, 0.5);
children[1].addHeight(0.5);
}
}
return Math.log(prob);
}
use of dr.evolution.tree.Clade in project beast-mcmc by beast-dev.
the class WeightedMultiplicativeBinary method calculateTreeProbabilityLogRecursive.
/**
* Calculates the probability of a given tree recursively.
*
* @param tree - the tree to be analyzed
* @param node - the node at which the subtree is rooted for which the
* probability has to be calculated
* @return estimated posterior probability in log
*/
private double calculateTreeProbabilityLogRecursive(Tree tree, NodeRef node) {
double prob = 0.0;
NodeRef leftChild = tree.getChild(node, 0);
NodeRef rightChild = tree.getChild(node, 1);
if (tree.isExternal(leftChild) && tree.isExternal(rightChild)) {
// both children are external nodes
return 0.0;
} else if (!tree.isExternal(leftChild) && !tree.isExternal(rightChild)) {
// both children are internal nodes
Clade leftSubclade = getClade(tree, leftChild);
Clade rightSubclade = getClade(tree, rightChild);
double sum = 0.0;
if (cladeProbabilities.containsKey(leftSubclade.getBits())) {
sum += (cladeProbabilities.get(leftSubclade.getBits()).getSampleCount() + EPSILON) / samples;
} else {
sum += EPSILON / samples;
}
if (cladeProbabilities.containsKey(rightSubclade.getBits())) {
sum += (cladeProbabilities.get(rightSubclade.getBits()).getSampleCount() + EPSILON) / samples;
} else {
sum += EPSILON / samples;
}
prob += Math.log(sum / 2.0);
prob += calculateTreeProbabilityLogRecursive(tree, leftChild);
prob += calculateTreeProbabilityLogRecursive(tree, rightChild);
return prob;
} else {
Clade leftSubclade = getClade(tree, leftChild);
Clade rightSubclade = getClade(tree, rightChild);
double sum = 0.0;
if (leftSubclade.getSize() > 1) {
if (cladeProbabilities.containsKey(leftSubclade.getBits())) {
sum += (cladeProbabilities.get(leftSubclade.getBits()).getSampleCount() + EPSILON) / samples;
} else {
sum += EPSILON / samples;
}
}
if (rightSubclade.getSize() > 1) {
if (cladeProbabilities.containsKey(rightSubclade.getBits())) {
sum += (cladeProbabilities.get(rightSubclade.getBits()).getSampleCount() + EPSILON) / samples;
} else {
sum += EPSILON / samples;
}
}
prob += Math.log(sum);
if (!tree.isExternal(leftChild)) {
prob += calculateTreeProbabilityLogRecursive(tree, leftChild);
}
if (!tree.isExternal(rightChild)) {
prob += calculateTreeProbabilityLogRecursive(tree, rightChild);
}
return prob;
}
}
use of dr.evolution.tree.Clade in project beast-mcmc by beast-dev.
the class WeightedMultiplicativeBinary method addTree.
/**
* increments the number of occurrences for all conditional clades
*
* @param tree - the tree to be added
*/
public void addTree(Tree tree, HashMap<String, Integer> taxonMap) {
samples++;
List<Clade> clades = new ArrayList<Clade>();
List<Clade> parentClades = new ArrayList<Clade>();
// get clades contained in the tree
getClades(tree, tree.getRoot(), parentClades, clades, taxonMap);
// increment the occurrences of the clade and the conditional clade
for (Clade c : clades) {
// increment the clade occurrences
if (cladeProbabilities.containsKey(c.getBits())) {
Clade tmp = cladeProbabilities.get(c.getBits());
tmp.addHeight(c.getHeight());
// frequency += cladeProbabilities.get(c);
} else {
// just to set the first value of the height value list
c.addHeight(c.getHeight());
cladeProbabilities.put(c.getBits(), c);
}
}
}
Aggregations