%matplotlib inline

import numpy as np, matplotlib as mpl, scipy as sp, networkx as nx
from matplotlib import pylab, mlab, pyplot as plt
from matplotlib.pylab import plot, scatter, contour, xlabel, ylabel, title, legend
from matplotlib.animation import FuncAnimation
from tqdm.notebook import tqdm, trange

from numpy import sqrt, pi, exp, log, floor, ceil, sin, cos
from numpy import linspace, arange, empty, zeros, ones, empty_like, zeros_like, ones_like
from numpy.linalg import norm

plt.rcParams.update({
    'image.cmap': 'coolwarm',
    'animation.html': 'jshtml',
    'animation.embed_limit': 40, # 40 mb
})
rng = np.random.default_rng()

# Choose a uniform grid as our Graph.
(N1, N2) = (10, 10)
β = 10

G = nx.grid_2d_graph(N1, N2)
pos = {(x, y): (y, -x) for x, y in G.nodes()}

# Number of realization of the Markov chain we will simulate
n_reals = 1000

x0 = rng.choice([-1, 1], size=(N1, N2, n_reals))
nx.draw(G, pos=pos, node_color=x0[:, :, 0])
_ = title('Uniformly random spin configuration')

Our goal is to simulate n_reals independent realizations of a MCMC algorithm sampling from the Gibbs distribution in the Ising model. In order to test our simulations, we also need to compute $\ln \pi_u$ in a given state, which we do with the function logπu.

The functions mcmc and logπu have been stripped from this notebook. In order for the notebook to run, you will have to implement these functions yourself. Here are the proto-types and doc-strings:

def mcmc(x0, n_iters, β):
    """Perform an MCMC iteration

    x0 -- initial spin configuartion (N1 × N2 × n_reals array)
    n_iters -- number of iterations to perform
    β -- inverse temperature

    Returns: an N1×N2×n_reals array of spin configurations after iteration
    """

def logπu(x, β):
    """Returns the log of the un-normalized Gibbs measure with inverse temperature β"""

Implementation notes and tips:

1. Do not compute $\pi_u$, and then take it's logarithm for the function logπu. Doing this will get you overflow errors.

2. DO NOT use the function logπu when computing the acceptance ratio in your implementation of mcmc. Since the proposed state and current state differ in only one spin, you can compute the acceptance ratio by examining the spin at just a few points. Do a pencil and paper calculation and figure out the exact formula. (If you use the function logπu, then your code will be about 100 times slower (when $N = 10$). It gets worse as $N$ increases.)

3. Feel free to use any algorithm you can think of: Metropolis Hastings with the proposal suggested in class; with a different proposal, Glauber dynamics, etc. To get a better feel for such MCMC simulations you may want to experiment with the proposal chain a little and see if something gives you better results.

Few things that are helpful during implementation:

1. v = tuple(rng.choice(G.nodes())) Chooses one random node of G

2. G.neighbors(v) generates an iterable of neighbors of the node v

Performance notes:

You don't have to optimize your code now (the numbers are small enough that it will still run within minutes with bad code). But if you're thinking about how to make your code faster here are a few tips:

1. Vectorize your code; don't run a for loop overall realizations -- that's really slow. Instead use the vectorized numpy functions to operate on all realizations using an array.

2. To be fast (run in a few seconds), you can "cheat" a tiny bit: The MCMC algorithms we've seen so far involve choosing a vertex $v$ uniformly at random, and then changing the spin at $v$ in some manner. To get n_reals independent runs of your Markov chain, you need to choose this vertex at random indpendently for every realization. But doing this makes the code hard to vectorize. Instead, if you choose the same vertex for all realizations, you can vectorize the code and get somethign that runs in a few seconds. It's not going to produce independent identically distributed runs of the Markov Chain; but it works well enough in practice.

Once you have implemented these functions, the rest of this notebook should run and produce the figures shown on the main page.

β = 1
y = x0.copy()
lπu = [logπu(y, β)]

iters_per_step = 10
for _ in trange(500):
    y = mcmc(y, iters_per_step, β)
    lπu.append(logπu(y, β))

lπu = np.array(lπu)

  0%|          | 0/500 [00:00<?, ?it/s]

iters = arange(lπu.shape[0]) * iters_per_step

μ = lπu.mean(axis=-1)
plot(iters, μ, label=r'$E\log(\pi_u)$')

σ = lπu.std(axis=-1)
plt.fill_between(iters, μ - σ, μ + σ, alpha=.3, label='StdDev')

for n in sorted(rng.integers(n_reals, size=3)):
    plot(iters, lπu[:, rng.integers(n_reals)], label=f'Trial {n}', alpha=.7)

lπ1 = logπu(ones((N1, N1)), β)
plt.hlines(lπ1, iters[0], iters[-1], linestyles='--', label=r'Max $\log \pi_u$', alpha=.5)

xlabel('# Iterations')
legend()
# plt.savefig( 'figures/ising-logpiu.svg' )

<matplotlib.legend.Legend at 0x7d9a79799940>

def M(y:np.ndarray):
    """Average magnetization"""
    return y.mean(axis=(0, 1))


iters_per_step = 100
ββ = np.logspace(-2, 0, num=7)
mag = np.empty((100, len(ββ)))

for i, β in tqdm(enumerate(ββ), total=len(ββ)):
    y = x0.copy()
    mag[0, i] = abs(M(y)).mean()
    for t in trange(1, mag.shape[0], leave=False):
        y = mcmc(y, iters_per_step, β)
        mag[t, i] = abs(M(y)).mean()

  0%|          | 0/7 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

  0%|          | 0/99 [00:00<?, ?it/s]

tt = np.arange( mag.shape[0] ) * iters_per_step
for i, β in enumerate( ββ ):
    plot( tt, mag[:, i], label=f'β={β:.2f}' )

plt.ylabel( 'Expected magnetization strength' )
plt.xlabel( '# MCMC iterations' )
plt.legend()

#plt.savefig( 'figures/ising-mag-transition.svg' )

<matplotlib.legend.Legend at 0x7d9a782e9810>

β=.5
iters_per_step = 10

y = np.zeros((N1, N2, 100))
yn = x0[:, :, 0:1].copy()
y[:, :, 0] = yn[:, :, 0]
for n in trange(1, y.shape[-1]):
    yn = mcmc(yn, iters_per_step, β)
    y[:, :, n] = yn[:, :, 0]

  0%|          | 0/99 [00:00<?, ?it/s]

n = -1
nx.draw(G, pos=pos, node_color=y[:, :, n])

fig, ax = plt.subplots()
def update(frame):
    ax.clear()
    ax.set_title(f'Iteration {frame*iters_per_step}')
    nx.draw(G, pos=pos, node_color=y[:, :, frame])


ani = FuncAnimation(fig, update, frames=y.shape[-1])
ani

Once Loop Reflect

# with open( 'figures/mcmc-ising.html.inc', 'w' ) as f: f.write(  ani.to_jshtml() )