diff --git a/src/main/java/ccd/algorithms/credibleSets/CredibleCCDComputer.java b/src/main/java/ccd/algorithms/credibleSets/CredibleCCDComputer.java
index 293b07b..32c5f71 100644
--- a/src/main/java/ccd/algorithms/credibleSets/CredibleCCDComputer.java
+++ b/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.
*/
diff --git a/src/main/java/ccd/model/AbstractCCD.java b/src/main/java/ccd/model/AbstractCCD.java
index 73a8bce..fd5876b 100644
--- a/src/main/java/ccd/model/AbstractCCD.java
+++ b/src/main/java/ccd/model/AbstractCCD.java
@@ -755,6 +755,13 @@ public double getEntropy() {
return -testro;
}
+ /**
+ * {@inheritDoc}
+ *
+ * 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.
+ *
+ *
{@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.
diff --git a/src/main/java/ccd/model/ITreeDistribution.java b/src/main/java/ccd/model/ITreeDistribution.java
index 3bdb59d..2cff739 100644
--- a/src/main/java/ccd/model/ITreeDistribution.java
+++ b/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 support, that is, the number of distinct
+ * tree topologies to which it assigns non-zero probability.
+ *
+ *
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.
+ *
+ *
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();
/**
diff --git a/src/main/java/ccd/model/KRegCCD.java b/src/main/java/ccd/model/KRegCCD.java
index 0018769..ad66285 100644
--- a/src/main/java/ccd/model/KRegCCD.java
+++ b/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 full support (every tree
@@ -645,6 +646,18 @@ public boolean containsTree(Tree tree) {
return true;
}
+ /**
+ * {@inheritDoc}
+ *
+ *
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
diff --git a/src/main/java/ccd/model/MRegCCD.java b/src/main/java/ccd/model/MRegCCD.java
index 2f11b7a..7d68761 100644
--- a/src/main/java/ccd/model/MRegCCD.java
+++ b/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}
+ *
+ *
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 bits = new HashMap<>();
computeBits(tree.getRoot(), bits);
diff --git a/src/main/java/ccd/model/UniformEscapeCCD.java b/src/main/java/ccd/model/UniformEscapeCCD.java
index ac1fb31..1570ffa 100644
--- a/src/main/java/ccd/model/UniformEscapeCCD.java
+++ b/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());
- }
}
diff --git a/src/test/java/ccd/model/NumberOfTreesContractTest.java b/src/test/java/ccd/model/NumberOfTreesContractTest.java
new file mode 100644
index 0000000..e3b46b1
--- /dev/null
+++ b/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 support, the
+ * number of topologies with non-zero probability, and so must agree with
+ * {@link ITreeDistribution#containsTree(Tree)}.
+ *
+ * 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 TAXA = Arrays.asList("A", "B", "C", "D");
+
+ private static Tree parse(String newick) {
+ return new TreeParser(TAXA, newick, 1, false);
+ }
+
+ private static List sampleTreeList() {
+ List 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));
+ }
+}