use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.
the class BirthDeathSSLikelihoodTest method setUp.
public void setUp() throws Exception {
super.setUp();
NewickImporter importer = new NewickImporter("((1:1.0,2:1.0):1.0,3:2.0);");
tree = (FlexibleTree) importer.importTree(null);
importer = new NewickImporter("((3 : 1.5, 4 : 0.5) : 1 , (1 : 2, 2 : 1) : 3);");
tree2 = (FlexibleTree) importer.importTree(null);
}
use of dr.evolution.io.NewickImporter 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.evolution.io.NewickImporter 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.evolution.io.NewickImporter 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
}
use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.
the class BeagleSeqSimTest method simulateAminoAcid.
// END: simulateThreePartitions
static void simulateAminoAcid() {
try {
System.out.println("Test case 4: simulateAminoAcid");
MathUtils.setSeed(666);
int sequenceLength = 10;
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 site model
GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
// create branch rate model
BranchRateModel branchRateModel = new DefaultBranchRateModel();
// create Frequency Model
Parameter freqs = new Parameter.Default(new double[] { 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05 });
FrequencyModel freqModel = new FrequencyModel(AminoAcids.INSTANCE, freqs);
// create substitution model
EmpiricalRateMatrix rateMatrix = Blosum62.INSTANCE;
EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, freqModel);
HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
// create partition
Partition partition1 = new //
Partition(//
treeModel, //
substitutionModel, //
siteRateModel, //
branchRateModel, //
freqModel, // from
0, // to
sequenceLength - 1, // every
1);
partitionsList.add(partition1);
// 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
}
Aggregations