diff --git a/R/cosine.R b/R/cosine.R index acf43df..c2f4021 100644 --- a/R/cosine.R +++ b/R/cosine.R @@ -67,8 +67,26 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') { #' partner sets grow quickly, \code{min_frac} prunes each layer to types that #' receive at least that fraction of effective input (per-type, applied at #' every hop including the final one); it defaults to a small non-zero value -#' when \code{nhops>0}. A grouping column (e.g. the default \code{group="type"}) -#' is required for \code{nhops>0}. +#' when \code{nhops>0}. +#' +#' Partner neurons lacking a group label (e.g. an untyped neuron when +#' \code{group="type"}) are handled differently depending on where they sit. In +#' an \emph{intermediate} layer they form their own singleton groups and are +#' thresholded per neuron, so connectivity still propagates through a poorly +#' typed layer to well typed neurons beyond it - which is often the whole point +#' of looking more than one hop away. In the \emph{final} layer they are +#' dropped, since an untyped neuron cannot serve as a shared feature; this +#' matches the behaviour for direct partners when \code{nhops=0}. +#' +#' \code{group=FALSE} works for multihop too: each neuron is then treated as +#' its own group, so the features are individual n-th order partner neurons +#' and \code{min_frac} becomes a per-neuron rather than a per-type cut. This is +#' the only option for datasets without metadata (e.g. FANC) and is useful for +#' columnar neurons, but as for \code{nhops=0} it only really makes sense +#' within a single dataset, since individual neurons do not correspond across +#' datasets. Note that per-neuron effective weights are smaller than their +#' per-type aggregates, so you may want a smaller \code{min_frac} than the +#' default. #' #' The \code{labRow} argument is most conveniently specified as a length 1 #' string to be interpolated by \code{\link[glue]{glue}}; this will happen in @@ -103,7 +121,8 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') { #' @param min_frac Per-type fractional threshold (default \code{0.005}) used #' when \code{nhops>0} to prune each layer (including the final one) to partner #' types receiving at least this fraction of effective input. A scalar or a -#' vector with one value per hop. Ignored when \code{nhops=0}. +#' vector with one value per hop. Becomes a per-neuron cut when +#' \code{group=FALSE} (see \bold{details}). Ignored when \code{nhops=0}. #' @param remove_query Whether to exclude the query neurons from the partners at #' every hop, both as intermediate interneurons and as final targets (default #' \code{FALSE}). Only relevant when \code{nhops>0}. @@ -387,9 +406,6 @@ multi_connection_table <- function(ids, partners=c("inputs", "outputs"), if(isTRUE(group)) group='type' partners=match.arg(partners, several.ok = T) - if(nhops>0 && !is.character(group)) - stop("Multihop connectivity (nhops>0) requires a grouping column ", - "(e.g. group='type') to define partner layers across hops.") kk=keys(ids) if(length(partners)>1) { l=sapply(partners, simplify = F, function(p) diff --git a/R/multihop.R b/R/multihop.R index d3fcf41..1170cb4 100644 --- a/R/multihop.R +++ b/R/multihop.R @@ -35,7 +35,7 @@ multihop_effective_matrix <- function(dskeys, partners, nhops, threshold, tbl <- cf_partners(frontier, threshold = threshold, partners = partners, MoreArgs = MoreArgs) if (is.null(tbl) || nrow(tbl) == 0) return(NULL) - if (!group %in% colnames(tbl)) + if (is.character(group) && !group %in% colnames(tbl)) stop("Grouping column `", group, "` not present in cf_partners result!") Mn <- coconat::colScaleM(.mh_adjacency(tbl)) @@ -50,17 +50,31 @@ multihop_effective_matrix <- function(dskeys, partners, nhops, threshold, else coconat::effective_connectivity(list(running, step), normalise = FALSE) - # type of each new-frontier neuron (partner side of this hop) - ntypes <- tbl[[group]][!duplicated(tbl[[newfrontier_key]])] - names(ntypes) <- tbl[[newfrontier_key]][!duplicated(tbl[[newfrontier_key]])] - ntypes <- ntypes[colnames(running)] + # group label of each new-frontier neuron (partner side of this hop). + # group=FALSE means each neuron is its own group, so pruning becomes a + # per-neuron rather than per-type cut and the terminal layer stays ungrouped. + ntypes <- if (isFALSE(group)) { + stats::setNames(colnames(running), colnames(running)) + } else { + nt <- tbl[[group]][!duplicated(tbl[[newfrontier_key]])] + names(nt) <- tbl[[newfrontier_key]][!duplicated(tbl[[newfrontier_key]])] + nt[colnames(running)] + } if (h <= nhops) { - # intermediate layer: prune frontier by type (selection only) - G <- coconat::grouping_matrix(colnames(running), ntypes) - grp <- running %*% G + # Intermediate layer: prune the frontier by group (selection only). With + # group=FALSE the grouping is the identity so we skip the multiplication. + # Neurons with no group label form their own singleton groups rather than + # being dropped, so connectivity still propagates through a poorly typed + # intermediate layer to well typed neurons beyond it; they are simply + # thresholded per neuron instead of per type. + ptypes <- ntypes + if (anyNA(ptypes)) + ptypes[is.na(ptypes)] <- colnames(running)[is.na(ptypes)] + grp <- if (isFALSE(group)) running + else running %*% coconat::grouping_matrix(colnames(running), ptypes) keep_types <- colnames(grp)[apply(as.matrix(grp), 2, max) >= mf[h]] - surviving <- colnames(running)[ntypes %in% keep_types] + surviving <- colnames(running)[ptypes %in% keep_types] if (remove_query) surviving <- setdiff(surviving, query_keys) running <- running[, surviving, drop = FALSE] @@ -77,16 +91,20 @@ multihop_effective_matrix <- function(dskeys, partners, nhops, threshold, running <- running[, keepcols, drop = FALSE] far_types <- far_types[keepcols] } - # group far neurons to type (per-neuron normalisation already done) - eff <- running %*% coconat::grouping_matrix(colnames(running), far_types) - # final per-type cut + # group far neurons to type (per-neuron normalisation already done). With + # group=FALSE the terminal layer is left at neuron resolution. + eff <- if (isFALSE(group)) running + else running %*% coconat::grouping_matrix(colnames(running), far_types) + # final cut eff <- Matrix::drop0(eff * (eff >= mf[nhops + 1L])) eff } -# Melt an effective query x target-type matrix into a cf_partners-like long -# table for one direction, with the query on the appropriate key column and the -# target type in `group`. Feeds straight into multi_cosine_matrix. +# Melt an effective query x target matrix into a cf_partners-like long table for +# one direction, with the query on the appropriate key column and the target in +# `group`. Feeds straight into multi_cosine_matrix. When group=FALSE the target +# is a neuron key, which multi_cosine_matrix reads from post_key/pre_key +# directly, so no grouping column is added. .mh_matrix2df <- function(eff, dataset, partners, group = "type") { if (is.null(eff) || length(eff) == 0 || sum(eff != 0) == 0) return(NULL) @@ -94,7 +112,8 @@ multihop_effective_matrix <- function(dskeys, partners, nhops, threshold, qk <- rownames(eff)[s$i] tp <- colnames(eff)[s$j] df <- data.frame(stringsAsFactors = FALSE, weight = s$x) - df[[group]] <- tp + if (!isFALSE(group)) + df[[group]] <- tp # query lives in pre_key for outputs, post_key for inputs if (partners == "outputs") { df$pre_key <- qk diff --git a/man/cf_cosine_plot.Rd b/man/cf_cosine_plot.Rd index f716a4a..c48bffd 100644 --- a/man/cf_cosine_plot.Rd +++ b/man/cf_cosine_plot.Rd @@ -62,7 +62,8 @@ partners; \code{1} is the \dQuote{one-hop} (2nd-order) case, etc. See \item{min_frac}{Per-type fractional threshold (default \code{0.005}) used when \code{nhops>0} to prune each layer (including the final one) to partner types receiving at least this fraction of effective input. A scalar or a -vector with one value per hop. Ignored when \code{nhops=0}.} +vector with one value per hop. Becomes a per-neuron cut when +\code{group=FALSE} (see \bold{details}). Ignored when \code{nhops=0}.} \item{remove_query}{Whether to exclude the query neurons from the partners at every hop, both as intermediate interneurons and as final targets (default @@ -169,8 +170,26 @@ examples. partner sets grow quickly, \code{min_frac} prunes each layer to types that receive at least that fraction of effective input (per-type, applied at every hop including the final one); it defaults to a small non-zero value - when \code{nhops>0}. A grouping column (e.g. the default \code{group="type"}) - is required for \code{nhops>0}. + when \code{nhops>0}. + + Partner neurons lacking a group label (e.g. an untyped neuron when + \code{group="type"}) are handled differently depending on where they sit. In + an \emph{intermediate} layer they form their own singleton groups and are + thresholded per neuron, so connectivity still propagates through a poorly + typed layer to well typed neurons beyond it - which is often the whole point + of looking more than one hop away. In the \emph{final} layer they are + dropped, since an untyped neuron cannot serve as a shared feature; this + matches the behaviour for direct partners when \code{nhops=0}. + + \code{group=FALSE} works for multihop too: each neuron is then treated as + its own group, so the features are individual n-th order partner neurons + and \code{min_frac} becomes a per-neuron rather than a per-type cut. This is + the only option for datasets without metadata (e.g. FANC) and is useful for + columnar neurons, but as for \code{nhops=0} it only really makes sense + within a single dataset, since individual neurons do not correspond across + datasets. Note that per-neuron effective weights are smaller than their + per-type aggregates, so you may want a smaller \code{min_frac} than the + default. The \code{labRow} argument is most conveniently specified as a length 1 string to be interpolated by \code{\link[glue]{glue}}; this will happen in diff --git a/tests/testthat/test-multihop.R b/tests/testthat/test-multihop.R index b4b2be4..4874a14 100644 --- a/tests/testthat/test-multihop.R +++ b/tests/testthat/test-multihop.R @@ -85,6 +85,84 @@ test_that("remove_query drops query neurons from the target set", { expect_false("Q" %in% colnames(eff1)) }) +test_that("untyped intermediate neurons still propagate connectivity", { + # q1 -> i1(type A), i2(UNTYPED) ; i1 -> t1(X) ; i2 -> t2(Y). The whole point + # of multihop is to reach a well typed layer through a poorly typed one, so + # the untyped interneuron must not silently discard the path to Y. + edges <- data.frame( + pre_key = c("fw:1","fw:1","fw:11","fw:12"), + post_key = c("fw:11","fw:12","fw:21","fw:22"), + weight = c(10,10,8,8), stringsAsFactors=FALSE) + mk <- function(typ) function(ids, threshold=1L, partners="outputs", + MoreArgs=list(), ...) { + e <- edges[edges$pre_key %in% as.character(ids) & edges$weight>=threshold,,drop=FALSE] + if(nrow(e)==0) return(NULL) + e$type <- unname(typ[e$post_key]); e$dataset <- "flywire"; e + } + + testthat::local_mocked_bindings( + cf_partners = mk(c("fw:11"="A","fw:12"=NA,"fw:21"="X","fw:22"="Y"))) + eff <- coconatfly:::multihop_effective_matrix( + "fw:1", "outputs", 1L, 1L, 0, "type") + expect_setequal(colnames(eff), c("X","Y")) +}) + +test_that("untyped terminal partners are dropped (as at nhops=0)", { + edges <- data.frame( + pre_key = c("fw:1","fw:1","fw:11","fw:12"), + post_key = c("fw:11","fw:12","fw:21","fw:22"), + weight = c(10,10,8,8), stringsAsFactors=FALSE) + typ <- c("fw:11"="A","fw:12"="B","fw:21"="X","fw:22"=NA) + testthat::local_mocked_bindings( + cf_partners = function(ids, threshold=1L, partners="outputs", + MoreArgs=list(), ...) { + e <- edges[edges$pre_key %in% as.character(ids),,drop=FALSE] + if(nrow(e)==0) return(NULL) + e$type <- unname(typ[e$post_key]); e$dataset <- "flywire"; e + }) + eff <- coconatfly:::multihop_effective_matrix( + "fw:1", "outputs", 1L, 1L, 0, "type") + # an untyped target cannot serve as a shared feature, so only X remains + expect_equal(colnames(eff), "X") +}) + +test_that("group=FALSE keeps targets at neuron resolution", { + testthat::local_mocked_bindings(cf_partners = mock_cf_partners) + # same network/values as the hand-computed one-hop test, but ungrouped: the + # single X neuron is fw:21 and the single Y neuron is fw:22 + eff <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), partners="outputs", nhops=1L, + threshold=1L, min_frac=0, group=FALSE) + m <- as.matrix(eff) + expect_setequal(colnames(m), c("fw:21","fw:22")) + m <- m[c("fw:1","fw:2"), c("fw:21","fw:22")] + expect_equal(unname(m), matrix(c(2/3, 2/9, 1/3, 7/9), nrow=2, byrow=TRUE)) +}) + +test_that("group=FALSE prunes the frontier per neuron", { + testthat::local_mocked_bindings(cf_partners = mock_cf_partners) + # per-neuron effective input fractions to the interneurons: i11 max 2/3, + # i12 max 1. min_frac 0.7 prunes i11, leaving only q2 -> i12 -> fw:22 + eff <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), partners="outputs", nhops=1L, + threshold=1L, min_frac=0.7, group=FALSE) + m <- as.matrix(eff) + expect_equal(colnames(m), "fw:22") + expect_equal(unname(m["fw:2","fw:22"]), 1) +}) + +test_that(".mh_matrix2df omits the grouping column when group=FALSE", { + m <- Matrix::Matrix(c(2/3, 1/3), nrow=1, + dimnames=list("fw:1", c("fw:21","fw:22"))) + df <- coconatfly:::.mh_matrix2df(m, dataset="flywire", partners="outputs", + group=FALSE) + expect_false("FALSE" %in% colnames(df)) + expect_false("type" %in% colnames(df)) + # multi_cosine_matrix reads the ungrouped target straight from post_key + expect_setequal(df$post_key, c("fw:21","fw:22")) + expect_equal(unique(df$pre_key), "fw:1") +}) + test_that("nhops=2 handles an interneuron dead end (row alignment)", { # q1 -> {i1(A), i2(B)} ; i1 -> j1(X) ; i2 has NO outputs (dead end) ; # j1 -> t1(Z). With alignment, i2's missing onward row must become zero.