diff --git a/.gitignore b/.gitignore index fcc5f9f..005f59e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,13 @@ # Build output /out +/bin # IntelliJ cruft /.idea /sampled-ancestors.iml + +# Eclipse files +.classpath +.project +/build +/build-lib \ No newline at end of file diff --git a/dist/SA.v2.0.1.zip b/dist/SA.v2.0.1.zip new file mode 100644 index 0000000..549c6c5 Binary files /dev/null and b/dist/SA.v2.0.1.zip differ diff --git a/examples/brachiopods.xml b/examples/brachiopods.xml new file mode 100644 index 0000000..48b093d --- /dev/null +++ b/examples/brachiopods.xml @@ -0,0 +1,236 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + beast.math.distributions.Uniform + beast.math.distributions.Exponential + beast.math.distributions.LogNormalDistributionModel + beast.math.distributions.Normal + beast.math.distributions.Beta + beast.math.distributions.Gamma + beast.math.distributions.LaplaceDistribution + beast.math.distributions.Prior + beast.math.distributions.InverseGamma + beast.math.distributions.OneOnX + + + + +Floweria_anomala=16.3335189128993, +Floweria_arctostriata=17.7354320870945, +Floweria_becraftensis=42.5012110778829, +Floweria_chemungensis=19.6528204575879, +Floweria_chemungensis=10.9193568470655, +Floweria_chemungensis=2.76948535852137, +Floweria_cornucopia=14.9861523577711, +Floweria_crassa=18.5540733357193, +Floweria_deformis=32.0874609937426, +Floweria_iowensis=5.62849222030491, +Floweria_lirella=18.3036514847772, +Floweria_magnacicatrix=12.0972767624771, +Floweria_pandora=32.0440664593131, +Floweria_perversa=30.5905422627926, +Floweria_perversa=19.7067617740249, +Floweria_prava=7.85325374931563, +Floweria_transversalis=9.78931124799419, +Schuchertella_lens=3.90088740836364, +Schuchertella_percha=0.146661693835654, +Xystostrophia_woolworthana=43.2312531741103, +Xystostrophia_umbraculum=23.362639985024, +Xystostrophia_umbraculum=15.6586777545745, +Xystostrophia_umbraculum=5.6815330258105, +Xystostrophia_umbraculum=0 + + + 0.1 + 0.5 + 0.5 + 500.0 + 1.0 + + + + 1.0 + + + + + + 0.0 + 1.0 + + + + + 0.2 + 1.25 + + + + + + + + + + + + + + + + + + + 1.0 + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/beast/app/simulators/SABDSimulator.java b/src/beast/app/simulators/SABDSimulator.java index d981a3d..23f84e1 100644 --- a/src/beast/app/simulators/SABDSimulator.java +++ b/src/beast/app/simulators/SABDSimulator.java @@ -5,7 +5,6 @@ import java.io.*; import java.util.*; -import java.util.concurrent.ExecutionException; /** * diff --git a/src/beast/app/simulators/SABDSkylineTreeSimulator.java b/src/beast/app/simulators/SABDSkylineTreeSimulator.java index 8eee81b..8f84f96 100644 --- a/src/beast/app/simulators/SABDSkylineTreeSimulator.java +++ b/src/beast/app/simulators/SABDSkylineTreeSimulator.java @@ -2,7 +2,6 @@ import beast.core.Description; import beast.evolution.tree.Node; -import beast.evolution.tree.ZeroBranchSANode; import beast.util.Randomizer; import java.io.File; @@ -60,7 +59,7 @@ public SABDSkylineTreeSimulator(double newLambda, double newMu, double newPsi, d public int simulate(PrintStream writer) { //create an initial node (origin of tree) - Node initial = new ZeroBranchSANode(); + Node initial = new Node(); initial.setNr(-1); initial.setHeight(0.0); ArrayList tipNodes = new ArrayList(); // an array of temporary tip nodes sorted by their heights @@ -160,7 +159,7 @@ public int simulate(PrintStream writer) { } //the unique node remains in the array is the root of the sampled tree - Node root = new ZeroBranchSANode(); + Node root = new Node(); for (Node node:children) { root=node; } @@ -214,9 +213,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter switch (typeOfEvent) { case BIRTH: { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(-1); left.setHeight(height); right.setHeight(height); @@ -238,9 +237,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter double remain = Randomizer.nextDouble(); if (r < remain) { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(sampleCount); left.setHeight(height); right.setHeight(height); @@ -434,7 +433,7 @@ private int countSA(Node node){ if (!node.isLeaf()) { return countSA(node.getLeft()) + countSA(node.getRight()); } else { - if (((ZeroBranchSANode)node).isDirectAncestor()) { + if (node.isDirectAncestor()) { return 1; } else return 0; } diff --git a/src/beast/app/simulators/SABDSkylineTreeSimulatorArbitraryRateChanges.java b/src/beast/app/simulators/SABDSkylineTreeSimulatorArbitraryRateChanges.java index d4085c5..535807c 100644 --- a/src/beast/app/simulators/SABDSkylineTreeSimulatorArbitraryRateChanges.java +++ b/src/beast/app/simulators/SABDSkylineTreeSimulatorArbitraryRateChanges.java @@ -2,7 +2,6 @@ import beast.core.Description; import beast.evolution.tree.Node; -import beast.evolution.tree.ZeroBranchSANode; import beast.util.Randomizer; import java.util.*; @@ -65,7 +64,7 @@ public SABDSkylineTreeSimulatorArbitraryRateChanges(double[] newLambda, double[] public int simulate() { //create an initial node (origin of tree) - Node initial = new ZeroBranchSANode(); + Node initial = new Node(); initial.setNr(-1); initial.setHeight(0.0); ArrayList tipNodes = new ArrayList(); // an array of nodes at the previous stage of simulation @@ -171,7 +170,7 @@ public int simulate() { } //the unique node remained in the array is the root of the sampled tree - Node root = new ZeroBranchSANode(); + Node root = new Node(); for (Node node:children) { root=node; } @@ -266,9 +265,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter node.setHeight(height); switch (typeOfEvent) { case BIRTH: { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(-1); left.setHeight(height); right.setHeight(height); @@ -284,9 +283,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter case SAMPLING: { double remain = Randomizer.nextDouble(); if (r[currentInterval] < remain) { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(sampleCount); left.setHeight(height); right.setHeight(height); diff --git a/src/beast/app/simulators/SABDTreeSimulator.java b/src/beast/app/simulators/SABDTreeSimulator.java index 4524429..066566f 100644 --- a/src/beast/app/simulators/SABDTreeSimulator.java +++ b/src/beast/app/simulators/SABDTreeSimulator.java @@ -1,7 +1,6 @@ package beast.app.simulators; import beast.evolution.tree.Node; -import beast.evolution.tree.ZeroBranchSANode; import beast.math.distributions.LogNormalDistributionModel; import beast.util.Randomizer; @@ -77,7 +76,7 @@ public SABDTreeSimulator(double newLambda, double newMu, double newPsi, double n */ public int simulate(PrintStream writer) { //create an initial node (origin of tree) - Node initial = new ZeroBranchSANode(); + Node initial = new Node(); initial.setNr(-1); initial.setHeight(0.0); ArrayList tipNodes = new ArrayList(); // an array of nodes at the previous stage of simulation @@ -145,7 +144,7 @@ public int simulate(PrintStream writer) { } //the unique node remained in the array is the root of the sampled tree - Node root = new ZeroBranchSANode(); + Node root = new Node(); for (Node node:children) { root=node; } @@ -170,7 +169,7 @@ public int simulate(PrintStream writer) { } public int simulateWithRho(PrintStream writer, double[] rootHeight) { - Node initial = new ZeroBranchSANode(); + Node initial = new Node(); initial.setNr(-1); initial.setHeight(0.0); ArrayList tipNodes = new ArrayList(); // an array of nodes at the previous stage of simulation @@ -261,7 +260,7 @@ public int simulateWithRho(PrintStream writer, double[] rootHeight) { } //the unique node remained in the array is the root of the sampled tree - Node root = new ZeroBranchSANode(); + Node root = new Node(); for (Node node:children) { root=node; } @@ -295,7 +294,7 @@ public int simulateWithRho(PrintStream writer, double[] rootHeight) { */ public int simulateWithRhoSamplingTime(PrintStream treeWriter, PrintStream writer, int[] leafCount) { //create an initial node (origin of tree) - Node initial = new ZeroBranchSANode(); + Node initial = new Node(); initial.setNr(-1); initial.setHeight(0.0); ArrayList tipNodes = new ArrayList(); // an array of nodes at the previous stage of simulation @@ -345,7 +344,7 @@ public int simulateWithRhoSamplingTime(PrintStream treeWriter, PrintStream write } //the unique node remained in the array is the root of the sampled tree - Node root = new ZeroBranchSANode(); + Node root = new Node(); for (Node node:children) { root=node; } @@ -426,9 +425,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter node.setHeight(height); switch (typeOfEvent) { case BIRTH: { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(-1); left.setHeight(height); right.setHeight(height); @@ -444,9 +443,9 @@ private ArrayList getNewNodes(Node node, int typeOfEvent, double timeInter case SAMPLING: { double remain = Randomizer.nextDouble(); if (r < remain) { - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(-1); - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(sampleCount); left.setHeight(height); right.setHeight(height); @@ -723,7 +722,7 @@ private int countSA(Node node){ if (!node.isLeaf()) { return countSA(node.getLeft()) + countSA(node.getRight()); } else { - if (((ZeroBranchSANode)node).isDirectAncestor()) { + if (node.isDirectAncestor()) { return 1; } else return 0; } @@ -749,7 +748,7 @@ private void printTraitsWithRhoSamplingTime(Node node, PrintStream writer){ private double printSAWithRhoSamplingTime(Node node, PrintStream writer){ if (node.isLeaf()){ - if (((ZeroBranchSANode)node).isDirectAncestor()) { + if (node.isDirectAncestor()) { writer.println(node.getNr() + "=1"); } else { writer.println(node.getNr() + "=0"); diff --git a/src/beast/app/tools/FullToExtantTreeConverter.java b/src/beast/app/tools/FullToExtantTreeConverter.java index 87ab3c7..2e9ce8b 100644 --- a/src/beast/app/tools/FullToExtantTreeConverter.java +++ b/src/beast/app/tools/FullToExtantTreeConverter.java @@ -64,7 +64,7 @@ private void printConvertedTrees(List trees, String outputFile) throws Exc tree.init(out); out.println(); - tree.log(0, out); + tree.log((long) 0, out); out.println(); for (int i=1; i< trees.size(); i++) { @@ -72,7 +72,7 @@ private void printConvertedTrees(List trees, String outputFile) throws Exc removeFossils(oldTree.getRoot(), oldTree, new ArrayList<>()); numberNodes(oldTree.getRoot(), new int[] {extantTaxa.size()}); tree = new Tree(oldTree.getRoot()); - tree.log(i, out); + tree.log((long) i, out); out.println(); } out.println("End;"); @@ -144,7 +144,7 @@ public void initAndValidate() { @Override public void run() throws Exception { - java.io.File file, file_out; + java.io.File file; String outputFile = ""; double customThreshold = -1.; // ArrayList taxa = null; diff --git a/src/beast/app/tools/NexusImporter.java b/src/beast/app/tools/NexusImporter.java index ed2511d..7531d16 100644 --- a/src/beast/app/tools/NexusImporter.java +++ b/src/beast/app/tools/NexusImporter.java @@ -3,7 +3,6 @@ import java.io.IOException; import java.io.Reader; import java.util.ArrayList; -import java.util.Map; /** * @author Alexandra Gavryushkina diff --git a/src/beast/app/tools/SATreeComparingAnalysis.java b/src/beast/app/tools/SATreeComparingAnalysis.java index 397f0a5..ab0d5ed 100644 --- a/src/beast/app/tools/SATreeComparingAnalysis.java +++ b/src/beast/app/tools/SATreeComparingAnalysis.java @@ -2,8 +2,6 @@ import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.evolution.tree.ZeroBranchSATree; import beast.util.NexusParser; import java.util.ArrayList; @@ -50,7 +48,7 @@ private SATreeComparingAnalysis.TreeSummary[] makeTreeSummaryForAllTrees(List dAPattern = new ArrayList(); for (int i=tree.getLeafNodeCount(); i< tree.getNodeCount(); i++){ - if (((ZeroBranchSANode)tree.getNode(i)).isFake()) { + if (tree.getNode(i).isFake()) { int descendantsCount = tree.getNode(i).getLeafNodeCount() - 1; if (descendantsCount > 0) { for(int j=dAPattern.size(); j dAPattern = new ArrayList(); for (int i=tree.getLeafNodeCount(); i< tree.getNodeCount(); i++){ - if (((ZeroBranchSANode)tree.getNode(i)).isFake()) { + if (tree.getNode(i).isFake()) { int descendantsCount = tree.getNode(i).getLeafNodeCount() - 1; if (descendantsCount > 0) { for(int j=dAPattern.size(); j a = new ArrayList (Arrays.asList(new Integer[]{1, 1})); diff --git a/src/beast/app/tools/SATreeToBinaryTreeConverter.java b/src/beast/app/tools/SATreeToBinaryTreeConverter.java index f269405..e5cd85f 100644 --- a/src/beast/app/tools/SATreeToBinaryTreeConverter.java +++ b/src/beast/app/tools/SATreeToBinaryTreeConverter.java @@ -1,6 +1,5 @@ package beast.app.tools; -import beast.evolution.tree.ZeroBranchSANode; import beast.util.SANexusParser; import java.io.BufferedOutputStream; @@ -19,14 +18,14 @@ * **/ public class SATreeToBinaryTreeConverter { - public static void perform(SampledAncestorTreeTrace trace, boolean useNumbers, String outputFile) throws Exception { + public static void perform(SANexusParser trace, boolean useNumbers, String outputFile) throws Exception { PrintStream out = new PrintStream(new BufferedOutputStream(new FileOutputStream(outputFile))); - trace.beastTrees.get(0).init(out); + trace.trees.get(0).init(out); out.println(); - for (int i=0; i< trace.beastTrees.size(); i++) { + for (int i=0; i< trace.trees.size(); i++) { out.print("tree STATE_" + i + " = "); - out.print(((ZeroBranchSANode)trace.beastTrees.get(i).getRoot()).toSortedNewickWithZeroBranches(new int[]{0})); + out.print(trace.trees.get(i).getRoot().toSortedNewick(new int[]{0}, false)); out.println(";"); } out.println("End;"); @@ -66,8 +65,7 @@ public static void main(String[] args) throws IOException, Exception { reader = new FileReader(file); SANexusParser parser = new SANexusParser(); parser.parseFile(file); - SampledAncestorTreeTrace trace = new SampledAncestorTreeTrace(parser); - perform(trace, useNumbers, outputFile); + perform(parser, useNumbers, outputFile); } catch (IOException e) { // diff --git a/src/beast/app/tools/SampledAncestorTreeAnalysis.java b/src/beast/app/tools/SampledAncestorTreeAnalysis.java deleted file mode 100644 index 10b550c..0000000 --- a/src/beast/app/tools/SampledAncestorTreeAnalysis.java +++ /dev/null @@ -1,446 +0,0 @@ -package beast.app.tools; - -import beast.core.util.ESS; -import beast.evolution.tree.Node; -import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.util.FrequencySet; -import beast.util.Randomizer; - -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collections; -import java.util.Comparator; - -/** - * @author Alexandra Gavryushkina - */ -@Deprecated -public class SampledAncestorTreeAnalysis { //TODO implement burn in - - SampledAncestorTreeTrace trace; - int percentCredSet; - int count; - - FrequencySet pairs = new FrequencySet(); - - public SampledAncestorTreeAnalysis(SampledAncestorTreeTrace newTrace, int percent) { - trace = newTrace; - percentCredSet = percent; - } - - /** - * @parameter useNumbers is true if the indices of sampled nodes should be use as labels - * and taxa names are used otherwise - */ - public void perform(boolean useNumbers) throws Exception { - countClades(true, true); -// countSampledAncestors(true); -// countSAFrequencies(true, false, 0.445); -// printTreeHeights(); -// removeFossils(); -// countTopologies(true); - } - - - public void countTreesWithDClades() throws Exception { - - int dCladeCount = 0; - Tree tree; - - for (int i =0; i < trace.treeCount; i++) { - tree = trace.beastTrees.get(i); - int j; - for (j=0; j clades = new FrequencySet(); - ArrayList tmp = new ArrayList(); - - int dCladeCount = 0; - - for (int i=0; i < trace.treeCount; i++) { - Tree tree = trace.beastTrees.get(i); - ArrayList dClades = extractAllDClades(tree.getRoot(),zeroBranchTrees); - tmp.addAll(dClades); - for (int j=0; j countSAFrequencies(boolean print, boolean useRanking, double cutoff) { - FrequencySet sampledAncestors = new FrequencySet(); - ArrayList tmp = new ArrayList(); - ArrayList predictedSA = new ArrayList(); - - int burnIn = trace.treeCount/10; - - for (int i=burnIn; i < trace.treeCount; i++) { - Tree tree = trace.beastTrees.get(i); - tmp.addAll(listSA(tree, useRanking)); - } - - for (int i=0; i < tmp.size(); i++) { - sampledAncestors.add(tmp.get(i)); - } - - if (print) { - System.out.println("There are " + trace.treeCount + " trees in the file. The first " + burnIn + " are removed as a burn-in. " + - "Then sampled ancestors are counted"); - System.out.println(); - System.out.println("Count \t Percent \t SA"); - System.out.println("flag1"); - for (int i =0; i < sampledAncestors.size(); i++) { - double percent = (double) (sampledAncestors.getFrequency(i) * 100)/(trace.treeCount - burnIn); - if (percent >= cutoff*100) { - predictedSA.add(sampledAncestors.get(i)); - } - System.out.format("%-10d %-10.2f", sampledAncestors.getFrequency(i), percent); - System.out.println(sampledAncestors.get(i)); - } - System.out.println("flag2"); - } - - int trueSACount; - boolean falseSADetected; - Tree randomTree; - int randomIndex; - do { - randomIndex = Randomizer.nextInt(trace.beastTrees.size()); - randomTree = trace.beastTrees.get(randomIndex); - trueSACount=0; - falseSADetected = false; - for (int i=0; i topologies = new FrequencySet(); - String[] trees; - if (useNumbers) { - trees = trace.shortTrees; - } else { - trees = trace.labeledTrees; - } - - for (int i=0; i < trees.length; i++) { - topologies.add(trees[i]); - } - - Double[] numericTrees = new Double[trees.length]; - - Arrays.fill(numericTrees, 0.0); - - for (int i=0; i< topologies.size(); i++) { - for (int j=0; j < trees.length; j++) { - if (trees[j].equals(topologies.get(i))) { - numericTrees[j] = (double)i+1; - } - } - } - - int nSampleInterval = 100; - // calc effective sample size - double ACT = ESS.ACT(numericTrees, nSampleInterval); - double ESS = numericTrees.length / (ACT / nSampleInterval); - - System.out.println("The number of trees in the file = " + trees.length + "."); - //System.out.println("ACT = " + ACT); - System.out.println("ESS = " + ESS); - System.out.println(); - - double sumPercentage = 0; - System.out.println(percentCredSet + "% credible set: "); - System.out.println(); - System.out.println("Count \t Percent \t Topology"); - System.out.println(); - int i; - - for (i =0; i < topologies.size(); i++) - if (sumPercentage < percentCredSet) { - double percent = (double) (topologies.getFrequency(i) * 100)/(trees.length); - System.out.format("%-10d %-10f ", topologies.getFrequency(i), percent); - System.out.println(topologies.get(i)); - sumPercentage += percent; - } else { - break; - } - System.out.println(); - System.out.println("Total \t" + Math.round(sumPercentage) + "%"); - System.out.println(); - System.out.println(percentCredSet + "% credible set has " + i + " trees."); - System.out.println(); - - /*System.out.println("Numeric tree representation"); - for(int j =0; j < 50; j++) - System.out.print(j + ", "); - System.out.println(); - for (int j =0; j < 50; j++) - System.out.print(numericTrees[j] + ", "); */ - - } - - private ArrayList listNodesUnder(Node node) { - - ArrayList tmp = new ArrayList(); - if (node.getChildCount() < 2) - tmp.add(node.getNr()+1); - if (node.getLeft() != null) - tmp.addAll(listNodesUnder(node.getLeft())); - if (node.getRight() != null) - tmp.addAll(listNodesUnder(node.getRight())); - return tmp; - } - - private ArrayList listNodesUnder(Node node, boolean useID) { - ArrayList tmp = new ArrayList(); - if (!node.isLeaf()) { - for (Node child : node.getChildren()) { - tmp.addAll(listNodesUnder(child, useID)); - } - } else tmp.add(node.getID()); - Collections.sort(tmp); - return tmp; - } - - private String extractDClade(Node node, boolean zeroBranchTrees) { - String tmp = new String(); - if (zeroBranchTrees) { - String ancestor = node.getID(); - tmp += ancestor + '<'; - if (((ZeroBranchSANode)node).isDirectAncestor()) { - ArrayList descendants = listNodesUnder(node.getParent(), true); - tmp+= descendants; - for (String des:descendants) { - if (!des.equals(ancestor)) { - pairs.add(ancestor + "<" + des); - } - } - } - } else { - String ancestor = Integer.toString(node.getNr() + 1); - tmp += ancestor + '<'; - if (node.getChildCount() == 1) { - Integer[] descendants = listNodesUnder(node.getLeft()).toArray(new Integer[0]); - Arrays.sort(descendants); - for (int i=0; i < descendants.length; i++){ - String pair = ancestor + '<' + descendants[i]; - pairs.add(pair); - } - tmp += Arrays.toString(descendants); - } - } - - - return tmp; - } - - public ArrayList extractAllDClades(Node node, boolean zeroBranchTrees) { - ArrayList tmp = new ArrayList(); - - if ((!zeroBranchTrees && node.getChildCount() < 2) || (zeroBranchTrees && node.isLeaf())) - tmp.add(extractDClade(node, zeroBranchTrees)); - if (node.getLeft() != null) - tmp.addAll(extractAllDClades(node.getLeft(), zeroBranchTrees)); - if (node.getRight() != null) - tmp.addAll(extractAllDClades(node.getRight(), zeroBranchTrees)); - - return tmp; - } - - /** - * retern the list of sampled ancestors, WARNING works only for zeroBranchTrees - * @param tree - * @return - */ - public ArrayList listSA(Tree tree, boolean useRanking){ - ArrayList sampledAncestors = new ArrayList(); - count = 0; - for (int i=0; i nodes = new ArrayList(tree.getExternalNodes()); - Comparator comp = new NodeComparator(); - - Collections.sort(nodes, comp); - return nodes.size() - nodes.indexOf(node); - } - - public void printTreeHeights(){ - Tree tree; - - System.out.print("heights <- c("); - for (int i =0; i < trace.treeCount-1; i++) { - tree = trace.beastTrees.get(i); - System.out.print(tree.getRoot().getHeight() + ","); - } - tree = trace.beastTrees.get(trace.treeCount-1); - System.out.println(tree.getRoot().getHeight() + ")"); - - System.out.print("lengths <- c("); - for (int i =0; i < trace.treeCount-1; i++) { - tree = trace.beastTrees.get(i); - double length = 0; - for (int j=0; j< tree.getNodeCount(); j++){ - length += tree.getNode(j).getLength(); - } - - System.out.print(length + ","); - } - tree = trace.beastTrees.get(trace.treeCount-1); - System.out.println(tree.getRoot().getHeight() + ")"); - } - - public void removeFossils(){ - Tree tree; - - for (int i =0; i < trace.treeCount; i++) { - tree = trace.beastTrees.get(i); - for (int j=0; j0.0000000005) || (fake.getRight().isLeaf() && fake.getRight().getHeight() > 0.0000000005))) { - Node parent = fake.getParent(); - Node otherChild= (fake.getLeft().isLeaf() && fake.getLeft().getHeight()>0.0000000005)?fake.getRight():fake.getLeft(); - parent.removeChild(fake); - parent.addChild(otherChild); - } - } - System.out.println("tree STATE_"+ i*1000 + " = "+tree.getRoot().toSortedNewick(new int[]{0}, false) + ";"); - } - - } - - - - -} \ No newline at end of file diff --git a/src/beast/app/tools/SampledAncestorTreeTrace.java b/src/beast/app/tools/SampledAncestorTreeTrace.java deleted file mode 100644 index 03ebeba..0000000 --- a/src/beast/app/tools/SampledAncestorTreeTrace.java +++ /dev/null @@ -1,213 +0,0 @@ -package beast.app.tools; - -import beast.evolution.tree.Node; -import beast.evolution.tree.Tree; -import beast.util.NexusParser; - -import java.util.Arrays; -import java.util.HashSet; -import java.util.List; -import java.util.Map; - -/** - * @author Alexandra Gavryushkina - */ -@Deprecated -public class SampledAncestorTreeTrace { - - public List beastTrees; - - public int treeCount; - - /** - * an array of short newick representation of trees from the parser - */ - public String[] trees; - - /** - * an array of newick strings with no lengths - */ - public String[] shortTrees; - - public int labelCount; - - public String[] labeledTrees; - - // private String[] rankedTrees; - - public HashSet labelSet; - - public SampledAncestorTreeTrace(NexusParser parser) throws Exception { - - beastTrees = parser.trees; - - String[] newick = new String[parser.trees.size()]; - for (int i =0 ; i < newick.length; i++) { - newick[i] = beastTrees.get(i).getRoot().toSortedNewick(new int[] {1}, false); - } - - trees = newick; - treeCount = trees.length; - - shortTrees = new String[treeCount]; - - for (int i =0; i < treeCount; i++) { - shortTrees[i] = convertToShortTree(trees[i]); - } - - labeledTrees = new String[treeCount]; - - //labelCount = parser.translationMap.keySet().size(); - -// Integer[] tmp = new Integer[labelCount]; -// for (int i=0; i< labelCount; i++) -// tmp[i] = i+1; -// -// labelSet = new HashSet(Arrays.asList(tmp)); - -// for (int i=0; i < treeCount; i++){ -// labeledTrees[i] = convertToLabeledTree(shortTrees[i], null); //parser.translationMap); -// } - - } - - public SampledAncestorTreeTrace(int newLabelCount) { - - labelCount = newLabelCount; - - Integer[] tmp = new Integer[labelCount]; - for (int i=0; i< labelCount; i++) - tmp[i] = i+1; - - labelSet = new HashSet(Arrays.asList(tmp)); - } - - public String convertToLabeledTree(String strIn, Map translationMap) { - - StringBuilder buf = new StringBuilder(); - boolean readingInt = false; - int begin=0; - int end=1; - String sInt; - - for (int i=0; i < strIn.length(); i++) { - if (!readingInt) { - if (strIn.charAt(i) == '(' || strIn.charAt(i) == ')' || strIn.charAt(i) == ',' ) { - buf.append(strIn.charAt(i)); - } else { - begin = i; - end = i+1; - readingInt = true; - } - } else { - if (strIn.charAt(i) == '(' || strIn.charAt(i) == ')' || strIn.charAt(i) == ',' ) { - String key = strIn.substring(begin, end); - buf.append(translationMap.get(key)); - buf.append(strIn.charAt(i)); - readingInt = false; - } else { - end++; - } - } - } - if(readingInt) { - String key = strIn.substring(begin, end); - buf.append(translationMap.get(key)); - } - - return buf.toString(); - } - - - public static String convertToShortTree(String strIn) { - - StringBuilder buf = new StringBuilder(); - boolean skipping = false; - boolean inComment = false; - - for (int i=0; i < strIn.length(); i++) { - if (!skipping && !inComment) { - - switch (strIn.charAt(i)) { - case ':': - if (strIn.charAt(i+1) != '(') { - skipping = true; - i += 1; - } else { - buf.append(strIn.charAt(i)); - } - break; - case '[': - inComment = true; - break; - default: - buf.append(strIn.charAt(i)); - } - } else { - if (!inComment && (strIn.charAt(i) == ',' || strIn.charAt(i) == ')')) { - buf.append(strIn.charAt(i)); - skipping = false; - } - if (inComment && strIn.charAt(i) == ']') { - inComment = false; - } - } - } - - return buf.toString(); - } - - private void shiftHeights(Node node, double shift) { - if (node.getLeft() != null) { - double tmp = node.getLeft().getHeight(); - node.getLeft().setHeight(tmp + shift); - shiftHeights(node.getLeft(), shift); - } - if (node.getRight() != null) { - double tmp = node.getRight().getHeight(); - node.getRight().setHeight(tmp + shift); - shiftHeights(node.getRight(), shift); - } - } - - private double readLength(String strTree, int[] currentPosition){ - - String tmp = new String(); - int start = currentPosition[0]; - - currentPosition[0] = strTree.length(); - - for (int i = start; i< strTree.length(); i++) { - if (strTree.charAt(i) != ')' && strTree.charAt(i) != ',') - tmp += strTree.charAt(i); - else { - currentPosition[0] = i; - break; - } - } - - return Double.parseDouble(tmp); - - } - - public void changeTaxaNumbers(Tree tree, Map newTranslationMap) { - for (int i=0; i < tree.getLeafNodeCount(); i++) { - String taxon = newTranslationMap.get(String.valueOf(i)); - int oldIndex = i; - for (int j=i; j < tree.getLeafNodeCount(); j++) { - if (tree.getNode(j).getID().equals(taxon)) { - oldIndex = j; - break; - } - } - if (i != oldIndex) { - tree.getNode(i).setNr(oldIndex); - tree.getNode(oldIndex).setNr(i); - // tree.initArray TODO implement this, but an array of nodes is private so... - } - - } - } - -} - diff --git a/src/beast/evolution/operators/ExchangeForZeroBranchSampledAncestorTrees.java b/src/beast/evolution/operators/ExchangeForZeroBranchSampledAncestorTrees.java deleted file mode 100644 index face304..0000000 --- a/src/beast/evolution/operators/ExchangeForZeroBranchSampledAncestorTrees.java +++ /dev/null @@ -1,135 +0,0 @@ -package beast.evolution.operators; - -/** - * @Alexandra Gavryushkina - */ - - -import beast.core.Description; -import beast.core.Input; -import beast.evolution.tree.Node; -import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.util.Randomizer; -/** - * @deprecated Use SAExchange instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -@Description("Implement Narrow and Wide Exchange for sampled ancestor trees." + - "Narrow move chooses a random internal node (not a fake node) with two non-leaf children." + - "Then it takes the older child of this node and exchange one of its children (or just a child" + - "if there is only one) with the younger child. Wide remains the same as for regular trees.") -public class ExchangeForZeroBranchSampledAncestorTrees extends Exchange { - - public double narrow(final Tree tree) { - - final int nodeCount = tree.getNodeCount(); - - //make sure that there are at least two distinct non-root nodes which are not direct ancestors. - if (nodeCount == 3 && ((ZeroBranchSANode)tree.getRoot()).isFake()) { - return Double.NEGATIVE_INFINITY; - } - - Node i; - do { - i = tree.getNode(Randomizer.nextInt(nodeCount)); - } while (i.isRoot() || i.getParent().isRoot() || ((ZeroBranchSANode)i).isDirectAncestor()); - - final Node iParent = i.getParent(); - final Node iGrandParent = iParent.getParent(); - Node iUncle = iGrandParent.getLeft(); - if (iUncle.getNr() == iParent.getNr()) { - iUncle = iGrandParent.getRight(); - //assert (iUncle.getNr() != iParent.getNr()); - } - //assert iUncle == getOtherChild(iGrandParent, iParent); - - //assert i.getHeight() <= iGrandParent.getHeight(); - - if (iUncle.getHeight() < iParent.getHeight()) { - exchangeNodes(i, iUncle, iParent, iGrandParent); - return 0.0; - } else { - // Couldn't find valid narrow move on this beast.tree!! - return Double.NEGATIVE_INFINITY; - } - -// final int nInternalNodes = tree.getInternalNodeCount(); //TODO look if I can implement this more efficient code for SA trees -// final int leafNodeCount = tree.getLeafNodeCount(); -// // make sure that the tree has at least two internal nodes -// if (nInternalNodes <= 1 ) { -// return Double.NEGATIVE_INFINITY; -// } -// -// //choose one of internal nodes that has at least one non-leaf node -// //(there is always at least one such node as long as the tree has at least 2 internal nodes) -// Node iGrandParent; -// do { -// iGrandParent = tree.getNode(leafNodeCount + Randomizer.nextInt(nInternalNodes)); -// } while (iGrandParent.getLeft().isLeaf() && iGrandParent.getRight().isLeaf()); -// -// Node iParent = iGrandParent.getLeft(); -// Node iUncle = iGrandParent.getRight(); -// if (iParent.getHeight() < iUncle.getHeight()) { -// iParent = iGrandParent.getRight(); -// iUncle = iGrandParent.getLeft(); -// } -// -// if( iParent.isLeaf() ) { -// return Double.NEGATIVE_INFINITY; -// } -// -// final Node i; -// -// if (iParent.isFake()) { -// if (iParent.getLeft().isDirectAncestor()) { -// i = iParent.getRight(); -// } else i = iParent.getLeft(); -// } else { -// i = (Randomizer.nextBoolean() ? iParent.getLeft() : iParent.getRight()); -// } -// -// exchangeNodes(i, iUncle, iParent, iGrandParent); -// -// return 0.0; - - } - - /** - * @param tree - */ - public double wide(final Tree tree) { - - final int nodeCount = tree.getNodeCount(); - - //make sure that there are at least two distinct non-root nodes which are not direct ancestors. - if (nodeCount == 3 && ((ZeroBranchSANode)tree.getRoot()).isFake()) { - return Double.NEGATIVE_INFINITY; - } - - Node i, j, iP, jP; - do { - i = tree.getNode(Randomizer.nextInt(nodeCount)); - } while (i.isRoot() || ((ZeroBranchSANode)i).isDirectAncestor()); - - do { - j = tree.getNode(Randomizer.nextInt(nodeCount)); - } while (j.getNr() == i.getNr() || j.isRoot() || ((ZeroBranchSANode)j).isDirectAncestor()); - - iP = i.getParent(); - jP = j.getParent(); - - if ((iP != jP) && (i != jP) && (j != iP) - && (j.getHeight() < iP.getHeight()) - && (i.getHeight() < jP.getHeight())) { - exchangeNodes(i, j, iP, jP); - - return 0.0; - } - - // Couldn't find valid wide move on this beast.tree! - return Double.NEGATIVE_INFINITY; - } - -} diff --git a/src/beast/evolution/operators/IntUniformWithExclusion.java b/src/beast/evolution/operators/IntUniformWithExclusion.java index 69d58cd..55d29ab 100644 --- a/src/beast/evolution/operators/IntUniformWithExclusion.java +++ b/src/beast/evolution/operators/IntUniformWithExclusion.java @@ -5,7 +5,6 @@ import beast.core.Operator; import beast.core.parameter.IntegerParameter; import beast.core.parameter.Parameter; -import beast.core.parameter.RealParameter; import beast.util.Randomizer; /** diff --git a/src/beast/evolution/operators/JumpToPoint.java b/src/beast/evolution/operators/JumpToPoint.java index 2257201..521065c 100644 --- a/src/beast/evolution/operators/JumpToPoint.java +++ b/src/beast/evolution/operators/JumpToPoint.java @@ -4,7 +4,6 @@ import beast.core.Operator; import beast.core.parameter.RealParameter; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSATree; import beast.util.Randomizer; /** @@ -36,7 +35,7 @@ public double proposal() { double r = rParameter.getValue(); - if (r != 1 && ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() == 0) { + if (r != 1 && tree.getDirectAncestorNodeCount() == 0) { rParameter.setValue(1.); return 0.0; } else if (r == 1) { diff --git a/src/beast/evolution/operators/SAScaleOperator.java b/src/beast/evolution/operators/SAScaleOperator.java index 76f7d74..494bb29 100644 --- a/src/beast/evolution/operators/SAScaleOperator.java +++ b/src/beast/evolution/operators/SAScaleOperator.java @@ -39,15 +39,7 @@ public double proposal() { return Double.NEGATIVE_INFINITY; } - if ((root).isFake() && scaleSNodes) { - Node directAncestor = root.getLeft(); - if (!(directAncestor).isDirectAncestor()) - directAncestor = root.getRight(); - root.setHeight(fNewHeight); - directAncestor.setHeight(fNewHeight); - } else { - root.setHeight(fNewHeight); - } + root.setHeight(fNewHeight); return -Math.log(scale); } else { diff --git a/src/beast/evolution/operators/SampledNodeDateRandomWalker.java b/src/beast/evolution/operators/SampledNodeDateRandomWalker.java index 63d8577..560f55c 100644 --- a/src/beast/evolution/operators/SampledNodeDateRandomWalker.java +++ b/src/beast/evolution/operators/SampledNodeDateRandomWalker.java @@ -1,28 +1,38 @@ package beast.evolution.operators; +import java.util.ArrayList; +import java.util.List; + import beast.core.Description; import beast.core.Input; import beast.evolution.tree.Node; import beast.evolution.tree.SamplingDate; import beast.evolution.tree.Tree; +import beast.evolution.tree.TreeWOffset; import beast.util.Randomizer; -import java.util.ArrayList; -import java.util.Collections; -import java.util.List; - @Description("Randomly select a sampled node and shifts the date of the node within a given window") public class SampledNodeDateRandomWalker extends TipDatesRandomWalker { + + public Input treeWOffsetInput = + new Input("treeWOffset", "Optional fully extinct tree", (TreeWOffset)null); public Input> samplingDatesInput = new Input<>("samplingDates", "List of sampling dates", new ArrayList()); boolean useNodeNumbers; List samplingDateTaxonNames = new ArrayList<>(); - + TreeWOffset combinedTree; @Override public void initAndValidate() { + combinedTree = treeWOffsetInput.get(); + if(combinedTree == null) { + combinedTree = new TreeWOffset(); + combinedTree.setInputValue("tree", treeInput.get()); + combinedTree.initAndValidate(); + } + windowSize = windowSizeInput.get(); useGaussian = useGaussianInput.get(); @@ -58,7 +68,7 @@ public void initAndValidate() { public double proposal() { // randomly select a leaf node - Tree tree = treeInput.get(); + Tree tree = combinedTree.getTree(); Node node; if (useNodeNumbers) { @@ -70,7 +80,7 @@ public double proposal() { node = tree.getNode(taxonIndices[i]); } - double value = node.getHeight(); + double value = combinedTree.getHeightOfNode(node.getNr()); if (value == 0.0) { return Double.NEGATIVE_INFINITY; @@ -96,14 +106,14 @@ public double proposal() { if ((node).isDirectAncestor()) { fake = node.getParent(); - lower = getOtherChild(fake, node).getHeight(); + lower = combinedTree.getHeightOfNode(getOtherChild(fake, node).getNr()); if (fake.getParent() != null) { - upper = fake.getParent().getHeight(); + upper = combinedTree.getHeightOfNode(fake.getParent().getNr()); } else upper = Double.POSITIVE_INFINITY; } else { //lower = Double.NEGATIVE_INFINITY; lower = 0.0; - upper = node.getParent().getHeight(); + upper = combinedTree.getHeightOfNode(node.getParent().getNr()); } if (newValue < lower || newValue > upper) { @@ -116,10 +126,10 @@ public double proposal() { } if (fake != null) { - fake.setHeight(newValue); + combinedTree.setHeightOfNode(fake.getNr(), newValue); } - node.setHeight(newValue); - + combinedTree.setHeightOfNode(node.getNr(), newValue); + return 0.0; } diff --git a/src/beast/evolution/operators/SampledNodeDateRandomWalkerForZeroBranchSATrees.java b/src/beast/evolution/operators/SampledNodeDateRandomWalkerForZeroBranchSATrees.java deleted file mode 100644 index 04b48f7..0000000 --- a/src/beast/evolution/operators/SampledNodeDateRandomWalkerForZeroBranchSATrees.java +++ /dev/null @@ -1,340 +0,0 @@ -package beast.evolution.operators; - -import beast.core.Description; -import beast.core.Distribution; -import beast.core.Input; -import beast.core.util.Log; -import beast.evolution.tree.Node; -import beast.evolution.tree.SamplingDate; -import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.math.distributions.*; -import beast.math.distributions.Uniform; -import beast.util.Randomizer; - -import java.io.*; -import java.util.ArrayList; -import java.util.Collections; -import java.util.List; - -/** - * @deprecated Use SampledNodeDateRandomWalker instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -@Description("Randomly select a sampled node and shifts the date of the node within a given window") -public class SampledNodeDateRandomWalkerForZeroBranchSATrees extends TipDatesRandomWalker { - - public Input> samplingDatesInput = new Input<>("samplingDates", - "List of sampling dates", new ArrayList<>()); - - public Input weightFileInput = new Input<>("weightFile", "tab-delimited file (no header) containing relative weights to apply to each sampling date (and optionally ESS values for sampling dates from previous run.) " + - "If the ESSs are included then they are assumed to be from a previous run using the given weights. They are used to calculate a new set of weights to even out the ESSs. " + - "Each new relative weights is calculated to be proportional to old weight divided by corresponding ESS. " + - "The weight file is only used if a taxon set is also provided. "); - - public Input weightOutFileInput = new Input<>("weightOutFile", "tab-delimited file (no header) containing weights applied to each sampling date", (String)null); - - public Input useWindowSizeWithSamplingDates = new Input<>("useWindowSizeWithSamplingDates", "if true, use window size even when sampling date range is specified.", Boolean.FALSE); - - boolean useNodeNumbers; - List samplingDateTaxonNames = new ArrayList<>(); - - double[] relativeWeights = null; - - @Override - public void initAndValidate() { - windowSize = windowSizeInput.get(); - useGaussian = useGaussianInput.get(); - - for (SamplingDate taxon : samplingDatesInput.get()) { - samplingDateTaxonNames.add(taxon.taxonInput.get().getID()); - } - - // determine taxon set to choose from - if (m_taxonsetInput.get() != null) { - useNodeNumbers = false; - List sTaxaNames = new ArrayList(); - for (String sTaxon : treeInput.get().getTaxaNames()) { - sTaxaNames.add(sTaxon); - } - - List set = m_taxonsetInput.get().asStringList(); - int nNrOfTaxa = set.size(); - taxonIndices = new int[nNrOfTaxa]; - int k = 0; - for (String sTaxon : set) { - int iTaxon = sTaxaNames.indexOf(sTaxon); - if (iTaxon < 0) { - throw new IllegalArgumentException("Cannot find taxon " + sTaxon + " in tree"); - } - taxonIndices[k++] = iTaxon; - } - - if (weightFileInput.get() != null) { - relativeWeights = new double[nNrOfTaxa]; - double[] ess = new double[nNrOfTaxa]; - int count = 0; - double essSum = 0.0; - double prevWeightSum = 0.0; - - double[] prevWeights = new double[nNrOfTaxa]; - - - try { - BufferedReader reader = new BufferedReader(new FileReader(weightFileInput.get())); - - String line = reader.readLine(); - while (line != null ) { - String[] parts = line.split("\t"); - int index = set.indexOf(parts[0]); - if (index < 0) { - Log.warning("taxon '" + parts[0] + "' in weight file '" + weightFileInput.get() + "' not found in taxa set."); - } else { - if (ess[index] > 0) { - throw new RuntimeException("Taxon '" + parts[0] + "' is duplicated in ESS file '" + weightFileInput.get() + "'"); - } - count += 1; - ess[index] = Double.parseDouble(parts[2]); - essSum += ess[index]; - prevWeights[index] = Double.parseDouble(parts[1]); - prevWeightSum += prevWeights[index]; - } - line = reader.readLine(); - } - - if (count < ess.length) { - double meanESS = essSum / count; - double meanPrevWeight = prevWeightSum / count; - Log.warning("Only " + count + " out of " + ess.length + " taxa found in ESS file. The remainder will be set to the mean weight and ESS."); - for (int i = 0; i < ess.length; i++) { - if (ess[i] == 0) { - ess[i] = meanESS; - prevWeights[i] = meanPrevWeight; - } - } - } - - double sum = 0.0; - for (int i = 0; i < relativeWeights.length; i++) { - relativeWeights[i] = prevWeights[i] / ess[i]; - sum += relativeWeights[i]; - } - - PrintWriter writer = null; - if (weightOutFileInput.get() != null) { - writer = new PrintWriter(new FileWriter(weightOutFileInput.get())); - } - - // normalize - for (int i = 0; i < relativeWeights.length; i++) { - relativeWeights[i] /= sum; - if (weightOutFileInput.get() != null) { - - if (writer != null) { - writer.println(set.get(i) + "\t" + relativeWeights[i]); - } else { - Log.info(set.get(i) + "\t" + relativeWeights[i]); - } - } - } - - if (writer != null) { - writer.flush(); - writer.close(); - } - - //cumulative - for (int i = 1; i < relativeWeights.length; i++) { - relativeWeights[i] = relativeWeights[i-1] + relativeWeights[i]; - } - //protect from rounding errors. - relativeWeights[relativeWeights.length-1] = 1.0; - - } catch (FileNotFoundException e) { - Log.warning("Weight file named '" + weightFileInput.get() + "' not found. Defaulting to equal weights."); - relativeWeights = null; - } catch (IOException e) { - Log.warning("IO exception reading file named '" + weightFileInput.get() + "'. Defaulting to equal weights."); - relativeWeights = null; - } - } - - } else { - useNodeNumbers = true; - } - - - // check that all nodes are within bounds. - List nodes = getNodesToOperateOn(); - boolean err = false; - for (Node node : nodes) { - double age = node.getHeight(); - boolean drawFromDistribution = samplingDateTaxonNames.contains(node.getID()); - if (drawFromDistribution) { - SamplingDate taxonSamplingDate = samplingDatesInput.get().get(samplingDateTaxonNames.indexOf(node.getID())); - double lower = taxonSamplingDate.getLower(); - double upper = taxonSamplingDate.getUpper(); - if (age > upper || age < lower) { - err = true; - Log.err("Node " + node.getID() + " has an age (" + age + ") outside the sampling date range (" + lower + ", " + upper + ")."); - } - } - } - if (err) { - throw new RuntimeException("Error: Stopping because nodes found out of sampling date range."); - } - } - - private List getNodesToOperateOn() { - Tree tree = treeInput.get(); - - ArrayList nodeList = new ArrayList<>(); - - if (useNodeNumbers) { - int leafNodeCount = tree.getLeafNodeCount(); - } else { - for (int i = 0; i < taxonIndices.length; i++) { - nodeList.add(tree.getNode(taxonIndices[i])); - } - } - return nodeList; - } - - @Override - public double proposal() { - - // randomly select a leaf node - Tree tree = treeInput.get(); - - Node node; - if (useNodeNumbers) { - int leafNodeCount = tree.getLeafNodeCount(); - int i = Randomizer.nextInt(leafNodeCount); - node = tree.getNode(i); - } else { - - if (relativeWeights != null) { - double X = Randomizer.nextDouble(); - int i = 0; - while (relativeWeights[i] < X) { - i += 1; - } - node = tree.getNode(taxonIndices[i]); - } else { - int i = Randomizer.nextInt(taxonIndices.length); - node = tree.getNode(taxonIndices[i]); - } - } - - double value = node.getHeight(); - - if (value == 0.0) { - return Double.NEGATIVE_INFINITY; - } - double newValue = value; - - boolean drawFromDistribution = (samplingDateTaxonNames.contains(node.getID()) && !useWindowSizeWithSamplingDates.get()); - if (drawFromDistribution) { - SamplingDate taxonSamplingDate = samplingDatesInput.get().get(samplingDateTaxonNames.indexOf(node.getID())); - double range = taxonSamplingDate.getUpper() - taxonSamplingDate.getLower(); - newValue = taxonSamplingDate.getLower() + Randomizer.nextDouble() * range; - } else { - if (useGaussian) { - newValue += Randomizer.nextGaussian() * windowSize; - } else { - newValue += Randomizer.nextDouble() * 2 * windowSize - windowSize; - } - if (samplingDatesInput.get() != null) { - SamplingDate taxonSamplingDate = samplingDatesInput.get().get(samplingDateTaxonNames.indexOf(node.getID())); - if (newValue < taxonSamplingDate.getLower() || newValue > taxonSamplingDate.getUpper()) { - return Double.NEGATIVE_INFINITY; - } - } - } - - Node fake = null; - double lower, upper; - - if (node.isDirectAncestor()) { - fake = node.getParent(); - lower = getOtherChild(fake, node).getHeight(); - if (fake.getParent() != null) { - upper = fake.getParent().getHeight(); - } else upper = Double.POSITIVE_INFINITY; - } else { - //lower = Double.NEGATIVE_INFINITY; - lower = 0.0; - upper = node.getParent().getHeight(); - } - - if (newValue < lower || newValue > upper) { - return Double.NEGATIVE_INFINITY; - } - - if (newValue == value) { - // this saves calculating the posterior - return Double.NEGATIVE_INFINITY; - } - - if (fake != null) { - fake.setHeight(newValue); - } - node.setHeight(newValue); - - if (newValue < 0) { - for (int i=0; i tipNodeHeights= new ArrayList(); - for (int i=0; i m_pScaleSNodes = new Input("scaleSampledNodes", "If it is true then sampled node dates are scaled (default false).", false); - - @Override //WARNING works with bifurcating (exactly 2 children) trees only - // sampled ancestors are assumed to be on zero branches - - public double proposal() { - - final double scale = getScaler(); - final boolean scaleSNodes = m_pScaleSNodes.get(); - - try { - - if (m_bIsTreeScaler) { - Tree tree = treeInput.get(this); - if (rootOnlyInput.get()) { - Node root = tree.getRoot(); - if (((ZeroBranchSANode)root).isFake() && !scaleSNodes) { - return Double.NEGATIVE_INFINITY; - } - double fNewHeight = root.getHeight() * scale; - - //make sure the new height doesn't make a parent younger than a child - double oldestChildHeight; - if (((ZeroBranchSANode)root).isFake()) { - oldestChildHeight = root.getNonDirectAncestorChild().getHeight(); - } else oldestChildHeight = Math.max(root.getLeft().getHeight(), root.getRight().getHeight()); - if (fNewHeight < oldestChildHeight) { - return Double.NEGATIVE_INFINITY; - } - - if (((ZeroBranchSANode)root).isFake() && scaleSNodes) { - Node directAncestor = root.getLeft(); - if (!((ZeroBranchSANode)directAncestor).isDirectAncestor()) - directAncestor = root.getRight(); - root.setHeight(fNewHeight); - directAncestor.setHeight(fNewHeight); - } else { - root.setHeight(fNewHeight); - } - - return -Math.log(scale); - } else { - // scale the beast.tree - final int nScaledDimensions = ((ZeroBranchSATree)tree).scale(scale, scaleSNodes); - return Math.log(scale) * (nScaledDimensions - 2); - } - } - return Double.NEGATIVE_INFINITY; - - } catch (Exception e) { - return Double.NEGATIVE_INFINITY; - } - } - -} diff --git a/src/beast/evolution/operators/TreeDimensionJump.java b/src/beast/evolution/operators/TreeDimensionJump.java index df5c595..13bc856 100644 --- a/src/beast/evolution/operators/TreeDimensionJump.java +++ b/src/beast/evolution/operators/TreeDimensionJump.java @@ -6,8 +6,6 @@ import beast.core.parameter.RealParameter; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.evolution.tree.ZeroBranchSATree; import beast.util.Randomizer; /** @@ -48,7 +46,7 @@ public double proposal() { Node leaf = tree.getNode(Randomizer.nextInt(leafNodeCount)); Node parent = leaf.getParent(); - if (((ZeroBranchSANode)leaf).isDirectAncestor()) { + if (leaf.isDirectAncestor()) { oldRange = (double) 1; if (parent.isRoot()) { final double randomNumber = Randomizer.nextExponential(1); @@ -84,7 +82,7 @@ public double proposal() { parent.setHeight(newHeight); //make sure that either there are no direct ancestors or r<1 - if ((rInput.get() != null) && (((ZeroBranchSATree)tree).getDirectAncestorNodeCount() > 0 && rInput.get().getValue() == 1)) { + if ((rInput.get() != null) && (tree.getDirectAncestorNodeCount() > 0 && rInput.get().getValue() == 1)) { return Double.NEGATIVE_INFINITY; } diff --git a/src/beast/evolution/operators/UniformForZeroBranchSATrees.java b/src/beast/evolution/operators/UniformForZeroBranchSATrees.java deleted file mode 100644 index 90fe2df..0000000 --- a/src/beast/evolution/operators/UniformForZeroBranchSATrees.java +++ /dev/null @@ -1,57 +0,0 @@ -package beast.evolution.operators; - -import beast.core.Description; -import beast.evolution.tree.Node; -import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.evolution.tree.ZeroBranchSATree; -import beast.util.Randomizer; - -/** - * @author Alexandra Gavryushkina - */ - -/** - * @deprecated Use SAUniform instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -@Description("Randomly selects true internal node (i.e. not the root and not a fake node) and move node height uniformly in interval " + - "restricted by the node's parent and children.") -public class UniformForZeroBranchSATrees extends TreeOperator { - - @Override - public void initAndValidate() { - } - - /** - * change the parameter and return the hastings ratio. - * - * @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted * - */ - @Override - public double proposal() { - final Tree tree = treeInput.get(this); - final int nNodeCount = tree.getNodeCount(); - int leafNodeCount = tree.getLeafNodeCount(); - - //make sure that there is at least one non-fake and non-root internal node - int fakeNodeCount = ((ZeroBranchSATree)tree).getDirectAncestorNodeCount(); - if (fakeNodeCount == leafNodeCount-1 || (fakeNodeCount == leafNodeCount-2 && !((ZeroBranchSANode)tree.getRoot()).isFake())) { - return Double.NEGATIVE_INFINITY; - } - - // randomly select internal node - Node node; - do { - node = tree.getNode(leafNodeCount + Randomizer.nextInt(nNodeCount / 2)); - } while (node.isRoot() || ((ZeroBranchSANode)node).isFake()); - final double fUpper = node.getParent().getHeight(); - final double fLower = Math.max(node.getLeft().getHeight(), node.getRight().getHeight()); - final double newValue = (Randomizer.nextDouble() * (fUpper - fLower)) + fLower; - node.setHeight(newValue); - - return 0.0; - } - -} diff --git a/src/beast/evolution/operators/WilsonBaldingForZeroBranchSampledAncestorTrees.java b/src/beast/evolution/operators/WilsonBaldingForZeroBranchSampledAncestorTrees.java deleted file mode 100644 index c3ba6d5..0000000 --- a/src/beast/evolution/operators/WilsonBaldingForZeroBranchSampledAncestorTrees.java +++ /dev/null @@ -1,193 +0,0 @@ -package beast.evolution.operators; - -import beast.core.Input; -import beast.core.parameter.RealParameter; -import beast.evolution.tree.Node; -import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.evolution.tree.ZeroBranchSATree; -import beast.util.Randomizer; - -/** - *@author Alexandra Gavryushkina - */ - -/** - * @deprecated Use SAWilsonBalding instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -public class WilsonBaldingForZeroBranchSampledAncestorTrees extends TreeOperator { - - public Input rInput = - new Input("removalProbability", "The probability of an individual to become noninfectious immediately after the sampling"); - - @Override - public void initAndValidate() { - } - - /** - * @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted * - */ - @Override - public double proposal() { - - Tree tree = treeInput.get(this); - - //double x0 = 10; - - double oldMinAge, newMinAge, newRange, oldRange, newAge, fHastingsRatio, DimensionCoefficient; - int newDimension, oldDimension; - - // choose a random node avoiding root and leaves that are direct ancestors - int nodeCount = tree.getNodeCount(); - Node i; - - do { - i = tree.getNode(Randomizer.nextInt(nodeCount)); - } while (i.isRoot() || (((ZeroBranchSANode)i).isDirectAncestor())); - - Node iP = i.getParent(); - Node CiP; - if (iP.getLeft().getNr() == i.getNr()) { - CiP = iP.getRight(); - } else { - CiP = iP.getLeft(); - } - - // make sure that there is at least one candidate edge to attach node iP to - if (iP.getParent() == null && CiP.getHeight() < i.getHeight()) { - return Double.NEGATIVE_INFINITY; - } - - // choose another random node to insert i above or to attach i to this node if it is a leaf - Node j; - Node jP; - - final int leafNodeCount = tree.getLeafNodeCount(); - - if (leafNodeCount != tree.getExternalNodes().size()) { - System.out.println("node counts are incorrect. NodeCount = " + nodeCount + " leafNodeCount = " + leafNodeCount + " exteranl node count = " + tree.getExternalNodes().size()); - } - - // make sure that the target branch or target leaf j is above the subtree being moved - - int nodeNumber; - double newParentHeight; - boolean attachingToLeaf; - boolean adjacentEdge; - //boolean adjacentLeaf; - do { - adjacentEdge = false; - //adjacentLeaf = false; - nodeNumber = Randomizer.nextInt(nodeCount + leafNodeCount); - if (nodeNumber < nodeCount) { - j = tree.getNode(nodeNumber); - jP = j.getParent(); - if (jP != null) - newParentHeight = jP.getHeight(); - else newParentHeight = Double.POSITIVE_INFINITY; - if (!((ZeroBranchSANode)CiP).isDirectAncestor()) - adjacentEdge = (CiP.getNr() == j.getNr() || iP.getNr() == j.getNr()); - attachingToLeaf = false; - } else { - j = tree.getExternalNodes().get(nodeNumber - nodeCount); - jP = j.getParent(); - newParentHeight = j.getHeight(); - attachingToLeaf = true; - //adjacentLeaf = (iP.getNr() == j.getNr()); - } - } while (((ZeroBranchSANode)j).isDirectAncestor() || (newParentHeight <= i.getHeight()) || (i.getNr() == j.getNr()) || adjacentEdge /*|| adjacentLeaf */); - - - if (attachingToLeaf && iP.getNr() == j.getNr()) { - System.out.println("Proposal failed because j = iP"); - return Double.NEGATIVE_INFINITY; - } - - if (jP != null && jP.getNr() == i.getNr()) { - System.out.println("Proposal failed because jP = i. Heights of i = " + i.getHeight() + " Height of jP = " + jP.getHeight()); - return Double.NEGATIVE_INFINITY; - } - - oldDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 1; - - //Hastings numerator calculation + newAge of iP - if (attachingToLeaf) { - newRange = 1; - newAge = j.getHeight(); - } else { - if (jP != null) { - newMinAge = Math.max(i.getHeight(), j.getHeight()); - newRange = jP.getHeight() - newMinAge; - newAge = newMinAge + (Randomizer.nextDouble() * newRange); - } else { - double randomNumberFromExponential; - randomNumberFromExponential = Randomizer.nextExponential(1); - //newRange = x0 - j.getHeight(); - //randomNumberFromExponential = Randomizer.nextDouble() * newRange; - newRange = Math.exp(randomNumberFromExponential); - newAge = j.getHeight() + randomNumberFromExponential; - } - } - - Node PiP = iP.getParent(); - - //Hastings denominator calculation - if (((ZeroBranchSANode)CiP).isDirectAncestor()) { - oldRange = 1; - } - else { - oldMinAge = Math.max(i.getHeight(), CiP.getHeight()); - if (PiP != null) { - oldRange = PiP.getHeight() - oldMinAge; - } else { - oldRange = Math.exp(iP.getHeight() - oldMinAge); - //oldRange = x0 - oldMinAge; - } - } - - //update - if (iP.getNr() != j.getNr() && CiP.getNr() != j.getNr()) { - iP.removeChild(CiP); //remove - - if (PiP != null) { - PiP.removeChild(iP); // remove - PiP.addChild(CiP); // add - PiP.makeDirty(Tree.IS_FILTHY); - CiP.makeDirty(Tree.IS_FILTHY); - } else { - CiP.setParent(null); // completely remove - tree.setRootOnly(CiP); - } - - if (jP != null) { - jP.removeChild(j); // remove - jP.addChild(iP); // add - jP.makeDirty(Tree.IS_FILTHY); - } else { - iP.setParent(null); // completely remove - tree.setRootOnly(iP); - } - iP.addChild(j); - iP.makeDirty(Tree.IS_FILTHY); - j.makeDirty(Tree.IS_FILTHY); - } - iP.setHeight(newAge); - - //make sure that either there are no direct ancestors or r<1 - if ((rInput.get() != null) && (((ZeroBranchSATree)tree).getDirectAncestorNodeCount() > 0 && rInput.get().getValue() == 1)) { - return Double.NEGATIVE_INFINITY; - } - - newDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 1; - DimensionCoefficient = (double) oldDimension / newDimension; - - fHastingsRatio = Math.abs(DimensionCoefficient * newRange / oldRange); - - return Math.log(fHastingsRatio); - - } - - -} diff --git a/src/beast/evolution/operators/WilsonBaldingWithRateCategories.java b/src/beast/evolution/operators/WilsonBaldingWithRateCategories.java index 2de3a6a..609603c 100644 --- a/src/beast/evolution/operators/WilsonBaldingWithRateCategories.java +++ b/src/beast/evolution/operators/WilsonBaldingWithRateCategories.java @@ -3,11 +3,8 @@ import beast.core.Description; import beast.core.Input; import beast.core.parameter.IntegerParameter; -import beast.core.parameter.RealParameter; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; -import beast.evolution.tree.ZeroBranchSATree; import beast.util.Randomizer; /** @@ -53,7 +50,7 @@ public double proposal() { do { i = tree.getNode(Randomizer.nextInt(nodeCount)); iP = i.getParent(); - } while (i.isRoot() || ((ZeroBranchSANode)i).isDirectAncestor() || iP.isRoot()); + } while (i.isRoot() || i.isDirectAncestor() || iP.isRoot()); Node CiP; @@ -83,7 +80,7 @@ public double proposal() { if (jP != null) newParentHeight = jP.getHeight(); else newParentHeight = Double.POSITIVE_INFINITY; - if (!((ZeroBranchSANode)CiP).isDirectAncestor()) + if (!CiP.isDirectAncestor()) adjacentEdge = (CiP.getNr() == j.getNr() || iP.getNr() == j.getNr()); attachingToLeaf = false; } else { @@ -93,7 +90,7 @@ public double proposal() { attachingToLeaf = true; //adjacentLeaf = (iP.getNr() == j.getNr()); } - } while (jP == null || ((ZeroBranchSANode)j).isDirectAncestor() || (newParentHeight <= i.getHeight()) || (i.getNr() == j.getNr()) || adjacentEdge /*|| adjacentLeaf */); + } while (jP == null || j.isDirectAncestor() || (newParentHeight <= i.getHeight()) || (i.getNr() == j.getNr()) || adjacentEdge /*|| adjacentLeaf */); if (attachingToLeaf && iP.getNr() == j.getNr()) { @@ -106,10 +103,10 @@ public double proposal() { return Double.NEGATIVE_INFINITY; } - if (((ZeroBranchSANode)tree.getRoot()).isFake()) { - oldDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 2; + if (tree.getRoot().isFake()) { + oldDimension = nodeCount - tree.getDirectAncestorNodeCount() - 2; } else { - oldDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 3; + oldDimension = nodeCount - tree.getDirectAncestorNodeCount() - 3; } int iPCategory = categoriesInput.get().getValue(iP.getNr()), @@ -135,7 +132,7 @@ public double proposal() { //Hastings denominator calculation - if (((ZeroBranchSANode)CiP).isDirectAncestor()) { + if (CiP.isDirectAncestor()) { oldRange = (double) 1; categoriesInput.get().setValue(CiP.getNr(), iPCategory); } @@ -160,10 +157,10 @@ public double proposal() { } iP.setHeight(newAge); - if (((ZeroBranchSANode)tree.getRoot()).isFake()) { - newDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 2; + if (tree.getRoot().isFake()) { + newDimension = nodeCount - tree.getDirectAncestorNodeCount() - 2; } else { - newDimension = nodeCount - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 3; + newDimension = nodeCount - tree.getDirectAncestorNodeCount() - 3; } DimensionCoefficient = (double) oldDimension / newDimension; diff --git a/src/beast/evolution/speciation/ParameterizedSABirthDeathModel.java b/src/beast/evolution/speciation/ParameterizedSABirthDeathModel.java index 91eb8ea..7db1687 100644 --- a/src/beast/evolution/speciation/ParameterizedSABirthDeathModel.java +++ b/src/beast/evolution/speciation/ParameterizedSABirthDeathModel.java @@ -54,7 +54,8 @@ public class ParameterizedSABirthDeathModel extends SpeciesTreeDistribution { private boolean lambdaExceedsMu = false; - public void initAndValidate() { + @Override + public void initAndValidate() { updateParameters(); diff --git a/src/beast/evolution/speciation/RateParameterization.java b/src/beast/evolution/speciation/RateParameterization.java index 215fdcd..6fc33a7 100644 --- a/src/beast/evolution/speciation/RateParameterization.java +++ b/src/beast/evolution/speciation/RateParameterization.java @@ -1,11 +1,8 @@ package beast.evolution.speciation; -import beast.core.BEASTInterface; import beast.core.Input; import beast.core.parameter.RealParameter; -import java.util.Set; - /** * Created by alexei on 7/09/15. */ diff --git a/src/beast/evolution/speciation/SABDSamplingThroughTimeModel.java b/src/beast/evolution/speciation/SABDSamplingThroughTimeModel.java deleted file mode 100644 index 1e4e72d..0000000 --- a/src/beast/evolution/speciation/SABDSamplingThroughTimeModel.java +++ /dev/null @@ -1,252 +0,0 @@ -package beast.evolution.speciation; - -import beast.core.Citation; -import beast.core.Description; -import beast.core.Input; -import beast.core.parameter.RealParameter; -import beast.evolution.tree.*; -import beast.util.ZeroBranchSATreeParser; - -/** - * @author Alexandra Gavryushkina - */ - -//The tree density is from: Tanja Stadler et al. "Estimating the Basic Reproductive Number from Viral Sequence Data" - -/** - * @deprecated Use SABirthDeathModel instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -@Description("Calculate tree density under Birth Death Sampling Through Time Model for Epidemics " + - "that is the BDM where an individual is sampled at a time with a constant rate psi" + - " and where an individual becomes noninfectious immediately after the sampling" + - "with a constant probability r") -@Citation(value = "Gavryushkina A, Welch D, Stadler T, Drummond AJ (2014) \n" + - "Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration. \n" + - "PLoS Comput Biol 10(12): e1003919. doi:10.1371/journal.pcbi.1003919", - year = 2014, firstAuthorSurname = "Gavryushkina", DOI="10.1371/journal.pcbi.1003919") -public class SABDSamplingThroughTimeModel extends SpeciesTreeDistribution { - - public Input originInput = - new Input("origin", "The origin of infection",(RealParameter)null); - - //'direct' parameters - public Input birthRateInput = - new Input("birthRate", "Birth rate"); - public Input deathRateInput = - new Input("deathRate", "Death rate"); - public Input samplingRateInput = - new Input("samplingRate", "Sampling rate per individual"); - - //transformed parameters: - public Input diversificationRateInput = - new Input("diversificationRate", "Net diversification rate. Birth rate - death rate", Input.Validate.XOR, birthRateInput); - public Input turnoverInput = - new Input("turnover", "Turnover. Death rate/birth rate", Input.Validate.XOR, deathRateInput); - public Input samplingProportionInput = - new Input("samplingProportion", "The probability of sampling prior to death. Sampling rate/(sampling rate + death rate)", Input.Validate.XOR, samplingRateInput); - - - // r parameter - public Input removalProbability = - new Input("removalProbability", "The probability that an individual is removed from the process after the sampling", Input.Validate.REQUIRED); - - public Input rhoProbability = - new Input("rho", "Probability of an individual to be sampled at present", (RealParameter)null); - - // if the tree likelihood is condition on sampling at least one individual then set to true one of the inputs: - public Input conditionOnSamplingInput = new Input("conditionOnSampling", "the tree " + - "likelihood is conditioned on sampling at least one individual", false); - public Input conditionOnRhoSamplingInput = new Input("conditionOnRhoSampling", "the tree " + - "likelihood is conditioned on sampling at least one individual in present", false); - - public Input conditionOnRootInput = new Input("conditionOnRoot", "the tree " + - "likelihood is conditioned on the root height otherwise on the time of origin", false); - - - protected double r; - protected double lambda; - protected double mu; - protected double psi; - protected double c1; - protected double c2; - protected double origin; - protected double rho; - protected boolean transform; //is true if the model is parametrised through transformed parameters - - public void initAndValidate() { - - if (originInput.get() == null && !conditionOnRootInput.get()) { - throw new IllegalArgumentException("Either specify origin input or set conditionOnRoot input to \"true\""); - } - - if (originInput.get() != null && conditionOnRootInput.get()){ - throw new IllegalArgumentException("Either don't specify origin input or set conditionOnRoot input to \"false\""); - } - - - if (conditionOnSamplingInput.get() && conditionOnRhoSamplingInput.get()){ - throw new IllegalArgumentException("Either set to \"true\" only one of conditionOnSampling and conditionOnRhoSampling inputs or don't specify both!"); - } - - if (birthRateInput.get() != null && deathRateInput.get() != null && samplingRateInput.get() != null) { - - transform = false; - //mu = deathRateInput.get().getValue(); - //psi = samplingRateInput.get().getValue(); - //lambda = birthRateInput.get().getValue(); - - } else if (diversificationRateInput.get() != null && turnoverInput.get() != null && samplingProportionInput.get() != null) { - - transform = true; - - } else { - throw new IllegalArgumentException("Either specify birthRate, deathRate and samplingRate OR specify diversificationRate, turnover and samplingProportion!"); - } - - if (treeInput.get() instanceof ZeroBranchSATreeParser && originInput.get() != null && originInput.get().getValue() < treeInput.get().getRoot().getHeight()){ - throw new IllegalArgumentException("Initial value of origin should be greater than initial root height"); - - } - - -// r = becomeNoninfectiousAfterSamplingProbability.get().getValue(); -// if (rhoProbability.get() != null ) { -// rho = rhoProbability.get().getValue(); -// } else { -// rho = 0.; -// } -// c1 = Math.sqrt((lambda - mu - psi) * (lambda - mu - psi) + 4 * lambda * psi); -// c2 = -(lambda - mu - 2*lambda*rho - psi) / c1; -// origin = originInput.get().getValue(); - } - - private double p0s(double t, double c1, double c2) { - double p0 = (lambda + mu + psi - c1 * ((1 + c2) - Math.exp(-c1 * t) * (1 - c2)) / ((1 + c2) + Math.exp(-c1 * t) * (1 - c2))) / (2 * lambda); - return r + (1 - r) * p0; - } - - private double oneMinusP0(double t, double c1, double c2) { - return 1 - (lambda + mu + psi - c1 * ((1 + c2) - Math.exp(-c1 * t) * (1 - c2)) / ((1 + c2) + Math.exp(-c1 * t) * (1 - c2))) / (2 * lambda); - } - - private double oneMinusP0Hat(double t, double c1, double c2) { - return rho*(lambda-mu)/(lambda*rho + (lambda*(1-rho) - mu)* Math.exp((mu-lambda) * t)) ; - } - - private double pS(double t) { - return psi*Math.exp((lambda - mu - psi*r) * t); - } - - private double q(double t, double c1, double c2) { - return Math.exp(c1 * t) * (1 + c2) * (1 + c2) + Math.exp(-c1 * t) * (1 - c2) * (1 - c2) + 2 * (1 - c2 * c2); - } - - private void transformParameters() { - double d = diversificationRateInput.get().getValue(); - double r_turnover = turnoverInput.get().getValue(); - double s = samplingProportionInput.get().getValue(); - lambda = d/(1-r_turnover); - mu = r_turnover*lambda; - psi = mu*s/(1-s); - } - - private void updateParameters() { - - if (transform) { - transformParameters(); - } else { - lambda = birthRateInput.get().getValue(); - mu = deathRateInput.get().getValue(); - psi = samplingRateInput.get().getValue(); - } - - r = removalProbability.get().getValue(); - if (rhoProbability.get() != null ) { - rho = rhoProbability.get().getValue(); - } else { - rho = 0.; - } - c1 = Math.sqrt((lambda - mu - psi) * (lambda - mu - psi) + 4 * lambda * psi); - c2 = -(lambda - mu - 2*lambda*rho - psi) / c1; - if (originInput.get() != null){ - origin = originInput.get().getValue(); - } else { - origin = Double.POSITIVE_INFINITY; - } - } - - @Override - public double calculateTreeLogLikelihood(TreeInterface tree) - { - int nodeCount = tree.getNodeCount(); - updateParameters(); - //double x0 = tree.getRoot().getHeight() + origToRootDistance; - double x0 = origin; - double x1=tree.getRoot().getHeight(); - - if (x0 < x1 || r==1 && ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() > 0) { - return Double.NEGATIVE_INFINITY; - } - - double logPost; - if (!conditionOnRootInput.get()){ - logPost = -Math.log(q(x0, c1, c2)); - } else { - if (((ZeroBranchSANode)tree.getRoot()).isFake()){ //when conditioning on the root we assume the process - //starts at the time of the first branching event and - //that means that the root can not be a sampled ancestor - return Double.NEGATIVE_INFINITY; - } else { - logPost = -Math.log(q(x1, c1, c2)); - } - } - - if (conditionOnSamplingInput.get()) { - logPost -= Math.log(oneMinusP0(x0, c1, c2)); - } - - if (conditionOnRhoSamplingInput.get()) { - if (conditionOnRootInput.get()) { - logPost -= Math.log(lambda*oneMinusP0Hat(x1, c1, c2)* oneMinusP0Hat(x1, c1, c2)); - } else { - logPost -= Math.log(oneMinusP0Hat(x0, c1, c2)); - } - } - - int internalNodeCount = tree.getLeafNodeCount() - ((ZeroBranchSATree)tree).getDirectAncestorNodeCount() - 1; - - logPost += internalNodeCount*Math.log(2); - - for (int i = 0; i < nodeCount; i++) { - if (tree.getNode(i).isLeaf()) { - if (!((ZeroBranchSANode)tree.getNode(i)).isDirectAncestor()) { - if (tree.getNode(i).getHeight() > 0.000000000005 || rho == 0.) { - logPost += Math.log(psi) + Math.log(q(tree.getNode(i).getHeight(), c1, c2)) + Math.log(p0s(tree.getNode(i).getHeight(), c1, c2)); - } else { - logPost += Math.log(4*rho); - } - } - } else { - if (((ZeroBranchSANode)tree.getNode(i)).isFake()) { - if (r == 1) { - System.out.println("r = 1 but there are sampled ancestors in the tree"); - System.exit(0); - } - logPost += Math.log(psi) + Math.log(1 - r); - } else { - logPost += Math.log(lambda) - Math.log(q(tree.getNode(i).getHeight(), c1, c2)); - } - } - } - - return logPost; - } - - @Override - protected boolean requiresRecalculation() { - return true; - } - -} diff --git a/src/beast/evolution/speciation/SABirthDeathModel.java b/src/beast/evolution/speciation/SABirthDeathModel.java index fb30490..79846b8 100644 --- a/src/beast/evolution/speciation/SABirthDeathModel.java +++ b/src/beast/evolution/speciation/SABirthDeathModel.java @@ -8,7 +8,9 @@ import beast.evolution.operators.*; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; +import beast.evolution.tree.TreeDistribution; import beast.evolution.tree.TreeInterface; +import beast.evolution.tree.TreeWOffset; import beast.math.distributions.Uniform; import java.util.List; @@ -27,9 +29,11 @@ "Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration. \n" + "PLoS Comput Biol 10(12): e1003919. doi:10.1371/journal.pcbi.1003919", year = 2014, firstAuthorSurname = "Gavryushkina", DOI="10.1371/journal.pcbi.1003919") -public class SABirthDeathModel extends SpeciesTreeDistribution { - +public class SABirthDeathModel extends TreeDistribution { + public Input treeWOffsetInput = + new Input("treeWOffset", "Optional fully extinct tree", (TreeWOffset)null); + //'direct' parameters public Input originInput = new Input("origin", "The time when the process started", (RealParameter)null); @@ -90,6 +94,8 @@ public class SABirthDeathModel extends SpeciesTreeDistribution { protected boolean lambdaExceedsMu = false; protected String taxonName; protected double taxonAge; + + TreeWOffset combinedTree; public void initAndValidate() { @@ -123,8 +129,17 @@ public void initAndValidate() { } else { throw new IllegalArgumentException("Either specify birthRate, deathRate and samplingRate OR specify diversificationRate, turnover and samplingProportion!"); } - - double rootHeight = treeInput.get().getRoot().getHeight(); + + combinedTree = treeWOffsetInput.get(); + if(combinedTree == null) { + combinedTree = new TreeWOffset(); + combinedTree.setInputValue("tree", treeInput.get()); + combinedTree.initAndValidate(); + } + else if(treeInput.get() != null) { + System.err.println("Both tree and treeWOffset specified as inputs, using treeWOffset."); + } + double rootHeight = combinedTree.getTree().getRoot().getHeight(); if (originSpecified && origin() < rootHeight){ throw new IllegalArgumentException("Initial value of origin (" + origin() + ") should be greater than initial root height (" +rootHeight + ")"); } @@ -328,8 +343,9 @@ protected void updateParameters() { } @Override - public double calculateTreeLogLikelihood(TreeInterface tree) - { + public double calculateLogP() { + Tree tree = combinedTree.getTree(); + int nodeCount = tree.getNodeCount(); updateParameters(); if (lambdaExceedsMu && lambda <= mu) { @@ -348,72 +364,72 @@ public double calculateTreeLogLikelihood(TreeInterface tree) if (taxonAge > origin) { return Double.NEGATIVE_INFINITY; } - double logPost = 0.0; + logP = 0.0; if (conditionOnSamplingInput.get()) { - logPost -= Math.log(oneMinusP0(x0, c1, c2)); + logP -= Math.log(oneMinusP0(x0, c1, c2)); } if (conditionOnRhoSamplingInput.get()) { - logPost -= Math.log(oneMinusP0Hat(x0, c1, c2)); + logP -= Math.log(oneMinusP0Hat(x0, c1, c2)); } if (SATaxonInput.get().getValue() == 0) { - logPost += Math.log(1 - oneMinusP0(taxonAge, c1, c2)); + logP += Math.log(1 - oneMinusP0(taxonAge, c1, c2)); } else { - logPost += Math.log(oneMinusP0(taxonAge, c1, c2)); + logP += Math.log(oneMinusP0(taxonAge, c1, c2)); } - return logPost; + return logP; } - double x1=tree.getRoot().getHeight(); + double x1 = combinedTree.getHeightOfNode(tree.getRoot().getNr()); - if (x0 < x1 || r==1 && ((Tree)tree).getDirectAncestorNodeCount() > 0) { + if (x0 < x1 || r==1 && tree.getDirectAncestorNodeCount() > 0) { return Double.NEGATIVE_INFINITY; } - - double logPost; + + logP = 0; if (!conditionOnRootInput.get()){ - logPost = -Math.log(q(x0, c1, c2)); + logP = -Math.log(q(x0, c1, c2)); } else { if (tree.getRoot().isFake()){ //when conditioning on the root we assume the process //starts at the time of the first branching event and //that means that the root can not be a sampled ancestor return Double.NEGATIVE_INFINITY; } else { - logPost = -Math.log(q(x1, c1, c2)); + logP = -Math.log(q(x1, c1, c2)); } } if (conditionOnSamplingInput.get()) { if (conditionOnRootInput.get()) { - logPost -= Math.log(lambda*oneMinusP0(x1, c1, c2)* oneMinusP0(x1, c1, c2)); + logP -= Math.log(lambda*oneMinusP0(x1, c1, c2)* oneMinusP0(x1, c1, c2)); } else { - logPost -= Math.log(oneMinusP0(x0, c1, c2)); + logP -= Math.log(oneMinusP0(x0, c1, c2)); } } if (conditionOnRhoSamplingInput.get()) { if (conditionOnRootInput.get()) { - logPost -= Math.log(lambda*oneMinusP0Hat(x1, c1, c2)* oneMinusP0Hat(x1, c1, c2)); + logP -= Math.log(lambda*oneMinusP0Hat(x1, c1, c2)* oneMinusP0Hat(x1, c1, c2)); } else { - logPost -= Math.log(oneMinusP0Hat(x0, c1, c2)); + logP -= Math.log(oneMinusP0Hat(x0, c1, c2)); } } - int internalNodeCount = tree.getLeafNodeCount() - ((Tree)tree).getDirectAncestorNodeCount() - 1; - - logPost += internalNodeCount*Math.log(2); + int internalNodeCount = tree.getLeafNodeCount() - tree.getDirectAncestorNodeCount() - 1; + logP += internalNodeCount*Math.log(2); + for (int i = 0; i < nodeCount; i++) { if (tree.getNode(i).isLeaf()) { if (!tree.getNode(i).isDirectAncestor()) { - if (tree.getNode(i).getHeight() > 0.000000000005 || rho == 0.) { - logPost += Math.log(psi) + Math.log(q(tree.getNode(i).getHeight(), c1, c2)) + Math.log(p0s(tree.getNode(i).getHeight(), c1, c2)); + if (combinedTree.getHeightOfNode(i) > 0.000000000005 || rho == 0.) { + logP += Math.log(psi) + Math.log(q(combinedTree.getHeightOfNode(i), c1, c2)) + Math.log(p0s(combinedTree.getHeightOfNode(i), c1, c2)); } else { - logPost += Math.log(4*rho); + logP += Math.log(4*rho); } } } else { @@ -422,14 +438,14 @@ public double calculateTreeLogLikelihood(TreeInterface tree) System.out.println("r = 1 but there are sampled ancestors in the tree"); System.exit(0); } - logPost += Math.log(psi) + Math.log(1 - r); + logP += Math.log(psi) + Math.log(1 - r); } else { - logPost += Math.log(lambda) - Math.log(q(tree.getNode(i).getHeight(), c1, c2)); + logP += Math.log(lambda) - Math.log(q(combinedTree.getHeightOfNode(i), c1, c2)); } } } - return logPost; + return logP; } @Override diff --git a/src/beast/evolution/tree/CladeConstraint.java b/src/beast/evolution/tree/CladeConstraint.java index 6ab8a93..491cb24 100644 --- a/src/beast/evolution/tree/CladeConstraint.java +++ b/src/beast/evolution/tree/CladeConstraint.java @@ -32,7 +32,6 @@ public class CladeConstraint extends Distribution { private double storedMrcaHeight; - @SuppressWarnings("unchecked") @Override public void initAndValidate() { tree = treeInput.get(); diff --git a/src/beast/evolution/tree/OffsetLogger.java b/src/beast/evolution/tree/OffsetLogger.java new file mode 100644 index 0000000..73c9082 --- /dev/null +++ b/src/beast/evolution/tree/OffsetLogger.java @@ -0,0 +1,37 @@ +package beast.evolution.tree; + +import java.io.PrintStream; + +import beast.core.BEASTObject; +import beast.core.Input; +import beast.core.Input.Validate; +import beast.core.Loggable; + +public class OffsetLogger extends BEASTObject implements Loggable { + + public Input treeWOffsetInput = + new Input("treeWOffset", "Fully extinct tree", Validate.REQUIRED); + + @Override + public void initAndValidate() {} + + @Override + public void init(PrintStream out) { + final Tree tree = treeWOffsetInput.get().getTree(); + if (getID() == null || getID().matches("\\s*")) { + out.print(tree.getID() + ".offset\t"); + } else { + out.print(getID() + "\t"); + } + } + + @Override + public void log(int sample, PrintStream out) { + final TreeWOffset tree = treeWOffsetInput.get(); + out.print(tree.getOffset() + "\t"); + } + + @Override + public void close(PrintStream out) {} + +} diff --git a/src/beast/evolution/tree/SALogger.java b/src/beast/evolution/tree/SALogger.java deleted file mode 100644 index c83d2c3..0000000 --- a/src/beast/evolution/tree/SALogger.java +++ /dev/null @@ -1,62 +0,0 @@ -package beast.evolution.tree; - -import beast.core.CalculationNode; -import beast.core.Function; -import beast.core.Input; -import beast.core.Loggable; - -import java.io.PrintStream; - -/** - * @author Alexandra Gavryushkina - */ - -/** - * @deprecated Use SampledAncestorLogger instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -public class SALogger extends CalculationNode implements Loggable, Function { - public Input treeInput = new Input("tree", "tree to report height for.", Input.Validate.REQUIRED); - - @Override - public void initAndValidate() { - // nothing to do - } - - @Override - public void init(PrintStream out) { - final Tree tree = treeInput.get(); - if (getID() == null || getID().matches("\\s*")) { - out.print(tree.getID() + ".SAcount\t"); - } else { - out.print(getID() + "\t"); - } - } - - @Override - public void log(long nSample, PrintStream out) { - final Tree tree = treeInput.get(); - out.print(((ZeroBranchSATree)tree).getDirectAncestorNodeCount() + "\t"); - } - - @Override - public void close(PrintStream out) { - // nothing to do - } - - @Override - public int getDimension() { - return 1; - } - - @Override - public double getArrayValue() { - return ((ZeroBranchSATree)treeInput.get()).getDirectAncestorNodeCount(); - } - - @Override - public double getArrayValue(int iDim) { - return ((ZeroBranchSATree)treeInput.get()).getDirectAncestorNodeCount(); - } - -} diff --git a/src/beast/evolution/tree/SATreeLogger.java b/src/beast/evolution/tree/SATreeLogger.java deleted file mode 100644 index 14c6416..0000000 --- a/src/beast/evolution/tree/SATreeLogger.java +++ /dev/null @@ -1,26 +0,0 @@ -package beast.evolution.tree; - -import beast.core.Input; -import beast.core.Logger; - -/** - * Alexandra Gavryushkina - */ -/** - * @deprecated This option is not supported. Starting from v.1.1.1 Tree class supports sampled ancestors and - * some of the classes that work with ZeroBranchSATree are not used anymore others are replaced. - * */ -@Deprecated -public class SATreeLogger extends Logger { - - public Input logWithZeroBranchesInput = new Input("logWithZeroBranches", "If it is true than sampled ancestors " + - "are logged as tips on zero branches otherwise they are logged as single child nodes", false); - - @Override - public void initAndValidate() { - super.initAndValidate(); - ZeroBranchSATree tree = (ZeroBranchSATree)loggersInput.get().get(0); - tree.logWithZeroBranches = logWithZeroBranchesInput.get(); - } - -} diff --git a/src/beast/evolution/tree/SampledAncestorLogger.java b/src/beast/evolution/tree/SampledAncestorLogger.java index 53f7bc0..dda7c20 100644 --- a/src/beast/evolution/tree/SampledAncestorLogger.java +++ b/src/beast/evolution/tree/SampledAncestorLogger.java @@ -1,17 +1,17 @@ package beast.evolution.tree; +import java.io.PrintStream; + import beast.core.CalculationNode; import beast.core.Function; import beast.core.Input; import beast.core.Loggable; -import java.io.PrintStream; - /** * @author Alexandra Gavryushkina */ public class SampledAncestorLogger extends CalculationNode implements Loggable, Function { - public Input treeInput = new Input("tree", "tree to report height for.", Input.Validate.REQUIRED); + public Input treeInput = new Input("tree", "tree to report SA count for.", Input.Validate.REQUIRED); @Override public void initAndValidate() { diff --git a/src/beast/evolution/tree/SamplingDate.java b/src/beast/evolution/tree/SamplingDate.java index 79c7ae6..21efc46 100644 --- a/src/beast/evolution/tree/SamplingDate.java +++ b/src/beast/evolution/tree/SamplingDate.java @@ -3,7 +3,6 @@ import beast.core.BEASTObject; import beast.core.Input; import beast.evolution.alignment.Taxon; -import beast.math.distributions.ParametricDistribution; import beast.util.Randomizer; /** @@ -18,7 +17,8 @@ public class SamplingDate extends BEASTObject { private double upper, lower; - public void initAndValidate() { + @Override + public void initAndValidate() { upper=Double.parseDouble(upperInput.get()); lower=Double.parseDouble(lowerInput.get()); if (upper < 0 || lower < 0 || upper < lower) { diff --git a/src/beast/evolution/tree/TreeWOffset.java b/src/beast/evolution/tree/TreeWOffset.java new file mode 100644 index 0000000..2cfcc69 --- /dev/null +++ b/src/beast/evolution/tree/TreeWOffset.java @@ -0,0 +1,109 @@ +package beast.evolution.tree; + +import beast.core.CalculationNode; +import beast.core.Input; +import beast.core.Input.Validate; + +public class TreeWOffset extends CalculationNode { + + public Input treeInput = new Input ("tree", "tree topology", Validate.REQUIRED); + public Input offsetInput = new Input("offset", "Starting offset to present, default 0", 0.0); + + Tree tree; + Double oldOffset, offset; + Double[] oldHeights, leaves_heights; + int oldMin_leaf, min_leaf; + private boolean stored = false; + + @Override + public void initAndValidate() { + tree = treeInput.get(); + offset = offsetInput.get(); + + leaves_heights = new Double[tree.getLeafNodeCount()]; + for(int i = 0; i < tree.getLeafNodeCount(); i++) { + double h = tree.getNode(i).getHeight(); + if(h == 0) min_leaf = i; + leaves_heights[i] = h + offset; + } + } + + public double getHeightOfNode(int nr) { + return tree.getNode(nr).getHeight() + offset; + } + + public void setHeightOfLeaf(int nr, double height) { + store(); + + double oldOffset = offset; + if(height < offset) { + offset = height; + min_leaf = nr; + } + leaves_heights[nr] = height; + // if we increased the height of the minimum leaf then which one it is may have changed + if(height > offset && nr == min_leaf) { + offset = height; + for(int i = 0; i < leaves_heights.length; i++) { + if(leaves_heights[i] < offset) { + min_leaf = i; + offset = leaves_heights[i]; + } + } + } + // if offset has changed then all heights need to be updated to keep true heights the same + if(Math.abs(oldOffset - offset) > 0) { + for(Node n : tree.getNodesAsArray()) n.setHeight(n.getHeight() + oldOffset - offset); + } + } + + public double getOffset() { + return offset; + } + + public Tree getTree() { + return tree; + } + + public void setHeightOfNode(int nr, double height) { + Node n = tree.getNode(nr); + if(n.isLeaf()) this.setHeightOfLeaf(nr, height); + n.setHeight(height - offset); + } + + // for testing purposes + public double getStoredHeightOfLeaf(int nr) { + return leaves_heights[nr]; + } + + /* This entire section is to work around the fact that this class is modified DURING the operator proposal of SampledNodeDateRandomWalker + * whereas Beast2 assumes that only state nodes are modified, and that calculation nodes are only updated afterwards + * So, using the default mechanic will store/restore the modified state not the original one + * There is probably a better way to do this + */ + + @Override + public void store() { + super.store(); + if(stored) return; + oldMin_leaf = min_leaf; + oldOffset = offset; + oldHeights = leaves_heights.clone(); + stored = true; + } + + @Override + public void restore() { + super.restore(); + min_leaf = oldMin_leaf; + offset = oldOffset; + leaves_heights = oldHeights.clone(); + stored = false; + } + + @Override + public void accept() { + super.accept(); + stored = false; + } +} diff --git a/src/beast/evolution/tree/ZeroBranchSANode.java b/src/beast/evolution/tree/ZeroBranchSANode.java deleted file mode 100644 index d14df53..0000000 --- a/src/beast/evolution/tree/ZeroBranchSANode.java +++ /dev/null @@ -1,439 +0,0 @@ -package beast.evolution.tree; - -import beast.util.Randomizer; - -import java.util.List; - -/** - * @author Alexandra Gavryushkina - */ - -/** - * @deprecated Use Node instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -public class ZeroBranchSANode extends Node { - - /** - * @return (deep) copy of node - */ - public Node copy() { - final Node node = new ZeroBranchSANode(); - node.height = height; - node.labelNr = labelNr; - node.metaDataString = metaDataString; - node.parent = null; - node.setID(ID); - - for (final Node child : getChildren()) { - node.addChild(child.copy()); - } - return node; - } // copy - - /** - * Writes a short newick format. - * Note: if bPrintInternalNodeNumbers=false this method suppresses internal - * nodes that have a zero branch length child (replacing them with the child in the - * newick string). - * - * @return beast.tree in Newick format, with length and meta data - * information. Unlike toNewick(), here Nodes are numbered, instead of - * using the node labels. - * All internal nodes are labelled if bPrintInternalNodeNumbers - * is set true. This is useful for example when storing a State to file - * so that it can be restored. - */ - @Override - public String toShortNewick(boolean bPrintInternalNodeNumbers) { - StringBuilder buf = new StringBuilder(); - if (getLeft() != null) { - buf.append("("); - if (isFake() && !bPrintInternalNodeNumbers) { - Node directAncestor; - Node otherChild; - if (((ZeroBranchSANode)getLeft()).isDirectAncestor()) { - directAncestor = getLeft(); - otherChild = getRight(); - } else { - directAncestor = getRight(); - otherChild = getLeft(); - } - buf.append(otherChild.toShortNewick(bPrintInternalNodeNumbers)); - buf.append(")"); - buf.append(directAncestor.getNr()); - } else { - buf.append(getLeft().toShortNewick(bPrintInternalNodeNumbers)); - if (getRight() != null) { - buf.append(','); - buf.append(getRight().toShortNewick(bPrintInternalNodeNumbers)); - } - buf.append(")"); - if (bPrintInternalNodeNumbers) { - buf.append(getNr()); - } - } - } else { - buf.append(getNr()); - } - buf.append(getNewickMetaData()); - buf.append(":").append(getLength()); - return buf.toString(); - } - - /** - * Note: this method suppresses internal nodes that have a zero branch length child - * (replacing them with the child in the newick string). - * - * @param iMaxNodeInClade - * @param printMetaData - * @return - */ - public String toSortedNewick(int[] iMaxNodeInClade, boolean printMetaData) { - StringBuilder buf = new StringBuilder(); - if (getLeft() != null) { - buf.append("("); - if (isFake()) { - Node directAncestor; - Node otherChild; - if (((ZeroBranchSANode)getLeft()).isDirectAncestor()) { - directAncestor = getLeft(); - otherChild = getRight(); - } else { - directAncestor = getRight(); - otherChild = getLeft(); - } - buf.append(otherChild.toSortedNewick(iMaxNodeInClade, printMetaData)); - buf.append(")"); - buf.append(directAncestor.getNr() + 1); - } else { - String sChild1 = getLeft().toSortedNewick(iMaxNodeInClade, printMetaData); - int iChild1 = iMaxNodeInClade[0]; - if (getRight() != null) { - String sChild2 = getRight().toSortedNewick(iMaxNodeInClade, printMetaData); - int iChild2 = iMaxNodeInClade[0]; - if (iChild1 > iChild2) { - buf.append(sChild2); - buf.append(","); - buf.append(sChild1); - } else { - buf.append(sChild1); - buf.append(","); - buf.append(sChild2); - iMaxNodeInClade[0] = iChild1; - } - } else { - buf.append(sChild1); - } - buf.append(")"); - } - } else { - iMaxNodeInClade[0] = labelNr; - buf.append(labelNr + 1); - } - - if (printMetaData) { - buf.append(getNewickMetaData()); - } - buf.append(":").append(getLength()); - return buf.toString(); - } - - public String toSortedNewickWithZeroBranches(int[] iMaxNodeInClade) { - StringBuilder buf = new StringBuilder(); - if (getLeft() != null) { - buf.append("("); - String sChild1 = ((ZeroBranchSANode)getLeft()).toSortedNewickWithZeroBranches(iMaxNodeInClade); - int iChild1 = iMaxNodeInClade[0]; - if (getRight() != null) { - String sChild2 = ((ZeroBranchSANode)getRight()).toSortedNewickWithZeroBranches(iMaxNodeInClade); - int iChild2 = iMaxNodeInClade[0]; - if (iChild1 > iChild2) { - buf.append(sChild2); - buf.append(","); - buf.append(sChild1); - } else { - buf.append(sChild1); - buf.append(","); - buf.append(sChild2); - iMaxNodeInClade[0] = iChild1; - } - } else { - buf.append(sChild1); - } - buf.append(")"); - if (getID() != null) { - buf.append(labelNr+1); - } - } else { - iMaxNodeInClade[0] = labelNr; - buf.append(labelNr + 1); - } - - buf.append(":").append(getLength()); - return buf.toString(); - } - - - /** //TODO make this method work for ZeroBranchSATree - * @return beast.tree in Newick format with taxon labels for labelled tip nodes - * and labeled (having non-null ID) internal nodes. - * If a tip node doesn't have an ID (taxon label) then node number (m_iLabel) is printed. - */ - public String toNewick() { - StringBuilder buf = new StringBuilder(); - if (getLeft() != null) { - buf.append("("); - buf.append(getLeft().toNewick()); - if (getRight() != null) { - buf.append(','); - buf.append(getRight().toNewick()); - } - buf.append(")"); - } else { - if (getID() == null) { - buf.append(labelNr); - } else { - buf.append(getID()); - } - } - buf.append(getNewickMetaData()); - buf.append(":").append(getLength()); - return buf.toString(); - } - - /** - * @param sLabels //TODO make this method work for ZeroBranchSATree - * @return beast.tree in long Newick format, with all length and meta data - * information, but with leafs labelled with their names - */ - public String toString(List sLabels) { - StringBuilder buf = new StringBuilder(); - if (getLeft() != null) { - buf.append("("); - buf.append(getLeft().toString(sLabels)); - if (getRight() != null) { - buf.append(','); - buf.append(getRight().toString(sLabels)); - } - buf.append(")"); - } else { - buf.append(sLabels.get(labelNr)); - } - if (metaDataString != null) { - buf.append('['); - buf.append(metaDataString); - buf.append(']'); - } - buf.append(":").append(getLength()); - return buf.toString(); - } - - public String toString() { - return toShortNewick(true); - } - - /** - * Scales this node and all its descendants (either all descendants, or only non-sampled descendants) - * - * @param fScale the scalar to multiply each scaled node age by - * @param scaleSNodes true if sampled nodes should be scaled as well as internal nodes, false if only non-sampled - * internal nodes should be scaled. - * @throws Exception throws exception if resulting tree would have negative branch lengths. - */ - public void scale(double fScale, boolean scaleSNodes) { - startEditing(); - isDirty |= Tree.IS_DIRTY; - if (scaleSNodes || (!isLeaf() && !isFake())) { - height *= fScale; - } - if (!isLeaf()) { - ((ZeroBranchSANode)getLeft()).scale(fScale, scaleSNodes); - if (getRight() != null) { - ((ZeroBranchSANode)getRight()).scale(fScale, scaleSNodes); - } - if (height < getLeft().height || height < getRight().height) { - throw new IllegalArgumentException("Scale gives negative branch length"); - } - } - } - - public void scaleAndSlide(double fScale, int[] splitCount) throws Exception { - startEditing(); - if (!isLeaf() && !isFake()) { - height *= fScale; - isDirty |= Tree.IS_DIRTY; - } - if (!isLeaf()) { - ((ZeroBranchSANode)getLeft()).scaleAndSlide(fScale, splitCount); - if (getRight() != null) { - ((ZeroBranchSANode)getRight()).scaleAndSlide(fScale, splitCount); - } - if (height < getLeft().height || height < getRight().height) { - if (isFake()){ // this can only be the case if we scale up - - if (fScale < 1.0) { - throw new Exception("Scale and slide doesn't work properly"); - } - - Node child = this.getNonDirectAncestorChild(); //child is the node that has been scaled up causing - //this parent sampled ancestor node to be lower than its child - - //this parent sampled ancestor node that has to be lowered down - Node sampledAncestor = this; - double saHeight = sampledAncestor.getHeight(); - - //apart from the direct parent of child there can be other sampled ancestor nodes lying on - //the branch going from child up to the first non-fake node - //search for such nodes and lower them down - do { - Node newChild = child; - - //search for the right place to put the sampled ancestor node - do { - if (((ZeroBranchSANode)newChild).isFake()) { - newChild = ((ZeroBranchSANode) newChild).getNonDirectAncestorChild(); - } else { - newChild = (Randomizer.nextInt(2) == 1)?newChild.getLeft():newChild.getRight(); - splitCount[0]++; //count how many levels the sa node goes down - } - } while (saHeight < newChild.getHeight()); - - //first remember the sampled ancestor node parent - Node newSampledAncestor = sampledAncestor.getParent(); - - //place the node above its newChild - placeSampledAncestorAboveItsNewChild(sampledAncestor, newChild); - - //update sampled ancestor node - sampledAncestor = newSampledAncestor; - if (sampledAncestor != null) { - saHeight = sampledAncestor.getHeight(); - } - - //continue until you find all closely ancestral sampled ancestor nodes (those are nodes - //that lies on the branch going from this node up to the fist non-fake node). - - } while (sampledAncestor != null && ((ZeroBranchSANode)sampledAncestor).isFake() && saHeight < child.getHeight()); - - } else { // this can only be the case if we scale down - - if (fScale > 1.0) { - throw new Exception("Scale and slide doesn't work properly"); - } - - Node sampledAncestor = getFakeChild(); - double saHeight = sampledAncestor.getHeight(); - - do { - Node newParent = this; - Node newChild; //newChild remembers on which branch to place the sampled ancestor node - //it will be a child of the sampled ancestor node - do { - newChild = newParent; - newParent = newParent.getParent(); - if (newParent!= null && !((ZeroBranchSANode)newParent).isFake()) { - splitCount[0]++; - } - } while (newParent != null && saHeight > newParent.getHeight()); - - - //first remember the sampled ancestor node child - Node newSampledAncestor = ((ZeroBranchSANode)sampledAncestor).getNonDirectAncestorChild(); - - - //place the node above its newChild - placeSampledAncestorAboveItsNewChild(sampledAncestor, newChild); - - //update sampled ancestor node - sampledAncestor = newSampledAncestor; - saHeight = sampledAncestor.getHeight(); - - - //continue until you find all closely descendant sampled ancestor nodes (those are nodes - //that lies on the branch going from this node down to the first non-fake node). - } while (((ZeroBranchSANode)sampledAncestor).isFake() && saHeight > this.getHeight()); - } - } - } - } - - - /** - * @return true if this leaf actually represent a direct ancestor - * (i.e. is on the end of a zero-length branch) - */ - public boolean isDirectAncestor() { - return (isLeaf() && !isRoot() && this.getParent().getHeight() == this.getHeight()); - } - - /** - * @return true if this is a "fake" internal node (i.e. one of its children is a direct ancestor) - */ - public boolean isFake() { - if (this.isLeaf()) - return false; - return (((ZeroBranchSANode)this.getLeft()).isDirectAncestor() || (this.getRight() != null && ((ZeroBranchSANode)this.getRight()).isDirectAncestor())); - } - - public Node getNonDirectAncestorChild(){ - if (((ZeroBranchSANode)this.getLeft()).isDirectAncestor()){ - return getRight(); - } - if (((ZeroBranchSANode)this.getRight()).isDirectAncestor()){ - return getLeft(); - } - return null; - } - - public Node getFakeChild(){ - - if (((ZeroBranchSANode)this.getLeft()).isFake()){ - return getLeft(); - } - if (((ZeroBranchSANode)this.getRight()).isFake()){ - return getRight(); - } - return null; - } - - public void placeSampledAncestorAboveItsNewChild(Node node, Node newChild){ - - //reattach node from the branch it belongs to - Node nodeChild = ((ZeroBranchSANode)node).getNonDirectAncestorChild(); - Node nodeParent = node.getParent(); - if (nodeParent != null) { - int nodePosition = nodeParent.getChildren().indexOf(node); - nodeParent.setChild(nodePosition, nodeChild); - nodeChild.setParent(nodeParent); - } else { - nodeChild.setParent(null); - } - node.removeChild(nodeChild); //you also need to set nodeChild as a new root but do this at the end - - //attach node to the branch above newChild - Node newParent = newChild.getParent(); - if (newParent != null){ - int newChildPosition = newParent.getChildren().indexOf(newChild); - newParent.setChild(newChildPosition, node); - node.setParent(newParent); - - } else { - node.setParent(null); - } - node.addChild(newChild); - if (newParent == null){ - this.m_tree.setRoot(node); - } - - //set the new root here because setting root automatically assigns the number of nodes below a new root - // to tree.nodeCount field and if you were to set it early the nodeCount would not be correct - if (nodeParent == null){ - this.m_tree.setRoot(nodeChild); - } - - - - } -} diff --git a/src/beast/evolution/tree/ZeroBranchSARandomTree.java b/src/beast/evolution/tree/ZeroBranchSARandomTree.java deleted file mode 100644 index acac942..0000000 --- a/src/beast/evolution/tree/ZeroBranchSARandomTree.java +++ /dev/null @@ -1,605 +0,0 @@ -/* - * CoalescentSimulator.java - * - * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut - * - * 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 beast.evolution.tree; - - -import beast.core.Description; -import beast.core.Input; -import beast.core.Input.Validate; -import beast.core.StateNode; -import beast.core.StateNodeInitialiser; -import beast.evolution.alignment.Alignment; -import beast.evolution.alignment.TaxonSet; -import beast.evolution.tree.coalescent.PopulationFunction; -import beast.util.HeapSort; -import beast.util.Randomizer; - -import java.util.ArrayList; -import java.util.BitSet; -import java.util.List; - - - -/** - * @deprecated Use RandomTree instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -@Description("This class provides the basic engine for coalescent simulation of a given demographic model over a given time period. ") -public class ZeroBranchSARandomTree extends ZeroBranchSATree implements StateNodeInitialiser { - public Input taxaInput = new Input("taxa", "set of taxa to initialise tree specified by alignment"); -// public Input m_taxonset = new Input("taxonset","set of taxa to initialise tree with specified by a taxonset", Validate.REQUIRED); - public Input populationFunctionInput = new Input("populationModel", "population function for generating coalescent???", Validate.REQUIRED); - public Input> cladeConstraintsInput = new Input>("cladeConstraint", "specifies monophyletic clades", new ArrayList()); - public Input rootHeightInput = new Input("rootHeight", "If specified the tree will be scaled to match the root height, if constraints allow this"); - - // total nr of taxa - int nrOfTaxa; - // list of bitset representation of the taxon sets - List taxonSets; - // strongMonophyletic inputs for clade constraints - List strongMonophyletic; - // the first m_nIsMonophyletic of the m_bTaxonSets are monophyletic, while the remainder are not - - int isMonophyletic; - - List taxonSetIDs; - - List[] children; - - // number of the next internal node, us ed when creating new internal nodes - int nextNodeNr; - - // used to indicate one of the clade constraints could not be met - protected class ConstraintViolatedException extends Exception { - private static final long serialVersionUID = 1L; - } - - @Override - public void initAndValidate() { - final List sTaxa; - if (taxaInput.get() != null) { - sTaxa = taxaInput.get().getTaxaNames(); - } else { - sTaxa = m_taxonset.get().asStringList(); - } - nrOfTaxa = sTaxa.size(); - - initStateNodes(); - super.initAndValidate(); - } - - @SuppressWarnings({"rawtypes", "unchecked"}) - private void swap(final List list, final int i, final int j) { - final Object tmp = list.get(i); - list.set(i, list.get(j)); - list.set(j, tmp); - } - - // taxonset intersection test - private boolean intersects(final BitSet bitSet, final BitSet bitSet2) { - for (int k = bitSet.nextSetBit(0); k >= 0; k = bitSet.nextSetBit(k + 1)) { - if (bitSet2.get(k)) { - return true; - } - } - return false; - } - - // taxonset subset test - private boolean isSubset(final BitSet bitSet, final BitSet bitSet2) { - boolean bIsSubset = true; - for (int k = bitSet.nextSetBit(0); bIsSubset && k >= 0; k = bitSet.nextSetBit(k + 1)) { - bIsSubset = bitSet2.get(k); - } - return bIsSubset; - } - - //@Override - public void initStateNodes() { - // find taxon sets we are dealing with - taxonSets = new ArrayList(); - strongMonophyletic = new ArrayList(); - taxonSetIDs = new ArrayList(); - isMonophyletic = 0; - - final List sTaxa; - if (taxaInput.get() != null) { - sTaxa = taxaInput.get().getTaxaNames(); - } else { - sTaxa = m_taxonset.get().asStringList(); - } - - - List calibrations = new ArrayList(); - calibrations.addAll(cladeConstraintsInput.get()); - - - for (final CladeConstraint constraint: calibrations) { - final TaxonSet taxonSet = constraint.taxonsetInInput.get(); - if (taxonSet != null) { - final BitSet bTaxa = new BitSet(nrOfTaxa); - if (taxonSet.asStringList() == null) { - taxonSet.initAndValidate(); - } - for (final String sTaxonID : taxonSet.asStringList()) { - final int iID = sTaxa.indexOf(sTaxonID); - if (iID < 0) { - throw new IllegalArgumentException("Taxon <" + sTaxonID + "> could not be found in list of taxa. Choose one of " + sTaxa.toArray(new String[0])); - } - bTaxa.set(iID); - } - - // add monophyletic constraint - taxonSets.add(isMonophyletic, bTaxa); - strongMonophyletic.add(isMonophyletic, constraint.isStronglyMonophyleticInput.get()); - taxonSetIDs.add(constraint.getID()); - isMonophyletic++; - } - } - - isMonophyletic = taxonSets.size(); - - - // sort constraints such that if taxon set i is subset of taxon set j, then i < j - for (int i = 0; i < isMonophyletic; i++) { - for (int j = i + 1; j < isMonophyletic; j++) { - final boolean bIntersects = intersects(taxonSets.get(i), taxonSets.get(j)); - if (bIntersects) { - final boolean bIsSubset = isSubset(taxonSets.get(j), taxonSets.get(i)); - final boolean bIsSubset2 = isSubset(taxonSets.get(i), taxonSets.get(j)); - // sanity check: make sure either - // o taxonset1 is subset of taxonset2 OR - // o taxonset1 is superset of taxonset2 OR - // o taxonset1 does not intersect taxonset2 - if (!(bIsSubset || bIsSubset2)) { - throw new IllegalArgumentException("333: Don't know how to generate a Random Tree for taxon sets that intersect, " + - "but are not inclusive. Taxonset " + taxonSetIDs.get(i) + " and " + taxonSetIDs.get(j)); - } - // swap i & j if b1 subset of b2 - if (bIsSubset) { - swap(taxonSets, i, j); - swap(strongMonophyletic, i, j); - swap(taxonSetIDs, i, j); - } - } - } - } - strongMonophyletic.add(false); - - // build tree of mono constraints such that i is parent of j => j is subset of i - children = new List[isMonophyletic + 1]; - for (int i = 0; i < isMonophyletic + 1; i++) { - children[i] = new ArrayList(); - } - for (int i = 0; i < isMonophyletic; i++) { - int j = i + 1; - while (j < isMonophyletic && !isSubset(taxonSets.get(i), taxonSets.get(j))) { - j++; - } - children[j].add(i); - } - - - final PopulationFunction popFunction = populationFunctionInput.get(); - - simulateTree(sTaxa, popFunction); - if (rootHeightInput.get() != null) { - scaleToFit(rootHeightInput.get() / root.getHeight(), root); - } - - nodeCount = 2 * sTaxa.size() - 1; - internalNodeCount = sTaxa.size() - 1; - leafNodeCount = sTaxa.size(); - initArrays(); - - if (m_initial.get() != null) { - m_initial.get().assignFromWithoutID(this); - } - } - - private void scaleToFit(double scale, Node node) { - if (!node.isLeaf()) { - double oldHeight = node.getHeight(); - node.height *= scale; - scaleToFit(scale, node.getLeft()); - scaleToFit(scale, node.getRight()); - if (node.height < Math.max(node.getLeft().getHeight(), node.getRight().getHeight())) { - // this can happen if a child node is constrained and the default tree is higher than desired - node.height = 1.0000001 * Math.max(node.getLeft().getHeight(), node.getRight().getHeight()); - } - } - } - - //@Override - public void getInitialisedStateNodes(final List stateNodes) { - stateNodes.add(m_initial.get()); - } - - /** - * Simulates a coalescent tree, given a taxon list. - * - * @param taxa the set of taxa to simulate a coalescent tree between - * @param demoFunction the demographic function to use - */ - public void simulateTree(final List taxa, final PopulationFunction demoFunction) { - if (taxa.size() == 0) - return; - - for (int attempts = 0; attempts < 1000; ++attempts) { - try { - nextNodeNr = nrOfTaxa; - final List candidates = new ArrayList(); - for (int i = 0; i < taxa.size(); i++) { - final Node node = new ZeroBranchSANode(); - node.setNr(i); - node.setID(taxa.get(i)); - node.setHeight(0.0); - candidates.add(node); - } - - if (m_initial.get() != null) - processCandidateTraits(candidates, m_initial.get().m_traitList.get()); - else - processCandidateTraits(candidates, m_traitList.get()); - - // TODO: deal with dated taxa - double fMostRecent = 0; - for (final Node node : candidates) { - fMostRecent = Math.max(fMostRecent, node.getHeight()); - } - - final List allCandidates = new ArrayList(); - allCandidates.addAll(candidates); - root = simulateCoalescent(isMonophyletic, allCandidates, candidates, demoFunction); - return; - } catch (ConstraintViolatedException e) { - // need to generate another tree - } - } - } - - /** - * Apply traits to a set of nodes. - * @param candidates List of nodes - * @param traitSets List of TraitSets to apply - */ - private void processCandidateTraits(List candidates, List traitSets) { - for (TraitSet traitSet : traitSets) { - for (Node node : candidates) - node.setMetaData(traitSet.getTraitName(), traitSet.getValue(node.getNr())); - } - } - - - private Node simulateCoalescent(final int iIsMonophyleticNode, final List allCandidates, final List candidates, final PopulationFunction demoFunction) - throws ConstraintViolatedException { - final List remainingCandidates = new ArrayList(); - final BitSet taxaDone = new BitSet(nrOfTaxa); - for (final int iMonoNode : children[iIsMonophyleticNode]) { - // create list of leaf nodes for this monophyletic MRCA - final List candidates2 = new ArrayList(); - final BitSet bTaxonSet = taxonSets.get(iMonoNode); - for (int k = bTaxonSet.nextSetBit(0); k >= 0; k = bTaxonSet.nextSetBit(k + 1)) { - candidates2.add(allCandidates.get(k)); - } - - final Node MRCA = simulateCoalescent(iMonoNode, allCandidates, candidates2, demoFunction); - remainingCandidates.add(MRCA); - taxaDone.or(bTaxonSet); - } - - for (final Node node : candidates) { - if (!taxaDone.get(node.getNr())) { - remainingCandidates.add(node); - } - } - - final Node MRCA = simulateCoalescent(strongMonophyletic.get(iIsMonophyleticNode), remainingCandidates, demoFunction); - return MRCA; - } - - /** - * @param nodes - * @param demographic - * @return the root node of the given array of nodes after simulation of the - * coalescent under the given demographic model. - * @throws beast.evolution.tree.RandomTree.ConstraintViolatedException - */ - public Node simulateCoalescent(boolean isStronglyMonophyletic, final List nodes, final PopulationFunction demographic) throws ConstraintViolatedException { - // sanity check - disjoint trees - - // if( ! Tree.Utils.allDisjoint(nodes) ) { - // throw new RuntimeException("non disjoint trees"); - // } - - if (nodes.size() == 0) { - throw new IllegalArgumentException("empty nodes set"); - } - - for (int attempts = 0; attempts < 1000; ++attempts) { - final List rootNode = simulateCoalescent(isStronglyMonophyletic, nodes, demographic, 0.0, Double.POSITIVE_INFINITY); - if (rootNode.size() == 1) { - return rootNode.get(0); - } - } - - throw new RuntimeException("failed to merge trees after 1000 tries!"); - } - - public List simulateCoalescent(boolean isStronglyMonophyletic, final List nodes, final PopulationFunction demographic, double currentHeight, - final double maxHeight) throws ConstraintViolatedException { - - Node extantNode = nodes.get(0); - - if (isStronglyMonophyletic) { - for (int i = 0; i < nodes.size(); i++) { - if (nodes.get(i).getHeight() == 0.0) { - extantNode = nodes.get(i); - nodes.remove(i); - break; - } - } -// int extantTaxaCount =0; -// for (int i = 0; i < nodes.size(); i++) { -// if (nodes.get(i).getHeight() == 0.0) { -// extantTaxaCount++; -// if (extantTaxaCount == 2) { -// extantNode = nodes.get(i); -// nodes.remove(i); -// break; -// } -// } -// } -// if (extantTaxaCount < 2) { -// throw new Exception("there is only one extant taxon in a clade that is specified as strongly monophyletic"); -// } - } - - // If only one node, return it - // continuing results in an infinite loop - if (nodes.size() == 1 && !isStronglyMonophyletic) { - return nodes; - } - - - - final double[] heights = new double[nodes.size()]; - for (int i = 0; i < nodes.size(); i++) { - heights[i] = nodes.get(i).getHeight(); - } - final int[] indices = new int[nodes.size()]; - HeapSort.sort(heights, indices); - - // node list - nodeList.clear(); - activeNodeCount = 0; - for (int i = 0; i < nodes.size(); i++) { - nodeList.add(nodes.get(indices[i])); - } - - if (nodes.size() != 1) { - setCurrentHeight(currentHeight); - - // get at least two tips - while (getActiveNodeCount() < 2) { - currentHeight = getMinimumInactiveHeight(); - setCurrentHeight(currentHeight); - } - - // simulate coalescent events - double nextCoalescentHeight = currentHeight - + PopulationFunction.Utils.getSimulatedInterval(demographic, getActiveNodeCount(), currentHeight); - - // while (nextCoalescentHeight < maxHeight && (getNodeCount() > 1)) { - while (nextCoalescentHeight < maxHeight && (nodeList.size() > 1)) { - - if (nextCoalescentHeight >= getMinimumInactiveHeight()) { - currentHeight = getMinimumInactiveHeight(); - setCurrentHeight(currentHeight); - } else { - currentHeight = coalesceTwoActiveNodes(currentHeight, nextCoalescentHeight); - } - - // if (getNodeCount() > 1) { - if (nodeList.size() > 1) { - // get at least two tips - while (getActiveNodeCount() < 2) { - currentHeight = getMinimumInactiveHeight(); - setCurrentHeight(currentHeight); - } - - // nextCoalescentHeight = currentHeight + - // DemographicFunction.Utils.getMedianInterval(demographic, - // getActiveNodeCount(), currentHeight); - nextCoalescentHeight = currentHeight - + PopulationFunction.Utils.getSimulatedInterval(demographic, getActiveNodeCount(), - currentHeight); - } - } - } - - if (isStronglyMonophyletic && nodeList.size() == 1) { - Node currentMRCA = nodeList.get(0); - double nextCoalescentHeight = currentMRCA.getHeight() + PopulationFunction.Utils.getSimulatedInterval(demographic, 2, - currentMRCA.getHeight()); - final Node newNode = new ZeroBranchSANode(); - newNode.setNr(nextNodeNr++); - newNode.setHeight(nextCoalescentHeight); - newNode.setLeft(extantNode); - extantNode.setParent(newNode); - newNode.setRight(currentMRCA); - currentMRCA.setParent(newNode); - nodeList.remove(0); - nodeList.add(newNode); - } - - return nodeList; - } - - /** - * @return the height of youngest inactive node. - */ - private double getMinimumInactiveHeight() { - if (activeNodeCount < nodeList.size()) { - return (nodeList.get(activeNodeCount)).getHeight(); - } else - return Double.POSITIVE_INFINITY; - } - - /** - * Set the current height. - * @param height - */ - private void setCurrentHeight(final double height) { - while (getMinimumInactiveHeight() <= height) { - activeNodeCount += 1; - } - } - - /** - * @return the number of active nodes (equate to lineages) - */ - private int getActiveNodeCount() { - return activeNodeCount; - } - - // - // /** - // * @return the total number of nodes both active and inactive - // */ - // private int getNodeCount() { - // return nodeList.size(); - // } - // - - /** - * Coalesce two nodes in the active list. This method removes the two - * (randomly selected) active nodes and replaces them with the new node at - * the top of the active list. - * @param fMinHeight - * @param height - * @return - */ - private double coalesceTwoActiveNodes(final double fMinHeight, double height) throws ConstraintViolatedException { - final int node1 = Randomizer.nextInt(activeNodeCount); - int node2 = node1; - while (node2 == node1) { - node2 = Randomizer.nextInt(activeNodeCount); - } - - final Node left = nodeList.get(node1); - final Node right = nodeList.get(node2); - - final Node newNode = new ZeroBranchSANode(); -// System.err.println(2 * m_taxa.get().getNrTaxa() - nodeList.size()); - newNode.setNr(nextNodeNr++); - newNode.setHeight(height); - newNode.setLeft(left); - left.setParent(newNode); - newNode.setRight(right); - right.setParent(newNode); - - nodeList.remove(left); - nodeList.remove(right); - - activeNodeCount -= 2; - - nodeList.add(activeNodeCount, newNode); - - activeNodeCount += 1; - - if (getMinimumInactiveHeight() < height) { - throw new RuntimeException( - "This should never happen! Somehow the current active node is older than the next inactive node!"); - } - return height; - } - - private Node findNodeWithExtantDescendantsOnBothSides(Node node, int[] nodeIndex) { - if (node.isLeaf()) { - nodeIndex[0] = -1; - } else { - if (CladeConstraint.hasExtantDescendant(node)) { - nodeIndex[0] = node.getNr(); - return node; - } else { - Node nodeOnTheLeft=findNodeWithExtantDescendantsOnBothSides(node.getLeft(), nodeIndex); - if (nodeIndex[0] > 0) return nodeOnTheLeft; - Node nodeOnTheRight = findNodeWithExtantDescendantsOnBothSides(node.getRight(), nodeIndex); - if (nodeIndex[0] > 0) return nodeOnTheRight; - } - } - return node; - } - - int traverse(final Node node, final BitSet MRCATaxonSet, final int nNrOfMRCATaxa, final int[] nTaxonCount) { - if (node.isLeaf()) { - nTaxonCount[0]++; - if (MRCATaxonSet.get(node.getNr())) { - return 1; - } else { - return 0; - } - } else { - int iTaxons = traverse(node.getLeft(), MRCATaxonSet, nNrOfMRCATaxa, nTaxonCount); - final int nLeftTaxa = nTaxonCount[0]; - nTaxonCount[0] = 0; - if (node.getRight() != null) { - iTaxons += traverse(node.getRight(), MRCATaxonSet, nNrOfMRCATaxa, nTaxonCount); - final int nRightTaxa = nTaxonCount[0]; - nTaxonCount[0] = nLeftTaxa + nRightTaxa; - } - if (iTaxons == nrOfTaxa + 127) { - iTaxons++; - } - if (iTaxons == nNrOfMRCATaxa) { - // we are at the MRCA, return magic nr - return nrOfTaxa + 127; - } - return iTaxons; - } - } - - - @Override - public String[] getTaxaNames() { - if (m_sTaxaNames == null) { - final List sTaxa; - if (taxaInput.get() != null) { - sTaxa = taxaInput.get().getTaxaNames(); - } else { - sTaxa = m_taxonset.get().asStringList(); - } - m_sTaxaNames = sTaxa.toArray(new String[sTaxa.size()]); - } - return m_sTaxaNames; - } - - final private ArrayList nodeList = new ArrayList(); - private int activeNodeCount = 0; - -} diff --git a/src/beast/evolution/tree/ZeroBranchSATree.java b/src/beast/evolution/tree/ZeroBranchSATree.java deleted file mode 100644 index 804f9d9..0000000 --- a/src/beast/evolution/tree/ZeroBranchSATree.java +++ /dev/null @@ -1,142 +0,0 @@ -package beast.evolution.tree; - -import beast.core.StateNodeInitialiser; - -import java.io.PrintStream; -import java.util.List; - -/** - *@author Alexandra Gavryushkina - */ - -/** - * @deprecated Use Tree instead. Starting from v.1.1.1 Tree class supports sampled ancestors and - * all the classes that contain ZeroBranch in their names (and some others) are replaced by similar classes. - */ -@Deprecated -public class ZeroBranchSATree extends Tree { - - public boolean logWithZeroBranches=false; - - @Override - public void initAndValidate() { - if (m_initial.get() != null && !(this instanceof StateNodeInitialiser)) { - final Tree other = m_initial.get(); - root = other.root.copy(); - nodeCount = other.nodeCount; - internalNodeCount = other.internalNodeCount; - leafNodeCount = other.leafNodeCount; - } - - if (nodeCount < 0) { - if (m_taxonset.get() != null) { - // make a caterpillar - final List sTaxa = m_taxonset.get().asStringList(); - Node left = newNode(); - left.labelNr = 0; - left.height = 0; - left.setID(sTaxa.get(0)); - for (int i = 1; i < sTaxa.size(); i++) { - Node right = newNode(); - right.labelNr = i; - right.height = 0; - right.setID(sTaxa.get(i)); - Node parent = newNode(); - parent.labelNr = sTaxa.size() + i - 1; - parent.height = i; - left.parent = parent; - parent.setLeft(left); - right.parent = parent; - parent.setRight(right); - left = parent; - } - root = left; - leafNodeCount = sTaxa.size(); - nodeCount = leafNodeCount * 2 - 1; - internalNodeCount = leafNodeCount - 1; - - } else { - // make dummy tree with a single root node - root = newNode(); - root.labelNr = 0; - root.height = 0; - root.m_tree = this; - nodeCount = 1; - internalNodeCount = 0; - leafNodeCount = 1; - } - } - - if (nodeCount >= 0) { - initArrays(); - } - - if (leafNodeCount < 0) { - leafNodeCount = getLeafNodeCount(); - } - if (internalNodeCount < 0) { - internalNodeCount = getInternalNodeCount(); - } - - - processTraits(m_traitList.get()); - - // Ensure tree is compatible with time trait. - if (timeTraitSet != null) - adjustTreeNodeHeights(root); - - - } - - @Override - public int scale(double fScale) { - return scale(fScale,false); - } - - public int scale(double fScale, boolean scaleSNodes) { - ((ZeroBranchSANode)root).scale(fScale, scaleSNodes); - if (scaleSNodes) { - return getNodeCount() - getDirectAncestorNodeCount(); - } else { - return getInternalNodeCount() - getDirectAncestorNodeCount(); - } - } - - public int getDirectAncestorNodeCount() { - int directAncestorNodeCount = 0; - for (int i = 0; i < leafNodeCount; i++) { - if (((ZeroBranchSANode)this.getNode(i)).isDirectAncestor()) { - directAncestorNodeCount += 1; - } - } - return directAncestorNodeCount; - } - - @Override - public void log(long nSample, PrintStream out) { - ZeroBranchSATree tree = (ZeroBranchSATree) getCurrent(); - out.print("tree STATE_" + nSample + " = "); - final int[] dummy = new int[1]; - final String sNewick; - if (logWithZeroBranches) { - sNewick = ((ZeroBranchSANode)tree.getRoot()).toSortedNewickWithZeroBranches(dummy); - } else { - sNewick = tree.getRoot().toSortedNewick(dummy); - } - out.print(sNewick); - out.print(";"); - } - - public double scaleAndSlide(double scale) throws Exception { - int[] splitCount = new int[]{0}; - ((ZeroBranchSANode)root).scaleAndSlide(scale, splitCount); - System.out.println("tree becomes " + root.toShortNewick(false)); - double scaleContribution = (getInternalNodeCount() - getDirectAncestorNodeCount())* Math.log(scale); - if (scale > 1.0) { - return scaleContribution - splitCount[0]*Math.log(2.0); - } else { - return scaleContribution + splitCount[0]*Math.log(2.0); - } - } - -} diff --git a/src/beast/util/ClusterZBSATree.java b/src/beast/util/ClusterZBSATree.java index 316c419..ec0fcb1 100644 --- a/src/beast/util/ClusterZBSATree.java +++ b/src/beast/util/ClusterZBSATree.java @@ -10,8 +10,6 @@ import beast.evolution.alignment.distance.JukesCantorDistance; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSATree; - import java.text.DecimalFormat; import java.text.DecimalFormatSymbols; import java.util.*; @@ -31,7 +29,7 @@ "
o adjusted complete link " + "
o neighborjoining " + "
o neighborjoining2 - corrects tree for tip data, unlike plain neighborjoining") -public class ClusterZBSATree extends ZeroBranchSATree implements StateNodeInitialiser { +public class ClusterZBSATree extends Tree implements StateNodeInitialiser { final static String M_SINGLE = "single"; final static String M_AVERAGE = "average"; final static String M_COMPLETE = "complete"; @@ -88,9 +86,9 @@ public void initAndValidate() { parent.setNr(taxaNames.size() + i - 1); parent.setHeight(i); left.setParent(parent); - parent.setLeft(left); + parent.addChild(left); right.setParent(parent); - parent.setRight(right); + parent.addChild(right); left = parent; } root = left; diff --git a/src/beast/util/ZeroBranchSATreeParser.java b/src/beast/util/ZeroBranchSATreeParser.java index 7907857..1f49b7d 100644 --- a/src/beast/util/ZeroBranchSATreeParser.java +++ b/src/beast/util/ZeroBranchSATreeParser.java @@ -8,8 +8,6 @@ import beast.evolution.alignment.TaxonSet; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSATree; - import java.util.ArrayList; import java.util.List; import java.util.Map; @@ -22,7 +20,7 @@ @Description("Create beast.tree by parsing from a specification of a beast.tree (which is a fake SA tree) in Newick" + " format (includes parsing of any meta data in the Newick string).") -public class ZeroBranchSATreeParser extends ZeroBranchSATree implements StateNodeInitialiser { +public class ZeroBranchSATreeParser extends Tree implements StateNodeInitialiser { /** * default beast.tree branch length, used when that info is not in the Newick beast.tree @@ -125,12 +123,9 @@ public void initAndValidate() { */ List m_bTaxonIndexInUse = new ArrayList(); - public ZeroBranchSATreeParser() throws Exception { - nodeTypeInput.setValue("beast.evolution.tree.ZeroBranchSANode", this); - } + public ZeroBranchSATreeParser() throws Exception {} public ZeroBranchSATreeParser(Alignment alignment, String newick) throws Exception { - nodeTypeInput.setValue("beast.evolution.tree.ZeroBranchSANode", this); dataInput.setValue(alignment, this); newickInput.setValue(newick, this); initAndValidate(); @@ -150,7 +145,6 @@ public ZeroBranchSATreeParser(Alignment alignment, String newick) throws Excepti public ZeroBranchSATreeParser(List taxaNames, String newick, int offset) throws Exception { - nodeTypeInput.setValue("beast.evolution.tree.ZeroBranchSANode", this); if (taxaNames == null) { isLabelledNewickInput.setValue(true, this); } else { @@ -195,7 +189,6 @@ public ZeroBranchSATreeParser(String newick, boolean allowSingleChildNodes, boolean isLabeled, int offset) throws Exception { - nodeTypeInput.setValue("beast.evolution.tree.ZeroBranchSANode", this); newickInput.setValue(newick, this); isLabelledNewickInput.setValue(isLabeled, this); allowSingleChildInput.setValue(allowSingleChildNodes, this); diff --git a/src/test/beast/evolution/operators/SampledNodeDateRandomWalkerTest.java b/src/test/beast/evolution/operators/SampledNodeDateRandomWalkerTest.java new file mode 100644 index 0000000..0ef9c9f --- /dev/null +++ b/src/test/beast/evolution/operators/SampledNodeDateRandomWalkerTest.java @@ -0,0 +1,77 @@ +package test.beast.evolution.operators; + +import org.junit.Test; + +import beast.evolution.alignment.Taxon; +import beast.evolution.alignment.TaxonSet; +import beast.evolution.operators.SampledNodeDateRandomWalker; +import beast.evolution.tree.SamplingDate; +import beast.evolution.tree.Tree; +import beast.evolution.tree.TreeWOffset; +import beast.util.Randomizer; +import beast.util.ZeroBranchSATreeParser; +import junit.framework.TestCase; + +public class SampledNodeDateRandomWalkerTest extends TestCase { + + @Test + public void testHeightWOffset() throws Exception { + + Tree tree = new ZeroBranchSATreeParser("((t1:1.5,t2:0.5):0.5)3:0.0", true, true, 1); + TaxonSet tx = new TaxonSet(), txs = new TaxonSet(); + tx.initByName("taxon", new Taxon("t1"), "taxon", new Taxon("t2")); + tree.setInputValue("taxonset", tx); + + TreeWOffset treewoffset = new TreeWOffset(); + treewoffset.initByName("tree", tree, "offset", 25.0); + + SamplingDate sd1 = new SamplingDate(); + sd1.initByName("taxon", tree.getTaxonset().taxonsetInput.get().get(0), "lower", "24.0", "upper", "30.0"); + + txs.initByName("taxon", new Taxon("t1")); + SampledNodeDateRandomWalker op = new SampledNodeDateRandomWalker(); + op.initByName("tree", tree, "treeWOffset", treewoffset, "samplingDates", sd1, + "windowSize", 10.0, "weight", 1.0, "taxonset", txs); + + Randomizer.setSeed(2); + // test that tip t2 doesn't move + for(int i = 0; i < 25; i++) { + op.proposal(); + assertEquals(treewoffset.getHeightOfNode(1), 26.0); + } + } + + @Test + public void testSwitchWOffset() throws Exception { + + Tree tree = new ZeroBranchSATreeParser("((t1:1.5,t2:0.5):0.5)3:0.0", true, true, 1); + TaxonSet tx = new TaxonSet(), txs = new TaxonSet(); + tx.initByName("taxon", new Taxon("t1"), "taxon", new Taxon("t2")); + tree.setInputValue("taxonset", tx); + + TreeWOffset treewoffset = new TreeWOffset(); + treewoffset.initByName("tree", tree, "offset", 25.0); + + SamplingDate sd1 = new SamplingDate(); + sd1.initByName("taxon", tree.getTaxonset().taxonsetInput.get().get(0), "lower", "24.0", "upper", "30.0"); + + txs.initByName("taxon", new Taxon("t1")); + SampledNodeDateRandomWalker op = new SampledNodeDateRandomWalker(); + op.initByName("tree", tree, "treeWOffset", treewoffset, "samplingDates", sd1, + "windowSize", 10.0, "weight", 1.0, "taxonset", txs); + + // test that switching leaves happens correctly + Randomizer.setSeed(2); + op.proposal(); + assertEquals(treewoffset.getOffset(), 26.0); + assertEquals(treewoffset.getTree().getNode(1).getHeight(), 0.0); + assertTrue(treewoffset.getTree().getNode(0).getHeight() > 0); + + treewoffset.restore(); + assertEquals(treewoffset.getOffset(), 25.0); + assertEquals(treewoffset.getStoredHeightOfLeaf(0), 25.0); + assertTrue(treewoffset.getStoredHeightOfLeaf(1) > 25.0); + + } + +} diff --git a/src/test/beast/evolution/operators/TreeDimensionJumpTest.java b/src/test/beast/evolution/operators/TreeDimensionJumpTest.java index d8f3baa..7722768 100644 --- a/src/test/beast/evolution/operators/TreeDimensionJumpTest.java +++ b/src/test/beast/evolution/operators/TreeDimensionJumpTest.java @@ -3,7 +3,6 @@ import beast.evolution.operators.TreeDimensionJump; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; -import beast.evolution.tree.ZeroBranchSANode; import junit.framework.TestCase; import org.junit.Test; @@ -20,20 +19,20 @@ public void testOperator1() throws Exception { int taxaSize = 3; // make a caterpillar - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(0); left.setHeight(0.0); for (int i = 1; i < taxaSize; i++) { - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(i); right.setHeight(i); - Node parent = new ZeroBranchSANode(); + Node parent = new Node(); parent.setNr(taxaSize + i - 1); parent.setHeight(i); left.setParent(parent); - parent.setLeft(left); + parent.addChild(left); right.setParent(parent); - parent.setRight(right); + parent.addChild(right); left = parent; } tree = new Tree(left); @@ -56,20 +55,20 @@ public void testOperator2() throws Exception { int taxaSize = 3; // make a caterpillar - Node left = new ZeroBranchSANode(); + Node left = new Node(); left.setNr(0); left.setHeight(0.0); for (int i = 1; i < taxaSize; i++) { - Node right = new ZeroBranchSANode(); + Node right = new Node(); right.setNr(i); right.setHeight(i); - Node parent = new ZeroBranchSANode(); + Node parent = new Node(); parent.setNr(taxaSize + i - 1); parent.setHeight(i); left.setParent(parent); - parent.setLeft(left); + parent.addChild(left); right.setParent(parent); - parent.setRight(right); + parent.addChild(right); left = parent; } left.setHeight(left.getRight().getHeight()+2); diff --git a/src/test/beast/evolution/operators/ZeroBranchSATreeScalerTest.java b/src/test/beast/evolution/operators/ZeroBranchSATreeScalerTest.java deleted file mode 100644 index a18bd85..0000000 --- a/src/test/beast/evolution/operators/ZeroBranchSATreeScalerTest.java +++ /dev/null @@ -1,56 +0,0 @@ -package test.beast.evolution.operators; - -import beast.evolution.operators.ScaleOperatorForZeroBranchSATrees; -import beast.util.ZeroBranchSATreeParser; -import junit.framework.TestCase; -import org.junit.Test; - - -import java.util.ArrayList; - -/** - *@author Alexandra Gavryushkina - */ -public class ZeroBranchSATreeScalerTest extends TestCase { - - @Test - public void testOperator() throws Exception { - - String newick = "((((0:0.5)2:0.5,1:0.7):0.5)4:0.8,3:1.5):0.0"; - - ArrayList taxa = new ArrayList(); - for (int i=0; i<5; i++) { - taxa.add(Integer.toString(i)); - } - ZeroBranchSATreeParser tree = new ZeroBranchSATreeParser(taxa, newick, 0); - - double oldTreeHeight = tree.getRoot().getHeight(); - double oldLeftLeftHeight = tree.getRoot().getLeft().getLeft().getHeight(); - - ScaleOperatorForZeroBranchSATrees operator = new ScaleOperatorForZeroBranchSATrees(); - operator.initByName("tree", tree, "weight", 1.0); - operator.initByName("scaleFactor", 0.95); - operator.initAndValidate(); - - - double HR; - do { - HR = operator.proposal(); - } while (HR == Double.NEGATIVE_INFINITY); - - System.out.println(tree.getRoot().toShortNewick(false)); - - double scaler = tree.getRoot().getHeight()/oldTreeHeight; - double newLefLeftHeight = tree.getRoot().getLeft().getLeft().getHeight(); - - assertEquals(oldLeftLeftHeight, newLefLeftHeight/scaler, 1e-5); - boolean sampledNodeScaled = false; - for (int i=0; i cladeConstraints = new ArrayList<>(); - - CladeConstraint clade = new CladeConstraint(); - clade.setID("constraint1"); - TaxonSet taxonSetIn = new TaxonSet(); - taxonSetIn.initByName("taxon", A, "taxon", B, "taxon", C); - TaxonSet taxonSetOut = new TaxonSet(); - taxonSetOut.initByName("taxon", D); - clade.initByName( - "tree", tree, - "taxonsetIn", taxonSetIn, - "taxonsetOut", taxonSetOut - ); - cladeConstraints.add(clade); - - clade = new CladeConstraint(); - clade.setID("constraint2"); - taxonSetIn = new TaxonSet(); - taxonSetIn.initByName("taxon", A, "taxon", B); - clade.initByName( - "tree", tree, - "taxonsetIn", taxonSetIn - ); - cladeConstraints.add(clade); - - clade = new CladeConstraint(); - clade.setID("constraint3"); - taxonSetIn = new TaxonSet(); - taxonSetIn.initByName("taxon", D, "taxon", E); - clade.initByName( - "tree", tree, - "taxonsetIn", taxonSetIn - ); - cladeConstraints.add(clade); - - ConstantPopulation popFunc = new ConstantPopulation(); - popFunc.initByName("popSize", new RealParameter("1.0")); - - ZeroBranchSARandomTree sARandomTree; - for (int i=0; i<1000; i++) { - sARandomTree = new ZeroBranchSARandomTree(); - sARandomTree.initByName( - "nodetype", "beast.evolution.tree.ZeroBranchSANode", - "initial", tree, - "taxonset", tree.getTaxonset(), - "populationModel", popFunc, - "cladeConstraint", cladeConstraints - ); - - String sortedNewickTopology = TreeUtils.sortedNewickTopology(sARandomTree.getRoot(), true); -// System.out.println(sortedNewickTopology); - assert sortedNewickTopology.contentEquals(treeTopology); - } - } - - @Test - public void testTreeTopology() throws Exception { -// Randomizer.setSeed(777); - - ZeroBranchSATree tree = new ZeroBranchSATree(); - tree.initByName("taxonset", taxonSetAll, "trait", traitSet); - - // cladeConstraints - List cladeConstraints = new ArrayList<>(); - - CladeConstraint clade = new CladeConstraint(); - clade.setID("constraint1"); - TaxonSet taxonSetIn = new TaxonSet(); - taxonSetIn.initByName("taxon", A, "taxon", B, "taxon", C); -// TaxonSet taxonSetOut = new TaxonSet(); -// taxonSetOut.initByName("taxon", D); - clade.initByName( - "tree", tree, - "taxonsetIn", taxonSetIn, - //"taxonsetOut", taxonSetOut, - "stronglyMonophyletic", true - ); - cladeConstraints.add(clade); - - clade = new CladeConstraint(); - clade.setID("constraint2"); - clade.initByName( - "tree", tree, - "taxonsetIn", taxonSetAll, - "stronglyMonophyletic", true - ); - cladeConstraints.add(clade); - - ConstantPopulation popFunc = new ConstantPopulation(); - popFunc.initByName("popSize", new RealParameter("1.0")); - - ZeroBranchSARandomTree sARandomTree; - for (int i=0; i<1000; i++) { - sARandomTree = new ZeroBranchSARandomTree(); - sARandomTree.initByName( - "nodetype", "beast.evolution.tree.ZeroBranchSANode", - "initial", tree, - "taxonset", tree.getTaxonset(), - "populationModel", popFunc, - "cladeConstraint", cladeConstraints - ); - - String sortedNewickTopology = TreeUtils.sortedNewickTopology(sARandomTree.getRoot(), true); - //System.out.println(sortedNewickTopology); - - boolean contain = false; - for (String topology:topologies) { - contain = (contain || sortedNewickTopology.contentEquals(topology)); - } - assert contain; - } - } -}