/*
* Coalescent.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.evolution.coalescent;
import dr.evolution.tree.Tree;
import dr.evolution.util.Units;
import dr.math.Binomial;
import dr.math.MultivariateFunction;
/**
* A likelihood function for the coalescent. Takes a tree and a demographic model.
*
* Parts of this class were derived from C++ code provided by Oliver Pybus.
*
* @version $Id: Coalescent.java,v 1.12 2005/05/24 20:25:55 rambaut Exp $
*
* @author Andrew Rambaut
* @author Alexei Drummond
*/
public class Coalescent implements MultivariateFunction, Units {
// PUBLIC STUFF
public Coalescent(Tree tree, DemographicFunction demographicFunction) {
this(new TreeIntervals(tree), demographicFunction);
}
public Coalescent(IntervalList intervals, DemographicFunction demographicFunction) {
this.intervals = intervals;
this.demographicFunction = demographicFunction;
}
/**
* Calculates the log likelihood of this set of coalescent intervals,
* given a demographic model.
*/
public double calculateLogLikelihood() {
return calculateLogLikelihood(intervals, demographicFunction);
}
/**
* Calculates the log likelihood of this set of coalescent intervals,
* given a demographic model.
*/
public static double calculateLogLikelihood(IntervalList intervals, DemographicFunction demographicFunction) {
return calculateLogLikelihood(intervals, demographicFunction, 0.0);
}
/**
* Calculates the log likelihood of this set of coalescent intervals,
* given a demographic model.
*/
public static double calculateLogLikelihood(IntervalList intervals, DemographicFunction demographicFunction, double threshold) {
double logL = 0.0;
double startTime = 0.0;
final int n = intervals.getIntervalCount();
for (int i = 0; i < n; i++) {
final double duration = intervals.getInterval(i);
final double finishTime = startTime + duration;
final double intervalArea = demographicFunction.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) {
final double demographicAtCoalPoint = demographicFunction.getDemographic(finishTime);
// if value at end is many orders of magnitude different than mean over interval reject the interval
// This is protection against cases where ridiculous infinitesimal population size at the end of a
// linear interval drive coalescent values to infinity.
if( duration == 0.0 || demographicAtCoalPoint * (intervalArea/duration) >= threshold ) {
// if( duration == 0.0 || demographicAtCoalPoint >= threshold * (duration/intervalArea) ) {
logL -= Math.log(demographicAtCoalPoint);
} else {
// remove this at some stage
// System.err.println("Warning: " + i + " " + demographicAtCoalPoint + " " + (intervalArea/duration) );
return Double.NEGATIVE_INFINITY;
}
}
startTime = finishTime;
}
return logL;
}
/**
* Calculates the log likelihood of this set of coalescent intervals,
* using an analytical integration over theta.
*/
public static double calculateAnalyticalLogLikelihood(IntervalList intervals) {
if (!intervals.isCoalescentOnly()) {
throw new IllegalArgumentException("Can only calculate analytical likelihood for pure coalescent intervals");
}
final double lambda = getLambda(intervals);
final int n = intervals.getSampleCount();
// assumes a 1/theta prior
//logLikelihood = Math.log(1.0/Math.pow(lambda,n));
// assumes a flat prior
return (1-n) * Math.log(lambda); // Math.log(1.0/Math.pow(lambda,n-1));
}
/**
* Returns a factor lambda such that the likelihood can be expressed as
* 1/theta^(n-1) * exp(-lambda/theta). This allows theta to be integrated
* out analytically. :-)
*/
private static double getLambda(IntervalList intervals) {
double lambda = 0.0;
for (int i= 0; i < intervals.getIntervalCount(); i++) {
lambda += (intervals.getInterval(i) * intervals.getLineageCount(i));
}
lambda /= 2;
return lambda;
}
// **************************************************************
// MultivariateFunction IMPLEMENTATION
// **************************************************************
public double evaluate(double[] argument) {
for (int i = 0; i < argument.length; i++) {
demographicFunction.setArgument(i, argument[i]);
}
return calculateLogLikelihood();
}
public int getNumArguments() {
return demographicFunction.getNumArguments();
}
public double getLowerBound(int n) {
return demographicFunction.getLowerBound(n);
}
public double getUpperBound(int n) {
return demographicFunction.getUpperBound(n);
}
// **************************************************************
// Units IMPLEMENTATION
// **************************************************************
/**
* Sets the units these coalescent intervals are
* measured in.
*/
public final void setUnits(Type u)
{
demographicFunction.setUnits(u);
}
/**
* Returns the units these coalescent intervals are
* measured in.
*/
public final Type getUnits()
{
return demographicFunction.getUnits();
}
public DemographicFunction getDemographicFunction(){
return demographicFunction;
}
public IntervalList getIntervals(){
return intervals;
}
/** The demographic function. */
DemographicFunction demographicFunction = null;
/** The intervals. */
IntervalList intervals = null;
}