'''
<title>367 Bozo Sort</title>

<meta name="solved-by"  content="Joe Tacheron">
<meta name="difficulty" content="55%">
<meta name="run time"   content="1s">
<meta name="answer"     content="48271207">

<h1>367 Bozo Sort</h1>

<h2>Problem Description</h2>
<pre>
Bozo sort, not to be confused with the slightly less efficient bogo sort, consists out of checking if the input sequence is sorted and if not swapping randomly two elements. This is repeated until eventually the sequence is sorted.

If we consider all permutations of the first 4 natural numbers as input the expectation value of the number of swaps, averaged over all 4! input sequences is 24.75.
The already sorted sequence takes 0 steps.

In this problem we consider the following variant on bozo sort.
If the sequence is not in order we pick three elements at random and shuffle these three elements randomly.
All 3!=6 permutations of those three elements are equally likely.
The already sorted sequence will take 0 steps.
If we consider all permutations of the first 4 natural numbers as input the expectation value of the number of shuffles, averaged over all 4! input sequences is 27.5.
Consider as input sequences the permutations of the first 11 natural numbers.
Averaged over all 11! input sequences, what is the expected number of shuffles this sorting algorithm will perform?

Give your answer rounded to the nearest integer.
</pre>

<h2>Solution Method</h2>
<pre>
First I'll define what I call a "swap list": a description of a permuted sequence given as a list of 2-swaps, 3-swaps, etc. needed to sort the sequence. E.g., the sequence 53421 has a swap list [2,3], indicating one 2-swap (the 1 & 5) and one 3-swap (the 2, 3, & 4); the swap list [] or [0] means the sequence is correctly sorted.

Note that each permutation is described by exactly one swap list, while each each swap list encompasses multiple permutations. In other words, permutations can be partitioned into "swap-classes." Also, note that after a Bozo swap, the transition probability distribution (from one swap-class to another) is the same for all permutations of the same swap-class.

So with those definitions out of the way, here's my solution method:

1. 	Create a representative permutation for each swap-class

2. 	Calculate its transition probability distribution.

3. 	Solve the "expected number of Bozo swaps" (ENBS) system of equations:

	S(0) = 0
	S(2) = 1 + P(2→2)×S(2) + P(2→2,2)×S(2,2) + ... + P(2→11)×S(11)
	...
	S(11) = 1 + P(11→2)×S(11) + P(11→2,2)×S(2,2) + ... + P(11→11)×S(11)

	where S(x) is the ENBS to go from swap-class [x] to swap-class [], and
	P(x→y) means the probability of transitioning from [x] to [y] in a single Bozo swap.

4. 	Calculate the average ENBS across every permutation. Since the ENBS is the same for every permutation in the same swap-class, you just need to know how large each swap-class is.

	answer = Σ(size(x)*S(x) for each swap-class x) / 11!

Solved in one second in python.
</pre>
'''

import math
import numpy as np
from collections import Counter

def bozo_swap(seq, i, j, k, type):
	'''Bozo swap 3 elements, 'type' selects how.'''
	I,J,K = ((i,j,k),(j,k,i),(k,i,j),(j,i,k),(i,k,j),(k,j,i))[type]
	seq[i], seq[j], seq[k] = seq[I], seq[J], seq[K]
	return seq

def make_swap_lists(length, _all = None, _parent = []):
	'''Create every possible swap list for a sequence of given length.'''
	if _all is None:
		_all = [tuple()]
	for s in range(2, length+1):
		if len(_parent) > 0 and s < _parent[-1]:
			continue
		if sum(_parent) + s > length:
			break
		swaps = _parent.copy()
		swaps.append(s)
		_all.append(tuple(swaps))
		make_swap_lists(length, _all, swaps)
	return _all

def representative_perm(length, swaps):
	'''Create an arbitrary permutation with the given length and swap list.'''
	seq = list(range(length))
	# an n-swap can be made by n-1 2-swaps
	i = 0
	for s in swaps:
		for j in range(s-1):
			seq[i+j], seq[i+j+1] = seq[i+j+1], seq[i+j]
		i += s
	return seq

def count_perms(length, swaps):
	'''Count the number of permutations with the given swap list.'''
	total = 1
	for s in swaps:
		total *= math.comb(length, s) * math.factorial(s-1) # choose s items & permute them
		length -= s # chosen s items cannot be chosen again
	# divide by repeated swaps
	# e.g. a [2,2] swap list double-counts; [3,3,3] overcounts by a factor of 3!
	for swap,count in Counter(swaps).items():
		total //= math.factorial(count)
	return total

def swap_list(seq):
	'''Identify the swap list from the given permutation.'''
	swaps = []
	for i in range(len(seq)):
		if seq[i] in [-1, i]:
			continue
		s = 0
		while seq[i] != -1:
			seq[i], i = -1, seq[i]
			s += 1
		swaps.append(s)
	swaps.sort()
	return tuple(swaps)

def transition_dist(length, swaps):
	'''
	Calculate the probability distribution for an arbitrary permutation of a given swap-class to transition into another swap-class after undergoing a Bozo swap.
	'''
	seq = representative_perm(length, swaps)
	dist = Counter()	
	# count swap-classes visited after a Bozo swap
	for i in range(len(seq)-2):
		for j in range(i+1, len(seq)-1):
			for k in range(j+1, len(seq)):
				for type in range(6):
					swapped = bozo_swap(seq.copy(), i, j, k, type)
					s = swap_list(swapped)
					dist[s] += 1
	# convert to list of probabilities
	total = sum(dist.values())
	dist = [dist[s]/total for s in make_swap_lists(len(seq))]
	return dist

def transition_matrix(length):
	'''Compile transition distributions into a matrix.'''
	swap_lists = make_swap_lists(length)	
	mat = []
	for swaps in swap_lists:
		mat.append(transition_dist(length, swaps))
	mat = np.matrix(mat)
	# The first row is special in that the transition probability from a
	# correctly ordered sequence to any other needs to be 0, including [] -> [].
	mat[0, :] = 0
	return mat

def expected_swap_count(length):
	'''Solve the expected swap count system of equations.'''
	A = transition_matrix(length)
	A = A - np.eye(A.shape[0])
	b = -1 * np.ones(A.shape[0])
	b[0] = 0 # 0 swaps needed for the sorted sequence
	expected_swap_count = np.linalg.solve(A, b)	
	# average "expected swap count" across all permutations
	class_size = [count_perms(length, swaps) for swaps in make_swap_lists(length)]
	average = expected_swap_count.dot(class_size) / math.factorial(length)
	return average

assert math.isclose(expected_swap_count(4), 27.5)
ans = round(expected_swap_count(11))
print(ans)