getNumberOfTrees: report the support, not the graph

src/main/java/ccd/algorithms/credibleSets/CredibleCCDComputer.java

@@ -370,18 +370,9 @@ private void writeCurrentResult(double remainingProbability) {
         }
     }
 
-    public static double logBigInteger(BigInteger val) {
-        int precision = Math.max((int) (Math.log(val.bitLength()) / Math.log(2)), 20); // Ensure sufficient precision
-        BigDecimal bigDecimalVal = new BigDecimal(val);
-        int scale = bigDecimalVal.scale();
-
-        // Scale value for improved precision
-        BigDecimal scaledValue = bigDecimalVal.movePointLeft(scale);
-
-        // Compute the logarithm using BigDecimal
-        double log2 = Math.log(scaledValue.doubleValue());
-
-        // Adjust the logarithm based on the scale
-        return log2 + scale * Math.log(10);
-    }
+    // NOTE: a logBigInteger implementation lived here and was removed. It was incorrect:
+    // new BigDecimal(aBigInteger).scale() is always 0, so movePointLeft(scale) was a no-op and
+    // the method reduced to Math.log(val.doubleValue()), which is infinite above ~1.8e308 --
+    // i.e. for essentially every tree count this class deals with.
+    // Use AbstractCCD.logBigInteger(BigInteger) instead.
     */

src/main/java/ccd/model/AbstractCCD.java

@@ -755,6 +755,13 @@ public double getEntropy() {
         return -testro;
     }
 
+    /**
+     * {@inheritDoc}
+     *
+     * <p>The default implementation returns the number of topologies represented by the CCD
+     * graph, which is the support for any model that assigns probability only within its graph.
+     * Full-support subclasses must override this; see {@link ITreeDistribution#getNumberOfTrees()}.
+     */
     @Override
     public BigInteger getNumberOfTrees() {
         if (numberOfTopologiesDirty) {
@@ -764,6 +771,48 @@ public BigInteger getNumberOfTrees() {
         return this.rootClade.getNumberOfTopologies();
     }
 
+    /**
+     * The number of rooted binary topologies on {@code n} labelled taxa, {@code (2n-3)!!}.
+     * Returns {@code 1} for {@code n <= 2}. This is the support size of any full-support model
+     * on {@code n} taxa.
+     *
+     * @param n number of taxa
+     * @return {@code (2n-3)!!} as a {@link BigInteger}
+     */
+    public static BigInteger numberOfRootedTopologies(int n) {
+        BigInteger result = BigInteger.ONE;
+        for (int k = 2 * n - 3; k > 1; k -= 2) {
+            result = result.multiply(BigInteger.valueOf(k));
+        }
+        return result;
+    }
+
+    /**
+     * Natural logarithm of a positive {@link BigInteger}, correct for values far outside
+     * {@code double} range.
+     *
+     * <p>{@code Math.log(value.doubleValue())} is not usable here: a {@code double} overflows to
+     * infinity above about 1.8e308, i.e. beyond roughly 1024 bits, and tree counts in this
+     * package routinely exceed that. Instead the value is shifted right so that at most 1000
+     * bits remain, comfortably inside {@code double} range, and the shift is added back as
+     * {@code shift * log 2}. Precision is unaffected, since {@code double} carries only 53
+     * mantissa bits either way.
+     *
+     * @param value a strictly positive value
+     * @return the natural logarithm of {@code value}
+     * @throws IllegalArgumentException if {@code value} is not positive
+     */
+    public static double logBigInteger(BigInteger value) {
+        if (value.signum() <= 0) {
+            throw new IllegalArgumentException("log of non-positive BigInteger: " + value);
+        }
+        int shift = value.bitLength() - 1000;
+        if (shift > 0) {
+            return Math.log(value.shiftRight(shift).doubleValue()) + shift * Math.log(2.0);
+        }
+        return Math.log(value.doubleValue());
+    }
+
     /**
      * Returns the AIC score of this CCD.
      * The number of parameters depends on the specific CCD.

src/main/java/ccd/model/ITreeDistribution.java

@@ -72,7 +72,25 @@ public interface ITreeDistribution {
      */
     public boolean containsTree(Tree tree);
 
-    /** @return the number of trees (topologies) in this distribution */
+    /**
+     * Returns the size of this distribution's <em>support</em>, that is, the number of distinct
+     * tree topologies to which it assigns non-zero probability.
+     *
+     * <p>For a CCD whose support is exactly the set of topologies its graph represents (CCD0,
+     * CCD1, regCCD) this equals the number of topologies of the graph. For a full-support model
+     * (KRegCCD, MRegCCD, {@link UniformEscapeCCD}) it is the number of rooted topologies on the
+     * taxon set, which is strictly larger: those models place probability outside their graph.
+     * Implementations must report the support, not the graph, so that
+     * {@code getNumberOfTrees()} and {@link #containsTree(Tree)} agree.
+     *
+     * <p>The result is a {@link BigInteger} because it overflows {@code long} at a handful of
+     * taxa and {@code double} not long after: the number of rooted topologies needs 840 bits on
+     * 129 taxa and 1093 bits on 160, and a {@code double} overflows past about 1024 bits. Take
+     * logarithms with {@link AbstractCCD#logBigInteger(BigInteger)}, not via
+     * {@code doubleValue()}, which is infinite from roughly 155 taxa upwards.
+     *
+     * @return the number of topologies with non-zero probability under this distribution
+     */
     public BigInteger getNumberOfTrees();
 
     /**

src/main/java/ccd/model/KRegCCD.java

@@ -12,6 +12,7 @@
 import java.util.concurrent.ConcurrentHashMap;
 import java.util.List;
 import java.util.Map;
+import java.math.BigInteger;
 
 /**
  * Remco's "blue-region" CCD regularisation with <em>full support</em> (every tree
@@ -645,6 +646,18 @@ public boolean containsTree(Tree tree) {
         return true;
     }
 
+    /**
+     * {@inheritDoc}
+     *
+     * <p>KRegCCD is full support, so its support is every rooted topology on the taxon set, not just
+     * the topologies of its graph. The inherited graph count would understate this by many orders
+     * of magnitude and would contradict {@link #containsTree(Tree)}, which is always {@code true}.
+     */
+    @Override
+    public BigInteger getNumberOfTrees() {
+        return numberOfRootedTopologies(getNumberOfLeaves());
+    }
+
     /**
      * Number of internal clades of {@code tree} that are not present in this CCD's observed clade
      * set (the novel-clade count = the total {@code m-2} over the tree's blue regions). This is the

src/main/java/ccd/model/MRegCCD.java

@@ -9,6 +9,7 @@
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
+import java.math.BigInteger;
 
 /**
  * MRegCCD -- the one-parameter "per-new-split" regularised CCD. It unifies RegCCD's split-expansion
@@ -182,6 +183,18 @@ public boolean containsTree(Tree tree) {
         return true;
     }
 
+    /**
+     * {@inheritDoc}
+     *
+     * <p>MRegCCD is full support, so its support is every rooted topology on the taxon set, not just
+     * the topologies of its graph. The inherited graph count would understate this by many orders
+     * of magnitude and would contradict {@link #containsTree(Tree)}, which is always {@code true}.
+     */
+    @Override
+    public BigInteger getNumberOfTrees() {
+        return numberOfRootedTopologies(getNumberOfLeaves());
+    }
+
     private double scoreTree(Tree tree, double scoreMu) {
         Map<Node, BitSet> bits = new HashMap<>();
         computeBits(tree.getRoot(), bits);

src/main/java/ccd/model/UniformEscapeCCD.java

@@ -241,25 +241,7 @@ public String toString() {
      * @return {@code (2n-3)!!} as a {@link BigInteger}
      */
     public static BigInteger numberOfRootedTopologies(int n) {
-        BigInteger result = BigInteger.ONE;
-        for (int k = 2 * n - 3; k > 1; k -= 2) {
-            result = result.multiply(BigInteger.valueOf(k));
-        }
-        return result;
+        return AbstractCCD.numberOfRootedTopologies(n);
     }
 
-    /**
-     * Natural logarithm of a {@link BigInteger}, robust to values far beyond {@code double} range
-     * (shifts the value into the mantissa range and corrects with the shift).
-     */
-    static double logBigInteger(BigInteger value) {
-        if (value.signum() <= 0) {
-            throw new IllegalArgumentException("log of non-positive BigInteger: " + value);
-        }
-        int shift = value.bitLength() - 1000; // keep the top ~1000 bits, well within double range
-        if (shift > 0) {
-            return Math.log(value.shiftRight(shift).doubleValue()) + shift * Math.log(2.0);
-        }
-        return Math.log(value.doubleValue());
-    }
 }

src/test/java/ccd/model/NumberOfTreesContractTest.java

@@ -0,0 +1,116 @@
+package ccd.model;
+
+import beast.base.evolution.tree.Tree;
+import beast.base.evolution.tree.TreeParser;
+import org.junit.jupiter.api.Test;
+
+import java.math.BigInteger;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+
+import static org.junit.jupiter.api.Assertions.assertEquals;
+import static org.junit.jupiter.api.Assertions.assertThrows;
+import static org.junit.jupiter.api.Assertions.assertTrue;
+
+/**
+ * {@link ITreeDistribution#getNumberOfTrees()} must report the size of the <em>support</em>, the
+ * number of topologies with non-zero probability, and so must agree with
+ * {@link ITreeDistribution#containsTree(Tree)}.
+ *
+ * <p>Uses the same four-taxon worked example as {@link UniformEscapeCCDTest}: two sampled trees,
+ * a backbone covering only those two, and 15 rooted topologies in total. A full-support model
+ * must therefore report 15, not 2 -- the failure mode this test exists to catch is a full-support
+ * model inheriting the graph-based count from {@link AbstractCCD}.
+ */
+public class NumberOfTreesContractTest {
+
+    private static final List<String> TAXA = Arrays.asList("A", "B", "C", "D");
+
+    private static Tree parse(String newick) {
+        return new TreeParser(TAXA, newick, 1, false);
+    }
+
+    private static List<Tree> sampleTreeList() {
+        List<Tree> trees = new ArrayList<>();
+        trees.add(parse("(((A:1,B:1):1,C:1):1,D:1):0;"));
+        trees.add(parse("(((D:1,C:1):1,B:1):1,A:1):0;"));
+        return trees;
+    }
+
+    @Test
+    public void ccd1ReportsItsGraph() {
+        CCD1 ccd1 = new CCD1(sampleTreeList(), 0.0);
+        assertEquals(BigInteger.TWO, ccd1.getNumberOfTrees(),
+                "CCD1's support is exactly the topologies its graph represents");
+    }
+
+    @Test
+    public void kRegCCDReportsFullSupport() {
+        KRegCCD kreg = new KRegCCD(sampleTreeList(), 0.0, 0.05, 0.4, 2);
+        assertEquals(BigInteger.valueOf(15), kreg.getNumberOfTrees(),
+                "KRegCCD is full support, so all 15 rooted topologies on four taxa");
+    }
+
+    @Test
+    public void mRegCCDReportsFullSupport() {
+        MRegCCD mreg = new MRegCCD(sampleTreeList(), 0.0, 0.05);
+        assertEquals(BigInteger.valueOf(15), mreg.getNumberOfTrees(),
+                "MRegCCD is full support, so all 15 rooted topologies on four taxa");
+    }
+
+    /** The count and containsTree describe the same set, so they must not disagree. */
+    @Test
+    public void fullSupportCountAgreesWithContainsTree() {
+        for (ITreeDistribution d : List.of(
+                new KRegCCD(sampleTreeList(), 0.0, 0.05, 0.4, 2),
+                new MRegCCD(sampleTreeList(), 0.0, 0.05),
+                new UniformEscapeCCD(sampleTreeList(), 0.0, 0.05))) {
+            boolean full = d.getNumberOfTrees()
+                    .equals(AbstractCCD.numberOfRootedTopologies(d.getNumberOfLeaves()));
+            assertTrue(full, d.getClass().getSimpleName()
+                    + " claims every tree via containsTree, so must count every tree");
+            // A topology absent from the two sampled trees, hence outside the backbone graph.
+            assertTrue(d.containsTree(parse("((A:1,B:1):1,(C:1,D:1):1):0;")),
+                    d.getClass().getSimpleName() + " should contain an off-graph topology");
+        }
+    }
+
+    @Test
+    public void numberOfRootedTopologiesMatchesDoubleFactorial() {
+        assertEquals(BigInteger.ONE, AbstractCCD.numberOfRootedTopologies(1));
+        assertEquals(BigInteger.ONE, AbstractCCD.numberOfRootedTopologies(2));
+        assertEquals(BigInteger.valueOf(3), AbstractCCD.numberOfRootedTopologies(3));
+        assertEquals(BigInteger.valueOf(15), AbstractCCD.numberOfRootedTopologies(4));
+        assertEquals(BigInteger.valueOf(105), AbstractCCD.numberOfRootedTopologies(5));
+    }
+
+    /**
+     * The reason {@link AbstractCCD#logBigInteger(BigInteger)} exists: the counts here overflow
+     * {@code double}, so the naive {@code Math.log(v.doubleValue())} returns infinity.
+     */
+    @Test
+    public void logBigIntegerSurvivesValuesBeyondDoubleRange() {
+        // The naive conversion survives 129 taxa (840 bits) but not 160 (1093 bits): a double
+        // overflows past ~1024 bits. Several data sets in routine use are on the wrong side of
+        // that boundary, which is why the helper exists.
+        BigInteger ok = AbstractCCD.numberOfRootedTopologies(129);
+        assertEquals(840, ok.bitLength(), "129 taxa needs 840 bits");
+        assertTrue(Double.isFinite(ok.doubleValue()), "129 taxa still fits a double");
+        assertEquals(582.13, AbstractCCD.logBigInteger(ok), 0.01);
+
+        BigInteger big = AbstractCCD.numberOfRootedTopologies(160);
+        assertEquals(1093, big.bitLength(), "160 taxa needs 1093 bits");
+        assertTrue(Double.isInfinite(big.doubleValue()), "the naive conversion overflows here");
+        assertTrue(Double.isFinite(AbstractCCD.logBigInteger(big)), "the helper does not");
+
+        // Exact on small values, and consistent with the identity log(a*b) = log a + log b.
+        assertEquals(Math.log(15), AbstractCCD.logBigInteger(BigInteger.valueOf(15)), 1e-12);
+        BigInteger a = AbstractCCD.numberOfRootedTopologies(60);
+        assertEquals(AbstractCCD.logBigInteger(a) + AbstractCCD.logBigInteger(big),
+                AbstractCCD.logBigInteger(a.multiply(big)), 1e-6);
+
+        assertThrows(IllegalArgumentException.class,
+                () -> AbstractCCD.logBigInteger(BigInteger.ZERO));
+    }
+}