use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.
the class MCMCTest method testMCMC.
public void testMCMC() {
// Sub model
// new double[]{0.25, 0.25, 0.25, 0.25});
Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
HKY hky = new HKY(kappa, f);
// siteModel
GammaSiteModel siteModel = new GammaSiteModel(hky);
Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
siteModel.setMutationRateParameter(mu);
// treeLikelihood
SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
// Operators
OperatorSchedule schedule = new SimpleOperatorSchedule();
MCMCOperator operator = new ScaleOperator(kappa, 0.5);
operator.setWeight(1.0);
schedule.addOperator(operator);
// Parameter rootParameter = treeModel.createNodeHeightsParameter(true, false, false);
// ScaleOperator scaleOperator = new ScaleOperator(rootParameter, 0.75, AdaptationMode.ADAPTATION_ON, 1.0);
Parameter rootHeight = ((DefaultTreeModel) treeModel).getRootHeightParameter();
rootHeight.setId(TREE_HEIGHT);
operator = new ScaleOperator(rootHeight, 0.5);
operator.setWeight(1.0);
schedule.addOperator(operator);
Parameter internalHeights = ((DefaultTreeModel) treeModel).createNodeHeightsParameter(false, true, false);
operator = new UniformOperator(internalHeights, 10.0);
schedule.addOperator(operator);
operator = new SubtreeSlideOperator(((DefaultTreeModel) treeModel), 1, 1, true, false, false, false, AdaptationMode.ADAPTATION_ON, AdaptableMCMCOperator.DEFAULT_ADAPTATION_TARGET);
schedule.addOperator(operator);
operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
// operator.doOperation();
schedule.addOperator(operator);
operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
// operator.doOperation();
schedule.addOperator(operator);
operator = new WilsonBalding(treeModel, 1.0);
// operator.doOperation();
schedule.addOperator(operator);
// Log
ArrayLogFormatter formatter = new ArrayLogFormatter(false);
MCLogger[] loggers = new MCLogger[2];
loggers[0] = new MCLogger(formatter, 1000, false);
loggers[0].add(treeLikelihood);
loggers[0].add(rootHeight);
loggers[0].add(kappa);
loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
loggers[1].add(treeLikelihood);
loggers[1].add(rootHeight);
loggers[1].add(kappa);
// MCMC
MCMC mcmc = new MCMC("mcmc1");
MCMCOptions options = new MCMCOptions(10000000);
mcmc.setShowOperatorAnalysis(true);
mcmc.init(options, treeLikelihood, schedule, loggers);
mcmc.run();
// time
System.out.println(mcmc.getTimer().toString());
// Tracer
List<Trace> traces = formatter.getTraces();
ArrayTraceList traceList = new ArrayTraceList("MCMCTest", traces, 0);
for (int i = 1; i < traces.size(); i++) {
traceList.analyseTrace(i);
}
// <expectation name="likelihood" value="-1815.75"/>
// <expectation name="treeModel.rootHeight" value="6.42048E-2"/>
// <expectation name="hky.kappa" value="32.8941"/>
TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.75);
TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
assertExpectation(TREE_HEIGHT, treeHeightStats, 6.42048E-2);
TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
assertExpectation(HKYParser.KAPPA, kappaStats, 32.8941);
}
use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.
the class TraceCorrelationAssert method createTreeModel.
private void createTreeModel(ConstantPopulation constant) {
CoalescentSimulator simulator = new CoalescentSimulator();
Tree tree = simulator.simulateTree(alignment, constant);
// treeModel
treeModel = new DefaultTreeModel(tree);
}
use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.
the class BASTADensityTester method createSpecifiedTree.
private static TreeModel createSpecifiedTree(String t) throws Exception {
NewickImporter importer = new NewickImporter(t);
Tree tree = importer.importTree(null);
// treeModel
return new DefaultTreeModel(tree);
}
use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.
the class TreeModelParser method parseXMLObject.
/**
* @return a tree object based on the XML element it was passed.
*/
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
Tree tree = (Tree) xo.getChild(Tree.class);
boolean fixHeights = xo.getAttribute(FIX_HEIGHTS, false);
boolean fixTree = xo.getAttribute(FIX_TREE, false);
DefaultTreeModel treeModel = new DefaultTreeModel(xo.getId(), tree, fixHeights, fixTree);
Logger.getLogger("dr.evomodel").info("\nCreating the tree model, '" + xo.getId() + "'");
for (int i = 0; i < xo.getChildCount(); i++) {
if (xo.getChild(i) instanceof XMLObject) {
XMLObject cxo = (XMLObject) xo.getChild(i);
if (cxo.getName().equals(ROOT_HEIGHT)) {
ParameterParser.replaceParameter(cxo, treeModel.getRootHeightParameter());
} else if (cxo.getName().equals(LEAF_HEIGHT)) {
String taxonName;
if (cxo.hasAttribute(TAXON)) {
taxonName = cxo.getStringAttribute(TAXON);
} else {
throw new XMLParseException("taxa element missing from leafHeight element in treeModel element");
}
int index = treeModel.getTaxonIndex(taxonName);
if (index == -1) {
throw new XMLParseException("taxon " + taxonName + " not found for leafHeight element in treeModel element");
}
NodeRef node = treeModel.getExternalNode(index);
Parameter newParameter = treeModel.getLeafHeightParameter(node);
ParameterParser.replaceParameter(cxo, newParameter);
Taxon taxon = treeModel.getTaxon(index);
setUncertaintyBounds(newParameter, taxon);
} else if (cxo.getName().equals(LEAF_HEIGHTS)) {
// get a set of leaf height parameters out as a compound parameter...
TaxonList taxa = (TaxonList) cxo.getChild(TaxonList.class);
Parameter offsetParameter = (Parameter) cxo.getChild(Parameter.class);
CompoundParameter leafHeights = new CompoundParameter("leafHeights");
for (Taxon taxon : taxa) {
int index = treeModel.getTaxonIndex(taxon);
if (index == -1) {
throw new XMLParseException("taxon " + taxon.getId() + " not found for leafHeight element in treeModel element");
}
NodeRef node = treeModel.getExternalNode(index);
Parameter newParameter = treeModel.getLeafHeightParameter(node);
leafHeights.addParameter(newParameter);
setUncertaintyBounds(newParameter, taxon);
}
ParameterParser.replaceParameter(cxo, leafHeights);
} else if (cxo.getName().equals(NODE_HEIGHTS)) {
boolean rootNode = cxo.getAttribute(ROOT_NODE, false);
boolean internalNodes = cxo.getAttribute(INTERNAL_NODES, false);
boolean leafNodes = cxo.getAttribute(LEAF_NODES, false);
if (!rootNode && !internalNodes && !leafNodes) {
throw new XMLParseException("one or more of root, internal or leaf nodes must be selected for the nodeHeights element");
}
ParameterParser.replaceParameter(cxo, treeModel.createNodeHeightsParameter(rootNode, internalNodes, leafNodes));
} else if (cxo.getName().equals(NODE_RATES)) {
boolean rootNode = cxo.getAttribute(ROOT_NODE, false);
boolean internalNodes = cxo.getAttribute(INTERNAL_NODES, false);
boolean leafNodes = cxo.getAttribute(LEAF_NODES, false);
double[] initialValues = null;
if (cxo.hasAttribute(INITIAL_VALUE)) {
initialValues = cxo.getDoubleArrayAttribute(INITIAL_VALUE);
}
if (!rootNode && !internalNodes && !leafNodes) {
throw new XMLParseException("one or more of root, internal or leaf nodes must be selected for the nodeRates element");
}
ParameterParser.replaceParameter(cxo, treeModel.createNodeRatesParameter(initialValues, rootNode, internalNodes, leafNodes));
} else if (cxo.getName().equals(NODE_TRAITS)) {
parseNodeTraits(cxo, treeModel);
} else if (cxo.getName().equals(LEAF_TRAIT)) {
String name = cxo.getAttribute(NAME, "trait");
String taxonName;
if (cxo.hasAttribute(TAXON)) {
taxonName = cxo.getStringAttribute(TAXON);
} else {
throw new XMLParseException("taxa element missing from leafTrait element in treeModel element");
}
int index = treeModel.getTaxonIndex(taxonName);
if (index == -1) {
throw new XMLParseException("taxon '" + taxonName + "' not found for leafTrait element in treeModel element");
}
NodeRef node = treeModel.getExternalNode(index);
Parameter parameter = treeModel.getNodeTraitParameter(node, name);
if (parameter == null)
throw new XMLParseException("trait '" + name + "' not found for leafTrait (taxon, " + taxonName + ") element in treeModel element");
ParameterParser.replaceParameter(cxo, parameter);
} else {
throw new XMLParseException("illegal child element in " + getParserName() + ": " + cxo.getName());
}
} else if (xo.getChild(i) instanceof Tree) {
// do nothing - already handled
} else {
throw new XMLParseException("illegal child element in " + getParserName() + ": " + xo.getChildName(i) + " " + xo.getChild(i));
}
}
// Logger.getLogger("dr.evomodel").info(" initial tree topology = " + TreeUtils.uniqueNewick(treeModel, treeModel.getRoot()));
Logger.getLogger("dr.evomodel").info(" taxon count = " + treeModel.getExternalNodeCount());
Logger.getLogger("dr.evomodel").info(" tree height = " + treeModel.getNodeHeight(treeModel.getRoot()));
return treeModel;
}
use of dr.evomodel.tree.DefaultTreeModel in project beast-mcmc by beast-dev.
the class MRCATraitStatisticParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
String name = xo.getAttribute(NAME, xo.getId());
String trait = xo.getStringAttribute(TRAIT);
DefaultTreeModel tree = (DefaultTreeModel) xo.getChild(DefaultTreeModel.class);
TaxonList taxa = (TaxonList) xo.getElementFirstChild(MRCA);
try {
return new MRCATraitStatistic(name, trait, tree, taxa);
} catch (TreeUtils.MissingTaxonException mte) {
throw new XMLParseException("Taxon, " + mte + ", in " + getParserName() + "was not found in the tree.");
}
}
Aggregations