Examples with Node beast.evolution.tree.Node used on opensource projects

Example 41 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodAmbiguities.

@Test
public void testLikelihoodAmbiguities() throws Exception {
    Locus locus = new Locus("locus", getAlignment());
    // ConversionGraph
    ConversionGraph acg = new ConversionGraph();
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
    acg.assignFrom(tree);
    acg.initByName("locus", locus);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", jc);
    // Likelihood
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    acg.setEverythingDirty(true);
    double logP = argLikelihood.calculateLogP();
    double logPtrue = argLikelihoodSlow.calculateLogP();
    double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    System.out.println("Test 1.  Truth: " + logPtrue + " Value: " + logP);
    assertTrue(relativeDiff < 1e-14);
    // Add a single recombination event
    Node node1 = acg.getExternalNodes().get(0);
    Node node2 = node1.getParent();
    double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
    double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    int startLocus = 100;
    int endLocus = 200;
    Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb1);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
    // Add another recombination event
    node1 = acg.getExternalNodes().get(0);
    node2 = acg.getNode(20);
    height1 = 0.75 * (node1.getHeight() + node1.getParent().getHeight());
    height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    startLocus = 250;
    endLocus = 300;
    Conversion recomb2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb2);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
}

Example 42 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class SimulatedAlignment method traverse.

/**
 * Traverse a marginal tree simulating a region of the sequence alignment
 * down it.
 *
 * @param node Node of the marginal tree
 * @param parentSequence Sequence at the parent node in the marginal tree
 * @param categories Mapping from sites to categories
 * @param transitionProbs
 * @param regionAlignment
 */
private void traverse(Node node, int[] parentSequence, int[] categories, double[][] transitionProbs, int[][] regionAlignment) {
    for (Node child : node.getChildren()) {
        // Calculate transition probabilities
        for (int i = 0; i < siteModel.getCategoryCount(); i++) {
            siteModel.getSubstitutionModel().getTransitionProbabilities(child, node.getHeight(), child.getHeight(), siteModel.getRateForCategory(i, child), transitionProbs[i]);
        }
        // Draw characters on child sequence
        int[] childSequence = new int[parentSequence.length];
        int nStates = dataType.getStateCount();
        double[] charProb = new double[nStates];
        for (int i = 0; i < childSequence.length; i++) {
            int category = categories[i];
            System.arraycopy(transitionProbs[category], parentSequence[i] * nStates, charProb, 0, nStates);
            childSequence[i] = Randomizer.randomChoicePDF(charProb);
        }
        if (child.isLeaf()) {
            System.arraycopy(childSequence, 0, regionAlignment[child.getNr()], 0, childSequence.length);
        } else {
            traverse(child, childSequence, categories, transitionProbs, regionAlignment);
        }
    }
}

Example 43 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGOperator method connectEdge.

/**
 * Connect edge <node, node.parent> above destEdgeBase, forming new
 * edge <destEdgeBase, node.parent> and <node.parent, destEdgeBase.parent>.
 * All conversions above destEdgeBase that are older than destTime
 * are transferred to the new edge above node.parent.
 *
 * Conversions on edge above node are not modified.
 *
 * @param node base of edge to attach
 * @param destEdgeBase base of edge to be bisected
 * @param destTime height at which bisection will occur
 */
protected void connectEdge(Node node, Node destEdgeBase, double destTime) {
    if (node.isRoot())
        throw new IllegalArgumentException("Programmer error: " + "root argument passed to connectEdge().");
    Node parent = node.getParent();
    if (destEdgeBase.isRoot()) {
        parent.addChild(destEdgeBase);
    } else {
        Node grandParent = destEdgeBase.getParent();
        grandParent.removeChild(destEdgeBase);
        grandParent.addChild(parent);
        parent.addChild(destEdgeBase);
    }
    parent.setHeight(destTime);
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode1() == destEdgeBase && conv.getHeight1() > destTime)
                conv.setNode1(parent);
            if (conv.getNode2() == destEdgeBase && conv.getHeight2() > destTime)
                conv.setNode2(parent);
        }
    }
}

Example 44 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGScaler method proposal.

@Override
public double proposal() {
    // Keep track of number of positively scaled elements minus
    // negatively scaled elements.
    int count = 0;
    // Choose scaling factor:
    double f = scaleParam + Randomizer.nextDouble() * (1.0 / scaleParam - scaleParam);
    // Scale clonal frame:
    if (rootOnly) {
        acg.getRoot().setHeight(acg.getRoot().getHeight() * f);
        count += 1;
    } else {
        for (Node node : acg.getInternalNodes()) {
            node.setHeight(node.getHeight() * f);
            count += 1;
        }
    }
    // Scale recombinant edges:
    for (Locus locus : acg.getLoci()) {
        for (Conversion recomb : acg.getConversions(locus)) {
            if (!rootOnly || recomb.getNode1().getParent().isRoot()) {
                recomb.setHeight1(recomb.getHeight1() * f);
                count += 1;
            }
            if (!rootOnly || recomb.getNode2().isRoot() || recomb.getNode2().getParent().isRoot()) {
                recomb.setHeight2(recomb.getHeight2() * f);
                count += 1;
            }
            if (recomb.getHeight1() < recomb.getNode1().getHeight() || recomb.getHeight2() < recomb.getNode2().getHeight() || recomb.getHeight1() > recomb.getHeight2())
                return Double.NEGATIVE_INFINITY;
        }
    }
    // Check for illegal node heights:
    if (rootOnly) {
        for (Node node : acg.getRoot().getChildren()) {
            if (node.getHeight() > node.getParent().getHeight())
                return Double.NEGATIVE_INFINITY;
        }
    } else {
        for (Node node : acg.getExternalNodes()) {
            if (node.getHeight() > node.getParent().getHeight())
                return Double.NEGATIVE_INFINITY;
        }
    }
    for (RealParameter param : parametersInput.get()) {
        try {
            param.startEditing(null);
            param.scale(f);
        } catch (Exception e) {
            // bounds.  Needs to change!
            return Double.NEGATIVE_INFINITY;
        }
        count += param.getDimension();
    }
    for (RealParameter paramInv : parametersInverseInput.get()) {
        try {
            paramInv.startEditing(null);
            paramInv.scale(1.0 / f);
        } catch (Exception e) {
            // bounds.  Needs to change!
            return Double.NEGATIVE_INFINITY;
        }
        count -= paramInv.getDimension();
    }
    // Return log of Hastings ratio:
    return (count - 2) * Math.log(f);
}

Example 45 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class AddRemoveDetour method removeDetour.

/**
 * Detour deletion variant of move.
 *
 * @return log HGF
 */
private double removeDetour() {
    double logHGF = 0.0;
    // Choose non-root detour edge at random
    Node detour = acg.getNode(Randomizer.nextInt(acg.getNodeCount() - 1));
    logHGF -= Math.log(1.0 / (acg.getNodeCount() - 1));
    List<Conversion> convApotentials = new ArrayList<>();
    List<Conversion> convBpotentials = new ArrayList<>();
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode2() == detour && conv.getNode1() != detour)
                convApotentials.add(conv);
            if (conv.getNode1() == detour && conv.getNode2() != detour)
                convBpotentials.add(conv);
        }
    }
    if (convApotentials.isEmpty() || convBpotentials.isEmpty())
        return Double.NEGATIVE_INFINITY;
    Conversion convA = convApotentials.get(Randomizer.nextInt(convApotentials.size()));
    Conversion convB = convBpotentials.get(Randomizer.nextInt(convBpotentials.size()));
    if (convA.getHeight2() > convB.getHeight1())
        return Double.NEGATIVE_INFINITY;
    logHGF -= Math.log(1.0 / (convApotentials.size() * convBpotentials.size()));
    double tLowerBound = convA.getHeight1();
    double tUpperBound = convB.getHeight2();
    Conversion conv = new Conversion();
    conv.setNode1(convA.getNode1());
    conv.setHeight1(convA.getHeight1());
    conv.setNode2(convB.getNode2());
    conv.setHeight2(convB.getHeight2());
    conv.setStartSite(convA.getStartSite());
    conv.setEndSite(convA.getEndSite());
    conv.setLocus(convA.getLocus());
    acg.deleteConversion(convA);
    acg.deleteConversion(convB);
    acg.addConversion(conv);
    logHGF += Math.log(1.0 / acg.getTotalConvCount()) + Math.log(1.0 / (acg.getNodeCount() - 1)) + Math.log(2.0) - 2.0 * Math.log(tUpperBound - tLowerBound) + getAffectedRegionProb(convB);
    return logHGF;
}