Fast method of getting all the descendants of a parent

Question

With the parent-child relationships data frame as below:

  parent_id child_id
1         1        2
2         2        3
3         3        4

The goal is to achieve the following, i.e. expanded version of previous data frame where all the descendants (children, grandchildren, etc.) are assigned to each parent (and incl. the parent/child itself):

   parent_id child_id
1          1        1
2          1        2
3          1        3
4          1        4
5          2        2
6          2        3
7          2        4
8          3        3
9          3        4
10         4        4

The question I have: What is the fastest possible way (or one of them) of achieving that in R?

I've already tried various methods - from a for loop, SQL recursion to using igraph (as described here). They are all rather slow, and some of them are also prone to crashing when dealing with a larger number of combinations.

Below are examples with sqldf and igraph, benchmarked on a slightly larger data frame than above.

library(sqldf)
library(purrr)
library(dplyr)
library(igraph)

df <- data.frame(parent_id = 1:1000L)
df$child_id <- df$parent_id + 1L

# SQL recursion

sqlQuery <- 'with recursive
             dfDescendants (parent_id, child_id)
             as
             (select parent_id, child_id from df
             union all
             select d.parent_id, s.child_id from dfDescendants d
             join df s
             on d.child_id = s.parent_id)
             select distinct parent_id, parent_id as child_id from dfDescendants
             union
             select distinct child_id as parent_id, child_id from dfDescendants
             union
             select * from dfDescendants;'

sqldf(sqlQuery)

# igraph with purrr

df_g = graph_from_data_frame(df, directed = TRUE)

map(V(df_g), ~ names(subcomponent(df_g, .x, mode = "out"))) %>% 
  map_df(~ data.frame(child_id = .x), .id = "parent_id")

Benchmark (excl. query creation in sqldf and conversion to graph in igraph):

set.seed(23423)

microbenchmark::microbenchmark(
  sqldf = sqldf(sqlQuery),
  tidyigraph = map(V(df_g), ~ names(subcomponent(df_g, .x, mode = "out"))) %>% 
    map_df(~ data.frame(child_id = .x), .id = "parent_id"),
  times = 5
)

#    Unit: seconds
#           expr      min       lq     mean   median       uq      max neval
#          sqldf 7.815179 8.002836 8.113392 8.084038 8.315207 8.349701     5
#     tidyigraph 5.784239 5.806539 5.883241 5.889171 5.964906 5.971350     5

Answer 1

We can use ego like below

library(igraph)
df <- data.frame(parent_id = 1:3, child_id = 2:4)
g <- graph_from_data_frame(df)

setNames(
  rev(
    stack(
      Map(
        names,
        setNames(
          ego(g,
            order = vcount(g),
            mode = "out"
          ),
          names(V(g))
        )
      )
    )
  ),
  names(df)
)

which gives

   parent_id child_id
1          1        1
2          1        2
3          1        3
4          1        4
5          2        2
6          2        3
7          2        4
8          3        3
9          3        4
10         4        4

---

Benchmarking

set.seed(23423)

microbenchmark::microbenchmark(
  sqldf = sqldf(sqlQuery),
  tidyigraph = map(V(df_g), ~ names(subcomponent(df_g, .x, mode = "out"))) %>%
    map_df(~ data.frame(child_id = .x), .id = "parent_id"),
  ego = setNames(
    rev(
      stack(
        Map(
          names,
          setNames(
            ego(df_g,
              order = vcount(df_g),
              mode = "out"
            ),
            names(V(df_g))
          )
        )
      )
    ),
    names(df)
  ),
  times = 5
)

shows

Unit: milliseconds
       expr       min       lq      mean    median         uq        max neval
      sqldf 7156.2753 9072.155 9402.6904 9518.2796 10206.3683 11060.3738     5
 tidyigraph 2483.9943 2623.558 3136.7490 2689.8388  2879.5688  5006.7853     5
        ego  182.5941  219.151  307.2481  253.2171   325.8721   555.4064     5

---

The code using pipes for readability :

g |>
  ego(order = vcount(g), mode = "out") |>
  setNames(names(V(g))) |>
  Map(f = names) |>
  stack() |>
  rev() |>
  setNames(names(df))

Answer 2

igraph is of course a good way to answer graph questions (see also Bioconductor's graph + RBGL packages), but I think this has an iterative solution in R.

It seems like a reasonable approach is to perform a depth-first graph traversal (I was expecting a fancier solution). This is actually quite easy to implement efficiently in R. Suppose vectors pid and cid describe the links between parent and child nodes in the graph (as in the data.frame in the question). Represent each node as a positive integer.

all_nodes <- unique(c(parent_id, child_id)  # all nodes
uid <- match(all_nodes, all_nodes)
pid <- match(parent_id, all_nodes)
cid <- match(child_id, all_nodes)

and form an edge list from each node to all its children.

edge_list <- unname(split(cid, factor(pid, levels = uid)))
edge_lengths <- lengths(edge_list)

The child nodes of the current child nodes is edge_list[cid], and the number of level 2 child nodes associated with each original parent node is rep(pid, edge_lengths[cid]). So the path from any node to any other reachable node is traversed by the simple iteration

while (length(pid)) {
    pid <- rep(pid, edge_lengths[cid])
    cid <- unlist(edge_list[cid])
}

@jblood94 points out that the traverse has to keep track of edges already visited. We can implement this efficiently (in time, not space!) by creating a logical vector of visited edges. We use a 'factory' pattern, where we create a function that retains state (the logical vector of nodes visited). The vector is indexed by the unique id (key) of the edge pid * n + cid. We're interested in keys that are not duplicated, and have not been already visited.

visitor <- function(uid, n_max = 3000) {
    n <- length(uid)
    if (n <= n_max) {
        ## over-allocate, to support `key = pid * n + cid`
        visited <- logical((n + 1L) * n) # FALSE on construction
    } else {
        stop("length(uid) greater than n_max = ", n_max)
    }
    function(pid, cid) {
        key <- pid * n + cid
        to_visit <- !(duplicated(key) | visited[key])
        visited[key[to_visit]] <<- TRUE  # update nodes that we will now visit
        to_visit
    }
}

Thus

> visit = visitor(1:10)
> visit(1:3, 2:4)
[1] TRUE TRUE TRUE
> visit(2:4, 3:5)
[1] FALSE FALSE  TRUE

Here is a more complete implementation of the entire solution, with additional book-keeping

visitor <- function(uid, n_max = 3000) {
    n <- length(uid)
    if (n <= n_max) {
        ## over-allocate, to support `key = pid * n + cid`
        visited <- logical((n + 1L) * n) # FALSE on construction
    } else {
        stop("length(uid) greater than n_max = ", n_max)
    }
    function(pid, cid) {
        key <- pid * n + cid
        to_visit <- !(duplicated(key) | visited[key])
        visited[key[to_visit]] <<- TRUE
        to_visit
    }
}

ancestor_descendant <- function(df) {
    ## encode parent and child to unique integer values
    ids <- unique(c(df$parent_id, df$child_id))
    uid <- match(ids, ids)
    pid <- match(df$parent_id, ids)
    cid <- match(df$child_id, ids)
    n <- length(uid)

    ## edge list of parent-offspring relationships, based on unique
    ## integer values; list is ordered by id, all ids are present, ids
    ## without children have zero-length elements. Use `unname()` so
    ## that edge_list is always indexed by integer
    edge_list <- unname(split(cid, factor(pid, levels = uid), drop = FALSE))
    edge_lengths <- lengths(edge_list)

    visit <- visitor(uid)
    keep <- visit(uid, uid) # all TRUE
    aid = did = list(uid) # results -- all uid's are there own ancestor / descendant
    i = 1L
   
    while (length(pid)) {
        ## only add new edges
        keep <- visit(pid, cid)
        ## record current generation ancestors and descendants
        pid <- pid[keep]
        cid <- cid[keep]
        i <- i + 1L
        aid[[i]] <- pid
        did[[i]] <- cid

        ## calculate next generation pid and cid.
        pid <- rep(pid, edge_lengths[cid])
        cid <- unlist(edge_list[cid])
    }
    ## decode results to original ids and clean up return value
    df <- data.frame(
        ancestor_id = ids[unlist(aid)],
        descendant_id = ids[unlist(did)]
    )
    df <- df[order(df$ancestor_id, df$descendant_id),]
    rownames(df) <- NULL
    df
}

This seems to be correct and performant, at least superficially

## Original example
df <- data.frame(parent_id = 1:1000L)
df$child_id <- df$parent_id + 1L
df = df[sample(nrow(df)),]
system.time(result <- ancestor_descendant(df))
##  user  system elapsed 
## 0.243   0.001   0.245 
dim(result)
## [1] 501501      2

## updated example from comments
df <- data.frame(parent_id = 1:1000L)
df$child_id <- df$parent_id + 1L
df <- rbind(df, data.frame(parent_id = 1000L, child_id = 1002L))
system.time(result <- ancestor_descendant(df))
##  user  system elapsed 
## 0.195   0.001   0.195 
dim(result)
## [1] 502502      2

## problematic case from @jblood94
df <- data.frame(
    parent_id=c(1, 1, 2),
    child_id = c(2, 3, 3)
)
ancestor_descendant(df)
##   ancestor_id descendant_id
## 1           1             1
## 2           1             2
## 3           1             3
## 4           2             2
## 5           2             3
## 6           3             3

## previously failed without filtering re-visited nodes
df <- data.frame(
    parent_id = rep(1:100, each = 2),
    child_id = c(2, rep(3:101, each = 2), 102)
)
system.time(result <- ancestor_descendant(df))
##  user  system elapsed 
## 0.005   0.000   0.006 
dim(result)
## [1] 5252    2