Examples with InvalidTreeException dr.evolution.tree.MutableTree.InvalidTreeException used on opensource projects

Example 1 with InvalidTreeException

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

the class ImportancePruneAndRegraft method importancePruneAndRegraft.

private double importancePruneAndRegraft() {
    final int nodeCount = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    NodeRef i;
    do {
        int indexI = MathUtils.nextInt(nodeCount);
        i = tree.getNode(indexI);
    } while (root == i || tree.getParent(i) == root);
    List<Integer> secondNodeIndices = new ArrayList<Integer>();
    List<Double> probabilities = new ArrayList<Double>();
    NodeRef j, iP, jP;
    iP = tree.getParent(i);
    double iParentHeight = tree.getNodeHeight(iP);
    double sum = 0.0;
    double backwardLikelihood = calculateTreeProbability(tree);
    int offset = (int) -backwardLikelihood;
    double backward = Math.exp(backwardLikelihood + offset);
    final NodeRef oldBrother = getOtherChild(tree, iP, i);
    final NodeRef oldGrandfather = tree.getParent(iP);
    tree.beginTreeEdit();
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((iP != jP) && (tree.getNodeHeight(j) < iParentHeight && iParentHeight < tree.getNodeHeight(jP))) {
                secondNodeIndices.add(n);
                pruneAndRegraft(tree, i, iP, j, jP);
                double prob = Math.exp(calculateTreeProbability(tree) + offset);
                probabilities.add(prob);
                sum += prob;
                pruneAndRegraft(tree, i, iP, oldBrother, oldGrandfather);
            }
        }
    }
    double ran = Math.random() * sum;
    int index = 0;
    while (ran > 0.0) {
        ran -= probabilities.get(index);
        index++;
    }
    index--;
    j = tree.getNode(secondNodeIndices.get(index));
    jP = tree.getParent(j);
    if (iP != jP) {
        pruneAndRegraft(tree, i, iP, j, jP);
        tree.pushTreeChangedEvent(i);
    }
    tree.endTreeEdit();
    // AR - not sure whether this check is necessary
    try {
        tree.checkTreeIsValid();
    } catch (InvalidTreeException e) {
        throw new RuntimeException(e.getMessage());
    }
    double forward = probabilities.get(index);
    // tree.pushTreeChangedEvent(jP);
    // tree.pushTreeChangedEvent(oldGrandfather);
    tree.pushTreeChangedEvent(i);
    double forwardProb = (forward / sum);
    double backwardProb = (backward / (sum - forward + backward));
    final double hastingsRatio = Math.log(backwardProb / forwardProb);
    return hastingsRatio;
}

Example 2 with InvalidTreeException

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

the class ImportanceSubtreeSwap method importanceExchange.

/**
     * WARNING: Assumes strictly bifurcating tree.
     *
     * @throws InvalidTreeException
     */
private double importanceExchange() {
    final int nodeCount = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    NodeRef i;
    int indexI;
    int indexJ;
    do {
        indexI = MathUtils.nextInt(nodeCount);
        i = tree.getNode(indexI);
    } while (root == i || (tree.getParent(i) == root && tree.getNodeHeight(i) > tree.getNodeHeight(getOtherChild(tree, tree.getParent(i), i))));
    List<Integer> secondNodeIndices = new ArrayList<Integer>();
    List<Double> probabilities = new ArrayList<Double>();
    NodeRef j, iP, jP;
    iP = tree.getParent(i);
    double sum = 0.0;
    double backward = calculateTreeProbability(tree);
    int offset = (int) -backward;
    backward = Math.exp(backward + offset);
    tree.beginTreeEdit();
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((iP != jP) && (i != jP) && (j != iP) && (tree.getNodeHeight(j) < tree.getNodeHeight(iP)) && (tree.getNodeHeight(i) < tree.getNodeHeight(jP))) {
                secondNodeIndices.add(n);
                swap(tree, tree.getNode(indexI), tree.getNode(n));
                double prob = Math.exp(calculateTreeProbability(tree) + offset);
                probabilities.add(prob);
                swap(tree, tree.getNode(indexI), tree.getNode(n));
                sum += prob;
            }
        }
    }
    double ran = Math.random() * sum;
    int index = 0;
    while (ran > 0.0) {
        ran -= probabilities.get(index);
        index++;
    }
    index--;
    j = tree.getNode(secondNodeIndices.get(index));
    jP = tree.getParent(j);
    // *******************************************
    // assuming we would have chosen j first
    double sumForward2 = 0.0;
    NodeRef k, kP;
    indexJ = secondNodeIndices.get(index);
    for (int n = 0; n < nodeCount; n++) {
        k = tree.getNode(n);
        if (k != root) {
            kP = tree.getParent(k);
            if ((jP != kP) && (j != kP) && (k != jP) && (tree.getNodeHeight(k) < tree.getNodeHeight(jP)) && (tree.getNodeHeight(j) < tree.getNodeHeight(kP))) {
                swap(tree, tree.getNode(indexJ), tree.getNode(n));
                double prob = Math.exp(calculateTreeProbability(tree) + offset);
                sumForward2 += prob;
                swap(tree, tree.getNode(indexJ), tree.getNode(n));
            }
        }
    }
    swap(tree, i, j);
    double forward = probabilities.get(index);
    iP = tree.getParent(i);
    double sumBackward = 0.0;
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((iP != jP) && (i != jP) && (j != iP) && (tree.getNodeHeight(j) < tree.getNodeHeight(iP)) && (tree.getNodeHeight(i) < tree.getNodeHeight(jP))) {
                swap(tree, tree.getNode(indexI), tree.getNode(n));
                double prob = Math.exp(calculateTreeProbability(tree) + offset);
                sumBackward += prob;
                swap(tree, tree.getNode(indexI), tree.getNode(n));
            }
        }
    }
    // *******************************************
    // assuming we would have chosen j first
    double sumBackward2 = 0.0;
    j = tree.getNode(secondNodeIndices.get(index));
    jP = tree.getParent(j);
    for (int n = 0; n < nodeCount; n++) {
        k = tree.getNode(n);
        if (k != root) {
            kP = tree.getParent(k);
            if ((jP != kP) && (j != kP) && (k != jP) && (tree.getNodeHeight(k) < tree.getNodeHeight(jP)) && (tree.getNodeHeight(j) < tree.getNodeHeight(kP))) {
                swap(tree, tree.getNode(indexJ), tree.getNode(n));
                double prob = Math.exp(calculateTreeProbability(tree) + offset);
                sumBackward2 += prob;
                swap(tree, tree.getNode(indexJ), tree.getNode(n));
            }
        }
    }
    tree.endTreeEdit();
    // AR - not sure whether this check is necessary
    try {
        tree.checkTreeIsValid();
    } catch (InvalidTreeException e) {
        throw new RuntimeException(e.getMessage());
    }
    double forwardProb = (forward / sum) + (forward / sumForward2);
    double backwardProb = (backward / sumBackward) + (backward / sumBackward2);
    double hastingsRatio = Math.log(backwardProb / forwardProb);
    return hastingsRatio;
}

Example 3 with InvalidTreeException

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

the class AbstractImportanceDistributionOperator method doImportanceDistributionOperation.

protected double doImportanceDistributionOperation(Likelihood likelihood) {
    final NodeRef root = tree.getRoot();
    BitSet all = new BitSet();
    all.set(0, (tree.getNodeCount() + 1) / 2);
    Clade rootClade = new Clade(all, tree.getNodeHeight(root));
    internalNodes.clear();
    fillInternalNodes(root);
    // remove the root
    internalNodes.poll();
    externalNodes.clear();
    fillExternalNodes(root);
    double prob;
    double back = probabilityEstimater.getTreeProbability(tree);
    try {
        tree.beginTreeEdit();
        List<Clade> originalClades = new ArrayList<Clade>();
        extractClades(tree, tree.getRoot(), originalClades, null);
        double[] originalNodeHeights = getAbsoluteNodeHeights(originalClades);
        Arrays.sort(originalNodeHeights);
        back += getChanceForNodeHeights(originalNodeHeights);
        prob = createTree(root, rootClade);
        assignDummyHeights(root);
        //			assignCladeHeights(tree.getRoot(), originalClades, null);
        //			double[] originalNodeHeights = getAbsoluteNodeHeights(originalClades);
        //			Arrays.sort(originalNodeHeights);
        //			prob += setMissingNodeHeights(tree.getChild(tree.getRoot(),0));
        //			prob += setMissingNodeHeights(tree.getChild(tree.getRoot(),1));
        prob += setNodeHeights(originalNodeHeights);
        //			List<Clade> newClades = new ArrayList<Clade>();
        //			extractClades(tree, tree.getRoot(), newClades, null);
        tree.endTreeEdit();
        tree.checkTreeIsValid();
    } catch (InvalidTreeException e) {
        throw new RuntimeException(e.getMessage());
    }
    tree.pushTreeChangedEvent(root);
    return back - prob;
}