From 871f125838a9e3d1da442a8be9ce36151a76d7e6 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 12:47:51 +0100 Subject: [PATCH 1/5] feat(concordance): add unit="trit" (Nelson-Ladiges) to QuartetConcordance Add `unit = c("quartet", "trit")` to QuartetConcordance(). "quartet" (default) is unchanged; "trit" scores in Nelson & Ladiges (1992) redundancy-corrected currency, where a split of sizes (k, t-k) carries only (k-1)(t-k-1) independent quartets (4-cycles of K_{k,t-k}; the N&L entailment is GF(2) cycle addition). Scoring uses coverage ("option b"): per state-pair, concordant trits A over the reported unit's own content (Wk for edges, Wc for chars), pooled by shared information M = min(Wc, Wk). Only a split *displayed* by the character scores 1; nested/crossing get strict partial credit. Multistate is handled in the same currency, not deferred: a quartet is decisive only within a state-pair, those K_{n_i,n_j} blocks are edge-disjoint, and the entailment never crosses them, so trits add over pairs (W = sum (n_i-1)(n_j-1), weight 4/(n_i n_j)). Verified by GF(2) rank of the decisive/concordant sets on 2/3/4-state characters. Missing / ambiguous / inapplicable tokens drop per character by construction. Validation (dev/benchmarks/frac-quart/): cell extraction cross-checked against the C++ kernel; A <= min(Wc, Wk) everywhere; spec ladder reproduced; package == an independent re-implementation on binary, multistate and missing-data inputs across all weight x return combos. Co-Authored-By: Claude Opus 4.8 --- R/Concordance.R | 147 ++++++++++- dev/benchmarks/frac-quart/multistate_trit.R | 88 +++++++ dev/benchmarks/frac-quart/validate_trit.R | 229 ++++++++++++++++++ dev/benchmarks/frac-quart/verify_conc_disc.R | 73 ++++++ .../frac-quart/verify_direction_dep.R | 47 ++++ dev/benchmarks/frac-quart/verify_trit_ranks.R | 59 +++++ dev/plans/frac-quart-weight-agreement.md | 222 +++++++++++++++++ inst/REFERENCES.bib | 13 + man/SiteConcordance.Rd | 35 ++- tests/testthat/test-Concordance.R | 100 ++++++++ 10 files changed, 1007 insertions(+), 6 deletions(-) create mode 100644 dev/benchmarks/frac-quart/multistate_trit.R create mode 100644 dev/benchmarks/frac-quart/validate_trit.R create mode 100644 dev/benchmarks/frac-quart/verify_conc_disc.R create mode 100644 dev/benchmarks/frac-quart/verify_direction_dep.R create mode 100644 dev/benchmarks/frac-quart/verify_trit_ranks.R create mode 100644 dev/plans/frac-quart-weight-agreement.md diff --git a/R/Concordance.R b/R/Concordance.R index 2ba797245..7c755c9f9 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -765,8 +765,31 @@ MutualClusteringConcordance <- function(tree, dataset) { #' `QuartetConcordance(return = "char")` returns a numeric vector giving the #' concordance index calculated at each site, averaged across all splits. #' +#' With `unit = "trit"` (see below) the same vectors are returned, but scored in +#' redundancy-corrected currency: values are typically lower, and reach 1 only +#' where a split is *displayed* by the character rather than merely concordant +#' with many of its quartets. +#' #' @param weight Logical specifying whether to weight sites according to the #' number of quartets they are decisive for. +#' @param unit Character specifying the currency in which quartets are counted: +#' - `"quartet"` (default): each resolved quartet counts once, so a character +#' is credited with the full combinatorial volume of quartets it resolves; +#' - `"trit"`: quartets are counted as the *independent* information they carry +#' \insertCite{Nelson1992}{TreeSearch}. Because the quartets resolved by a +#' split of sizes \eqn{(k, t - k)} are logically redundant -- +#' \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only +#' \eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}} +#' resolved quartets are independent. Concordance is then measured against +#' this reduced (redundancy-corrected) content, so that only a character +#' whose split is *identical* to the tree split scores full marks; nested +#' (compatible) characters receive genuine partial support, and crossing +#' (incompatible) characters score lower still. A multistate character is +#' scored in the same currency: trits are counted independently within each +#' pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of +#' \eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary +#' characters remain directly comparable. +#' @references \insertAllCited{} #' @importFrom ape keep.tip #' @importFrom cli cli_progress_bar cli_progress_update #' @importFrom utils combn @@ -776,7 +799,8 @@ QuartetConcordance <- function( tree, dataset = NULL, weight = TRUE, - return = "edge" + return = "edge", + unit = c("quartet", "trit") ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -785,6 +809,7 @@ QuartetConcordance <- function( if (!inherits(dataset, "phyDat")) { stop("`dataset` must be a phyDat object.") } + unit <- match.arg(unit) tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { warning("No overlap between tree labels and dataset.") @@ -816,14 +841,19 @@ QuartetConcordance <- function( dimnames = dimnames(characters) ) + options <- c("character", "site", "default") + return <- options[[pmatch(tolower(trimws(return)), options, + nomatch = length(options))]] + + if (unit == "trit") { + # Return: + return(.TritConcordance(logiSplits, charInt, weight, return, splits)) + } + raw_counts <- quartet_concordance(logiSplits, charInt) num <- raw_counts$concordant den <- raw_counts$decisive - options <- c("character", "site", "default") - return <- options[[pmatch(tolower(trimws(return)), options, - nomatch = length(options))]] - if (return == "default") { if (isTRUE(weight)) { @@ -869,6 +899,113 @@ QuartetConcordance <- function( } } +# Nelson-Ladiges fractional ("trit") currency for QuartetConcordance(). +# +# The quartets a bipartition of sizes (k, t - k) resolves are the 4-cycles of +# the complete bipartite graph K_{k, t-k}; the Nelson-Ladiges entailment +# (ab,cd) + (ab,ce) -> (ab,de) is GF(2) cycle addition, so only the cyclomatic +# number (k - 1)(t - k - 1) of them are independent ("trits"). +# +# A multistate character is the disjoint union of its state-pairs: a quartet is +# decisive only when two taxa share one state and two share another, so every +# decisive quartet lives in the K_{n_i, n_j} block between two states. Those +# blocks are edge-disjoint and the entailment never crosses them (it forces the +# shared "c, d, e" taxa into a single state), so trits add over state pairs: +# W_char = sum_{i 0, m * aij / wk, 0) + numChar[, ci] <- numChar[, ci] + ifelse(wc > 0, m * aij / wc, 0) + wcTot[ci] <- wcTot[ci] + wc[1] # wc is constant across splits + } + } + } + + informative <- wcTot > 0 + + if (return == "default") { + # edge: one value per split + if (isTRUE(weight)) { + denom <- rowSums(denM) + ret <- ifelse(denom == 0, NA_real_, rowSums(numEdge) / denom) + } else { + # Mean per-site quality over informative characters; uninformative + # characters carry no trits and are dropped, as in the quartet path. + sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) + ret <- if (any(informative)) { + rowMeans(sEdge[, informative, drop = FALSE], na.rm = TRUE) + } else { + rep(NA_real_, nSplit) + } + ret[is.nan(ret)] <- NA_real_ + } + setNames(ret, names(splits)) + } else { + # char: one value per character + if (isTRUE(weight)) { + denom <- colSums(denM) + ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + } else { + sChar <- ifelse(denM > 0, numChar / denM, NA_real_) + vapply(seq_len(nChar), function(ci) { + if (informative[ci]) { + m <- mean(sChar[, ci], na.rm = TRUE) + if (is.nan(m)) NA_real_ else m + } else { + NA_real_ + } + }, double(1)) + } + } +} + .ExpectedMICache <- new.env(hash = TRUE, parent = emptyenv()) # @param a must be a vector of length <= 2 diff --git a/dev/benchmarks/frac-quart/multistate_trit.R b/dev/benchmarks/frac-quart/multistate_trit.R new file mode 100644 index 000000000..d11e38e1d --- /dev/null +++ b/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 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)) +} diff --git a/dev/benchmarks/frac-quart/validate_trit.R b/dev/benchmarks/frac-quart/validate_trit.R new file mode 100644 index 000000000..b0f9142e9 --- /dev/null +++ b/dev/benchmarks/frac-quart/validate_trit.R @@ -0,0 +1,229 @@ +# Validation of the unit = "trit" path in QuartetConcordance(). +# +# Ordering follows the advisor's priority: +# 1. Cell extraction cross-checked against the verified C++ kernel (conc/dec). +# 2. A_ij <= min(Wc, Wk) asserted on every real (split, char, pair). +# 3. Spec ladder (t = 8) reproduced from the per-pair formula. +# 4. An INDEPENDENT re-implementation of the measure matches the package on +# binary (congreveLamsdell), multistate, and missing-data inputs. +# 5. `return` parsing quirk mirrored between quartet and trit paths. +# 6. congreveLamsdell quartet -> trit edge-score movement + regime census. +# +# Run: Rscript dev/benchmarks/frac-quart/validate_trit.R + +suppressMessages(pkgload::load_all(".", quiet = TRUE)) +library("TreeTools", quietly = TRUE) + +ok <- function(label, cond) { + cat(sprintf("[%s] %s\n", if (isTRUE(cond)) "PASS" else "FAIL", label)) + if (!isTRUE(cond)) stop("FAILED: ", label) +} + +# ---- Rebuild the (logiSplits, charInt) the function feeds to its kernels ---- +build_inputs <- function(tree, dataset) { + tipLabels <- intersect(TipLabels(tree), names(dataset)) + dataset <- dataset[tipLabels, drop = FALSE] + splits <- as.Splits(tree, dataset) + logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]), + logical(NTip(dataset))) + contrast <- attr(dataset, "contrast") + charLevels <- attr(dataset, "allLevels") + isInapp <- charLevels == "-" + isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1 + isGrouping <- !isAmbig & !isInapp + groupingCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which) + levelToInt <- rep(NA_integer_, length(charLevels)) + levelToInt[isGrouping] <- as.integer(groupingCols) + characters <- PhyDatToMatrix(dataset) + charInt <- array(levelToInt[match(characters, charLevels)], + dim = dim(characters), dimnames = dimnames(characters)) + list(logiSplits = logiSplits, charInt = charInt, splits = splits) +} + +ch2 <- function(x) x * (x - 1) / 2 +posm <- function(z) { z[z < 0] <- 0; z } + +# Per (split, char) cells summed over state-pairs, returned as component lists +# so the same crosstab feeds both the kernel cross-check and the trit measure. +pair_cells <- function(logiSplits, charInt) { + nSplit <- ncol(logiSplits); nChar <- ncol(charInt) + out <- vector("list", nChar) + for (ci in seq_len(nChar)) { + col <- charInt[, ci]; scored <- !is.na(col) + states <- sort(unique(col[scored])) + prs <- list() + if (length(states) >= 2) { + for (a in seq_len(length(states) - 1L)) for (b in seq(a + 1L, length(states))) { + inI <- scored & col == states[a]; inJ <- scored & col == states[b] + prs[[length(prs) + 1L]] <- list( + aI = colSums(logiSplits & inI), bI = colSums(!logiSplits & inI), + aJ = colSums(logiSplits & inJ), bJ = colSums(!logiSplits & inJ)) + } + } + out[[ci]] <- prs + } + out +} + +# ---------- 1. Cell extraction vs the verified kernel ---------- +{ + data("congreveLamsdellMatrices", package = "TreeSearch") + dataset <- congreveLamsdellMatrices[[1]] + tree <- TreeSearch::referenceTree + inp <- build_inputs(tree, dataset) + cells <- pair_cells(inp$logiSplits, inp$charInt) + nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt) + reconC <- reconD <- matrix(0, nSplit, nChar) + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + reconC[, ci] <- reconC[, ci] + ch2(p$aI) * ch2(p$bJ) + ch2(p$bI) * ch2(p$aJ) + reconD[, ci] <- reconD[, ci] + p$aI * p$bI * p$aJ * p$bJ + } + reconD <- reconD + reconC + kern <- TreeSearch:::quartet_concordance(inp$logiSplits, inp$charInt) + ok("cell -> concordant matches kernel", all(abs(reconC - kern$concordant) < 1e-9)) + ok("cell -> decisive matches kernel", all(abs(reconD - kern$decisive) < 1e-9)) + + # ---------- 2. A_ij <= min(Wc, Wk) on every pair ---------- + worst <- -Inf + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1) + nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ + wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1) + worst <- max(worst, max(aij - pmin(wc, wk))) + } + ok("A_ij <= min(Wc, Wk) everywhere", worst <= 1e-9) +} + +# ---------- 3. Spec ladder (t = 8) from the per-pair formula ---------- +pairQ <- function(p, q, r, s) { # edge quality A / Wk for one 2x2 table + A <- posm(p - 1) * posm(s - 1) + posm(q - 1) * posm(r - 1) + m <- p + r; t <- p + q + r + s + Wk <- posm(m - 1) * posm(t - m - 1) + A / Wk +} +ladder <- rbind( + identical = c(4, 0, 0, 4), + nested = c(3, 0, 2, 3), + mildCrossing = c(3, 1, 1, 3), + maxCrossing = c(2, 2, 2, 2)) +Qlad <- apply(ladder, 1, function(x) pairQ(x[1], x[2], x[3], x[4])) +expected <- c(identical = 1, nested = 0.5, mildCrossing = 4 / 9, maxCrossing = 2 / 9) +ok("ladder reproduces spec table", all(abs(Qlad - expected) < 1e-9)) +ok("only identical scores 1", sum(abs(Qlad - 1) < 1e-9) == 1L) +ok("crossing < nested", Qlad["mildCrossing"] < Qlad["nested"] && + Qlad["maxCrossing"] < Qlad["mildCrossing"]) + +# ---------- 4. Independent re-implementation vs the package ---------- +# Deliberately different code shape: per-pair matrices, explicit pooling. +ref_trit <- function(tree, dataset, weight = TRUE, return = "edge") { + inp <- build_inputs(tree, dataset) + cells <- pair_cells(inp$logiSplits, inp$charInt) + nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt) + numE <- numC <- den <- matrix(0, nSplit, nChar); wcTot <- numeric(nChar) + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1) + nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ + wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1) + m <- pmin(wc, wk) + den[, ci] <- den[, ci] + m + numE[, ci] <- numE[, ci] + ifelse(wk > 0, m * aij / wk, 0) + numC[, ci] <- numC[, ci] + ifelse(wc > 0, m * aij / wc, 0) + wcTot[ci] <- wcTot[ci] + wc[1] + } + info <- wcTot > 0 + ret <- pmatch(tolower(trimws(return)), c("character", "site", "default"), + nomatch = 3L) + if (ret == 3L) { # edge + if (isTRUE(weight)) { + d <- rowSums(den); v <- ifelse(d == 0, NA_real_, rowSums(numE) / d) + } else { + sE <- ifelse(den > 0, numE / den, NA_real_) + v <- if (any(info)) rowMeans(sE[, info, drop = FALSE], na.rm = TRUE) + else rep(NA_real_, nSplit) + v[is.nan(v)] <- NA_real_ + } + setNames(v, names(inp$splits)) + } else { # char + if (isTRUE(weight)) { + d <- colSums(den); ifelse(d == 0, NA_real_, colSums(numC) / d) + } else { + sC <- ifelse(den > 0, numC / den, NA_real_) + vapply(seq_len(nChar), function(ci) if (info[ci]) { + mm <- mean(sC[, ci], na.rm = TRUE); if (is.nan(mm)) NA_real_ else mm + } else NA_real_, double(1)) + } + } +} + +same <- function(a, b) all((is.na(a) & is.na(b)) | (abs(a - b) < 1e-9), na.rm = FALSE) + +data("congreveLamsdellMatrices", package = "TreeSearch") +binDat <- congreveLamsdellMatrices[[1]] +refTree <- TreeSearch::referenceTree +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(refTree, binDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(refTree, binDat, weight = w, return = r) + ok(sprintf("binary pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# Multistate toy (3- and 4-state characters, complete) +msMat <- matrix(c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states + 0, 0, 0, 1, 1, 2, 3, 3, # 4 states + 0, 0, 1, 1, 2, 2, 0, 1), 8, # 3 states + dimnames = list(paste0("t", 1:8), NULL)) +msDat <- MatrixToPhyDat(msMat) +msTree <- BalancedTree(8) +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(msTree, msDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(msTree, msDat, weight = w, return = r) + ok(sprintf("multi pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# Missing / ambiguous / inapplicable toy +naMat <- matrix(c(0, 0, 1, 1, "?", "?", 0, 1, + 0, 1, "-", 1, 0, "(01)", 1, 0, + 0, 0, 0, 1, 1, 2, "?", 2), 8, + dimnames = list(paste0("t", 1:8), NULL)) +naDat <- MatrixToPhyDat(naMat) +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(msTree, naDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(msTree, naDat, weight = w, return = r) + ok(sprintf("missing pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# ---------- 5. `return` parsing mirrors the quartet path ---------- +e1 <- QuartetConcordance(refTree, binDat, return = "edge", unit = "trit") +e2 <- QuartetConcordance(refTree, binDat, return = "default", unit = "trit") +c1 <- QuartetConcordance(refTree, binDat, return = "char", unit = "trit") +c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "trit") +c3 <- QuartetConcordance(refTree, binDat, return = "site", unit = "trit") +ok("return edge == default", same(e1, e2)) +ok("return char == character == site", same(c1, c2) && same(c1, c3)) +ok("edge length == n splits, char length == n chars", + length(e1) == length(inp <- as.Splits(refTree)) && + length(c1) == sum(attr(binDat, "weight"))) + +# ---------- 6. congreveLamsdell quartet -> trit movement + census ---------- +qEdge <- QuartetConcordance(refTree, binDat, unit = "quartet") +tEdge <- QuartetConcordance(refTree, binDat, unit = "trit") +cat("\nEdge concordance, quartet vs trit (congreveLamsdell[[1]]):\n") +print(round(rbind(quartet = qEdge, trit = tEdge, delta = tEdge - qEdge), 3)) +cat(sprintf("\nmean quartet = %.3f mean trit = %.3f (trit is stricter)\n", + mean(qEdge, na.rm = TRUE), mean(tEdge, na.rm = TRUE))) + +# Regime census per edge: identical / nested / crossing state-pairs +inpc <- build_inputs(refTree, binDat) +cells <- pair_cells(inpc$logiSplits, inpc$charInt) +regime <- matrix(0L, ncol(inpc$logiSplits), 3, + dimnames = list(names(inpc$splits), + c("identical", "nested", "crossing"))) +for (ci in seq_along(cells)) for (p in cells[[ci]]) { + empt <- (p$aI == 0) + (p$bI == 0) + (p$aJ == 0) + (p$bJ == 0) + regime[, "identical"] <- regime[, "identical"] + (empt == 2) + regime[, "nested"] <- regime[, "nested"] + (empt == 1) + regime[, "crossing"] <- regime[, "crossing"] + (empt == 0) +} +cat("\nState-pair regime census per edge:\n") +print(regime) + +cat("\nAll checks passed.\n") diff --git a/dev/benchmarks/frac-quart/verify_conc_disc.R b/dev/benchmarks/frac-quart/verify_conc_disc.R new file mode 100644 index 000000000..c60fd14c9 --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_conc_disc.R @@ -0,0 +1,73 @@ +# Brute-force check of concordant/discordant/independent-quartet formulae +# for two binary bipartitions with 2x2 taxon table {p,q,r,s}. +choose2 <- function(x) x * (x - 1) / 2 + +# quartet resolution induced by a 0/1 labelling of 4 taxa: returns the +# unordered pair-of-pairs if 2-2, else NA +res <- function(lab) { + if (sum(lab) != 2) return(NA_character_) # not 2-2 + paste(sort(which(lab == lab[order(lab)][1]))[1:2], collapse = "") # placeholder +} +# cleaner: canonical partition key for a 2-2 split of positions 1:4 +key22 <- function(side) { + if (length(unique(side)) != 2 || sum(side == side[1]) != 2) return(NA_character_) + a <- sort(which(side == side[1])) + paste(a, collapse = "") # the pair on 'first' side identifies the 2|2 split up to complement +} + +set.seed(1) +ok <- TRUE +for (rep in 1:2000) { + t <- sample(6:11, 1) + charState <- sample(0:1, t, replace = TRUE) # I = state1 + splitSide <- sample(0:1, t, replace = TRUE) # A = side1 + p <- sum(charState == 1 & splitSide == 1) + q <- sum(charState == 1 & splitSide == 0) + r <- sum(charState == 0 & splitSide == 1) + s <- sum(charState == 0 & splitSide == 0) + + conc <- 0; disc <- 0; decisive <- 0 + combs <- combn(t, 4) + for (j in seq_len(ncol(combs))) { + idx <- combs[, j] + ck <- key22(charState[idx]); sk <- key22(splitSide[idx]) + if (!is.na(ck)) decisive <- decisive + 1 + if (!is.na(ck) && !is.na(sk)) { + # same 2|2 split iff the pair-key matches OR is complementary + same <- (ck == sk) || (ck == paste(setdiff(1:4, as.integer(strsplit(ck,"")[[1]])), collapse="")) + if (same) conc <- conc + 1 else disc <- disc + 1 + } + } + fConc <- choose2(p)*choose2(s) + choose2(q)*choose2(r) + fDisc <- p*q*r*s + n <- p + q + fDecisive <- choose2(n) * choose2(t - n) + if (conc != fConc || disc != fDisc || decisive != fDecisive) { + cat("MISMATCH pqrs=", p,q,r,s, + " conc", conc, fConc, " disc", disc, fDisc, + " dec", decisive, fDecisive, "\n"); ok <- FALSE; break + } + # fractional (independent) agreement with floors + pf <- function(x) pmax(x - 1, 0) + indConc <- pf(p)*pf(s) + pf(q)*pf(r) + indDec <- pf(n)*pf(t - n) # char's own independent quartets + if (indDec > 0) { + FQ <- indConc / indDec + if (FQ < -1e-9 || FQ > 1 + 1e-9) { cat("FQ out of range", FQ, p,q,r,s, "\n"); ok <- FALSE; break } + } + # self-agreement: char == split => q=r=0 => FQ==1 + if (q == 0 && r == 0 && indDec > 0 && abs(FQ - 1) > 1e-9) { + cat("self-agree != 1", FQ, p,q,r,s, "\n"); ok <- FALSE; break + } +} +cat(if (ok) "ALL CHECKS PASS\n" else "FAILED\n") + +# ind_decisive - ind_concordant identity (advisor: = pr + qs - 1 without floors) +set.seed(2); id_ok <- TRUE +for (i in 1:5000) { + p<-sample(0:6,1);q<-sample(0:6,1);r<-sample(0:6,1);s<-sample(0:6,1) + n<-p+q; tt<-p+q+r+s + lhs <- (n-1)*(tt-n-1) - ((p-1)*(s-1) + (q-1)*(r-1)) + if (lhs != p*r + q*s - 1) { cat("identity fail\n"); id_ok<-FALSE; break } +} +cat(if (id_ok) "IDENTITY (no-floor) ind_dec - ind_conc == pr+qs-1: CONFIRMED\n" else "IDENTITY FAIL\n") diff --git a/dev/benchmarks/frac-quart/verify_direction_dep.R b/dev/benchmarks/frac-quart/verify_direction_dep.R new file mode 100644 index 000000000..cf3fcad4d --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_direction_dep.R @@ -0,0 +1,47 @@ +source("gfrank.R") # reuse gf2rank (prints table; ignore) +# Represent the SAME decisive quartet set two ways: +# - by the CHARACTER's pairing (4-cycle in K_{n,l}, I x O) -> Delta_char +# - by the SPLIT's pairing (4-cycle in K_{m,t-m}, A x B) -> Delta_split +# Concordant quartets pair the same either way; discordant differ. +delta_both <- function(p, q, r, s) { + n <- p + q; l <- r + s; m <- p + r; b <- q + s # split sizes + # taxon coords: state (I=1..n grouped p then q), side; build 4 cells + # cell membership: give each taxon (charState in {I,O}, side in {A,B}) + taxa <- rbind( + cbind(I=1, A=1, id=seq_len(p)), # p: I&A + cbind(I=1, A=0, id=seq_len(q)), # q: I&B + cbind(I=0, A=1, id=seq_len(r)), # r: O&A + cbind(I=0, A=0, id=seq_len(s))) # s: O&B + N <- nrow(taxa) + # index taxa within char-graph: I-vertices and O-vertices + Iidx <- which(taxa[,"I"]==1); Oidx <- which(taxa[,"I"]==0) + Aidx <- which(taxa[,"A"]==1); Bidx <- which(taxa[,"A"]==0) + vc <- match(seq_len(N), Iidx); vo <- match(seq_len(N), Oidx) # char graph coord + va <- match(seq_len(N), Aidx); vb <- match(seq_len(N), Bidx) # split graph coord + edgeC <- function(iA, iB) (iA-1)*length(Oidx) + iB # char edge id + edgeS <- function(a, b) (a-1)*length(Bidx) + b # split edge id + rowsC <- list(); rowsS <- list() + quads <- combn(N, 4) + for (j in seq_len(ncol(quads))) { + t4 <- quads[, j] + isI <- taxa[t4,"I"]==1; isA <- taxa[t4,"A"]==1 + if (sum(isI)!=2 || sum(isA)!=2) next # decisive = both resolve 2-2 + Is <- t4[isI]; Os <- t4[!isI]; As <- t4[isA]; Bs <- t4[!isA] + # char pairing: (Is)|(Os) + ec <- c(edgeC(vc[Is[1]],vo[Os[1]]), edgeC(vc[Is[1]],vo[Os[2]]), + edgeC(vc[Is[2]],vo[Os[1]]), edgeC(vc[Is[2]],vo[Os[2]])) + es <- c(edgeS(va[As[1]],vb[Bs[1]]), edgeS(va[As[1]],vb[Bs[2]]), + edgeS(va[As[2]],vb[Bs[1]]), edgeS(va[As[2]],vb[Bs[2]])) + vC <- integer(length(Iidx)*length(Oidx)); vC[ec] <- 1; rowsC[[length(rowsC)+1]] <- vC + vS <- integer(length(Aidx)*length(Bidx)); vS[es] <- 1; rowsS[[length(rowsS)+1]] <- vS + } + mk <- function(L, w) if (length(L)) do.call(rbind, L) else matrix(0,0,w) + c(Dchar = gf2rank(mk(rowsC, length(Iidx)*length(Oidx))), + Dsplit= gf2rank(mk(rowsS, length(Aidx)*length(Bidx))), + char_wt = (n-1)*(l-1), split_wt = (m-1)*(b-1)) +} +for (cfg in list(c(3,1,2,2), c(2,3,4,1), c(5,2,3,4), c(4,1,2,3))) { + d <- delta_both(cfg[1],cfg[2],cfg[3],cfg[4]) + cat(sprintf("p,q,r,s=%-10s Delta_char=%d (charWt=%d) Delta_split=%d (splitWt=%d)\n", + paste(cfg,collapse=","), d["Dchar"], d["char_wt"], d["Dsplit"], d["split_wt"])) +} diff --git a/dev/benchmarks/frac-quart/verify_trit_ranks.R b/dev/benchmarks/frac-quart/verify_trit_ranks.R new file mode 100644 index 000000000..f1830671c --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_trit_ranks.R @@ -0,0 +1,59 @@ +# Independent-quartet counts as GF(2) ranks in the character's bipartite +# cycle space. A character quartet {i1,i2 | o1,o2} = a 4-cycle of K_{n,l} +# (n = I-taxa, l = O-taxa); entailment = cycle addition (Nelson-Ladiges). +gf2rank <- function(M) { # rank over GF(2) + if (!nrow(M) || !ncol(M)) return(0L) + M <- M %% 2; r <- 0L; ncol <- ncol(M); row <- 1L + for (col in seq_len(ncol)) { + piv <- which(M[row:nrow(M), col] == 1) + if (!length(piv)) next + piv <- piv[1] + row - 1L + tmp <- M[row, ]; M[row, ] <- M[piv, ]; M[piv, ] <- tmp + others <- which(M[, col] == 1); others <- others[others != row] + for (o in others) M[o, ] <- (M[o, ] + M[row, ]) %% 2 + row <- row + 1L; r <- r + 1L + if (row > nrow(M)) break + } + r +} + +# taxa: I-vertices 1..n (in A iff <=p), O-vertices 1..l (in A iff <=r) +quartet_vec <- function(i1, i2, o1, o2, l) { + v <- integer(0) # edge id = (i-1)*l + o + e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2) + e +} +counts <- function(p, q, r, s) { + n <- p + q; l <- r + s; nEdge <- n * l + inA_I <- function(i) i <= p + inA_O <- function(o) o <= r + rows_conc <- list(); rows_disc <- list(); rows_dec <- list() + for (i1 in 1:(n-1)) for (i2 in (i1+1):n) for (o1 in 1:(l-1)) for (o2 in (o1+1):l) { + nA <- inA_I(i1) + inA_I(i2) + inA_O(o1) + inA_O(o2) + if (nA != 2) next # split doesn't resolve -> not decisive + e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2) + v <- integer(nEdge); v[e] <- 1 + # concordant iff the two I-taxa are on the same split side + conc <- (inA_I(i1) == inA_I(i2)) + rows_dec[[length(rows_dec)+1]] <- v + if (conc) rows_conc[[length(rows_conc)+1]] <- v + else rows_disc[[length(rows_disc)+1]] <- v + } + mk <- function(L) if (length(L)) do.call(rbind, L) else matrix(0, 0, nEdge) + list(conc = gf2rank(mk(rows_conc)), + disc = gf2rank(mk(rows_disc)), + dec = gf2rank(mk(rows_dec)), + n_conc = length(rows_conc), n_disc = length(rows_disc)) +} +pf <- function(x) pmax(x - 1, 0) +cat(sprintf("%-14s %5s %5s %5s | %-22s %-10s\n", + "p,q,r,s", "conc", "disc", "dec", "A=(p1)(s1)+(q1)(r1)", "n-1 l-1")) +for (cfg in list(c(2,1,1,2), c(3,2,2,3), c(3,1,2,2), c(2,2,2,2), + c(4,1,1,4), c(3,3,1,1), c(2,3,4,1), c(5,2,3,4), + c(1,2,3,4), c(4,0,0,4), c(3,0,2,3), c(2,1,0,3))) { + p<-cfg[1];q<-cfg[2];r<-cfg[3];s<-cfg[4]; n<-p+q; l<-r+s + x <- counts(p,q,r,s) + A <- pf(p)*pf(s) + pf(q)*pf(r) + cat(sprintf("%-14s %5d %5d %5d | A=%-3d disc? dec vs %-4d| char=(%d)(%d)=%d\n", + paste(cfg,collapse=","), x$conc, x$disc, x$dec, A, A + x$disc, n-1, l-1, (n-1)*(l-1))) +} diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md new file mode 100644 index 000000000..a7ff1b578 --- /dev/null +++ b/dev/plans/frac-quart-weight-agreement.md @@ -0,0 +1,222 @@ +# Fractional (Nelson–Ladiges) quartet weighting for `QuartetConcordance()` + +Branch: `frac-quart`. Goal: port Nelson & Ladiges (1992, *Syst. Biol.* 41:490–494) +fractional weighting of three-item statements into the **unrooted quartet** +concordance measure, so that logically redundant quartets are counted as the +information they actually carry (trits) rather than as raw combinatorial volume. + +This file is the settled design. Everything below was derived and, where noted, +verified by brute-force enumeration / GF(2) rank (scratch scripts, this session). + +--- + +## 1. Decisions locked in + +- **Unrooted quartets only.** No polarity / no "informative state" asymmetry. + N&L's rooted `2/n` symmetrises to the quartet weight below. +- **New arg `unit = c("quartet", "trit")`** on `QuartetConcordance()`. It + *replaces* any `fractional` flag. `"quartet"` = current behaviour (raw counts); + `"trit"` = N&L fractional (independent-quartet) currency. +- **Normalization = coverage by the *reported* unit** (option "b"): only a + character-split **identical** to the tree split scores full marks; nesting and + crossing get partial credit. See §4. +- **`weight` (TRUE/FALSE) is retained and orthogonal to `unit`.** It controls the + *pooling amount*, not the per-pair quality. See §5. +- **`return` (edge/char/tree) retained;** it selects which unit is "reported" and + therefore which weight normalises. See §4–5. +- This slots into the package's existing **quality × amount** decomposition + (`ConcordanceTable()` / `QACol()`): quality = per-pair coverage ratio, amount = + N&L fractional weight (trit content). + +--- + +## 2. Currency: trits, quartets = 4-cycles + +For `t` taxa, a binary character (no missing data) is a bipartition of sizes +`(k, t−k)`. The quartets it resolves are the **4-cycles of the complete bipartite +graph `K_{k, t−k}`** (two taxa each side = one 4-cycle). Nelson–Ladiges +entailment `(ab,cd) + (ab,ce) → (ab,de)` is exactly GF(2) cycle addition, so the +independent-quartet count is the cyclomatic number: + +``` +independent quartets of a split of sizes (k, t−k) = (k−1)(t−k−1) # "trits" +raw quartets it resolves = C(k,2)·C(t−k,2) +per-quartet fractional weight = 4 / (k(t−k)) # = indep/raw +``` + +`(k−1)(t−k−1)` is the "fractional weight" `W` of the whole character/split. It is +`0` for a constant (`k=0`) or autapomorphy (`k=1`) — so **uninformative characters +self-zero in trit currency**, no special-casing needed. + +Smallest N&L case (t=5, k=2): `4/(2·3)·3 = 2` trits from 3 raw quartets. ✓ + +**Multistate (solved, not deferred).** A quartet is decisive only when two taxa +share one state and two share another, so every decisive quartet lives in the +`K_{n_i, n_j}` edge-block between two states `i, j`. Those blocks are +edge-disjoint, and the N&L entailment never crosses them (it forces the shared +`c,d,e` taxa into one state), so **trits add over state pairs**: + +``` +W_char = Σ_{i 0` | genuine test; conflict lowers score | + +Nesting is a *passed test*, not homoplasy: no evidence against ⇒ evidence for +(cf. viviparity → Mammalia despite the platypus). It must get partial positive +credit — not exclusion. + +--- + +## 4. The measure (option b): coverage by the reported unit + +Per pair, quality = concordant trits over the **reported unit's** own weight: + +``` +edge (report split k): Q(k,c) = A(k,c) / W_k, W_k = (m−1)(t−m−1) +char (report char c): Q(k,c) = A(k,c) / W_c, W_c = (n−1)(t−n−1) +``` + +`A ≤ W_reported` always (⇒ `Q ∈ [0,1]`), and **per state-pair** +`Q_ij = 1 ⟺ A_ij = W_reported ⟺ identical sub-bipartition`. So for a **binary** +character, only an identical character-split scores full marks; nesting and +crossing are partial, with conflict strictly lowering the score. Verified ladder +(t=8, edge): + +| (p,q,r,s) | relationship | A | W_k | Q | +|---|---|---|---|---| +| 4,0,0,4 | identical | 9 | 9 | 1.00 | +| 3,0,2,3 | nested | 4 | 8 | 0.50 | +| 3,1,1,3 | mild crossing | 4 | 9 | 0.44 | +| 2,2,2,2 | max crossing | 2 | 9 | 0.22 | + +**Multistate: `Q = 1 ⟺ the split is displayed by the character`** (§5 pools the +per-pair `Q_ij` by `M = min(W_c, W_k)`). When every state block lies wholly on +one side of the split, each split-relevant pair scores `Q_ij = 1` and the one +split-orthogonal pair has `W_k^{ij} = 0 ⇒ M = 0` and drops out, so the pooled +`Q = 1` — even though `A_tot < W_c,tot` and the multistate character is *not* +"identical" to the (binary) split. This is the correct generalisation of +"identical" (for a binary character, displayed = identical), and the intended +behaviour: a reproductive character coded `{oviparous | viviparous | ovovivip.}` +that cleanly refines Mammalia fully supports the Mammalia split. A state block +that *straddles* the split makes its pairs crossing (`Q_ij < 1`), lowering the +pooled score. Verified: `t7,t8,t9`-refining char scores `1.0` on every displayed +split, `0.667` on splits that cut a block (see `test-Concordance.R`). + +Raw-currency (`unit="quartet"`) analogue: `Q = conc / (raw quartets of the +reported unit)` — same shape, `conc` and `C(m,2)C(t−m,2)` instead of `A`, `W_k`. + +--- + +## 5. Pooling (`weight` acts on the amount, not the quality) + +`W_k` is constant across characters for a fixed edge, so `weight` must act on the +**amount** (pooling weight), else it goes inert: + +``` +edge FQ(k) = Σ_c M(k,c)·Q(k,c) / Σ_c M(k,c) +char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) +``` + +- `M(k,c)` = shared information = **`min(W_c, W_k)`** (the hBest analogue used by + `ClusteringConcordance`). Preserves "only identical = full marks" (ceiling needs + `A = W_c = W_k`). +- `weight = FALSE` ⇒ `M ≡ 1` (per-pair-equal mean of `Q`). +- `unit = "quartet"` ⇒ `M` = raw shared-quartet count; `Q` uses raw conc. + +**Quality × amount map to the existing table:** `Q` is `QACol`'s *quality* axis; +`W_c` (trit content) is the *amount* axis. The trit port reuses +`ConcordanceTable()` rather than parallelling it. + +--- + +## 6. Build plan + +1. **R prototype first** (binary, no missing). Compute cells `{p,q,r,s}` per + (split, char) from the split logical vectors + character integer vectors; + derive `A`, `W_c`, `W_k`, `Q`, `M`; pool per §5. Keep the existing C++ kernel + for the raw path; the trit path can start in R since cells are a cheap crosstab. +2. **Validate**: (i) brute-force check `A` and the ladder on toy cells; + (ii) run on `congreveLamsdellMatrices` — report edge-score movement + `quartet → trit`, the identical/nested/crossing census per edge, and a + quality/amount table; (iii) confirm `unit="quartet"` reproduces current + `QuartetConcordance` bit-for-bit. +3. **Kernel** (later, perf only): the R per-state-pair path is already correct + and validated for binary + multistate + missing data, so a C++ port of `A` + is now a *performance* optimisation (large matrices), not a correctness + prerequisite. +4. **Tests**: extend `tests/testthat/test-QuartetConcordance*` — `unit` switch, + identical→1, nested partial, crossing < nested, uninformative→dropped, + `unit="quartet"` == legacy. + +## 7. Deferred / open + +- **Multistate**: ~~deferred~~ **DONE.** Trits add over state pairs + (`W_char = Σ_{i= 0 & qt <= 1, na.rm = TRUE)) + + # unit = "quartet" is the default and is byte-identical to omitting `unit`. + expect_identical(QuartetConcordance(tree, identChar), + QuartetConcordance(tree, identChar, unit = "quartet")) + # Coverage normalisation makes trit no laxer than quartet. + expect_true(mean(qt, na.rm = TRUE) <= + mean(QuartetConcordance(tree, identChar, unit = "quartet"), + na.rm = TRUE)) + + # Uninformative characters (constant / autapomorphy) carry no trits -> NA. + autap <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 0, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) + expect_equal(unname(QuartetConcordance(tree, autap, unit = "trit")), + rep(NA_real_, 5)) + + # `unit` is validated. + expect_error(QuartetConcordance(tree, identChar, unit = "trits"), + "should be one of") + + # `return` aliases mirror the quartet path. + expect_equal(QuartetConcordance(tree, identChar, return = "edge", unit = "trit"), + QuartetConcordance(tree, identChar, return = "default", unit = "trit")) + cA <- QuartetConcordance(tree, identChar, return = "char", unit = "trit") + expect_equal(cA, + QuartetConcordance(tree, identChar, return = "character", unit = "trit")) + expect_equal(cA, + QuartetConcordance(tree, identChar, return = "site", unit = "trit")) +}) + +test_that("QuartetConcordance() unit = 'trit' gives nested partial credit", { + # {t1..t5 | t6,t7,t8} split; state {t1,t2,t4} nested within side A but not a + # clade, giving cells (p,q,r,s) = (3,0,2,3) -> A/Wk = 4/8 = 0.5 exactly. + tree <- ape::read.tree(text = "(((((t1,t2),t3),t4),t5),(t6,(t7,t8)));") + char <- MatrixToPhyDat(matrix( + c(0, 0, 1, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) + qt <- QuartetConcordance(tree, char, unit = "trit") + + # The split isolating {t6,t7,t8} scores exactly the nested value 0.5. + sp <- as.Splits(tree) + tips <- TipLabels(sp) + member <- vapply(seq_along(sp), function(i) { + side <- tips[as.logical(sp[[i]])] + setequal(side, c("t6", "t7", "t8")) || + setequal(setdiff(tips, side), c("t6", "t7", "t8")) + }, logical(1)) + expect_equal(unname(qt[member]), 0.5) +}) + +test_that("QuartetConcordance() unit = 'trit' supports multistate", { + # Multistate no longer errors: trits sum over state-pairs (same currency). + tree <- BalancedTree(8) + ms <- MatrixToPhyDat(matrix( + c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states + 0, 0, 0, 1, 1, 2, 3, 3), 8, # 4 states + dimnames = list(paste0("t", 1:8), NULL))) + qt <- QuartetConcordance(tree, ms, unit = "trit") + expect_length(qt, 5L) + expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE)) + expect_false(anyNA(qt)) # both characters are informative on this tree +}) + +test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { + # A multistate character need not be *identical* to a split to score 1: it + # scores full marks for every split its own tree displays (each state block + # wholly on one side; the split-orthogonal state-pair drops via M = 0). This + # generalises "only identical scores 1" from binary to multistate -- a cleanly + # refining reproductive character fully supports the clade it refines. + tree <- ape::read.tree(text = "((((t1,t2),t3),((t4,t5),t6)),((t7,t8),t9));") + char <- MatrixToPhyDat(matrix( + c(0, 0, 0, 1, 1, 1, 2, 2, 2), 9, + dimnames = list(paste0("t", 1:9), NULL))) + qt <- QuartetConcordance(tree, char, unit = "trit") + sp <- as.Splits(tree) + tips <- TipLabels(sp) + atSplit <- function(members) { + i <- which(vapply(seq_along(sp), function(j) { + side <- tips[as.logical(sp[[j]])] + setequal(side, members) || setequal(setdiff(tips, side), members) + }, logical(1))) + unname(qt[i]) + } + # Splits the character displays -> exactly 1 + expect_equal(atSplit(c("t1", "t2", "t3")), 1) + expect_equal(atSplit(c("t4", "t5", "t6")), 1) + expect_equal(atSplit(c("t1", "t2", "t3", "t4", "t5", "t6")), 1) + # A split cutting through a state block is only partially supported + expect_true(atSplit(c("t1", "t2")) < 1) +}) + test_that(".Rezero() works", { expect_equal(TreeSearch:::.Rezero(seq(0, 1, by = 0.1), 0.1), -1:9 / 9) }) From b5226a06fbaf67a2f752ae40d4897fd09f993fa4 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 13:19:18 +0100 Subject: [PATCH 2/5] feat(concordance): add normalize= chance correction to unit="trit" Add a random-expectation baseline to QuartetConcordance(unit = "trit"), mirroring ClusteringConcordance()'s normalize interface: normalize = FALSE (default) no correction; published measure unchanged normalize = TRUE exact expected concordance under a fixed-marginal null normalize = Monte-Carlo tip-shuffle estimate of that expectation Null model: reassign each character's tokens across the leaves at random, holding state counts and split sizes fixed (the multivariate-hypergeometric confusion table; the same null as ClusteringConcordance's miRand and Consistency's ExpectedLength). This is the only coherent null for the trit currency -- it re-zeros each (split, char, state-pair) against its own combinatorial floor and extends unchanged to multistate + missing data. The flat 1/3 SCF baseline is a raw-quartet heuristic that does not survive the trit floors; it is documented as the rejected alternative. Implementation: only A (and, for multistate pairs where a pair's side-A count is random, wk and M = min(wc,wk)) vary under the null, so the exact estimator accumulates E[m], E[m*A/wk], E[m*A/wc] from the hypergeometric pmf (cached on (n_i,n_j,M,t)) and re-zeros the pooled edge/char ratio against the pooled expectation, with the same weight aggregation as the observed score. .Rezero guards the z->1 blow-up (-> NA) and returns values unclamped; .Rezero(1,z) = 1 so a displayed split still scores exactly 1 (ceiling invariant preserved), while conflicting characters go negative. Chance correction for the raw unit="quartet" currency is deferred (it needs the per-state-pair cells the C++ kernel does not expose, and a decision on its null) -- normalize errors there for now. Validated: exact vs MC tip-shuffle agree within MC error at the cell level (dev/benchmarks/frac-quart/chance_baseline.R) and end-to-end (test-Concordance.R). Design settled in dev/plans/frac-quart-weight-agreement.md section 7. Co-Authored-By: Claude Opus 4.8 --- R/Concordance.R | 307 ++++++++++++++++---- dev/benchmarks/frac-quart/chance_baseline.R | 86 ++++++ dev/plans/frac-quart-weight-agreement.md | 51 +++- man/ConcordanceTable.Rd | 30 +- man/SiteConcordance.Rd | 33 ++- tests/testthat/test-Concordance.R | 82 ++++++ 6 files changed, 516 insertions(+), 73 deletions(-) create mode 100644 dev/benchmarks/frac-quart/chance_baseline.R diff --git a/R/Concordance.R b/R/Concordance.R index 7c755c9f9..7db8c95cd 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -87,15 +87,30 @@ NULL #' #' Matching is case‑insensitive and partial. #' -#' @param normalize Controls how the *expected* mutual information (the zero -#' point of the scale) is determined. -#' - `FALSE`: no chance correction; MI is scaled only by its maximum. -#' - `TRUE`: subtract the analytical expected MI for random association. -#' - ``: subtract an empirical expected MI estimated from that -#' number of random trees. -#' -#' In all cases, 1 corresponds to the maximal attainable MI for the pair -#' (`hBest`), and 0 corresponds to the chosen expectation. +#' @param normalize Controls the zero point of the concordance scale by +#' subtracting the value expected under a chance (fixed-marginal) null, in which +#' each character's tokens are reassigned at random across the leaves while its +#' state frequencies and the split sizes are held fixed. +#' - `FALSE`: no chance correction; the measure is scaled only by its maximum, +#' so 1 marks a perfect match and 0 the measure's own (typically positive) +#' floor. +#' - `TRUE`: subtract the expected value under the null, so that 0 marks random +#' expectation and negative values fall below it. `ClusteringConcordance()` +#' uses an analytical approximation to the expected mutual information (fast +#' and generally accurate for large trees, ~200+ taxa, but neglecting +#' correlation between splits); `QuartetConcordance(unit = "trit")` uses the +#' exact hypergeometric expectation. +#' - a positive integer `n`: estimate that expectation by Monte Carlo. +#' `ClusteringConcordance()` fits each character to `n` random trees (and +#' returns Monte-Carlo standard errors, more accurate for small trees where +#' the analytical approximation is biased); `QuartetConcordance(unit = +#' "trit")` averages over `n` random reassignments of each character's +#' tokens. +#' +#' In all cases 1 corresponds to the maximum attainable value. For +#' `QuartetConcordance()`, chance correction is implemented for `unit = "trit"` +#' only, and the returned values are unclamped (they may fall below \eqn{-1}); +#' clamp to \eqn{[-1, 1]} before plotting with [QCol()] / [QACol()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each @@ -800,7 +815,8 @@ QuartetConcordance <- function( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit") + unit = c("quartet", "trit"), + normalize = FALSE ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -810,6 +826,21 @@ QuartetConcordance <- function( stop("`dataset` must be a phyDat object.") } unit <- match.arg(unit) + if (!isFALSE(normalize)) { + if (!isTRUE(normalize)) { + if (!is.numeric(normalize) || length(normalize) != 1L || + is.na(normalize) || normalize < 1) { + stop("`normalize` must be FALSE, TRUE, or a positive integer.") + } + } + if (unit != "trit") { + # The chance baseline needs the per-state-pair cell counts, which the raw + # quartet path computes in the C++ kernel but does not expose; re-baselining + # the raw currency (and choosing its null; see the 1/3 note in the design + # spec) is deferred to the manuscript-figure review gate. + stop("`normalize` chance-correction is implemented for `unit = \"trit\"` only.") + } + } tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { warning("No overlap between tree labels and dataset.") @@ -847,7 +878,8 @@ QuartetConcordance <- function( if (unit == "trit") { # Return: - return(.TritConcordance(logiSplits, charInt, weight, return, splits)) + return(.TritConcordance(logiSplits, charInt, weight, return, splits, + normalize)) } raw_counts <- quartet_concordance(logiSplits, charInt) @@ -924,7 +956,8 @@ QuartetConcordance <- function( # split so both `return`s pool by the same amount. Ambiguous / inapplicable / # absent tokens drop out per character (treated as "?", as in the quartet path), # giving each character its own effective taxon count. -.TritConcordance <- function(logiSplits, charInt, weight, return, splits) { +.TritConcordance <- function(logiSplits, charInt, weight, return, splits, + normalize = FALSE) { nSplit <- ncol(logiSplits) nChar <- ncol(charInt) pos <- function(z) { @@ -932,39 +965,32 @@ QuartetConcordance <- function( z } - # Contributions to the pools, summed over each character's state-pairs. + # Observed contributions to the pools, summed over each character's + # state-pairs; the chance baseline (if requested) uses the same accumulators + # for the expected pools, so the re-zero divides like against like. numEdge <- numChar <- denM <- matrix(0, nSplit, nChar) wcTot <- numeric(nChar) # character trit content (split-free) + doNorm <- !isFALSE(normalize) + if (doNorm) { + baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar) + } for (ci in seq_len(nChar)) { col <- charInt[, ci] - scored <- !is.na(col) - states <- sort(unique(col[scored])) - if (length(states) < 2L) { - next # constant / autapomorphic / empty - } - for (a in seq_len(length(states) - 1L)) { - for (b in seq(a + 1L, length(states))) { - inI <- scored & col == states[a] - inJ <- scored & col == states[b] - aI <- colSums(logiSplits & inI) # state i, side A - bI <- colSums(!logiSplits & inI) # state i, side B - aJ <- colSums(logiSplits & inJ) # state j, side A - bJ <- colSums(!logiSplits & inJ) # state j, side B - nI <- aI + bI # n_i (constant across splits) - nJ <- aJ + bJ # n_j - mA <- aI + aJ # taxa of this pair on side A - tP <- nI + nJ # taxa scored in this pair - # Concordant trits; the (x - 1)_+ floors stop self-agreement exceeding 1. - aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1) - wc <- pos(nI - 1) * pos(nJ - 1) # pair's character content - wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content - m <- pmin(wc, wk) # shared information - denM[, ci] <- denM[, ci] + m - numEdge[, ci] <- numEdge[, ci] + ifelse(wk > 0, m * aij / wk, 0) - numChar[, ci] <- numChar[, ci] + ifelse(wc > 0, m * aij / wc, 0) - wcTot[ci] <- wcTot[ci] + wc[1] # wc is constant across splits + obs <- .CharTritContrib(col, logiSplits, pos) + numEdge[, ci] <- obs[["numEdge"]] + numChar[, ci] <- obs[["numChar"]] + denM[, ci] <- obs[["denM"]] + wcTot[ci] <- obs[["wcTot"]] + if (doNorm && obs[["wcTot"]] > 0) { + base <- if (isTRUE(normalize)) { + .CharTritExpect(col, logiSplits, pos) # exact hypergeometric + } else { + .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo tip-shuffle } + baseNumEdge[, ci] <- base[["numEdge"]] + baseNumChar[, ci] <- base[["numChar"]] + baseDenM[, ci] <- base[["denM"]] } } @@ -975,37 +1001,212 @@ QuartetConcordance <- function( if (isTRUE(weight)) { denom <- rowSums(denM) ret <- ifelse(denom == 0, NA_real_, rowSums(numEdge) / denom) + if (doNorm) { + bDen <- rowSums(baseDenM) + base <- ifelse(bDen == 0, NA_real_, rowSums(baseNumEdge) / bDen) + ret <- .RezeroGuarded(ret, base) + } } else { # Mean per-site quality over informative characters; uninformative # characters carry no trits and are dropped, as in the quartet path. sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) - ret <- if (any(informative)) { - rowMeans(sEdge[, informative, drop = FALSE], na.rm = TRUE) - } else { - rep(NA_real_, nSplit) + ret <- .RowMeanInformative(sEdge, informative, nSplit) + if (doNorm) { + sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) + base <- .RowMeanInformative(sBase, informative, nSplit) + ret <- .RezeroGuarded(ret, base) } - ret[is.nan(ret)] <- NA_real_ } setNames(ret, names(splits)) } else { # char: one value per character if (isTRUE(weight)) { denom <- colSums(denM) - ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + ret <- ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + if (doNorm) { + bDen <- colSums(baseDenM) + base <- ifelse(bDen == 0, NA_real_, colSums(baseNumChar) / bDen) + ret <- .RezeroGuarded(ret, base) + } + ret } else { sChar <- ifelse(denM > 0, numChar / denM, NA_real_) - vapply(seq_len(nChar), function(ci) { - if (informative[ci]) { - m <- mean(sChar[, ci], na.rm = TRUE) - if (is.nan(m)) NA_real_ else m - } else { - NA_real_ - } - }, double(1)) + ret <- .ColMeanInformative(sChar, informative, nChar) + if (doNorm) { + sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_) + base <- .ColMeanInformative(sBase, informative, nChar) + ret <- .RezeroGuarded(ret, base) + } + ret } } } +# Observed per-character trit contributions, summed over the character's +# state-pairs, returned as per-split vectors. Shared by the observed pass and +# the Monte-Carlo baseline so the two use byte-identical arithmetic. +.CharTritContrib <- function(col, logiSplits, pos) { + nSplit <- ncol(logiSplits) + numEdge <- numChar <- denM <- numeric(nSplit) + wcTot <- 0 + scored <- !is.na(col) + states <- sort(unique(col[scored])) + nStates <- length(states) + if (nStates >= 2L) { + for (a in seq_len(nStates - 1L)) { + for (b in seq(a + 1L, nStates)) { + inI <- scored & col == states[a] + inJ <- scored & col == states[b] + aI <- colSums(logiSplits & inI) # state i, side A + bI <- colSums(!logiSplits & inI) # state i, side B + aJ <- colSums(logiSplits & inJ) # state j, side A + bJ <- colSums(!logiSplits & inJ) # state j, side B + nI <- aI + bI # n_i (constant across splits) + nJ <- aJ + bJ # n_j + mA <- aI + aJ # taxa of this pair on side A + tP <- nI + nJ # taxa scored in this pair + # Concordant trits; the (x - 1)_+ floors stop self-agreement exceeding 1. + aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1) + wc <- pos(nI - 1) * pos(nJ - 1) # pair's character content + wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content + m <- pmin(wc, wk) # shared information + denM <- denM + m + numEdge <- numEdge + ifelse(wk > 0, m * aij / wk, 0) + numChar <- numChar + ifelse(wc > 0, m * aij / wc, 0) + wcTot <- wcTot + wc[1] # wc is constant across splits + } + } + } + list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot) +} + +# Exact expected trit contributions under the fixed-marginal (hypergeometric) +# null: each split's scored side-A size `M` and each pair's state counts +# `n_i, n_j` are held fixed while the character's tokens are reassigned at +# random across the `t` scored leaves. Only `A`, and (for multistate pairs) +# `wk` and `m`, vary; so we accumulate E[m], E[m * A / wk] and E[m * A / wc] +# from the exact pmf, matching `.CharTritContrib`'s pool definitions. +.CharTritExpect <- function(col, logiSplits, pos) { + nSplit <- ncol(logiSplits) + numEdge <- numChar <- denM <- numeric(nSplit) + scored <- !is.na(col) + tc <- sum(scored) + mSideA <- colSums(logiSplits & scored) # M per split + uM <- unique(mSideA) + idx <- match(mSideA, uM) + states <- sort(unique(col[scored])) + nStates <- length(states) + if (nStates >= 2L) { + cnt <- tabulate(match(col[scored], states), nStates) # state counts + for (a in seq_len(nStates - 1L)) { + for (b in seq(a + 1L, nStates)) { + et <- vapply(uM, function(M) .ExpectedTrit(cnt[a], cnt[b], M, tc), + double(3)) + denM <- denM + et[1, idx] + numEdge <- numEdge + et[2, idx] + numChar <- numChar + et[3, idx] + } + } + } + list(numEdge = numEdge, numChar = numChar, denM = denM) +} + +.ExpectedTritCache <- new.env(hash = TRUE, parent = emptyenv()) + +# Exact E[m], E[m * A / wk], E[m * A / wc] for one state-pair, over the +# trivariate hypergeometric placement of `nI` state-i and `nJ` state-j leaves +# (and `t - nI - nJ` other leaves) among the `M` leaves on side A of a split. +.ExpectedTrit <- function(nI, nJ, M, t) { + key <- paste(nI, nJ, M, t, sep = ",") + cached <- .ExpectedTritCache[[key]] + if (!is.null(cached)) { + return(cached) + } + pos1 <- function(z) if (z < 0) 0 else z + nOther <- t - nI - nJ + tP <- nI + nJ # taxa in this pair + wc <- pos1(nI - 1) * pos1(nJ - 1) # pair's character content (fixed) + logDen <- lchoose(t, M) + acc <- c(0, 0, 0) # E[m], E[m*A/wk], E[m*A/wc] + for (p in 0:min(nI, M)) { + for (r in 0:min(nJ, M - p)) { + oA <- M - p - r # other-state leaves on side A + 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 # this pair's taxa on side A (random) + A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1) + wk <- pos1(mA - 1) * pos1(tP - mA - 1) # pair's split content (random) + m <- min(wc, wk) + acc[1] <- acc[1] + prob * m + if (wk > 0) acc[2] <- acc[2] + prob * m * A / wk + if (wc > 0) acc[3] <- acc[3] + prob * m * A / wc + } + } + .ExpectedTritCache[[key]] <- acc + acc +} + +# Monte-Carlo baseline: average `.CharTritContrib` over `nRelabel` random +# reassignments of the character's tokens across its scored leaves (the split, +# the scored set and the state counts are held fixed). Averaging the pools +# (rather than per-split ratios) mirrors the exact ratio-of-expectations path. +.CharTritMC <- function(col, logiSplits, pos, nRelabel) { + nSplit <- ncol(logiSplits) + scored <- !is.na(col) + tokens <- col[scored] + accEdge <- accChar <- accDen <- numeric(nSplit) + for (i in seq_len(nRelabel)) { + shuffled <- col + shuffled[scored] <- sample(tokens) + cc <- .CharTritContrib(shuffled, logiSplits, pos) + accEdge <- accEdge + cc[["numEdge"]] + accChar <- accChar + cc[["numChar"]] + accDen <- accDen + cc[["denM"]] + } + list(numEdge = accEdge / nRelabel, + numChar = accChar / nRelabel, + denM = accDen / nRelabel) +} + +# Mean of per-split quality across informative characters, matching the +# quartet path's treatment of uninformative characters (dropped). +.RowMeanInformative <- function(mat, informative, nSplit) { + ret <- if (any(informative)) { + rowMeans(mat[, informative, drop = FALSE], na.rm = TRUE) + } else { + rep(NA_real_, nSplit) + } + ret[is.nan(ret)] <- NA_real_ + ret +} + +.ColMeanInformative <- function(mat, informative, nChar) { + vapply(seq_len(nChar), function(ci) { + if (informative[ci]) { + m <- mean(mat[, ci], na.rm = TRUE) + if (is.nan(m)) NA_real_ else m + } else { + NA_real_ + } + }, double(1)) +} + +# Re-zero to a chance baseline, guarding the degenerate `zero -> 1` case +# (returns NA rather than dividing by ~0). Values are returned UNCLAMPED, as +# in `ClusteringConcordance`; clamping to [-1, 1] is deferred to plotting. +.RezeroGuarded <- function(value, zero) { + denom <- 1 - zero + ifelse(is.na(value) | is.na(zero) | denom < sqrt(.Machine[["double.eps"]]), + NA_real_, (value - zero) / denom) +} + .ExpectedMICache <- new.env(hash = TRUE, parent = emptyenv()) # @param a must be a vector of length <= 2 diff --git a/dev/benchmarks/frac-quart/chance_baseline.R b/dev/benchmarks/frac-quart/chance_baseline.R new file mode 100644 index 000000000..93b9d9690 --- /dev/null +++ b/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") diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md index a7ff1b578..6840a9142 100644 --- a/dev/plans/frac-quart-weight-agreement.md +++ b/dev/plans/frac-quart-weight-agreement.md @@ -212,9 +212,54 @@ char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) toy with `?`/`-`/ambiguous tokens (`validate_trit.R`), but the validation corpus (congreveLamsdell) has no missing data, so the semantics are proven self-consistent (`Q∈[0,1]`, pkg==independent ref) rather than benchmarked. -- **Random-expectation baseline**: a `.Rezero`/`miRand` analogue (crossing floor - ≈ 0.22 in §4 is not zero) — decide whether trit concordance needs chance - correction. Deferred. +- **Random-expectation baseline**: ~~deferred~~ **DONE** (`unit = "trit"`). + New arg `normalize = FALSE / TRUE / ` on `QuartetConcordance()`, mirroring + `ClusteringConcordance()`. The crossing floor (≈ 0.22 in §4) is not zero, so + without correction a maximally conflicting character still scores positive and + the floor is split-size dependent; `normalize` re-zeros against the + concordance *expected by chance*, exactly as `ClusteringConcordance()` re-zeros + MI against `miRand`. + - **Null model (chosen): fixed-marginal randomization** — reassign each + character's tokens across the leaves at random, holding its state counts and + the split sizes fixed (the multivariate-hypergeometric confusion table; the + same null as `ClusteringConcordance`'s `miRand` and `Consistency`'s + `ExpectedLength` tip-shuffle). This is the *only* coherent null for the trit + currency: it re-zeros each `(split, char, state-pair)` against its own + combinatorial floor and extends unchanged to multistate and missing data. + - **Rejected alternative: the flat Minh-style `1/3`** (three quartet + resolutions ⇒ expected concordance `1/3`). It is a *raw-quartet* heuristic + that ignores marginals and does **not** extend to floored trit counts (the + trit floor is `0.22`, not `1/3`), so it cannot re-zero the trit measure + coherently. Its only remaining use would be re-baselining the *raw* `quartet` + unit — see the raw-unit note below. + - **Two estimators** (mirroring `ClusteringConcordance`): `normalize = TRUE` + computes the **exact** expected pools from the hypergeometric pmf (a small + double sum, cached on `(n_i, n_j, M, t)`); `normalize = ` estimates them + by **Monte-Carlo** tip-shuffle. Validated to agree within MC error + (`dev/benchmarks/frac-quart/chance_baseline.R` at the cell level; + `test-Concordance.R` end-to-end). + - **What varies under the null.** Only `A` (and, for *multistate* pairs where a + pair's side-A count `mA` is itself random, `wk` and `M = min(w_c, w_k)`) vary; + for *binary* characters `mA = M`, `tP = t` so `wk` is fixed and only `A` + varies (the clean `ClusteringConcordance` case). The exact estimator therefore + accumulates `E[m]`, `E[m·A/w_k]`, `E[m·A/w_c]` and re-zeros the *pooled* + edge/char ratio against the pooled expectation, with the **same** `weight` + aggregation as the observed score. + - **Guarantees / guards.** `.Rezero(1, z) = 1`, so a **displayed** split still + scores exactly `1` (the ceiling invariant of §4 survives chance correction — + only the floor is lifted); conflicting characters go **negative** ("below + random", as in `ClusteringConcordance`). Values are returned **unclamped** + (may fall below `−1`); clamp to `[−1, 1]` only at plot time (`QCol`/`QACol`), + never inside the measure. The `z → 1` blow-up is guarded (→ `NA`). + `normalize` defaults to `FALSE`, so the published measure is unchanged unless + the user opts in. + - **Raw `quartet` unit: deferred.** `normalize` currently errors for + `unit = "quartet"`. Re-baselining the raw currency needs the per-state-pair + cells (computed inside the C++ kernel but not exposed) and a decision on its + null (fixed-marginal vs the flat `1/3`); both are folded into the C++ port + (task 2) and the manuscript-figure review (task 3), which will compare + trit-with vs trit-without correction and decide whether the raw unit needs it + too. - **`M = min(W_c,W_k)`** (settled): the pooling amount is the *shared* information, so it must be symmetric between character and split (both `return`s pool by the same `M`); `M = W_k` would make the char return weight a diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 9d6a0831b..cc617bca8 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -36,17 +36,31 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.} \item{ylab}{Character giving a label for the y axis.} -\item{normalize}{Controls how the \emph{expected} mutual information (the zero -point of the scale) is determined. +\item{normalize}{Controls the zero point of the concordance scale by +subtracting the value expected under a chance (fixed-marginal) null, in which +each character's tokens are reassigned at random across the leaves while its +state frequencies and the split sizes are held fixed. \itemize{ -\item \code{FALSE}: no chance correction; MI is scaled only by its maximum. -\item \code{TRUE}: subtract the analytical expected MI for random association. -\item \verb{}: subtract an empirical expected MI estimated from that -number of random trees. +\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, +so 1 marks a perfect match and 0 the measure's own (typically positive) +floor. +\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random +expectation and negative values fall below it. \code{ClusteringConcordance()} +uses an analytical approximation to the expected mutual information (fast +and generally accurate for large trees, ~200+ taxa, but neglecting +correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the +exact hypergeometric expectation. +\item a positive integer \code{n}: estimate that expectation by Monte Carlo. +\code{ClusteringConcordance()} fits each character to \code{n} random trees (and +returns Monte-Carlo standard errors, more accurate for small trees where +the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's +tokens. } -In all cases, 1 corresponds to the maximal attainable MI for the pair -(\code{hBest}), and 0 corresponds to the chosen expectation.} +In all cases 1 corresponds to the maximum attainable value. For +\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} +only, and the returned values are unclamped (they may fall below \eqn{-1}); +clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index f61cddce8..3dd1a5604 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -18,7 +18,8 @@ QuartetConcordance( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit") + unit = c("quartet", "trit"), + normalize = FALSE ) PhylogeneticConcordance(tree, dataset) @@ -43,17 +44,31 @@ split–character pair. Matching is case‑insensitive and partial.} -\item{normalize}{Controls how the \emph{expected} mutual information (the zero -point of the scale) is determined. +\item{normalize}{Controls the zero point of the concordance scale by +subtracting the value expected under a chance (fixed-marginal) null, in which +each character's tokens are reassigned at random across the leaves while its +state frequencies and the split sizes are held fixed. \itemize{ -\item \code{FALSE}: no chance correction; MI is scaled only by its maximum. -\item \code{TRUE}: subtract the analytical expected MI for random association. -\item \verb{}: subtract an empirical expected MI estimated from that -number of random trees. +\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, +so 1 marks a perfect match and 0 the measure's own (typically positive) +floor. +\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random +expectation and negative values fall below it. \code{ClusteringConcordance()} +uses an analytical approximation to the expected mutual information (fast +and generally accurate for large trees, ~200+ taxa, but neglecting +correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the +exact hypergeometric expectation. +\item a positive integer \code{n}: estimate that expectation by Monte Carlo. +\code{ClusteringConcordance()} fits each character to \code{n} random trees (and +returns Monte-Carlo standard errors, more accurate for small trees where +the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's +tokens. } -In all cases, 1 corresponds to the maximal attainable MI for the pair -(\code{hBest}), and 0 corresponds to the chosen expectation.} +In all cases 1 corresponds to the maximum attainable value. For +\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} +only, and the returned values are unclamped (they may fall below \eqn{-1}); +clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the number of quartets they are decisive for.} diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index 88e32f625..83c7adcd2 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -233,6 +233,88 @@ test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { expect_true(atSplit(c("t1", "t2")) < 1) }) +test_that("QuartetConcordance() unit = 'trit' normalize contract", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + ident <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, + dimnames = list(paste0("t", 1:9), NULL))) + + # normalize = FALSE is the default and leaves the measure untouched. + expect_identical( + QuartetConcordance(tree, ident, unit = "trit", normalize = FALSE), + QuartetConcordance(tree, ident, unit = "trit")) + + # Invalid normalize is rejected. + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = 0), + "positive integer") + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = -3), + "positive integer") + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"), + "positive integer") + + # Chance correction is implemented for the trit currency only, for now. + expect_error( + QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE), + "trit") + + # An identical (displayed) split still scores exactly 1 after re-zeroing: + # .Rezero(1, z) == 1, so chance correction lifts the floor without moving the + # ceiling. + sp <- as.Splits(tree) + tips <- TipLabels(sp) + col789 <- which(vapply(seq_along(sp), function(j) { + side <- tips[as.logical(sp[[j]])] + setequal(side, c("t7", "t8", "t9")) || + setequal(setdiff(tips, side), c("t7", "t8", "t9")) + }, logical(1))) + corrected <- QuartetConcordance(tree, ident, unit = "trit", normalize = TRUE) + expect_equal(unname(corrected[col789]), 1) +}) + +test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + dat <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} + 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested: {t8,t9} + 0, 0, 1, 1, 0, 0, 1, 1, 0), 9, # crossing + dimnames = list(paste0("t", 1:9), NULL))) + + # The exact hypergeometric baseline and the Monte-Carlo tip-shuffle baseline + # target the same null, so they agree to within MC error (edge and char). + exact <- QuartetConcordance(tree, dat, unit = "trit", normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "trit", normalize = 6000) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.03) + + exactC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", + normalize = TRUE) + set.seed(2) + mcC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", + normalize = 6000) + expect_lt(max(abs(exactC - mcC), na.rm = TRUE), 0.03) +}) + +test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + support <- MatrixToPhyDat(matrix( # agrees with the tree: displays {t7,t8,t9} + c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, + dimnames = list(paste0("t", 1:9), NULL))) + conflict <- MatrixToPhyDat(matrix( # crosses the two halves of the tree + c(0, 0, 1, 1, 0, 0, 1, 1, 0), 9, + dimnames = list(paste0("t", 1:9), NULL))) + + sChar <- QuartetConcordance(tree, support, return = "char", unit = "trit", + normalize = TRUE) + cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "trit") + cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "trit", + normalize = TRUE) + # A character that agrees with the tree scores above random expectation; + # a crossing character is pulled below it, and below its uncorrected score. + expect_gt(sChar, 0) + expect_lt(cChar, cNone) + expect_lt(cChar, 0) +}) + test_that(".Rezero() works", { expect_equal(TreeSearch:::.Rezero(seq(0, 1, by = 0.1), 0.1), -1:9 / 9) }) From 1318e8f4313200b3c0604d10744833de51c0dc18 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 13:24:55 +0100 Subject: [PATCH 3/5] fix(concordance): match cells in weight=FALSE trit re-zero; test multistate Two gaps in the unit="trit" chance correction, both in the multistate path: - weight=FALSE re-zero averaged observed and baseline over DIFFERENT (split, char) cells: observed drops a cell when denM==0 (na.rm), but the baseline kept it whenever baseDenM>0 (a degenerate multistate pair can have observed wk=0 yet E[m]>0). Mask both to NA on the union of NA cells so the edge/char row-means cover the same population before .Rezero. weight=TRUE (ratio-of-sums) is unaffected. - The exact-vs-MC end-to-end test used only binary characters, so the random-wk multistate regime -- the reason E[A] was elaborated into E[m*A/wk] + E[m] -- had no oracle. Add a 3-state column and cover edge/char x weight in {TRUE, FALSE}, asserting exact and MC agree on which cells are NA (guards the cell-matching fix) and within MC error elsewhere. Co-Authored-By: Claude Opus 4.8 --- R/Concordance.R | 25 +++++++++++++++++++------ tests/testthat/test-Concordance.R | 31 ++++++++++++++++++------------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 7db8c95cd..7e3542ba6 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -1010,11 +1010,19 @@ QuartetConcordance <- function( # Mean per-site quality over informative characters; uninformative # characters carry no trits and are dropped, as in the quartet path. sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) - ret <- .RowMeanInformative(sEdge, informative, nSplit) if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) - base <- .RowMeanInformative(sBase, informative, nSplit) - ret <- .RezeroGuarded(ret, base) + # Average observed and baseline over the SAME (split, char) cells: a + # degenerate multistate pair can be dropped from one but not the other + # (observed wk = 0 yet E[m] > 0), which would otherwise re-zero mismatched + # populations. + naMask <- is.na(sEdge) | is.na(sBase) + sEdge[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + ret <- .RezeroGuarded(.RowMeanInformative(sEdge, informative, nSplit), + .RowMeanInformative(sBase, informative, nSplit)) + } else { + ret <- .RowMeanInformative(sEdge, informative, nSplit) } } setNames(ret, names(splits)) @@ -1031,11 +1039,16 @@ QuartetConcordance <- function( ret } else { sChar <- ifelse(denM > 0, numChar / denM, NA_real_) - ret <- .ColMeanInformative(sChar, informative, nChar) if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_) - base <- .ColMeanInformative(sBase, informative, nChar) - ret <- .RezeroGuarded(ret, base) + # Match cells, as in the edge return above. + naMask <- is.na(sChar) | is.na(sBase) + sChar[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + ret <- .RezeroGuarded(.ColMeanInformative(sChar, informative, nChar), + .ColMeanInformative(sBase, informative, nChar)) + } else { + ret <- .ColMeanInformative(sChar, informative, nChar) } ret } diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index 83c7adcd2..c1e100f16 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -276,22 +276,27 @@ test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { dat <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested: {t8,t9} - 0, 0, 1, 1, 0, 0, 1, 1, 0), 9, # crossing + 0, 0, 1, 1, 0, 0, 1, 1, 0, # crossing (binary) + 0, 0, 0, 1, 1, 1, 2, 2, 2), 9, # 3-state: exercises the random-wk path dimnames = list(paste0("t", 1:9), NULL))) # The exact hypergeometric baseline and the Monte-Carlo tip-shuffle baseline - # target the same null, so they agree to within MC error (edge and char). - exact <- QuartetConcordance(tree, dat, unit = "trit", normalize = TRUE) - set.seed(1) - mc <- QuartetConcordance(tree, dat, unit = "trit", normalize = 6000) - expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.03) - - exactC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", - normalize = TRUE) - set.seed(2) - mcC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", - normalize = 6000) - expect_lt(max(abs(exactC - mcC), na.rm = TRUE), 0.03) + # target the same null, so they agree to within MC error -- for edge and char, + # under both weightings, and crucially through the MULTISTATE pair path where a + # pair's wk and M = min(w_c, w_k) are themselves random under the null. + for (ret in c("edge", "char")) { + for (w in c(TRUE, FALSE)) { + exact <- QuartetConcordance(tree, dat, unit = "trit", return = ret, + weight = w, normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "trit", return = ret, + weight = w, normalize = 8000) + # Same cells resolve to NA under exact and MC (guards the cell-matching in + # the weight = FALSE path), and the finite values agree within MC error. + expect_identical(is.na(exact), is.na(mc)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } }) test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { From 2ac6a6c3c4679f3bda063f2405c9496c44267dc8 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 14:51:08 +0100 Subject: [PATCH 4/5] feat(concordance): enable normalize= chance correction for unit="quartet" Extend the chance-correction option to the raw quartet currency, using the same fixed-marginal (hypergeometric) null as the trit path. Per state-pair the raw concordant/decisive counts are polynomials in the 2x2 cells (no floors), so the exact expectation E[conc], E[dec] is a clean sum over the same trivariate-hypergeometric pmf (.ExpectedQuartet / .QuartetExpect); MC (.QuartetMC) reshuffles tokens and re-scores through the C++ kernel. The pooled conc/dec ratio is re-zeroed against E[conc]/E[dec] (weighted) or with the same NA-cell matching as the trit weight=FALSE path (unweighted). Removes the guard that errored on unit="quartet", normalize != FALSE. .Rezero(1, z) = 1, so a displayed split still scores exactly 1; conflicting characters go negative. normalize = FALSE is byte-identical to the published raw measure. Validated exact vs MC tip-shuffle across return x weight x binary/multistate (test-Concordance.R); all Concordance tests green. Design note in dev/plans/frac-quart-weight-agreement.md section 7 updated (raw unit DONE). This is the R reference for the forthcoming C++ port of the same expectation. Co-Authored-By: Claude Opus 4.8 --- R/Concordance.R | 198 +++++++++++++++++++---- dev/plans/frac-quart-weight-agreement.md | 19 ++- man/ConcordanceTable.Rd | 14 +- man/SiteConcordance.Rd | 14 +- tests/testthat/test-Concordance.R | 41 ++++- 5 files changed, 230 insertions(+), 56 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 7e3542ba6..ea2f353b7 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -98,19 +98,18 @@ NULL #' expectation and negative values fall below it. `ClusteringConcordance()` #' uses an analytical approximation to the expected mutual information (fast #' and generally accurate for large trees, ~200+ taxa, but neglecting -#' correlation between splits); `QuartetConcordance(unit = "trit")` uses the -#' exact hypergeometric expectation. +#' correlation between splits); `QuartetConcordance()` uses the exact +#' hypergeometric expectation (for either `unit`). #' - a positive integer `n`: estimate that expectation by Monte Carlo. #' `ClusteringConcordance()` fits each character to `n` random trees (and #' returns Monte-Carlo standard errors, more accurate for small trees where -#' the analytical approximation is biased); `QuartetConcordance(unit = -#' "trit")` averages over `n` random reassignments of each character's -#' tokens. +#' the analytical approximation is biased); `QuartetConcordance()` averages +#' over `n` random reassignments of each character's tokens. #' #' In all cases 1 corresponds to the maximum attainable value. For -#' `QuartetConcordance()`, chance correction is implemented for `unit = "trit"` -#' only, and the returned values are unclamped (they may fall below \eqn{-1}); -#' clamp to \eqn{[-1, 1]} before plotting with [QCol()] / [QACol()]. +#' `QuartetConcordance()` (either `unit`), chance-corrected values are returned +#' unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +#' plotting with [QCol()] / [QACol()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each @@ -833,13 +832,6 @@ QuartetConcordance <- function( stop("`normalize` must be FALSE, TRUE, or a positive integer.") } } - if (unit != "trit") { - # The chance baseline needs the per-state-pair cell counts, which the raw - # quartet path computes in the C++ kernel but does not expose; re-baselining - # the raw currency (and choosing its null; see the 1/3 note in the design - # spec) is deferred to the manuscript-figure review gate. - stop("`normalize` chance-correction is implemented for `unit = \"trit\"` only.") - } } tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { @@ -887,6 +879,23 @@ QuartetConcordance <- function( num <- raw_counts$concordant den <- raw_counts$decisive + # Chance correction (option): re-zero the observed concordant/decisive ratio + # against the ratio expected under the same fixed-marginal null used for the + # trit currency. Only `conc` and `dec` vary under the null (the split sizes + # and state counts are fixed), so we need E[conc] and E[dec] per (split, char), + # computed exactly from the hypergeometric pmf (no floors here, unlike trits) + # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio. + doNorm <- !isFALSE(normalize) + if (doNorm) { + base <- if (isTRUE(normalize)) { + .QuartetExpect(charInt, logiSplits) + } else { + .QuartetMC(charInt, logiSplits, normalize) + } + eNum <- base[["concordant"]] + eDen <- base[["decisive"]] + } + if (return == "default") { if (isTRUE(weight)) { # Sum numerator and denominator across sites (columns), then divide @@ -898,37 +907,164 @@ QuartetConcordance <- function( NA_real_, split_sums_num / split_sums_den ) + if (doNorm) { + bDen <- rowSums(eDen) + base <- ifelse(bDen == 0, NA_real_, rowSums(eNum) / bDen) + ret <- .RezeroGuarded(ret, base) + } } else { # Mean of ratios per site # Avoid division by zero (0/0 -> NaN -> NA handled by na.rm) ratios <- num / den # Replace NaN/Inf with NA for rowMeans calculation ratios[!is.finite(ratios)] <- NA - ret <- rowMeans(ratios, na.rm = TRUE) + if (doNorm) { + bRatios <- eNum / eDen + bRatios[!is.finite(bRatios)] <- NA + # Average observed and baseline over the same (split, char) cells. + naMask <- is.na(ratios) | is.na(bRatios) + ratios[naMask] <- NA + bRatios[naMask] <- NA + ret <- .RezeroGuarded(rowMeans(ratios, na.rm = TRUE), + rowMeans(bRatios, na.rm = TRUE)) + } else { + ret <- rowMeans(ratios, na.rm = TRUE) + } } setNames(ret, names(splits)) } else { # return = "char" - p <- num / den if (isTRUE(weight)) { - vapply( - seq_len(dim(num)[[2]]), - function(i) { - weighted.mean(num[, i] / den[, i], den[, i]) - }, - double(1) - ) + if (doNorm) { + obs <- vapply(seq_len(ncol(num)), function(i) { + d <- sum(den[, i]); if (d == 0) NA_real_ else sum(num[, i]) / d + }, double(1)) + b <- vapply(seq_len(ncol(eNum)), function(i) { + d <- sum(eDen[, i]); if (d == 0) NA_real_ else sum(eNum[, i]) / d + }, double(1)) + .RezeroGuarded(obs, b) + } else { + vapply( + seq_len(dim(num)[[2]]), + function(i) { + weighted.mean(num[, i] / den[, i], den[, i]) + }, + double(1) + ) + } } else { - vapply( - seq_len(dim(num)[[2]]), - function(i) { - mean(num[den[, i] > 0, i] / den[den[, i] > 0, i]) - }, - double(1) - ) + if (doNorm) { + sObs <- ifelse(den > 0, num / den, NA_real_) + sBase <- ifelse(eDen > 0, eNum / eDen, NA_real_) + naMask <- is.na(sObs) | is.na(sBase) + sObs[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + obs <- vapply(seq_len(ncol(sObs)), function(i) { + m <- mean(sObs[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m + }, double(1)) + b <- vapply(seq_len(ncol(sBase)), function(i) { + m <- mean(sBase[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m + }, double(1)) + .RezeroGuarded(obs, b) + } else { + vapply( + seq_len(dim(num)[[2]]), + function(i) { + mean(num[den[, i] > 0, i] / den[den[, i] > 0, i]) + }, + double(1) + ) + } + } + } +} + +# Exact / Monte-Carlo expected raw concordant + decisive quartet counts under +# the fixed-marginal null (the raw-currency analogue of `.CharTritExpect` / +# `.CharTritMC`). Per state-pair, conc = C(p,2)C(s,2) + C(q,2)C(r,2) and +# dec = conc + p*q*r*s are polynomials in the cells (no floors), so their exact +# expectation is a clean sum over the same trivariate-hypergeometric pmf. +.ExpectedQuartetCache <- new.env(hash = TRUE, parent = emptyenv()) + +.ExpectedQuartet <- function(nI, nJ, M, t) { + key <- paste(nI, nJ, M, t, sep = ",") + cached <- .ExpectedQuartetCache[[key]] + if (!is.null(cached)) { + return(cached) + } + choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2 + nOther <- t - nI - nJ + logDen <- lchoose(t, M) + eConc <- 0 + eDec <- 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 + conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r) + eConc <- eConc + prob * conc + eDec <- eDec + prob * (conc + p * q * r * s) + } + } + res <- c(eConc, eDec) + .ExpectedQuartetCache[[key]] <- res + res +} + +.QuartetExpect <- function(charInt, logiSplits) { + nSplit <- ncol(logiSplits) + nChar <- ncol(charInt) + eConc <- eDec <- 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) .ExpectedQuartet(cnt[a], cnt[b], M, tc), + double(2)) + eConc[, ci] <- eConc[, ci] + e[1, idx] + eDec[, ci] <- eDec[, ci] + e[2, idx] + } + } + } + list(concordant = eConc, decisive = eDec) +} + +.QuartetMC <- function(charInt, logiSplits, nRelabel) { + nSplit <- ncol(logiSplits) + nChar <- ncol(charInt) + accConc <- accDec <- matrix(0, nSplit, nChar) + for (s in seq_len(nRelabel)) { + shuffled <- charInt + for (ci in seq_len(nChar)) { + col <- charInt[, ci] + scored <- !is.na(col) + shuffled[scored, ci] <- sample(col[scored]) } + kc <- quartet_concordance(logiSplits, shuffled) + accConc <- accConc + kc[["concordant"]] + accDec <- accDec + kc[["decisive"]] } + list(concordant = accConc / nRelabel, decisive = accDec / nRelabel) } # Nelson-Ladiges fractional ("trit") currency for QuartetConcordance(). diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md index 6840a9142..7fe19b3d4 100644 --- a/dev/plans/frac-quart-weight-agreement.md +++ b/dev/plans/frac-quart-weight-agreement.md @@ -253,13 +253,18 @@ char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) never inside the measure. The `z → 1` blow-up is guarded (→ `NA`). `normalize` defaults to `FALSE`, so the published measure is unchanged unless the user opts in. - - **Raw `quartet` unit: deferred.** `normalize` currently errors for - `unit = "quartet"`. Re-baselining the raw currency needs the per-state-pair - cells (computed inside the C++ kernel but not exposed) and a decision on its - null (fixed-marginal vs the flat `1/3`); both are folded into the C++ port - (task 2) and the manuscript-figure review (task 3), which will compare - trit-with vs trit-without correction and decide whether the raw unit needs it - too. + - **Raw `quartet` unit: DONE.** `normalize` now works for `unit = "quartet"` + too, using the *same* fixed-marginal null. Per state-pair `conc` and `dec` + are polynomials in the cells (no floors), so the exact expectation + `E[conc]`, `E[dec]` is a clean sum over the same trivariate-hypergeometric + pmf as `.ExpectedTrit` (`.ExpectedQuartet` / `.QuartetExpect`); MC via + `.QuartetMC` reshuffles tokens and re-scores through the kernel. The pooled + `conc/dec` ratio is re-zeroed against `E[conc]/E[dec]` (weighted) or with + the same cell-matching as trits (unweighted). Validated exact≈MC across + edge/char × weight × binary/multistate (`test-Concordance.R`), and + `normalize = FALSE` is byte-identical to the published raw measure. The flat + `1/3` baseline stays rejected (it does not adapt to marginals; the + fixed-marginal null generalises it). - **`M = min(W_c,W_k)`** (settled): the pooling amount is the *shared* information, so it must be symmetric between character and split (both `return`s pool by the same `M`); `M = W_k` would make the char return weight a diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index cc617bca8..b420da41c 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -48,19 +48,19 @@ floor. expectation and negative values fall below it. \code{ClusteringConcordance()} uses an analytical approximation to the expected mutual information (fast and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the -exact hypergeometric expectation. +correlation between splits); \code{QuartetConcordance()} uses the exact +hypergeometric expectation (for either \code{unit}). \item a positive integer \code{n}: estimate that expectation by Monte Carlo. \code{ClusteringConcordance()} fits each character to \code{n} random trees (and returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's -tokens. +the analytical approximation is biased); \code{QuartetConcordance()} averages +over \code{n} random reassignments of each character's tokens. } In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} -only, and the returned values are unclamped (they may fall below \eqn{-1}); -clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned +unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 3dd1a5604..f8b3c6cc0 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -56,19 +56,19 @@ floor. expectation and negative values fall below it. \code{ClusteringConcordance()} uses an analytical approximation to the expected mutual information (fast and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the -exact hypergeometric expectation. +correlation between splits); \code{QuartetConcordance()} uses the exact +hypergeometric expectation (for either \code{unit}). \item a positive integer \code{n}: estimate that expectation by Monte Carlo. \code{ClusteringConcordance()} fits each character to \code{n} random trees (and returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's -tokens. +the analytical approximation is biased); \code{QuartetConcordance()} averages +over \code{n} random reassignments of each character's tokens. } In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} -only, and the returned values are unclamped (they may fall below \eqn{-1}); -clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned +unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the number of quartets they are decisive for.} diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index c1e100f16..b455f6012 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -252,10 +252,10 @@ test_that("QuartetConcordance() unit = 'trit' normalize contract", { expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"), "positive integer") - # Chance correction is implemented for the trit currency only, for now. - expect_error( - QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE), - "trit") + # Chance correction also works for the raw quartet currency. + qc <- QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE) + expect_type(qc, "double") + expect_false(anyNA(qc)) # An identical (displayed) split still scores exactly 1 after re-zeroing: # .Rezero(1, z) == 1, so chance correction lifts the floor without moving the @@ -299,6 +299,39 @@ test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { } }) +test_that("QuartetConcordance() unit = 'quartet' exact baseline matches MC", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + dat <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} + 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested + 0, 0, 1, 1, 0, 0, 1, 1, 0, # crossing (binary) + 0, 0, 0, 1, 1, 1, 2, 2, 2), 9, # 3-state: multistate decisive path + dimnames = list(paste0("t", 1:9), NULL))) + + # The raw-currency chance baseline (exact E[conc]/E[dec] from the + # hypergeometric pmf) matches the Monte-Carlo tip-shuffle baseline within MC + # error, across edge/char and both weightings, incl. the multistate path. + for (ret in c("edge", "char")) { + for (w in c(TRUE, FALSE)) { + exact <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, normalize = 8000) + expect_identical(is.na(exact), is.na(mc)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } + + # normalize = FALSE leaves the published raw measure untouched. + for (ret in c("edge", "char")) for (w in c(TRUE, FALSE)) { + expect_identical( + QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w), + QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w, + normalize = FALSE)) + } +}) + test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") support <- MatrixToPhyDat(matrix( # agrees with the tree: displays {t7,t8,t9} From 8461fc6bd447ff04c36974b636d7294037439b1d Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 15:14:44 +0100 Subject: [PATCH 5/5] perf(concordance): port chance-baseline expectation to C++ Move the exact hypergeometric expectation for QuartetConcordance(normalize = TRUE) from R to C++ (src/concordance_expect.cpp): `quartet_expect` returns E[conc], E[dec]; `trit_expect` returns E[m], E[m*A/wk], E[m*A/wc]. Both sum the per-state-pair trivariate-hypergeometric double sum, hoisting the log-choose terms out of the M/p/r loops (the dominant cost). The R exact helpers (.ExpectedTrit/.CharTritExpect/.ExpectedQuartet and the per-char .QuartetExpect loop) are replaced by thin wrappers / a single all-characters call in .TritConcordance; the Monte-Carlo (normalize = ) path stays in R. New Rcpp exports are registered in the manual R_CallMethodDef table (src/TreeSearch-init.c; the package uses R_useDynamicSymbols(FALSE)). Speedup at split-support scale (48 tips, 96 chars): exact normalize=TRUE quartet 2.88s -> 0.14s (~20x), trit 3.48s -> 0.45s (~8x); a full 1000-matrix run drops from ~100 min to ~10 min. Validated bit-for-bit (max |C++ - R| ~1e-14) against a self-contained R reference across binary/multistate/missing data (dev/benchmarks/frac-quart/cpp_expect_parity.R); the exact-vs-MC end-to-end tests in test-Concordance.R now exercise the C++ path and pass. Co-Authored-By: Claude Opus 4.8 --- R/Concordance.R | 167 ++------------ R/RcppExports.R | 8 + dev/benchmarks/frac-quart/cpp_expect_parity.R | 108 +++++++++ src/RcppExports.cpp | 24 ++ src/TreeSearch-init.c | 4 + src/concordance_expect.cpp | 215 ++++++++++++++++++ 6 files changed, 382 insertions(+), 144 deletions(-) create mode 100644 dev/benchmarks/frac-quart/cpp_expect_parity.R create mode 100644 src/concordance_expect.cpp diff --git a/R/Concordance.R b/R/Concordance.R index ea2f353b7..51a484d22 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -980,73 +980,15 @@ QuartetConcordance <- function( } } -# Exact / Monte-Carlo expected raw concordant + decisive quartet counts under -# the fixed-marginal null (the raw-currency analogue of `.CharTritExpect` / -# `.CharTritMC`). Per state-pair, conc = C(p,2)C(s,2) + C(q,2)C(r,2) and -# dec = conc + p*q*r*s are polynomials in the cells (no floors), so their exact -# expectation is a clean sum over the same trivariate-hypergeometric pmf. -.ExpectedQuartetCache <- new.env(hash = TRUE, parent = emptyenv()) - -.ExpectedQuartet <- function(nI, nJ, M, t) { - key <- paste(nI, nJ, M, t, sep = ",") - cached <- .ExpectedQuartetCache[[key]] - if (!is.null(cached)) { - return(cached) - } - choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2 - nOther <- t - nI - nJ - logDen <- lchoose(t, M) - eConc <- 0 - eDec <- 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 - conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r) - eConc <- eConc + prob * conc - eDec <- eDec + prob * (conc + p * q * r * s) - } - } - res <- c(eConc, eDec) - .ExpectedQuartetCache[[key]] <- res - res -} - +# Expected raw concordant + decisive quartet counts under the fixed-marginal +# null. The exact expectation (per state-pair, an exact sum over the +# trivariate-hypergeometric pmf) is computed in C++ (`quartet_expect` in +# src/concordance_expect.cpp) for speed; `.QuartetExpect` is a thin R wrapper. +# The `.QuartetMC` Monte-Carlo baseline reshuffles tokens and re-scores through +# the C++ `quartet_concordance` kernel. (The C++ path is validated bit-for-bit +# against the earlier R implementation; see dev/benchmarks/frac-quart.) .QuartetExpect <- function(charInt, logiSplits) { - nSplit <- ncol(logiSplits) - nChar <- ncol(charInt) - eConc <- eDec <- 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) .ExpectedQuartet(cnt[a], cnt[b], M, tc), - double(2)) - eConc[, ci] <- eConc[, ci] + e[1, idx] - eDec[, ci] <- eDec[, ci] + e[2, idx] - } - } - } - list(concordant = eConc, decisive = eDec) + quartet_expect(logiSplits, charInt) } .QuartetMC <- function(charInt, logiSplits, nRelabel) { @@ -1109,6 +1051,15 @@ QuartetConcordance <- function( doNorm <- !isFALSE(normalize) if (doNorm) { baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar) + if (isTRUE(normalize)) { + # Exact hypergeometric expectation for all characters at once (C++; + # `trit_expect` in src/concordance_expect.cpp). Uninformative characters + # yield zero pools, matching the observed loop's `wcTot > 0` guard. + be <- trit_expect(logiSplits, charInt) + baseNumEdge <- be[["numEdge"]] + baseNumChar <- be[["numChar"]] + baseDenM <- be[["denM"]] + } } for (ci in seq_len(nChar)) { @@ -1118,12 +1069,8 @@ QuartetConcordance <- function( numChar[, ci] <- obs[["numChar"]] denM[, ci] <- obs[["denM"]] wcTot[ci] <- obs[["wcTot"]] - if (doNorm && obs[["wcTot"]] > 0) { - base <- if (isTRUE(normalize)) { - .CharTritExpect(col, logiSplits, pos) # exact hypergeometric - } else { - .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo tip-shuffle - } + if (doNorm && !isTRUE(normalize) && obs[["wcTot"]] > 0) { + base <- .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo shuffle baseNumEdge[, ci] <- base[["numEdge"]] baseNumChar[, ci] <- base[["numChar"]] baseDenM[, ci] <- base[["denM"]] @@ -1229,78 +1176,10 @@ QuartetConcordance <- function( list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot) } -# Exact expected trit contributions under the fixed-marginal (hypergeometric) -# null: each split's scored side-A size `M` and each pair's state counts -# `n_i, n_j` are held fixed while the character's tokens are reassigned at -# random across the `t` scored leaves. Only `A`, and (for multistate pairs) -# `wk` and `m`, vary; so we accumulate E[m], E[m * A / wk] and E[m * A / wc] -# from the exact pmf, matching `.CharTritContrib`'s pool definitions. -.CharTritExpect <- function(col, logiSplits, pos) { - nSplit <- ncol(logiSplits) - numEdge <- numChar <- denM <- numeric(nSplit) - scored <- !is.na(col) - tc <- sum(scored) - mSideA <- colSums(logiSplits & scored) # M per split - uM <- unique(mSideA) - idx <- match(mSideA, uM) - states <- sort(unique(col[scored])) - nStates <- length(states) - if (nStates >= 2L) { - cnt <- tabulate(match(col[scored], states), nStates) # state counts - for (a in seq_len(nStates - 1L)) { - for (b in seq(a + 1L, nStates)) { - et <- vapply(uM, function(M) .ExpectedTrit(cnt[a], cnt[b], M, tc), - double(3)) - denM <- denM + et[1, idx] - numEdge <- numEdge + et[2, idx] - numChar <- numChar + et[3, idx] - } - } - } - list(numEdge = numEdge, numChar = numChar, denM = denM) -} - -.ExpectedTritCache <- new.env(hash = TRUE, parent = emptyenv()) - -# Exact E[m], E[m * A / wk], E[m * A / wc] for one state-pair, over the -# trivariate hypergeometric placement of `nI` state-i and `nJ` state-j leaves -# (and `t - nI - nJ` other leaves) among the `M` leaves on side A of a split. -.ExpectedTrit <- function(nI, nJ, M, t) { - key <- paste(nI, nJ, M, t, sep = ",") - cached <- .ExpectedTritCache[[key]] - if (!is.null(cached)) { - return(cached) - } - pos1 <- function(z) if (z < 0) 0 else z - nOther <- t - nI - nJ - tP <- nI + nJ # taxa in this pair - wc <- pos1(nI - 1) * pos1(nJ - 1) # pair's character content (fixed) - logDen <- lchoose(t, M) - acc <- c(0, 0, 0) # E[m], E[m*A/wk], E[m*A/wc] - for (p in 0:min(nI, M)) { - for (r in 0:min(nJ, M - p)) { - oA <- M - p - r # other-state leaves on side A - 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 # this pair's taxa on side A (random) - A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1) - wk <- pos1(mA - 1) * pos1(tP - mA - 1) # pair's split content (random) - m <- min(wc, wk) - acc[1] <- acc[1] + prob * m - if (wk > 0) acc[2] <- acc[2] + prob * m * A / wk - if (wc > 0) acc[3] <- acc[3] + prob * m * A / wc - } - } - .ExpectedTritCache[[key]] <- acc - acc -} +# The exact trit expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++ +# (`trit_expect` in src/concordance_expect.cpp) for all characters at once, and +# consumed directly by `.TritConcordance`; there is no per-character R exact +# helper. The Monte-Carlo baseline below is the opt-in `normalize = ` path. # Monte-Carlo baseline: average `.CharTritContrib` over `nRelabel` random # reassignments of the character's tokens across its scored leaves (the split, diff --git a/R/RcppExports.R b/R/RcppExports.R index afce48989..4d29eae1d 100644 --- a/R/RcppExports.R +++ b/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 diff --git a/dev/benchmarks/frac-quart/cpp_expect_parity.R b/dev/benchmarks/frac-quart/cpp_expect_parity.R new file mode 100644 index 000000000..0e8afba7f --- /dev/null +++ b/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") diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dc52e9b59..f46299b2e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -31,6 +31,30 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// quartet_expect +List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters); +RcppExport SEXP _TreeSearch_quartet_expect(SEXP splitsSEXP, SEXP charactersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP); + rcpp_result_gen = Rcpp::wrap(quartet_expect(splits, characters)); + return rcpp_result_gen; +END_RCPP +} +// trit_expect +List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters); +RcppExport SEXP _TreeSearch_trit_expect(SEXP splitsSEXP, SEXP charactersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP); + rcpp_result_gen = Rcpp::wrap(trit_expect(splits, characters)); + return rcpp_result_gen; +END_RCPP +} // expected_mi double expected_mi(const IntegerVector& ni, const IntegerVector& nj); RcppExport SEXP _TreeSearch_expected_mi(SEXP niSEXP, SEXP njSEXP) { diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c index 7b166a616..890fd4255 100644 --- a/src/TreeSearch-init.c +++ b/src/TreeSearch-init.c @@ -21,6 +21,8 @@ extern SEXP _TreeSearch_mi_key(SEXP, SEXP); // extern SEXP _TreeSearch_astar_search_r(SEXP, SEXP, SEXP); extern SEXP _TreeSearch_quartet_concordance(SEXP, SEXP); +extern SEXP _TreeSearch_quartet_expect(SEXP, SEXP); +extern SEXP _TreeSearch_trit_expect(SEXP, SEXP); extern SEXP _TreeSearch_ts_fitch_score(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_na_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); @@ -75,6 +77,8 @@ static const R_CallMethodDef callMethods[] = { // {"_TreeSearch_astar_search_r", (DL_FUNC) &_TreeSearch_astar_search_r, 3}, {"_TreeSearch_quartet_concordance",(DL_FUNC) &_TreeSearch_quartet_concordance, 2}, + {"_TreeSearch_quartet_expect", (DL_FUNC) &_TreeSearch_quartet_expect, 2}, + {"_TreeSearch_trit_expect", (DL_FUNC) &_TreeSearch_trit_expect, 2}, {"_TreeSearch_ts_fitch_score", (DL_FUNC) &_TreeSearch_ts_fitch_score, 12}, {"_TreeSearch_ts_na_char_steps", (DL_FUNC) &_TreeSearch_ts_na_char_steps, 5}, {"_TreeSearch_ts_char_steps", (DL_FUNC) &_TreeSearch_ts_char_steps, 5}, diff --git a/src/concordance_expect.cpp b/src/concordance_expect.cpp new file mode 100644 index 000000000..4faefac8e --- /dev/null +++ b/src/concordance_expect.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include + +using namespace Rcpp; + +// Expected concordance under the fixed-marginal (hypergeometric) null used by +// QuartetConcordance(normalize = TRUE). For one character each token is +// reassigned at random across the scored leaves, holding the state counts and +// the split sizes fixed. For a state-pair (i, j) with counts nI, nJ over t +// scored leaves, of which M are on side A of a split, the 2x2 cell vector +// (p, q, r, s) is trivariate-hypergeometric: +// P(p, r) = C(nI,p) C(nJ,r) C(t-nI-nJ, M-p-r) / C(t, M) +// with q = nI - p, s = nJ - r, and mA = p + r on side A. These kernels sum the +// per-pair contributions exactly (double sum over p, r), mirroring the R +// reference `.ExpectedQuartet`/`.QuartetExpect` and `.ExpectedTrit`/ +// `.CharTritExpect` in R/Concordance.R. The p-outer / r-inner accumulation +// order matches the R code so the results agree to machine precision. + +// C(z, 2) with the same z < 2 -> 0 guard as the R `choose2`. +static inline double choose2(double z) { + return z < 2.0 ? 0.0 : z * (z - 1.0) / 2.0; +} + +// (z - 1) floored at 0, matching the R `pos(z - 1)`. +static inline double posm1(double z) { + return z < 1.0 ? 0.0 : z - 1.0; +} + +// Per-character scored state counts + per-split side-A size, shared by both +// kernels. `stateOf[t]` is the compact 0-based state index of a scored leaf, or +// -1 if the leaf is NA for this character. +struct CharCounts { + int tc; // number of scored leaves + std::vector cnt; // count of each compact state + std::vector mSideA; // per split: scored leaves on side A +}; + +static CharCounts char_counts(const LogicalMatrix& splits, + const IntegerMatrix& characters, int c) { + const int n_taxa = splits.nrow(); + const int n_splits = splits.ncol(); + CharCounts out; + out.tc = 0; + out.mSideA.assign(n_splits, 0); + + // Map raw state values to compact 0-based indices (sorted ascending, as R's + // sort(unique(...)) does; the pair loop is symmetric so order only affects + // which is "i" vs "j", not the sums). + std::vector raw; + raw.reserve(n_taxa); + std::vector stateOf(n_taxa, -1); + for (int t = 0; t < n_taxa; ++t) { + int st = characters(t, c); + if (IntegerVector::is_na(st)) continue; + ++out.tc; + raw.push_back(st); + } + std::sort(raw.begin(), raw.end()); + raw.erase(std::unique(raw.begin(), raw.end()), raw.end()); + out.cnt.assign(raw.size(), 0); + + for (int t = 0; t < n_taxa; ++t) { + int st = characters(t, c); + if (IntegerVector::is_na(st)) continue; + int idx = int(std::lower_bound(raw.begin(), raw.end(), st) - raw.begin()); + stateOf[t] = idx; + ++out.cnt[idx]; + for (int s = 0; s < n_splits; ++s) { + if (splits(t, s)) ++out.mSideA[s]; + } + } + return out; +} + +// Exact E[conc], E[dec] for one state-pair over all needed side-A sizes. +// Returns, for each M in 0..tc, the pair's (Econc, Edec); callers pick M[split]. +static void expected_quartet_by_M(int nI, int nJ, int tc, + std::vector& eConc, + std::vector& eDec) { + const int nOther = tc - nI - nJ; + eConc.assign(tc + 1, 0.0); + eDec.assign(tc + 1, 0.0); + // Hoist the log-choose terms out of the M / p / r loops (the dominant cost). + std::vector lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1); + for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p); + for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r); + for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o); + for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M); + for (int M = 0; M <= tc; ++M) { + const double lt = lcT[M]; + double aConc = 0.0, aDec = 0.0; + const int pMax = std::min(nI, M); + for (int p = 0; p <= pMax; ++p) { + const int q = nI - p; + const double c_p = choose2(p), c_q = choose2(q), lp = lcI[p]; + const int rMax = std::min(nJ, M - p); + for (int r = 0; r <= rMax; ++r) { + const int oA = M - p - r; + if (oA < 0 || oA > nOther) continue; + const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt); + if (prob <= 0.0) continue; + const int s = nJ - r; + const double conc = c_p * choose2(s) + c_q * choose2(r); + aConc += prob * conc; + aDec += prob * (conc + double(p) * q * r * s); + } + } + eConc[M] = aConc; + eDec[M] = aDec; + } +} + +// [[Rcpp::export]] +List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters) { + const int n_splits = splits.ncol(); + const int n_chars = characters.ncol(); + NumericMatrix eConc(n_splits, n_chars), eDec(n_splits, n_chars); + + for (int c = 0; c < n_chars; ++c) { + const CharCounts cc = char_counts(splits, characters, c); + const int nStates = int(cc.cnt.size()); + if (nStates < 2) continue; + std::vector ec, ed; + for (int a = 0; a < nStates - 1; ++a) { + for (int b = a + 1; b < nStates; ++b) { + expected_quartet_by_M(cc.cnt[a], cc.cnt[b], cc.tc, ec, ed); + for (int s = 0; s < n_splits; ++s) { + const int M = cc.mSideA[s]; + eConc(s, c) += ec[M]; + eDec(s, c) += ed[M]; + } + } + } + } + return List::create(_["concordant"] = eConc, _["decisive"] = eDec); +} + +// Exact E[m], E[m*A/wk], E[m*A/wc] for one trit state-pair over all side-A +// sizes. wc = (nI-1)+ (nJ-1)+ is fixed; wk = (mA-1)+ (tP-mA-1)+ and +// m = min(wc, wk) vary with mA = p + r. Mirrors R `.ExpectedTrit`. +static void expected_trit_by_M(int nI, int nJ, int tc, + std::vector& eM, + std::vector& eMAwk, + std::vector& eMAwc) { + const int nOther = tc - nI - nJ; + const int tP = nI + nJ; + const double wc = posm1(nI) * posm1(nJ); + eM.assign(tc + 1, 0.0); + eMAwk.assign(tc + 1, 0.0); + eMAwc.assign(tc + 1, 0.0); + // Hoist the log-choose terms out of the M / p / r loops (the dominant cost). + std::vector lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1); + for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p); + for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r); + for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o); + for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M); + for (int M = 0; M <= tc; ++M) { + const double lt = lcT[M]; + double aM = 0.0, aWk = 0.0, aWc = 0.0; + const int pMax = std::min(nI, M); + for (int p = 0; p <= pMax; ++p) { + const int q = nI - p; + const double lp = lcI[p]; + const int rMax = std::min(nJ, M - p); + for (int r = 0; r <= rMax; ++r) { + const int oA = M - p - r; + if (oA < 0 || oA > nOther) continue; + const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt); + if (prob <= 0.0) continue; + const int s = nJ - r; + const int mA = p + r; + const double A = posm1(p) * posm1(s) + posm1(q) * posm1(r); + const double wk = posm1(mA) * posm1(tP - mA); + const double m = std::min(wc, wk); + aM += prob * m; + if (wk > 0.0) aWk += prob * m * A / wk; + if (wc > 0.0) aWc += prob * m * A / wc; + } + } + eM[M] = aM; + eMAwk[M] = aWk; + eMAwc[M] = aWc; + } +} + +// [[Rcpp::export]] +List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters) { + const int n_splits = splits.ncol(); + const int n_chars = characters.ncol(); + NumericMatrix numEdge(n_splits, n_chars); + NumericMatrix numChar(n_splits, n_chars); + NumericMatrix denM(n_splits, n_chars); + + for (int c = 0; c < n_chars; ++c) { + const CharCounts cc = char_counts(splits, characters, c); + const int nStates = int(cc.cnt.size()); + if (nStates < 2) continue; + std::vector eM, eWk, eWc; + for (int a = 0; a < nStates - 1; ++a) { + for (int b = a + 1; b < nStates; ++b) { + expected_trit_by_M(cc.cnt[a], cc.cnt[b], cc.tc, eM, eWk, eWc); + for (int s = 0; s < n_splits; ++s) { + const int M = cc.mSideA[s]; + denM(s, c) += eM[M]; + numEdge(s, c) += eWk[M]; + numChar(s, c) += eWc[M]; + } + } + } + } + return List::create(_["numEdge"] = numEdge, _["numChar"] = numChar, + _["denM"] = denM); +}