use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.
the class NarrowExchangeTest method testNarrow.
/**
* Test method for {@link dr.evomodel.operators.ExchangeOperator#narrow()}.
*/
public void testNarrow() throws IOException, ImportException {
// probability of picking B node is 1/(2n-4) = 1/6
// probability of swapping it with C is 1/1
// total = 1/6
System.out.println("Test 1: Forward");
String treeMatch = "((((A,C),B),D),E);";
int count = 0;
int reps = 100000;
for (int i = 0; i < reps; i++) {
TreeModel treeModel = new TreeModel("treeModel", tree5);
ExchangeOperator operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1);
operator.doOperation();
String tree = TreeUtils.newickNoLengths(treeModel);
if (tree.equals(treeMatch)) {
count += 1;
}
}
double p_1 = (double) count / (double) reps;
System.out.println("Number of proposals:\t" + count);
System.out.println("Number of tries:\t" + reps);
System.out.println("Number of ratio:\t" + p_1);
System.out.println("Number of expected ratio:\t" + 1.0 / 6.0);
assertExpectation(1.0 / 6.0, p_1, reps);
}
use of dr.evomodel.tree.TreeModel 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.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.
the class ImportancePruneAndRegraftTestProblem method testDoOperation.
/**
* Test method for {@link SimpleMCMCOperator#doOperation()}.
* @throws ImportException
* @throws IOException
*/
public void testDoOperation() throws IOException, ImportException {
// probability of picking (A,B) node is 1/(2n-3) = 1/7
// probability of swapping with D is 1/2
// total = 1/14
//probability of picking {D} node is 1/(2n-3) = 1/7
//probability of picking {A,B} is 1/5
// total = 1/35
//total = 1/14 + 1/35 = 7/70 = 0.1
System.out.println("Test 1: Forward");
String treeMatch = "(((D,C),(A,B)),E);";
int count = 0;
int reps = 1000000;
for (int i = 0; i < reps; i++) {
TreeModel treeModel = new TreeModel("treeModel", tree5);
ImportancePruneAndRegraft operator = new ImportancePruneAndRegraft(treeModel, 1.0, 0);
operator.doOperation();
String tree = TreeUtils.newickNoLengths(treeModel);
if (tree.equals(treeMatch)) {
count += 1;
}
}
double p_1 = (double) count / (double) reps;
System.out.println("Number of proposals:\t" + count);
System.out.println("Number of tries:\t" + reps);
System.out.println("Number of ratio:\t" + p_1);
System.out.println("Number of expected ratio:\t" + 0.1);
assertExpectation(0.1, p_1, reps);
// lets see what the backward probability is for the hastings ratio
// (((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0) -> ((((A,B),C),D),E)
// probability of picking (A,B) node is 1/(2n-3) = 1/7
// probability of swapping with D is 1/3
// total = 1/21
//probability of picking {D} node is 1/(2n-2) = 1/7
//probability of picking {A,B} is 1/4
// total = 1/28
//total = 1/21 + 1/28 = 7/84 = 0.08333333
System.out.println("Test 2: Backward");
treeMatch = "((((A,B),C),D),E);";
NewickImporter importer = new NewickImporter("(((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0);");
FlexibleTree tree5_2 = (FlexibleTree) importer.importTree(null);
count = 0;
for (int i = 0; i < reps; i++) {
TreeModel treeModel = new TreeModel("treeModel", tree5_2);
ImportancePruneAndRegraft operator = new ImportancePruneAndRegraft(treeModel, 1.0, 1);
operator.doOperation();
String tree = TreeUtils.newickNoLengths(treeModel);
if (tree.equals(treeMatch)) {
count += 1;
}
}
double p_2 = (double) count / (double) reps;
System.out.println("Number of proposals:\t" + count);
System.out.println("Number of tries:\t" + reps);
System.out.println("Number of ratio:\t" + p_2);
System.out.println("Number of expected ratio:\t" + 0.0833333);
assertExpectation(0.0833333, p_2, reps);
}
use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.
the class ImportanceSubtreeSwapTestProblem method testDoOperation.
/**
* Test method for {@link SimpleMCMCOperator#doOperation()}.
* @throws ImportException
* @throws IOException
*/
public void testDoOperation() throws IOException, ImportException {
// probability of picking (A,B) node is 1/(2n-3) = 1/7
// probability of swapping with D is 1/2
// total = 1/14
//probability of picking {D} node is 1/(2n-3) = 1/7
//probability of picking {A,B} is 1/5
// total = 1/35
//total = 1/14 + 1/35 = 7/70 = 0.1
System.out.println("Test 1: Forward");
String treeMatch = "(((D,C),(A,B)),E);";
int count = 0;
int reps = 100000;
for (int i = 0; i < reps; i++) {
TreeModel treeModel = new TreeModel("treeModel", tree5);
ImportanceSubtreeSwap operator = new ImportanceSubtreeSwap(treeModel, 1.0, 0);
operator.doOperation();
String tree = TreeUtils.newickNoLengths(treeModel);
if (tree.equals(treeMatch)) {
count += 1;
}
}
double p_1 = (double) count / (double) reps;
System.out.println("Number of proposals:\t" + count);
System.out.println("Number of tries:\t" + reps);
System.out.println("Number of ratio:\t" + p_1);
System.out.println("Number of expected ratio:\t" + 0.1);
assertExpectation(0.1, p_1, reps);
// lets see what the backward probability is for the hastings ratio
// (((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0) -> ((((A,B),C),D),E)
// probability of picking (A,B) node is 1/(2n-3) = 1/7
// probability of swapping with D is 1/3
// total = 1/21
//probability of picking {D} node is 1/(2n-2) = 1/7
//probability of picking {A,B} is 1/4
// total = 1/28
//total = 1/21 + 1/28 = 7/84 = 0.08333333
System.out.println("Test 2: Backward");
treeMatch = "((((A,B),C),D),E);";
NewickImporter importer = new NewickImporter("(((D:2.0,C:2.0):1.0,(A:1.0,B:1.0):2.0):1.0,E:4.0);");
FlexibleTree tree5_2 = (FlexibleTree) importer.importTree(null);
count = 0;
for (int i = 0; i < reps; i++) {
TreeModel treeModel = new TreeModel("treeModel", tree5_2);
ImportanceSubtreeSwap operator = new ImportanceSubtreeSwap(treeModel, 1.0, 0);
operator.doOperation();
String tree = TreeUtils.newickNoLengths(treeModel);
if (tree.equals(treeMatch)) {
count += 1;
}
}
double p_2 = (double) count / (double) reps;
System.out.println("Number of proposals:\t" + count);
System.out.println("Number of tries:\t" + reps);
System.out.println("Number of ratio:\t" + p_2);
System.out.println("Number of expected ratio:\t" + 0.0833333);
assertExpectation(0.0833333, p_2, reps);
}
use of dr.evomodel.tree.TreeModel 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