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

Example 1 with TreeIntervals

use of beast.evolution.tree.coalescent.TreeIntervals in project beast2 by CompEvol.

the class WilsonBaldingTest method topologyDistribution.

/**
 * Test topology distribution.
 * @throws Exception
 */
@Test
public void topologyDistribution() throws Exception {
    // Fix seed: will hopefully ensure success of test unless something
    // goes terribly wrong.
    Randomizer.setSeed(42);
    // Assemble model:
    ConstantPopulation constantPop = new ConstantPopulation();
    constantPop.initByName("popSize", new RealParameter("10000.0"));
    List<Object> alignmentInitArgs = new ArrayList<Object>();
    for (int i = 0; i < 4; i++) {
        Sequence thisSeq = new Sequence();
        thisSeq.initByName("taxon", String.valueOf(i), "value", "?");
        alignmentInitArgs.add("sequence");
        alignmentInitArgs.add(thisSeq);
    }
    Alignment alignment = new Alignment();
    alignment.initByName(alignmentInitArgs.toArray());
    Tree tree = new RandomTree();
    tree.initByName("taxa", alignment, "populationModel", constantPop);
    TreeIntervals treeIntervals = new TreeIntervals();
    treeIntervals.initByName("tree", tree);
    Coalescent coalescentDistrib = new Coalescent();
    coalescentDistrib.initByName("treeIntervals", treeIntervals, "populationModel", constantPop);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", tree);
    // Set up operator:
    WilsonBalding wilsonBalding = new WilsonBalding();
    wilsonBalding.initByName("weight", "1", "tree", tree);
    // Set up logger:
    TreeReport treeReport = new TreeReport();
    treeReport.initByName("logEvery", "100", "burnin", "200000", "credibleSetPercentage", "95.0", "log", tree, "silent", true);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "2000000", "state", state, "distribution", coalescentDistrib, "operator", wilsonBalding, "logger", treeReport);
    // Run MCMC:
    mcmc.run();
    // Obtain analysis results:
    TreeTraceAnalysis analysis = treeReport.getAnalysis();
    Map<String, Integer> topologyCounts = analysis.getTopologyCounts();
    int totalTreesUsed = analysis.getNTrees();
    // Test topology distribution against ideal:
    double tol = 0.005;
    for (int i = 0; i < topologies.length; i++) {
        double thisProb = topologyCounts.get(topologies[i]) / (double) totalTreesUsed;
        boolean withinTol = (thisProb > probs[i] - tol && thisProb < probs[i] + tol);
        Assert.assertTrue(withinTol);
        System.err.format("Topology %s rel. freq. %.3f", topologies[i], thisProb);
        if (withinTol)
            System.err.println(" (Within tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
        else
            System.err.println(" (FAILURE: outside tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
    }
}

Example 2 with TreeIntervals

use of beast.evolution.tree.coalescent.TreeIntervals in project beast2 by CompEvol.

the class RandomTreeTest method testCoalescentTimes.

@Test
public void testCoalescentTimes() throws Exception {
    Randomizer.setSeed(53);
    int Nleaves = 10;
    int Niter = 5000;
    // (Serially sampled) coalescent time means and variances
    // estimated from 50000 trees simulated using MASTER
    double[] coalTimeMeansTruth = { 1.754662, 2.833337, 3.843532, 4.850805, 5.849542, 6.847016, 7.8482, 8.855137, 10.15442 };
    double[] coalTimeVarsTruth = { 0.2751625, 0.2727121, 0.2685172, 0.2705117, 0.2678611, 0.2671793, 0.2686952, 0.2828477, 1.076874 };
    // Assemble BEASTObjects needed by RandomTree
    StringBuilder traitSB = new StringBuilder();
    List<Sequence> seqList = new ArrayList<Sequence>();
    for (int i = 0; i < Nleaves; i++) {
        String taxonID = "t " + i;
        seqList.add(new Sequence(taxonID, "?"));
        if (i > 0)
            traitSB.append(",");
        traitSB.append(taxonID).append("=").append(i);
    }
    Alignment alignment = new Alignment(seqList, "nucleotide");
    TaxonSet taxonSet = new TaxonSet(alignment);
    TraitSet timeTrait = new TraitSet();
    timeTrait.initByName("traitname", "date-backward", "taxa", taxonSet, "value", traitSB.toString());
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    // Create RandomTree and TreeInterval instances
    RandomTree tree = new RandomTree();
    TreeIntervals intervals = new TreeIntervals();
    // Estimate coalescence time moments
    double[] coalTimeMeans = new double[Nleaves - 1];
    double[] coalTimeVars = new double[Nleaves - 1];
    double[] coalTimes = new double[Nleaves - 1];
    for (int i = 0; i < Niter; i++) {
        tree.initByName("taxa", alignment, "populationModel", popFunc, "trait", timeTrait);
        intervals.initByName("tree", tree);
        intervals.getCoalescentTimes(coalTimes);
        for (int j = 0; j < Nleaves - 1; j++) {
            coalTimeMeans[j] += coalTimes[j];
            coalTimeVars[j] += coalTimes[j] * coalTimes[j];
        }
    }
    // Normalise means and variances
    for (int j = 0; j < Nleaves - 1; j++) {
        coalTimeMeans[j] /= Niter;
        coalTimeVars[j] /= Niter;
        coalTimeVars[j] -= coalTimeMeans[j] * coalTimeMeans[j];
    }
    // Test means and variances against independently estimated values
    for (int j = 0; j < Nleaves - 1; j++) {
        assert (relError(coalTimeMeans[j], coalTimeMeansTruth[j]) < 5e-3);
        assert (relError(coalTimeVars[j], coalTimeVarsTruth[j]) < 5e-2);
    }
}

Example 3 with TreeIntervals

use of beast.evolution.tree.coalescent.TreeIntervals in project beast2 by CompEvol.

the class CoalescentTest method testConstantPopulation.

public void testConstantPopulation() throws Exception {
    // *********** 3 taxon **********
    Tree tree = getTree(data, trees[0]);
    TreeIntervals treeIntervals = new TreeIntervals();
    treeIntervals.initByName("tree", tree);
    ConstantPopulation cp = new ConstantPopulation();
    cp.initByName("popSize", Double.toString(pop));
    Coalescent coal = new Coalescent();
    coal.initByName("treeIntervals", treeIntervals, "populationModel", cp);
    double logL = coal.calculateLogP();
    assertEquals(logL, -(4 / pop) - 2 * Math.log(pop), PRECISION);
// *********** 4 taxon **********
// tree = getTree(data, trees[1]);
// treeIntervals = new TreeIntervals();
// treeIntervals.initByName("tree", tree);
// 
// cp = new ConstantPopulation();
// cp.initByName("popSize", Double.toString(pop));
// 
// coal = new Coalescent();
// coal.initByName("treeIntervals", treeIntervals, "populationModel", cp);
// 
// logL = coal.calculateLogP();
// 
// assertEquals(logL, -(4 / pop) - 2 * Math.log(pop), PRECISION);
}

Example 4 with TreeIntervals

use of beast.evolution.tree.coalescent.TreeIntervals in project beast2 by CompEvol.

the class BayesianSkylineTest method testSkyline.

@Test
public void testSkyline() throws Exception {
    // RealParameter popSize = new RealParameter("1.0", 0.0, 10.0, 2);
    // IntegerParameter groupSize = new IntegerParameter("2", 1, 4, 2);
    // popSize.setValue(1, 2.0);
    Tree tree = new Tree("(((1:1,2:1):2.5,(3:1.5,4:1.5):2):2,5:5.5);");
    TreeIntervals intervals = new TreeIntervals(tree);
    BayesianSkyline skyline = new BayesianSkyline();
    // skyline.init(popSize, groupSize, intervals);
    skyline.initByName("popSizes", "1.0 2.0", "groupSizes", "2 2", "treeIntervals", intervals);
    assertEquals(skyline.getPopSize(0.01), 1.0);
    assertEquals(skyline.getPopSize(1.49), 1.0);
    assertEquals(skyline.getPopSize(1.51), 2.0);
    assertEquals(skyline.getPopSize(5.51), 2.0);
}