Support group=FALSE for multihop connectivity

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)

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,24 +91,29 @@ 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)
   s <- Matrix::summary(Matrix::drop0(eff))
   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

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.