Find the partitions labels for placing n distinct objects into k non-distinct boxes.

stirling.R

################################################################################
# file:   stirling.R
# author: Peter DeWitt
# date:   27 Nov 2013
#
# purpose: 
#
# 1) stirling2(n, k) returns the Stirling number of the second kind.  This is
# the number of ways to place n descript objects into k non-nondescript boxes,
# non-empty boxes.
#
# 2) bell number: B_n = sum_{k=0}^{n} stirling2(n, k)
#
# 3) stirlingGroups(n) returns a matrix with the group labels for enumerating
# the bell(n) partitions.
################################################################################

stirling2 <- function(n, k) { 
  if (k > n ) return(0)
  idx <- 0:k
  return(1 / factorial(k) * sum((-1)^(k-idx) * choose(k, idx) * idx^n))
}

bell <- function(n) { 
  out <- lapply(0:n, stirling2, n = n)
  out <- sum(do.call(c, out))
  return(out)
}

# test 
stirling2(4L, 2L)
stirling2(10, 4)

# test 
bell(5)
bell(6)

# define the next partition.
#
# start with k <- rep(0, n), the label vector
# and m <- k, the max value for the column
#
# Algorithm is thanks to 
# Orlov, Michael. "Efficient generation of set partitions." Engineering and
# Computer Sciences, University of Ulm, Tech. Rep (2002).
partitions <- function(k, m) { 
  n <- length(k) 
  i <- n

  while(i >= 2) { 
    if (k[i] <= m[i-1]) {
      k[i] <- k[i] + 1
      m[i] <- max(m[i], k[i])

      if (i < n) {
        k[(i+1):n] <- k[1]
        m[(i+1):n] <- m[i]
      }
      break
    }
    i <- i - 1
  } 
  return(list(k = k, m = m))
}

# test the partitions function for finding the bell(4) partitions for 4 distinct
# objects into non-empty boxes.
partitions(rep(0, 4), rep(0, 4))
partitions(c(0, 0, 0, 1), c(0, 0, 0, 1))
partitions(c(0, 0, 1, 0), c(0, 0, 1, 1))
partitions(c(0, 0, 1, 1), c(0, 0, 1, 1))
partitions(c(0, 0, 1, 2), c(0, 0, 1, 2))
partitions(c(0, 1, 0, 0), c(0, 1, 1, 1))
partitions(c(0, 1, 0, 1), c(0, 1, 1, 1))
partitions(c(0, 1, 0, 2), c(0, 1, 1, 2))
partitions(c(0, 1, 1, 0), c(0, 1, 1, 1))
partitions(c(0, 1, 1, 1), c(0, 1, 1, 1))
partitions(c(0, 1, 1, 2), c(0, 1, 1, 2))
partitions(c(0, 1, 2, 0), c(0, 1, 2, 2))
partitions(c(0, 1, 2, 1), c(0, 1, 2, 2))
partitions(c(0, 1, 2, 2), c(0, 1, 2, 2))
partitions(c(0, 1, 2, 3), c(0, 1, 2, 2))

# stringGroups, provide the number of objects.  Due to the possible very large
# size of bell(n), save the results to a file.
stirlingGroups <- function(n, file = NULL) { 
  if (is.null(file)) stop("Specify the file to store the partitions in.")

  out <- list(k = rep(0, n), m = rep(0, n))
  write(out$k, ncolumns = n, file = file, append = FALSE)

  while(!all(out$k == 0:(n-1))) { 
    out <- partitions(out$k, out$m)
    write(out$k, ncolumns = n, file = file, append = TRUE)
  }
}

# test
stirlingGroups(7, file = "~/Desktop/temp.txt")
bell(7)

################################################################################
# end of file
################################################################################