use of dr.evolution.alignment.SimpleSiteList in project beast-mcmc by beast-dev.
the class BASTADensityTester method main.
public static void main(String[] args) {
MathUtils.setSeed(666);
// turn off logging to avoid screen noise...
Logger logger = Logger.getLogger("dr");
logger.setUseParentHandlers(false);
System.out.println("EXAMPLE 1: 4 taxa with 2 demes");
SimpleAlignment alignment = createAlignment(sequences, Nucleotides.INSTANCE);
TreeModel treeModel;
try {
treeModel = createSpecifiedTree("(((AB192965Japan2004:0.42343888910376215,AF189018Indonesia2005:1.4234388891037622)5:9.725099517053918,AF071228Spain1997:3.14853840615768)6:1.7275747971782511,AF105975Portugal1995:2.8761132033359313)");
} catch (Exception e) {
throw new RuntimeException("Unable to parse Newick tree");
}
List<String> states = new ArrayList<String>();
states.add("Asia");
states.add("West_Medit");
GeneralDataType dataType = new GeneralDataType(states);
Taxa taxonList = new Taxa();
Taxon taxonOne = new Taxon(sequences[0][0]);
taxonOne.setAttribute(DEME, states.get(0));
taxonList.addTaxon(taxonOne);
Taxon taxonTwo = new Taxon(sequences[0][1]);
taxonTwo.setAttribute(DEME, states.get(1));
taxonList.addTaxon(taxonTwo);
Taxon taxonThree = new Taxon(sequences[0][2]);
taxonThree.setAttribute(DEME, states.get(1));
taxonList.addTaxon(taxonThree);
Taxon taxonFour = new Taxon(sequences[0][3]);
taxonFour.setAttribute(DEME, states.get(0));
taxonList.addTaxon(taxonFour);
int[] pattern = new int[sequences[0].length];
pattern[0] = 0;
pattern[1] = 1;
pattern[2] = 1;
pattern[3] = 0;
SimpleSiteList patterns = new SimpleSiteList(dataType, taxonList);
patterns.addPattern(pattern);
double[] freqValues = new double[2];
freqValues[0] = 0.5;
freqValues[1] = 0.5;
FrequencyModel freqModel = new FrequencyModel(dataType, freqValues);
Parameter ratesParameter = new Parameter.Default("migration rates", 2);
ratesParameter.setParameterValue(0, 1.0);
ratesParameter.setParameterValue(1, 1.0);
Parameter popSizesParameter = new Parameter.Default("population sizes", 2);
popSizesParameter.setParameterValue(0, 1.0);
popSizesParameter.setParameterValue(1, 1.0);
// GeneralSubstitutionModel migrationModel = new GeneralSubstitutionModel("migrationModel", dataType, null, ratesParameter, 0);
SVSComplexSubstitutionModel migrationModel = new SVSComplexSubstitutionModel("migrationModel", dataType, freqModel, ratesParameter, null);
TaxonList includeSubtree = null;
List<TaxonList> excludeSubtrees = new ArrayList<TaxonList>();
int subIntervals = 2;
BranchRateModel branchRateModel = new DefaultBranchRateModel();
StructuredCoalescentLikelihood structured = null;
try {
structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
} catch (TreeUtils.MissingTaxonException missing) {
System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
}
System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
System.out.println("EXAMPLE 2: 5 taxa (2 identical sampling dates) with 2 demes");
try {
treeModel = createSpecifiedTree("((((AB192965Japan2004:0.48972300366304955,AJ842306France2004:0.48972300366304955):1.2052625264106087,AF189018Indonesia2005:2.6949855300736583):6.86618775860274,AF071228Spain1997:1.5611732886763985):1.4313107562768952,AF105975Portugal1995:0.9924840449532937)");
} catch (Exception e) {
throw new RuntimeException("Unable to parse Newick tree");
}
Taxon taxonFive = new Taxon(sequencesTwo[0][4]);
taxonFive.setAttribute(DEME, states.get(1));
taxonList.addTaxon(taxonFive);
pattern = new int[sequencesTwo[0].length];
pattern[0] = 0;
pattern[1] = 1;
pattern[2] = 1;
pattern[3] = 0;
pattern[4] = 1;
patterns = new SimpleSiteList(dataType, taxonList);
patterns.addPattern(pattern);
try {
structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
} catch (TreeUtils.MissingTaxonException missing) {
System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
}
System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
System.out.println("EXAMPLE 3: 6 taxa (2 identical sampling dates) with 3 demes");
try {
treeModel = createSpecifiedTree("(((AB192965Japan2004:4.509764060335954,((AF189018Indonesia2005:1.496797606742874,AJ842306France2004:0.4967976067428741):0.3617175987505199,FM212660Cameroon2007:3.858515205493394):3.6512488548425597):4.957323993357177,AF071228Spain1997:2.4670880536931303):0.9760399779205056,AF105975Portugal1995:1.443128031613636)");
} catch (Exception e) {
throw new RuntimeException("Unable to parse Newick tree");
}
states = new ArrayList<String>();
states.add("Asia");
states.add("West_Medit");
states.add("Africa");
dataType = new GeneralDataType(states);
taxonOne.setAttribute(DEME, states.get(0));
taxonTwo.setAttribute(DEME, states.get(1));
taxonThree.setAttribute(DEME, states.get(1));
taxonFour.setAttribute(DEME, states.get(0));
taxonFive.setAttribute(DEME, states.get(1));
Taxon taxonSix = new Taxon(sequencesThree[0][5]);
taxonSix.setAttribute(DEME, states.get(2));
taxonList.addTaxon(taxonSix);
pattern = new int[sequencesThree[0].length];
pattern[0] = 0;
pattern[1] = 1;
pattern[2] = 1;
pattern[3] = 0;
pattern[4] = 1;
pattern[5] = 2;
patterns = new SimpleSiteList(dataType, taxonList);
patterns.addPattern(pattern);
freqValues = new double[3];
freqValues[0] = 1.0 / 3.0;
freqValues[1] = 1.0 / 3.0;
freqValues[2] = 1 - freqValues[0] - freqValues[1];
freqModel = new FrequencyModel(dataType, freqValues);
ratesParameter = new Parameter.Default("migration rates", 6);
ratesParameter.setParameterValue(0, 1.0);
ratesParameter.setParameterValue(1, 1.0);
ratesParameter.setParameterValue(2, 1.0);
ratesParameter.setParameterValue(3, 1.0);
ratesParameter.setParameterValue(4, 1.0);
ratesParameter.setParameterValue(5, 1.0);
popSizesParameter = new Parameter.Default("population sizes", 3);
popSizesParameter.setParameterValue(0, 1.0);
popSizesParameter.setParameterValue(1, 1.0);
popSizesParameter.setParameterValue(2, 1.0);
// GeneralSubstitutionModel migrationModel = new GeneralSubstitutionModel("migrationModel", dataType, null, ratesParameter, 0);
migrationModel = new SVSComplexSubstitutionModel("migrationModel", dataType, freqModel, ratesParameter, null);
try {
structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
} catch (TreeUtils.MissingTaxonException missing) {
System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
}
System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
System.out.println("EXAMPLE 4: 6 taxa (2 identical sampling dates) with 3 demes and simple branch lengths");
try {
treeModel = createSpecifiedTree("(((AB192965Japan2004:4.5,((AF189018Indonesia2005:1.5,AJ842306France2004:0.5):0.25,FM212660Cameroon2007:3.75):3.75):5.0,AF071228Spain1997:2.5):1.0,AF105975Portugal1995:1.5)");
} catch (Exception e) {
throw new RuntimeException("Unable to parse Newick tree");
}
try {
structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
} catch (TreeUtils.MissingTaxonException missing) {
System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
}
System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
System.out.println("EXAMPLE 5: 6 taxa (2 identical sampling dates) with 3 demes and simple (but different than in EXAMPLE 4) branch lengths");
try {
treeModel = createSpecifiedTree("(((((AF189018Indonesia2005:1.5,AJ842306France2004:0.5):0.75,AB192965Japan2004:1.25):3.25,FM212660Cameroon2007:7.5):5.0,AF071228Spain1997:2.5):1.0,AF105975Portugal1995:1.5)");
} catch (Exception e) {
throw new RuntimeException("Unable to parse Newick tree");
}
try {
structured = new StructuredCoalescentLikelihood(treeModel, branchRateModel, popSizesParameter, patterns, null, "", migrationModel, subIntervals, includeSubtree, excludeSubtrees, false);
} catch (TreeUtils.MissingTaxonException missing) {
System.out.println("Error thrown in test class dr.evomodel.coalescent.basta.SCLikelihoodTester: " + missing);
}
System.out.println("Structured coalescent lnL = " + structured.calculateLogLikelihood() + "\n");
}
use of dr.evolution.alignment.SimpleSiteList in project beast-mcmc by beast-dev.
the class AttributePatternsParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
String attributeName = xo.getStringAttribute(ATTRIBUTE);
String secondaryAttributeName = xo.getAttribute(SECONDARY_ATTRIBUTE, (String) null);
TaxonList taxa = (TaxonList) xo.getChild(TaxonList.class);
DataType dataType = DataTypeUtils.getDataType(xo);
if (dataType == null) {
throw new XMLParseException("dataType expected for attributePatterns element");
}
// using a SimpleSiteList rather than Patterns to allow ancestral reconstruction
SimpleSiteList patterns = new SimpleSiteList(dataType, taxa);
int[] pattern = new int[taxa.getTaxonCount()];
boolean attributeFound = true;
for (int i = 0; i < taxa.getTaxonCount(); i++) {
Taxon taxon = taxa.getTaxon(i);
if (secondaryAttributeName == null || secondaryAttributeName.isEmpty()) {
Object value = taxon.getAttribute(attributeName);
if (value != null) {
int state = dataType.getState(value.toString());
if (state < 0) {
throw new XMLParseException("State for attribute, " + attributeName + ", in taxon, " + taxon.getId() + ", is unknown: " + value.toString());
}
pattern[i] = state;
} else {
attributeFound = false;
}
} else {
Object value1 = taxon.getAttribute(attributeName);
Object value2 = taxon.getAttribute(secondaryAttributeName);
if (value1 != null && value2 != null) {
String code = value1.toString() + CompositeDataTypeParser.COMPOSITE_STATE_SEPARATOR + value2.toString();
int state = dataType.getState(code);
if (state < 0) {
throw new XMLParseException("State for attributes, " + attributeName + " & " + secondaryAttributeName + ", in taxon, " + taxon.getId() + ", is unknown: " + code);
}
pattern[i] = state;
} else {
attributeFound = false;
}
}
}
if (!attributeFound) {
throw new XMLParseException("The attribute, " + attributeName + " was not found in all taxa. Check the name of the attribute.");
}
patterns.addPattern(pattern);
if (xo.hasAttribute(XMLParser.ID)) {
Logger.getLogger("dr.evoxml").info("Read attribute patterns, '" + xo.getId() + "' for attribute, " + attributeName);
} else {
Logger.getLogger("dr.evoxml").info("Read attribute patterns for attribute, " + attributeName);
}
return patterns;
}
Aggregations