Search in sources :

Example 11 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class BalancedBeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(BeagleTreeLikelihoodParser.USE_AMBIGUITIES, false);
    /*int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
        if (instanceCount < 1) {
            instanceCount = 1;
        }

        String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
        if (ic != null && ic.length() > 0) {
            instanceCount = Integer.parseInt(ic);
        }*/
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    // if (xo.getChild(TipStatesModel.class) != null) {
    // throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    // }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    if (xo.hasAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME)) {
        // scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    boolean delayScaling = true;
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        // Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    /*if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
            return createTreeLikelihood(
                    patternList,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    tipStatesModel,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo
            );
        }*/
    // first run a test for instanceCount == 1
    System.err.println("\nTesting instanceCount == 1");
    Likelihood baseLikelihood = createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    double start = System.nanoTime();
    for (int i = 0; i < TEST_RUNS; i++) {
        baseLikelihood.makeDirty();
        baseLikelihood.getLogLikelihood();
    }
    double end = System.nanoTime();
    double baseResult = end - start;
    System.err.println("Evaluation took: " + baseResult);
    if (!(patternList instanceof SitePatterns)) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with BEAUti-selected codon partitioning.");
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    // List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    List<Likelihood> likelihoods = null;
    CompoundLikelihood compound = null;
    int instanceCount = 2;
    boolean optimal = false;
    while (optimal == false) {
        System.err.println("\nCreating instanceCount == " + instanceCount);
        likelihoods = new ArrayList<Likelihood>();
        for (int i = 0; i < instanceCount; i++) {
            Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }
        // construct compoundLikelihood
        compound = new CompoundLikelihood(instanceCount, likelihoods);
        // test timings
        System.err.println("\nTesting instanceCount == " + instanceCount);
        start = System.nanoTime();
        for (int i = 0; i < TEST_RUNS; i++) {
            compound.makeDirty();
            compound.getLogLikelihood();
        }
        end = System.nanoTime();
        double newResult = end - start;
        System.err.println("Evaluation took: " + newResult);
        if (baseResult / newResult > TEST_CUTOFF) {
            instanceCount++;
            baseResult = newResult;
        } else {
            optimal = true;
            instanceCount--;
            System.err.println("\nCreating final BeagleTreeLikelihood with instanceCount: " + instanceCount);
            likelihoods = new ArrayList<Likelihood>();
            for (int i = 0; i < instanceCount; i++) {
                Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
                AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
                treeLikelihood.setId(xo.getId() + "_" + instanceCount);
                likelihoods.add(treeLikelihood);
            }
            // construct compoundLikelihood
            compound = new CompoundLikelihood(instanceCount, likelihoods);
        }
    }
    return compound;
/*for (int i = 0; i < instanceCount; i++) {

            Patterns subPatterns = new Patterns((SitePatterns)patternList, 0, 0, 1, i, instanceCount);

            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(
                    subPatterns,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    null,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }

        return new CompoundLikelihood(likelihoods);*/
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Set(java.util.Set) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) BeagleTreeLikelihood(dr.evomodel.treelikelihood.BeagleTreeLikelihood) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PatternList(dr.evolution.alignment.PatternList) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) TreeModel(dr.evomodel.tree.TreeModel) Patterns(dr.evolution.alignment.Patterns) SitePatterns(dr.evolution.alignment.SitePatterns) TreeUtils(dr.evolution.tree.TreeUtils) SitePatterns(dr.evolution.alignment.SitePatterns) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Parameter(dr.inference.model.Parameter)

Example 12 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    if (instanceCount < 1) {
        instanceCount = 1;
    }
    String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    if (ic != null && ic.length() > 0) {
        instanceCount = Integer.parseInt(ic);
    }
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchModel instanceof EpochBranchModel && rootFreqModel != null) {
        EpochBranchModel epochBranchModel = (EpochBranchModel) branchModel;
        epochBranchModel.setRootFrequencyModel(rootFreqModel);
    }
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    // if (xo.getChild(TipStatesModel.class) != null) {
    // throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    // }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    boolean delayScaling = true;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        // Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    int beagleThreadCount = -1;
    if (System.getProperty(BEAGLE_THREAD_COUNT) != null) {
        beagleThreadCount = Integer.parseInt(System.getProperty(BEAGLE_THREAD_COUNT));
    }
    if (beagleThreadCount == -1) {
        // the default is -1 threads (automatic thread pool size) but an XML attribute can override it
        int threadCount = xo.getAttribute(THREADS, -1);
        if (System.getProperty(THREAD_COUNT) != null) {
            threadCount = Integer.parseInt(System.getProperty(THREAD_COUNT));
        }
        // Todo: allow for different number of threads per beagle instance according to pattern counts
        if (threadCount >= 0) {
            System.setProperty(BEAGLE_THREAD_COUNT, Integer.toString(threadCount / instanceCount));
        }
    }
    if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
        return createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    for (int i = 0; i < instanceCount; i++) {
        Patterns subPatterns = new Patterns(patternList, i, instanceCount);
        AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
        treeLikelihood.setId(xo.getId() + "_" + instanceCount);
        likelihoods.add(treeLikelihood);
    }
    return new CompoundLikelihood(likelihoods);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Set(java.util.Set) CompoundLikelihood(dr.inference.model.CompoundLikelihood) Likelihood(dr.inference.model.Likelihood) BeagleTreeLikelihood(dr.evomodel.treelikelihood.BeagleTreeLikelihood) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) AbstractTreeLikelihood(dr.evomodel.treelikelihood.AbstractTreeLikelihood) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) BranchModel(dr.evomodel.branchmodel.BranchModel) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) MutableTreeModel(dr.evolution.tree.MutableTreeModel) Patterns(dr.evolution.alignment.Patterns) TreeUtils(dr.evolution.tree.TreeUtils) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) TaxonList(dr.evolution.util.TaxonList) CompoundLikelihood(dr.inference.model.CompoundLikelihood) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Parameter(dr.inference.model.Parameter)

Example 13 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class TreeDataLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    boolean usePreOrder = xo.getAttribute(USE_PREORDER, false);
    boolean branchRateDerivative = xo.getAttribute(BRANCHRATE_DERIVATIVE, usePreOrder);
    boolean branchInfinitesimalDerivative = xo.getAttribute(BRANCHINFINITESIMAL_DERIVATIVE, false);
    if (usePreOrder != (branchRateDerivative || branchInfinitesimalDerivative)) {
        throw new RuntimeException("Need to specify derivative types.");
    }
    PreOrderSettings settings = new PreOrderSettings(usePreOrder, branchRateDerivative, branchInfinitesimalDerivative);
    // TreeDataLikelihood doesn't currently support Instances defined from the command line
    // int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    // if (instanceCount < 1) {
    // instanceCount = 1;
    // }
    // 
    // String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    // if (ic != null && ic.length() > 0) {
    // instanceCount = Integer.parseInt(ic);
    // }
    List<PatternList> patternLists = new ArrayList<PatternList>();
    List<SiteRateModel> siteRateModels = new ArrayList<SiteRateModel>();
    List<BranchModel> branchModels = new ArrayList<BranchModel>();
    boolean hasSinglePartition = false;
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    if (patternList != null) {
        hasSinglePartition = true;
        patternLists.add(patternList);
        GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
        siteRateModels.add(siteRateModel);
        FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
        BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
        if (branchModel == null) {
            SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
            if (substitutionModel == null) {
                substitutionModel = siteRateModel.getSubstitutionModel();
            }
            if (substitutionModel == null) {
                throw new XMLParseException("No substitution model available for partition in DataTreeLikelihood: " + xo.getId());
            }
            branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
        }
        branchModels.add(branchModel);
    }
    int k = 0;
    for (int i = 0; i < xo.getChildCount(); i++) {
        if (xo.getChildName(i).equals(PARTITION)) {
            if (hasSinglePartition) {
                throw new XMLParseException("Either a single set of patterns should be given or multiple 'partitions' elements within DataTreeLikelihood: " + xo.getId());
            }
            k += 1;
            XMLObject cxo = (XMLObject) xo.getChild(i);
            patternList = (PatternList) cxo.getChild(PatternList.class);
            patternLists.add(patternList);
            GammaSiteRateModel siteRateModel = (GammaSiteRateModel) cxo.getChild(GammaSiteRateModel.class);
            siteRateModels.add(siteRateModel);
            FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
            BranchModel branchModel = (BranchModel) cxo.getChild(BranchModel.class);
            if (branchModel == null) {
                SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
                if (substitutionModel == null) {
                    substitutionModel = siteRateModel.getSubstitutionModel();
                }
                if (substitutionModel == null) {
                    throw new XMLParseException("No substitution model available for partition " + k + " in DataTreeLikelihood: " + xo.getId());
                }
                branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
            }
            branchModels.add(branchModel);
            BranchRateModel branchRateModel = (BranchRateModel) cxo.getChild(BranchRateModel.class);
            if (branchRateModel != null) {
                throw new XMLParseException("Partitions are not currently allowed their own BranchRateModel in TreeDataLikelihood object '" + xo.getId());
            }
        }
    }
    if (patternLists.size() == 0) {
        throw new XMLParseException("Either a single set of patterns should be given or multiple 'partitions' elements within DataTreeLikelihood: " + xo.getId());
    }
    Tree treeModel = (Tree) xo.getChild(Tree.class);
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchRateModel == null) {
        throw new XMLParseException("BranchRateModel missing from TreeDataLikelihood object '" + xo.getId());
    // branchRateModel = new DefaultBranchRateModel();
    }
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    final boolean preferGPU = xo.getAttribute(PREFER_GPU, false);
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "TreeDataLikelihood object '" + xo.getId());
    }
    final boolean delayScaling = xo.getAttribute(DELAY_SCALING, true);
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    return createTreeDataLikelihood(patternLists, branchModels, siteRateModels, treeModel, branchRateModel, null, useAmbiguities, preferGPU, scalingScheme, delayScaling, settings);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SiteRateModel(dr.evomodel.siteratemodel.SiteRateModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) PreOrderSettings(dr.evomodel.treedatalikelihood.PreOrderSettings) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Tree(dr.evolution.tree.Tree)

Example 14 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class MultiPartitionDataLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    if (instanceCount < 1) {
        instanceCount = 1;
    }
    String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    if (ic != null && ic.length() > 0) {
        instanceCount = Integer.parseInt(ic);
    }
    if (DEBUG) {
        System.out.println("instanceCount: " + instanceCount);
    }
    List<PatternList> patternLists = xo.getAllChildren(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    List<SiteRateModel> siteRateModels = xo.getAllChildren(SiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    List<BranchModel> branchModels = xo.getAllChildren(BranchModel.class);
    if (branchModels.size() == 0) {
        if (DEBUG) {
            System.out.println("branchModels == null");
        }
        branchModels = new ArrayList<BranchModel>();
        List<SubstitutionModel> substitutionModels = xo.getAllChildren(SubstitutionModel.class);
        if (substitutionModels.size() == 0) {
            if (DEBUG) {
                System.out.println("substitutionModels == null");
            }
            for (SiteRateModel siteRateModel : siteRateModels) {
                SubstitutionModel substitutionModel = ((GammaSiteRateModel) siteRateModel).getSubstitutionModel();
                if (substitutionModel == null) {
                    throw new XMLParseException("No substitution model available for TreeDataLikelihood: " + xo.getId());
                }
                branchModels.add(new HomogeneousBranchModel(substitutionModel, rootFreqModel));
            }
        }
        if (DEBUG) {
            System.out.println("branchModels size: " + branchModels.size());
        }
        for (BranchModel branchModel : branchModels) {
            System.out.println("  " + branchModel.getId() + "  " + branchModel.getModelName());
        }
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchRateModel == null) {
        branchRateModel = new DefaultBranchRateModel();
    }
    if (DEBUG) {
        System.out.println("BranchRateModel: " + branchRateModel.getId());
    }
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    boolean delayScaling = true;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    if (instanceCount == 1) {
        if (DEBUG) {
            System.out.println("instanceCount == 1");
        }
        return createTreeDataLikelihood(patternLists, treeModel, branchModels, siteRateModels, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, xo);
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    List<PatternList> patternInstanceLists = new ArrayList<PatternList>();
    for (int j = 0; j < patternLists.size(); j++) {
        for (int i = 0; i < instanceCount; i++) {
            patternInstanceLists.add(new Patterns(patternLists.get(j), i, instanceCount));
        }
    }
    return createTreeDataLikelihood(patternLists, treeModel, branchModels, siteRateModels, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, xo);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) PartialsRescalingScheme(dr.evomodel.treelikelihood.PartialsRescalingScheme) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) SiteRateModel(dr.evomodel.siteratemodel.SiteRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) TipStatesModel(dr.evomodel.tipstatesmodel.TipStatesModel) TreeModel(dr.evomodel.tree.TreeModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) Patterns(dr.evolution.alignment.Patterns)

Example 15 with SubstitutionModel

use of dr.evomodel.substmodel.SubstitutionModel in project beast-mcmc by beast-dev.

the class MutationDeathModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Parameter dummyFreqParameter;
    Parameter delParam = (Parameter) xo.getChild(Parameter.class);
    Logger.getLogger("dr.evomodel").info("Creating MutationDeath substitution model.\n\tInitial death rate is " + delParam.getParameterValue(0));
    MutationDeathType dT = (MutationDeathType) xo.getChild(MutationDeathType.class);
    SubstitutionModel evoModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
    if (evoModel == null) {
        // Assuming pure survival model
        Logger.getLogger("dr.evomodel").info("\tSubstitutionModel not provided assuming pure death/survival model.");
        dummyFreqParameter = new Parameter.Default(new double[] { 1.0, 0.0 });
    } else {
        dummyFreqParameter = new Parameter.Default(dT.getStateCount());
        double[] freqs = evoModel.getFrequencyModel().getFrequencies();
        for (int i = 0; i < freqs.length; ++i) {
            dummyFreqParameter.setParameterValueQuietly(i, freqs[i]);
        }
        dummyFreqParameter.setParameterValueQuietly(dT.getStateCount() - 1, 0.0);
    }
    FrequencyModel dummyFrequencies = new FrequencyModel(dT, dummyFreqParameter);
    Parameter mutationRate;
    if (xo.hasChildNamed(MUTATION_RATE)) {
        mutationRate = (Parameter) xo.getElementFirstChild(MUTATION_RATE);
    } else {
        mutationRate = new Parameter.Default(new double[] { 1.0 });
    }
    Logger.getLogger("dr.evomodel").info("\tInitial mutation rate is " + mutationRate.getParameterValue(0));
    return new MutationDeathModel(delParam, dT, evoModel, dummyFrequencies, mutationRate);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Parameter(dr.inference.model.Parameter) MutationDeathModel(dr.evomodel.substmodel.MutationDeathModel) MutationDeathType(dr.evolution.datatype.MutationDeathType) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel)

Aggregations

SubstitutionModel (dr.evomodel.substmodel.SubstitutionModel)25 TreeModel (dr.evomodel.tree.TreeModel)13 Parameter (dr.inference.model.Parameter)13 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)11 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)9 ArrayList (java.util.ArrayList)9 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)8 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)8 BranchModel (dr.evomodel.branchmodel.BranchModel)7 XMLObject (dr.xml.XMLObject)6 PatternList (dr.evolution.alignment.PatternList)5 Partition (dr.app.beagle.tools.Partition)4 NodeRef (dr.evolution.tree.NodeRef)4 Tree (dr.evolution.tree.Tree)4 TaxonList (dr.evolution.util.TaxonList)4 EpochBranchModel (dr.evomodel.branchmodel.EpochBranchModel)4 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)4 TipStatesModel (dr.evomodel.tipstatesmodel.TipStatesModel)4 BeagleTreeLikelihood (dr.evomodel.treelikelihood.BeagleTreeLikelihood)4 PartialsRescalingScheme (dr.evomodel.treelikelihood.PartialsRescalingScheme)4