Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
53 changes: 32 additions & 21 deletions src/expected_mi.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -158,34 +158,45 @@ std::string mi_key(IntegerVector ni, IntegerVector nj) {
Rcpp::stop("ni must be a vector of length 2.");
}

std::vector<uint16_t> ni_vals = {static_cast<uint16_t>(ni[0]),
static_cast<uint16_t>(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<uint32_t> ni_vals = {static_cast<uint32_t>(ni[0]),
static_cast<uint32_t>(ni[1])};
std::sort(ni_vals.begin(), ni_vals.end());
std::vector<uint16_t> nj_vals;

std::vector<uint32_t> nj_vals;
nj_vals.reserve(nj.size());
for (int val : nj) {
nj_vals.push_back(static_cast<uint16_t>(val));
nj_vals.push_back(static_cast<uint32_t>(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;
}
31 changes: 31 additions & 0 deletions tests/testthat/test-expected-mi.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading