TVDTemperingLift.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 21:56, 16 March 2022 Werner (Talk | contribs) ← Previous diff |
Current revision Werner (Talk | contribs) |
||
| 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. | + | 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. | ||
| + | |||
| - | # | ||
| - | # 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)) | |
| - | PReplica = np.zeros((4 * n, 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]: | + | i = Table.index((x, Dir, 0)) |
| - | i = Table.index((x, Dir, 0)) | + | j = Table.index((x, Dir, 1)) |
| - | j = Table.index((x, Dir, 1)) | + | PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, Dir, 1)] / PiStat[(x, Dir, 0)]) |
| - | PReplica[i, j] = ReplicaChange * min(1.0, PiStat[(x, Dir, 1)] / PiStat[(x, Dir, 0)]) | + | PReplica[i, i] = 1.0 - PReplica[i, j] |
| - | 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, i] = ReplicaChange * min(1.0, PiStat[(x, Dir, 0)] / PiStat[(x, Dir, 1)]) | + | PReplica[j, j] = 1.0 - PReplica[j, i] |
| - | PReplica[j, j] = 1.0 - PReplica[j, i] | + | PResampling = np.zeros((4 * n, 4 * n)) |
| - | print('PReplica') | + | for x in range(1, n + 1): |
| - | print(PReplica) | + | for Dir in [-1, 1]: |
| - | P = PTrans @ PReplica | + | i = Table.index((x, Dir, 0)) |
| - | eigvals, eigvecsl, eigvecsr = la.eig(P, left=True) | + | j = Table.index((x, -Dir, 0)) |
| - | eigvals = eigvals.real | + | k = Table.index((x, Dir, 1)) |
| - | print('Eigenvalues') | + | l = Table.index((x, -Dir, 1)) |
| - | print(eigvals) | + | PResampling[i, j] = 1.0 / (4.0 * n) |
| - | gap = 1000.0 | + | PResampling[i, i] = 1 - 1. / (4.0 * n) |
| - | for k in range(4 * n): | + | PResampling[k, l] = 1.0 / (4.0 * n) |
| - | if abs(eigvals[k] - 1.0) < 1.e-5: max_val = k | + | PResampling[k, k] = 1 - 1. / (4.0 * n) |
| - | else: | + | |
| - | gap = min(gap, abs(eigvals[k] - 1.0)) | + | P = PTrans @ PReplica @ PResampling |
| - | print('gap', gap, ReplicaChange) | + | Pit = np.zeros([4 * n]) |
| - | print('maxval', max_val) | + | Pit[0] = 1.0 |
| - | PiStat = eigvecsl[:, max_val] | + | xvalues = [] |
| - | PiStat = PiStat.real | + | yvalues = [] |
| - | PiStatSum = PiStat.sum() | + | iter = 0 |
| - | PiStat0 = [] | + | while True: |
| - | PiStat1 = [] | + | iter += 1 |
| - | for k in range(4 * n): | + | Pit = Pit @ P |
| - | if k % 4 == 0: | + | TVD = sum(np.absolute(Pi - Pit) / 2.0) |
| - | PiStat0.append(4.0 * PiStat[k] / PiStatSum) | + | xvalues.append(iter / float(n)) |
| - | elif k % 4 == 2: | + | yvalues.append(TVD) |
| - | PiStat1.append(4.0 * PiStat[k] / PiStatSum) | + | if TVD < 0.1: break |
| - | XVec = [k for k in range(1, n + 1)] | + | pylab.plot(xvalues,yvalues, label='$n =$ '+str(n)) |
| - | print(PiStat0) | + | pylab.legend(loc='upper right') |
| - | pylab.plot(XVec, PiStat0) | + | pylab.xlabel("$t/ n$ (rescaled time) ") |
| - | pylab.plot(XVec, PiStat1) | + | pylab.ylabel("TVD") |
| - | pylab.xlabel("$x$ (position) ") | + | pylab.title("TVD lift tempering on the path graph of $n$ sites") |
| - | pylab.title("Simulated tempering (Lifted MCMC), Left EV transition matrix") | + | |
| pylab.show() | pylab.show() | ||
Current revision
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()
