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

Example 1 with IntervalList

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

the class CoalescentLikelihood method calculateLogLikelihood.

protected double calculateLogLikelihood(PopulationSizeFunction populationSizeFunction) {
    double logL = 0.0;
    IntervalList intervals = getIntervalList();
    final int n = intervals.getIntervalCount();
    if (n == 0) {
        return 0.0;
    }
    double absoluteStartTime = intervals.getStartTime();
    demographicModel.setTimeOffset(absoluteStartTime);
    DemographicFunction demographicFunction = demographicModel.getDemographicFunction();
    double startTime = 0;
    for (int i = 0; i < n; i++) {
        final double duration = intervals.getInterval(i);
        final double finishTime = startTime + duration;
        final double intervalArea = populationSizeFunction.getIntegral(startTime, finishTime);
        if (intervalArea == 0 && duration != 0) {
            return Double.NEGATIVE_INFINITY;
        }
        final int lineageCount = intervals.getLineageCount(i);
        final double kChoose2 = Binomial.choose2(lineageCount);
        // common part
        logL += -kChoose2 * intervalArea;
        if (intervals.getIntervalType(i) == IntervalType.COALESCENT) {
            logL -= populationSizeFunction.getLogDemographic(finishTime);
        }
        startTime = finishTime;
    }
    return logL;
}

Example 2 with IntervalList

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

the class GMRFSkygridLikelihood method calculateLogCoalescentLikelihood.

protected double calculateLogCoalescentLikelihood() {
    if (!intervalsKnown) {
        // intervalsKnown -> false when handleModelChanged event occurs in super.
        wrapSetupIntervals();
        for (IntervalList intervalList : intervalsList) {
            intervalList.calculateIntervals();
        }
        setupSufficientStatistics();
        intervalsKnown = true;
    }
    // Matrix operations taken from block update sampler to calculate data likelihood and field prior
    double currentLike = 0;
    double[] currentGamma = popSizeParameter.getParameterValues();
    for (int i = 0; i < fieldLength; i++) {
        currentLike += -numCoalEvents[i] * currentGamma[i] + ploidySums[i] - sufficientStatistics[i] * Math.exp(-currentGamma[i]);
    }
    return currentLike;
}

Example 3 with IntervalList

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

the class GMRFSkygridLikelihood method handleModelChangedEvent.

protected void handleModelChangedEvent(Model model, Object object, int index) {
    if (model instanceof IntervalList) {
        IntervalList intervalList = (IntervalList) model;
        int tn = intervalsList.indexOf(intervalList);
        if (tn >= 0) {
            intervalsKnown = false;
            likelihoodKnown = false;
        } else {
            throw new RuntimeException("Unknown tree modified in GMRFSkygridLikelihood");
        }
    } else {
        throw new RuntimeException("Unknown object modified in GMRFSkygridLikelihood");
    }
}

Example 4 with IntervalList

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

the class SmoothSkygridLikelihood method computeSufficientStatistics.

private void computeSufficientStatistics() {
    for (IntervalList intervals : intervalList) {
        // Find first grid-interval for tree
        double previousIntervalTimeOnTime = intervals.getStartTime();
        int currentGridIndex = 0;
        while (previousIntervalTimeOnTime > gridPoints[currentGridIndex]) {
            ++currentGridIndex;
        }
        // Interate over tree-intervals in time-order
        for (int j = 0; j < intervals.getIntervalCount() - 1; ++j) {
            double currentIntervalTimeOnTree = intervals.getIntervalTime(j);
            int currentIntervalLineageCount = intervals.getLineageCount(j);
            while (intervals.getIntervalTime(j + 1) > gridPoints[currentGridIndex]) {
            // Do somethimg until end of internval
            }
        // Do something with last bit of time
        }
    // Handle last tree-interval here
    }
    throw new RuntimeException("Not yet implemented");
}

Example 5 with IntervalList

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

the class CoalescentLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(MODEL);
    PopulationSizeModel populationSizeModel = (PopulationSizeModel) cxo.getChild(PopulationSizeModel.class);
    DemographicModel demographicModel = null;
    if (populationSizeModel == null) {
        demographicModel = (DemographicModel) cxo.getChild(DemographicModel.class);
    }
    List<TreeModel> trees = new ArrayList<TreeModel>();
    List<Double> popFactors = new ArrayList<Double>();
    MultiLociTreeSet treesSet = (demographicModel instanceof MultiLociTreeSet ? (MultiLociTreeSet) demographicModel : null);
    IntervalList intervalList = null;
    for (int k = 0; k < xo.getChildCount(); ++k) {
        final Object child = xo.getChild(k);
        if (child instanceof XMLObject) {
            cxo = (XMLObject) child;
            if (cxo.getName().equals(POPULATION_TREE)) {
                final TreeModel t = (TreeModel) cxo.getChild(TreeModel.class);
                assert t != null;
                trees.add(t);
                popFactors.add(cxo.getAttribute(POPULATION_FACTOR, 1.0));
            } else if (cxo.getName().equals(INTERVALS)) {
                intervalList = (IntervalList) cxo.getChild(IntervalList.class);
            }
        }
    // in the future we may have arbitrary multi-loci element
    // else if( child instanceof MultiLociTreeSet )  {
    // treesSet = (MultiLociTreeSet)child;
    // }
    }
    TreeModel treeModel = null;
    if (trees.size() == 1 && popFactors.get(0) == 1.0) {
        treeModel = trees.get(0);
    } else if (trees.size() > 1) {
        treesSet = new MultiLociTreeSet.Default(trees, popFactors);
    } else if (!(trees.size() == 0 && treesSet != null) && intervalList == null) {
        throw new XMLParseException("Incorrectly constructed likelihood element");
    }
    if (treeModel != null) {
        TaxonList includeSubtree = null;
        if (xo.hasChildNamed(INCLUDE)) {
            includeSubtree = (TaxonList) xo.getElementFirstChild(INCLUDE);
        }
        List<TaxonList> excludeSubtrees = new ArrayList<TaxonList>();
        if (xo.hasChildNamed(EXCLUDE)) {
            cxo = xo.getChild(EXCLUDE);
            for (int i = 0; i < cxo.getChildCount(); i++) {
                excludeSubtrees.add((TaxonList) cxo.getChild(i));
            }
        }
        try {
            intervalList = new TreeIntervals(treeModel, includeSubtree, excludeSubtrees);
        // TreeIntervals now deals with all the interval stuff
        // return new CoalescentLikelihood(treeModel, includeSubtree, excludeSubtrees, demoModel);
        } catch (TreeUtils.MissingTaxonException mte) {
            throw new XMLParseException("treeModel missing a taxon from taxon list in " + getParserName() + " element");
        }
    }
    if (intervalList != null) {
        if (demographicModel != null) {
            return new CoalescentLikelihood(intervalList, demographicModel);
        }
        return new CoalescentLikelihood(intervalList, populationSizeModel);
    } else {
        if (xo.hasChildNamed(INCLUDE) || xo.hasChildNamed(EXCLUDE)) {
            throw new XMLParseException("Include/Exclude taxa not supported for multi locus sets");
        }
        // a base - and modifing it will probsbly result in a bigger mess.
        return new OldAbstractCoalescentLikelihood(treesSet, demographicModel);
    }
}