package dr.evomodel.speciation; import dr.evolution.util.Units; import dr.inference.model.Likelihood; /** * Computes coalescent log-likelihood of a set of gene trees embedded inside a * allopolyploid species network. * * @author Graham Jones * Date: 03/06/2011 */ /* * AlloppMSCoalescent, AlloppSpeciesNetworkModel, AlloppLeggedTree, * AlloppDiploidHistory, and AlloppSpeciesBindings collaborate closely. * * AlloppSpeciesNetworkModel contains an AlloppDiploidHistory and an array of * AlloppLeggedTree's for representing the tetraploid trees. It also contains a * multiply labelled tree as an alternative representation. * * AlloppSpeciesBindings represents the species-individual-genome structure * of the data, and contains the gene trees. * * This module, AlloppMSCoalescent, uses AlloppSpeciesNetworkModel and * AlloppSpeciesBindings to calculate P(g_i|S) = probability that gene tree * g_i fits into network S. * * apsp.geneTreeFitsInNetwork calls * FitsInNetwork (inside a gene tree) calls * subTreeFitsInNetwork (recursive, inside a gene tree) calls * CoalescenceIsCompatible to test a single coalescence against mullab tree * * apsp.geneTreeLogLikelihood calls * TreeLogLikelihood (inside a gene tree) calls * clearCoalescences (in mullab tree) and * recordCoalescence (in mullab tree) and * recordLineageCounts (in mullab tree) and * to set up data in mullab tree nodes, then * mullabTreeLogLikelihood which calls * mullabSubTreeLogLikelihood (recursive, in mullab tree) calls * limbLogLike calls * limbLinPopIntegral * */ public class AlloppMSCoalescent extends Likelihood.Abstract implements Units { private final AlloppSpeciesNetworkModel asnetwork; private final AlloppSpeciesBindings apsp; public AlloppMSCoalescent(AlloppSpeciesBindings apspecies, AlloppSpeciesNetworkModel apspnetwork) { super(apspnetwork); apsp = apspecies; asnetwork = apspnetwork; asnetwork.addModelListener(this); apsp.addModelListeners(this); } @Override protected double calculateLogLikelihood() { for (int i = 0; i < apsp.numberOfGeneTrees(); i++) { if (!apsp.geneTreeFitsInNetwork(i, asnetwork)) { return Double.NEGATIVE_INFINITY; } } // grjtodo-oneday JH has compatible flags for efficiency. I'm checking // every time. double logl = 0; for(int i = 0; i < apsp.numberOfGeneTrees(); i++) { final double v = apsp.geneTreeLogLikelihood(i, asnetwork); assert ! Double.isNaN(v); logl += v; } return logl; } @Override protected boolean getLikelihoodKnown() { return false; } public Type getUnits() { return asnetwork.getUnits(); } public void setUnits(Type units) { // TODO Auto-generated method stub // one day may allow units other than substitutions } }