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

Example 41 with NodeRef

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

the class PrecisionMatrixGibbsOperator method incrementOuterProduct.

//    private void checkDiagonals(double[][] S) {
//        for (int i = 0; i < S.length; ++i) {
//            if (S[i][i] < 0.0) {
//                System.err.println("ERROR diag(S)\n" + new Matrix(S));
//                System.exit(-1);
//            }
//        }
//    }
private void incrementOuterProduct(double[][] S, NodeRef node) {
    if (!treeModel.isRoot(node)) {
        NodeRef parent = treeModel.getParent(node);
        double[] parentTrait = treeModel.getMultivariateNodeTrait(parent, traitName);
        double[] childTrait = treeModel.getMultivariateNodeTrait(node, traitName);
        double time = traitModel.getRescaledBranchLengthForPrecision(node);
        if (time > 0) {
            double sqrtTime = Math.sqrt(time);
            double[] delta = new double[dim];
            for (int i = 0; i < dim; i++) delta[i] = (childTrait[i] - parentTrait[i]) / sqrtTime;
            for (int i = 0; i < dim; i++) {
                // symmetric matrix,
                for (int j = i; j < dim; j++) S[j][i] = S[i][j] += delta[i] * delta[j];
            }
            // This assumes a *single* observation per tip
            numberObservations += 1;
        }
    }
    // recurse down tree
    for (int i = 0; i < treeModel.getChildCount(node); i++) incrementOuterProduct(S, treeModel.getChild(node, i));
}

Example 42 with NodeRef

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

the class RLCNarrowExchangeOperator method narrow.

/**
     * WARNING: Assumes strictly bifurcating tree.
     */
public void narrow() {
    final int nNodes = tree.getNodeCount();
    final NodeRef root = tree.getRoot();
    for (int tries = 0; tries < MAX_TRIES; ++tries) {
        NodeRef i = tree.getNode(MathUtils.nextInt(nNodes));
        while (root == i || tree.getParent(i) == root) {
            i = tree.getNode(MathUtils.nextInt(nNodes));
        }
        final NodeRef iParent = tree.getParent(i);
        final NodeRef iGrandParent = tree.getParent(iParent);
        NodeRef iUncle = tree.getChild(iGrandParent, 0);
        if (iUncle == iParent) {
            iUncle = tree.getChild(iGrandParent, 1);
        }
        assert tree.getNodeHeight(i) < tree.getNodeHeight(iGrandParent);
        if (tree.getNodeHeight(iUncle) < tree.getNodeHeight(iParent)) {
            NodeRef iSister = tree.getChild(iParent, 0);
            if (iSister == i) {
                iSister = tree.getChild(iParent, 1);
            }
            eupdate(i, iUncle, iParent, iGrandParent, iSister);
            tree.pushTreeChangedEvent(iParent);
            tree.pushTreeChangedEvent(iGrandParent);
            return;
        }
    }
    throw new RuntimeException("Couldn't find valid narrow move on this tree!!");
}

Example 43 with NodeRef

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

the class RLCNarrowExchangeOperator method eupdate.

/* exchange subtrees whose root are i and j */
private void eupdate(NodeRef i, NodeRef j, NodeRef iP, NodeRef jP, NodeRef iS) {
    tree.beginTreeEdit();
    tree.removeChild(iP, i);
    tree.removeChild(jP, j);
    tree.addChild(jP, i);
    tree.addChild(iP, j);
    tree.endTreeEdit();
    List<NodeRef> nodes = new ArrayList<NodeRef>();
    nodes.add(i);
    nodes.add(iP);
    nodes.add(j);
    nodes.add(jP);
    nodes.add(iS);
    NodeRef a = nodes.remove(MathUtils.nextInt(nodes.size()));
    NodeRef b = nodes.get(MathUtils.nextInt(nodes.size()));
    //swap traits in these two nodes
    double changedA = tree.getNodeTrait(a, "trait");
    double changedB = tree.getNodeTrait(a, "trait");
    tree.setNodeTrait(a, "trait", changedB);
    tree.setNodeTrait(b, "trait", changedA);
}

Example 44 with NodeRef

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

the class RandomWalkIntegerSetSizeWeightedOperator method computeSampleWeights.

private void computeSampleWeights() {
    TreeModel tree = msatSampleTreeModel.getTreeModel();
    int intNodeCount = tree.getInternalNodeCount();
    int extNodeCount = tree.getExternalNodeCount();
    weights = new double[intNodeCount];
    for (int i = 0; i < intNodeCount; i++) {
        NodeRef node = tree.getNode(i + extNodeCount);
        int lcState = msatSampleTreeModel.getNodeValue(tree.getChild(node, 0));
        int rcState = msatSampleTreeModel.getNodeValue(tree.getChild(node, 1));
        weights[i] = Math.abs(lcState - rcState) + baseSetSize;
    }
}

Example 45 with NodeRef

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

the class GibbsPruneAndRegraft method gibbsProposal.

private double gibbsProposal(Likelihood likelihood) {
    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, 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))) {
                secondNodeIndices.add(n);
                pruneAndRegraft(tree, i, iP, j, jP);
                double prob = Math.exp(calculateTreeLikelihood(likelihood, tree) + offset);
                probabilities.add(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);
    double forward = probabilities.get(index);
    double forwardProb = (forward / sum);
    double backwardProb = (backward / (sum - forward + backward));
    final double hastingsRatio = Math.log(backwardProb / forwardProb);
    return hastingsRatio;
}