Anchor the expected-MI recurrence at the mode; guard expected_mi() and quartet_concordance() inputs

NEWS.md

@@ -491,6 +491,32 @@
   both the `qmApp` (T-302) and `qm` (commit e8b318c3) scalar-unwrap paths,
   confirming all deltas are non-negative and match independent computation.
 
+- `ClusteringConcordance(normalize = TRUE)` now chance-corrects large trees,
+  which it previously left uncorrected while still describing the result as
+  corrected.  The expected mutual information that sets the zero point was
+  accumulated by a recurrence over the hypergeometric distribution of cell
+  overlaps, seeded at the smallest overlap the marginals allow.  That
+  probability sinks below the smallest representable double once a character
+  scores about 1080 tips, and the recurrence being multiplicative, every later
+  term then stayed zero: `expected_mi()` returned exactly 0 where the seed
+  vanished for every block of the character, and a silently truncated sum --
+  as little as a quarter of the true value -- where it vanished for some.  The recurrence is now anchored at the
+  mode of the distribution, whose probability is the largest of at most
+  `N + 1` values summing to one and so is always representable.  The threshold
+  is a property of the tips each character scores rather than of the tree, and
+  only marginals close to even reach it: of 600 random partitions, none below
+  1200 items was affected, 13 of 60 at 1500 items, and 30 of 60 at 3000.
+  Values that were already correct are unchanged, to of order 1e-11 relative.
+
+- `expected_mi()` now checks that `ni` gives exactly two block sizes, as its
+  documentation always required.  A shorter vector was read past its end, and
+  the arbitrary values that produced could index the log-factorial lookup
+  table out of bounds and crash the session.
+
+- `QuartetConcordance()`'s counting kernel now rejects a negative character
+  state code rather than indexing its count buffers out of bounds.  State
+  codes generated by the package are always positive, so no result changes.
+
 # TreeSearch 2.0.0
 
 ## Breaking changes

src/expected_mi.cpp

@@ -6,27 +6,39 @@
 #include <Rcpp.h>
 using namespace Rcpp;
 
-#define MAX_FACTORIAL_LOOKUP 8192
-static double log2_factorial_table[MAX_FACTORIAL_LOOKUP + 1];
-static const double LOG2_E = 1.4426950408889634;
+namespace {
 
-__attribute__((constructor))
-void initialize_factorial_cache() {
-  log2_factorial_table[0] = 0.0;
-  for (int i = 1; i <= MAX_FACTORIAL_LOOKUP; i++) {
-    log2_factorial_table[i] = log2_factorial_table[i - 1] + std::log2(i);
-  }
+constexpr int MAX_FACTORIAL_LOOKUP = 8192;
+constexpr double LOG2_E = 1.4426950408889634;
+
+// Block-scope static: C++17 guarantees the initialization runs exactly once
+// even if several threads reach it together.
+const std::vector<double>& log2_factorial_table() {
+  static const std::vector<double> table = []() {
+    std::vector<double> t(MAX_FACTORIAL_LOOKUP + 1);
+    t[0] = 0.0;
+    for (int i = 1; i <= MAX_FACTORIAL_LOOKUP; ++i) {
+      t[i] = t[i - 1] + std::log2(i);
+    }
+    return t;
+  }();
+  return table;
 }
 
 // Fast lookup with bounds checking
 inline double l2factorial(int n) {
+  if (n < 0) {
+    Rcpp::stop("Factorial undefined for negative arguments.");
+  }
   if (n <= MAX_FACTORIAL_LOOKUP) {
-    return log2_factorial_table[n];
+    return log2_factorial_table()[n];
   } else {
     return lgamma(n + 1) * LOG2_E;
   }
 }
 
+} // namespace
+
 //' Expected mutual information between two partitions
 //'
 //' Computes the mutual information expected purely by chance between two
@@ -54,6 +66,9 @@ inline double l2factorial(int n) {
 //' @export
 // [[Rcpp::export]]
 double expected_mi(const IntegerVector &ni, const IntegerVector &nj) {
+  if (ni.size() != 2) {
+    Rcpp::stop("ni must be a vector of length 2.");
+  }
   // ni and nj are vectors listing the number of entitites in each cluster
   // ni = {a, N-a}; nj = counts of character states
   const int a = ni[0];
@@ -77,36 +92,63 @@ double expected_mi(const IntegerVector &ni, const IntegerVector &nj) {
     if (kmin > kmax) continue;
 
     const double log2mj = std::log2(static_cast<double>(mj));
-
-    // compute P(K=kmin)
-    double log2P = (l2factorial(mj) - l2factorial(kmin) - l2factorial(mj - kmin))
-      + (l2factorial(N - mj) - l2factorial(a - kmin) - l2factorial(N - mj - (a - kmin)))
-      - log2_denom;
-      double Pk = std::pow(2.0, log2P);
-
-      for (int k = kmin; k <= kmax; ++k) {
-        if (Pk > 0.0) {
-          // contribution from inside the split
-          if (k > 0) {
-            double mi_in = std::log2(static_cast<double>(k)) + log2N - (log2a + log2mj);
-            emi += (static_cast<double>(k) * invN) * mi_in * Pk;
-          }
-          // contribution from outside the split
-          int kout = mj - k;
-          if (kout > 0) {
-            double mi_out = std::log2(static_cast<double>(kout)) + log2N - (log2Na + log2mj);
-            emi += (static_cast<double>(kout) * invN) * mi_out * Pk;
-          }
-        }
-        // Update P(k) → P(k+1)
-        if (k < kmax) {
-          double numer = static_cast<double>((mj - k) * (a - k));
-          double denom = static_cast<double>((k + 1) * (N - mj - a + k + 1));
-          Pk *= numer / denom;
-        }
+
+    // Mutual information contributed by an overlap of k, per unit probability
+    const auto cell_mi = [&](int k) {
+      double contribution = 0.0;
+      // contribution from inside the split
+      if (k > 0) {
+        double mi_in = std::log2(static_cast<double>(k)) + log2N - (log2a + log2mj);
+        contribution += (static_cast<double>(k) * invN) * mi_in;
       }
+      // contribution from outside the split
+      const int kout = mj - k;
+      if (kout > 0) {
+        double mi_out = std::log2(static_cast<double>(kout)) + log2N - (log2Na + log2mj);
+        contribution += (static_cast<double>(kout) * invN) * mi_out;
+      }
+      return contribution;
+    };
+
+    // Anchor the recurrence at the mode of the hypergeometric.  P(K = kmode)
+    // is the largest of at most N + 1 probabilities summing to one, so it is
+    // always representable; P(K = kmin) is not — at N = 1200 it is around
+    // 2^-1197, and a recurrence seeded with the zero it underflows to stays
+    // zero for every remaining k.
+    const int kmode = std::min(kmax, std::max(kmin, static_cast<int>(
+      (static_cast<double>(mj) + 1.0) * (static_cast<double>(a) + 1.0) /
+        (static_cast<double>(N) + 2.0))));
+
+    const double log2Pmode =
+      (l2factorial(mj) - l2factorial(kmode) - l2factorial(mj - kmode))
+      + (l2factorial(N - mj) - l2factorial(a - kmode)
+           - l2factorial(N - mj - (a - kmode)))
+      - log2_denom;
+    const double Pmode = std::exp2(log2Pmode);
+
+    emi += cell_mi(kmode) * Pmode;
+
+    // Walk down: P(k - 1) = P(k) * k(N - mj - a + k) / ((mj - k + 1)(a - k + 1))
+    double Pk = Pmode;
+    for (int k = kmode; k > kmin; --k) {
+      Pk *= (static_cast<double>(k) * (N - mj - a + k)) /
+        (static_cast<double>(mj - k + 1) * (a - k + 1));
+      // The distribution is unimodal, so once the tail underflows every
+      // remaining term is likewise negligible.
+      if (!(Pk > 0.0)) break;
+      emi += cell_mi(k - 1) * Pk;
+    }
+
+    // Walk up: P(k + 1) = P(k) * (mj - k)(a - k) / ((k + 1)(N - mj - a + k + 1))
+    Pk = Pmode;
+    for (int k = kmode; k < kmax; ++k) {
+      Pk *= (static_cast<double>(mj - k) * (a - k)) /
+        (static_cast<double>(k + 1) * (N - mj - a + k + 1));
+      if (!(Pk > 0.0)) break;
+      emi += cell_mi(k + 1) * Pk;
+    }
   }
-  
+
   return emi;
 }
 

src/quartet_concordance.cpp

@@ -24,7 +24,13 @@ List quartet_concordance(const LogicalMatrix splits, const IntegerMatrix charact
     for (int t = 0; t < n_taxa; ++t) {
       int state = characters(t, c);
       char_col[t] = state;
-      if (!IntegerVector::is_na(state) && state > max_state) max_state = state;
+      if (!IntegerVector::is_na(state)) {
+        // State codes index n0 / n1 directly
+        if (state < 0) {
+          Rcpp::stop("`characters` must contain non-negative state codes.");
+        }
+        if (state > max_state) max_state = state;
+      }
     }
     // Hoist resize outside split loop: only reallocate when a new character
     // has states beyond the current buffer capacity.

tests/testthat/test-expected-mi.R

@@ -0,0 +1,114 @@
+# Tier 1: arithmetic only, no search, whole file under a second.
+#
+# An independent reference for the expected mutual information under the
+# hypergeometric null.  lchoose() works in log space, so unlike the C++
+# recurrence it cannot underflow at the tails of the distribution.
+ReferenceEmi <- function(ni, nj) {
+  a <- ni[[1]]
+  n <- sum(ni)
+  emi <- 0
+  for (mj in nj) {
+    k <- max(0, a + mj - n):min(a, mj)
+    logP <- lchoose(mj, k) + lchoose(n - mj, a - k) - lchoose(n, a)
+    p <- exp(logP)
+    kOut <- mj - k
+    emi <- emi +
+      sum(p * ifelse(k > 0, (k / n) * log2(k * n / (a * mj)), 0)) +
+      sum(p * ifelse(kOut > 0, (kOut / n) * log2(kOut * n / ((n - a) * mj)), 0))
+  }
+  emi
+}
+
+test_that("expected_mi() is correct for large balanced partitions", {
+  # P(K = kmin) is around 2^-1197 at N = 1200; a recurrence seeded there
+  # returns exactly zero for every k.
+  expect_equal(expected_mi(c(550L, 550L), c(550L, 550L)),
+               ReferenceEmi(c(550L, 550L), c(550L, 550L)), tolerance = 1e-9)
+  expect_equal(expected_mi(c(600L, 600L), c(600L, 600L)),
+               ReferenceEmi(c(600L, 600L), c(600L, 600L)), tolerance = 1e-9)
+  expect_equal(expected_mi(c(1000L, 1000L), c(1000L, 1000L)),
+               ReferenceEmi(c(1000L, 1000L), c(1000L, 1000L)), tolerance = 1e-9)
+
+  # `nj` as ClusteringConcordance() supplies it: a tabulate() over states,
+  # which need not be two
+  expect_equal(expected_mi(c(600L, 600L), c(300L, 300L, 300L, 300L)),
+               ReferenceEmi(c(600L, 600L), c(300L, 300L, 300L, 300L)),
+               tolerance = 1e-9)
+  expect_equal(expected_mi(c(437L, 1063L), c(211L, 396L, 893L)),
+               ReferenceEmi(c(437L, 1063L), c(211L, 396L, 893L)),
+               tolerance = 1e-9)
+
+  # The value the caller that motivated the fix actually receives
+  expect_equal(TreeSearch:::.ExpectedMI(c(600L, 600L), c(600L, 600L)),
+               ReferenceEmi(c(600L, 600L), c(600L, 600L)), tolerance = 1e-9)
+
+  # Chance-corrected mutual information is positive and decreases with N
+  balanced <- vapply(c(500L, 1000L, 1100L, 1200L, 2000L), function(n) {
+    expected_mi(c(n %/% 2L, n %/% 2L), c(n %/% 2L, n %/% 2L))
+  }, double(1))
+  expect_true(all(balanced > 0))
+  expect_true(all(diff(balanced) < 0))
+})
+
+test_that("expected_mi() is unchanged for small partitions", {
+  # Values produced before the recurrence was re-anchored at the mode,
+  # in the regime where seeding it at kmin was safe.
+  expect_equal(expected_mi(c(3L, 4L), c(2L, 5L)),
+               0.15383715015513183, tolerance = 1e-10)
+  expect_equal(expected_mi(c(9L, 11L), c(4L, 7L, 9L)),
+               0.084112593221791668, tolerance = 1e-10)
+  expect_equal(expected_mi(c(50L, 50L), c(50L, 50L)),
+               0.007323652940324857, tolerance = 1e-10)
+  expect_equal(expected_mi(c(37L, 163L), c(11L, 60L, 129L)),
+               0.0077921375666311615, tolerance = 1e-10)
+  expect_equal(expected_mi(c(500L, 500L), c(500L, 500L)),
+               0.00072243147032421289, tolerance = 1e-10)
+  expect_equal(expected_mi(c(400L, 800L), c(300L, 400L, 500L)),
+               0.0012047105794341356, tolerance = 1e-10)
+  expect_equal(expected_mi(c(1L, 6L), c(3L, 4L)), 0.15809905413668374,
+               tolerance = 1e-10)
+  expect_equal(expected_mi(c(0L, 7L), c(3L, 4L)), 0)
+  expect_equal(expected_mi(c(7L, 0L), c(3L, 4L)), 0)
+})
+
+test_that("expected_mi() rejects an `ni` that is not a pair", {
+  expect_error(expected_mi(3L, c(2L, 5L)), "length 2")
+  expect_error(expected_mi(integer(0), c(2L, 5L)), "length 2")
+  expect_error(expected_mi(c(1L, 2L, 4L), c(2L, 5L)), "length 2")
+})
+
+test_that("expected_mi() agrees across the factorial lookup boundary", {
+  # N exceeds the 8192-entry log-factorial table, so l2factorial() must
+  # return matching values from the table and from its lgamma() fallback.
+  expect_equal(expected_mi(c(4500L, 4500L), c(4500L, 4500L)),
+               ReferenceEmi(c(4500L, 4500L), c(4500L, 4500L)),
+               tolerance = 1e-8)
+  expect_equal(expected_mi(c(3000L, 7000L), c(4096L, 5904L)),
+               ReferenceEmi(c(3000L, 7000L), c(4096L, 5904L)),
+               tolerance = 1e-8)
+})
+
+test_that("quartet_concordance() rejects negative state codes", {
+  splits <- matrix(c(TRUE, TRUE, FALSE, FALSE), ncol = 1)
+  characters <- matrix(c(1L, 1L, 2L, 2L), ncol = 1)
+  counts <- TreeSearch:::quartet_concordance(splits, characters)
+  expect_equal(counts[["concordant"]], matrix(1))
+  expect_equal(counts[["decisive"]], matrix(1))
+
+  negative <- matrix(c(1L, -1L, 2L, 2L), ncol = 1)
+  expect_error(TreeSearch:::quartet_concordance(splits, negative),
+               "non-negative")
+
+  # NA marks the absence of a state, and is not a negative code: a taxon
+  # scored NA counts as if it were not in the matrix at all.  Six taxa give
+  # a quartet count that a dropped taxon could change, unlike three.
+  sixSplits <- matrix(c(TRUE, TRUE, TRUE, FALSE, FALSE, FALSE), ncol = 1)
+  sixChars <- matrix(c(1L, 1L, 2L, 2L, 1L, 2L), ncol = 1)
+  sixCounts <- TreeSearch:::quartet_concordance(sixSplits, sixChars)
+  expect_equal(sixCounts[["concordant"]], matrix(1))
+  expect_equal(sixCounts[["decisive"]], matrix(5))
+  expect_equal(
+    TreeSearch:::quartet_concordance(rbind(sixSplits, TRUE),
+                                     rbind(sixChars, NA_integer_)),
+    sixCounts)
+})