Examples with ConstantPopulation beast.evolution.tree.coalescent.ConstantPopulation used on opensource projects

Example 6 with ConstantPopulation

use of beast.evolution.tree.coalescent.ConstantPopulation in project bacter by tgvaughan.

the class SimulatedAlignmentTest method test.

@Test
public void test() throws Exception {
    // Randomizer.setSeed(26);
    Randomizer.setSeed(7);
    Locus locus = new Locus("locus", 100000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / 100000, "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.0"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    for (Region region : acg.getRegions(locus)) System.out.println(new MarginalTree(acg, region));
    // (Should be enough info here for precise agreement)
    for (Region region : acg.getRegions(locus)) {
        Alignment margAlign = createMarginalAlignment(alignment, region);
        ClusterTree upgmaTree = new ClusterTree();
        upgmaTree.initByName("clusterType", "upgma", "taxa", margAlign);
        MarginalTree marginalTree = new MarginalTree(acg, region);
        // outfMarg.println(marginalTree.getRoot() + ";");
        // outfMarg.flush();
        // 
        // outfUPGMA.println(upgmaTree.getRoot() + ";");
        // outfUPGMA.flush();
        assertTrue(topologiesEquivalent(marginalTree.getRoot(), upgmaTree.getRoot()));
    }
// outfMarg.close();
// outfUPGMA.close();
}

Example 7 with ConstantPopulation

use of beast.evolution.tree.coalescent.ConstantPopulation in project bacter by tgvaughan.

the class AddRemoveConversionTest method testHR.

/**
 * Tests that probability density of forward move calculated
 * by drawNewConversion() matches probability density of backward
 * move calculated by getConversionProb().
 *
 * @throws Exception
 */
@Test
public void testHR() 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);
    AddRemoveConversion operator = new AddRemoveConversion();
    // Loop until a valid proposal is made
    double logP1;
    List<Conversion> oldConversions;
    do {
        operator.initByName("weight", 1.0, "acg", acg, "delta", new RealParameter("50.0"), "populationModel", popFunc);
        oldConversions = Lists.newArrayList(acg.getConversions(locus));
        logP1 = operator.drawNewConversion();
    } while (Double.isInfinite(logP1));
    System.out.println("logP1 = " + logP1);
    // Identify new recomination
    Conversion newRecomb = null;
    for (Conversion recomb : acg.getConversions(locus)) {
        if (!oldConversions.contains(recomb))
            newRecomb = recomb;
    }
    assertNotNull(newRecomb);
    double logP2 = operator.getConversionProb(newRecomb);
    System.out.println("logP2 = " + logP2);
    assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}

Example 8 with ConstantPopulation

use of beast.evolution.tree.coalescent.ConstantPopulation in project bacter by tgvaughan.

the class AddRemoveConversion method main.

public static void main(String[] args) throws Exception {
    ConversionGraph acg = new ConversionGraph();
    ConstantPopulation popFunc = new ConstantPopulation();
    AddRemoveConversion operator = new AddRemoveConversion();
    operator.initByName("weight", 1.0, "acg", acg, "populationModel", popFunc, "rho", new RealParameter(Double.toString(1.0 / 10000.0)), "delta", new RealParameter("50.0"));
    popFunc.initByName("popSize", new RealParameter("1.0"));
    TaxonSet taxonSet = new TaxonSet();
    taxonSet.taxonsetInput.setValue(new Taxon("t1"), taxonSet);
    taxonSet.taxonsetInput.setValue(new Taxon("t2"), taxonSet);
    Locus locus = new Locus("locus", 10000);
    try (PrintStream ps = new PrintStream("out.txt")) {
        for (int i = 0; i < 100000; i++) {
            acg.initByName("locus", locus, "taxonset", taxonSet, "fromString", "(0:1.0,1:1.0)2:0.0;");
            operator.drawNewConversion();
            ps.println(acg.getConversions(locus).get(0).getStartSite() + " " + acg.getConversions(locus).get(0).getEndSite());
        }
    }
}

Example 9 with ConstantPopulation

use of beast.evolution.tree.coalescent.ConstantPopulation 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 10 with ConstantPopulation

use of beast.evolution.tree.coalescent.ConstantPopulation 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);
}