# TVDTemperingLift.py

(Difference between revisions)
 Revision as of 21:56, 16 March 2022Werner (Talk | contribs)← Previous diff Revision as of 10:26, 8 September 2022Werner (Talk | contribs) Next diff → Line 1: Line 1: - This is the program TVDTemperingLift.py, which allows one to compute the TVD as a function of time for the simulated tempering problem with the flat and the V-shaped probability distribution. - - # - # Simulated tempering for the V-shaped stationary distribution --- Lifting - # algorithm - # import random import random import pylab import pylab import numpy as np import numpy as np - import scipy.linalg as la + for n in [10, 20, 40, 80, 160, 320]: - ReplicaChange = 0.1 + ReplicaChange = 0.1 - n = 100 + const = 4.0 / n ** 2 - const = 4.0 / n ** 2 + PiStat = {} - PiStat = {} + Table = [] - Table = [] + - + for x in range(1, n + 1): - for x in range(1, n + 1): + Table.append((x, -1, 0)) - Table.append((x, -1, 0)) + Table.append((x, 1, 0)) - Table.append((x, 1, 0)) + Table.append((x, -1, 1)) - Table.append((x, -1, 1)) + Table.append((x, 1, 1)) - Table.append((x, 1, 1)) + # # # factor of 1/4 because the total must be normalized, with two liftings and # factor of 1/4 because the total must be normalized, with two liftings and # two replicas # 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, 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[(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, 0)] = 0.0 + PiStat[(0, 1, 0)] = 0.0 - PiStat[(0, -1, 1)] = 0.0 + PiStat[(0, -1, 1)] = 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, 0)] = 0.0 + PiStat[(n + 1, 1, 0)] = 0.0 - PiStat[(n + 1, -1, 1)] = 0.0 + PiStat[(n + 1, -1, 1)] = 0.0 - PiStat[(n + 1, 1, 1)] = 0.0 + PiStat[(n + 1, 1, 1)] = 0.0 - Position = (1, 1, 0) + Position = (1, 1, 0) - print(Table) + PTrans = np.zeros((4 * n, 4 * n)) - PTrans = np.zeros((4 * n, 4 * n)) + Pi = np.zeros([4 * n]) - for x in range(1, n + 1): + for x in range(1, n + 1): - for Dir in [-1, 1]: + for Dir in [-1, 1]: - for Rep in [0, 1]: + for Rep in [0, 1]: - i = Table.index((x, Dir, Rep)) + i = Table.index((x, Dir, Rep)) - k = Table.index((x, -Dir, Rep)) + Pi[i] = PiStat[(x, Dir, Rep)] - if PiStat[(x + Dir, Dir, Rep)] > 0.0: + k = Table.index((x, -Dir, Rep)) - j = Table.index((x + Dir, Dir, Rep)) + if PiStat[(x + Dir, Dir, Rep)] > 0.0: - PTrans[i, j] = min(1.0, PiStat[(x + Dir, Dir, Rep)] / PiStat[(x, Dir, Rep)]) + j = Table.index((x + Dir, Dir, Rep)) - PTrans[i, k] = 1.0 - PTrans[i, j] + PTrans[i, j] = min(1.0, PiStat[(x + Dir, Dir, Rep)] / PiStat[(x, Dir, Rep)]) - else: + PTrans[i, k] = 1.0 - PTrans[i, j] - PTrans[i, k] = 1.0 + else: - print('PTrans') + PTrans[i, k] = 1.0 - print(PTrans) + + 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) - PReplica = np.zeros((4 * n, 4 * n)) + P = PTrans @ PReplica @ PResampling - for x in range(1, n + 1): + Pit = np.zeros([4 * n]) - for Dir in [-1, 1]: + Pit[0] = 1.0 - i = Table.index((x, Dir, 0)) + xvalues = [] - j = Table.index((x, Dir, 1)) + yvalues = [] - PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, Dir, 1)] / PiStat[(x, Dir, 0)]) + iter = 0 - PReplica[i, i] = 1.0 - PReplica[i, j] + while True: - PReplica[j, i] = ReplicaChange * min(1.0, PiStat[(x, Dir, 0)] / PiStat[(x, Dir, 1)]) + iter += 1 - PReplica[j, j] = 1.0 - PReplica[j, i] + Pit = Pit @ P - print('PReplica') + TVD = sum(np.absolute(Pi - Pit) / 2.0) - print(PReplica) + xvalues.append(iter / float(n)) - P = PTrans @ PReplica + yvalues.append(TVD) - eigvals, eigvecsl, eigvecsr = la.eig(P, left=True) + if TVD < 0.1: break - eigvals = eigvals.real + pylab.plot(xvalues,yvalues, label='\$n =\$ '+str(n)) - print('Eigenvalues') + pylab.legend(loc='upper right') - print(eigvals) + pylab.xlabel("\$t/ n\$ (rescaled time) ") - gap = 1000.0 + pylab.ylabel("TVD") - for k in range(4 * n): + pylab.title("TVD lift tempering on the path graph of \$n\$ sites") - if abs(eigvals[k] - 1.0) < 1.e-5: max_val = k + - else: + - gap = min(gap, abs(eigvals[k] - 1.0)) + - print('gap', gap, ReplicaChange) + - print('maxval', max_val) + - PiStat = eigvecsl[:, max_val] + - PiStat = PiStat.real + - PiStatSum = PiStat.sum() + - PiStat0 = [] + - PiStat1 = [] + - for k in range(4 * n): + - if k % 4 == 0: + - PiStat0.append(4.0 * PiStat[k] / PiStatSum) + - elif k % 4 == 2: + - PiStat1.append(4.0 * PiStat[k] / PiStatSum) + - XVec = [k for k in range(1, n + 1)] + - print(PiStat0) + - pylab.plot(XVec, PiStat0) + - pylab.plot(XVec, PiStat1) + - pylab.xlabel("\$x\$ (position) ") + - pylab.title("Simulated tempering (Lifted MCMC), Left EV transition matrix") + pylab.show() pylab.show()

## Revision as of 10:26, 8 September 2022

```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()
```