[R] Slow enumeration of vectors with restrictions

[Bryan Keller]

I have a vector T and a scalar m such that 1 <= m <= sum(T).  I also have a vector U that satisfies three properties...
 1. length(T)==length(U)
 2. sum(U)==m
 3. for (i in 1:length(T)) U[i]<=T[i]

The function "nextu" (given below) calculates the "next" vector subject to the same restrictions.  The recursion A(T,m) is given below.  It calculates the total number of different vectors there are for given T and m (in the example here it is 20). 

For example, if 

T <- c(2,4,3,1)
m <- 4

#first I generate the "first" value of U 
#to be plugged into function "nextu"

U.first <- numeric(c)
i<-1
while(sum(U.first) < m) {
if(T[i] > U.first[i]) {
U.first[i]<-(U.first[i]+1)
} else {
i<-(i+1)
}
}

#then I create a place to store the 
#vectors I will create with "nextu"

np <- A(T,m)
c <- length(T)
mat <- matrix(numeric(np*c),np,c)

#then I run nextu starting with U.first
#and store results in mat

i <- 1
U <- U.first
while (i <= np) {
	for (j in 1:c) {
		mat[i,j] <- U[j]
	}
	U <- nextu(U)
	i <- i+1
}


In the example I gave above we get the following matrix...

      [,1] [,2] [,3] [,4]
 [1,]    2    2    0    0
 [2,]    1    3    0    0
 [3,]    0    4    0    0
 [4,]    2    1    1    0
 [5,]    1    2    1    0
 [6,]    0    3    1    0
 [7,]    2    0    2    0
 [8,]    1    1    2    0
 [9,]    0    2    2    0
[10,]    1    0    3    0
[11,]    0    1    3    0
[12,]    2    1    0    1
[13,]    1    2    0    1
[14,]    0    3    0    1
[15,]    2    0    1    1
[16,]    1    1    1    1
[17,]    0    2    1    1
[18,]    1    0    2    1
[19,]    0    1    2    1
[20,]    0    0    3    1

Although this works perfectly, it is way too slow to be useful in a program I am writing.  I’m wondering if anyone has any ideas about how to speed up the process.

If it helps, one can think of reversing the vectors and seeing them as numbers where addition is carried out in base T[i]+1 using carries from the (i-1) coordinate.

#function "nextu" creates the next
#sequential ordering of U that satisfies
#the restrictions given by T and m

nextu <- function(U) {
s <- -1
i <- 1

while (i <= c) {
	if (U[i] == 0) {
		i <- i+1
		} else {
	s <- s+U[i]
	i <- i+1
	break
	}
}

while (i <=c) {
	if (U[i] == T[i]) {
		s <- s+U[i]
		i <- i+1
		} else {
	k <<- i
	U[k] <- (U[k] + 1)
	i <- 1
	break	
	}
}

U[i] <- min(s,T[i])
s <- s-U[i]
i <- i+1

while (i < k) {
	U[i] <- min(s,T[i])
	s <- s-U[i]
	i <- i+1
	}
U
}
#end of function "nextu"

#function "A" implements a recursion to 
#get the total number of nextu calculations
#needed
#Thanks to Martin Morgan and Bill Dunlap for
#help speeding A(T,m) up more than 1000 times!

A <- function(T, m) {
    C <- function(lt, m) {
        if (lt == 1L) {
            R <- as.numeric(0 <= m & m <= T[1])
        } else if (lt == 2L) {
            mu <- m - (0:T[2L])
            R <- sum(mu <= T[1L] & 0L <= mu)
        } else {
            R <- 0
            lt1 <- lt-1L
            for (mu in m-(0:T[lt])) {
                if (is.na(memo[lt1, offset + mu]))
                    memo[lt1, offset + mu] <<- C(lt1, mu)
                R <- R + memo[lt1, offset + mu]
            }
        }
        R
    }
    T <- rev(sort(T))
    m <- as.integer(m)
    offset <- sum(T[-1]) - m + 1L
    nrow <- length(T) - 1L
    memo <- matrix(rep(NA_real_, nrow * (offset + m)), nrow=nrow)
    C(length(T), m)
}
#end of function A



-------------
Bryan Keller, Doctoral Student/Project Assistant
Educational Psychology - Quantitative Methods
The University of Wisconsin - Madison