Examples with DefaultBranchRateModel dr.evomodel.branchratemodel.DefaultBranchRateModel used on opensource projects

Example 6 with DefaultBranchRateModel

use of dr.evomodel.branchratemodel.DefaultBranchRateModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateRandomBranchAssignment.

// END: main
static void simulateRandomBranchAssignment() {
    try {
        System.out.println("Test case I dunno which: simulate random branch assignments");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        File treeFile = new File("/home/filip/Dropbox/BeagleSequenceSimulator/SimTree/SimTree.figtree");
        Tree tree = Utils.importTreeFromFile(treeFile);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create base subst model
        Parameter omegaParameter = new Parameter.Default("omega", 1, 1.0);
        Parameter kappaParameter = new Parameter.Default("kappa", 1, 1.0);
        GY94CodonModel baseSubModel = new GY94CodonModel(Codons.UNIVERSAL, omegaParameter, kappaParameter, freqModel);
        RandomBranchModel substitutionModel = new RandomBranchModel(treeModel, baseSubModel, 0.25, false, -1);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        //			Sequence ancestralSequence = new Sequence();
        //			ancestralSequence.appendSequenceString("TCAAGTGAGG");
        //			partition1.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        SimpleAlignment alignment = simulator.simulate(simulateInPar, true);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
// END: try-catch
}

Example 7 with DefaultBranchRateModel

use of dr.evomodel.branchratemodel.DefaultBranchRateModel in project beast-mcmc by beast-dev.

the class ContinuousDataLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    MultivariateDiffusionModel diffusionModel = (MultivariateDiffusionModel) xo.getChild(MultivariateDiffusionModel.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TreeTraitParserUtilities utilities = new TreeTraitParserUtilities();
    String traitName = TreeTraitParserUtilities.DEFAULT_TRAIT_NAME;
    TreeTraitParserUtilities.TraitsAndMissingIndices returnValue = utilities.parseTraitsFromTaxonAttributes(xo, traitName, treeModel, true);
    CompoundParameter traitParameter = returnValue.traitParameter;
    List<Integer> missingIndices = returnValue.missingIndices;
    Parameter sampleMissingParameter = returnValue.sampleMissingParameter;
    traitName = returnValue.traitName;
    final int dim = diffusionModel.getPrecisionmatrix().length;
    PrecisionType precisionType = PrecisionType.SCALAR;
    if (missingIndices.size() > 0 && !xo.getAttribute(FORCE_COMPLETELY_MISSING, false)) {
        precisionType = PrecisionType.FULL;
    }
    System.err.println("Using precisionType == " + precisionType + " for data model.");
    ContinuousTraitDataModel dataModel = new ContinuousTraitDataModel(traitName, traitParameter, missingIndices, dim, precisionType);
    ConjugateRootTraitPrior rootPrior = ConjugateRootTraitPrior.parseConjugateRootTraitPrior(xo, dim);
    boolean useTreeLength = xo.getAttribute(USE_TREE_LENGTH, false);
    boolean scaleByTime = xo.getAttribute(SCALE_BY_TIME, false);
    if (rateModel == null) {
        rateModel = new DefaultBranchRateModel();
    }
    ContinuousRateTransformation rateTransformation = new ContinuousRateTransformation.Default(treeModel, scaleByTime, useTreeLength);
    ContinuousDataLikelihoodDelegate delegate = new ContinuousDataLikelihoodDelegate(treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel);
    TreeDataLikelihood treeDataLikelihood = new TreeDataLikelihood(delegate, treeModel, rateModel);
    boolean reconstructTraits = xo.getAttribute(RECONSTRUCT_TRAITS, true);
    if (reconstructTraits) {
        if (missingIndices.size() == 0) {
            ProcessSimulationDelegate simulationDelegate = new ProcessSimulationDelegate.ConditionalOnTipsRealizedDelegate(traitName, treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel, delegate);
            TreeTraitProvider traitProvider = new ProcessSimulation(traitName, treeDataLikelihood, simulationDelegate);
            treeDataLikelihood.addTraits(traitProvider.getTreeTraits());
        } else {
            ProcessSimulationDelegate simulationDelegate = delegate.getPrecisionType() == PrecisionType.SCALAR ? new ProcessSimulationDelegate.ConditionalOnTipsRealizedDelegate(traitName, treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel, delegate) : new ProcessSimulationDelegate.MultivariateConditionalOnTipsRealizedDelegate(traitName, treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel, delegate);
            TreeTraitProvider traitProvider = new ProcessSimulation(traitName, treeDataLikelihood, simulationDelegate);
            treeDataLikelihood.addTraits(traitProvider.getTreeTraits());
            ProcessSimulationDelegate fullConditionalDelegate = new ProcessSimulationDelegate.TipRealizedValuesViaFullConditionalDelegate(traitName, treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel, delegate);
            treeDataLikelihood.addTraits(new ProcessSimulation(("fc." + traitName), treeDataLikelihood, fullConditionalDelegate).getTreeTraits());
        //                String partialTraitName = getPartiallyMissingTraitName(traitName);
        //
        //                ProcessSimulationDelegate parialSimulationDelegate = new ProcessSimulationDelegate.ConditionalOnPartiallyMissingTipsDelegate(partialTraitName,
        //                        treeModel, diffusionModel, dataModel, rootPrior, rateTransformation, rateModel, delegate);
        //
        //                TreeTraitProvider partialTraitProvider = new ProcessSimulation(partialTraitName,
        //                        treeDataLikelihood, parialSimulationDelegate);
        //
        //                treeDataLikelihood.addTraits(partialTraitProvider.getTreeTraits());
        }
    }
    return treeDataLikelihood;
}

Example 8 with DefaultBranchRateModel

use of dr.evomodel.branchratemodel.DefaultBranchRateModel 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 == null) {
        if (DEBUG) {
            System.out.println("branchModels == null");
        }
        branchModels = new ArrayList<BranchModel>();
        List<SubstitutionModel> substitutionModels = xo.getAllChildren(SubstitutionModel.class);
        if (substitutionModels == null) {
            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);
}

Example 9 with DefaultBranchRateModel

use of dr.evomodel.branchratemodel.DefaultBranchRateModel 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);
    // 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);
        }
    }
    if (patternLists.size() == 0) {
        throw new XMLParseException("Either a single set of patterns should be given or multiple 'partitions' elements within DataTreeLikelihood: " + xo.getId());
    }
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (branchRateModel == null) {
        branchRateModel = new DefaultBranchRateModel();
    }
    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 " + "BeagleDataLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    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, scalingScheme, delayScaling);
}

Example 10 with DefaultBranchRateModel

use of dr.evomodel.branchratemodel.DefaultBranchRateModel in project beast-mcmc by beast-dev.

the class CompleteHistorySimulatorTest method testHKYSimulation.

public void testHKYSimulation() {
    Parameter kappa = new Parameter.Default(1, 2.0);
    double[] pi = { 0.45, 0.05, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    int stateCount = hky.getDataType().getStateCount();
    Parameter mu = new Parameter.Default(1, 0.5);
    Parameter alpha = new Parameter.Default(1, 0.5);
    GammaSiteRateModel siteModel = new GammaSiteRateModel("gammaModel", mu, alpha, 4, null);
    siteModel.setSubstitutionModel(hky);
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    double analyticResult = TreeUtils.getTreeLength(tree, tree.getRoot()) * mu.getParameterValue(0);
    int nSites = 200;
    double[] register1 = new double[stateCount * stateCount];
    double[] register2 = new double[stateCount * stateCount];
    // Count all jumps
    MarkovJumpsCore.fillRegistrationMatrix(register1, stateCount);
    // Move some jumps from 1 to 2
    register1[1 * stateCount + 2] = 0;
    register2[1 * stateCount + 2] = 1;
    register1[1 * stateCount + 3] = 0;
    register2[1 * stateCount + 3] = 1;
    register1[2 * stateCount + 3] = 0;
    register2[2 * stateCount + 3] = 1;
    runSimulation(tree, siteModel, branchRateModel, nSites, new double[][] { register1, register2 }, analyticResult);
}