Examples with Tree dr.evolution.tree.Tree used on opensource projects

Example 6 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class WeightedMultiplicativeBinary method analyzeTrace.

/**
     * Actually analyzes the trace given the burn-in. Each tree from the trace
     * is read and the conditional clade frequencies incremented.
     *
     * @param verbose if true then progress is logged to stdout
     */
public void analyzeTrace(boolean verbose) {
    if (verbose) {
        if (traces.length > 1)
            System.out.println("Combining " + traces.length + " traces.");
    }
    // get first tree to extract the taxon
    Tree tree = getTree(0);
    // read every tree from the trace
    for (TreeTrace trace : traces) {
        // do some output stuff
        int treeCount = trace.getTreeCount(burnin * trace.getStepSize());
        double stepSize = treeCount / 60.0;
        int counter = 1;
        if (verbose) {
            System.out.println("Analyzing " + treeCount + " trees...");
            System.out.println("0              25             50             75            100");
            System.out.println("|--------------|--------------|--------------|--------------|");
            System.out.print("*");
        }
        for (int i = 1; i < treeCount; i++) {
            // get the next tree
            tree = trace.getTree(i, burnin * trace.getStepSize());
            // add the tree and its clades to the frequencies
            addTree(tree);
            // some more output stuff
            if (i >= (int) Math.round(counter * stepSize) && counter <= 60) {
                if (verbose) {
                    System.out.print("*");
                    System.out.flush();
                }
                counter += 1;
            }
        }
        if (verbose) {
            System.out.println("*");
        }
    }
}

Example 7 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class OperatorAssert method irreducibilityTester.

private void irreducibilityTester(Tree tree, int numLabelledTopologies, int chainLength, int sampleTreeEvery) throws IOException, Importer.ImportException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    TreeModel treeModel = new TreeModel("treeModel", tree);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    OperatorSchedule schedule = getOperatorSchedule(treeModel);
    Parameter b = new Parameter.Default("b", 2.0, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    MCLogger[] loggers = new MCLogger[2];
    //        loggers[0] = new MCLogger(new ArrayLogFormatter(false), 100, false);
    //        loggers[0].add(likelihood);
    //        loggers[0].add(rootHeight);
    //        loggers[0].add(tls);
    loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[0].add(likelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tls);
    File file = new File("yule.trees");
    file.deleteOnExit();
    FileOutputStream out = new FileOutputStream(file);
    loggers[1] = new TreeLogger(treeModel, new TabDelimitedFormatter(out), sampleTreeEvery, true, true, false);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, likelihood, schedule, loggers);
    mcmc.run();
    out.flush();
    out.close();
    Set<String> uniqueTrees = new HashSet<String>();
    HashMap<String, Integer> topologies = new HashMap<String, Integer>();
    HashMap<String, HashMap<String, Integer>> treeCounts = new HashMap<String, HashMap<String, Integer>>();
    NexusImporter importer = new NexusImporter(new FileReader(file));
    int sampleSize = 0;
    while (importer.hasTree()) {
        sampleSize++;
        Tree t = importer.importNextTree();
        String uniqueNewick = TreeUtils.uniqueNewick(t, t.getRoot());
        String topology = uniqueNewick.replaceAll("\\w+", "X");
        if (!uniqueTrees.contains(uniqueNewick)) {
            uniqueTrees.add(uniqueNewick);
        }
        HashMap<String, Integer> counts;
        if (topologies.containsKey(topology)) {
            topologies.put(topology, topologies.get(topology) + 1);
            counts = treeCounts.get(topology);
        } else {
            topologies.put(topology, 1);
            counts = new HashMap<String, Integer>();
            treeCounts.put(topology, counts);
        }
        if (counts.containsKey(uniqueNewick)) {
            counts.put(uniqueNewick, counts.get(uniqueNewick) + 1);
        } else {
            counts.put(uniqueNewick, 1);
        }
    }
    TestCase.assertEquals(numLabelledTopologies, uniqueTrees.size());
    TestCase.assertEquals(sampleSize, chainLength / sampleTreeEvery + 1);
    Set<String> keys = topologies.keySet();
    double ep = 1.0 / topologies.size();
    for (String topology : keys) {
        double ap = ((double) topologies.get(topology)) / (sampleSize);
        //          	assertExpectation(ep, ap, sampleSize);
        HashMap<String, Integer> counts = treeCounts.get(topology);
        Set<String> trees = counts.keySet();
        double MSE = 0;
        double ep1 = 1.0 / counts.size();
        for (String t : trees) {
            double ap1 = ((double) counts.get(t)) / (topologies.get(topology));
            //              	assertExpectation(ep1, ap1, topologies.get(topology));
            MSE += (ep1 - ap1) * (ep1 - ap1);
        }
        MSE /= counts.size();
        System.out.println("The Mean Square Error for the topolgy " + topology + " is " + MSE);
    }
}

Example 8 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class TestCalibratedYuleModel method calibration.

private Tree calibration(final TaxonList taxa, DemographicModel demoModel) throws Exception {
    dr.evolution.coalescent.CoalescentSimulator simulator = new dr.evolution.coalescent.CoalescentSimulator();
    DemographicFunction demoFunction = demoModel.getDemographicFunction();
    SimpleNode[] firstHalfNodes = new SimpleNode[taxa.getTaxonCount() / 2];
    SimpleNode[] secondHalfNodes = new SimpleNode[taxa.getTaxonCount() - taxa.getTaxonCount() / 2];
    for (int i = 0; i < firstHalfNodes.length; i++) {
        firstHalfNodes[i] = new SimpleNode();
        firstHalfNodes[i].setTaxon(taxa.getTaxon(i));
    }
    for (int i = 0; i < secondHalfNodes.length; i++) {
        secondHalfNodes[i] = new SimpleNode();
        secondHalfNodes[i].setTaxon(taxa.getTaxon(i + taxa.getTaxonCount() / 2));
    }
    SimpleNode firstHalfRootNode = simulator.simulateCoalescent(firstHalfNodes, demoFunction);
    SimpleNode[] restNodes = simulator.simulateCoalescent(secondHalfNodes, demoFunction, 0, firstHalfRootNode.getHeight());
    SimpleNode[] together = new SimpleNode[restNodes.length + 1];
    for (int i = 0; i < restNodes.length; i++) {
        together[i + 1] = restNodes[i];
    }
    together[0] = firstHalfRootNode;
    SimpleNode root = simulator.simulateCoalescent(together, demoFunction);
    Tree tree = new SimpleTree(root);
    return tree;
}

Example 9 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class TestCalibratedYuleModel method createTreeModel.

protected TreeModel createTreeModel(int treeSize) throws Exception {
    taxa = new Taxa();
    for (int i = 0; i < treeSize; i++) {
        taxa.addTaxon(new Taxon("T" + Integer.toString(i)));
    }
    //System.out.println("taxaSubSet_size = " + taxaSubSet.getTaxonCount());
    Parameter popSize = new Parameter.Default(treeSize);
    popSize.setId(ConstantPopulationModelParser.POPULATION_SIZE);
    ConstantPopulationModel startingTree = new ConstantPopulationModel(popSize, Units.Type.YEARS);
    Tree tree = calibration(taxa, startingTree);
    //treeModel
    return new TreeModel(tree);
}

Example 10 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateTwoPartitions.

// END: simulateOnePartition
static void simulateTwoPartitions() {
    try {
        System.out.println("Test case 3: simulateTwoPartitions");
        MathUtils.setSeed(666);
        int sequenceLength = 11;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        3, // every
        1);
        // create partition
        Partition Partition = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        4, // to
        sequenceLength - 1, // every
        1);
        Sequence ancestralSequence = new Sequence();
        ancestralSequence.appendSequenceString("TCAAGTG");
        Partition.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        partitionsList.add(Partition);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}