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

Example 1 with ExponentialBSPGrowth

use of dr.evolution.coalescent.ExponentialBSPGrowth 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;
}

Example 2 with ExponentialBSPGrowth

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

the class ExponentialSkythingLikelihood 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;
    ExponentialBSPGrowth eg = new ExponentialBSPGrowth(Type.YEARS);
    double startGroupPopSize = 0;
    double endGroupPopSize;
    for (int j = intervalCount - 1; j >= 0; j--) {
        if (j == intervalCount - 1) {
            endGroupPopSize = startingPopSize.getParameterValue(0);
        } else {
            endGroupPopSize = startGroupPopSize;
        }
        double startTime = currentTime;
        double endTime = currentTime + intervals[j];
        startGroupPopSize = endGroupPopSize * Math.exp(slopeParameter.getParameterValue(j) * (endTime - startTime));
        eg.setupN1(endGroupPopSize, slopeParameter.getParameterValue(j), 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];
    }
    return logL;
}