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

Example 6 with DemographicFunction

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

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

the class SpeciesTreeSimplePrior method calculateLogLikelihood.

protected double calculateLogLikelihood() {
    double ll = 0;
    final NodeRef root = sTree.getRoot();
    double l = sTree.getNodeHeight(root) + ((VDdemographicFunction) sTree.getNodeDemographic(root)).naturalLimit();
    for (int nn = 0; nn < sTree.getNodeCount(); ++nn) {
        final NodeRef n = sTree.getNode(nn);
        final DemographicFunction demog = sTree.getNodeDemographic(n);
        if (sTree.isExternal(n)) {
            ll += tips.logPdf(demog.getDemographic(0));
        }
        final double branch = sTree.isRoot(n) ? ((VDdemographicFunction) demog).naturalLimit() : sTree.getBranchLength(n);
        final double avg = branch / demog.getIntegral(0, branch);
        final double p = demog.getDemographic(branch);
        final double s = sigma.getParameterValue(0) * Math.sqrt((1 - Math.exp(-branch / l)) / d1);
        //final double vx = (new LogNormalDistribution(-0.5 * s * s, s)).logPdf(p / avg);
        double z = p / avg;
        final double v2 = Math.log(z);
        final double v1 = v2 / s + s / 2;
        final double v = -Math.log(s) - f2 - v2 - v1 * v1 / 2;
        ll += v;
    //ll += dist.logPdf(p/avg);
    //
    //            double[] pa = {demog.getDemographic(0), demog.getDemographic(branch)};
    //            for( double p : pa ) {
    //               ll += dist.logPdf(p/avg);
    //            }
    }
    return ll;
}

Example 8 with DemographicFunction

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

the class CoalescentLikelihood method calculateLogLikelihood.

// **************************************************************
// Likelihood IMPLEMENTATION
// **************************************************************
/**
	 * Calculates the log likelihood of this set of coalescent intervals,
	 * given a demographic model.
	 */
public double calculateLogLikelihood() {
    DemographicFunction demoFunction = demoModel.getDemographicFunction();
    //double lnL =  Coalescent.calculateLogLikelihood(getIntervals(), demoFunction);
    double lnL = Coalescent.calculateLogLikelihood(getIntervals(), demoFunction, demoFunction.getThreshold());
    if (Double.isNaN(lnL) || Double.isInfinite(lnL)) {
        Logger.getLogger("warning").severe("CoalescentLikelihood for " + demoModel.getId() + " is " + Double.toString(lnL));
    }
    return lnL;
}

Example 9 with DemographicFunction

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

the class MultiSpeciesCoalescent method treeLogLikelihood.

private double treeLogLikelihood(SpeciesBindings.GeneTreeInfo geneTree, NodeRef node, int[] info, double popFactor) {
    // number of lineages remaining at node
    int nLineages;
    // location in coalescent list (optimization)
    int indexInClist = 0;
    // accumulated log-likelihood inBranchh from node to it's parent
    double like = 0;
    final double t0 = spTree.getNodeHeight(node);
    final SpeciesBindings.CoalInfo[] cList = geneTree.getCoalInfo();
    if (verbose) {
        if (spTree.isRoot(node)) {
            System.err.println("gtree:" + geneTree.tree.getId());
            System.err.println("t0 " + t0);
            for (int k = 0; k < cList.length; ++k) {
                System.err.println(k + " " + cList[k].ctime + " " + cList[k].sinfo[0] + " " + cList[k].sinfo[1]);
            }
        }
    }
    if (spTree.isExternal(node)) {
        nLineages = geneTree.nLineages(spTree.speciesIndex(node));
        indexInClist = 0;
    } else {
        //assert spTree.getChildCount(node) == 2;
        nLineages = 0;
        for (int nc = 0; nc < 2; ++nc) {
            final NodeRef child = spTree.getChild(node, nc);
            like += treeLogLikelihood(geneTree, child, info, popFactor);
            nLineages += info[0];
            indexInClist = Math.max(indexInClist, info[1]);
        }
        // root of species tree
        if (indexInClist >= cList.length) {
            int k = 1;
        }
        assert indexInClist < cList.length;
        while (cList[indexInClist].ctime < t0) {
            ++indexInClist;
        }
    }
    final boolean isRoot = spTree.isRoot(node);
    // Upper limit
    // use of (t0 + spTree.getBranchLength(node)) caused problem since there was a tiny difference
    // between those (supposedly equal) values. we should track where the discrepancy comes from.
    final double stopTime = isRoot ? Double.MAX_VALUE : spTree.getNodeHeight(spTree.getParent(node));
    // demographic function is 0 based (relative to node height)
    // time away from node
    double lastTime = 0.0;
    // demographic function across branch
    DemographicFunction demog = spTree.getNodeDemographic(node);
    if (popFactor > 0) {
        demog = new ScaledDemographic(demog, popFactor);
    }
    //        if(false) {
    //            final double duration = isRoot ? cList[cList.length - 1].ctime - t0 : stopTime;
    //            double demographicAtCoalPoint = demog.getDemographic(duration);
    //            double intervalArea = demog.getIntegral(0, duration);
    //            if( demographicAtCoalPoint < 1e-12 * (duration/intervalArea) ) {
    //               return Double.NEGATIVE_INFINITY;
    //            }
    //        }
    // Species sharing this branch
    FixedBitSet subspeciesSet = spTree.spSet(node);
    if (verbose) {
        System.err.println(TreeUtils.uniqueNewick(spTree, node) + " nl " + nLineages + " " + subspeciesSet + " t0 - st " + t0 + " - " + stopTime);
    }
    while (nLineages > 1) {
        //            }
        assert (indexInClist < cList.length);
        final double nextT = cList[indexInClist].ctime;
        // while rare they can be equal
        if (nextT >= stopTime) {
            break;
        }
        if (nonEmptyIntersection(cList[indexInClist].sinfo, subspeciesSet)) {
            final double time = nextT - t0;
            if (time > 0) {
                final double interval = demog.getIntegral(lastTime, time);
                lastTime = time;
                final int nLineageOver2 = (nLineages * (nLineages - 1)) / 2;
                like -= nLineageOver2 * interval;
                final double pop = demog.getDemographic(time);
                assert (pop > 0);
                like -= Math.log(pop);
            }
            --nLineages;
        }
        ++indexInClist;
    }
    if (nLineages > 1) {
        // add term for No coalescent until root
        final double interval = demog.getIntegral(lastTime, stopTime - t0);
        final int nLineageOver2 = (nLineages * (nLineages - 1)) / 2;
        like -= nLineageOver2 * interval;
    }
    info[0] = nLineages;
    info[1] = indexInClist;
    if (verbose) {
        System.err.println(TreeUtils.uniqueNewick(spTree, node) + " stopTime " + stopTime + " nl " + nLineages + " icl " + indexInClist);
    }
    return like;
}

Example 10 with DemographicFunction

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