Examples with AncestralStateBeagleTreeLikelihood dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood used on opensource projects

Example 1 with AncestralStateBeagleTreeLikelihood

use of dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood in project beast-mcmc by beast-dev.

the class AncestralStateTreeLikelihoodParser method createTreeLikelihood.

protected BeagleTreeLikelihood createTreeLikelihood(//
PatternList patternList, //
TreeModel treeModel, //
BranchModel branchModel, //
GammaSiteRateModel siteRateModel, //
BranchRateModel branchRateModel, //
TipStatesModel tipStatesModel, //
boolean useAmbiguities, //
PartialsRescalingScheme scalingScheme, boolean delayScaling, Map<//
Set<String>, Parameter> partialsRestrictions, //
XMLObject xo) throws XMLParseException {
    //		System.err.println("XML object: " + xo.toString());
    DataType dataType = branchModel.getRootSubstitutionModel().getDataType();
    // default tag is RECONSTRUCTION_TAG
    String tag = xo.getAttribute(RECONSTRUCTION_TAG_NAME, RECONSTRUCTION_TAG);
    boolean useMAP = xo.getAttribute(MAP_RECONSTRUCTION, false);
    boolean useMarginalLogLikelihood = xo.getAttribute(MARGINAL_LIKELIHOOD, true);
    if (patternList.areUnique()) {
        throw new XMLParseException("Ancestral state reconstruction cannot be used with compressed (unique) patterns.");
    }
    return new // Current just returns a OldBeagleTreeLikelihood
    AncestralStateBeagleTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, dataType, tag, useMAP, useMarginalLogLikelihood);
}

Example 2 with AncestralStateBeagleTreeLikelihood

use of dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood in project beast-mcmc by beast-dev.

the class CodonPartitionedRobustCountingParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    AncestralStateBeagleTreeLikelihood[] partition = new AncestralStateBeagleTreeLikelihood[3];
    String[] labels = new String[] { FIRST, SECOND, THIRD };
    int patternCount = -1;
    BranchRateModel testBranchRateModel = null;
    for (int i = 0; i < 3; i++) {
        partition[i] = (AncestralStateBeagleTreeLikelihood) xo.getChild(labels[i]).getChild(AncestralStateBeagleTreeLikelihood.class);
        if (i == 0) {
            patternCount = partition[i].getPatternCount();
        } else {
            if (partition[i].getPatternCount() != patternCount) {
                throw new XMLParseException("Codon-partitioned robust counting requires all partitions to have the same length." + " Make sure that partitions include all unique sites and do not strip gaps.");
            }
        }
        // Ensure that siteRateModel has one category
        if (partition[i].getSiteRateModel().getCategoryCount() > 1) {
            throw new XMLParseException("Robust counting currently only implemented for single category models");
        }
        // Ensure that branchRateModel is the same across all partitions
        if (testBranchRateModel == null) {
            testBranchRateModel = partition[i].getBranchRateModel();
        } else if (testBranchRateModel != partition[i].getBranchRateModel()) {
            throw new XMLParseException("Robust counting currently requires the same branch rate model for all partitions");
        }
    }
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    Codons codons = Codons.UNIVERSAL;
    if (xo.hasAttribute(GeneticCode.GENETIC_CODE)) {
        String codeStr = xo.getStringAttribute(GeneticCode.GENETIC_CODE);
        codons = Codons.findByName(codeStr);
    }
    String labelingString = (String) xo.getAttribute(LABELING);
    CodonLabeling codonLabeling = CodonLabeling.parseFromString(labelingString);
    if (codonLabeling == null) {
        throw new XMLParseException("Unrecognized codon labeling '" + labelingString + "'");
    }
    String branchFormatString = xo.getAttribute(BRANCH_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat branchFormat = StratifiedTraitOutputFormat.parseFromString(branchFormatString);
    if (branchFormat == null) {
        throw new XMLParseException("Unrecognized branch output format '" + branchFormat + "'");
    }
    String logFormatString = xo.getAttribute(LOG_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat logFormat = StratifiedTraitOutputFormat.parseFromString(logFormatString);
    if (logFormat == null) {
        throw new XMLParseException("Unrecognized log output format '" + branchFormat + "'");
    }
    boolean useUniformization = xo.getAttribute(USE_UNIFORMIZATION, false);
    boolean includeExternalBranches = xo.getAttribute(INCLUDE_EXTERNAL, true);
    boolean includeInternalBranches = xo.getAttribute(INCLUDE_INTERNAL, true);
    boolean doUnconditionedPerBranch = xo.getAttribute(DO_UNCONDITIONED_PER_BRANCH, false);
    boolean averageRates = xo.getAttribute(AVERAGE_RATES, true);
    boolean saveCompleteHistory = xo.getAttribute(SAVE_HISTORY, false);
    boolean useNewNeutralModel = xo.getAttribute(USE_NEW_NEUTRAL_MODEL, false);
    String prefix = xo.hasAttribute(PREFIX) ? xo.getStringAttribute(PREFIX) : null;
    return new CodonPartitionedRobustCounting(xo.getId(), tree, partition, codons, codonLabeling, useUniformization, includeExternalBranches, includeInternalBranches, doUnconditionedPerBranch, saveCompleteHistory, averageRates, useNewNeutralModel, branchFormat, logFormat, prefix);
}

Example 3 with AncestralStateBeagleTreeLikelihood

use of dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood in project beast-mcmc by beast-dev.

the class AncestralSequenceAnnotator method sampleTree.

private FlexibleTree sampleTree(Tree tree, PatternList alignment, GammaSiteRateModel siteModel, BranchRateModel rateModel) {
    FlexibleTree flexTree = new FlexibleTree(tree, true);
    flexTree.adoptTreeModelOrdering();
    FlexibleTree finalTree = new FlexibleTree(tree);
    finalTree.adoptTreeModelOrdering();
    TreeModel treeModel = new DefaultTreeModel(tree);
    // Turn off noisy logging by TreeLikelihood constructor
    Logger logger = Logger.getLogger("dr.evomodel");
    boolean useParentHandlers = logger.getUseParentHandlers();
    logger.setUseParentHandlers(false);
    // AncestralStateTreeLikelihood likelihood = new AncestralStateTreeLikelihood(
    // alignment,
    // treeModel,
    // siteModel,
    // rateModel,
    // false, true,
    // alignment.getDataType(),
    // TAG,
    // false);
    AncestralStateBeagleTreeLikelihood likelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, new HomogeneousBranchModel(siteModel.getSubstitutionModel()), siteModel, rateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, alignment.getDataType(), TAG, false, true);
    // PatternList patternList, TreeModel treeModel,
    // BranchSiteModel branchSiteModel,
    // SiteRateModel siteRateModel,
    // BranchRateModel branchRateModel,
    // boolean useAmbiguities,
    // PartialsRescalingScheme scalingScheme,
    // Map<Set<String>, Parameter> partialsRestrictions,
    // final DataType dataType,
    // final String tag,
    // SubstitutionModel substModel,
    // boolean useMAP,
    // boolean returnML) {
    // PatternList patternList, TreeModel treeModel,
    // SiteModel siteModel, BranchRateModel branchRateModel,
    // boolean useAmbiguities, boolean storePartials,
    // final DataType dataType,
    // final String tag,
    // boolean forceRescaling,
    // boolean useMAP,
    // boolean returnML) {
    logger.setUseParentHandlers(useParentHandlers);
    // Sample internal nodes
    likelihood.makeDirty();
    double logLikelihood = likelihood.getLogLikelihood();
    System.out.println("The new and old Likelihood (this value should be roughly the same, debug?): " + logLikelihood + ", " + Double.parseDouble(tree.getAttribute(LIKELIHOOD).toString()));
    if (Double.parseDouble(tree.getAttribute(LIKELIHOOD).toString()) != logLikelihood) {
    /* Newly written check, not sure if this is correct. May need to round values at least */
    // throw new RuntimeException("The values of likelihood are not identical");
    }
    // System.err.printf("New logLikelihood = %4.1f\n", logLikelihood);
    flexTree.setAttribute(LIKELIHOOD, logLikelihood);
    TreeTrait ancestralStates = likelihood.getTreeTrait(TAG);
    for (int i = 0; i < treeModel.getNodeCount(); i++) {
        NodeRef node = treeModel.getNode(i);
        String sample = ancestralStates.getTraitString(treeModel, node);
        String oldSeq = (String) flexTree.getNodeAttribute(flexTree.getNode(i), SEQ_STRING);
        if (oldSeq != null) {
            char[] seq = (sample.substring(1, sample.length() - 1)).toCharArray();
            int length = oldSeq.length();
            for (int j = 0; j < length; j++) {
                if (oldSeq.charAt(j) == GAP)
                    seq[j] = GAP;
            }
            String newSeq = new String(seq);
            // if( newSeq.contains("MMMMMMM") ) {
            // System.err.println("bad = "+newSeq);
            // System.exit(-1);
            // }
            finalTree.setNodeAttribute(finalTree.getNode(i), NEW_SEQ, newSeq);
        }
    // Taxon taxon = finalTree.getNodeTaxon(finalTree.getNode(i));
    // System.err.println("node: "+(taxon == null ? "null" : taxon.getId())+" "+
    // finalTree.getNodeAttribute(finalTree.getNode(i),NEW_SEQ));
    }
    return finalTree;
}

Example 4 with AncestralStateBeagleTreeLikelihood

use of dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood in project beast-mcmc by beast-dev.

the class TipStateSwapOperatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    AncestralStateBeagleTreeLikelihood treeLikelihood = (AncestralStateBeagleTreeLikelihood) xo.getChild(AncestralStateBeagleTreeLikelihood.class);
    final double weight = xo.getDoubleAttribute("weight");
    final boolean uniformRandomization = xo.getAttribute(UNIFORM_RANDOMIZATION, false);
    return new TipStateSwapOperator(treeLikelihood, weight, uniformRandomization);
}

Example 5 with AncestralStateBeagleTreeLikelihood

use of dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood in project beast-mcmc by beast-dev.

the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new DefaultTreeModel("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);
}