Combinatorics (was: Merging lists has made my brain hurt.)

[Thorsten Kampe]

* Alex Martelli
> Maybe it's an artefact of having worked too much in combinatorics
> (where double-counting is a perennial worry:-), [...]

I would like you to have a look at my "combinatorics" program. It's a
straight translation from the original RPL (mixture between Lisp and
Forth) code. RPL has no dictionaries just lists so the data
representation could be ineffective.

The code is loosely commented so here's what it does: it deals with
the whole field of combinatorics, that is permutation, variation and
combination with and without repetition.

The resulting list can be quite large, so for convenience the main
"if-branch" computes the length of the resulting list rather than the
list itself. If you're just interested in the count, you can always
give the actual list /or/ the size of the list as argument except for
permutation with repetition where you have to give the actual list.

'P' means 'permutation', 'C' combination and 'V' variation
'#' means return the /count/ otherwise the actual list
'R' means 'with repetition' otherwise it's without
combinatoric([1, 2, 3, 4], '#P')        -> 24  # permutation without repetition
combinatoric(4, '#P')                   -> 24
combinatoric([1, 2, 3, 3], '#PR')       -> 12  # permutation with repetition
combinatoric([1, 2, 3, 4, 5], 2, '#C')  -> 10  # combination without repetition
combinatoric(5, 2, '#C')                -> 10
combinatoric([1, 2, 3, 4, 5], 2, '#CR') -> 15  # combination with repetition
combinatoric(5, 2, '#CR')               -> 15
combinatoric([1, 2, 3, 4, 5], 2, '#V')  -> 20  # variation without repetition
combinatoric(5, 2, '#V')  -> 20
combinatoric([1, 2, 3, 4, 5], 2, '#VR') -> 25  # variation with repetition
combinatoric(5, 2, '#VR')               -> 25

combinatoric([1, 2, 3, 4], 2, 'C')  -> [[2, 3], [3, 4], [1, 3], [2, 4], [1, 4], [1, 2]]
combinatoric([1, 2, 3, 3], 2, 'VR') -> [[1, 1], [1, 2], [1, 3], [1, 3], [2, 1], [2, 2],
                                        [2, 3], [2, 3], [3, 1], [3, 2], [3, 3], [3, 3],
                                        [3, 1], [3, 2], [3, 3], [3, 3]]

                                        
One comment on the permutation algorithm: GvR wrote a nice
demonstration (http://www.python.org/doc/current/ref/indentation.html)
that gives the permutations lexicographically sorted but is about 8
times slower for range(9) than my version.

I treat the argument list strictly as a tuple (where repetition
matters) so for combination and variation I'm actually generating the
lists for the index set range(len(seq)) and finally convert this
back to the original objects.

The program uses some standard external routines:
comb, fact, perm, quotient_set, set and cartes listed beneath the code
for combinatoric.

Thorsten

#v+
def combinatoric(seq, k, modus=''):
    # combination: order is irrelevant, variation: order is relevant
    #from function import comb, fact, perm
    if '#' in modus or '#' in str(k):
        if isinstance(seq, list) and k != '#PR':
            seq = len(seq)
        if modus == '#C':     # combination without repetition
            return comb(seq, k)
        elif modus == '#CR':  # combination with repetition
            return comb(seq+k-1, k)
        elif k == '#P':       # permutation without repetition
            return fact(seq)
        elif k == '#PR':      # permutation with repetition
            return reduce(lambda x, y: x / fact(len(y)),
                          [fact(len(seq))] + quotient_set(seq, lambda x: x).values())
        elif modus == '#V':   # variation without repetition
            return perm(seq, k)
        elif modus == '#VR':  # variation with repetition
            return seq ** k
    else:
        if k == 'P':     # permutation without repetition
            if len(seq) <= 1:
                return [seq]
            else:
                perms = []
                for permutation in combinatoric(seq[:-1], 'P'):
                    for i in range(len(seq)):
                        perms.append(permutation[:])
                        perms[-1].insert(i, seq[-1])
                return perms
        elif k == 'PR':  # permutation with repetition
            return set(combinatoric(seq, 'P'))
        else:            # combination/variation with/without repetition
            if k == 1:
                return [[item] for item in seq]
            else:
                # auxiliary functions
                def setord(seq):  # remove sets with duplicates (order is relevant)
                    def sort(seq):
                        seq.sort()
                        return seq
                    return set([sort(item) for item in seq])

                def setrep(seq):  # remove sets with duplicates (repetition is relevant)
                    def delrep(seq):  # remove duplicates while maintaining order
                        result = []
                        for item in seq:
                            if item not in result:
                                result.append(item)
                        return result
                    return [item for item in seq if item == delrep(item)]


                cartesmodus = 'pair'
                # treat 'seq' as a tuple not as a set; use its index set instead
                result = range(len(seq))
                for i in range(k-1):
                     cartesprod = cartes(result, range(len(seq)), cartesmodus)
                     # filter the cartesian product by order or repetition
                     if modus == 'C':     # combination without repetition
                         result = setrep(setord(cartesprod))
                     elif modus == 'CR':  # combination with repetition
                         result = setord(cartesprod)
                     elif modus == 'V':   # variation without repetition
                         result = setrep(cartesprod)
                     elif modus == 'VR':  # variation with repetition
                         result = cartesprod
                     cartesmodus = 'triple'
                # convert the index sets back to the actual objects
                return [[seq[index] for index in indices] for indices in result]
#v-

#v+
def comb(n, k):
    return fact(n) / (fact(k) * fact(n-k))

def fact(integer):
    factorial = 1
    for factor in range(2, integer+1):
        factorial *= factor
    return factorial

def perm(n, k):
    return fact(n) / fact(n-k)

def cartes(seq0, seq1, modus='pair'):
    """ return the Cartesian Product of two sequences """
    if  modus == 'pair':
        return [[item0, item1] for item0 in seq0 for item1 in seq1]
    elif modus == 'triple':
        return [item0 + [item1] for item0 in seq0 for item1 in seq1]    

def set(seq):
    """ make seq a true set by removing duplicates """
    return [item[0] for item in quotient_set(seq, lambda x: x).values()]

def quotient_set(seq, func):
    """ partition seq into equivalence classes """
    quotient_set = {}
    for item in seq:
        quotient_set.setdefault(repr(func(item)),[]).append(item)
    return quotient_set
#v-