Examples with TaxonSet beast.evolution.alignment.TaxonSet used on opensource projects

Example 21 with TaxonSet

use of beast.evolution.alignment.TaxonSet 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 22 with TaxonSet

use of beast.evolution.alignment.TaxonSet 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 23 with TaxonSet

use of beast.evolution.alignment.TaxonSet 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 24 with TaxonSet

use of beast.evolution.alignment.TaxonSet 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 25 with TaxonSet

use of beast.evolution.alignment.TaxonSet in project beast2 by CompEvol.

the class Tree method getTaxaNames.

/**
 * @returns an array of taxon names in order of their node numbers.
 * Note that in general no special order of taxa is assumed, just the order
 * assumed in this tree. Consider using tree.m_taxonset.get().asStringList()
 * instead.
 */
public String[] getTaxaNames() {
    if (m_sTaxaNames == null || (m_sTaxaNames.length == 1 && m_sTaxaNames[0] == null) || m_sTaxaNames.length == 0) {
        // take taxa from tree if one exists
        if (root != null) {
            m_sTaxaNames = new String[getNodeCount()];
            collectTaxaNames(getRoot());
            List<String> taxaNames = new ArrayList<>();
            for (int i = 0; i < m_sTaxaNames.length && i < getLeafNodeCount(); i++) {
                String name = m_sTaxaNames[i];
                if (name != null) {
                    taxaNames.add(name);
                }
            }
            m_sTaxaNames = taxaNames.toArray(new String[] {});
        } else {
            // no tree? use taxon set.
            final TaxonSet taxonSet = m_taxonset.get();
            if (taxonSet != null) {
                final List<String> txs = taxonSet.asStringList();
                m_sTaxaNames = txs.toArray(new String[txs.size()]);
            } else {
                Log.err("No taxa specified");
            }
        }
    }
    // sanity check
    if (m_sTaxaNames.length == 1 && m_sTaxaNames[0] == null) {
        Log.warning("WARNING: tree interrogated for taxa, but the tree was not initialised properly. To fix this, specify the taxonset input");
    }
    return m_sTaxaNames;
}