diff --git a/NEWS.md b/NEWS.md index c4372e120..7005da83e 100644 --- a/NEWS.md +++ b/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 diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index 3e02b6e00..94e768290 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -6,27 +6,39 @@ #include 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& log2_factorial_table() { + static const std::vector table = []() { + std::vector 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(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(k)) + log2N - (log2a + log2mj); - emi += (static_cast(k) * invN) * mi_in * Pk; - } - // contribution from outside the split - int kout = mj - k; - if (kout > 0) { - double mi_out = std::log2(static_cast(kout)) + log2N - (log2Na + log2mj); - emi += (static_cast(kout) * invN) * mi_out * Pk; - } - } - // Update P(k) → P(k+1) - if (k < kmax) { - double numer = static_cast((mj - k) * (a - k)); - double denom = static_cast((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(k)) + log2N - (log2a + log2mj); + contribution += (static_cast(k) * invN) * mi_in; } + // contribution from outside the split + const int kout = mj - k; + if (kout > 0) { + double mi_out = std::log2(static_cast(kout)) + log2N - (log2Na + log2mj); + contribution += (static_cast(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( + (static_cast(mj) + 1.0) * (static_cast(a) + 1.0) / + (static_cast(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(k) * (N - mj - a + k)) / + (static_cast(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(mj - k) * (a - k)) / + (static_cast(k + 1) * (N - mj - a + k + 1)); + if (!(Pk > 0.0)) break; + emi += cell_mi(k + 1) * Pk; + } } - + return emi; } diff --git a/src/quartet_concordance.cpp b/src/quartet_concordance.cpp index 0f338be6f..107e6ea55 100644 --- a/src/quartet_concordance.cpp +++ b/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. diff --git a/tests/testthat/test-expected-mi.R b/tests/testthat/test-expected-mi.R new file mode 100644 index 000000000..31e1ac409 --- /dev/null +++ b/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) +})