TVDTemperingLift.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 21:56, 16 March 2022 Werner (Talk | contribs) ← Previous diff |
Revision as of 10:26, 8 September 2022 Werner (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()
