Top to random TVD mix rel.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 00:08, 6 June 2024 Werner (Talk | contribs) ← Previous diff |
Current revision Werner (Talk | contribs) |
||
| Line 1: | Line 1: | ||
| + | ==Context== | ||
| + | This page is part of my [[BegRohu_Lectures_2024|2024 Beg Rohu Lectures]] on "The second Markov chain revolution" at the [https://www.ipht.fr/Meetings/BegRohu2024/index.html Summer School] "Concepts and Methods of Statistical Physics" (3 - 15 June 2024). | ||
| + | |||
| + | ==Python program== | ||
| + | |||
| import math | import math | ||
| import copy | import copy | ||
Current revision
[edit]
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).
[edit]
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()
