TVDTemperingLift.py

From Werner KRAUTH

Jump to: navigation, search

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