[FEAT] fcount(drop = FALSE) and collap(drop = FALSE)

[egoipse]

Hi,

I'm a great fan of collapse and fastverse.

Sometimes, however, I miss some features that would make collapse's functions more powerful.

One of them is the capability to keep factor levels in collap or other function even when some of those levels are not present in final output.

A similar request was made in another issue.

The tidyverse turtle allows to keep levels with group_by's argument drop = FALSE.

Same functionality would be great in collap, fcount, and fgroup_by function to avoid tidyverse pretty slower equivalents.

[SebKrantz]

Hi @egoipse. Thanks for your kind words. I believe this issue has already been requested in a similar form here. The problem is that collapse's algorithms are designed for multivariate data, and making such an univariate exemption would be difficult. What does dplyr do on multivariate grouping data if one variable is a factor? But perhaps @claude has some ideas on how we could efficiently integrate this?

[claude]

Claude finished @SebKrantz's task in 2m 21s —— View job

---

Analysis: drop = FALSE for fcount/collap/fgroup_by

• Gather context on the request
• Explore GRP system and factor handling
• Examine fcount and collap implementations
• Research how dplyr handles drop = FALSE for multivariate grouping
• Propose integration approaches

---

Here's my analysis of the current architecture and feasible integration paths.

What the codebase already has

GRP.factor (R/GRP.R:335) already supports drop = FALSE — and it is already the default (drop = FALSE keeps all levels). When drop = TRUE, it calls fdroplevels to discard unobserved levels. So single-factor grouping via GRP.factor already preserves empty levels. The issue is that GRP.default (which handles list/data.frame inputs, i.e. all multivariate cases and most of the named API) never dispatches to GRP.factor — it uses radix sort which naturally only creates groups for observed combinations.

How dplyr handles multivariate drop = FALSE

For group_by(df, f1, f2, .drop = FALSE), dplyr generates the full Cartesian product of all factor levels — including unobserved combinations. So if f1 has levels A,B,C and f2 has levels X,Y, you get 6 groups, regardless of which appear in the data.

The key insight for a collapse-friendly implementation

The Cartesian product approach maps naturally to integer arithmetic over factor codes, which fits well with collapse's integer group.id system:

# For factors f1 (L1 levels) and f2 (L2 levels):
# group.id = (as.integer(f1) - 1L) * L2 + as.integer(f2)
# Total groups = L1 * L2 (unobserved combinations have group.size = 0)

This means no special-casing is needed in the C layer — group.id is still an integer vector of the same length as the data rows, pointing into a potentially larger-than-observed group.sizes vector. The key change is entirely in R.

Proposed approaches (ordered by scope)

---

Option 1 — fcount(drop = FALSE) only (targeted, low complexity)

For fcount, drop = FALSE is a pure post-processing step: compute counts normally, then expand to include all factor level combinations with N = 0.

fcount <- function(x, ..., w = NULL, name = "N", add = FALSE, sort = FALSE, 
                   decreasing = FALSE, drop = TRUE) {
  # ... existing code to build g ...
  res <- fcount_core(x, g, w, name, add)
  if(!drop && !add) res <- .expand_fcount(res, g, name)
  res
}

.expand_fcount <- function(res, g, name) {
  factor_cols <- which(vapply(g$groups, is.factor, TRUE))
  if(!length(factor_cols)) return(res)
  # Build full grid of all factor level combinations
  full_grid <- do.call(expand.grid, c(
    lapply(g$groups[factor_cols], levels),
    list(stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)
  ))
  names(full_grid) <- names(g$groups)[factor_cols]
  # Convert non-factor group vars to character for matching
  # Left-join res into full_grid, fill unmatched N with 0
  # (simple merge or fjoin-based approach)
}

This is self-contained, doesn't touch GRP internals, and covers the most common use case (fcount is the primary request). The tricky part is the join when some grouping variables are non-factors, but for the factor-only case it's clean.

---

Option 2 — GRP.default(drop = FALSE) (general, medium complexity)

Add drop = FALSE to GRP.default. When any factor columns exist in the input:

1. Extract factor levels for each factor column; non-factor columns only have observed levels
2. Compute the Cartesian product of all factor level combinations as the canonical group list
3. Compute group.id for each data row via integer arithmetic on factor codes
4. Set group.sizes[i] = 0 for unobserved group indices

GRP.default <- function(X, ..., drop = TRUE) {
  # ... existing logic for non-factor case (drop = TRUE behavior) ...
  
  # New path for drop = FALSE with factor columns:
  if(!drop && is.list(X)) {
    factor_idx <- which(vapply(unclass(X), is.factor, TRUE))
    if(length(factor_idx)) {
      # All factor cols: use their levels; others: use observed values
      level_sets <- lapply(seq_along(X), function(i) {
        if(i %in% factor_idx) levels(X[[i]]) else unique(X[[i]])
      })
      # Cartesian product group table
      groups <- do.call(expand.grid, c(level_sets, list(stringsAsFactors = FALSE)))
      N.groups <- nrow(groups)
      # group.id via positional arithmetic
      strides <- cumprod(c(1L, vapply(level_sets[-length(level_sets)], length, 1L)))
      codes <- lapply(seq_along(X), function(i) {
        if(i %in% factor_idx) as.integer(X[[i]]) 
        else match(X[[i]], level_sets[[i]])
      })
      group.id <- Reduce("+", Map(function(c, s) (c - 1L) * s, codes, strides)) + 1L
      group.sizes <- tabulate(group.id, N.groups)
      # ... build and return GRP object ...
    }
  }
}

This enables fgroup_by(.drop = FALSE), which would make the feature available across all aggregation functions automatically.

---

Option 3 — Separate complete() utility (non-invasive)

Add a complete_groups() function that post-processes any GRP-based result:

complete_groups <- function(x, fill = list(N = 0L))

This is like tidyr::complete() and doesn't touch existing functions at all, but puts the burden on the user.

---

Recommendation

For minimum viable implementation, I'd suggest Option 1 (fcount only) as a first step — it addresses the most common use case, is self-contained, and the post-processing overhead is negligible compared to the counting step.

For the general case, Option 2 is architecturally cleaner since drop = FALSE on the GRP object propagates automatically to collap, fsummarise, and all other aggregation functions. The integer-arithmetic approach for group.id means no changes to C code are needed — the GRP object produced is structurally identical to what the existing C layer already handles (including group.size = 0 groups, which it already handles for factor inputs via GRP.factor).

The one design question for Option 2: should drop = FALSE silently have no effect when no factors are present, or should it error/warn? dplyr silently does nothing, which seems right.

Sebastian's concern about "univariate exemption"

I think this concern is addressed by the integer-arithmetic approach — it handles both univariate and multivariate factor grouping uniformly, with no special paths based on the number of grouping variables. The logic is: "if any column is a factor and drop = FALSE, expand to Cartesian product of their levels." Non-factor columns are treated as having only their observed levels (same as current behavior).

[egoipse]

Thanks for the answer @SebKrantz

I believe this issue has already been requested in a similar form here

Sure! This is the exact issue I linked in my first post.

What does dplyr do on multivariate grouping data if one variable is a factor?

I think @claude explanation answer your question. The goal is to keep levels of a factor even when the manipulation data output has no data for one or more of them. Exactly what drop = F does in dplyr::group_by

[SebKrantz]

Thanks @egoipse. @claude, if we do this, we need a thorough solution in GRP.default(), GRP.factor() already has a drop = FALSE argument. The solution needs to be implemented largely in C and have minimal overhead. For the multivariate case with the cartesian product, I think something similar to the multi_match() C-level helper for joins with multiple = TRUE should give inspiration:

collapse/src/join.c

Line 329 in 015fdaa

SEXP multi_match(SEXP m, SEXP g) {

. Please draft a plan where we largely use C or existing helper functions implemented in C to add this efficiently.

[claude]

Claude finished @SebKrantz's task in 4m 25s —— View job

---

Implementation Plan: drop = FALSE in GRP.default() via C

• Gather context and explore codebase
• Read multi_match(), GRP.default(), GRP.factor(), fcrosscolon(), fwtabulate()
• Verify aggregation C functions handle zero-count groups (confirmed: no changes needed)
• Draft full implementation plan

---

Key Findings

What already works: GRP.factor already supports drop = FALSE — it uses factor integer codes directly as group.id and fwtabulate to get group.sizes (which includes zeros for absent levels). This is the template.

What doesn't exist yet: GRP.default (the list/data.frame path) never dispatches to GRP.factor. It uses radixsort or hash, both of which only produce observed groups.

Critical insight — no changes needed in aggregation C code: Every aggregation function (e.g. fsum_double_g_impl in src/fsum.c) allocates output using g[[1L]] (N.groups) and writes via g[[2L]][i] (group.id). Slots that receive zero data points keep their initialization value (NA_REAL for na.rm=TRUE, 0.0 otherwise). This is already correct for empty factor levels.

The multi_match inspiration: multi_match (src/join.c:329) builds a positional expansion array using stride arithmetic over cumulative group sizes. The same pattern — computing per-column level indices for the full Cartesian product — is what we need to build the groups table for drop = FALSE.

The fcrosscolon shortcut: src/small_helper.c already has fcrosscolon(x, ng, y, ckna) which does in-place x[i] += (y[i] - 1) * ng with NA propagation. This is precisely the stride multiplication needed to combine per-column group IDs into a single Cartesian group ID. We can reuse this logic directly.

---

Implementation Plan

Step 1 — New C function GRP_factor_cross_C in src/small_helper.c

This is the main new piece of C code. It takes per-column group ID vectors and returns the combined group.id, group.sizes, and per-column level indices for every group in the Cartesian product.

/*
 * GRP_factor_cross_C: Build combined group.id for drop=FALSE grouping.
 *
 * ids      - R list of nc integer vectors, each length n, 1-based group codes
 *            per column (NAs permitted). For factor columns these are as.integer()
 *            codes; for non-factor columns they come from C_group.
 * nlevels  - integer vector of length nc: number of groups per column
 *            (= nlevels for factors, N.groups from C_group for non-factors)
 * ret_grp  - logical: whether to also return the level-index table
 *
 * Returns a list:
 *   [[1]] group.id   : integer[n],   1-based combined group ID (NA for any NA column)
 *   [[2]] group.sizes: integer[Ng],  Ng = product(nlevels); 0 for unobserved combos
 *   [[3]] N.groups   : scalar int,   Ng
 *   [[4]] level.idx  : list of nc integer[Ng] vectors giving per-column 1-based
 *                      level index for each of the Ng Cartesian combinations,
 *                      or R_NilValue if ret_grp = FALSE
 */
SEXP GRP_factor_cross_C(SEXP ids, SEXP nlevels, SEXP ret_grp) {

  int nc = length(ids), n = length(VECTOR_ELT(ids, 0));
  int *pnl = INTEGER(nlevels), retg = asLogical(ret_grp);

  // --- Overflow-safe product of nlevels ---
  long long Ng_ll = 1;
  for(int j = 0; j < nc; ++j) {
    Ng_ll *= pnl[j];
    if(Ng_ll > INT_MAX)
      error("Total number of factor level combinations (%lld) exceeds INT_MAX. "
            "Consider using drop = TRUE.", Ng_ll);
  }
  int Ng = (int)Ng_ll;

  // --- Build combined group.id via stride arithmetic ---
  // group.id[i] = code_0[i] + (code_1[i]-1)*nl_0 + (code_2[i]-1)*nl_0*nl_1 + ...
  // This is identical to what fcrosscolon does, but non-destructive:
  SEXP gid = PROTECT(allocVector(INTSXP, n));
  int *pgid = INTEGER(gid);
  int *p0 = INTEGER(VECTOR_ELT(ids, 0));
  for(int i = 0; i < n; ++i) pgid[i] = p0[i];  // copy first column codes

  int stride = pnl[0];
  for(int j = 1; j < nc; ++j) {
    int *pj = INTEGER(VECTOR_ELT(ids, j));
    for(int i = 0; i < n; ++i) {
      if(pgid[i] == NA_INTEGER) continue;       // already NA: propagate
      if(pj[i]  == NA_INTEGER) { pgid[i] = NA_INTEGER; continue; }
      pgid[i] += (pj[i] - 1) * stride;
    }
    stride *= pnl[j];
    // stride == Ng after the last iteration
  }

  // --- Tabulate group.sizes (like fwtabulate, narm = TRUE) ---
  SEXP gs = PROTECT(allocVector(INTSXP, Ng));
  int *pgs = INTEGER(gs);
  memset(pgs, 0, sizeof(int) * Ng);
  --pgs;  // 1-based offset
  for(int i = 0; i < n; ++i)
    if(pgid[i] != NA_INTEGER) ++pgs[pgid[i]];

  // --- Per-column level indices for every Cartesian combination ---
  // For group g (0-based): col j's level = (g / stride_j) % nl_j + 1  (1-based)
  // This mirrors the multi_match pattern of building a positional expansion:
  // instead of expanding rows, we enumerate all level combinations.
  SEXP lidx = R_NilValue;
  if(retg) {
    PROTECT(lidx = allocVector(VECSXP, nc));
    int stride_j = 1;
    for(int j = 0; j < nc; ++j) {
      SEXP col = PROTECT(allocVector(INTSXP, Ng));
      int *pcol = INTEGER(col);
      for(int g = 0; g < Ng; ++g)
        pcol[g] = (g / stride_j) % pnl[j] + 1;   // 1-based
      SET_VECTOR_ELT(lidx, j, col);
      UNPROTECT(1);
      stride_j *= pnl[j];
    }
  }

  // --- Assemble result ---
  SEXP res = PROTECT(allocVector(VECSXP, 4));
  SET_VECTOR_ELT(res, 0, gid);
  SET_VECTOR_ELT(res, 1, gs);
  SET_VECTOR_ELT(res, 2, ScalarInteger(Ng));
  SET_VECTOR_ELT(res, 3, lidx);
  UNPROTECT(retg ? 4 : 3);
  return res;
}

Register in src/ExportSymbols.c as {"C_GRP_factor_cross", (DL_FUNC) &GRP_factor_cross_C, 3} and add the declaration to src/collapse_c.h.

---

Step 2 — GRP.default() R changes (R/GRP.R)

Add drop = TRUE parameter and a new code path triggered when drop = FALSE and any column is a factor:

GRP.default <- function(X, by = NULL, sort = .op[["sort"]], decreasing = FALSE,
                        na.last = TRUE, return.groups = TRUE, return.order = sort,
                        method = "auto", drop = TRUE, call = TRUE, ...) {

  # ... existing setup (resolving by, namby) unchanged ...

  # New drop=FALSE path for list input with factor columns
  if(!drop && is.list(X)) {
    cols <- if(is.null(by)) unclass(X) else .subset(unclass(X), by)
    factor_idx <- vapply(cols, is.factor, TRUE)

    if(any(factor_idx)) {
      # Per-column group IDs and level sets
      ids <- vector("list", length(cols))
      levsets <- vector("list", length(cols))
      nlevels_vec <- integer(length(cols))

      for(j in seq_along(cols)) {
        col <- cols[[j]]
        if(factor_idx[[j]]) {
          col <- addNA2(col)               # NA -> extra level, matches GRP.factor
          ids[[j]] <- unattrib(col)        # integer codes, 1-based
          levsets[[j]] <- attr(col, "levels")
          nlevels_vec[[j]] <- length(levsets[[j]])
        } else {
          # Non-factor: hash into observed groups only
          g <- .Call(C_group, col, FALSE, FALSE)
          ids[[j]] <- unattrib(g)
          # Observed unique values (in hash encounter order)
          levsets[[j]] <- .Call(C_funique, col)  # or use group.starts to subset
          nlevels_vec[[j]] <- attr(g, "N.groups")
        }
      }

      # Combined group.id, group.sizes, level indices via C
      res <- .Call(C_GRP_factor_cross, ids, nlevels_vec, return.groups)
      gid  <- res[[1L]]
      gs   <- res[[2L]]
      Ng   <- res[[3L]]
      lidx <- res[[4L]]   # list of nc int vectors, or NULL

      # Build groups data.frame from level indices
      groups <- NULL
      if(return.groups) {
        groups <- vector("list", length(cols))
        names(groups) <- namby
        for(j in seq_along(cols)) {
          lev <- levsets[[j]]
          # level index vector (1-based) -> subset actual level values
          groups[[j]] <- lev[lidx[[j]]]
          if(factor_idx[[j]]) {
            # Re-attach factor class so print/GRPnames works
            attr(groups[[j]], "levels") <- lev
            oldClass(groups[[j]]) <- "factor"
          }
        }
        # Wrap as data.frame-like list (what GRP.default normally returns via subsetDT)
        class(groups) <- "data.frame"
        attr(groups, "row.names") <- .set_row_names(Ng)
      }

      return(`oldClass<-`(list(
        N.groups   = Ng,
        group.id   = gid,
        group.sizes = gs,
        groups     = groups,
        group.vars = namby,
        ordered    = c(ordered = NA, sorted = NA),  # hash-like: no sort guarantee
        order      = NULL,
        group.starts = NULL,   # consistent with GRP.factor for drop=FALSE
        call       = if(call) match.call() else NULL
      ), "GRP"))
    }
    # fall through to existing path when no factors present (drop=FALSE is a no-op)
  }

  # ... existing radix/hash path unchanged ...
}

---

Step 3 — fgroup_by() change (R/GRP.R)

fgroup_by <- function(.X, ..., sort = .op[["sort"]], decreasing = FALSE, na.last = TRUE,
                      return.groups = TRUE, return.order = sort, method = "auto",
                      .drop = TRUE) {
  # ... existing body ...
  attr(.X, "groups") <- GRP.default(e, NULL, sort, decreasing, na.last,
                                    return.groups, return.order, method,
                                    drop = .drop, FALSE)
  # ...
}

---

Why this is efficient

Concern	How addressed
Combined group.id	Single O(n·nc) loop with integer stride arithmetic — identical cost to fcrosscolon
group.sizes	Single O(n) tabulation pass, no allocation beyond the output vector
Groups table	O(Ng·nc) positional index computation (Ng ≪ n in practice); mirrors multi_match's expansion loop
Aggregation C code	Zero changes — fsum_double_g_impl etc. already handle zero-count groups via output-vector initialization
Existing callers	Zero overhead — drop = TRUE default preserves current codepath entirely

---

What the multi_match analogy teaches us

multi_match builds per-element positional mappings using cumulative stride offsets (cgs[k] + cnt[k]++). The level index loop in GRP_factor_cross_C does the same in reverse: instead of expanding match positions, it decomposes a flat group index into per-column level positions via (g / stride_j) % nl_j. Both are O(Ng·nc) positional enumeration loops — no sorting, no hashing.

---

Edge cases

• No factor columns + drop=FALSE: falls through to existing radix/hash path silently (dplyr behavior)
• NA levels: handled identically to GRP.factor via addNA2 — NA becomes a valid extra level
• Overflow: product of all nlevels is checked against INT_MAX before allocation
• sort parameter: set ordered = c(ordered=NA, sorted=NA) for the drop=FALSE case (group IDs are determined by factor code order, not radixsort)
• group.starts: set to NULL (consistent with GRP.factor)

[SebKrantz]

Thanks. This provides good inspiration, and I will use it to craft something @egoipse. I may take some time to do it, though. The implementation here works but does not satisfy my demands in terms of parsimony and performance for internal collapse code.

[egoipse]

I'm pretty impress. Thanks for the disposition.

There's no rush on my part. Mine was just an enhancement suggestion.

But it's impossible not to appreciate your answer to a suggestion. That's why collapse improves every time and is a great asset for R and free software community.

Thanks again @SebKrantz