Top to random TVD mix rel.py

From Werner KRAUTH

import math
import copy
import matplotlib.pyplot as plt

def factorial(n):
    return 1 if n == 0 else (0 if n == 0 else factorial(n - 1) * n)

for N in [8]:
    print(N, 'N')
    XValues = []
    YValues = []
#
#   computing the TVD
#
    Pit = {}
    Conf = [k for k in range(N)]
    Pit[tuple(Conf)] = 1.0
    NConfs = factorial(N)
    weight = 1.0 / NConfs
    DiameterFlag = True
    for step in range(40000):
        TVD = (NConfs - len(Pit)) * weight / 2.0
        for k in Pit:
            TVD += abs(Pit[k] - weight) / 2.0
#        print('============', N, step, len(Pit) / NConfs, TVD, '==============')
        if len(Pit) == NConfs and DiameterFlag == True:
            DiameterFlag = False
            print(N, step, 'N, Diameter')
        XValues.append(step)
        YValues.append(TVD)
        if TVD < 0.0000001: break
        if TVD < 1.0 / 5.4: tmix = step
        NewPit = {}
        for OldConf, OldPi in Pit.items():
            OldConf = list(OldConf)
            a = OldConf.pop(0)
            for k in range(N):
                Conf = OldConf[0:k] + [a] + OldConf[k:N-1]
                if tuple(Conf) in NewPit.keys():
                    NewPit[tuple(Conf)] += OldPi / N
                else:
                    NewPit[tuple(Conf)] = OldPi / N
        Pit = copy.deepcopy(NewPit)
    plt.semilogy(XValues, YValues, label = '$N= $' + str(N))
    tau_rel = N  / 2.3
    for i in range(len(XValues)):
        X = XValues[i]
        YValues[i] =  (1 - 2.0 / N) ** X
    plt.semilogy(XValues, YValues, label = 'Rel $N= $' + str(N))
    for i in range(len(XValues)):
        X = XValues[i]
        YValues[i] =  (1.0 / 2.0 ) ** (X / tmix)
    plt.semilogy(XValues, YValues, label = 'Mix N= ' + str(N))

plt.title('TVD of Top-to-random shuffle')
plt.xlabel('$t$ (non-rescaled time)')
plt.ylabel('TVD')
plt.legend(loc='upper center')
plt.show()