# Go deh!

Mainly Tech projects on Python and Electronic Design Automation.

## Tuesday, July 31, 2012

### A permutation journey

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

$\det(A) = \sum_\sigma\sgn(\sigma)\prod_{i=1}^n M_{i,\sigma_i}$
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):
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))

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.

1. I just answered this Stack overflow question with a revised version of l_perm that is now an iterative generator function.

2. As a player of Go, your post titles drive me nuts. What does 'deh' mean?

1. Go deh is Jamaican patois meaning 'go there' literally, but also used as an expression meaning great! Keep on doing the great job that you are doing - something like that although I find it hard to express the feelings that go with the term when spoken between those of Jamaican ancestry.
.