Examples with GammaSiteModel dr.oldevomodel.sitemodel.GammaSiteModel used on opensource projects

Example 16 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class GammaSiteModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    SubstitutionModel substitutionModel = (SubstitutionModel) xo.getElementFirstChild(SUBSTITUTION_MODEL);
    String msg = "";
    Parameter muParam = null;
    if (xo.hasChildNamed(SUBSTITUTION_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(SUBSTITUTION_RATE);
        msg += "\n  with initial substitution rate = " + muParam.getParameterValue(0);
    } else if (xo.hasChildNamed(MUTATION_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(MUTATION_RATE);
        msg += "\n  with initial substitution rate = " + muParam.getParameterValue(0);
    } else if (xo.hasChildNamed(RELATIVE_RATE)) {
        muParam = (Parameter) xo.getElementFirstChild(RELATIVE_RATE);
        msg += "\n  with initial relative rate = " + muParam.getParameterValue(0);
    }
    Parameter shapeParam = null;
    int catCount = 4;
    if (xo.hasChildNamed(GAMMA_SHAPE)) {
        final XMLObject cxo = xo.getChild(GAMMA_SHAPE);
        catCount = cxo.getIntegerAttribute(GAMMA_CATEGORIES);
        shapeParam = (Parameter) cxo.getChild(Parameter.class);
        msg += "\n  " + catCount + " category discrete gamma with initial shape = " + shapeParam.getParameterValue(0);
    }
    Parameter invarParam = null;
    if (xo.hasChildNamed(PROPORTION_INVARIANT)) {
        invarParam = (Parameter) xo.getElementFirstChild(PROPORTION_INVARIANT);
        msg += "\n  initial proportion of invariant sites = " + invarParam.getParameterValue(0);
    }
    Logger.getLogger("dr.evomodel").info("Creating site model." + (msg.length() > 0 ? msg : ""));
    return new GammaSiteModel(substitutionModel, muParam, shapeParam, catCount, invarParam);
}

Example 17 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class MsatFullLikelihoodTest method setUpExample3.

private void setUpExample3() throws Exception {
    //taxa
    ArrayList<Taxon> taxonList3 = new ArrayList<Taxon>();
    Collections.addAll(taxonList3, new Taxon("taxon1"), new Taxon("taxon2"), new Taxon("taxon3"), new Taxon("taxon4"), new Taxon("taxon5"), new Taxon("taxon6"), new Taxon("taxon7"));
    Taxa taxa3 = new Taxa(taxonList3);
    //msat datatype
    Microsatellite msat = new Microsatellite(1, 4);
    Patterns msatPatterns = new Patterns(msat, taxa3);
    //pattern in the correct code form.
    msatPatterns.addPattern(new int[] { 0, 3, 1, 2, 3, 0, 1 });
    //create tree
    NewickImporter importer = new NewickImporter("(((taxon1:0.3,taxon2:0.3):0.6,taxon3:0.9):0.9,((taxon4:0.5,taxon5:0.5):0.3,(taxon6:0.7,taxon7:0.7):0.1):1.0);");
    Tree tree = importer.importTree(null);
    //treeModel
    TreeModel treeModel = new TreeModel(tree);
    //msatsubstModel
    AsymmetricQuadraticModel aqm3 = new AsymmetricQuadraticModel(msat, null);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(aqm3);
    //treeLikelihood
    treeLikelihood3 = new TreeLikelihood(msatPatterns, treeModel, siteModel, null, null, false, false, true, false, false);
}

Example 18 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class SequenceLikelihoodTest method computeSitePatternLikelihoods.

protected double[] computeSitePatternLikelihoods(SitePatterns patterns) {
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 29.739445, 0, 100);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    //treeLikelihood
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, null, null, false, false, true, false, false);
    return treeLikelihood.getPatternLogLikelihoods();
}

Example 19 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class MCMCTest method testMCMC.

public void testMCMC() {
    // Sub model
    //new double[]{0.25, 0.25, 0.25, 0.25});
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    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(kappa, 0.5);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    //        Parameter rootParameter = treeModel.createNodeHeightsParameter(true, false, false);
    //        ScaleOperator scaleOperator = new ScaleOperator(rootParameter, 0.75, CoercionMode.COERCION_ON, 1.0);
    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(kappa);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(kappa);
    // 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("MCMCTest", 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="hky.kappa" value="32.8941"/>
    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, 6.42048E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 32.8941);
}

Example 20 with GammaSiteModel

use of dr.oldevomodel.sitemodel.GammaSiteModel in project beast-mcmc by beast-dev.

the class PMDTestProblem method testPMD.

public void testPMD() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 496432.69917113904, 0, Double.POSITIVE_INFINITY);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter rateParameter = new Parameter.Default(StrictClockBranchRates.RATE, 4.0E-7, 0, 100.0);
    StrictClockBranchRates branchRateModel = new StrictClockBranchRates(rateParameter);
    // Sub model
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 1.0E-8, Double.POSITIVE_INFINITY);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    // SequenceErrorModel
    Parameter ageRelatedRateParameter = new Parameter.Default(SequenceErrorModelParser.AGE_RELATED_RATE, 4.0E-7, 0, 100.0);
    TipStatesModel aDNADamageModel = new SequenceErrorModel(null, null, SequenceErrorModel.ErrorType.TRANSITIONS_ONLY, null, ageRelatedRateParameter, null);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, aDNADamageModel, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(rateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter allInternalHeights = treeModel.createNodeHeightsParameter(true, true, false);
    operator = new UpDownOperator(new Scalable[] { new Scalable.Default(rateParameter) }, new Scalable[] { new Scalable.Default(allInternalHeights) }, 0.75, 3.0, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(ageRelatedRateParameter, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 15.0, 49643.2699171139, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    // ??? correct?
    operator = new DeltaExchangeOperator(freqs, new int[] { 1, 1, 1, 1 }, 0.01, 1.0, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    //CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(rateParameter);
    loggers[0].add(ageRelatedRateParameter);
    loggers[0].add(popSize);
    loggers[0].add(kappa);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(rateParameter);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("PMDTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="clock.rate" value="1.5E-7"/>
    //        <expectation name="errorModel.ageRate" value="0.7E-7"/>
    //        <expectation name="hky.kappa" value="10"/>
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 10);
    TraceCorrelation rateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(StrictClockBranchRates.RATE));
    assertExpectation(StrictClockBranchRates.RATE, rateStats, 1.5E-7);
    TraceCorrelation ageRateStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(SequenceErrorModelParser.AGE_RELATED_RATE));
    assertExpectation(SequenceErrorModelParser.AGE_RELATED_RATE, ageRateStats, 0.7E-7);
}