From b0a04614d5040cd94a2ddef64e9af019abfacc1f Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:24:50 +0100 Subject: [PATCH 1/2] Widen expected-MI cache keys to the full integer range mi_key() narrowed each block size to uint16_t, so partitions whose block sizes differed by a multiple of 65536 shared a .ExpectedMICache entry and the second was served the first one's expected mutual information. uint32_t spans the whole of int, so the encoding is now injective and the collision is impossible rather than unlikely. The two sort invariances the key relies on hold: expected_mi() agrees to 2.8e-13 when the blocks of ni are swapped and to 1.5e-15 when nj is permuted, over 2000 random partitions. Fixes #132 Co-Authored-By: Claude Opus 5 --- NEWS.md | 8 ++++++ src/expected_mi.cpp | 43 +++++++++++++++++-------------- tests/testthat/test-expected-mi.R | 23 +++++++++++++++++ 3 files changed, 54 insertions(+), 20 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7005da83e..df9e27296 100644 --- a/NEWS.md +++ b/NEWS.md @@ -517,6 +517,14 @@ state code rather than indexing its count buffers out of bounds. State codes generated by the package are always positive, so no result changes. +- 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..87315f092 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -158,34 +158,37 @@ 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 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]; + key.reserve((2 + nj_vals.size()) * 8); + + const auto append_hex = [&](uint32_t v) { + for (int shift = 28; shift >= 0; shift -= 4) { + key += hex[(v >> shift) & 0xF]; + } + }; + for (uint32_t v : ni_vals) { + append_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) { + append_hex(v); } - + return key; } diff --git a/tests/testthat/test-expected-mi.R b/tests/testthat/test-expected-mi.R index 31e1ac409..b8b02b427 100644 --- a/tests/testthat/test-expected-mi.R +++ b/tests/testthat/test-expected-mi.R @@ -88,6 +88,29 @@ 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))) + + # 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) From e0e4f4284943845d50bc3f3476f5f357e9687c4a Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Fri, 7 Aug 2026 14:12:45 +0100 Subject: [PATCH 2/2] Write the wider cache key through a pointer, not by appending Doubling the hex width made mi_key() measurably slower in isolation: at 2e6 calls per arm, appending to a reserved string costs 250 ns against the 16-bit original's 215 ns at three blocks, rising to 318 vs 239 at eight. Sizing the string once and writing through a pointer drops the per-character capacity check, and is 1-4% faster than the 16-bit original at every block count from 2 to 16 despite emitting twice the characters. The emitted key is unchanged; a new assertion pins it to fixed-width hex of the sorted values, using 0x12345678 so a mis-shifted nibble fails. Co-Authored-By: Claude Opus 5 --- src/expected_mi.cpp | 28 ++++++++++++++++++---------- tests/testthat/test-expected-mi.R | 8 ++++++++ 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index 87315f092..c759a5547 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -173,21 +173,29 @@ std::string mi_key(IntegerVector ni, IntegerVector nj) { } std::sort(nj_vals.begin(), nj_vals.end()); - // Encode each value as 8 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()) * 8); - - const auto append_hex = [&](uint32_t v) { - for (int shift = 28; shift >= 0; shift -= 4) { - key += hex[(v >> shift) & 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) { - append_hex(v); + write_hex(v); } for (uint32_t v : nj_vals) { - append_hex(v); + write_hex(v); } return key; diff --git a/tests/testthat/test-expected-mi.R b/tests/testthat/test-expected-mi.R index b8b02b427..f7f995a1c 100644 --- a/tests/testthat/test-expected-mi.R +++ b/tests/testthat/test-expected-mi.R @@ -96,6 +96,14 @@ test_that("mi_key() distinguishes block sizes above 65535", { 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) {