Listing all integer compositions in R

Viewed 130

What I look for is basically an R-version of the answer to this question: Generating all permutation of numbers that sums up to N. First of all the answer uses java, which I have a really hard time reading. Second of all the code uses "deque", which I cant figure out a way to implement in R. I have found several algorithms to do this, but they have all been written in programming languages using structures not available in R such as deques, heaps or list-comprehensions.

What I actually need is a way of finding all the vectors v of length N-1 where: sum(v * 1:(N-1)) == N and I think I can manage that myself if only I find a way of obtaining all the ordered integer partitions.

As an example for N = 4 all the ordered integer partitions using numbers 1 to N-1 are:

1+1+1+1

1+1+2

1+3

2+2

What I effectively need is output of the either form:

c(1,1,1,1)

c(1,1,2)

c(1,3)

c(2,2)

Or of the form:

c(4,0,0)

c(2,1,0)

c(1,0,1)

c(0,2,0)

since I should be able to convert the former format to the latter by myself. Any hint as to how to approach this problem using R would be greatly appreciated. The latter format is excactly the vectors v such that sum(v * 1:3) is 4.

EDIT: My own attempt:

rek = function(mat, id1, id2){
  if(id1 + id2 != length(mat) + 1){ #If next state not absorbing
    mat[id1] = mat[id1] - 1
    mat[id2] = mat[id2] - 1
    mat[id1+id2] = mat[id1+id2] + 1
    out = mat
    
    id = which(mat > 0)
    for(i in id){
      for(j in id[id>=i]){
        if(j == i & mat[i] == 1){
          next
        }
        out = rbind(out, rek(mat,i,j))
      }
    }
    return(out)
  }
}
start = c(n, rep(0, n-2))
states = rbind(start, rek(start, 1, 1))
states = states[!duplicated(states), ] #only unique states.

This is incredibly inefficient. E. g. when n = 11, my states has over 120,000 rows prior to removing duplicates, which leaves only 55 rows.

EDIT 2: Using the parts() function described below I came up with:

temp = partitions::parts(n)
temp = t(temp)
for(i in 1:length(temp[,1])){
  row = temp[i,]
  if(any(row>(n-1))){#if absorbing state
    next
  }
  counts = plyr::count(row[row>0])
  newrow = rep(0,n-1)
  id = counts$x
  numbs = counts$freq
  newrow[id] = numbs
  states = rbind(states, newrow)
}
states = states[-1,]#removing the first row, added manually

which excactly gives me the vectors v such that sum(v * 1:(N-1)) is N.

If anyone is interested, this is to be used within coalescent theory, as a way to describe the possible relations between N individuals omitting when all are related. As an example with N = 4:

(4, 0, 0) -- No individuals are related

(2, 1, 0) -- Two individuals are related, the rest are not

(0, 2, 0) -- The individuals are pair-wise related

(1, 0, 1) -- Three individuals are related, the other individual is not.

1 Answers

Hope parts from package partitions could help

library(partitions)
N <- 4
res <- unique(lapply(asplit(parts(N),2),function(x) sort(x[x>0])))[-1]

which gives

> res
[[1]]
[1] 1 3

[[2]]
[1] 2 2

[[3]]
[1] 1 1 2

[[4]]
[1] 1 1 1 1

If you would like to write a custom base R function, here is a recursive version

f <- function(n, vhead = n, v = c()) {
  if (n == 0) return(list(v))
  unlist(lapply(seq_len(min(n, vhead)), function(k) f(n - k, k, c(k,v))), recursive = FALSE)
}

then we can run

res <- Filter(function(x) length(x)>1,f(N))
Related