use of dr.evomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class FrequencyModelParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
DataType dataType = DataTypeUtils.getDataType(xo);
Parameter freqsParam = (Parameter) xo.getElementFirstChild(FREQUENCIES);
double[] frequencies = null;
for (int i = 0; i < xo.getChildCount(); i++) {
Object obj = xo.getChild(i);
if (obj instanceof PatternList) {
PatternList patternList = (PatternList) obj;
if (xo.getAttribute(COMPRESS, false) && (patternList.getDataType() instanceof HiddenDataType)) {
double[] hiddenFrequencies = patternList.getStateFrequencies();
int hiddenCount = ((HiddenDataType) patternList.getDataType()).getHiddenClassCount();
int baseStateCount = hiddenFrequencies.length / hiddenCount;
frequencies = new double[baseStateCount];
for (int j = 0; j < baseStateCount; ++j) {
for (int k = 0; k < hiddenCount; ++k) {
frequencies[j] += hiddenFrequencies[j + k * baseStateCount];
}
}
} else {
// TODO
if (xo.hasAttribute(COMPOSITION)) {
String type = xo.getStringAttribute(COMPOSITION);
if (type.equalsIgnoreCase(FREQ_3x4)) {
frequencies = getEmpirical3x4Freqs(patternList);
}
} else {
frequencies = patternList.getStateFrequencies();
}
// END: composition check
}
break;
}
// END: patternList check
}
StringBuilder sb = new StringBuilder("\nCreating state frequencies model '" + freqsParam.getParameterName() + "': ");
if (frequencies != null) {
if (freqsParam.getDimension() != frequencies.length) {
throw new XMLParseException("dimension of frequency parameter and number of sequence states don't match.");
}
for (int j = 0; j < frequencies.length; j++) {
freqsParam.setParameterValue(j, frequencies[j]);
}
sb.append("Using empirical frequencies from data ");
} else {
sb.append("Initial frequencies ");
}
sb.append("= {");
if (xo.getAttribute(NORMALIZE, false)) {
double sum = 0;
for (int j = 0; j < freqsParam.getDimension(); j++) sum += freqsParam.getParameterValue(j);
for (int j = 0; j < freqsParam.getDimension(); j++) {
if (sum != 0)
freqsParam.setParameterValue(j, freqsParam.getParameterValue(j) / sum);
else
freqsParam.setParameterValue(j, 1.0 / freqsParam.getDimension());
}
}
NumberFormat format = NumberFormat.getNumberInstance();
format.setMaximumFractionDigits(5);
sb.append(format.format(freqsParam.getParameterValue(0)));
for (int j = 1; j < freqsParam.getDimension(); j++) {
sb.append(", ");
sb.append(format.format(freqsParam.getParameterValue(j)));
}
sb.append("}");
Logger.getLogger("dr.evomodel").info(sb.toString());
return new FrequencyModel(dataType, freqsParam);
}
use of dr.evomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class GTRParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
FrequencyModel freqModel = (FrequencyModel) xo.getElementFirstChild(FREQUENCIES);
Parameter rates = null;
if (xo.hasChildNamed(RATES)) {
rates = (Parameter) xo.getElementFirstChild(RATES);
rates.setDimensionNames(new String[] { rates.getId() + A_TO_C, rates.getId() + A_TO_G, rates.getId() + A_TO_T, rates.getId() + C_TO_G, rates.getId() + C_TO_T, rates.getId() + G_TO_T });
return new GTR(rates, freqModel);
} else {
Variable<Double> rateACVariable = null;
if (xo.hasChildNamed(A_TO_C)) {
rateACVariable = (Variable<Double>) xo.getElementFirstChild(A_TO_C);
}
Variable<Double> rateAGVariable = null;
if (xo.hasChildNamed(A_TO_G)) {
rateAGVariable = (Variable<Double>) xo.getElementFirstChild(A_TO_G);
}
Variable<Double> rateATVariable = null;
if (xo.hasChildNamed(A_TO_T)) {
rateATVariable = (Variable<Double>) xo.getElementFirstChild(A_TO_T);
}
Variable<Double> rateCGVariable = null;
if (xo.hasChildNamed(C_TO_G)) {
rateCGVariable = (Variable<Double>) xo.getElementFirstChild(C_TO_G);
}
Variable<Double> rateCTVariable = null;
if (xo.hasChildNamed(C_TO_T)) {
rateCTVariable = (Variable<Double>) xo.getElementFirstChild(C_TO_T);
}
Variable<Double> rateGTVariable = null;
if (xo.hasChildNamed(G_TO_T)) {
rateGTVariable = (Variable<Double>) xo.getElementFirstChild(G_TO_T);
}
int countNull = 0;
if (rateACVariable == null)
countNull++;
if (rateAGVariable == null)
countNull++;
if (rateATVariable == null)
countNull++;
if (rateCGVariable == null)
countNull++;
if (rateCTVariable == null)
countNull++;
if (rateGTVariable == null)
countNull++;
if (countNull != 1)
throw new XMLParseException("Only five parameters may be specified in GTR, leave exactly one out, the others will be specifed relative to the one left out.");
return new GTR(rateACVariable, rateAGVariable, rateATVariable, rateCGVariable, rateCTVariable, rateGTVariable, freqModel);
}
}
use of dr.evomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class OldGLMSubstitutionModelParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
DataType dataType = DataTypeUtils.getDataType(xo);
if (dataType == null)
dataType = (DataType) xo.getChild(DataType.class);
int rateCount = (dataType.getStateCount() - 1) * dataType.getStateCount();
LogLinearModel glm = (LogLinearModel) xo.getChild(GeneralizedLinearModel.class);
int length = glm.getXBeta().length;
if (length != rateCount) {
throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (rateCount) + " dimensions. However GLM dimension is " + length);
}
XMLObject cxo = xo.getChild(dr.oldevomodelxml.substmodel.ComplexSubstitutionModelParser.ROOT_FREQUENCIES);
FrequencyModel rootFreq = (FrequencyModel) cxo.getChild(FrequencyModel.class);
if (dataType != rootFreq.getDataType()) {
throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its rootFrequencyModel.");
}
return new OldGLMSubstitutionModel(xo.getId(), dataType, rootFreq, glm);
}
use of dr.evomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class MarkovJumpsSubstitutionModelTest method testMarkovJumpsCounts.
public void testMarkovJumpsCounts() {
HKY substModel = new HKY(2.0, new FrequencyModel(Nucleotides.INSTANCE, // A,C,G,T
new double[] { 0.3, 0.2, 0.25, 0.25 }));
int states = substModel.getDataType().getStateCount();
MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(substModel, MarkovJumpsType.COUNTS);
double[] r = new double[states * states];
double[] q = new double[states * states];
double[] j = new double[states * states];
double[] c = new double[states * states];
double[] p = new double[states * states];
double time = 1.0;
// A
int from = 0;
// C
int to = 1;
MarkovJumpsCore.fillRegistrationMatrix(r, from, to, states, 1.0);
markovjumps.setRegistration(r);
substModel.getInfinitesimalMatrix(q);
substModel.getTransitionProbabilities(time, p);
markovjumps.computeJointStatMarkovJumps(time, j);
markovjumps.computeCondStatMarkovJumps(time, c);
MarkovJumpsCore.makeComparableToRPackage(q);
System.out.println("Q = " + new Vector(q));
MarkovJumpsCore.makeComparableToRPackage(p);
System.out.println("P = " + new Vector(p));
System.out.println("Counts:");
MarkovJumpsCore.makeComparableToRPackage(r);
System.out.println("R = " + new Vector(r));
MarkovJumpsCore.makeComparableToRPackage(j);
System.out.println("J = " + new Vector(j));
assertEquals(rMarkovJumpsJ, j, tolerance);
MarkovJumpsCore.makeComparableToRPackage(c);
System.err.println("C = " + new Vector(c));
assertEquals(rMarkovJumpsC, c, tolerance);
}
use of dr.evomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.
the class MarkovJumpsSubstitutionModelTest method testMarkovJumpsReward.
public void testMarkovJumpsReward() {
HKY substModel = new HKY(2.0, new FrequencyModel(Nucleotides.INSTANCE, // A,C,G,T
new double[] { 0.3, 0.2, 0.25, 0.25 }));
int states = substModel.getDataType().getStateCount();
double[] j = new double[states * states];
double[] c = new double[states * states];
double time = 1.0;
MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(substModel, MarkovJumpsType.REWARDS);
double[] rewards = { 1.0, 1.0, 1.0, 1.0 };
markovjumps.setRegistration(rewards);
markovjumps.computeJointStatMarkovJumps(time, j);
markovjumps.computeCondStatMarkovJumps(time, c);
System.out.println("Rewards:");
MarkovJumpsCore.makeComparableToRPackage(rewards);
System.out.println("R = " + new Vector(rewards));
MarkovJumpsCore.makeComparableToRPackage(j);
System.out.println("J = " + new Vector(j));
assertEquals(rMarkovRewardsJ, j, tolerance);
MarkovJumpsCore.makeComparableToRPackage(c);
System.out.println("C = " + new Vector(c));
assertEquals(rMarkovRewardsC, c, tolerance);
}
Aggregations