[Numpy-discussion] Extracting all the possible combinations of a grid

[Gael Varoquaux]

On Sat, Sep 22, 2007 at 10:35:16AM +0200, Gael Varoquaux wrote:
> I would go for the "generate_fourplets" solution if I had a way to
> calculate the binomial coefficient without overflows.

Sorry, premature optimisation is the root of all evil, but turning ones
brain on early is good.

"""
##############################################################################
# Some routines for calculation of binomial coefficients
def gcd(m,n):
    while n:
        m,n=n,m%n
    return m
    
def binom_(n,k):
        if k==0:
            return 1
        else:
            g = gcd(n,k)
            return binomial(n-1, k-1)/(k/g)*(n/g)
        
def binomial(n,k): 
    if k > n/2: # Limit recursion depth
        k=n-k 
    return binom_(n,k) 
"""

This is surprisingly fast (surprising for me, at least).

Using this and the C code I have, I can generate the quadruplets of 200
integers quite quickly:

In [5]: %timeit b = [1 for i in  generate_quadruplets(200)]
10 loops, best of 3: 1.61 s per loop

With generate_quadruplets given by:

"""
##############################################################################
def generate_quadruplets(size):
    """ Returns an iterator on tables listing all the possible unique 
    combinations of four integers below size. """
    
    C_code = """
    int index = 0;
    for (int j=0; j<i+1; j++) {
        for (int k=0; k<j+1; k++) {
            for (int l=0; l<k+1; l++) {
                quadruplets(index, 0) = i;
                quadruplets(index, 1) = j;
                quadruplets(index, 2) = k;
                quadruplets(index, 3) = l;
                index++ ;
            }
        }
    }
    """
   
    for i in xrange(size):
        multiset_coef = binomial(i+3, 3)
        quadruplets = empty((multiset_coef, 4), dtype=uint32)
        inline(C_code, ['quadruplets', 'i'],
				type_converters=converters.blitz)

        yield quadruplets
"""

This fits my needs.

Maybe I should add this on the cookbook. It is a bit specific, but it
shows a serie of interesting problems, I think.

I also think the binomial and gcd should go in scipy (I could not find
them, but maybe they are already there). Maybe in a new module, as I
don't really see where this fits in.

Gaël