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

Example 11 with GammaSiteModel

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

the class Coevolve method runModel2.

private void runModel2(PatternList patternList, PrintWriter pw, Tree tree, SubstitutionModel substModel, final Parameter betaParameter) {
    final Parameter muParameter = new Parameter.Default(1.0);
    muParameter.setId("mu");
    SiteModel siteModel = new GammaSiteModel(substModel, muParameter, null, 1, null);
    TreeModel treeModel = new TreeModel(tree);
    final TreeLikelihood treeLikelihood = new TreeLikelihood(patternList, treeModel, siteModel, null, null, false, false, true, false, false);
    treeLikelihood.setId("likelihood");
    MultivariateFunction function = new MultivariateFunction() {

        public double evaluate(double[] argument) {
            betaParameter.setParameterValue(0, argument[0]);
            betaParameter.setParameterValue(1, argument[1]);
            muParameter.setParameterValue(0, argument[2]);
            double lnL = -treeLikelihood.getLogLikelihood();
            // System.err.println("" + argument[0] + "\t" + argument[1] + "\t" + argument[2] + "\t" + lnL);
            return lnL;
        }

        public int getNumArguments() {
            return 3;
        }

        public double getLowerBound(int n) {
            return 0.0;
        }

        public double getUpperBound(int n) {
            return 100.0;
        }
    };
    MultivariateMinimum optimizer = new ConjugateGradientSearch();
    double lnL = optimizer.findMinimum(function, new double[] { 1.0, 1.0, 1.0 }, 6, 6);
    pw.write(betaParameter.getParameterValue(0) + "\t");
    pw.write(betaParameter.getParameterValue(1) + "\t");
    pw.write(muParameter.getParameterValue(0) + "\t");
    pw.write(lnL + "\n");
    pw.flush();
    System.out.println("" + betaParameter.getParameterValue(0) + "\t" + betaParameter.getParameterValue(1) + "\t" + muParameter.getParameterValue(0) + "\t" + lnL);
}

Example 12 with GammaSiteModel

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

the class RandomLocalClockTestProblem method testRandomLocalClock.

public void testRandomLocalClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 0.077, 0, Double.POSITIVE_INFINITY);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter ratesParameter = new Parameter.Default(RandomLocalClockModelParser.RATES, 10, 1);
    Parameter rateIndicatorParameter = new Parameter.Default(RandomLocalClockModelParser.RATE_INDICATORS, 10, 1);
    Parameter meanRateParameter = new Parameter.Default(RandomLocalClockModelParser.CLOCK_RATE, 1, 1.0);
    RandomLocalClockModel branchRateModel = new RandomLocalClockModel(treeModel, meanRateParameter, rateIndicatorParameter, ratesParameter, false, 0.5);
    SumStatistic rateChanges = new SumStatistic("rateChangeCount", true);
    rateChanges.addStatistic(rateIndicatorParameter);
    RateStatistic meanRate = new RateStatistic("meanRate", treeModel, branchRateModel, true, true, RateStatisticParser.MEAN);
    RateStatistic coefficientOfVariation = new RateStatistic(RateStatisticParser.COEFFICIENT_OF_VARIATION, treeModel, branchRateModel, true, true, RateStatisticParser.COEFFICIENT_OF_VARIATION);
    RateCovarianceStatistic covariance = new RateCovarianceStatistic("covariance", treeModel, branchRateModel);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 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, branchRateModel, null, 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(ratesParameter, 0.75);
    operator.setWeight(10.0);
    schedule.addOperator(operator);
    operator = new BitFlipOperator(rateIndicatorParameter, 15.0, true);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 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, 0.0077, 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);
    //CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    DistributionLikelihood likelihood3 = new DistributionLikelihood(new GammaDistribution(0.5, 2.0), 0.0);
    likelihood3.addData(ratesParameter);
    DistributionLikelihood likelihood4 = new DistributionLikelihood(new PoissonDistribution(1.0), 0.0);
    likelihood4.addData(rateChanges);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    likelihoods.add(likelihood3);
    likelihoods.add(likelihood4);
    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(prior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(kappa);
    //        loggers[0].add(meanRate);
    loggers[0].add(rateChanges);
    loggers[0].add(coefficientOfVariation);
    loggers[0].add(covariance);
    loggers[0].add(popSize);
    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(meanRate);
    loggers[1].add(rateChanges);
    // 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("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="posterior" value="-1818.26"/>
    //        <expectation name="prior" value="-2.70143"/>
    //        <expectation name="likelihood" value="-1815.56"/>
    //        <expectation name="treeModel.rootHeight" value="6.363E-2"/>
    //        <expectation name="constant.popSize" value="9.67405E-2"/>
    //        <expectation name="hky.kappa" value="30.0394"/>
    //        <expectation name="coefficientOfVariation" value="7.02408E-2"/>
    //        covariance 0.47952
    //        <expectation name="rateChangeCount" value="0.40786"/>
    //        <expectation name="coalescent" value="7.29521"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -1818.26);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.PRIOR));
    assertExpectation(CompoundLikelihoodParser.PRIOR, likelihoodStats, -2.70143);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.56);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 6.363E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 30.0394);
    TraceCorrelation rateChangeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("rateChangeCount"));
    assertExpectation("rateChangeCount", rateChangeStats, 0.40786);
    TraceCorrelation coefficientOfVariationStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(RateStatisticParser.COEFFICIENT_OF_VARIATION));
    assertExpectation(RateStatisticParser.COEFFICIENT_OF_VARIATION, coefficientOfVariationStats, 7.02408E-2);
    TraceCorrelation covarianceStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("covariance"));
    assertExpectation("covariance", covarianceStats, 0.47952);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 9.67405E-2);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, 7.29521);
}

Example 13 with GammaSiteModel

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

the class ALSSiteModelParser 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 14 with GammaSiteModel

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

the class SeqGen method main.

public static void main(String[] argv) {
    String treeFileName = argv[0];
    String outputFileStem = argv[1];
    int length = 500;
    double[] frequencies = new double[] { 0.25, 0.25, 0.25, 0.25 };
    double kappa = 10.0;
    double alpha = 0.5;
    double substitutionRate = argv.length < 3 ? 1.0E-3 : Double.parseDouble(argv[2]);
    int categoryCount = argv.length < 4 ? 8 : Integer.parseInt(argv[3]);
    //1.56E-6;
    double damageRate = argv.length < 5 ? 0 : Double.parseDouble(argv[4]);
    System.out.println("substitutionRate = " + substitutionRate + "; categoryCount = " + categoryCount + "; damageRate = " + damageRate);
    FrequencyModel freqModel = new FrequencyModel(dr.evolution.datatype.Nucleotides.INSTANCE, frequencies);
    HKY hkyModel = new HKY(kappa, freqModel);
    SiteModel siteModel = null;
    if (categoryCount > 1) {
        siteModel = new GammaSiteModel(hkyModel, alpha, categoryCount);
    } else {
        // no rate heterogeneity
        siteModel = new GammaSiteModel(hkyModel);
    }
    List<Tree> trees = new ArrayList<Tree>();
    FileReader reader = null;
    try {
        reader = new FileReader(treeFileName);
        //            TreeImporter importer = new NexusImporter(reader);
        TreeImporter importer = new NewickImporter(reader);
        while (importer.hasTree()) {
            Tree tree = importer.importNextTree();
            trees.add(tree);
            System.out.println("tree height = " + tree.getNodeHeight(tree.getRoot()) + "; leave nodes = " + tree.getExternalNodeCount());
        }
    } catch (FileNotFoundException e) {
        e.printStackTrace();
        return;
    } catch (Importer.ImportException e) {
        e.printStackTrace();
        return;
    } catch (IOException e) {
        e.printStackTrace();
        return;
    }
    SeqGen seqGen = new SeqGen(length, substitutionRate, freqModel, hkyModel, siteModel, damageRate);
    int i = 1;
    for (Tree tree : trees) {
        Alignment alignment = seqGen.simulate(tree);
        FileWriter writer = null;
        try {
            //                writer = new FileWriter(outputFileStem + (i < 10 ? "00" : (i < 100 ? "0" : "")) + i + ".nex");
            //                NexusExporter exporter = new NexusExporter(writer);
            //
            //                exporter.exportAlignment(alignment);
            //
            //                writer.close();
            String outputFileName = outputFileStem + "-" + substitutionRate + ".fasta";
            writer = new FileWriter(outputFileName);
            BufferedWriter bf = new BufferedWriter(writer);
            FastaExporter exporter = new FastaExporter(bf);
            exporter.exportSequences(alignment.getSequenceList());
            bf.close();
            System.out.println("Write " + i + "th sequence file : " + outputFileName);
            i++;
        } catch (IOException e) {
            e.printStackTrace();
            return;
        }
    }
}

Example 15 with GammaSiteModel

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

the class MicrosatelliteSimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Microsatellite msatDataType = (Microsatellite) xo.getChild(Microsatellite.class);
    Taxa taxa = (Taxa) xo.getChild(Taxa.class);
    Tree tree = (Tree) xo.getChild(Tree.class);
    MicrosatelliteModel msatModel = (MicrosatelliteModel) xo.getChild(MicrosatelliteModel.class);
    BranchRateModel brModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (brModel == null) {
        brModel = new DefaultBranchRateModel();
    }
    MicrosatelliteSimulator msatSim = new MicrosatelliteSimulator(msatDataType, taxa, tree, new GammaSiteModel(msatModel), brModel);
    Patterns patterns = msatSim.simulateMsatPattern();
    String msatPatId = xo.getAttribute("id", "simMsatPat");
    patterns.setId(msatPatId);
    MicrosatellitePatternParser.printDetails(patterns);
    MicrosatellitePatternParser.printMicrosatContent(patterns);
    return patterns;
}