A permutation journey

Rosetta Code had a new task that had this expression for calculating the determinant of a square matrix (Leibniz's formula):

 

 You can read the wikipedia article to find out more about sigma and sgn or this article which took me step-by-step through to understanding it.

Although the Leibniz formula is inefficient for calculation in general I decided to use it, which meant I needed to create permutations of n items,  σ, and the sign of each permutation

Swapping signs

The sign of a permutation is defined as the parity of the number of swaps between any two items of a permutation needed to generate that permutation from the initial permutation.
So for example, given a starting perm of [0,1,2] the permutation [2,0,1] can be achieved by first swapping the 2 and the 0 to form [2,1,0], then swapping the 1 and 0 to form [2,0,1. That is two swaps in total. The sign of the permutation is given as -1 if an odd number of swaps are used in generating it otherwise as +1 if an even (or no) swaps are needed. (The sign of the permutation will not change even if other swaps, or redundant swaps are used to generate it from the same starting perm).

Pythons itertools.permutations function was of unknown sign, but I picked up on the  Steinhaus–Johnson–Trotter algorithm that generates permutations in such a way that it guarantees that successive permutations only differ by the swapping of two items - this meant that the sign would start as +1 and alternate in sign between successive items.

The Steinhaus–Johnson–Trotter Python program

I decided to create a new Rosetta code task for generating both permutations and their sign with a view for their subsequent use in calculating determinants:

from operator import itemgetter
 
DEBUG = False # like the built-in __debug__
 
def spermutations(n):
    """permutations by swapping. Yields: perm, sign"""
    sign = 1
    p = [[i, 0 if i == 0 else -1] # [num, direction]
         for i in range(n)]
 
    if DEBUG: print ' #', p
    yield tuple(pp[0] for pp in p), sign
 
    while any(pp[1] for pp in p): # moving
        i1, (n1, d1) = max(((i, pp) for i, pp in enumerate(p) if pp[1]),
                           key=itemgetter(1))
        sign *= -1
        if d1 == -1:
            # Swap down
            i2 = i1 - 1
            p[i1], p[i2] = p[i2], p[i1]
            # If this causes the chosen element to reach the First or last
            # position within the permutation, or if the next element in the
            # same direction is larger than the chosen element:
            if i2 == 0 or p[i2 - 1][0] > n1:
                # The direction of the chosen element is set to zero
                p[i2][1] = 0
        elif d1 == 1:
            # Swap up
            i2 = i1 + 1
            p[i1], p[i2] = p[i2], p[i1]
            # If this causes the chosen element to reach the first or Last
            # position within the permutation, or if the next element in the
            # same direction is larger than the chosen element:
            if i2 == n - 1 or p[i2 + 1][0] > n1:
                # The direction of the chosen element is set to zero
                p[i2][1] = 0
        if DEBUG: print ' #', p
        yield tuple(pp[0] for pp in p), sign
 
        for i3, pp in enumerate(p):
            n3, d3 = pp
            if n3 > n1:
                pp[1] = 1 if i3 < i2 else -1
                if DEBUG: print ' # Set Moving'
 
 
if __name__ == '__main__':
    from itertools import permutations
 
    for n in (3, 4):
        print '\nPermutations and sign of %i items' % n
        sp = set()
        for i in spermutations(n):
            sp.add(i[0])
            print('Perm: %r Sign: %2i' % i)
            #if DEBUG: raw_input('?')
        # Test
        p = set(permutations(range(n)))
        assert sp == p, 'Two methods of generating permutations do not agree'

Determinants using Leibniz's formula

I then had all I needed to create the determinant (and permanent) of a matrix code:

from itertools import permutations
from operator import mul
from math import fsum
from spermutations import spermutations
 
def prod(lst):
    return reduce(mul, lst, 1)
 
def perm(a):
    n = len(a); 
    r = range(n); 
    s = list(permutations(r))
    return fsum(prod(a[i][sigma[i]] for i in r) for sigma in s)
 
def det(a):
    n = len(a); 
    r = range(n); 
    s = list(spermutations(n))
    return fsum(sign*prod(a[i][sigma[i]] for i in r) for sigma, sign in s)
 
 
if __name__ == '__main__':
    from pprint import pprint as pp
 
    for a in ( 
            [
             [1, 2], 
             [3, 4]], 
 
            [
             [1, 2, 3, 4],
             [4, 5, 6, 7],
             [7, 8, 9, 10],
             [10, 11, 12, 13]],        
 
            [
             [ 0,  1,  2,  3,  4],
             [ 5,  6,  7,  8,  9],
             [10, 11, 12, 13, 14],
             [15, 16, 17, 18, 19],
             [20, 21, 22, 23, 24]],
        ):
        print('')
        pp(a)
        print('Perm: %s Det: %s' % (perm(a), det(a)))

Sign() of the times

I was still thinking about how to generate the sign of Pythons in-built permutation generator, and a little thought made me think that somewhere there was a sort routine that sorted by swapping two elements. If you count the swaps then you can calculate the sign if the number of swaps is even or odd.

This lead me to Selection Sort. 

def selectionSort(lst):
    for i in range(0,len(lst)-1):
        mn = min(range(i,len(lst)), key=lst.__getitem__)
        lst[i],lst[mn] = lst[mn],lst[i]
    return lst

Unfortunately the sorter does unnecessary swaps so I made a slight modification to create a routine that generated a sign of any permutation of the integers 0..n:

def perm_parity(lst):
    '''\
    Given a permutation of the digits 0..N in order as a list, 
    returns its parity (or sign): +1 for even parity; -1 for odd.
    '''
    parity = 1
    for i in range(0,len(lst)-1):
        if lst[i] != i:
            parity *= -1
            mn = min(range(i,len(lst)), key=lst.__getitem__)
            lst[i],lst[mn] = lst[mn],lst[i]
    return parity    

if __name__ == '__main__':
    from itertools import permutations

    for p in permutations(range(3)):
        l = list(p)
        print "%2i %r" % (perm_parity(l), p)

I published the sign generator on ActiveState as I didn't think it was enough for a Rosetta Code task.

(If you run the code you can see that their is in fact a pattern to the signs generated by successive members Pythons permutations).

Patterns

If you look at a few of the series generated by the Steinhaus–Johnson–Trotter Python program you can see patterns that point to a recursive method of generating the same sequences. I coded that up too.

Serendipity doo dah

A few days later I came across a Stack-exchange question that needed to generate all permutations when there were duplicate items in the initial list. (It goes on to need the permutations in a specific order, but I wasn't thinking of that bit at first).

Of coure I started by generating possible perm then using a set to remove duplicates, but it gets very (very) wasteful, so I started looking around for an algorithm to generate what are sometimes called lexicographically correct permutations. for example the l-perms of [1,1,1,2,2,2] are the twenty terms:

1 1 1 2 2 2
1 1 2 1 2 2
1 2 1 1 2 2
2 1 1 1 2 2
2 1 1 2 1 2
1 2 1 2 1 2
1 1 2 2 1 2
1 2 2 1 1 2
2 1 2 1 1 2
2 2 1 1 1 2
2 2 1 1 2 1
2 1 2 1 2 1
1 2 2 1 2 1
1 1 2 2 2 1
1 2 1 2 2 1
2 1 1 2 2 1
2 1 2 2 1 1
1 2 2 2 1 1
2 2 1 2 1 1
2 2 2 1 1 1

Which is a lot less than 6! = 720 terms of the standard permutation.

I spent an hour playing with my recursive permutation generator until I hit on the idea that if you had to do permutations of 123XX then you could think of each (identical) X as being X0 and X1, and that valid l-perms are the perms when X0 and X1 appear in only one order.

Now the recursive algorithm for perms introduces each item of the initial perm one-by-one. If you restrict each X as it is introduced to never be introduced 'above' and X that is already in the partial perm item that is being expanded, then, in effect, X0 and X1 can never 'cross' and so will stay in the same 'order'.

I made the change and it seems to work for me.

Lexicographic permutation generator in Python

def l_perm(items):
    if not items:
        return [[]]
    else:
        dir = 1
        new_items = []
        this = [items.pop()]
        for item in l_perm(items):
            lenitem = len(item)
            try:
                # Never insert 'this' above any other 'this' in the item 
                maxinsert = item.index(this[0])
            except ValueError:
                maxinsert = lenitem
            if dir == 1:
                # step down
                new_items += [item[:i] + this + item[i:] for i in range(lenitem, -1, -1)
                              if i <= maxinsert]
            else:    
                # step up
                new_items += [item[:i] + this + item[i:] for i in range(lenitem + 1)
                              if i <= maxinsert]
            dir *= -1
        return new_items
 
 
if __name__ == '__main__':
    n = [1, 1, 1, 2, 2, 2]
    print '\nLexicograpic Permutations of %i items: %r' % (len(n), n)
    for i, x in enumerate(l_perm(n)):
        print('%3i %r' % (i, x))

End bit.

I'm tired and I want to go to bed so turning l-perm into a RC task or sticking it on Activestate can wait. I am wondering if their is anyone else interested in this as I just pottered around the less deep bits.


- Paddy.