diff --git a/beast-base/src/main/java/beast/base/evolution/RateStatistic.java b/beast-base/src/main/java/beast/base/evolution/RateStatistic.java index 1880606..e529226 100644 --- a/beast-base/src/main/java/beast/base/evolution/RateStatistic.java +++ b/beast-base/src/main/java/beast/base/evolution/RateStatistic.java @@ -26,13 +26,7 @@ package beast.base.evolution; -import java.io.PrintStream; - -import beast.base.core.BEASTObject; -import beast.base.core.Description; -import beast.base.core.Function; -import beast.base.core.Input; -import beast.base.core.Loggable; +import beast.base.core.*; import beast.base.core.Input.Validate; import beast.base.evolution.branchratemodel.BranchRateModel; import beast.base.evolution.likelihood.GenericTreeLikelihood; @@ -40,9 +34,13 @@ import beast.base.evolution.tree.Tree; import beast.base.util.DiscreteStatistics; +import java.io.PrintStream; - +/** + * @deprecated use {@link beast.base.spec.evolution.RateStatistic} instead + */ +@Deprecated @Description("A statistic that tracks the mean, variance and coefficent of variation of rates. " + "It has three dimensions, one for each statistic.") public class RateStatistic extends BEASTObject implements Loggable, Function { diff --git a/beast-base/src/main/java/beast/base/spec/evolution/RateStatistic.java b/beast-base/src/main/java/beast/base/spec/evolution/RateStatistic.java new file mode 100644 index 0000000..9ce7383 --- /dev/null +++ b/beast-base/src/main/java/beast/base/spec/evolution/RateStatistic.java @@ -0,0 +1,221 @@ +/* + * RateStatistic.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.base.spec.evolution; + + +import beast.base.core.BEASTObject; +import beast.base.core.Description; +import beast.base.core.Input; +import beast.base.core.Input.Validate; +import beast.base.core.Loggable; +import beast.base.evolution.tree.Node; +import beast.base.evolution.tree.Tree; +import beast.base.spec.domain.NonNegativeReal; +import beast.base.spec.evolution.branchratemodel.Base; +import beast.base.spec.evolution.likelihood.GenericTreeLikelihood; +import beast.base.spec.type.RealVector; +import beast.base.util.DiscreteStatistics; + +import java.io.PrintStream; +import java.util.Arrays; +import java.util.List; +import java.util.stream.Collectors; + + +@Description("A statistic that tracks the mean, variance and coefficent of variation of rates. " + + "It has three dimensions, one for each statistic.") +public class RateStatistic extends BEASTObject implements Loggable, RealVector { + + final public Input likelihoodInput = new Input<>("treeLikelihood", "TreeLikelihood containing branch rate model that provides rates for a tree"); + final public Input branchRateModelInput = new Input<>("branchratemodel", "model that provides rates for a tree", Validate.XOR, likelihoodInput); + final public Input treeInput = new Input<>("tree", "tree for which the rates apply", Validate.REQUIRED); + final public Input internalInput = new Input<>("internal", "consider internal nodes, default true", true); + final public Input externalInput = new Input<>("external", "consider external nodes, default true", true); + + private Tree tree = null; + private Base branchRateModel = null; + private boolean internal = true; + private boolean external = true; + + // array index + final static int MEAN = 0; + final static int VARIANCE = 1; + final static int COEFFICIENT_OF_VARIATION = 2; + + @Override + public void initAndValidate() { + tree = treeInput.get(); + branchRateModel = branchRateModelInput.get(); + if (branchRateModel == null) { + branchRateModel = likelihoodInput.get().branchRateModelInput.get(); + } + this.internal = internalInput.get(); + this.external = externalInput.get(); + } + + /** + * calculate the three statistics from scratch * + */ + public double[] calcValues() { + int length = 0; + int offset = 0; + + final int nrOfLeafs = tree.getLeafNodeCount(); + + if (external) { + length += nrOfLeafs; + } + if (internal) { + length += tree.getInternalNodeCount() - 1; + } + + final double[] rates = new double[length]; + // need those only for mean + final double[] branchLengths = new double[length]; + + final Node[] nodes = tree.getNodesAsArray(); + + /** handle leaf nodes **/ + if (external) { + for (int i = 0; i < nrOfLeafs; i++) { + final Node child = nodes[i]; + final Node parent = child.getParent(); + branchLengths[i] = parent.getHeight() - child.getHeight(); + rates[i] = branchRateModel.getRateForBranch(child); + } + offset = nrOfLeafs; + } + + /** handle internal nodes **/ + if (internal) { + final int n = tree.getNodeCount(); + int k = offset; + for (int i = nrOfLeafs; i < n; i++) { + final Node child = nodes[i]; + if (!child.isRoot()) { + final Node parent = child.getParent(); + branchLengths[k] = parent.getHeight() - child.getHeight(); + rates[k] = branchRateModel.getRateForBranch(child); + k++; + } + } + } + + final double[] values = new double[3]; + double totalWeightedRate = 0.0; + double totalTreeLength = 0.0; + for (int i = 0; i < rates.length; i++) { + totalWeightedRate += rates[i] * branchLengths[i]; + totalTreeLength += branchLengths[i]; + } + values[MEAN] = totalWeightedRate / totalTreeLength; + // compute mean/variance once and reuse: DiscreteStatistics.variance(rates) alone + // would recompute mean(rates) internally, duplicating both passes done here + final double mean = DiscreteStatistics.mean(rates); + values[VARIANCE] = DiscreteStatistics.variance(rates, mean); + values[COEFFICIENT_OF_VARIATION] = Math.sqrt(values[VARIANCE]) / mean; + return values; + } + + + /** + * Valuable implementation * + */ + + @Deprecated + public int getDimension() { + return size(); + } + + @Deprecated + public double getArrayValue() { + return get(0); + } + + @Deprecated + public double getArrayValue(final int dim) { + return get(dim); + } + + // new API + + @Override + public int size() { + return 3; + } + + @Override + public double get(int i) { + if (i < 0 || i >= size()) { + throw new IllegalArgumentException(); + } + return calcValues()[i]; + } + + @Override + public List getElements() { + return Arrays.stream(calcValues()).boxed().collect(Collectors.toList()); + } + + /** + * Non-negative rather than positive: the variance and the coefficient of + * variation are exactly 0 whenever every branch carries the same rate, which + * is the normal state of affairs under a strict clock. {@code PositiveReal} + * excludes 0 ({@code lowerInclusive()} is false), so declaring it here makes + * {@link #isValid()} report false for a perfectly ordinary analysis. + */ + @Override + public NonNegativeReal getDomain() { + return NonNegativeReal.INSTANCE; + } + + /** + * Loggable implementation * + */ + + @Override + public void init(final PrintStream out) { + String id = getID(); + if (id == null) { + id = ""; + } + out.print(id + ".mean\t" + id + ".variance\t" + id + ".coefficientOfVariation\t"); + } + + + @Override + public void log(final long sample, final PrintStream out) { + final double[] values = calcValues(); + out.print(values[0] + "\t" + values[1] + "\t" + values[2] + "\t"); + } + + + @Override + public void close(final PrintStream out) { + // nothing to do + } + +} diff --git a/beast-base/src/main/java/module-info.java b/beast-base/src/main/java/module-info.java index 49ef4f2..aa657e1 100644 --- a/beast-base/src/main/java/module-info.java +++ b/beast-base/src/main/java/module-info.java @@ -281,6 +281,7 @@ beast.base.spec.inference.distribution.TruncatedReal, beast.base.spec.inference.distribution.TruncatedInt, beast.base.spec.FunctionOfTensor, + beast.base.spec.evolution.RateStatistic, beast.base.spec.evolution.Sum, beast.base.spec.evolution.IntSum, beast.base.spec.evolution.TreeWithMetaDataLogger, diff --git a/beast-base/src/test/java/beast/base/spec/evolution/RateStatisticTest.java b/beast-base/src/test/java/beast/base/spec/evolution/RateStatisticTest.java new file mode 100644 index 0000000..4e402a5 --- /dev/null +++ b/beast-base/src/test/java/beast/base/spec/evolution/RateStatisticTest.java @@ -0,0 +1,122 @@ +package beast.base.spec.evolution; + +import beast.base.evolution.tree.Node; +import beast.base.evolution.tree.Tree; +import beast.base.spec.evolution.branchratemodel.Base; +import beast.base.spec.evolution.branchratemodel.StrictClockModel; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * Tests for the spec {@link RateStatistic}, in particular the domain it declares. + * + *

The statistic exposes three elements — mean, variance, coefficient of + * variation — through one {@code RealVector} domain, so that domain has to admit + * every value all three can legitimately take. Variance and coefficient of + * variation are exactly 0 whenever the branches share a rate, which is the + * ordinary state of affairs under a strict clock, so the domain has to include + * zero.

+ */ +public class RateStatisticTest { + + /** Rates keyed off node number, so each branch can differ. */ + private static class FixedRates extends Base { + private final double[] rates; + FixedRates(double... rates) { this.rates = rates; } + @Override public void initAndValidate() {} + @Override public double getRateForBranch(Node node) { return rates[node.getNr()]; } + } + + @Test + public void strictClockLeavesVarianceAtZeroAndStillValidates() throws Exception { + // Every branch shares a rate, so variance and coefficient of variation are + // exactly 0. PositiveReal excludes 0, so declaring it here made isValid() + // report false for a perfectly ordinary strict-clock analysis. + RateStatistic rs = statisticOn(new StrictClockModel()); + + double[] values = rs.calcValues(); + assertEquals(0.0, values[1], 0.0, "strict clock should give zero variance"); + assertEquals(0.0, values[2], 0.0, "strict clock should give zero coefficient of variation"); + + assertTrue(rs.isValid(0.0), "the declared domain must admit zero"); + assertTrue(rs.isValid(), "a strict-clock RateStatistic must validate"); + } + + @Test + public void domainStillRejectsNegatives() throws Exception { + // Widening to NonNegativeReal must not weaken the domain to plain Real: + // none of the three statistics can be negative. + RateStatistic rs = statisticOn(new StrictClockModel()); + assertFalse(rs.isValid(-1.0), "negative values are not in the domain"); + } + + @Test + public void varyingRatesProduceTheDocumentedStatistics() throws Exception { + // 3 tips: A, B at height 0, C at height 0; P at 1 joins A and B; root at 2. + // Node numbering: A=0, B=1, C=2, P=3, root=4. The root has no branch, so + // four branches are measured: A, B, C (external) and P (internal). + // branch lengths: A=1, B=1, C=2, P=1 + // rates: A=1, B=2, C=3, P=4 (root's entry is never read) + RateStatistic rs = statisticOn(new FixedRates(1.0, 2.0, 3.0, 4.0, 99.0)); + + double[] values = rs.calcValues(); + // mean is branch-length weighted: (1*1 + 2*1 + 3*2 + 4*1) / (1+1+2+1) + assertEquals((1.0 + 2.0 + 6.0 + 4.0) / 5.0, values[0], 1e-12, + "mean should be weighted by branch length"); + // variance and cv are unweighted over the same four rates + double unweightedMean = (1.0 + 2.0 + 3.0 + 4.0) / 4.0; + double sumSq = 0.0; + for (double r : new double[] { 1.0, 2.0, 3.0, 4.0 }) { + sumSq += (r - unweightedMean) * (r - unweightedMean); + } + double expectedVariance = sumSq / 3.0; // DiscreteStatistics uses count - 1 + assertEquals(expectedVariance, values[1], 1e-12); + assertEquals(Math.sqrt(expectedVariance) / unweightedMean, values[2], 1e-12); + + assertTrue(rs.isValid(), "varying rates should validate too"); + } + + @Test + public void getRejectsOutOfRangeIndices() throws Exception { + // The legacy getArrayValue(dim) guarded with `dim > 3`, so dim == 3 fell + // through to an ArrayIndexOutOfBoundsException and negatives went unchecked. + RateStatistic rs = statisticOn(new StrictClockModel()); + assertEquals(3, rs.size()); + assertThrows(IllegalArgumentException.class, () -> rs.get(3)); + assertThrows(IllegalArgumentException.class, () -> rs.get(-1)); + } + + /** Three-tip ultrametric tree wired to the given branch rate model. */ + private RateStatistic statisticOn(Base branchRateModel) throws Exception { + Node a = leaf("A", 0); + Node b = leaf("B", 1); + Node c = leaf("C", 2); + Node p = internal(3, 1.0, a, b); + Tree tree = new Tree(internal(4, 2.0, p, c)); + + branchRateModel.initAndValidate(); + RateStatistic rs = new RateStatistic(); + rs.initByName("tree", tree, "branchratemodel", branchRateModel); + return rs; + } + + private Node leaf(String id, int nr) { + Node n = new Node(id); + n.setNr(nr); + n.setHeight(0.0); + return n; + } + + private Node internal(int nr, double height, Node left, Node right) { + Node n = new Node(); + n.setNr(nr); + n.setHeight(height); + n.addChild(left); + n.addChild(right); + return n; + } +} diff --git a/beast-base/src/test/resources/beast.base/examples/beast2vs1/testUCRelaxedClockLogNormal.xml b/beast-base/src/test/resources/beast.base/examples/beast2vs1/testUCRelaxedClockLogNormal.xml index 0e88114..9f8e477 100644 --- a/beast-base/src/test/resources/beast.base/examples/beast2vs1/testUCRelaxedClockLogNormal.xml +++ b/beast-base/src/test/resources/beast.base/examples/beast2vs1/testUCRelaxedClockLogNormal.xml @@ -142,7 +142,7 @@ - +