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\):

  1. Compute energy at current density.

  2. Move randomly chosen agent randomly left or right.

  3. Compute second energy.

  4. Compare the two energies:

  5. If second energy is lower, accept move.

  6. \(\beta\) is a parameter which determines how likely the simulation is to move from a ‘less favourable’ situation to a ‘more favourable’ one.

  7. Compute \(P_0=e^{-\beta (E_1 - E_0)}\) and \(P_1\) a random number between 0 and 1,

  8. If \(P_0 > P_1\), do the move anyway.

  9. 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:

  1. A function to generate a random change: random_agent(), random_direction()

  2. A function to compute the energy before the change and after it: energy()

  3. A function to determine the probability of a change given the energy difference (1 if decreases, otherwise based on exponential): change_density()

  4. A function to determine whether to execute a change or not by drawing a random numberaccept_change()

  5. A method to iterate the above procedure: step()

Exercise: Fill in the gaps for testing

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

  2. We have provide the production code for change_density() and accept_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())