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

Example 1 with ConstantPopulation

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

the class VariableSkylineLikelihood method getLogLikelihood.

// **************************************************************
// Likelihood IMPLEMENTATION
// **************************************************************
/**
     * Calculates the log likelihood of this set of coalescent intervals,
     * given a demographic model.
     */
public synchronized double getLogLikelihood() {
    insureValid();
    double logL = 0.0;
    double currentTime = 0.0;
    // index of current interval in the population function
    int groupIndex = 0;
    // how many intervals precessed in the current group - tell when to switch to next group
    int subIndex = 0;
    ConstantPopulation cp = new ConstantPopulation(Units.Type.YEARS);
    final double[] pop = new double[1];
    if (type == Type.LINEAR) {
        cp = new ConstantPopulation(Units.Type.YEARS) {

            public double getDemographic(double t) {
                return pop[0];
            }
        };
    }
    for (int j = 0; j < intervalCount; j++) {
        // set the population size to the size of the middle of the current interval
        if (type == Type.LINEAR) {
            pop[0] = getPopSize(groupIndex, currentTime + intervals[j]);
        }
        cp.setN0(getPopSize(groupIndex, currentTime + (intervals[j] / 2.0)));
        final CoalescentEventType iType = getIntervalType(j);
        if (iType == CoalescentEventType.COALESCENT) {
            subIndex += 1;
            if (subIndex >= groupSizes.get(groupIndex)) {
                groupIndex += 1;
                subIndex = 0;
            }
        }
        logL += calculateIntervalLikelihood(cp, intervals[j], currentTime, lineageCounts[j], iType);
        if (logL > 0 && j > 1) {
            System.out.println(logL);
        }
        // insert zero-length coalescent intervals
        final int diff = getCoalescentEvents(j) - 1;
        for (int k = 0; k < diff; k++) {
            // not clear, seems wrong - how can population size change in 0 time?
            pop[0] = getPopSize(groupIndex, currentTime);
            cp.setN0(pop[0]);
            logL += calculateIntervalLikelihood(cp, 0.0, currentTime, lineageCounts[j] - k - 1, CoalescentEventType.COALESCENT);
            subIndex += 1;
            if (subIndex >= groupSizes.get(groupIndex)) {
                groupIndex += 1;
                subIndex = 0;
            }
        }
        currentTime += intervals[j];
    }
    return logL;
}

Example 2 with ConstantPopulation

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

the class BayesianSkylineGibbsOperator method doOperation.

public double doOperation() {
    if (!bayesianSkylineLikelihood.getIntervalsKnown())
        bayesianSkylineLikelihood.setupIntervals();
    assert (populationSizeParameter != null);
    assert (groupSizeParameter != null);
    // Now enter a loop similar to the likelihood calculation, except that
    // we now just collect
    // waiting times, and exponents of the population size parameters, for
    // each group.
    int groupIndex = 0;
    int subIndex = 0;
    int[] groupSizes = bayesianSkylineLikelihood.getGroupSizes();
    double[] groupEnds = bayesianSkylineLikelihood.getGroupHeights();
    ConstantPopulation cp = new ConstantPopulation(Units.Type.YEARS);
    double exponent = 0.0;
    double gammaRate = 0.0;
    double currentTime = 0.0;
    double[] populationSizes = new double[groupSizes.length];
    double[] exponents = new double[groupSizes.length];
    double[] gammaRates = new double[groupSizes.length];
    // distributions describing the posteriors on N
    for (int j = 0; j < bayesianSkylineLikelihood.getIntervalCount(); j++) {
        // set the population size to the size of the middle of the current
        // interval
        double curpopsize = bayesianSkylineLikelihood.getPopSize(groupIndex, currentTime + (bayesianSkylineLikelihood.getInterval(j) / 2.0), groupEnds);
        populationSizes[groupIndex] = curpopsize;
        cp.setN0(curpopsize);
        exponent += bayesianSkylineLikelihood.calculateIntervalShapeParameter(cp, bayesianSkylineLikelihood.getInterval(j), currentTime, bayesianSkylineLikelihood.getLineageCount(j), bayesianSkylineLikelihood.getIntervalType(j));
        gammaRate += bayesianSkylineLikelihood.calculateIntervalRateParameter(cp, bayesianSkylineLikelihood.getInterval(j), currentTime, bayesianSkylineLikelihood.getLineageCount(j), bayesianSkylineLikelihood.getIntervalType(j)) * curpopsize;
        if (bayesianSkylineLikelihood.getIntervalType(j) == OldAbstractCoalescentLikelihood.CoalescentEventType.COALESCENT) {
            subIndex += 1;
            if (subIndex >= groupSizes[groupIndex]) {
                exponents[groupIndex] = exponent;
                gammaRates[groupIndex] = gammaRate;
                // Reset accumulators
                exponent = 0.0;
                gammaRate = 0.0;
                // Next group
                groupIndex += 1;
                subIndex = 0;
            }
        }
        // insert zero-length coalescent intervals (for multifurcating trees
        int diff = bayesianSkylineLikelihood.getCoalescentEvents(j) - 1;
        for (int k = 0; k < diff; k++) {
            curpopsize = bayesianSkylineLikelihood.getPopSize(groupIndex, currentTime + (bayesianSkylineLikelihood.getInterval(j) / 2.0), groupEnds);
            populationSizes[groupIndex] = curpopsize;
            cp.setN0(curpopsize);
            exponent += bayesianSkylineLikelihood.calculateIntervalShapeParameter(cp, bayesianSkylineLikelihood.getInterval(j), currentTime, bayesianSkylineLikelihood.getLineageCount(j), bayesianSkylineLikelihood.getIntervalType(j));
            gammaRate += bayesianSkylineLikelihood.calculateIntervalRateParameter(cp, bayesianSkylineLikelihood.getInterval(j), currentTime, bayesianSkylineLikelihood.getLineageCount(j), bayesianSkylineLikelihood.getIntervalType(j)) * curpopsize;
            subIndex += 1;
            if (subIndex >= groupSizes[groupIndex]) {
                exponents[groupIndex] = exponent;
                gammaRates[groupIndex] = gammaRate;
                // Reset accumulators
                exponent = 0.0;
                gammaRate = 0.0;
                // Next group
                groupIndex += 1;
                subIndex = 0;
            }
        }
        currentTime += bayesianSkylineLikelihood.getInterval(j);
    }
    // Next, sample new population sizes
    double hastingsRatio = 0.0;
    int numGibbsMoves = iterations;
    boolean[] tabu = new boolean[groupSizes.length];
    for (int i = 0; i < groupSizes.length; i++) tabu[i] = false;
    while (numGibbsMoves > 0) {
        int randomCoord = dr.math.MathUtils.nextInt(groupSizes.length);
        // neighbours haven't been sampled
        if (!tabu[randomCoord]) {
            hastingsRatio += getSample(exponents, gammaRates, populationSizes, randomCoord);
            // make current index tabu, but un-tabuize its neighbours
            tabu[randomCoord] = true;
            if (randomCoord != 0)
                tabu[randomCoord - 1] = false;
            if (randomCoord != groupSizes.length - 1)
                tabu[randomCoord + 1] = false;
        }
        numGibbsMoves -= 1;
    }
    // store the new array of population sizes
    for (int j = 0; j < groupIndex; j++) {
        populationSizeParameter.setParameterValue(j, populationSizes[j]);
    }
    return hastingsRatio;
}

Example 3 with ConstantPopulation

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

the class SkylineLikelihood method calculateLogLikelihood.

// **************************************************************
// Likelihood IMPLEMENTATION
// **************************************************************
/**
	 * Calculates the log likelihood of this set of coalescent intervals,
	 * given a demographic model.
	 */
public double calculateLogLikelihood() {
    if (!intervalsKnown)
        setupIntervals();
    double logL = 0.0;
    double currentTime = 0.0;
    int popIndex = 0;
    ConstantPopulation cp = new ConstantPopulation(Units.Type.YEARS);
    for (int j = 0; j < intervalCount; j++) {
        cp.setN0(popSizeParameter.getParameterValue(popIndex));
        if (getIntervalType(j) == CoalescentEventType.COALESCENT) {
            popIndex += 1;
        }
        logL += calculateIntervalLikelihood(cp, intervals[j], currentTime, lineageCounts[j], getIntervalType(j));
        // insert zero-length coalescent intervals
        int diff = getCoalescentEvents(j) - 1;
        for (int k = 0; k < diff; k++) {
            cp.setN0(popSizeParameter.getParameterValue(popIndex));
            logL += calculateIntervalLikelihood(cp, 0.0, currentTime, lineageCounts[j] - k - 1, CoalescentEventType.COALESCENT);
            popIndex += 1;
        }
        currentTime += intervals[j];
    }
    return logL;
}

Example 4 with ConstantPopulation

use of dr.evolution.coalescent.ConstantPopulation 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 5 with ConstantPopulation

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

the class BayesianSkylineLikelihood method getLogLikelihood.

// **************************************************************
// Likelihood IMPLEMENTATION
// **************************************************************
/**
     * Calculates the log likelihood of this set of coalescent intervals,
     * given a demographic model.
     */
public double getLogLikelihood() {
    setupIntervals();
    double logL = 0.0;
    double currentTime = 0.0;
    int groupIndex = 0;
    int[] groupSizes = getGroupSizes();
    double[] groupEnds = getGroupHeights();
    int subIndex = 0;
    if (type == EXPONENTIAL_TYPE) {
        ExponentialBSPGrowth eg = new ExponentialBSPGrowth(Units.Type.YEARS);
        for (int j = 0; j < intervalCount; j++) {
            double startGroupPopSize = popSizeParameter.getParameterValue(groupIndex);
            double endGroupPopSize = popSizeParameter.getParameterValue(groupIndex + 1);
            double startTime = currentTime;
            double endTime = currentTime + intervals[j];
            eg.setup(startGroupPopSize, endGroupPopSize, endTime - startTime);
            if (getIntervalType(j) == CoalescentEventType.COALESCENT) {
                subIndex += 1;
                if (subIndex >= groupSizes[groupIndex]) {
                    groupIndex += 1;
                    subIndex = 0;
                }
            }
            logL += calculateIntervalLikelihood(eg, intervals[j], currentTime, lineageCounts[j], getIntervalType(j));
            // insert zero-length coalescent intervals
            int diff = getCoalescentEvents(j) - 1;
            for (int k = 0; k < diff; k++) {
                eg.setup(startGroupPopSize, startGroupPopSize, endTime - startTime);
                logL += calculateIntervalLikelihood(eg, 0.0, currentTime, lineageCounts[j] - k - 1, CoalescentEventType.COALESCENT);
                subIndex += 1;
                if (subIndex >= groupSizes[groupIndex]) {
                    groupIndex += 1;
                    subIndex = 0;
                }
            }
            currentTime += intervals[j];
        }
    } else {
        ConstantPopulation cp = new ConstantPopulation(Units.Type.YEARS);
        for (int j = 0; j < intervalCount; j++) {
            // set the population size to the size of the middle of the current interval
            final double ps = getPopSize(groupIndex, currentTime + (intervals[j] / 2.0), groupEnds);
            cp.setN0(ps);
            if (getIntervalType(j) == CoalescentEventType.COALESCENT) {
                subIndex += 1;
                if (subIndex >= groupSizes[groupIndex]) {
                    groupIndex += 1;
                    subIndex = 0;
                }
            }
            logL += calculateIntervalLikelihood(cp, intervals[j], currentTime, lineageCounts[j], getIntervalType(j));
            // insert zero-length coalescent intervals
            int diff = getCoalescentEvents(j) - 1;
            for (int k = 0; k < diff; k++) {
                cp.setN0(getPopSize(groupIndex, currentTime, groupEnds));
                logL += calculateIntervalLikelihood(cp, 0.0, currentTime, lineageCounts[j] - k - 1, CoalescentEventType.COALESCENT);
                subIndex += 1;
                if (subIndex >= groupSizes[groupIndex]) {
                    groupIndex += 1;
                    subIndex = 0;
                }
            }
            currentTime += intervals[j];
        }
    }
    return logL;
}