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

Example 46 with NodeRef

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

the class GibbsSubtreeSwap method prunedWide.

/**
     * WARNING: Assumes strictly bifurcating tree.
     *
     * @throws InvalidTreeException
     */
public double prunedWide(Likelihood likelihood) {
    final int nodeCount = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    NodeRef i;
    int indexI;
    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, jP;
    NodeRef iP = tree.getParent(i);
    double heightIP = tree.getNodeHeight(iP);
    double heightI = tree.getNodeHeight(i);
    double sum = 0.0;
    double backward = calculateTreeLikelihood(likelihood, tree);
    int offset = (int) -backward;
    backward = Math.exp(backward + offset);
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((iP != jP) && (tree.getNodeHeight(j) < heightIP) && (heightI < tree.getNodeHeight(jP)) && getNodeDistance(iP, jP) <= MAX_DISTANCE) {
                secondNodeIndices.add(n);
                swap(tree, i, j, iP, jP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                probabilities.add(prob);
                swap(tree, i, j, jP, iP);
                sum += prob;
            }
        }
    }
    if (sum <= 1E-100) {
        // neglected
        throw new RuntimeException("Couldn't find another proposal with a decent likelihood.");
    }
    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);
    double heightJP = tree.getNodeHeight(jP);
    double heightJ = tree.getNodeHeight(j);
    // *******************************************
    // assuming we would have chosen j first
    double sumForward2 = 0.0;
    NodeRef k, kP;
    for (int n = 0; n < nodeCount; n++) {
        k = tree.getNode(n);
        if (k != root) {
            kP = tree.getParent(k);
            if ((jP != kP) && (tree.getNodeHeight(k) < heightJP) && (heightJ < tree.getNodeHeight(kP)) && getNodeDistance(kP, jP) <= MAX_DISTANCE) {
                swap(tree, j, k, jP, kP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                sumForward2 += prob;
                swap(tree, j, k, kP, jP);
            }
        }
    }
    swap(tree, i, j, iP, jP);
    double forward = probabilities.get(index);
    iP = jP;
    heightIP = heightJP;
    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) && (tree.getNodeHeight(j) < heightIP) && (heightI < tree.getNodeHeight(jP)) && getNodeDistance(iP, jP) <= MAX_DISTANCE) {
                swap(tree, i, j, iP, jP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                sumBackward += prob;
                swap(tree, i, j, jP, iP);
            }
        }
    }
    // *******************************************
    // assuming we would have chosen j first
    double sumBackward2 = 0.0;
    j = tree.getNode(secondNodeIndices.get(index));
    jP = tree.getParent(j);
    heightJP = tree.getNodeHeight(jP);
    heightJ = tree.getNodeHeight(j);
    for (int n = 0; n < nodeCount; n++) {
        k = tree.getNode(n);
        if (k != root) {
            kP = tree.getParent(k);
            if ((jP != kP) && (tree.getNodeHeight(k) < heightJP) && (heightJ < tree.getNodeHeight(kP)) && getNodeDistance(kP, jP) <= MAX_DISTANCE) {
                swap(tree, j, k, jP, kP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                sumBackward2 += prob;
                swap(tree, j, k, kP, jP);
            }
        }
    }
    double forwardProb = (forward / sum) + (forward / sumForward2);
    double backwardProb = (backward / sumBackward) + (backward / sumBackward2);
    double hastingsRatio = Math.log(backwardProb / forwardProb);
    return hastingsRatio;
}

Example 47 with NodeRef

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

the class ImportanceNarrowExchange method getNode.

private int getNode() {
    traverseTree(tree.getRoot());
    totalWeight = 0;
    for (int k = 0; k < tree.getInternalNodeCount(); ++k) {
        final NodeRef node = tree.getInternalNode(k);
        final double w = nodeWeight(node);
        weights[node.getNumber()] = w;
        if (DEBUG > 5 && w > 0) {
            System.out.println("" + w + " " + TreeUtils.uniqueNewick(tree, node));
        }
        totalWeight += w;
    }
    double r = MathUtils.nextDouble() * totalWeight;
    for (int k = 0; k < tree.getInternalNodeCount(); ++k) {
        final NodeRef node = tree.getInternalNode(k);
        final int nodeIndex = node.getNumber();
        r -= weights[nodeIndex];
        if (r < 0) {
            if (DEBUG > 0) {
                System.out.println("" + weights[nodeIndex] + "/" + totalWeight + " " + TreeUtils.uniqueNewick(tree, node));
            }
            return k;
        }
    }
    //assert false;
    return -1;
}

Example 48 with NodeRef

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

the class ImportanceNarrowExchange method nodeWeight.

private double nodeWeight(final NodeRef node) {
    final NodeRef ch0 = tree.getChild(node, 0);
    final NodeRef ch1 = tree.getChild(node, 1);
    if (tree.isExternal(ch0) && tree.isExternal(ch1)) {
        return 0;
    }
    final boolean leftSubtree = tree.getNodeHeight(ch0) < tree.getNodeHeight(ch1);
    final int st0 = nodeCounts[(leftSubtree ? ch0 : ch1).getNumber()];
    final int st1 = nodeCounts[tree.getChild(leftSubtree ? ch1 : ch0, 0).getNumber()];
    final int st2 = nodeCounts[tree.getChild(leftSubtree ? ch1 : ch0, 1).getNumber()];
    final double w = (epsilon + st0) * (epsilon + st1) + (epsilon + st0) * (epsilon + st2) + (epsilon + st1) * (epsilon + st2) - 3 * epsilon * epsilon;
    return w;
}

Example 49 with NodeRef

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

the class RateExchangeOperator method doOperation.

/**
     * Do a probablistic subtree slide move.
     *
     * @return the log-transformed hastings ratio
     */
public double doOperation() {
    NodeRef node0 = tree.getInternalNode(MathUtils.nextInt(tree.getInternalNodeCount()));
    NodeRef node1 = tree.getChild(node0, 0);
    NodeRef node2 = tree.getChild(node0, 1);
    if (swapRates) {
        if (swapAtRoot) {
            double[] rates = new double[] { tree.getNodeRate(node0), tree.getNodeRate(node1), tree.getNodeRate(node2) };
            int r1 = MathUtils.nextInt(3);
            tree.setNodeRate(node0, rates[r1]);
            // swap down the top trait
            rates[r1] = rates[2];
            int r2 = MathUtils.nextInt(2);
            tree.setNodeRate(node1, rates[r2]);
            // swap down the top trait
            rates[r2] = rates[1];
            tree.setNodeRate(node2, rates[0]);
        } else {
            // just swap the two child rates...
            double tmp = tree.getNodeRate(node1);
            tree.setNodeRate(node1, tree.getNodeRate(node2));
            tree.setNodeRate(node2, tmp);
        }
    }
    if (swapTraits) {
        if (swapAtRoot) {
            double[] traits = new double[] { tree.getNodeTrait(node0, TRAIT), tree.getNodeTrait(node1, TRAIT), tree.getNodeTrait(node2, TRAIT) };
            int r1 = MathUtils.nextInt(3);
            tree.setNodeTrait(node0, TRAIT, traits[r1]);
            // swap down the top trait
            traits[r1] = traits[2];
            int r2 = MathUtils.nextInt(2);
            tree.setNodeTrait(node1, TRAIT, traits[r2]);
            // swap down the top trait
            traits[r2] = traits[1];
            tree.setNodeTrait(node2, TRAIT, traits[0]);
        } else {
            // just swap the two child traits...
            double tmp = tree.getNodeTrait(node1, TRAIT);
            tree.setNodeTrait(node1, TRAIT, tree.getNodeTrait(node2, TRAIT));
            tree.setNodeTrait(node2, TRAIT, tmp);
        }
    }
    // If the node is not the root, do a uniform pick of its height
    if (!tree.isRoot(node0) && moveHeight) {
        double lower = tree.getNodeHeightLower(node0);
        double upper = tree.getNodeHeightUpper(node0);
        double newValue = (MathUtils.nextDouble() * (upper - lower)) + lower;
        tree.setNodeHeight(node0, newValue);
    }
    return 0.0;
}

Example 50 with NodeRef

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

the class RateSampleOperator method sampleNode.

void sampleNode(NodeRef parent) {
    int nbChildren = tree.getChildCount(parent);
    if (nbChildren > 0) {
        final int c = MathUtils.nextInt(nbChildren);
        final NodeRef node = tree.getChild(parent, c);
        rateEvolution.sampleRate(node);
    }
}