In [1]:
%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()
In [2]:
(N1, N2) = (10, 10)
β = 10

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

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' )
No description has been provided for this image

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 β"""

Few notes:

  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.)

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

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

In [7]:
β = 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]
In [8]:
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' )
Out[8]:
<matplotlib.legend.Legend at 0x799259b6dd00>
No description has been provided for this image
In [9]:
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]
In [10]:
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' )
Out[10]:
<matplotlib.legend.Legend at 0x799259a47f50>
No description has been provided for this image
In [11]:
β=.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]
In [12]:
n = -1
nx.draw( G, pos=pos, node_color=y[:, :, n] )
No description has been provided for this image
In [13]:
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
Out[13]:
No description has been provided for this image
No description has been provided for this image
In [14]:
#with open( 'figures/mcmc-ising.html.inc', 'w' ) as f: f.write(  ani.to_jshtml() )
In [ ]: