Search in sources :

Example 51 with FrequencyModel

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);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) PatternList(dr.evolution.alignment.PatternList) HiddenDataType(dr.evolution.datatype.HiddenDataType) DataType(dr.evolution.datatype.DataType) Parameter(dr.inference.model.Parameter) HiddenDataType(dr.evolution.datatype.HiddenDataType) NumberFormat(java.text.NumberFormat)

Example 52 with FrequencyModel

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);
    }
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) GTR(dr.evomodel.substmodel.nucleotide.GTR) Parameter(dr.inference.model.Parameter)

Example 53 with FrequencyModel

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);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) GeneralizedLinearModel(dr.inference.distribution.GeneralizedLinearModel) DataType(dr.evolution.datatype.DataType) LogLinearModel(dr.inference.distribution.LogLinearModel) OldGLMSubstitutionModel(dr.evomodel.substmodel.OldGLMSubstitutionModel)

Example 54 with FrequencyModel

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);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) HKY(dr.evomodel.substmodel.nucleotide.HKY) Vector(dr.math.matrixAlgebra.Vector) MarkovJumpsSubstitutionModel(dr.evomodel.substmodel.MarkovJumpsSubstitutionModel)

Example 55 with FrequencyModel

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);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) HKY(dr.evomodel.substmodel.nucleotide.HKY) Vector(dr.math.matrixAlgebra.Vector) MarkovJumpsSubstitutionModel(dr.evomodel.substmodel.MarkovJumpsSubstitutionModel)

Aggregations

FrequencyModel (dr.evomodel.substmodel.FrequencyModel)57 Parameter (dr.inference.model.Parameter)42 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)23 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)22 HKY (dr.evomodel.substmodel.nucleotide.HKY)21 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)19 TreeModel (dr.evomodel.tree.TreeModel)19 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)17 ArrayList (java.util.ArrayList)14 BranchModel (dr.evomodel.branchmodel.BranchModel)12 Partition (dr.app.beagle.tools.Partition)11 SubstitutionModel (dr.evomodel.substmodel.SubstitutionModel)11 BeagleSequenceSimulator (dr.app.beagle.tools.BeagleSequenceSimulator)10 DataType (dr.evolution.datatype.DataType)10 NewickImporter (dr.evolution.io.NewickImporter)9 Tree (dr.evolution.tree.Tree)9 Vector (dr.math.matrixAlgebra.Vector)9 PatternList (dr.evolution.alignment.PatternList)7 ImportException (dr.evolution.io.Importer.ImportException)7 BeagleTreeLikelihood (dr.evomodel.treelikelihood.BeagleTreeLikelihood)7