Search in sources :

Example 21 with Clade

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;
}
Also used : NodeRef(dr.evolution.tree.NodeRef) BitSet(java.util.BitSet) Clade(dr.evolution.tree.Clade)

Example 22 with Clade

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;
}
Also used : Clade(dr.evolution.tree.Clade)

Example 23 with Clade

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);
}
Also used : Clade(dr.evolution.tree.Clade)

Example 24 with Clade

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;
    }
}
Also used : NodeRef(dr.evolution.tree.NodeRef) Clade(dr.evolution.tree.Clade)

Example 25 with Clade

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);
        }
    }
}
Also used : Clade(dr.evolution.tree.Clade)

Aggregations

Clade (dr.evolution.tree.Clade)25 NodeRef (dr.evolution.tree.NodeRef)13 BitSet (java.util.BitSet)7 InvalidTreeException (dr.evolution.tree.MutableTree.InvalidTreeException)1