{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Recap example: Monte-Carlo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem: Implement and test a simple Monte-Carlo algorithm\n", "\n", "Given an input function (energy) and starting point (density) and a temperature $T$: \n", "\n", "1. Compute energy at current density.\n", "1. Move randomly chosen agent randomly left or right.\n", "1. Compute second energy.\n", "1. Compare the two energies:\n", "1. If second energy is lower, accept move.\n", "1. $\\beta$ is a parameter which determines how likely\n", " the simulation is to move from a 'less favourable' situation to a 'more favourable' one.\n", "1. Compute $P_0=e^{-\\beta (E_1 - E_0)}$ and $P_1$ a random number between 0 and 1,\n", "1. If $P_0 > P_1$, do the move anyway.\n", "1. Repeat." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* the algorithm should work for (m)any energy function(s).\n", "* there should be separate tests for separate steps! What constitutes a step?\n", "* tests for the Monte-Carlo should not depend on other parts of code.\n", "* Use [matplotlib](http://matplotlib.org/) to plot density at each iteration, and make an animation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NB.** If you are using the Windows command prompt, you will have to replace all %%bash directives in this notebook with %%cmd " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise - Partial Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have given you a partial solution below. In the solution we have broken the problem down into pieces:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. A function to generate a random change: random_agent(), random_direction()\n", "1. A function to compute the energy before the change and after it: energy()\n", "1. A function to determine the probability of a change given the energy difference (1 if decreases, otherwise based on exponential): change_density()\n", "1. A function to determine whether to execute a change or not by drawing a random numberaccept_change()\n", "1. A method to iterate the above procedure: step()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**_Exercise_: Fill in the gaps for testing**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "1. We have provide the production code for change_density() and accept_change(). Add tests to ensure these functions work as expected. Consider these senarios:\n", " - change_density(): density is change by a particle hopping left or right? Do all positions have an equal chance of moving?\n", " - accept_change() will move be accepted when second energy is lower?" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "rm -rf DiffusionExercise\n", "mkdir DiffusionExercise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Windows:** You will need to run the following instead\n", " \n", "cmd\n", "%%cmd\n", "rmdir /s DiffusionExercise\n", "mkdir DiffusionExercise\n", "" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing DiffusionExercise/MonteCarlo.py\n" ] } ], "source": [ "%%writefile DiffusionExercise/MonteCarlo.py\n", "import matplotlib.pyplot as plt\n", "from numpy import sum, array\n", "from numpy.random import randint, choice\n", "\n", "\n", "class MonteCarlo:\n", " \"\"\"A simple Monte Carlo implementation\"\"\"\n", "\n", " def __init__(self, energy, density, temperature=1, itermax=1000):\n", " from numpy import any, array\n", "\n", " # ADD CODE HERE\n", " # Sanitise your inputs\n", " \n", " density = array(density)\n", " self.itermax = itermax\n", "\n", " self.current_energy = energy(density)\n", " self.temperature = temperature\n", " self.density = density\n", "\n", " def random_direction(self):\n", " return choice([-1, 1])\n", "\n", " def random_agent(self, density):\n", " # Particle index\n", " particle = randint(sum(density))\n", " current = 0\n", " for location, n in enumerate(density):\n", " current += n\n", " if current > particle:\n", " break\n", " return location\n", "\n", " def change_density(self, density):\n", " \"\"\"Move one particle left or right.\"\"\"\n", "\n", " location = self.random_agent(density)\n", "\n", " # Move direction\n", " if density[location] - 1 < 0:\n", " return array(density)\n", " if location == 0:\n", " direction = 1\n", " elif location == len(density) - 1:\n", " direction = -1\n", " else:\n", " direction = self.random_direction()\n", "\n", " # Now make change\n", " result = array(density)\n", " result[location] -= 1\n", " result[location + direction] += 1\n", " return result\n", "\n", " def accept_change(self, prior, successor):\n", " \"\"\"Returns true if should accept change.\"\"\"\n", " from numpy import exp\n", " from numpy.random import uniform\n", "\n", " if successor <= prior:\n", " return True\n", " else:\n", " return exp(-(successor - prior) / self.temperature) > uniform()\n", "\n", " def step(self):\n", " iteration = 0\n", " while iteration < self.itermax:\n", " new_density = self.change_density(self.density)\n", " new_energy = energy(new_density)\n", "\n", " accept = self.accept_change(self.current_energy, new_energy)\n", " if accept:\n", " self.density, self.current_energy = new_density, new_energy\n", " iteration += 1\n", "\n", " return self.current_energy, self.density\n", "\n", "\n", "def energy(density, coefficient=1):\n", " \"\"\"Energy associated with the diffusion model\n", " :Parameters:\n", " density: array of positive integers\n", " Number of particles at each position i in the array/geometry\n", " \"\"\"\n", " from numpy import array, any, sum\n", "\n", " # Make sure input is an array\n", " density = array(density)\n", "\n", " # of the right kind (integer). Unless it is zero length, in which case type does not matter.\n", " if density.dtype.kind != \"i\" and len(density) > 0:\n", " raise TypeError(\"Density should be an array of *integers*.\")\n", " # and the right values (positive or null)\n", " if any(density < 0):\n", " raise ValueError(\"Density should be an array\" + \"of *positive* integers.\")\n", " if density.ndim != 1:\n", " raise ValueError(\n", " \"Density should be an a *1-dimensional*\" + \"array of positive integers.\"\n", " )\n", "\n", " return coefficient * 0.5 * sum(density * (density - 1))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "