Examples with ExponentialGrowth dr.evolution.coalescent.ExponentialGrowth used on opensource projects

Example 1 with ExponentialGrowth

use of dr.evolution.coalescent.ExponentialGrowth in project beast-mcmc by beast-dev.

the class TransmissionTreeToVirusTree method main.

public static void main(String[] args) {
    ModelType model = ModelType.CONSTANT;
    double startNe = 1;
    double growthRate = 0;
    double t50 = 0;
    Arguments arguments = new Arguments(new Arguments.Option[] { new Arguments.StringOption(DEMOGRAPHIC_MODEL, demographics, false, "The type of within-host" + " demographic function to use, default = constant"), new Arguments.RealOption(STARTING_POPULATION_SIZE, "The effective population size at time zero" + " (used in all models), default = 1"), new Arguments.RealOption(GROWTH_RATE, "The effective population size growth rate (used in" + " exponential and logistic models), default = 0"), new Arguments.RealOption(T50, "The time point, relative to the time of infection in backwards " + "time, at which the population is equal to half its final asymptotic value, in the " + "logistic model default = 0") });
    try {
        arguments.parseArguments(args);
    } catch (Arguments.ArgumentException ae) {
        System.out.println(ae);
        printUsage(arguments);
        System.exit(1);
    }
    if (arguments.hasOption(HELP)) {
        printUsage(arguments);
        System.exit(0);
    }
    if (arguments.hasOption(DEMOGRAPHIC_MODEL)) {
        String modelString = arguments.getStringOption(DEMOGRAPHIC_MODEL);
        if (modelString.toLowerCase().startsWith("c")) {
            model = ModelType.CONSTANT;
        } else if (modelString.toLowerCase().startsWith("e")) {
            model = ModelType.EXPONENTIAL;
        } else if (modelString.toLowerCase().startsWith("l")) {
            model = ModelType.LOGISTIC;
        } else {
            progressStream.print("Unrecognised demographic model type");
            System.exit(1);
        }
    }
    if (arguments.hasOption(STARTING_POPULATION_SIZE)) {
        startNe = arguments.getRealOption(STARTING_POPULATION_SIZE);
    }
    if (arguments.hasOption(GROWTH_RATE) && model != ModelType.CONSTANT) {
        growthRate = arguments.getRealOption(GROWTH_RATE);
    }
    if (arguments.hasOption(T50) && model == ModelType.LOGISTIC) {
        t50 = arguments.getRealOption(T50);
    }
    DemographicFunction demoFunction = null;
    switch(model) {
        case CONSTANT:
            {
                demoFunction = new ConstantPopulation(Units.Type.YEARS);
                ((ConstantPopulation) demoFunction).setN0(startNe);
            }
        case EXPONENTIAL:
            {
                demoFunction = new ExponentialGrowth(Units.Type.YEARS);
                ((ExponentialGrowth) demoFunction).setN0(startNe);
                ((ExponentialGrowth) demoFunction).setGrowthRate(growthRate);
            }
        case LOGISTIC:
            {
                demoFunction = new LogisticGrowthN0(Units.Type.YEARS);
                ((LogisticGrowthN0) demoFunction).setN0(startNe);
                ((LogisticGrowthN0) demoFunction).setGrowthRate(growthRate);
                ((LogisticGrowthN0) demoFunction).setT50(t50);
            }
    }
    final String[] args2 = arguments.getLeftoverArguments();
    if (args2.length != 3) {
        printUsage(arguments);
        System.exit(1);
    }
    String infectionsFileName = args2[0];
    String samplesFileName = args2[1];
    String outputFileRoot = args2[2];
    TransmissionTreeToVirusTree instance = new TransmissionTreeToVirusTree(samplesFileName, infectionsFileName, demoFunction, outputFileRoot);
    try {
        instance.run();
    } catch (IOException e) {
        e.printStackTrace();
    }
}

Example 2 with ExponentialGrowth

use of dr.evolution.coalescent.ExponentialGrowth in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateTopology.

// END: annotateTree
static void simulateTopology() {
    try {
        System.out.println("Test case 1: simulateTopology");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        // set demographic function
        ExponentialGrowth exponentialGrowth = new ExponentialGrowth(Units.Type.YEARS);
        exponentialGrowth.setN0(10);
        exponentialGrowth.setGrowthRate(0.5);
        Taxa taxa = new Taxa();
        for (Taxon taxon : tree.asList()) {
            double absoluteHeight = Utils.getAbsoluteTaxonHeight(taxon, tree);
            taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, absoluteHeight);
            // taxon.setAttribute("date", new Date(absoluteHeight,
            // Units.Type.YEARS, true));
            taxa.addTaxon(taxon);
        }
        // END: taxon loop
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        TreeModel treeModel = new TreeModel(topologySimulator.simulateTree(taxa, exponentialGrowth));
        System.out.println(treeModel.toString());
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 3 with ExponentialGrowth

use of dr.evolution.coalescent.ExponentialGrowth in project beast-mcmc by beast-dev.

the class PartitionData method createDemographicFunction.

public DemographicFunction createDemographicFunction() {
    DemographicFunction demographicFunction = null;
    if (this.demographicModelIndex == 0) {
    // No model
    // do nothing
    } else if (this.demographicModelIndex == 1) {
        // Constant Population
        demographicFunction = new ConstantPopulation(Units.Type.YEARS);
        ((ConstantPopulation) demographicFunction).setN0(demographicParameterValues[0]);
    } else if (this.demographicModelIndex == 2) {
        // Exponential Growth (Growth Rate)
        demographicFunction = new ExponentialGrowth(Units.Type.YEARS);
        ((ExponentialGrowth) demographicFunction).setN0(demographicParameterValues[1]);
        ((ExponentialGrowth) demographicFunction).setGrowthRate(demographicParameterValues[2]);
    } else if (this.demographicModelIndex == 3) {
        // Exponential Growth (Doubling Time)
        demographicFunction = new ExponentialGrowth(Units.Type.YEARS);
        ((ExponentialGrowth) demographicFunction).setN0(demographicParameterValues[3]);
        ((ExponentialGrowth) demographicFunction).setDoublingTime(demographicParameterValues[4]);
    //		} else if (this.demographicModelIndex == 4) {// Logistic Growth (Growth Rate)
    //			
    //			demographicFunction = new LogisticGrowth(Units.Type.YEARS);
    //            ((LogisticGrowth) demographicFunction).setN0(demographicParameterValues[5]);
    //            ((LogisticGrowth) demographicFunction).setGrowthRate(demographicParameterValues[6]);
    //            ((LogisticGrowth) demographicFunction).setTime50(demographicParameterValues[7]);
    //			
    //		} else if (this.demographicModelIndex == 5) {// Logistic Growth (Doubling Time)
    //			
    //			demographicFunction = new LogisticGrowth(Units.Type.YEARS);
    //            ((LogisticGrowth) demographicFunction).setN0(demographicParameterValues[8]);
    //            ((LogisticGrowth) demographicFunction).setDoublingTime(demographicParameterValues[9]);
    //            ((LogisticGrowth) demographicFunction).setTime50(demographicParameterValues[10]);
    //			
    //		} else if (this.demographicModelIndex == 6) {// Expansion (Growth Rate)
    //			
    //			demographicFunction = new Expansion(Units.Type.YEARS);
    //            ((Expansion) demographicFunction).setN0(demographicParameterValues[11]);
    //            ((Expansion) demographicFunction).setProportion(demographicParameterValues[12]);
    //            ((Expansion) demographicFunction).setGrowthRate(demographicParameterValues[13]);
    //			
    //		} else if (this.demographicModelIndex == 7) {// Expansion (Doubling Time)
    //			
    //			demographicFunction = new Expansion(Units.Type.YEARS);
    //            ((Expansion) demographicFunction).setN0(demographicParameterValues[14]);
    //            ((Expansion) demographicFunction).setProportion(demographicParameterValues[15]);
    //            ((Expansion) demographicFunction).setDoublingTime(demographicParameterValues[16]);
    } else {
        System.out.println("Not yet implemented");
    }
    return demographicFunction;
}

Example 4 with ExponentialGrowth

use of dr.evolution.coalescent.ExponentialGrowth in project beast-mcmc by beast-dev.

the class MultiEpochExponentialTest method testThreeExponential.

public void testThreeExponential() {
    Units.Type units = Units.Type.YEARS;
    ExponentialGrowth e = new ExponentialGrowth(units);
    e.setN0(N0);
    e.setGrowthRate(rates3[0]);
    MultiEpochExponential mee = new MultiEpochExponential(units, 3);
    mee.setN0(N0);
    for (int i = 0; i < rates3.length; ++i) {
        mee.setGrowthRate(i, rates3[i]);
    }
    for (int i = 0; i < transitionTimes3.length; ++i) {
        mee.setTransitionTime(i, transitionTimes3[i]);
    }
    double start = 0.0;
    double finish = 1.0;
    for (; finish < 20.0; finish += 1.0) {
        double eeInt = e.getIntegral(start, finish);
        double meeIntN = mee.getNumericalIntegral(start, finish);
        double meeIntA = mee.getAnalyticIntegral(start, finish);
        //            System.err.println(finish + ": " + eeInt + " " + meeIntN + " " + meeIntA);
        assertEquals(eeInt, meeIntN, tolerance1);
        assertEquals(meeIntN, meeIntA, tolerance2);
    }
    start = 0.5;
    finish = 1.0;
    for (; finish < 20.0; finish += 1.0) {
        double eeInt = e.getIntegral(start, finish);
        double meeIntN = mee.getNumericalIntegral(start, finish);
        double meeIntA = mee.getAnalyticIntegral(start, finish);
        //            System.err.println(finish + ": " + eeInt + " " + meeIntN + " " + meeIntA);
        assertEquals(eeInt, meeIntN, tolerance1);
        assertEquals(meeIntN, meeIntA, tolerance2);
    }
    start = 11.0;
    finish = 11.0;
    for (; finish < 20.0; finish += 1.0) {
        double eeInt = e.getIntegral(start, finish);
        double meeIntN = mee.getNumericalIntegral(start, finish);
        double meeIntA = mee.getAnalyticIntegral(start, finish);
        //            System.err.println(finish + ": " + eeInt + " " + meeIntN + " " + meeIntA);
        assertEquals(eeInt, meeIntN, tolerance1);
        assertEquals(meeIntN, meeIntA, tolerance2);
    }
}