use of dr.evomodel.siteratemodel.GammaSiteRateModel 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);*/
}
use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.
the class BeagleOperationParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
SitePatterns patternList = (SitePatterns) xo.getChild(PatternList.class);
TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
Alignment alignment = (Alignment) xo.getChild(Alignment.class);
GammaSiteRateModel substitutionModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
PrintWriter branch = null, operation = null;
if (xo.hasAttribute(BRANCH_FILE_NAME)) {
branch = XMLParser.getFilePrintWriter(xo, OPERATION_REPORT, BRANCH_FILE_NAME);
}
if (xo.hasAttribute(OPERATION_FILE_NAME)) {
operation = XMLParser.getFilePrintWriter(xo, OPERATION_REPORT, OPERATION_FILE_NAME);
}
return new BeagleOperationReport(treeModel, patternList, rateModel, substitutionModel, alignment, branch, operation);
}
use of dr.evomodel.siteratemodel.GammaSiteRateModel 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);
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;
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");
}
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);
}
use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.
the class BeagleTreeLikelihood method main.
public static void main(String[] args) {
try {
MathUtils.setSeed(666);
System.out.println("Test case 1: simulateOnePartition");
int sequenceLength = 1000;
ArrayList<Partition> partitionsList = new ArrayList<Partition>();
// create tree
NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
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);
Parameter kappa2 = new Parameter.Default(1, 1);
HKY hky1 = new HKY(kappa1, freqModel);
HKY hky2 = new HKY(kappa2, freqModel);
HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
substitutionModels.add(hky1);
substitutionModels.add(hky2);
List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
freqModels.add(freqModel);
Parameter epochTimes = new Parameter.Default(1, 20);
// create branch rate model
Parameter rate = new Parameter.Default(1, 0.001);
BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
// create site model
GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
// create partition
Partition partition1 = new //
Partition(//
treeModel, //
homogenousBranchSubstitutionModel, //
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);
BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
} catch (Exception e) {
e.printStackTrace();
System.exit(-1);
}
// END: try-catch block
}
use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.
the class ProductChainSubstitutionModelTest method setUpTwoStatesUnequalRate.
private void setUpTwoStatesUnequalRate() {
FrequencyModel freqModel0 = new FrequencyModel(TwoStates.INSTANCE, new double[] { 1.0 / 3.0, 2.0 / 3.0 });
FrequencyModel freqModel1 = new FrequencyModel(TwoStates.INSTANCE, new double[] { 1.0 / 4.0, 3.0 / 4.0 });
GeneralSubstitutionModel substModel0 = new GeneralSubstitutionModel("model0", TwoStates.INSTANCE, freqModel0, new Parameter.Default(new double[] { 1 }), 1);
GeneralSubstitutionModel substModel1 = new GeneralSubstitutionModel("model1", TwoStates.INSTANCE, freqModel1, new Parameter.Default(new double[] { 1 }), 1);
baseModels = new ArrayList<SubstitutionModel>();
baseModels.add(substModel0);
baseModels.add(substModel1);
SiteRateModel rateModel0 = new GammaSiteRateModel("rate0", new Parameter.Default(new double[] { 0.5 }), null, -1, null);
SiteRateModel rateModel1 = new GammaSiteRateModel("rate0", // Runs twice as fast
new Parameter.Default(new double[] { 2 }), null, -1, null);
List<SiteRateModel> rateModels = new ArrayList<SiteRateModel>();
rateModels.add(rateModel0);
rateModels.add(rateModel1);
productChainModel = new ProductChainSubstitutionModel("productChain", baseModels, rateModels);
stateCount = 2 * 2;
/*
model0 = as.eigen.two.state(0.5 * 1.5, 0.5 * 0.75) # rate = 0.5
model1 = as.eigen.two.state(2 * 2.0, 2 * 2.0 / 3.0) # rate = 2.0
pc = ind.two.eigen(model0, model1)
pc$rate.matrix
matexp(pc, 0.5)
*/
markovJumpsInfinitesimalResult = new double[] { -4.750000, 4.000000, 0.750000, 0.000000, 1.333333, -2.083333, 0.000000, 0.750000, 0.375000, 0.000000, -4.375000, 4.000000, 0.000000, 0.375000, 1.333333, -1.708333 };
markovJumpsEigenValues = new double[] { -6.458333, -5.333333, -1.125000, 0.000000 };
markovJumpsProbs = new double[] { 0.21546324, 0.4977253, 0.08664935, 0.2001621, 0.16590844, 0.5472801, 0.06672070, 0.2200907, 0.04332467, 0.1000811, 0.25878791, 0.5978064, 0.03336035, 0.1100454, 0.19926879, 0.6573255 };
}
Aggregations