diff --git a/NEWS.md b/NEWS.md index fd41a64ce..249c3ed92 100644 --- a/NEWS.md +++ b/NEWS.md @@ -597,6 +597,14 @@ met during search, and `MultiRatchet()` no longer errors when a starting tree already meets `stopAtScore`. +- The cache behind `ClusteringConcordance(normalize = TRUE)` keyed partitions + on block sizes narrowed to 16 bits, so two partitions whose block sizes + differed by a multiple of 65536 shared an entry and the second was given the + first one's expected mutual information. Reaching this needed a tree of at + least 65536 tips, so no published result is affected; keys now span the full + range of an integer, which rules the collision out rather than making it + unlikely. + # TreeSearch 2.0.0 ## Breaking changes diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index 94e768290..c759a5547 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -158,34 +158,45 @@ std::string mi_key(IntegerVector ni, IntegerVector nj) { Rcpp::stop("ni must be a vector of length 2."); } - std::vector ni_vals = {static_cast(ni[0]), - static_cast(ni[1])}; + // 32 bits spans the whole of `int`, so distinct block sizes always give + // distinct keys. A narrower code aliases: encoded in 16 bits, block sizes + // differing by a multiple of 65536 shared a key, and the cache then served + // one partition's expected mutual information for the other's. + std::vector ni_vals = {static_cast(ni[0]), + static_cast(ni[1])}; std::sort(ni_vals.begin(), ni_vals.end()); - - std::vector nj_vals; + + std::vector nj_vals; nj_vals.reserve(nj.size()); for (int val : nj) { - nj_vals.push_back(static_cast(val)); + nj_vals.push_back(static_cast(val)); } std::sort(nj_vals.begin(), nj_vals.end()); - - // Encode each uint16_t as 4 hex characters — no R allocation needed + + // Encode each value as 8 hex characters — no R allocation needed. Sizing + // the string up front and writing through a pointer beats appending to a + // reserved string, which re-checks capacity on every character. static const char hex[] = "0123456789abcdef"; - std::string key; - key.reserve((2 + nj_vals.size()) * 4); - - for (uint16_t v : ni_vals) { - key += hex[(v >> 12) & 0xF]; - key += hex[(v >> 8) & 0xF]; - key += hex[(v >> 4) & 0xF]; - key += hex[(v) & 0xF]; + std::string key((2 + nj_vals.size()) * 8, '0'); + char *out = &key[0]; + + const auto write_hex = [&out](uint32_t v) { + out[0] = hex[(v >> 28) & 0xF]; + out[1] = hex[(v >> 24) & 0xF]; + out[2] = hex[(v >> 20) & 0xF]; + out[3] = hex[(v >> 16) & 0xF]; + out[4] = hex[(v >> 12) & 0xF]; + out[5] = hex[(v >> 8) & 0xF]; + out[6] = hex[(v >> 4) & 0xF]; + out[7] = hex[(v) & 0xF]; + out += 8; + }; + for (uint32_t v : ni_vals) { + write_hex(v); } - for (uint16_t v : nj_vals) { - key += hex[(v >> 12) & 0xF]; - key += hex[(v >> 8) & 0xF]; - key += hex[(v >> 4) & 0xF]; - key += hex[(v) & 0xF]; + for (uint32_t v : nj_vals) { + write_hex(v); } - + return key; } diff --git a/tests/testthat/test-expected-mi.R b/tests/testthat/test-expected-mi.R index 31e1ac409..f7f995a1c 100644 --- a/tests/testthat/test-expected-mi.R +++ b/tests/testthat/test-expected-mi.R @@ -88,6 +88,37 @@ test_that("expected_mi() agrees across the factorial lookup boundary", { tolerance = 1e-8) }) +test_that("mi_key() distinguishes block sizes above 65535", { + # Sorting is only sound because expected_mi() is invariant under both + # canonicalizations the key applies. + expect_equal(expected_mi(c(3L, 61L), c(30L, 31L)), + expected_mi(c(61L, 3L), c(31L, 30L))) + expect_identical(TreeSearch:::mi_key(c(3L, 61L), c(30L, 31L)), + TreeSearch:::mi_key(c(61L, 3L), c(31L, 30L))) + + # Fixed-width hex of the sorted values. 0x12345678 has eight distinct + # nibbles, so an emit that mis-shifts one is caught here; injectivity + # alone would not notice. + expect_identical( + TreeSearch:::mi_key(c(65597L, 3L), c(305419896L, 30L)), + paste(sprintf("%08x", c(3L, 65597L, 30L, 305419896L)), collapse = "") + ) + + # Block sizes differing by a multiple of 65536 must not share a key + aliases <- c(60L, 61L, 65596L, 65597L, 131133L) + keys <- vapply(aliases, function(n) { + TreeSearch:::mi_key(c(3L, n), c(30L, 31L)) + }, character(1)) + expect_equal(anyDuplicated(keys), 0L) + + # The cached value must belong to the partition asked for. Populate the + # small key first, so a colliding key would return it. + expect_equal(TreeSearch:::.ExpectedMI(c(3L, 61L), c(30L, 31L)), + expected_mi(c(3L, 61L), c(30L, 31L))) + expect_equal(TreeSearch:::.ExpectedMI(c(3L, 65597L), c(30L, 31L)), + expected_mi(c(3L, 65597L), c(30L, 31L))) +}) + 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)