TVDTemperingLift.py
From Werner KRAUTH
This is the final program TVDTemperingLift.py to be discussed in my fourth lecture at the summerschool .... at the International Center for Theoretical Physics in Trieste, in September 2022.
In previous lecture, we had started the discussion of the diffusion on the path graph P_n, in other words the one-dimensional linear lattice with n vertices and without periodic boundary conditions.
import random import pylab import numpy as np for n in [10, 20, 40, 80, 160, 320]: ReplicaChange = 0.1 const = 4.0 / n ** 2 PiStat = {} Table = [] for x in range(1, n + 1): Table.append((x, -1, 0)) Table.append((x, 1, 0)) Table.append((x, -1, 1)) Table.append((x, 1, 1)) # # factor of 1/4 because the total must be normalized, with two liftings and # two replicas # PiStat[(x, -1, 0)] = 1.0 / float(n) / 4.0 PiStat[(x, 1, 0)] = 1.0 / float(n) / 4.0 PiStat[(x, -1, 1)] = const * abs( (n + 1) / 2 - x) / 4.0 PiStat[(x, 1, 1)] = const * abs( (n + 1) / 2 - x) / 4.0 PiStat[(0, -1, 0)] = 0.0 PiStat[(0, 1, 0)] = 0.0 PiStat[(0, -1, 1)] = 0.0 PiStat[(0, 1, 1)] = 0.0 PiStat[(n + 1, -1, 0)] = 0.0 PiStat[(n + 1, 1, 0)] = 0.0 PiStat[(n + 1, -1, 1)] = 0.0 PiStat[(n + 1, 1, 1)] = 0.0 Position = (1, 1, 0) PTrans = np.zeros((4 * n, 4 * n)) Pi = np.zeros([4 * n]) for x in range(1, n + 1): for Dir in [-1, 1]: for Rep in [0, 1]: i = Table.index((x, Dir, Rep)) Pi[i] = PiStat[(x, Dir, Rep)] k = Table.index((x, -Dir, Rep)) if PiStat[(x + Dir, Dir, Rep)] > 0.0: j = Table.index((x + Dir, Dir, Rep)) PTrans[i, j] = min(1.0, PiStat[(x + Dir, Dir, Rep)] / PiStat[(x, Dir, Rep)]) PTrans[i, k] = 1.0 - PTrans[i, j] else: PTrans[i, k] = 1.0 PReplica = np.zeros((4 * n, 4 * n)) for x in range(1, n + 1): for Dir in [-1, 1]: i = Table.index((x, Dir, 0)) j = Table.index((x, Dir, 1)) PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, Dir, 1)] / PiStat[(x, Dir, 0)]) PReplica[i, i] = 1.0 - PReplica[i, j] PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, Dir, 0)] / PiStat[(x, Dir, 1)]) PReplica[j, j] = 1.0 - PReplica[j, i] PResampling = np.zeros((4 * n, 4 * n)) for x in range(1, n + 1): for Dir in [-1, 1]: i = Table.index((x, Dir, 0)) j = Table.index((x, -Dir, 0)) k = Table.index((x, Dir, 1)) l = Table.index((x, -Dir, 1)) PResampling[i, j] = 1.0 / (4.0 * n) PResampling[i, i] = 1 - 1. / (4.0 * n) PResampling[k, l] = 1.0 / (4.0 * n) PResampling[k, k] = 1 - 1. / (4.0 * n) P = PTrans @ PReplica @ PResampling Pit = np.zeros([4 * n]) Pit[0] = 1.0 xvalues = [] yvalues = [] iter = 0 while True: iter += 1 Pit = Pit @ P TVD = sum(np.absolute(Pi - Pit) / 2.0) xvalues.append(iter / float(n)) yvalues.append(TVD) if TVD < 0.1: break pylab.plot(xvalues,yvalues, label='$n =$ '+str(n)) pylab.legend(loc='upper right') pylab.xlabel("$t/ n$ (rescaled time) ") pylab.ylabel("TVD") pylab.title("TVD lift tempering on the path graph of $n$ sites") pylab.show()