Search in sources :

Example 36 with DefaultTreeModel

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);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) ExchangeOperator(dr.evomodel.operators.ExchangeOperator) MCMC(dr.inference.mcmc.MCMC) SubtreeSlideOperator(dr.evomodel.operators.SubtreeSlideOperator) MCMCOptions(dr.inference.mcmc.MCMCOptions) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) WilsonBalding(dr.evomodel.operators.WilsonBalding) TraceCorrelation(dr.inference.trace.TraceCorrelation) SitePatterns(dr.evolution.alignment.SitePatterns) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) Trace(dr.inference.trace.Trace) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) ArrayTraceList(dr.inference.trace.ArrayTraceList) HKY(dr.oldevomodel.substmodel.HKY) Parameter(dr.inference.model.Parameter) MCLogger(dr.inference.loggers.MCLogger)

Example 37 with DefaultTreeModel

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);
}
Also used : SimpleTree(dr.evolution.tree.SimpleTree) Tree(dr.evolution.tree.Tree) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) CoalescentSimulator(dr.evolution.coalescent.CoalescentSimulator)

Example 38 with DefaultTreeModel

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);
}
Also used : NewickImporter(dr.evolution.io.NewickImporter) Tree(dr.evolution.tree.Tree) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel)

Example 39 with DefaultTreeModel

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;
}
Also used : TaxonList(dr.evolution.util.TaxonList) Taxon(dr.evolution.util.Taxon) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) CompoundParameter(dr.inference.model.CompoundParameter) NodeRef(dr.evolution.tree.NodeRef) Tree(dr.evolution.tree.Tree) CompoundParameter(dr.inference.model.CompoundParameter) Parameter(dr.inference.model.Parameter)

Example 40 with DefaultTreeModel

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.");
    }
}
Also used : MRCATraitStatistic(dr.evomodel.tree.MRCATraitStatistic) TaxonList(dr.evolution.util.TaxonList) DefaultTreeModel(dr.evomodel.tree.DefaultTreeModel) TreeUtils(dr.evolution.tree.TreeUtils)

Aggregations

DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)46 Tree (dr.evolution.tree.Tree)21 NewickImporter (dr.evolution.io.NewickImporter)19 Parameter (dr.inference.model.Parameter)19 ArrayList (java.util.ArrayList)14 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)10 TreeModel (dr.evomodel.tree.TreeModel)10 GammaSiteModel (dr.oldevomodel.sitemodel.GammaSiteModel)10 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)9 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)9 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)9 TreeLikelihood (dr.oldevomodel.treelikelihood.TreeLikelihood)9 BeagleSequenceSimulator (dr.app.beagle.tools.BeagleSequenceSimulator)8 Partition (dr.app.beagle.tools.Partition)8 ImportException (dr.evolution.io.Importer.ImportException)8 Taxa (dr.evolution.util.Taxa)8 Taxon (dr.evolution.util.Taxon)8 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)8 IOException (java.io.IOException)8 ExchangeOperator (dr.evomodel.operators.ExchangeOperator)7