# TVDMetroLift.py

```#
# TVD for the Lifted Metropolis algorithm on the path graph.
#
import random
import pylab
import numpy as np
model = 'Flat'
model = 'VShape'
for n in [10, 20, 40, 80, 160, 320]:
const = 4.0 / n ** 2
PiStat = {}
Table = []
for x in range(1, n + 1):
Table.append((x, -1))
Table.append((x,  1))
if model == 'Flat':
PiStat[(x, -1)] = 1.0 / float(n) / 2.0
PiStat[(x, +1)] = 1.0 / float(n) / 2.0
elif model == 'VShape':
PiStat[(x, -1)] = const * abs( (n + 1) / 2 - x) / 2.0
PiStat[(x,  1)] = const * abs( (n + 1) / 2 - x) / 2.0
PiStat[(0, -1)] = 0.0
PiStat[(0,  1)] = 0.0
PiStat[(n + 1, -1)] = 0.0
PiStat[(n + 1,  1)] = 0.0
PTrans   = np.zeros((2 * n, 2 * n))
Pi = np.zeros([2 * n])
for x in range(1, n + 1):
for Sigma in [-1, 1]:
i = Table.index((x, Sigma))
Pi[i] = PiStat[(x, Sigma)]
k = Table.index((x, -Sigma))
if PiStat[(x + Sigma, Sigma)] > 0.0:
j = Table.index((x + Sigma, Sigma))
PTrans[i, j] = min(1.0, PiStat[(x + Sigma, Sigma)] / PiStat[(x, Sigma)])
PTrans[i, k] = 1.0 - PTrans[i, j]
else:
PTrans[i, k] = 1.0
PResampling = np.zeros((2 * n, 2 * n))
for x in range(1, n + 1):
for Sigma in [-1, 1]:
i = Table.index((x, Sigma))
j = Table.index((x, -Sigma))
PResampling[i, j] = 1.0 / n
PResampling[i, i] = 1.0 - PResampling[i,j]

P = PTrans @ PResampling
Pit = np.zeros([2 * n])
Pit[0] = 1.0
xvalues = []
yvalues = []
iter = 0
while True:
iter += 1
Pit = np.array(Pit)
Pit = Pit@P
TVD = sum(np.absolute(Pi - Pit) / 2.0)
xvalues.append(iter / float(n ** 2))
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 for the lifted Metropolis algorithm on the path graph of \$n\$ sites")
pylab.show()
```

~