Top to random TVD mix rel.py

From Werner KRAUTH

Context

This page is part of my 2024 Beg Rohu Lectures on "The second Markov chain revolution" at the Summer School "Concepts and Methods of Statistical Physics" (3 - 15 June 2024).

Python program

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()