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
}
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;
}
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);
}
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);
}
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);
}
Aggregations