Examples with GeneralSubstitutionModel dr.evomodel.substmodel.GeneralSubstitutionModel used on opensource projects

Example 1 with GeneralSubstitutionModel

use of dr.evomodel.substmodel.GeneralSubstitutionModel in project beast-mcmc by beast-dev.

the class GeneralSubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Parameter ratesParameter = null;
    FrequencyModel freqModel = null;
    if (xo.hasChildNamed(FREQUENCIES)) {
        XMLObject cxo = xo.getChild(FREQUENCIES);
        freqModel = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    }
    DataType dataType = DataTypeUtils.getDataType(xo);
    if (dataType == null)
        dataType = (DataType) xo.getChild(DataType.class);
    if (dataType == null)
        dataType = freqModel.getDataType();
    if (dataType != freqModel.getDataType()) {
        throw new XMLParseException("Data type of " + getParserName() + " element does not match that of its frequencyModel.");
    }
    XMLObject cxo = xo.getChild(RATES);
    ratesParameter = (Parameter) cxo.getChild(Parameter.class);
    int states = dataType.getStateCount();
    Logger.getLogger("dr.evomodel").info("  General Substitution Model (stateCount=" + states + ")");
    boolean hasRelativeRates = cxo.hasChildNamed(RELATIVE_TO) || (cxo.hasAttribute(RELATIVE_TO) && cxo.getIntegerAttribute(RELATIVE_TO) > 0);
    int nonReversibleRateCount = ((dataType.getStateCount() - 1) * dataType.getStateCount());
    int reversibleRateCount = (nonReversibleRateCount / 2);
    boolean isNonReversible = ratesParameter.getDimension() == nonReversibleRateCount;
    boolean hasIndicator = xo.hasChildNamed(INDICATOR);
    if (!hasRelativeRates) {
        Parameter indicatorParameter = null;
        if (ratesParameter.getDimension() != reversibleRateCount && ratesParameter.getDimension() != nonReversibleRateCount) {
            throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (reversibleRateCount) + " dimensions for reversible model or " + nonReversibleRateCount + " dimensions for non-reversible. " + "However parameter dimension is " + ratesParameter.getDimension());
        }
        if (hasIndicator) {
            // this is using BSSVS
            cxo = xo.getChild(INDICATOR);
            indicatorParameter = (Parameter) cxo.getChild(Parameter.class);
            if (indicatorParameter.getDimension() != ratesParameter.getDimension()) {
                throw new XMLParseException("Rates and indicator parameters in " + getParserName() + " element must be the same dimension.");
            }
            boolean randomize = xo.getAttribute(ComplexSubstitutionModelParser.RANDOMIZE, false);
            if (randomize) {
                BayesianStochasticSearchVariableSelection.Utils.randomize(indicatorParameter, dataType.getStateCount(), !isNonReversible);
            }
        }
        if (isNonReversible) {
            //                if (xo.hasChildNamed(ROOT_FREQ)) {
            //                    cxo = xo.getChild(ROOT_FREQ);
            //                    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.");
            //                    }
            //
            //                    Logger.getLogger("dr.evomodel").info("  Using BSSVS Complex Substitution Model");
            //                    return new SVSComplexSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, indicatorParameter);
            //
            //                } else {
            //                    throw new XMLParseException("Non-reversible model missing " + ROOT_FREQ + " element");
            //                }
            Logger.getLogger("dr.evomodel").info("  Using BSSVS Complex Substitution Model");
            return new SVSComplexSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, indicatorParameter);
        } else {
            Logger.getLogger("dr.evomodel").info("  Using BSSVS General Substitution Model");
            return new SVSGeneralSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, indicatorParameter);
        }
    } else {
        if (ratesParameter.getDimension() != reversibleRateCount - 1) {
            throw new XMLParseException("Rates parameter in " + getParserName() + " element should have " + (reversibleRateCount - 1) + " dimensions. However parameter dimension is " + ratesParameter.getDimension());
        }
        int relativeTo = 0;
        if (hasRelativeRates) {
            relativeTo = cxo.getIntegerAttribute(RELATIVE_TO) - 1;
        }
        if (relativeTo < 0 || relativeTo >= reversibleRateCount) {
            throw new XMLParseException(RELATIVE_TO + " must be 1 or greater");
        } else {
            int t = relativeTo;
            int s = states - 1;
            int row = 0;
            while (t >= s) {
                t -= s;
                s -= 1;
                row += 1;
            }
            int col = t + row + 1;
            Logger.getLogger("dr.evomodel").info("  Rates relative to " + dataType.getCode(row) + "<->" + dataType.getCode(col));
        }
        if (ratesParameter == null) {
            if (reversibleRateCount == 1) {
            // simplest model for binary traits...
            } else {
                throw new XMLParseException("No rates parameter found in " + getParserName());
            }
        }
        return new GeneralSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, relativeTo);
    }
}

Example 2 with GeneralSubstitutionModel

use of dr.evomodel.substmodel.GeneralSubstitutionModel in project beast-mcmc by beast-dev.

the class HiddenMarkovRatesTest method testHiddenRates.

public void testHiddenRates() {
    final int index = 0;
    double[] freqs = getBinaryFreqs(index);
    FrequencyModel binaryFreqModel = new FrequencyModel(TwoStates.INSTANCE, freqs);
    int relativeTo = 0;
    Parameter ratesParameter = new Parameter.Default(0);
    GeneralSubstitutionModel binaryModel = new GeneralSubstitutionModel("binary", TwoStates.INSTANCE, binaryFreqModel, ratesParameter, relativeTo);
    UniformizedSubstitutionModel uSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS);
    uSM.setSaveCompleteHistory(true);
    double[] rewardRegister = new double[] { 0.0, 1.0 };
    uSM.setRegistration(rewardRegister);
    final double[] hkyFreqs = getHKYFreqs(index);
    FrequencyModel hkyFreqModel = new FrequencyModel(Nucleotides.INSTANCE, hkyFreqs);
    final double kappa = getKappa(index);
    final HKY hky = new HKY(kappa, hkyFreqModel);
    final double length = getLength(index);
    double[] resultCompleteHistory = new double[16];
    final int replicates = getNumberReplicates(index);
    double result = 0.0;
    for (int r = 0; r < replicates; ++r) {
        result += oneCompleteHistoryReplicate(resultCompleteHistory, hky, uSM, length);
    }
    result /= replicates;
    normalize(resultCompleteHistory, replicates);
    System.out.println("Averaged probabilities");
    System.out.println(result);
    System.out.println(new Vector(resultCompleteHistory));
    System.out.println();
    double[] intermediate = new double[16];
    hky.getTransitionProbabilities(result, intermediate);
    System.out.println("Intermediate using above average reward");
    System.out.println(result);
    System.out.println(new Vector(intermediate));
    System.out.println();
    double[] resultExpected = new double[16];
    UniformizedSubstitutionModel expectedUSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS, replicates);
    expectedUSM.setRegistration(rewardRegister);
    result = oneCompleteHistoryReplicate(resultExpected, hky, expectedUSM, length);
    System.out.println("Averaged reward");
    System.out.println(result);
    System.out.println(new Vector(resultExpected));
    System.out.println();
    double[] originalProbs = new double[16];
    hky.getTransitionProbabilities(length, originalProbs);
    System.out.println("Original probabilities");
    System.out.println(new Vector(originalProbs));
    System.out.println();
}

Example 3 with GeneralSubstitutionModel

use of dr.evomodel.substmodel.GeneralSubstitutionModel in project beast-mcmc by beast-dev.

the class BinarySubstitutionModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Parameter ratesParameter;
    XMLObject cxo = xo.getChild(dr.oldevomodelxml.substmodel.GeneralSubstitutionModelParser.FREQUENCIES);
    FrequencyModel freqModel = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    DataType dataType = freqModel.getDataType();
    if (dataType != TwoStates.INSTANCE)
        throw new XMLParseException("Frequency model must have binary (two state) data type.");
    int relativeTo = 0;
    ratesParameter = new Parameter.Default(0);
    return new GeneralSubstitutionModel(getParserName(), dataType, freqModel, ratesParameter, relativeTo);
}

Example 4 with GeneralSubstitutionModel

use of dr.evomodel.substmodel.GeneralSubstitutionModel in project beast-mcmc by beast-dev.

the class LewisMkSubstitutionModelParser method parseXMLObject.

//public static XMLObjectParser PARSER=new LewisMkSubstitutionModelParser();
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(FREQUENCIES);
    FrequencyModel freqModel = (FrequencyModel) cxo.getChild(FrequencyModel.class);
    DataType dataType = freqModel.getDataType();
    int k = dataType.getStateCount();
    System.err.println("Number of states " + k);
    Parameter ratesParameter;
    if (xo.hasAttribute(TOTAL_ORDER) && xo.getBooleanAttribute(TOTAL_ORDER)) {
        //TOTAL ORDERING OF THE STATES BASED ON DATATYPE
        ratesParameter = new Parameter.Default(k * (k - 1) / 2, 0);
        int j = k - 1;
        for (int i = 0; i < (k - 1) * k / 2; i = i + j + 1) {
            ratesParameter.setParameterValue(i, 1);
            j -= 1;
        }
    } else if (xo.hasChildNamed(ORDER)) {
        // USER-SPECIFIED ORDERING OF THE STATES
        ratesParameter = new Parameter.Default(k * (k - 1) / 2, 0);
        for (int i = 0; i < xo.getChildCount(); ++i) {
            if (xo.getChildName(i).equals(ORDER)) {
                cxo = (XMLObject) xo.getChild(i);
                if (cxo.getName().equals(ORDER)) {
                    int from = dataType.getState(cxo.getStringAttribute(STATE).charAt(0));
                    int to = dataType.getState(cxo.getStringAttribute(ADJACENT).charAt(0));
                    if (from > to) {
                        //SWAP: from should have the smaller state number
                        to += from;
                        from = to - from;
                        to -= from;
                    }
                    int ratesIndex = (from * (2 * k - 3) - from * from) / 2 + to - 1;
                    ratesParameter.setParameterValue(ratesIndex, 1);
                }
            }
        }
    } else {
        ratesParameter = new Parameter.Default(k * (k - 1) / 2, 1);
    }
    System.err.println(ratesParameter.toString());
    System.err.println("Infinitesimal matrix:");
    for (int i = 0; i < k; ++i) {
        for (int j = 0; j < k; ++j) {
            int from, to;
            if (i < j) {
                from = i;
                to = j;
            } else {
                from = j;
                to = i;
            }
            //This is right now!!! Thanks, Marc!
            int ratesIndex = (from * (2 * k - 3) - from * from) / 2 + to - 1;
            if (i != j)
                System.err.print(Double.toString(ratesParameter.getValue(ratesIndex)) + "\t(" + ratesIndex + ")\t");
            else
                System.err.print("-\t\t");
        }
        //newline
        System.err.println("");
    }
    System.err.println("");
    if (!checkConnected(ratesParameter.getValues(), k)) {
        throw (new XMLParseException("The state transitions form a disconnected graph! This model is not suited for this case."));
    }
    return new GeneralSubstitutionModel(LEWIS_MK_MODEL, dataType, freqModel, ratesParameter, -1);
}