first stab at integrated prior; works somewhat

src/dr/app/beast/development_parsers.properties

@@ -125,6 +125,8 @@ dr.evomodelxml.coalescent.PowerLawGrowthModelParser
 dr.evomodelxml.coalescent.PeakAndDeclineModelParser
 dr.evomodelxml.coalescent.AsymptoticGrowthModelParser
 
+# INTEGRATED COALESCENT 
+dr.evomodelxml.coalescent.IGCoalescentLikelihoodParser
 
 # UNIFORM INTERNAL NODE HEIGHT PRIOR
 dr.evomodelxml.operators.FunkyPriorMixerOperatorParser

src/dr/evomodel/coalescent/IGCoalescentLikelihood.java

@@ -0,0 +1,137 @@
+/*
+ * CoalescentLikelihood.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.evomodel.coalescent;
+
+import dr.evolution.coalescent.DemographicFunction;
+import dr.evolution.coalescent.IntervalList;
+import dr.evolution.coalescent.IntervalType;
+import dr.evolution.tree.Tree;
+import dr.evolution.tree.TreeUtils;
+import dr.evolution.util.TaxonList;
+import dr.evolution.util.Units;
+import dr.evomodelxml.coalescent.CoalescentLikelihoodParser;
+import dr.math.Binomial;
+import dr.math.LogTricks;
+
+import java.util.List;
+
+
+/**
+ * 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: NewCoalescentLikelihood.java,v 1.6 2005/05/24 20:25:57 rambaut Exp $
+ *
+ * @author Andrew Rambaut
+ * @author Alexei Drummond
+ * @author Luiz Max Carvalho
+ */
+public final class IGCoalescentLikelihood extends AbstractCoalescentLikelihood implements Units {
+
+
+
+	// PUBLIC STUFF
+	public IGCoalescentLikelihood(Tree tree,
+			TaxonList includeSubtree,
+			List<TaxonList> excludeSubtrees,
+			double alpha, double beta) throws TreeUtils.MissingTaxonException {
+
+		super(CoalescentLikelihoodParser.COALESCENT_LIKELIHOOD, tree, includeSubtree, excludeSubtrees);
+
+		this.alpha = alpha;
+		this.beta = beta;
+	}
+
+
+	/**
+	 * Calculates the log likelihood of this set of coalescent intervals,
+	 * given doubles alpha and beta
+	 */
+	public double calculateLogLikelihood() {
+		double lnL =  calculateLogLikelihood(getIntervals(), alpha, beta);
+		return lnL;
+	}
+
+	public static double calculateLogLikelihood(IntervalList intervals, double alpha, double beta) {
+
+		double logL = 0.0;
+		double ldens = 0.0;		
+		final int n = intervals.getIntervalCount();
+		for (int i = 0; i < n; i++) {
+
+			final double duration = intervals.getInterval(i);
+			final int lineageCount = intervals.getLineageCount(i);
+			final double kChoose2 = lineageCount*(lineageCount-1.0)/2.0;
+
+			if (intervals.getIntervalType(i) == IntervalType.COALESCENT) {
+				ldens = marginalLogDensity(duration, kChoose2, alpha, beta);
+				logL +=  ldens;
+			}else {
+				// sampling interval
+				if(lineageCount == 1.0 || duration == 0.0) {
+					ldens = 0.0;
+					logL += ldens;
+				}else {
+					ldens = invGammaNoCoalescentlogProb(duration, kChoose2, alpha, beta);
+					logL += ldens;
+				}
+			}
+			//		System.err.println("i=" + i + " k=" + lineageCount + " delta= " + duration + " isCoal=" + intervals.getIntervalType(i) +  " ldens= " + ldens + "\n");
+		}
+//				System.err.println("logLik=" + logL + "\n");
+		return logL;
+	}
+
+	public static double marginalLogDensity(double t, double coeff, double alpha, double beta) {
+		double lnum = Math.log(coeff) + Math.log(alpha) + alpha * Math.log(beta);
+		double ldenom = (alpha + 1) * Math.log(coeff * t + beta);
+		return(lnum - ldenom);
+	}
+
+	public static double invGammaNoCoalescentlogProb(double x, double coeff, double alpha, double beta) {
+		//Gives (log) Pr(no coalescence | time = x, coeff, alpha, beta)
+		// Pr(no coal | time =x, coeff, Ne) = exp(-x * coeff/Ne). Then integrate against prior on Ne.
+		double logL = alpha * ( Math.log(beta)-Math.log(coeff*x + beta) );
+		return logL;
+	}
+
+	@Override
+	public Type getUnits() {
+		// TODO Auto-generated method stub
+		return null;
+	}
+
+
+	@Override
+	public void setUnits(Type units) {
+		// TODO Auto-generated method stub			
+	}
+
+	// PRIVATE STUFF 
+	private double alpha;
+	private double beta;
+}

src/dr/evomodelxml/coalescent/IGCoalescentLikelihoodParser.java

@@ -0,0 +1,165 @@
+/*
+ * IGCoalescentLikelihoodParser.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.evomodelxml.coalescent;
+
+import dr.evolution.tree.TreeUtils;
+import dr.evolution.util.Taxa;
+import dr.evolution.util.TaxonList;
+import dr.evomodel.coalescent.IGCoalescentLikelihood;
+import dr.evomodel.coalescent.DemographicModel;
+import dr.evomodel.coalescent.MultiLociTreeSet;
+import dr.evomodel.coalescent.OldAbstractCoalescentLikelihood;
+import dr.evomodel.tree.TreeModel;
+import dr.xml.*;
+
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ */
+public class IGCoalescentLikelihoodParser extends AbstractXMLObjectParser {
+
+    public static final String IGCOALESCENT_LIKELIHOOD = "coalescentLikelihoodIG";
+    public static final String MODEL = "model";
+    public static final String POPULATION_TREE = "populationTree";
+    public static final String POPULATION_FACTOR = "factor";
+
+    public static final String ALPHA = "alpha";
+    public static final String BETA = "beta";
+
+    public static final String INCLUDE = "include";
+    public static final String EXCLUDE = "exclude";
+
+    public String getParserName() {
+        return IGCOALESCENT_LIKELIHOOD;
+    }
+
+    public Object parseXMLObject(XMLObject xo) throws XMLParseException {
+
+        XMLObject cxo = xo.getChild(MODEL);
+        DemographicModel demoModel = null;
+
+        double alpha = xo.getDoubleAttribute(ALPHA);
+        double beta = xo.getDoubleAttribute(BETA);
+
+        List<TreeModel> trees = new ArrayList<TreeModel>();
+        List<Double> popFactors = new ArrayList<Double>();
+        MultiLociTreeSet treesSet = demoModel instanceof MultiLociTreeSet ? (MultiLociTreeSet) demoModel : 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));
+                }
+            }
+//                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)) {
+            throw new XMLParseException("Incorrectly constructed likelihood element");
+        }
+
+        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));
+            }
+        }
+
+        if (treeModel != null) {
+            try {
+                return new IGCoalescentLikelihood(treeModel, includeSubtree, excludeSubtrees, alpha, beta);
+            } catch (TreeUtils.MissingTaxonException mte) {
+                throw new XMLParseException("treeModel missing a taxon from taxon list in " + getParserName() + " element");
+            }
+        } else {
+            if (includeSubtree != null || excludeSubtrees.size() > 0) {
+                throw new XMLParseException("Include/Exclude taxa not supported for multi locus sets");
+            }
+            // Use old code for multi locus sets.
+            // This is a little unfortunate but the current code is using AbstractCoalescentLikelihood as
+            // a base - and modifing it will probsbly result in a bigger mess.
+            return new OldAbstractCoalescentLikelihood(treesSet, demoModel);
+        }
+    }
+
+    //************************************************************************
+    // AbstractXMLObjectParser implementation
+    //************************************************************************
+
+    public String getParserDescription() {
+        return "This element represents the likelihood of the tree given the demographic function.";
+    }
+
+    public Class getReturnType() {
+        return IGCoalescentLikelihood.class;
+    }
+
+    public XMLSyntaxRule[] getSyntaxRules() {
+        return rules;
+    }
+
+    private final XMLSyntaxRule[] rules = {
+            new ElementRule(MODEL, new XMLSyntaxRule[]{
+                    new ElementRule(DemographicModel.class)
+            }, "The demographic model which describes the effective population size over time"),
+
+            new ElementRule(POPULATION_TREE, new XMLSyntaxRule[]{
+                    AttributeRule.newDoubleRule(POPULATION_FACTOR, true),
+                    new ElementRule(TreeModel.class)
+            }, "Tree(s) to compute likelihood for", 0, Integer.MAX_VALUE),
+
+            new ElementRule(INCLUDE, new XMLSyntaxRule[]{
+                    new ElementRule(Taxa.class)
+            }, "An optional subset of taxa on which to calculate the likelihood (should be monophyletic)", true),
+
+            new ElementRule(EXCLUDE, new XMLSyntaxRule[]{
+                    new ElementRule(Taxa.class, 1, Integer.MAX_VALUE)
+            }, "One or more subsets of taxa which should be excluded from calculate the likelihood (should be monophyletic)", true)
+    };
+}

src/dr/math/distributions/InverseGammaDistribution.java

@@ -41,13 +41,13 @@ public class InverseGammaDistribution implements Distribution {
     private double shape, scale;
 
     private final double factor;
-    private final double logFacor;
+    private final double logFactor;
 
     public InverseGammaDistribution(double shape, double scale) {
         this.shape = shape;
         this.scale = scale;
         this.factor = Math.pow(scale, shape) / Math.exp(GammaFunction.lnGamma(shape));
-        this.logFacor = shape * Math.log(scale) - GammaFunction.lnGamma(shape);
+        this.logFactor = shape * Math.log(scale) - GammaFunction.lnGamma(shape);
     }
 
     public double getShape() {
@@ -71,7 +71,7 @@ public double pdf(double x) {
     }
 
     public double logPdf(double x) {
-        return logPdf(x, shape, scale, logFacor);
+        return logPdf(x, shape, scale, logFactor);
     }
 
     public double cdf(double x) {
@@ -138,11 +138,11 @@ public static double pdf(double x, double shape, double scale, double factor) {
      * @param scale scale parameter
      * @return log pdf value
      */
-    public static double logPdf(double x, double shape, double scale, double factor) {
+    public static double logPdf(double x, double shape, double scale, double lfactor) {
         if (x <= 0)
             return Double.NEGATIVE_INFINITY;
 
-        return  factor + shape*Math.log(scale) - (shape + 1)*Math.log(x) - (scale/x) - GammaFunction.lnGamma(shape);
+        return  lfactor + shape*Math.log(scale) - (shape + 1)*Math.log(x) - (scale/x) - GammaFunction.lnGamma(shape);
     }
 
     /**