Rate statistic

beast-base/src/main/java/beast/base/evolution/RateStatistic.java

@@ -26,23 +26,21 @@
 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;
 import beast.base.evolution.tree.Node;
 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 {

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<NonNegativeReal> {
+
+    final public Input<GenericTreeLikelihood> likelihoodInput = new Input<>("treeLikelihood", "TreeLikelihood containing branch rate model that provides rates for a tree");
+    final public Input<Base> branchRateModelInput = new Input<>("branchratemodel", "model that provides rates for a tree", Validate.XOR, likelihoodInput);
+    final public Input<Tree> treeInput = new Input<>("tree", "tree for which the rates apply", Validate.REQUIRED);
+    final public Input<Boolean> internalInput = new Input<>("internal", "consider internal nodes, default true", true);
+    final public Input<Boolean> 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<Double> 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
+    }
+
+}

beast-base/src/main/java/module-info.java

@@ -148,7 +148,7 @@
         beast.base.inference.operator.kernel.Transform.NoTransform,
         beast.base.inference.operator.kernel.Transform.NoTransformMultivariable,
         beast.base.inference.operator.kernel.Transform.PowerTransform,
         beast.base.evolution.sitemodel.SiteModel,
         beast.base.evolution.speciation.BirthDeathGernhard08Model,
         beast.base.evolution.speciation.CalibratedBirthDeathModel,
         beast.base.evolution.speciation.CalibratedYuleInitialTree,
@@ -166,7 +166,7 @@
         beast.base.evolution.substitutionmodel.Blosum62,
         beast.base.evolution.substitutionmodel.CPREV,
         beast.base.evolution.substitutionmodel.Dayhoff,
         beast.base.evolution.substitutionmodel.Frequencies,
         beast.base.evolution.substitutionmodel.GTR,
         beast.base.evolution.substitutionmodel.GeneralSubstitutionModel,
         beast.base.evolution.substitutionmodel.ComplexSubstitutionModel,
@@ -225,7 +225,7 @@
         beast.base.inference.operator.SwapOperator,
         beast.base.inference.operator.UniformOperator,
         beast.base.inference.operator.UpDownOperator,
         beast.base.inference.operator.kernel.BactrianDeltaExchangeOperator,
         beast.base.inference.operator.kernel.BactrianIntervalOperator,
         beast.base.inference.operator.kernel.BactrianRandomWalkOperator,
         beast.base.inference.operator.kernel.BactrianUpDownOperator,
@@ -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,

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.
+ *
+ * <p>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.</p>
+ */
+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;
+    }
+}

beast-base/src/test/resources/beast.base/examples/beast2vs1/testUCRelaxedClockLogNormal.xml

@@ -142,7 +142,7 @@
             <log idref="hky.frequencies" />
             <log idref="M" />
             <log idref="S" />
-            <log id='rate' spec='beast.base.evolution.RateStatistic' tree='@tree' branchratemodel='@branchRates'/>
+            <log id='rate' spec='beast.base.spec.evolution.RateStatistic' tree='@tree' branchratemodel='@branchRates'/>
             <log idref="rateCategories" /> 
         </logger>
         <!-- <logger logEvery="10000" fileName="test.$(seed).trees">

version.xml

@@ -235,6 +235,7 @@
 
 
                 <provider classname="beast.base.spec.FunctionOfTensor"/>
+                <provider classname="beast.base.spec.evolution.RateStatistic"/>
                 <provider classname="beast.base.spec.evolution.Sum"/>
                 <provider classname="beast.base.spec.evolution.IntSum"/>
                 <provider classname="beast.base.spec.evolution.TreeWithMetaDataLogger"/>