Examples with Locus bacter.Locus used on opensource projects

Example 16 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodUsingSimulatedData.

@Test
public void testLikelihoodUsingSimulatedData() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    // Calculate likelihood:
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logP = argLikelihood.calculateLogP();
    // Compare product of likelihoods of "marginal alignments" with
    // likelihood computed using RGL.
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logPprime = argLikelihoodSlow.calculateLogP();
    double relError = 2.0 * Math.abs(logP - logPprime) / Math.abs(logP + logPprime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP, logPprime, relError);
    assertTrue(relError < 1e-13);
}

Example 17 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGLikelihoodTest method testBeagleLikelihood.

@Test
public void testBeagleLikelihood() 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
    ACGLikelihoodBeagle argLikelihood = new ACGLikelihoodBeagle();
    argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    acg.setEverythingDirty(true);
    try {
        double logP = argLikelihood.calculateLogP();
        double logPtrue = argLikelihoodSlow.calculateLogP();
        double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + 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);
    } catch (RuntimeException ex) {
        System.err.println("Beagle library not found: skipping beagle likelihood test.");
    }
}

Example 18 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodFixedData.

@Test
public void testLikelihoodFixedData() 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);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    acg.setEverythingDirty(true);
    double logP = argLikelihood.calculateLogP();
    double logPtrue = argLikelihoodSlow.calculateLogP();
    double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + 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 19 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class AddRemoveConversionTest method testProbability.

/**
 * Tests whether probability of proposing a conversion lines up with
 * conversion probability found in ACGCoalescent.
 * @throws java.lang.Exception
 */
@Test
public void testProbability() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    SimulatedACG acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
    RealParameter rho = new RealParameter(Double.toString(1.0 / locus.getSiteCount()));
    RealParameter delta = new RealParameter("50.0");
    AddRemoveConversion operator = new AddRemoveConversion();
    operator.initByName("weight", 1.0, "acg", acg, "delta", delta, "populationModel", popFunc);
    ACGCoalescent coal = new ACGCoalescent();
    coal.initByName("tree", acg, "populationModel", popFunc, "rho", rho, "delta", delta);
    double logP1 = 0.0;
    double logP2 = 0.0;
    for (Conversion conv : acg.getConversions(locus)) {
        logP1 += operator.getConversionProb(conv);
        logP2 += coal.calculateConversionLogP(conv);
    }
    System.out.println("logP1 = " + logP1);
    System.out.println("logP2 = " + logP2);
    assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}

Example 20 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGCoalescent method calculateLogP.

@Override
public double calculateLogP() {
    // Check whether conversion count exceeds bounds.
    if (acg.getTotalConvCount() < lowerCCBoundInput.get() || acg.getTotalConvCount() > upperCCBoundInput.get())
        return Double.NEGATIVE_INFINITY;
    logP = calculateClonalFrameLogP();
    double poissonMean = rhoInput.get().getValue() * acg.getClonalFrameLength() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
    // Probability of conversion count:
    if (poissonMean > 0.0) {
        logP += -poissonMean + acg.getTotalConvCount() * Math.log(poissonMean);
    // - GammaFunction.lnGamma(acg.getConvCount()+1);
    } else {
        if (acg.getTotalConvCount() > 0)
            logP = Double.NEGATIVE_INFINITY;
    }
    for (Locus locus : acg.getLoci()) for (Conversion conv : acg.getConversions(locus)) logP += calculateConversionLogP(conv);
    if (lowerCCBoundInput.get() > 0 || upperCCBoundInput.get() < Integer.MAX_VALUE) {
        try {
            logP -= new PoissonDistributionImpl(poissonMean).cumulativeProbability(lowerCCBoundInput.get(), upperCCBoundInput.get());
        } catch (MathException e) {
            throw new RuntimeException("Error computing modification to ARG " + "prior density required by conversion number constraint.");
        }
    }
    return logP;
}