use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.
the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.
public void testJointLikelihood() {
TreeModel treeModel = new TreeModel("treeModel", tree);
Sequence[] sequence = new Sequence[3];
sequence[0] = new Sequence(new Taxon("0"), "A");
sequence[1] = new Sequence(new Taxon("1"), "C");
sequence[2] = new Sequence(new Taxon("2"), "C");
Taxa taxa = new Taxa();
for (Sequence s : sequence) {
taxa.addTaxon(s.getTaxon());
}
SimpleAlignment alignment = new SimpleAlignment();
for (Sequence s : sequence) {
alignment.addSequence(s);
}
Parameter mu = new Parameter.Default(1, 1.0);
Parameter kappa = new Parameter.Default(1, 1.0);
double[] pi = { 0.25, 0.25, 0.25, 0.25 };
Parameter freqs = new Parameter.Default(pi);
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
HKY hky = new HKY(kappa, f);
GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", mu, null, -1, null);
siteRateModel.setSubstitutionModel(hky);
BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel());
BranchRateModel branchRateModel = null;
AncestralStateBeagleTreeLikelihood treeLikelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, hky.getDataType(), "stateTag", // useMap = true
true, false);
double logLike = treeLikelihood.getLogLikelihood();
StringBuffer buffer = new StringBuffer();
// Tree.Utils.newick(treeModel, treeModel.getRoot(), false, Tree.BranchLengthType.LENGTHS_AS_TIME,
// null, null, new NodeAttributeProvider[]{treeLikelihood}, null, null, buffer);
TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
System.out.println(buffer);
System.out.println("t_CA(2) = " + t(false, 2.0));
System.out.println("t_CC(1) = " + t(true, 1.0));
double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
assertEquals(logLike, Math.log(trueValue), 1e-6);
}
use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.
the class PhylipSequentialImporter method importAlignment.
/**
* importAlignment.
*/
public Alignment importAlignment() throws IOException, Importer.ImportException {
SimpleAlignment alignment = null;
try {
int taxonCount = readInteger();
int siteCount = readInteger();
String firstSeq = null;
for (int i = 0; i < taxonCount; i++) {
StringBuffer name = new StringBuffer();
char ch = read();
int n = 0;
while (!Character.isWhitespace(ch) && (maxNameLength < 1 || n < maxNameLength)) {
name.append(ch);
ch = read();
n++;
}
StringBuffer seq = new StringBuffer(siteCount);
readSequence(seq, dataType, "", siteCount, "-", "?", ".", firstSeq);
if (firstSeq == null) {
firstSeq = seq.toString();
}
if (alignment == null) {
alignment = new SimpleAlignment();
}
alignment.addSequence(new Sequence(new Taxon(name.toString()), seq.toString()));
}
} catch (EOFException e) {
}
return alignment;
}
use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.
the class DataLikelihoodTester2 method createAlignment.
private static SimpleAlignment createAlignment(Object[][] taxa_sequence, DataType dataType) {
SimpleAlignment alignment = new SimpleAlignment();
alignment.setDataType(dataType);
// alignment.setDataType(Nucleotides.INSTANCE);
// 6, 17
Taxon[] taxa = new Taxon[taxa_sequence[0].length];
System.out.println("Taxon len = " + taxa_sequence[0].length);
System.out.println("Alignment len = " + taxa_sequence[1].length);
if (taxa_sequence.length > 2)
System.out.println("Date len = " + taxa_sequence[2].length);
for (int i = 0; i < taxa_sequence[0].length; i++) {
taxa[i] = new Taxon(taxa_sequence[0][i].toString());
if (taxa_sequence.length > 2) {
Date date = new Date((Double) taxa_sequence[2][i], Units.Type.YEARS, (Boolean) taxa_sequence[3][0]);
taxa[i].setDate(date);
}
//taxonList.addTaxon(taxon);
Sequence sequence = new Sequence(taxa_sequence[1][i].toString());
sequence.setTaxon(taxa[i]);
sequence.setDataType(dataType);
alignment.addSequence(sequence);
}
return alignment;
}
use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.
the class SequenceLikelihoodTest method createAlignmentWithAllUniquePatterns.
protected void createAlignmentWithAllUniquePatterns(Object[][] taxa_sequence, DataType dataType) {
alignment = new SimpleAlignment();
alignment.setDataType(dataType);
int nTaxa = taxa_sequence[0].length;
String[] allUniquePatterns = createAllUniquePatterns(nTaxa, dataType);
taxa_sequence[1] = allUniquePatterns;
// 6, 17
taxa = new Taxon[nTaxa];
System.out.println("Taxon len = " + taxa_sequence[0].length);
System.out.println("Alignment len = " + taxa_sequence[1].length);
if (taxa_sequence.length > 2)
System.out.println("Date len = " + taxa_sequence[2].length);
for (int i = 0; i < taxa_sequence[0].length; i++) {
taxa[i] = new Taxon(taxa_sequence[0][i].toString());
if (taxa_sequence.length > 2) {
Date date = new Date((Double) taxa_sequence[2][i], Units.Type.YEARS, (Boolean) taxa_sequence[3][0]);
taxa[i].setDate(date);
}
//taxonList.addTaxon(taxon);
Sequence sequence = new Sequence(taxa_sequence[1][i].toString());
sequence.setTaxon(taxa[i]);
sequence.setDataType(dataType);
alignment.addSequence(sequence);
}
System.out.println("Sequence pattern count = " + alignment.getPatternCount());
}
use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.
the class AncestralStateTreeLikelihoodTest method testJointLikelihood.
public void testJointLikelihood() {
TreeModel treeModel = new TreeModel("treeModel", tree);
Sequence[] sequence = new Sequence[3];
sequence[0] = new Sequence(new Taxon("0"), "A");
sequence[1] = new Sequence(new Taxon("1"), "C");
sequence[2] = new Sequence(new Taxon("2"), "C");
Taxa taxa = new Taxa();
for (Sequence s : sequence) {
taxa.addTaxon(s.getTaxon());
}
SimpleAlignment alignment = new SimpleAlignment();
for (Sequence s : sequence) {
alignment.addSequence(s);
}
Parameter mu = new Parameter.Default(1, 1.0);
Parameter kappa = new Parameter.Default(1, 1.0);
double[] pi = { 0.25, 0.25, 0.25, 0.25 };
Parameter freqs = new Parameter.Default(pi);
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
HKY hky = new HKY(kappa, f);
AncestralStateTreeLikelihood treeLikelihood = new AncestralStateTreeLikelihood(alignment, treeModel, new GammaSiteModel(hky), new StrictClockBranchRates(mu), false, true, Nucleotides.INSTANCE, "state", false, // useMap = true
true, false);
double logLike = treeLikelihood.getLogLikelihood();
StringBuffer buffer = new StringBuffer();
TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
System.out.println(buffer);
System.out.println("t_CA(2) = " + t(false, 2.0));
System.out.println("t_CC(1) = " + t(true, 1.0));
double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
assertEquals(logLike, Math.log(trueValue), 1e-6);
}
Aggregations