feat(concordance): Nelson-Ladiges "trit" currency and chance correction for QuartetConcordance

R/RcppExports.R

@@ -20,6 +20,14 @@ MaddisonSlatkin_clear_cache <- function() {
     invisible(.Call(`_TreeSearch_MaddisonSlatkin_clear_cache`))
 }
 
+quartet_expect <- function(splits, characters) {
+    .Call(`_TreeSearch_quartet_expect`, splits, characters)
+}
+
+trit_expect <- function(splits, characters) {
+    .Call(`_TreeSearch_trit_expect`, splits, characters)
+}
+
 #' Expected mutual information between two partitions
 #'
 #' Computes the mutual information expected purely by chance between two

dev/benchmarks/frac-quart/chance_baseline.R

@@ -0,0 +1,86 @@
+# Validate the exact hypergeometric chance baseline used by
+# QuartetConcordance(unit = "trit", normalize = TRUE) against a Monte-Carlo
+# tip-shuffle oracle, at the level of a single (split, character, state-pair).
+#
+# Under the fixed-marginal null the character's tokens are reassigned across the
+# t scored leaves at random, holding the pair's state counts (n_i, n_j) and the
+# split's scored side-A size M fixed.  The cell vector (p,q,r,s) is then
+# trivariate-hypergeometric.  We need the EXPECTED per-pair pool contributions
+#   m       = min(w_c, w_k)          (shared-information pooling weight)
+#   m*A/w_k (edge-numerator contribution)
+#   m*A/w_c (char-numerator contribution)
+# because w_k = (mA-1)+(tP-mA-1)+ -- and hence m -- are THEMSELVES random for
+# multistate pairs (mA = p + r varies).  Binary characters are the special case
+# mA == M, tP == t, so w_k is fixed and only A varies (the clean
+# ClusteringConcordance case).  This mirrors `.ExpectedTrit()` in R/Concordance.R.
+
+pos <- function(z) { z[z < 0] <- 0; z }
+
+# Observed per-pair triple for one realised cell vector.
+pairTriple <- function(p, q, r, s) {
+  nI <- p + q; nJ <- r + s
+  mA <- p + r; tP <- nI + nJ
+  A  <- pos(p - 1) * pos(s - 1) + pos(q - 1) * pos(r - 1)
+  wc <- pos(nI - 1) * pos(nJ - 1)
+  wk <- pos(mA - 1) * pos(tP - mA - 1)
+  m  <- min(wc, wk)
+  c(m    = m,
+    mAwk = if (wk > 0) m * A / wk else 0,
+    mAwc = if (wc > 0) m * A / wc else 0)
+}
+
+# Exact expectation over the trivariate hypergeometric (the R kernel's method).
+expectedTriple <- function(nI, nJ, M, t) {
+  nOther <- t - nI - nJ
+  stopifnot(nOther >= 0, M >= 0, M <= t)
+  logC <- function(n, k) if (k < 0 || k > n) -Inf else lchoose(n, k)
+  denom <- lchoose(t, M)
+  acc <- c(m = 0, mAwk = 0, mAwc = 0)
+  for (p in 0:min(nI, M)) {
+    for (r in 0:min(nJ, M - p)) {
+      oA <- M - p - r
+      if (oA < 0 || oA > nOther) next
+      prob <- exp(logC(nI, p) + logC(nJ, r) + logC(nOther, oA) - denom)
+      if (prob <= 0) next
+      acc <- acc + prob * pairTriple(p, nI - p, r, nJ - r)
+    }
+  }
+  acc
+}
+
+# Monte-Carlo oracle: shuffle state labels across t leaves, M fixed on side A.
+mcTriple <- function(nI, nJ, M, t, nShuf = 2e5) {
+  sideA <- c(rep(TRUE, M), rep(FALSE, t - M))
+  acc <- c(m = 0, mAwk = 0, mAwc = 0)
+  for (i in seq_len(nShuf)) {
+    lab <- sample(rep(1:3, c(nI, nJ, t - nI - nJ)))   # 1 = i, 2 = j, 3 = other
+    p <- sum(lab == 1 & sideA); q <- sum(lab == 1 & !sideA)
+    r <- sum(lab == 2 & sideA); s <- sum(lab == 2 & !sideA)
+    acc <- acc + pairTriple(p, q, r, s)
+  }
+  acc / nShuf
+}
+
+set.seed(1)
+cases <- list(
+  c(nI = 4, nJ = 4, M = 4, t = 8),   # binary-like (nOther = 0): w_k fixed
+  c(nI = 3, nJ = 5, M = 4, t = 8),
+  c(nI = 3, nJ = 3, M = 4, t = 9),   # multistate (nOther > 0): w_k random
+  c(nI = 4, nJ = 2, M = 5, t = 10),
+  c(nI = 5, nJ = 4, M = 6, t = 12),
+  c(nI = 2, nJ = 2, M = 3, t = 7),
+  c(nI = 6, nJ = 3, M = 5, t = 15)
+)
+
+worst <- 0
+for (cs in cases) {
+  ex <- expectedTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"])
+  mc <- mcTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"])
+  worst <- max(worst, max(abs(ex - mc)))
+  cat(sprintf("nI=%d nJ=%d M=%2d t=%2d  exact=(%.4f,%.4f,%.4f)  mc=(%.4f,%.4f,%.4f)\n",
+              cs["nI"], cs["nJ"], cs["M"], cs["t"],
+              ex[1], ex[2], ex[3], mc[1], mc[2], mc[3]))
+}
+cat(sprintf("\nWorst exact-vs-MC discrepancy: %.5f\n", worst))
+stopifnot(worst < 0.05)
+cat("PASS: exact hypergeometric baseline matches the MC tip-shuffle oracle.\n")

dev/benchmarks/frac-quart/cpp_expect_parity.R

@@ -0,0 +1,108 @@
+# Parity guard for the C++ chance-baseline kernels `quartet_expect` and
+# `trit_expect` (src/concordance_expect.cpp), which compute the exact expected
+# concordant/decisive quartet counts (raw) and the expected trit pools
+# (E[m], E[m*A/wk], E[m*A/wc]) under the fixed-marginal (hypergeometric) null.
+#
+# The kernels replaced an earlier pure-R implementation; this script keeps a
+# self-contained R reference (the same trivariate-hypergeometric double sum) and
+# checks the C++ output matches it to machine precision across random binary,
+# multistate and missing-data matrices.  Run after any change to the kernels.
+
+suppressMessages(pkgload::load_all(
+  "C:/Users/pjjg18/GitHub/worktrees/TreeSearch/frac-quart-cpp",
+  quiet = TRUE, recompile = TRUE))
+suppressMessages(library(TreeTools))
+
+# Build logiSplits + charInt exactly as QuartetConcordance() does internally.
+build_inputs <- function(tree, dat) {
+  tl <- intersect(TipLabels(tree), names(dat))
+  dat <- dat[tl, drop = FALSE]
+  splits <- as.Splits(tree, dat)
+  logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]),
+                       logical(NTip(dat)))
+  contrast <- attr(dat, "contrast"); lv <- attr(dat, "allLevels")
+  isInapp <- lv == "-"
+  isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1
+  isGrouping <- !isAmbig & !isInapp
+  gCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which)
+  l2i <- rep(NA_integer_, length(lv)); l2i[isGrouping] <- as.integer(gCols)
+  ch <- PhyDatToMatrix(dat)
+  charInt <- array(l2i[match(ch, lv)], dim = dim(ch), dimnames = dimnames(ch))
+  list(logiSplits = logiSplits, charInt = charInt)
+}
+
+# ---- self-contained R reference (trivariate-hypergeometric double sum) ----
+pos1 <- function(z) if (z < 0) 0 else z
+choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2
+
+ref_pair <- function(nI, nJ, M, t) {           # returns c(Econc, Edec, Em, Ewk, Ewc)
+  nOther <- t - nI - nJ; tP <- nI + nJ
+  wc <- pos1(nI - 1) * pos1(nJ - 1)
+  logDen <- lchoose(t, M)
+  acc <- c(0, 0, 0, 0, 0)
+  for (p in 0:min(nI, M)) for (r in 0:min(nJ, M - p)) {
+    oA <- M - p - r
+    if (oA < 0 || oA > nOther) next
+    prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen)
+    if (prob <= 0) next
+    q <- nI - p; s <- nJ - r; mA <- p + r
+    conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r)
+    dec <- conc + p * q * r * s
+    A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1)
+    wk <- pos1(mA - 1) * pos1(tP - mA - 1); m <- min(wc, wk)
+    acc <- acc + prob * c(conc, dec, m, if (wk > 0) m * A / wk else 0,
+                          if (wc > 0) m * A / wc else 0)
+  }
+  acc
+}
+
+ref_expect <- function(charInt, logiSplits) {
+  nSplit <- ncol(logiSplits); nChar <- ncol(charInt)
+  out <- list(eConc = matrix(0, nSplit, nChar), eDec = matrix(0, nSplit, nChar),
+              eM = matrix(0, nSplit, nChar), eWk = matrix(0, nSplit, nChar),
+              eWc = matrix(0, nSplit, nChar))
+  for (ci in seq_len(nChar)) {
+    col <- charInt[, ci]; scored <- !is.na(col)
+    states <- sort(unique(col[scored])); nStates <- length(states)
+    if (nStates < 2L) next
+    tc <- sum(scored); mSideA <- colSums(logiSplits & scored)
+    uM <- unique(mSideA); idx <- match(mSideA, uM)
+    cnt <- tabulate(match(col[scored], states), nStates)
+    for (a in seq_len(nStates - 1L)) for (b in seq(a + 1L, nStates)) {
+      e <- vapply(uM, function(M) ref_pair(cnt[a], cnt[b], M, tc), double(5))
+      out$eConc[, ci] <- out$eConc[, ci] + e[1, idx]
+      out$eDec[, ci]  <- out$eDec[, ci]  + e[2, idx]
+      out$eM[, ci]    <- out$eM[, ci]    + e[3, idx]
+      out$eWk[, ci]   <- out$eWk[, ci]   + e[4, idx]
+      out$eWc[, ci]   <- out$eWc[, ci]   + e[5, idx]
+    }
+  }
+  out
+}
+
+set.seed(7)
+worst <- 0
+for (rep in 1:10) {
+  nTip <- sample(9:44, 1)
+  tree <- RandomTree(nTip, root = FALSE)
+  m <- replicate(sample(5:30, 1), {
+    ns <- sample(2:4, 1)
+    x <- as.character(sample(0:(ns - 1), nTip, replace = TRUE))
+    if (runif(1) < 0.3) x[sample(nTip, sample(1:3, 1))] <- "?"
+    x
+  })
+  rownames(m) <- paste0("t", seq_len(nTip))
+  io <- build_inputs(tree, MatrixToPhyDat(m))
+  ls <- io$logiSplits; ci <- io$charInt
+  R <- ref_expect(ci, ls)
+  cq <- quartet_expect(ls, ci); ct <- trit_expect(ls, ci)
+  d <- max(abs(cq$concordant - R$eConc), abs(cq$decisive - R$eDec),
+           abs(ct$numEdge - R$eWk), abs(ct$numChar - R$eWc),
+           abs(ct$denM - R$eM))
+  worst <- max(worst, d)
+  cat(sprintf("rep %2d: nTip=%2d nChar=%2d  max|C++ - R| = %.2e\n",
+              rep, nTip, ncol(ci), d))
+}
+cat(sprintf("\nWorst |C++ - R| over all kernels: %.3e\n", worst))
+stopifnot(worst < 1e-9)
+cat("PASS: C++ expectation kernels match the R reference to machine precision.\n")

dev/benchmarks/frac-quart/multistate_trit.R

@@ -0,0 +1,88 @@
+# Verify: for a multistate character, GF(2) rank of its decisive quartets
+# equals Sum_{i<j}(n_i-1)(n_j-1), and rank of concordant quartets (vs a split)
+# equals Sum_{i<j} A_ij.  Blocks are edge-disjoint => trits add across pairs.
+
+gf2rank <- function(M) {                      # rank over GF(2)
+  if (nrow(M) == 0 || ncol(M) == 0) return(0L)
+  M <- M %% 2L; r <- 0L; nc <- ncol(M)
+  for (col in seq_len(nc)) {
+    piv <- which(M[, col] == 1L)
+    piv <- piv[piv > r]
+    if (!length(piv)) next
+    p <- piv[1]; r <- r + 1L
+    M[c(r, p), ] <- M[c(p, r), ]
+    others <- which(M[, col] == 1L); others <- others[others != r]
+    if (length(others)) M[others, ] <- sweep(M[others, , drop = FALSE], 2, M[r, ], `+`) %% 2L
+    if (r == nrow(M)) break
+  }
+  r
+}
+
+# edge id for unordered pair (a,b) among t taxa
+eid <- function(a, b, t) { if (a > b) { tmp <- a; a <- b; b <- tmp }; (a - 1) * t + b }
+
+# All decisive quartets of `char` (state vector), each as its char-induced
+# 4-cycle edge vector over C(t,2)... encoded on t*t (sparse-safe) columns.
+quartet_rows <- function(char, split = NULL, concordant_only = FALSE) {
+  t <- length(char)
+  qs <- combn(t, 4)
+  rows <- list()
+  for (k in seq_len(ncol(qs))) {
+    q <- qs[, k]; st <- char[q]
+    tb <- table(st)
+    if (length(tb) == 2 && all(tb == 2)) {          # "xx yy": decisive
+      states <- as.integer(names(tb))
+      g1 <- q[st == states[1]]; g2 <- q[st == states[2]]   # the two same-state pairs
+      if (concordant_only) {
+        s1 <- split[g1]; s2 <- split[g2]
+        # concordant iff each same-state pair lies wholly on one split side,
+        # and the two pairs are on opposite sides
+        ok <- length(unique(s1)) == 1 && length(unique(s2)) == 1 && s1[1] != s2[1]
+        if (!ok) next
+      }
+      v <- integer(t * t)
+      for (a in g1) for (b in g2) v[eid(a, b, t)] <- 1L    # cross edges = 4-cycle
+      rows[[length(rows) + 1]] <- v
+    }
+  }
+  if (!length(rows)) return(matrix(0L, 0, t * t))
+  do.call(rbind, rows)
+}
+
+Aij <- function(p, q, r, s) max(p - 1, 0) * max(s - 1, 0) + max(q - 1, 0) * max(r - 1, 0)
+
+set.seed(1)
+cat(sprintf("%-28s %6s %6s %6s\n", "case", "rankD", "SumW", "match"))
+for (trial in 1:6) {
+  t <- sample(6:9, 1)
+  r <- sample(2:4, 1)
+  char <- sample(seq_len(r), t, replace = TRUE)
+  while (length(unique(char)) < 2 || any(tabulate(char) == 1 & tabulate(char) > 0 & FALSE)) char <- sample(seq_len(r), t, replace = TRUE)
+  ns <- as.integer(tabulate(char)); ns <- ns[ns > 0]
+  sumW <- sum(outer(seq_along(ns), seq_along(ns), Vectorize(function(i, j)
+    if (i < j) (ns[i] - 1) * (ns[j] - 1) else 0)))
+  rD <- gf2rank(quartet_rows(char))
+  cat(sprintf("t=%d states=%-14s %6d %6d %6s\n",
+              t, paste(ns, collapse = ","), rD, sumW, rD == sumW))
+}
+
+cat("\n-- concordant vs split --\n")
+cat(sprintf("%-34s %6s %6s %6s\n", "case", "rankC", "SumAij", "match"))
+for (trial in 1:6) {
+  t <- sample(6:9, 1); r <- sample(2:4, 1)
+  char <- sample(seq_len(r), t, replace = TRUE)
+  while (length(unique(char)) < 2) char <- sample(seq_len(r), t, replace = TRUE)
+  split <- sample(c(TRUE, FALSE), t, replace = TRUE)
+  states <- sort(unique(char))
+  sumA <- 0
+  for (ii in seq_along(states)) for (jj in seq_along(states)) if (ii < jj) {
+    i <- states[ii]; j <- states[jj]
+    p <- sum(char == i & split); qc <- sum(char == i & !split)
+    rr <- sum(char == j & split); ss <- sum(char == j & !split)
+    sumA <- sumA + Aij(p, qc, rr, ss)
+  }
+  rC <- gf2rank(quartet_rows(char, split, concordant_only = TRUE))
+  cat(sprintf("t=%d states=%-12s split=%-2d/%-2d %6d %6d %6s\n",
+              t, paste(as.integer(tabulate(char))[tabulate(char) > 0], collapse = ","),
+              sum(split), sum(!split), rC, sumA, rC == sumA))
+}