package test.dr.evomodel.treelikelihood; import dr.evolution.alignment.SitePatterns; import dr.evolution.datatype.Nucleotides; import dr.evolution.tree.SimpleNode; import dr.evolution.tree.SimpleTree; import dr.evolution.tree.Tree; import dr.evolution.util.Units; import dr.evomodel.tree.TreeModel; import dr.math.matrixAlgebra.Vector; import java.util.ArrayList; import java.util.List; /** * @author Marc A. Suchard */ public class NormalizedSequenceLikelihoodTest extends SequenceLikelihoodTest { public NormalizedSequenceLikelihoodTest(String name) { super(name); } public void setUp() throws Exception { format.setMaximumFractionDigits(5); int numTaxa = TWO_TAXON_SEQUENCE[0].length; createAlignmentWithAllUniquePatterns(TWO_TAXON_SEQUENCE, Nucleotides.INSTANCE); treeModel = createSillyTreeModel(numTaxa); } protected TreeModel createSillyTreeModel(int numTaxa) { SimpleNode[] nodes = new SimpleNode[2 * numTaxa - 1]; for (int n = 0; n < 2 * numTaxa - 1; n++) { nodes[n] = new SimpleNode(); } nodes[0].setTaxon(taxa[0]); nodes[1].setTaxon(taxa[1]); nodes[2].setHeight(1); nodes[2].addChild(nodes[0]); nodes[2].addChild(nodes[1]); SimpleNode root = nodes[2]; if (numTaxa == 3) { nodes[3].setTaxon(taxa[2]); nodes[4].setHeight(2); nodes[4].addChild(nodes[2]); nodes[4].addChild(nodes[3]); root = nodes[4]; } Tree tree = new SimpleTree(root); tree.setUnits(Units.Type.YEARS); return new TreeModel(tree); //treeModel } public void testAllPossibleAlignments() { SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true); tryAllPossibleAlignments(3, patterns); } private void tryAllPossibleAlignments(int length, SitePatterns patterns) { System.out.println("Trying all possible " + length + " site alignments"); double[] patternLogLikeihoods = computeSitePatternLikelihoods(patterns); System.out.println("Site logLikelihoods: " + new Vector(patternLogLikeihoods)); double total = 0.0; List<Double> allAlignmentLogProbabilities = new ArrayList<Double>(); recursivelyComputeAlignmentLikelihood(allAlignmentLogProbabilities, patternLogLikeihoods, length, 0, 0); System.out.println("Total possible alignments: " + allAlignmentLogProbabilities.size()); for (Double x : allAlignmentLogProbabilities) { total += Math.exp(x); } System.out.println("Total probability = " + total); assertEquals("uncorrected", 1.0, total, tolerance); } private void recursivelyComputeAlignmentLikelihood(List<Double> finalLogProbabilities, double[] patternLogLikelihoods, int alignmentLength, int level, double logProbability) { if (level < alignmentLength) { for (int i = 0; i < patternLogLikelihoods.length; i++) { double thisLogProb = logProbability + patternLogLikelihoods[i]; recursivelyComputeAlignmentLikelihood(finalLogProbabilities, patternLogLikelihoods, alignmentLength, level + 1, thisLogProb); } } else { finalLogProbabilities.add(logProbability); } } protected static final String[][] TWO_TAXON_SEQUENCE = {{"human", "chimp", "gorrila"}, {"","", ""}}; }