Top to random TVD.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 [3, 4, 5, 6, 7, 8, 9, 10]:
    print(N, 'N')
    XValues = []
    YValues = []
#
# lifted TASEP: sample a random initial active particle
#
    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):
#    print(Pit)
        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 / N ** 2')
        XValues.append(step / N / math.log(N))
        YValues.append(TVD)
        if TVD < 0.02: break
        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.plot(XValues, YValues, label = '$N= $' + str(N))
plt.title('TVD of Top-to-random shuffle')
plt.xlabel('$t /( N \log N)$ (rescaled time)')
plt.ylabel('TVD')
plt.legend(loc='upper center')
plt.show()