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("*");
}
}
}
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);
}
}
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;
}
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);
}
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
}
Aggregations