Examples with SubstitutionModel dr.evomodel.substmodel.SubstitutionModel used on opensource projects

Example 16 with SubstitutionModel

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

the class MarkovModulatedSubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (!(dataType instanceof HiddenDataType)) {
        throw new XMLParseException("Must construct " + MARKOV_MODULATED_MODEL + " with hidden data types. " + "You may need to provide the `-universal` extension to your hidden code type.");
    }
    Parameter switchingRates = (Parameter) xo.getElementFirstChild(SWITCHING_RATES);
    List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Object cxo = xo.getChild(i);
        if (cxo instanceof SubstitutionModel) {
            substModels.add((SubstitutionModel) cxo);
        }
    }
    boolean geometricRates = xo.getAttribute(GEOMETRIC_RATES, false);
    Parameter rateScalar = xo.hasChildNamed(RATE_SCALAR) ? (Parameter) xo.getChild(RATE_SCALAR).getChild(Parameter.class) : null;
    SiteRateModel siteRateModel = (SiteRateModel) xo.getChild(SiteRateModel.class);
    if (siteRateModel != null) {
        if (siteRateModel.getCategoryCount() != substModels.size() && substModels.size() % siteRateModel.getCategoryCount() != 0) {
            throw new XMLParseException("Number of gamma categories must equal number of substitution models in " + xo.getId());
        }
    }
    MarkovModulatedSubstitutionModel mmsm = new MarkovModulatedSubstitutionModel(xo.getId(), substModels, switchingRates, dataType, null, rateScalar, geometricRates, siteRateModel);
    if (xo.getAttribute(RENORMALIZE, false)) {
        mmsm.setNormalization(true);
    }
    return mmsm;
}

Example 17 with SubstitutionModel

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

the class CTMCScalePriorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    TaxonList taxa = (TaxonList) xo.getChild(TaxonList.class);
    Parameter ctmcScale = (Parameter) xo.getElementFirstChild(SCALEPARAMETER);
    boolean reciprocal = xo.getAttribute(RECIPROCAL, false);
    boolean trial = xo.getAttribute(TRIAL, false);
    SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
    Logger.getLogger("dr.evolution").info("Creating CTMC Scale Reference Prior model.");
    if (taxa != null) {
        Logger.getLogger("dr.evolution").info("Acting on subtree of size " + taxa.getTaxonCount());
    }
    return new CTMCScalePrior(MODEL_NAME, ctmcScale, treeModel, taxa, reciprocal, substitutionModel, trial);
}

Example 18 with SubstitutionModel

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

the class BeagleBranchLikelihood method main.

// END: LikelihoodColumn class
// ////////////
// ---TEST---//
// ////////////
public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("((SimSeq1:22.0,SimSeq2:22.0):12.0,(SimSeq3:23.1,SimSeq4:23.1):10.899999999999999);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 1.000);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogeneousBranchModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        System.out.println(alignment);
        BeagleTreeLikelihood btl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("BTL(homogeneous) = " + btl.getLogLikelihood());
        BeagleBranchLikelihood bbl = new BeagleBranchLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, freqModel, branchRateModel);
        int branchIndex = 4;
        System.out.println(bbl.getBranchLogLikelihood(branchIndex));
        bbl.finalizeBeagle();
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    } catch (Throwable e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 19 with SubstitutionModel

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

the class PartitionParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    int from = 0;
    int to = -1;
    int every = xo.getAttribute(EVERY, 1);
    DataType dataType = null;
    if (xo.hasAttribute(FROM)) {
        from = xo.getIntegerAttribute(FROM) - 1;
        if (from < 0) {
            throw new XMLParseException("Illegal 'from' attribute in patterns element");
        }
    }
    if (xo.hasAttribute(TO)) {
        to = xo.getIntegerAttribute(TO) - 1;
        if (to < 0 || to < from) {
            throw new XMLParseException("Illegal 'to' attribute in patterns element");
        }
    }
    if (every <= 0) {
        throw new XMLParseException("Illegal 'every' attribute in patterns element");
    }
    if (xo.hasAttribute(DATA_TYPE)) {
        dataType = DataTypeUtils.getDataType(xo);
    }
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    // FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    FrequencyModel freqModel;
    List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Object cxo = xo.getChild(i);
        if (cxo instanceof FrequencyModel) {
            freqModels.add((FrequencyModel) cxo);
        }
    }
    if (freqModels.size() == 1) {
        freqModel = freqModels.get(0);
    } else {
        double[] freqParameter = new double[freqModels.size() * freqModels.get(0).getFrequencyCount()];
        int index = 0;
        for (int i = 0; i < freqModels.size(); i++) {
            for (int j = 0; j < freqModels.get(i).getFrequencyCount(); j++) {
                freqParameter[index] = (freqModels.get(i).getFrequency(j)) / freqModels.size();
                index++;
            }
        }
        freqModel = new FrequencyModel(dataType, freqParameter);
    }
    Sequence rootSequence = (Sequence) xo.getChild(Sequence.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (rateModel == null) {
        rateModel = new DefaultBranchRateModel();
    }
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        branchModel = new HomogeneousBranchModel(substitutionModel);
    }
    Partition partition = new Partition(tree, branchModel, siteModel, rateModel, freqModel, from, to, every, dataType);
    if (rootSequence != null) {
        partition.setRootSequence(rootSequence);
    }
    return partition;
}

Example 20 with SubstitutionModel

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

the class LineageSpecificBranchModel method main.

// END: acceptState
public static void main(String[] args) {
    try {
        // the seed of the BEAST
        MathUtils.setSeed(666);
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        TreeModel tree = new DefaultTreeModel(importer.importTree(null));
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { // 
        0.0163936, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344 });
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        MG94HKYCodonModel mg94 = new MG94K80CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel, new CodonOptions());
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new // 
        Partition(// 
        tree, // 
        substitutionModel, // 
        siteRateModel, // 
        branchRateModel, // 
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
        for (int i = 0; i < 2; i++) {
            // alpha = new Parameter.Default(1, 10 );
            // beta = new Parameter.Default(1, 5 );
            // mg94 = new MG94HKYCodonModel(Codons.UNIVERSAL, alpha, beta,
            // freqModel);
            substModels.add(mg94);
        }
        Parameter uCategories = new Parameter.Default(2, 0);
        // CountableBranchCategoryProvider provider = new CountableBranchCategoryProvider.IndependentBranchCategoryModel(tree, uCategories);
        LineageSpecificBranchModel branchSpecific = new // provider,
        LineageSpecificBranchModel(// provider,
        tree, // provider,
        freqModel, // provider,
        substModels, uCategories);
        BeagleTreeLikelihood like = new // 
        BeagleTreeLikelihood(// 
        convert, // 
        tree, // 
        branchSpecific, // 
        siteRateModel, // 
        branchRateModel, // 
        null, // 
        false, PartialsRescalingScheme.DEFAULT, true);
        BeagleTreeLikelihood gold = new // 
        BeagleTreeLikelihood(// 
        convert, // 
        tree, // 
        substitutionModel, // 
        siteRateModel, // 
        branchRateModel, // 
        null, // 
        false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("likelihood (gold) = " + gold.getLogLikelihood());
        System.out.println("likelihood = " + like.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
    }
}