From 1b67e9acba572d55c9e8839d40d3463983cb5b48 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Tue, 23 Mar 2021 16:52:58 +0100 Subject: [PATCH 1/8] Move log computation earlier In linguistics, we tend to measure times in years, so heights are in the order of magnitude of 100s to 1000s, not the 0.1My often found in biology. In order to avoid overflow errors, the log should be taken ASAP, and in fact this does not make the computation more expensive. --- .../operators/LeafToSampledAncestorJump.java | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/beast/evolution/operators/LeafToSampledAncestorJump.java b/src/beast/evolution/operators/LeafToSampledAncestorJump.java index e9f7642..ad742f1 100644 --- a/src/beast/evolution/operators/LeafToSampledAncestorJump.java +++ b/src/beast/evolution/operators/LeafToSampledAncestorJump.java @@ -32,7 +32,7 @@ public void initAndValidate() { @Override public double proposal() { - double newHeight, newRange, oldRange; + double newHeight, logNewRange, logOldRange; int categoryCount = 1; if (categoriesInput.get() != null) { @@ -47,14 +47,15 @@ public double proposal() { Node parent = leaf.getParent(); if (leaf.isDirectAncestor()) { - oldRange = (double) 1; + logOldRange = (double) 0; if (parent.isRoot()) { final double randomNumber = Randomizer.nextExponential(1); newHeight = parent.getHeight() + randomNumber; - newRange = Math.exp(randomNumber); + logNewRange = randomNumber; } else { - newRange = parent.getParent().getHeight() - parent.getHeight(); + double newRange = parent.getParent().getHeight() - parent.getHeight(); newHeight = parent.getHeight() + Randomizer.nextDouble() * newRange; + logNewRange = Math.log(newRange); } if (categoriesInput.get() != null) { @@ -63,15 +64,15 @@ public double proposal() { categoriesInput.get().setValue(index, newValue); } } else { - newRange = (double) 1; + logNewRange = (double) 0; //make sure that the branch where a new sampled node to appear is not above that sampled node if (getOtherChild(parent, leaf).getHeight() >= leaf.getHeight()) { return Double.NEGATIVE_INFINITY; } if (parent.isRoot()) { - oldRange = Math.exp(parent.getHeight() - leaf.getHeight()); + logOldRange = parent.getHeight() - leaf.getHeight(); } else { - oldRange = parent.getParent().getHeight() - leaf.getHeight(); + logOldRange = Math.log(parent.getParent().getHeight() - leaf.getHeight()); } newHeight = leaf.getHeight(); if (categoriesInput.get() != null) { @@ -86,6 +87,6 @@ public double proposal() { return Double.NEGATIVE_INFINITY; } - return Math.log(newRange/oldRange); + return logNewRange - logOldRange; } } From 5295a3a99cc77cd72a1796d7108cec05b684a698 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Tue, 23 Mar 2021 22:17:01 +0100 Subject: [PATCH 2/8] Add option to sample a subset of leaves --- .../operators/LeafToSampledAncestorJump.java | 34 +++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/src/beast/evolution/operators/LeafToSampledAncestorJump.java b/src/beast/evolution/operators/LeafToSampledAncestorJump.java index ad742f1..a98e989 100644 --- a/src/beast/evolution/operators/LeafToSampledAncestorJump.java +++ b/src/beast/evolution/operators/LeafToSampledAncestorJump.java @@ -1,9 +1,13 @@ package beast.evolution.operators; +import java.util.ArrayList; +import java.util.List; + import beast.core.Description; import beast.core.Input; import beast.core.parameter.IntegerParameter; import beast.core.parameter.RealParameter; +import beast.evolution.alignment.Taxon; import beast.evolution.tree.Node; import beast.evolution.tree.Tree; import beast.util.Randomizer; @@ -24,9 +28,35 @@ public class LeafToSampledAncestorJump extends TreeOperator { public Input rInput = new Input("removalProbability", "The probability of an individual to be removed from the process immediately after the sampling"); + public Input> sampledTaxa = + new Input>("sampledTaxa", "Taxa that this operator should be allowed to let jump between sampled ancestor and leaf. Default: All non-recent leaves."); + protected List validLeaves; + @Override public void initAndValidate() { + if (sampledTaxa.get() == null) { + validLeaves = new ArrayList(treeInput.get().getLeafNodeCount()); + for (Node leaf: treeInput.get().getExternalNodes()) { + if (leaf.getHeight() > 1e-6) { + validLeaves.add(leaf.getNr()); + } + } + } else { + List taxa = sampledTaxa.get(); + List taxaNames = new ArrayList(taxa.size()); + for (Taxon taxon: taxa) { + taxaNames.add(taxon.getID()); + } + validLeaves = new ArrayList(taxa.size()); + Integer i = 0; + for (String leaf: treeInput.get().getTaxaNames()) { + if (taxaNames.contains(leaf)) { + validLeaves.add(i); + } + i += 1; + } + } } @Override @@ -41,9 +71,9 @@ public double proposal() { Tree tree = treeInput.get(); - int leafNodeCount = tree.getLeafNodeCount(); + int leafNodeCount = validLeaves.size(); - Node leaf = tree.getNode(Randomizer.nextInt(leafNodeCount)); + Node leaf = tree.getNode(validLeaves.get(Randomizer.nextInt(leafNodeCount))); Node parent = leaf.getParent(); if (leaf.isDirectAncestor()) { From c1a359846975f0da41060b068a8978ac76e53738 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Tue, 23 Mar 2021 23:28:12 +0100 Subject: [PATCH 3/8] Fix initialization --- examples/bears-not-all-leaves.xml | 267 ++++++++++++++++++ .../operators/LeafToSampledAncestorJump.java | 12 +- 2 files changed, 276 insertions(+), 3 deletions(-) create mode 100644 examples/bears-not-all-leaves.xml diff --git a/examples/bears-not-all-leaves.xml b/examples/bears-not-all-leaves.xml new file mode 100644 index 0000000..2a5d225 --- /dev/null +++ b/examples/bears-not-all-leaves.xml @@ -0,0 +1,267 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +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 + + + + + + + Canis_lupus_0.0=0.0, +Phoca_largha_0.0=0.0, +Ailuropoda_melanoleuca_0.0=0.0, +Tremarctos_ornatus_0.0=0.0, +Helarctos_malayanus_0.0=0.0, +Melursus_ursinus_0.0=0.0, +Ursus_americanus_0.0=0.0, +Ursus_arctos_0.0=0.0, +Ursus_thibetanus_0.0=0.0, +Ursus_maritimus_0.0=0.0, +Hesperocyon_gregarius_39.07=39.07, +Caedocyon_tedfordi_25.88=25.88, +Osbornodon_sesnoni_30.08=30.08, +Cormocyon_copei_26.14=26.14, +Borophagus_diversidens_4.284=4.284, +Enaliarctos_tedfordi_27.11=27.11, +Proneotherium_repenningi_17.92=17.92, +Leptophoca_lenis_14.99=14.99, +Acrophoca_6.695=6.695, +Phoca_vitulina_0.805=0.805, +Parictis_montanus_36.6=36.6, +Zaragocyon_daamsi_21.86=21.86, +Ballusia_elmensis_14.01=14.01, +Ursavus_primaevus_14.41=14.41, +Ursavus_brevihinus_16.2=16.2, +Indarctos_vireti_8.68=8.68, +Indarctos_arctoides_9.54=9.54, +Indarctos_punjabiensis_4.996=4.996, +Ailurarctos_lufengensis_7.652=7.652, +Agriarctos_spp_5.006=5.006, +Kretzoiarctos_beatrix_11.69=11.69, +Arctodus_simus_0.487=0.487, +Ursus_abstrusus_4.27=4.27, +Ursus_spelaeus_0.054=0.054 + + + + + + 1.0 + 0.1 + 1 + 100.0 + 1.0 + 0.5 + 0.5 + + + + + 1.0 + + + + + + + 0.0 + 1.0 + + + + + + + + + + + + + + + + + + + 0.5396 + 0.3819 + + + + + + + + + + + + + + + + + 1.0 + 1.0 + 0.0 + + + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/beast/evolution/operators/LeafToSampledAncestorJump.java b/src/beast/evolution/operators/LeafToSampledAncestorJump.java index a98e989..d768656 100644 --- a/src/beast/evolution/operators/LeafToSampledAncestorJump.java +++ b/src/beast/evolution/operators/LeafToSampledAncestorJump.java @@ -1,6 +1,7 @@ package beast.evolution.operators; import java.util.ArrayList; +import java.util.Arrays; import java.util.List; import beast.core.Description; @@ -29,13 +30,16 @@ public class LeafToSampledAncestorJump extends TreeOperator { public Input rInput = new Input("removalProbability", "The probability of an individual to be removed from the process immediately after the sampling"); public Input> sampledTaxa = - new Input>("sampledTaxa", "Taxa that this operator should be allowed to let jump between sampled ancestor and leaf. Default: All non-recent leaves."); + new Input>( + "sampledTaxa", + "Taxa that this operator should be allowed to let jump between sampled ancestor and leaf. Default: All non-recent leaves.", + new ArrayList<>()); - protected List validLeaves; + protected List validLeaves = new ArrayList(); @Override public void initAndValidate() { - if (sampledTaxa.get() == null) { + if (sampledTaxa.get().size() == 0) { validLeaves = new ArrayList(treeInput.get().getLeafNodeCount()); for (Node leaf: treeInput.get().getExternalNodes()) { if (leaf.getHeight() > 1e-6) { @@ -57,6 +61,8 @@ public void initAndValidate() { i += 1; } } + // System.out.println("Nodes to be jumped:"); + // System.out.println(Arrays.toString(validLeaves.toArray())); } @Override From 5c93ce386935347b259c51683f899235d0edb2ab Mon Sep 17 00:00:00 2001 From: Gereon Kaiping Date: Wed, 21 Jul 2021 17:33:45 +0200 Subject: [PATCH 4/8] Modify tree, then set root. This should fix one of the symptoms of https://github.com/CompEvol/beast2/issues/984 and also https://github.com/CompEvol/sampled-ancestors/issues/16 --- src/beast/evolution/operators/SAWilsonBalding.java | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/beast/evolution/operators/SAWilsonBalding.java b/src/beast/evolution/operators/SAWilsonBalding.java index 7b50fb3..9cbc86f 100644 --- a/src/beast/evolution/operators/SAWilsonBalding.java +++ b/src/beast/evolution/operators/SAWilsonBalding.java @@ -157,11 +157,12 @@ public double proposal() { jP.removeChild(j); // remove jP.addChild(iP); // add jP.makeDirty(Tree.IS_FILTHY); + iP.addChild(j); } else { iP.setParent(null); // completely remove - tree.setRootOnly(iP); + iP.addChild(j); + tree.setRoot(iP); } - iP.addChild(j); iP.makeDirty(Tree.IS_FILTHY); j.makeDirty(Tree.IS_FILTHY); } From db0fe3a6bb61c8e5a6bcf3d820233f10d9d918e1 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Wed, 21 Jul 2021 19:05:20 +0200 Subject: [PATCH 5/8] Add random tree that can start with sampled ancestors --- .../evolution/tree/RandomTreeWithSA.java | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 src/beast/evolution/tree/RandomTreeWithSA.java diff --git a/src/beast/evolution/tree/RandomTreeWithSA.java b/src/beast/evolution/tree/RandomTreeWithSA.java new file mode 100644 index 0000000..534b6cd --- /dev/null +++ b/src/beast/evolution/tree/RandomTreeWithSA.java @@ -0,0 +1,60 @@ +package beast.evolution.tree; + +import java.util.List; +import java.util.Set; + +import beast.core.Input; +import beast.core.util.Log; +import beast.evolution.alignment.TaxonSet; + +public class RandomTreeWithSA extends RandomTree { + + final public Input startAsSampledAncestors = new Input<>("sampledAncestor", + "Taxa to start as sampled ancestors"); + + @Override + public void initAndValidate() { + if (!taxaInput.get().getTaxaNames().containsAll(startAsSampledAncestors.get().getTaxaNames())) { + Set remaining = startAsSampledAncestors.get().getTaxaNames(); + remaining.removeAll(taxaInput.get().getTaxaNames()); + Log.warning("WARNING: Not all sampledAncestor taxa were part of the alignment taxa. Missing taxa " + + remaining.toString() + " will not be forced to be sampled ancestors."); + } + super.initAndValidate(); + } + + @Override + public void initStateNodes() { + super.initStateNodes(); + + List sa = startAsSampledAncestors.get().asStringList(); + + for (Node leaf : root.getAllLeafNodes()) { + + if (!sa.contains(this.getTaxonId(leaf))) { + continue; + } + Node parent = leaf.getParent(); + + double newHeight = leaf.getHeight(); + parent.setHeight(newHeight); + // setHeight does not check temporal order of nodes, so we nudge all children + // that are older than their parent down a bit. This can mess with other + // assumptions. + for (Node node : parent.getAllChildNodesAndSelf()) { + if (node.getHeight() > node.parent.getHeight()) { + if (sa.contains(getTaxonId(node))) { + // Keep sampled ancestors as sampled ancestors. Yes, this is a bit redundant, + // but at least it should be correct. The alternative is to recursively walk the + // tree in a somewhat erratic order. + node.setHeight(node.parent.getHeight()); + } else { + // Put everything else down a notch. + node.setHeight(node.parent.getHeight() - 1e-9); + } + } + } + assert leaf.isDirectAncestor(); + } + } +} From 59e958352dc00d2dedc93bf7a67fb904cc5b8b25 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Wed, 21 Jul 2021 22:57:48 +0200 Subject: [PATCH 6/8] Fix to corner cases --- .../evolution/tree/RandomTreeWithSA.java | 24 +++++++++++++++---- ...rators.SampledNodeDateRandomWalkerTest.txt | 5 ++++ ...lution.operators.TreeDimensionJumpTest.txt | 13 ++++++++++ ...ution.speciation.SABirthDeathModelTest.txt | 6 +++++ ...t.beast.evolution.tree.TreeWOffsetTest.txt | 4 ++++ 5 files changed, 48 insertions(+), 4 deletions(-) create mode 100644 test-reports/TEST-test.beast.evolution.operators.SampledNodeDateRandomWalkerTest.txt create mode 100644 test-reports/TEST-test.beast.evolution.operators.TreeDimensionJumpTest.txt create mode 100644 test-reports/TEST-test.beast.evolution.speciation.SABirthDeathModelTest.txt create mode 100644 test-reports/TEST-test.beast.evolution.tree.TreeWOffsetTest.txt diff --git a/src/beast/evolution/tree/RandomTreeWithSA.java b/src/beast/evolution/tree/RandomTreeWithSA.java index 534b6cd..6037e37 100644 --- a/src/beast/evolution/tree/RandomTreeWithSA.java +++ b/src/beast/evolution/tree/RandomTreeWithSA.java @@ -1,20 +1,28 @@ package beast.evolution.tree; +import java.util.LinkedHashSet; import java.util.List; import java.util.Set; import beast.core.Input; +import beast.core.Input.Validate; import beast.core.util.Log; import beast.evolution.alignment.TaxonSet; public class RandomTreeWithSA extends RandomTree { final public Input startAsSampledAncestors = new Input<>("sampledAncestor", - "Taxa to start as sampled ancestors"); + "Taxa to start as sampled ancestors", Validate.REQUIRED); @Override public void initAndValidate() { - if (!taxaInput.get().getTaxaNames().containsAll(startAsSampledAncestors.get().getTaxaNames())) { + taxa = new LinkedHashSet<>(); + if (taxaInput.get() != null) { + taxa.addAll(taxaInput.get().getTaxaNames()); + } else { + taxa.addAll(m_taxonset.get().asStringList()); + } + if (!taxa.containsAll(startAsSampledAncestors.get().getTaxaNames())) { Set remaining = startAsSampledAncestors.get().getTaxaNames(); remaining.removeAll(taxaInput.get().getTaxaNames()); Log.warning("WARNING: Not all sampledAncestor taxa were part of the alignment taxa. Missing taxa " @@ -31,7 +39,7 @@ public void initStateNodes() { for (Node leaf : root.getAllLeafNodes()) { - if (!sa.contains(this.getTaxonId(leaf))) { + if (!sa.contains(getTaxonId(leaf))) { continue; } Node parent = leaf.getParent(); @@ -43,7 +51,11 @@ public void initStateNodes() { // assumptions. for (Node node : parent.getAllChildNodesAndSelf()) { if (node.getHeight() > node.parent.getHeight()) { - if (sa.contains(getTaxonId(node))) { + if (!node.isLeaf()) { + // Not a leaf, so cannot be a sampled ancestor which needs to be handled with + // care + node.setHeight(node.parent.getHeight() - 1e-9); + } else if (sa.contains(getTaxonId(node))) { // Keep sampled ancestors as sampled ancestors. Yes, this is a bit redundant, // but at least it should be correct. The alternative is to recursively walk the // tree in a somewhat erratic order. @@ -56,5 +68,9 @@ public void initStateNodes() { } assert leaf.isDirectAncestor(); } + + if (m_initial.get() != null) { + m_initial.get().assignFromWithoutID(this); + } } } diff --git a/test-reports/TEST-test.beast.evolution.operators.SampledNodeDateRandomWalkerTest.txt b/test-reports/TEST-test.beast.evolution.operators.SampledNodeDateRandomWalkerTest.txt new file mode 100644 index 0000000..bf03ecd --- /dev/null +++ b/test-reports/TEST-test.beast.evolution.operators.SampledNodeDateRandomWalkerTest.txt @@ -0,0 +1,5 @@ +Testsuite: test.beast.evolution.operators.SampledNodeDateRandomWalkerTest +Tests run: 2, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.141 sec + +Testcase: testSwitchWOffset took 0.083 sec +Testcase: testHeightWOffset took 0.025 sec diff --git a/test-reports/TEST-test.beast.evolution.operators.TreeDimensionJumpTest.txt b/test-reports/TEST-test.beast.evolution.operators.TreeDimensionJumpTest.txt new file mode 100644 index 0000000..e5b5eca --- /dev/null +++ b/test-reports/TEST-test.beast.evolution.operators.TreeDimensionJumpTest.txt @@ -0,0 +1,13 @@ +Testsuite: test.beast.evolution.operators.TreeDimensionJumpTest +Tests run: 2, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.079 sec +------------- Standard Output --------------- +Tree was = ((0:1.0,1:0.0):3.0,2:2.0):0.0 +Proposed tree = ((0:1.5583203577630607,1:0.5583203577630607):2.4416796422369393,2:2.0):0.0 +Log Hastings ratio = 1.0986122886681098 +Tree was = ((0:1.0,1:0.0):1.0,2:0.0):0.0 +Proposed tree = ((0:1.0,1:0.0):1.0,2:0.0):0.0 +Log Hastings ratio = -Infinity +------------- ---------------- --------------- + +Testcase: testOperator2 took 0.051 sec +Testcase: testOperator1 took 0.005 sec diff --git a/test-reports/TEST-test.beast.evolution.speciation.SABirthDeathModelTest.txt b/test-reports/TEST-test.beast.evolution.speciation.SABirthDeathModelTest.txt new file mode 100644 index 0000000..1301faf --- /dev/null +++ b/test-reports/TEST-test.beast.evolution.speciation.SABirthDeathModelTest.txt @@ -0,0 +1,6 @@ +Testsuite: test.beast.evolution.speciation.SABirthDeathModelTest +Tests run: 3, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.15 sec + +Testcase: testLikelihoodCalculationSimple took 0.095 sec +Testcase: testLikelihoodCalculation1 took 0.011 sec +Testcase: testLikelihoodCalculation2 took 0.026 sec diff --git a/test-reports/TEST-test.beast.evolution.tree.TreeWOffsetTest.txt b/test-reports/TEST-test.beast.evolution.tree.TreeWOffsetTest.txt new file mode 100644 index 0000000..7496784 --- /dev/null +++ b/test-reports/TEST-test.beast.evolution.tree.TreeWOffsetTest.txt @@ -0,0 +1,4 @@ +Testsuite: test.beast.evolution.tree.TreeWOffsetTest +Tests run: 1, Failures: 0, Errors: 0, Skipped: 0, Time elapsed: 0.067 sec + +Testcase: testTree took 0.049 sec From 6b2b5e142d01595c310943719a14561d123b907d Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Wed, 21 Jul 2021 23:29:02 +0200 Subject: [PATCH 7/8] Add default value that reduces RTwSA to RT --- .../evolution/tree/RandomTreeWithSA.java | 24 ++++++++++++++----- 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/beast/evolution/tree/RandomTreeWithSA.java b/src/beast/evolution/tree/RandomTreeWithSA.java index 6037e37..73614a6 100644 --- a/src/beast/evolution/tree/RandomTreeWithSA.java +++ b/src/beast/evolution/tree/RandomTreeWithSA.java @@ -1,5 +1,7 @@ package beast.evolution.tree; +import java.util.ArrayList; +import java.util.HashSet; import java.util.LinkedHashSet; import java.util.List; import java.util.Set; @@ -7,12 +9,13 @@ import beast.core.Input; import beast.core.Input.Validate; import beast.core.util.Log; +import beast.evolution.alignment.Taxon; import beast.evolution.alignment.TaxonSet; public class RandomTreeWithSA extends RandomTree { - final public Input startAsSampledAncestors = new Input<>("sampledAncestor", - "Taxa to start as sampled ancestors", Validate.REQUIRED); + final public Input> startAsSampledAncestors = new Input<>("sampledAncestor", + "Taxa to start as sampled ancestors", new ArrayList()); @Override public void initAndValidate() { @@ -22,9 +25,14 @@ public void initAndValidate() { } else { taxa.addAll(m_taxonset.get().asStringList()); } - if (!taxa.containsAll(startAsSampledAncestors.get().getTaxaNames())) { - Set remaining = startAsSampledAncestors.get().getTaxaNames(); - remaining.removeAll(taxaInput.get().getTaxaNames()); + + Set startTaxa = new LinkedHashSet<>(); + for (Taxon taxon: startAsSampledAncestors.get()) { + startTaxa.add(taxon.getID()); + } + if (!taxa.containsAll(startTaxa)) { + Set remaining = startTaxa; + remaining.removeAll(taxa); Log.warning("WARNING: Not all sampledAncestor taxa were part of the alignment taxa. Missing taxa " + remaining.toString() + " will not be forced to be sampled ancestors."); } @@ -35,7 +43,11 @@ public void initAndValidate() { public void initStateNodes() { super.initStateNodes(); - List sa = startAsSampledAncestors.get().asStringList(); + List sa = new ArrayList(); + + for (Taxon taxon: startAsSampledAncestors.get()) { + sa.add(taxon.getID()); + } for (Node leaf : root.getAllLeafNodes()) { From 0b2d58e98d33fe0fee52771bafa9d87d5cc77212 Mon Sep 17 00:00:00 2001 From: "Gereon A. Kaiping" Date: Thu, 22 Jul 2021 11:12:59 +0200 Subject: [PATCH 8/8] Keep Sampled Ancestors as such. The SAWilsonBalding operator was permitted to turn sampled ancestors into normal tips. For an analysis where certain tips should be kept sampled ancestors, add a Boolean input (default: false) that aborts the operator proposal when it is about to do this change. --- .../evolution/operators/SAWilsonBalding.java | 352 +++++++++--------- 1 file changed, 180 insertions(+), 172 deletions(-) diff --git a/src/beast/evolution/operators/SAWilsonBalding.java b/src/beast/evolution/operators/SAWilsonBalding.java index 9cbc86f..7228be1 100644 --- a/src/beast/evolution/operators/SAWilsonBalding.java +++ b/src/beast/evolution/operators/SAWilsonBalding.java @@ -7,180 +7,188 @@ import beast.util.Randomizer; /** - *@author Alexandra Gavryushkina + * @author Alexandra Gavryushkina */ public class SAWilsonBalding 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() || 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 (!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 (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 - 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 (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); - iP.addChild(j); - } else { - iP.setParent(null); // completely remove - iP.addChild(j); - tree.setRoot(iP); - } - 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) && (tree.getDirectAncestorNodeCount() > 0 && rInput.get().getValue() == 1)) { - return Double.NEGATIVE_INFINITY; - } - - newDimension = nodeCount - tree.getDirectAncestorNodeCount() - 1; - DimensionCoefficient = (double) oldDimension / newDimension; - - fHastingsRatio = Math.abs(DimensionCoefficient * newRange / oldRange); - - return Math.log(fHastingsRatio); - - } - + public Input rInput = new Input("removalProbability", + "The probability of an individual to become noninfectious immediately after the sampling"); + public Input keepSAInput = new Input("keepSA", + "Avoid moving the descendant of a sampled ancestor away, turning it into a non-SA tip", false); + + @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; + + // There have been issues with this assumption in this package, so let's assert + // it explicitly. + assert tree.getNode(nodeCount - 1).isRoot(); + + do { + i = tree.getNode(Randomizer.nextInt(nodeCount - 1)); + } while (i.isDirectAncestor()); + + Node iP = i.getParent(); + Node CiP = getOtherChild(iP, i); + if (keepSAInput.get() && CiP.isDirectAncestor()) { + return Double.NEGATIVE_INFINITY; + } + + // 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 + " external 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 (!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 (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 - 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 (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); + iP.addChild(j); + } else { + iP.setParent(null); // completely remove + iP.addChild(j); + tree.setRoot(iP); + } + 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) && (tree.getDirectAncestorNodeCount() > 0 && rInput.get().getValue() == 1)) { + return Double.NEGATIVE_INFINITY; + } + + newDimension = nodeCount - tree.getDirectAncestorNodeCount() - 1; + DimensionCoefficient = (double) oldDimension / newDimension; + + fHastingsRatio = Math.abs(DimensionCoefficient * newRange / oldRange); + + return Math.log(fHastingsRatio); + + } }