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

Example 36 with NodeRef

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

the class ExchangeOperator method getRandomNode.

private NodeRef getRandomNode(NodeRef[] nodes, NodeRef ref) {
    calcDistances(nodes, ref);
    double sum = 0;
    for (double distance : distances) {
        sum += 1.0 / distance;
    }
    double randomValue = MathUtils.nextDouble() * sum;
    NodeRef n = null;
    for (int i = 0; i < distances.length; i++) {
        randomValue -= 1.0 / distances[i];
        if (randomValue <= 0) {
            n = nodes[i];
            break;
        }
    }
    return n;
}

Example 37 with NodeRef

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

the class ExchangeOperator method wide.

/**
     * WARNING: Assumes strictly bifurcating tree.
     */
public boolean wide() {
    final int nodeCount = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    NodeRef i = root;
    while (root == i) {
        i = tree.getNode(MathUtils.nextInt(nodeCount));
    }
    NodeRef j = i;
    while (j == i || j == root) {
        j = tree.getNode(MathUtils.nextInt(nodeCount));
    }
    final NodeRef iP = tree.getParent(i);
    final NodeRef jP = tree.getParent(j);
    if ((iP != jP) && (i != jP) && (j != iP) && (tree.getNodeHeight(j) < tree.getNodeHeight(iP)) && (tree.getNodeHeight(i) < tree.getNodeHeight(jP))) {
        exchangeNodes(tree, i, j, iP, jP);
        // System.out.println("tries = " + tries+1);
        return true;
    }
    return false;
}

Example 38 with NodeRef

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

the class GibbsPruneAndRegraft method prunedGibbsProposal.

private double prunedGibbsProposal(Likelihood likelihood) {
    final int nodeCount = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    for (int i = 0; i < nodeCount; i++) {
        scores[i] = Double.NEGATIVE_INFINITY;
    }
    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, jP;
    final NodeRef iP = tree.getParent(i);
    final double heightIP = tree.getNodeHeight(iP);
    double sum = 0.0;
    double backwardLikelihood = calculateTreeLikelihood(likelihood, tree);
    int offset = (int) -backwardLikelihood;
    double backward = Math.exp(backwardLikelihood + offset);
    final NodeRef oldBrother = getOtherChild(tree, iP, i);
    final NodeRef oldGrandfather = tree.getParent(iP);
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((i != j) && (tree.getNodeHeight(j) < heightIP) && (heightIP < tree.getNodeHeight(jP)) && getNodeDistance(iP, jP) <= MAX_DISTANCE) {
                secondNodeIndices.add(n);
                pruneAndRegraft(tree, i, iP, j, jP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                probabilities.add(prob);
                scores[n] = prob;
                sum += prob;
                pruneAndRegraft(tree, i, iP, oldBrother, oldGrandfather);
            }
        }
    }
    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);
    pruneAndRegraft(tree, i, iP, j, jP);
    // now simulate the backward move
    double sumBackward = 0.0;
    final NodeRef newBrother = j;
    final NodeRef newGrandfather = jP;
    for (int n = 0; n < nodeCount; n++) {
        j = tree.getNode(n);
        if (j != root) {
            jP = tree.getParent(j);
            if ((i != j) && (tree.getNodeHeight(j) < heightIP) && (heightIP < tree.getNodeHeight(jP)) && getNodeDistance(iP, jP) <= MAX_DISTANCE) {
                if (scores[n] != Double.NEGATIVE_INFINITY) {
                    sumBackward += scores[n];
                } else {
                    pruneAndRegraft(tree, i, iP, j, jP);
                    double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                    sumBackward += prob;
                    pruneAndRegraft(tree, i, iP, newBrother, newGrandfather);
                    evaluate(likelihood, 1.0);
                }
            }
        }
    }
    double forward = probabilities.get(index);
    final double forwardProb = (forward / sum);
    final double backwardProb = (backward / (sumBackward));
    final double hastingsRatio = Math.log(backwardProb / forwardProb);
    return hastingsRatio;
}

Example 39 with NodeRef

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

the class NNI method doOperation.

/*
     * (non-Javadoc)
     * 
     * @see dr.inference.operators.SimpleMCMCOperator#doOperation()
     */
@Override
public double doOperation() {
    final int nNodes = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    NodeRef i;
    // get a random node where neither you or your father is the root
    do {
        i = tree.getNode(MathUtils.nextInt(nNodes));
    } while (root == i || tree.getParent(i) == root);
    // get parent node
    final NodeRef iParent = tree.getParent(i);
    // get parent of parent -> grant parent :)
    final NodeRef iGrandParent = tree.getParent(iParent);
    // get left child of grant parent -> uncle
    NodeRef iUncle = tree.getChild(iGrandParent, 0);
    // check if uncle == father
    if (iUncle == iParent) {
        // if so take right child -> sibling of father
        iUncle = tree.getChild(iGrandParent, 1);
    }
    // change the height of my father to be randomly between my uncle's
    // heights and my grandfather's height
    // this is necessary for the hastings ratio to do also if the uncle is
    // younger anyway
    final double heightGrandfather = tree.getNodeHeight(iGrandParent);
    final double heightUncle = tree.getNodeHeight(iUncle);
    final double minHeightFather = Math.max(heightUncle, tree.getNodeHeight(getOtherChild(tree, iParent, i)));
    final double heightI = tree.getNodeHeight(i);
    final double minHeightReverse = Math.max(heightI, tree.getNodeHeight(getOtherChild(tree, iParent, i)));
    double ran;
    do {
        ran = Math.random();
    } while (ran == 0.0 || ran == 1.0);
    // now calculate the new height for father between the height of the
    // uncle and the grandparent
    final double newHeightFather = minHeightFather + (ran * (heightGrandfather - minHeightFather));
    // set the new height for the father
    tree.setNodeHeight(iParent, newHeightFather);
    // double prForward = 1 / (heightGrandfather - minHeightFather);
    // double prBackward = 1 / (heightGrandfather - minHeightReverse);
    // hastings ratio = backward Prob / forward Prob
    final double hastingsRatio = Math.log((heightGrandfather - minHeightFather) / (heightGrandfather - minHeightReverse));
    // now change the nodes
    exchangeNodes(tree, i, iUncle, iParent, iGrandParent);
    tree.pushTreeChangedEvent(iParent);
    tree.pushTreeChangedEvent(iGrandParent);
    return hastingsRatio;
}

Example 40 with NodeRef

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

the class OldLatentLiabilityGibbsOperator method doOperation.

//    private boolean nodeGeoSpatialPriorExists(NodeRef node) {
//        return nodeGeoSpatialPrior != null && nodeGeoSpatialPrior.containsKey(treeModel.getNodeTaxon(node));
//    }
//
//    private boolean nodeMVNPriorExists(NodeRef node) {
//        return nodeMVNPrior != null && nodeMVNPrior.containsKey(treeModel.getNodeTaxon(node));
//    }
public double doOperation() {
    traitModel.redrawAncestralStates();
    NodeRef node = treeModel.getNode(MathUtils.nextInt(treeModel.getExternalNodeCount()));
    int tip = node.getNumber();
    // Draw truncated MVN using rejection sampling
    do {
    // Nothing
    } while (!liabilityLikelihood.validTraitForTip(tip));
    return 0;
}