Recap example: Monte-Carlo
Contents
Recap example: Monte-Carlo#
Problem: Implement and test a simple Monte-Carlo algorithm#
Given an input function (energy) and starting point (density) and a temperature \(T\):
Compute energy at current density.
Move randomly chosen agent randomly left or right.
Compute second energy.
Compare the two energies:
If second energy is lower, accept move.
\(\beta\) is a parameter which determines how likely the simulation is to move from a ‘less favourable’ situation to a ‘more favourable’ one.
Compute \(P_0=e^{-\beta (E_1 - E_0)}\) and \(P_1\) a random number between 0 and 1,
If \(P_0 > P_1\), do the move anyway.
Repeat.
the algorithm should work for (m)any energy function(s).
there should be separate tests for separate steps! What constitutes a step?
tests for the Monte-Carlo should not depend on other parts of code.
Use matplotlib to plot density at each iteration, and make an animation
NB. If you are using the Windows command prompt, you will have to replace all %%bash directives in this notebook with %%cmd
Exercise - Partial Solution#
We have given you a partial solution below. In the solution we have broken the problem down into pieces:
A function to generate a random change:
random_agent(),random_direction()A function to compute the energy before the change and after it:
energy()A function to determine the probability of a change given the energy difference (1 if decreases, otherwise based on exponential):
change_density()A function to determine whether to execute a change or not by drawing a random number
accept_change()A method to iterate the above procedure:
step()
Exercise: Fill in the gaps for testing
Input insanity: We have provided some unittests can provide a range of possible valid and invalid inputs. Add the production code necessary to handle these senarios elegantly.
We have provide the production code for
change_density()andaccept_change(). Add tests to ensure these functions work as expected. Consider these senarios:
change_density(): density is change by a particle hopping left or right? Do all positions have an equal chance of moving?accept_change()will move be accepted when second energy is lower?
%%bash
rm -rf DiffusionExercise
mkdir DiffusionExercise
Windows: You will need to run the following instead
%%cmd
rmdir /s DiffusionExercise
mkdir DiffusionExercise
%%writefile DiffusionExercise/MonteCarlo.py
import matplotlib.pyplot as plt
from numpy import sum, array
from numpy.random import randint, choice
class MonteCarlo:
"""A simple Monte Carlo implementation"""
def __init__(self, energy, density, temperature=1, itermax=1000):
from numpy import any, array
# ADD CODE HERE
# Sanitise your inputs
density = array(density)
self.itermax = itermax
self.current_energy = energy(density)
self.temperature = temperature
self.density = density
def random_direction(self):
return choice([-1, 1])
def random_agent(self, density):
# Particle index
particle = randint(sum(density))
current = 0
for location, n in enumerate(density):
current += n
if current > particle:
break
return location
def change_density(self, density):
"""Move one particle left or right."""
location = self.random_agent(density)
# Move direction
if density[location] - 1 < 0:
return array(density)
if location == 0:
direction = 1
elif location == len(density) - 1:
direction = -1
else:
direction = self.random_direction()
# Now make change
result = array(density)
result[location] -= 1
result[location + direction] += 1
return result
def accept_change(self, prior, successor):
"""Returns true if should accept change."""
from numpy import exp
from numpy.random import uniform
if successor <= prior:
return True
else:
return exp(-(successor - prior) / self.temperature) > uniform()
def step(self):
iteration = 0
while iteration < self.itermax:
new_density = self.change_density(self.density)
new_energy = energy(new_density)
accept = self.accept_change(self.current_energy, new_energy)
if accept:
self.density, self.current_energy = new_density, new_energy
iteration += 1
return self.current_energy, self.density
def energy(density, coefficient=1):
"""Energy associated with the diffusion model
:Parameters:
density: array of positive integers
Number of particles at each position i in the array/geometry
"""
from numpy import array, any, sum
# Make sure input is an array
density = array(density)
# of the right kind (integer). Unless it is zero length, in which case type does not matter.
if density.dtype.kind != "i" and len(density) > 0:
raise TypeError("Density should be an array of *integers*.")
# and the right values (positive or null)
if any(density < 0):
raise ValueError("Density should be an array" + "of *positive* integers.")
if density.ndim != 1:
raise ValueError(
"Density should be an a *1-dimensional*" + "array of positive integers."
)
return coefficient * 0.5 * sum(density * (density - 1))
Writing DiffusionExercise/MonteCarlo.py
import sys
sys.path.append("DiffusionExercise")
import numpy as np
from IPython.display import HTML
from matplotlib import animation
from matplotlib import pyplot as plt
from MonteCarlo import MonteCarlo, energy
Temperature = 0.1
density = [np.sin(i) for i in np.linspace(0.1, 3, 100)]
density = np.array(density) * 100
density = density.astype(int)
fig = plt.figure()
ax = plt.axes(xlim=(-1, len(density)), ylim=(0, np.max(density) + 1))
image = ax.scatter(range(len(density)), density)
txt_energy = plt.text(0, 100, "Energy = 0")
plt.xlabel("Temperature = 0.1")
plt.ylabel("Energy Density")
mc = MonteCarlo(energy, density, temperature=Temperature)
def simulate(step):
energy, density = mc.step()
image.set_offsets(np.vstack((range(len(density)), density)).T)
txt_energy.set_text(f"Energy = {energy}")
anim = animation.FuncAnimation(fig, simulate, frames=200, interval=50)
HTML(anim.to_jshtml())