Examples with FrequencyModel dr.oldevomodel.substmodel.FrequencyModel used on opensource projects

Example 6 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class EmpiricalCodonModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Codons codons = Codons.UNIVERSAL;
    if (xo.hasAttribute(GeneticCode.GENETIC_CODE)) {
        String codeStr = xo.getStringAttribute(GeneticCode.GENETIC_CODE);
        if (codeStr.equals(GeneticCode.UNIVERSAL.getName())) {
            codons = Codons.UNIVERSAL;
        } else if (codeStr.equals(GeneticCode.VERTEBRATE_MT.getName())) {
            codons = Codons.VERTEBRATE_MT;
        } else if (codeStr.equals(GeneticCode.YEAST.getName())) {
            codons = Codons.YEAST;
        } else if (codeStr.equals(GeneticCode.MOLD_PROTOZOAN_MT.getName())) {
            codons = Codons.MOLD_PROTOZOAN_MT;
        } else if (codeStr.equals(GeneticCode.INVERTEBRATE_MT.getName())) {
            codons = Codons.INVERTEBRATE_MT;
        } else if (codeStr.equals(GeneticCode.CILIATE.getName())) {
            codons = Codons.CILIATE;
        } else if (codeStr.equals(GeneticCode.ECHINODERM_MT.getName())) {
            codons = Codons.ECHINODERM_MT;
        } else if (codeStr.equals(GeneticCode.EUPLOTID_NUC.getName())) {
            codons = Codons.EUPLOTID_NUC;
        } else if (codeStr.equals(GeneticCode.BACTERIAL.getName())) {
            codons = Codons.BACTERIAL;
        } else if (codeStr.equals(GeneticCode.ALT_YEAST.getName())) {
            codons = Codons.ALT_YEAST;
        } else if (codeStr.equals(GeneticCode.ASCIDIAN_MT.getName())) {
            codons = Codons.ASCIDIAN_MT;
        } else if (codeStr.equals(GeneticCode.FLATWORM_MT.getName())) {
            codons = Codons.FLATWORM_MT;
        } else if (codeStr.equals(GeneticCode.BLEPHARISMA_NUC.getName())) {
            codons = Codons.BLEPHARISMA_NUC;
        } else if (codeStr.equals(GeneticCode.NO_STOPS.getName())) {
            codons = Codons.NO_STOPS;
        }
    }
    Parameter omegaParam = (Parameter) xo.getElementFirstChild(OMEGA);
    Parameter kappaParam = null;
    Parameter mntParam = null;
    if (xo.hasChildNamed(KAPPATSTV)) {
        kappaParam = (Parameter) xo.getElementFirstChild(KAPPATSTV);
        if (kappaParam.getDimension() != 2 && kappaParam.getDimension() != 9) {
            throw new XMLParseException("If you use the kappa parameter, you need to enter exactly\n" + "two values for ts and tv or nine values\n" + "according to the Kosiol ECM+F+omega+9k model");
        }
    } else {
        mntParam = (Parameter) xo.getElementFirstChild(MULTI_NT_CHANGE);
    }
    String dirString = xo.getStringAttribute(ECM_DATA_DIR);
    String freqString = xo.getStringAttribute(ECM_FREQ_MATRIX);
    String matString = xo.getStringAttribute(ECM_DATA_MATRIX);
    EmpiricalCodonRateMatrix rateMat = new EmpiricalCodonRateMatrix(EMPIRICAL_RATE_MATRIX, codons, dirString, freqString, matString);
    // get frequencies from XML, from frequency csv file or estimate from data
    FrequencyModel freqModel = null;
    if (xo.getChild(FrequencyModel.class) != null) {
        freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    } else {
        freqModel = createNewFreqModel(codons, rateMat);
    }
    return new EmpiricalCodonModel(codons, omegaParam, kappaParam, mntParam, rateMat, freqModel);
}

Example 7 with FrequencyModel

use of dr.oldevomodel.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) {
            frequencies = ((PatternList) obj).getStateFrequencies();
            break;
        }
    }
    StringBuilder sb = new StringBuilder("Creating 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("= {");
    double sum = 0;
    for (int j = 0; j < freqsParam.getDimension(); j++) {
        sum += freqsParam.getParameterValue(j);
    }
    if (xo.getAttribute(NORMALIZE, false)) {
        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());
        }
        sum = 1.0;
    }
    if (Math.abs(sum - 1.0) > 1e-8) {
        throw new XMLParseException("Frequencies do not sum to 1 (they sum to " + sum + ")");
    }
    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);
}

Example 8 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class AsymQuadModelParser method parseXMLObject.

//AbstractXMLObjectParser implementation
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Microsatellite microsatellite = (Microsatellite) xo.getChild(Microsatellite.class);
    Parameter expanConst = processModelParameter(xo, EXPANSION_CONSTANT);
    Parameter expanLin = processModelParameter(xo, EXPANSION_LIN);
    Parameter expanQuad = processModelParameter(xo, EXPANSION_QUAD);
    Parameter contractConst = processModelParameter(xo, CONTRACTION_CONSTANT);
    Parameter contractLin = processModelParameter(xo, CONTRACTION_LIN);
    Parameter contractQuad = processModelParameter(xo, CONTRACTION_QUAD);
    //get FrequencyModel
    FrequencyModel freqModel = null;
    if (xo.hasChildNamed(FrequencyModelParser.FREQUENCIES)) {
        freqModel = (FrequencyModel) xo.getElementFirstChild(FrequencyModelParser.FREQUENCIES);
    }
    boolean isSubmodel = xo.getAttribute(IS_SUBMODEL, false);
    return new AsymmetricQuadraticModel(microsatellite, freqModel, expanConst, expanLin, expanQuad, contractConst, contractLin, contractQuad, isSubmodel);
}

Example 9 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel 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(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(dataType, freqModel, ratesParameter, relativeTo);
}

Example 10 with FrequencyModel

use of dr.oldevomodel.substmodel.FrequencyModel in project beast-mcmc by beast-dev.

the class GeneralSubstitutionModelTest method testGeneralSubstitutionModel.

public void testGeneralSubstitutionModel() {
    // Sub model
    FrequencyModel freqModel = new FrequencyModel(dataType, alignment.getStateFrequencies());
    // dimension="5" value="1.0"
    Parameter ratesPara = new Parameter.Default(GeneralSubstitutionModelParser.RATES, 5, 1.0);
    // relativeTo="5"
    GeneralSubstitutionModel generalSubstitutionModel = new GeneralSubstitutionModel(dataType, freqModel, ratesPara, 4);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(generalSubstitutionModel);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(ratesPara, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 10.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 1, 1, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 1.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(ratesPara);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(ratesPara);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(10000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, treeLikelihood, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("GeneralSubstitutionModelTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="likelihood" value="-1815.75"/>
    //		<expectation name="treeModel.rootHeight" value="6.42048E-2"/>
    //		<expectation name="rateAC" value="6.08986E-2"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.75);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 0.0640787258170083);
    TraceCorrelation rateACStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(GeneralSubstitutionModelParser.RATES + "1"));
    assertExpectation(GeneralSubstitutionModelParser.RATES + "1", rateACStats, 0.061071756742081366);
}