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

Example 1 with ScaledDemographic

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

the class OldAbstractCoalescentLikelihood method calculateLogLikelihood.

/**
     * @return the log likelihood of this set of coalescent intervals,
     *         given a demographic model
     */
public double calculateLogLikelihood() {
    if (treesSet != null) {
        final int nTrees = treesSet.nLoci();
        final DemographicFunction demogFunction = demoModel.getDemographicFunction();
        double logLike = 0.0;
        for (int nt = 0; nt < nTrees; ++nt) {
            final double popFactor = treesSet.getPopulationFactor(nt);
            DemographicFunction df = popFactor != 1.0 ? new ScaledDemographic(demogFunction, popFactor) : demogFunction;
            logLike += Coalescent.calculateLogLikelihood(treesSet.getTreeIntervals(nt), df);
        }
        return logLike;
    }
    if (!intervalsKnown)
        setupIntervals();
    if (demoModel == null)
        return calculateAnalyticalLogLikelihood();
    double logL = 0.0;
    double currentTime = 0.0;
    DemographicFunction demoFunction = demoModel.getDemographicFunction();
    for (int j = 0; j < intervalCount; j++) {
        logL += calculateIntervalLikelihood(demoFunction, intervals[j], currentTime, lineageCounts[j], getIntervalType(j));
        // insert zero-length coalescent intervals
        final int diff = getCoalescentEvents(j) - 1;
        for (int k = 0; k < diff; k++) {
            logL += calculateIntervalLikelihood(demoFunction, 0.0, currentTime, lineageCounts[j] - k - 1, CoalescentEventType.COALESCENT);
        }
        currentTime += intervals[j];
    }
    return logL;
}

Example 2 with ScaledDemographic

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

the class MulMSCoalescent method treeLogLikelihood.

private double treeLogLikelihood(MulSpeciesBindings.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 MulSpeciesBindings.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) {
            System.err.println("ERROR in evomodel.speciation.MulMSCoalescent.treeLogLikelihood()");
        }
        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 3 with ScaledDemographic

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