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

Example 51 with NodeRef

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

the class RateSampleOperator method sampleSubtree.

void sampleSubtree(NodeRef parent) {
    int nbChildren = tree.getChildCount(parent);
    for (int c = 0; c < nbChildren; c++) {
        final NodeRef node = tree.getChild(parent, c);
        rateEvolution.sampleRate(node);
        sampleSubtree(node);
    }
}

Example 52 with NodeRef

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

the class RateVarianceScaleOperator method doOperation.

/**
     * scale the rates of a subtree and return the hastings ratio.
     */
public final double doOperation() {
    final double scale = (scaleFactor + (MathUtils.nextDouble() * ((1.0 / scaleFactor) - scaleFactor)));
    //Scale the variance
    double oldValue = variance.getParameterValue(0);
    double newValue = scale * oldValue;
    double logq = -Math.log(scale);
    final Bounds<Double> bounds = variance.getBounds();
    if (newValue < bounds.getLowerLimit(0) || newValue > bounds.getUpperLimit(0)) {
        //            throw new OperatorFailedException("proposed value outside boundaries");
        return Double.NEGATIVE_INFINITY;
    }
    variance.setParameterValue(0, newValue);
    //Scale the rates of the tree accordingly
    NodeRef root = tree.getRoot();
    final int index = root.getNumber();
    List<NodeRef> listNode = new ArrayList<NodeRef>();
    getSubtree(listNode, tree.getNode(index));
    final double rateScale = Math.sqrt(scale);
    for (NodeRef node : listNode) {
        oldValue = tree.getNodeRate(node);
        newValue = oldValue * rateScale;
        tree.setNodeRate(node, newValue);
    }
    //  According to the hastings ratio in the scale Operator
    logq += (listNode.size() - 2) * Math.log(rateScale);
    return logq;
}

Example 53 with NodeRef

use of dr.evolution.tree.NodeRef 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 54 with NodeRef

use of dr.evolution.tree.NodeRef 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 55 with NodeRef

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

the class ImportanceSubtreeSwap method swap.

/* exchange subtrees whose root are i and j */
private void swap(TreeModel tree, NodeRef i, NodeRef j) {
    NodeRef iP = tree.getParent(i);
    NodeRef jP = tree.getParent(j);
    tree.removeChild(iP, i);
    tree.removeChild(jP, j);
    tree.addChild(jP, i);
    tree.addChild(iP, j);
}