From Werner KRAUTH
import numpy as np
import itertools
import scipy.linalg as la
def factorial(n):
return 1 if n == 0 else (0 if n == 0 else factorial(n - 1) * n)
for N in [2, 3, 4, 5, 6, 7]:
FacN = factorial(N)
ConfCopy = [0] * N
print(N, 'N', FacN)
#
# Setup of transition matrix
#
P = np.zeros((FacN, FacN))
Conf = [k for k in range(N)]
ConfList = list(itertools.permutations(Conf))
for Conf in ConfList:
i = ConfList.index(tuple(Conf))
ConfCopy[:] = Conf
a = ConfCopy.pop(0)
for k in range(N):
TargetConf = ConfCopy[0:k] + [a] + ConfCopy[k:N - 1]
j = ConfList.index(tuple(TargetConf))
P[i][j] = 1.0 / float(N)
eigvals, eigvecsl, eigvecsr = la.eig(P, left=True)
eigvals.sort()
stats = [0] * (N+1)
for a in eigvals:
index = int(N * a.real + 0.5)
stats[index] += 1
print(stats)