Calculating permutations

A common problem in statistics asks you to compute (and evaluate) all the permutations of an input sequence.

Ask yourself: Do I understand the problem?

Input is an array of N objects. We don't know N yet, but obviously the algorithm, whatever it is, shouldn't depend on a particular N.
Without loss of generality, let's assume all values are chars. We shall check at the end what happens if our values are different things - like integers, files, complex numbers or words in a document.
Output is a list of all permutations (all n! of them)

Setting up a test harness

As usual, our first step is to prepare a testbed. What we want to do is verify that our solution is correct, and the easiest way to verify that is to compare it to a known solution. Now, the Python itertools contains a permutations() function that returns all known permutations. There is a slight problem, however: we don't know if the order of the permutations returned by our code is identical to the one returned by the builtin method. Therefor we'll sort both results and compare them sorted. Our test strategy therefor will be this:

Calculate all permutations using our method
Sort our result
Calculate all permutations using the itertools method
Sort that result
Compare the two results and fail if they are not equal

This is the code:

#! -*- Encoding: Latin-1 -*-

def selftest():
    A = list("pyalda")

    received_permutations = permutation(A)
    received_permutations.sort()

    expected_permutations = [p for p in itertools.permutations(A)]
    expected_permutations.sort()

    assert received_permutations == expected_permutations

if __name__ == "__main__":
    selftest()

Using brute-force

The basic strategy to calculate all permutations is to evaluate all possible combinations of the input alphabet. What is the input alphabet?

Input is an array of N objects. They happen to be chars, but that is just an artifact of our test data.
The input alphabet are the individual objects in the input array

What are possible combinations?

A possible combination must use all items of the input alphabet, and
A possible combination must use each item only once

What sort of memory do we need?

We need an output list that will recieve all the permutations
For each current permutation:
- We need an array of N items to represent the current permutation
- We need to check condition (2) by using some temporary memory that tells us if a item is already used in the current permutation.
- For condition (1), we need to enumerate all items in the input alphabet.

That sounds doable, but there is one catch: how do we ensure we catch all possible combinations? It turns out that recursion will save the day:

For input position 0: choose all possible values (which are all values in the input alphabet), once
For input position 1: choose all possible values (except the ones already in use)
...
For input position n-1: choose all possible values (except the ones already in use)
For input position n: we've reached a possible combination = a permutation, store it

This logic can be represented pretty easily in Python:

def permutation(A):
    result = []
    N = len(A)

    # Condition (2): ensure that each item in the input alphabet is used only once
    used_in_current_permutation = [False] * N

    # temporary memory for the current permutation
    current_permutation = [None] * N

    # recursive logic, see above
    def recurse(n):
        if n == N:
            # we've reached the last position: we have a possible permutation
            result.append( tuple(current_permutation) )

        else:
            # try all input values still available and see what happens if we use them at the current position
            for index, possible_digit in enumerate(A):
                if not used_in_current_permutation[index]:
                    used_in_current_permutation[index] = True
                    current_permutation[n] = possible_digit
                    recurse(n+1)
                    used_in_current_permutation[index] = False

    recurse(0)
    return result

Permuting the elements of an array without extra memory

Now, in the very fine Elements of Programming Interviews problem 3.10 asks an interesting question: how to apply a given permutation to an array of input values - without additional memory. For me it turned out to be a much more interesting question than calculating all permutations (but then maybe this is so because I already knew the former).

Ask yourself: Do I understand the problem?

Input is
1. an array of objects, and
2. an array of permutation indices
Output is the very same input array, with the permutation applied

Setting up a test harness

As usual, our first step is to prepare a testbed. This time we'll take it easy: we'll provide a test array and a test permutation, and check out what happens. (A more complex testbed will compare different results and timings). This is the code:

#! -*- Encoding: Latin-1 -*-

def selftest():
    A = list('ABCDEF')
    P = [5,1,3,2,4,0]

    print apply_permutation(A, P)

if __name__ == "__main__":
    selftest()

Applying the permutation if additional memory is available

If you can allocate an output array, this is a trivial task:

def apply_permutation(A, P):
    output = [None] * len(A)
    for index in range(len(A)):
        output[index] = A[P[index]]

    return output

Actually, this being Python you might as well use list comprehension

def apply_permutation(A, P):
    return [A[P[index]] for index in range(len(A))]

Or, better yet:

def apply_permutation(A, P):
    return [A[p] for p in P]

OK, that was fun, but it doesn't solve the question posed: how to do this without additional memory.

Cheating using Python

Take a look at the input provided. You have two arrays:

A is the input array, which is to be modified inplace
P is the array of permutation indices

The key thing to notice is that we don't need P later, so we can just use P as scratch space

ONE WORD OF CAUTION: This will only work in Python. Or, it will only work for other programming languages if both A and P are arrays of the same type. This is so because in "type-safe" languages like C/C++ or Java you cannot store objects in int-arrays, like you can in Python.

OK, let's try it:

def apply_permutation(A, P):
    assert len(A) == len(P)
    
    N = len(A)
    # apply the permutation, storing the result in P (overwriting the data at P)
    for index, item in enumerate(A):
        P[index] = A[P[index]]

    # copy P back to A
    for index, item in enumerate(A):
        A[index] = P[index]

Now, we have followed the letter of the original problem, but not the spirit: for two reasons:

we've implicitly used O(N) space
we've destroyed the information in P

We'll turn to Problem 1 later, but also Problem 2 is interesting in its own right: how can we restore the original permutation?

Whenever we write something to P, we read a copy of A at a certain position k
Because we're mapping a permutation, we can be sure that A[k] isn't used again, so we can use that space to store the information of what we just did.
In the final step, we don't simply copy P back to A, we swap P with A

This is what the code looks like:

def apply_permutation(A, P):
    N = len(A)

    # Step 1: apply the permutation, storing the result in P (overwriting the data at P)
    for index, item in enumerate(A):
        p = P[index] 
        P[index] = A[p]

        # NEW: keep information on P in A
        A[p] = index
    
    # Step 2: swap items of P and A
    for index, item in enumerate(A):
        A[index], P[index] = P[index], A[index]

Now, there is a slight problem that can best be seen by looking at a table:

    
        
            Step A P
        
        
            Initial Input ['C', 'A', 'B'] [2,0,1]
        
        
            Step 1 [1, 2, 0] ['C', 'A', 'B']
        
        
            Step 2 ['C', 'A', 'B'] [1, 2, 0]

Step	A	P
Initial Input	['C', 'A', 'B']	[2,0,1]
Step 1	[1, 2, 0]	['C', 'A', 'B']
Step 2	['C', 'A', 'B']	[1, 2, 0]

We have successfully applied the permutation, and we've stored some information regarding the permutation in P - but not the original one!. How can we restore that?

Ask yourself: what does the information in the original P=[2,0,1] mean?

It means that
- A[index=0] should be A[P[index=0]=2],
- A[index=1] should be A[P[index=1]=0], and
- A[index=2] should be A[P[index=2]=1]

Ask yourself: what does the information in the modified P=[1,2,0] mean?

It means that
- P[index=0] should be P[P[index=0]=1] = 2,
- P[index=1] should be P[P[index=1]=2] = 0, and
- P[index=2] should be P[P[index=2]=0] = 1

Good! There is only a slight problem with that formulation:

    
        
            Index P before step What should we do P after step
        
        
            0 [1,2,0] P[index=0] should be P[P[index=0]=1] = 2 [2,2,0]
        
        
            1 [2,2,0] P[index=1] should be P[P[index=1]=2] = 0 [2,0,0]
        
        
            2 [2,0,0] P[index=2] should be P[P[index=2]=0] = 0 [2,0,0]

Index	P before step	What should we do	P after step
0	[1,2,0]	P[index=0] should be P[P[index=0]=1] = 2	[2,2,0]
1	[2,2,0]	P[index=1] should be P[P[index=1]=2] = 0	[2,0,0]
2	[2,0,0]	P[index=2] should be P[P[index=2]=0] = 0	[2,0,0]

Actually, everything was lost already after the first step: because the information that there was a 1 in the array was lost.

But do not despair: there is a simple trick to fix this: we can do two passes: in the first one we restore the information but store it somewhere else, in the second pass we use that updated information. The only question: where do we store it?

Well, this is where some bit of math comes to our rescue: we can treat P[i] as a number in the N-based number system :). Let's try to make this idea more clear with an example:

Input is an array of 10 values
So by definition all P[i] will be in the range 0..9 - like decimal digits
Assume we want to store digit x in P[i]: then we simply set P[i] = P[i] + 10*x
To restore the original P[i], we use P[i] mod 10, and to find the x we use P[i] div 10

This idea trivially generalizes to any sort of number system: and our number system has N digits, so it is a N-based number system. Let's look at the code for this:

def apply_permutation(A, P):
    N = len(A)

    # Step 1: apply the permutation, storing the result in P (overwriting the data at P)
    for index, item in enumerate(A):
        p = P[index] 
        P[index] = A[p]
        A[p] = index
    
    # Step 2: swap items of P and A
    for index, item in enumerate(A):
        A[index], P[index] = P[index], A[index]

    # Step 3: encode information on the target value in the N-based number system
    for index, p in enumerate(P):
        l = p % N
        P[l] = (P[l] % N) + index*N

    # Step 4: reset P to its original values
    for index, p in enumerate(P):
        P[index] = p / N

    return A

Cheating while preserving type-safety

Now, let's start this slowly.

If we have a permutation, we can shift the elements with only constant additional memory - a single temporary variable, plus an index. This is best explained by looking at an example:

    
        
            Input: Array A A B C Temp Value Temp Index
        
        
            Input: Permutation P 2 0 1    
        
        
            Expected Output A: C A B    
        
        
            Step 1: Start with P[0] C B C A[0] = 'A' 0
        
        
            Step 2: Find P[x] = Temp Index = 0 C A C A[1] = 'B' 1
        
        
            Step 3: Find P[x] = Temp Index = 0 C A B A[2] = 'C' 2
        
        
            Step 4: Done (Temp Index = Initial Index) C A B

Code that implements this shift looks like this:

def apply_permutation(A, P):
    for current_index in range(0, len(A)):

        # initial value: A[x] = A[P[x]]
        temp_index = current_index
        temp_value = A[current_index]
        A[current_index] = A[P[current_index]]
        while True:

            # find next position in loop
            for k in range(0, len(A)):
                if P[k] == temp_index:
                    A[k], temp_value = temp_value, A[k]
                    temp_index = k
                    break

            # we're done if the last item of the permutation loop has been found
            if temp_index == current_index:
                break
    return A

There is only a slight problem with this code: it doesn't work. You can see why here:

    
        
            A: 'A' 'B' 'C'
        
        
            P: 2 0 1
        
        
            current_index = 0 'C' 'A' 'B'
        
        
            current_index = 1 'B' 'C' 'A'
        
        
            current_index = 2 'A' 'B' 'C'

So what happened here is that the permutation was applied multiple times, because the code didn't detect that this bit of code had already been seen. What we need is a way to store the information "permutation P[x] has already been applied" that works without either destroying the information in P, or using additional memory. A simple way to do so is by negating the number. If it is negative, we've seen it, if it is positive, we haven't. Thus our improved code will look something like this:

def apply_permutation(A, P):
    N = len(A)
    for current_index in range(N):
        if P[current_index] >= 0:
            temp_index = current_index
            temp_value = A[current_index]
            A[current_index] = A[P[current_index]]

            # mark negative = visited
            P[current_index] -= N
            found = True
            while found:
                found = False
                for k in range(0, len(A)):
                    if P[k] == temp_index:

                        # mark negative = visited
                        P[k] -= N
                        A[k], temp_value = temp_value, A[k]
                        temp_index = k
                        found = True
    
    # restore P
    for k in range(N):
        if P[k] < 0:
            P[k] += N

    return A

This code correctly applies permutation P to A, and restores the information in P.

Note: this code is O(N²), because it looks for the update index in the permutation. EPI contains a solution that is right - but nontrivially so. It is right because it does not use the normal permutation, but its cyclic representation. That is kind of cheating: getting the cyclic representation out of a given permutation in non-cyclic form requires O(N) space and O(N²) time, so it kind of defeats the purpose. It is true that once given, the cyclic representation can be applied in O(N) time (if you use the negative-numbers trick from above), but it feels ... not right. Anyway, more information on getting the cycling representation (and a general good read) is Cyclesort - a curious little sorting algorithm

Input: Array A	A	B	C	Temp Value	Temp Index
Input: Permutation P	2	0	1
Expected Output A:	C	A	B
Step 1: Start with P[0]	C	B	C	A[0] = 'A'	0
Step 2: Find P[x] = Temp Index = 0	C	A	C	A[1] = 'B'	1
Step 3: Find P[x] = Temp Index = 0	C	A	B	A[2] = 'C'	2
Step 4: Done (Temp Index = Initial Index)	C	A	B

A:	'A'	'B'	'C'
P:	2	0	1
current_index = 0	'C'	'A'	'B'
current_index = 1	'B'	'C'	'A'
current_index = 2	'A'	'B'	'C'