Lesson_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Monte Carlo</h2> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes physical systems consist of many parameters that are not well known, or that evolve over time in unpredictable ways.  Magnetic spins, diffuse light scattering, and electron gases are all examples of systems that can better be modeled using Monte Carlo methods, so named because the randomness involved <a href=\"https://en.wikipedia.org/wiki/Monte_Carlo_method#History\">reminded Nicholas Metropolis</a> of gambling at the famous casino in Monte Carlo, in Monaco. (Metropolis and his colleagues were working at Los Alamos shortly after the end of World War II, and modeling neutron transport through materials.) \n",
    "\n",
    "In lab today we will explore these methods by modeling how light scatters through biological tissue. As you might imagine, understanding this is important for understanding different kinds of medical imaging.  Light scattering is a natural use of Monte Carlo methods, because in many cases, when an atom or molecule absorbs a photon, a short time later it re-emits a photon in a random direction, retaining no memory of the initial direction of the photon.  We can model this behavior by generating random numbers and using them to specify the direction and distance that the photon travels.  In some cases, there *is* some memory of the photon's initial direction (e.g. it is more likely to scatter forward than backward).  As we will see, we can also model such cases by specifying a non-uniform distribution (with whatever functional form suits the particular situation) for the new direction; this is called the *phase function* of the scattering. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random numbers\n",
    "\n",
    "Some physical processes are random (e.g. radioactive decay), and early calculations with Monte Carlo simulations used tables of these truly random numbers.  However, using tabulated random numbers is tedious, and it's useful if we can generate new ones with a computer. \n",
    "\n",
    "However, nothing a computer does is truly random, so these are more technically termed \"pseudo-random\" numbers.  For most purposes these are plenty random enough for anything we want to do, but if your analysis depends heavily on them, it is good to be aware of the strength (i.e. randomness) of the random number generator in the language you're using.  In particular, some algorithms will always give you the same pseudo-random sequence unless you are careful to alter the \"seed\" given to the algorithm.  (This isn't a problem with the Python routines we'll be using today.) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1\n",
    "To explore Python's random number generation, use `numpy.random.random()` to generate 5000 `x` values and 5000 `y` values.  (As with other specific routines we'll mention here, you can find the documentation for this routine with a little web searching, to understand what it expects as inputs and what it returns.)  With these values:\n",
    "\n",
    "* Plot them in a labeled x-y plot (as points, not lines!).  They should look randomly distributed. \n",
    "* Plot a histogram of their values (using `hist`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 2\n",
    "\n",
    "Using your `x` and `y` values from above, do a simple approximation of $\\pi$, as follows.  Select only those points that lie at a distance less than or equal to 1 away from the origin, i.e. that lie in the first quadrant of the unit circle.  With a little thought, you should be able to figure out how to use the number of values you find to make a calculation of the value of $\\pi$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a little hint about a `numpy` feature that will be useful here, and below.  Rather than looping over a numpy array one element at a time to find elements with certain values, you can do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This lets the code work in Python 2 or 3:\n",
    "from __future__ import division, print_function\n",
    "\n",
    "# Make a list of integers for our example:\n",
    "values = np.arange(10)\n",
    "\n",
    "# Assign a new variable to the output of a logical test:\n",
    "less_than_five = (values < 5)\n",
    "\n",
    "# Note that this produces an array of boolean values:\n",
    "print(\"Contents of array 'less_than_five':\")\n",
    "print(less_than_five)\n",
    "\n",
    "# Now we can use this boolean array as an *index*, to \n",
    "# pull out only the desired elements:\n",
    "new_values = values[less_than_five]\n",
    "print()  # Just a blank line.\n",
    "print(\"Result of array 'values' indexed by 'less_than_five' array:\")\n",
    "print(new_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now do the above exercise for simulating $\\pi$, using the indexing trick just outlined.  If your value for $\\pi$ is relatively inaccurate, increase the number of random values you use until you reliably get it correct to at least two or three decimal places.  Don't worry - you'll find that generating even millions of random numbers is relatively quick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications to light scattering\n",
    "\n",
    "OK, now let's do something more interesting than calculating $\\pi$!  We'll take an example of a light source (e.g. a laser) being shone through some tissue.  If there were no scattering or absorption, the image of the light coming out would be the same as it was going in, i.e. a relatively compact point source.  (\"No scattering or absorption\" is another way of saying \"transparent\", so just think of shining a light through a pane of glass.)  However, if there *is* some scattering, the image will change. (Think about the same example, but with a frosted pane of glass rather than a clear one.)  Being able to model how the image changes as it passes through tissue helps diagnosticians understand what the tissue properties are, since, e.g., a tumor might have different light scattering or absorption properties than healthy tissue. \n",
    "\n",
    "Read the first four pages of the PDF posted in Moodle.  We will be following the mathematical framework spelled out there to model light scattering through a slab of tissue.  You can stop when you get to the part labeled \"Lab Assignment\".  We'll be doing exercises similar to those, but breaking it down a bit more into steps.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 3\n",
    "\n",
    "Look at equation 1 in the background document, and use it to generate a `numpy` array of step sizes.  Write it as a function that takes as input the number of steps $n$, and the desired mean step size, and returns a `numpy` array of steps.  Hint: `np.log` gives the natural logarithm, as desired in this equation; the base-10 log is written `np.log10`.  \n",
    "\n",
    "After writing your function, write code to call your function to generate 5000 steps with a mean step size of 0.01 (a typical mean free path for a visible-light photon in tissue), and plot a histogram of the step values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 4\n",
    "\n",
    "Implement equations 5 and 6 as two functions, to return values of $\\theta$ and $\\phi$.  Your function for $\\theta$ should take as inputs the value of $g$, the anisotropy parameter, and $n$, the number of steps desired, and return an $n$-element `numpy` array of $\\theta$ values.  Alternatively, it's fine for your function to return $\\cos \\theta$ rather than $\\theta$ - in fact that's more efficient, since you need $\\cos \\theta$ in equations 7-9 anyway.  Just be sure to keep track of which one you have!  Also write a separate function that returns $n$ values of $\\phi$.  (Note that it's fine for $n$ to be 1; a `numpy` array of length 1 can be treated as a scalar in your code with no problems, i.e. any `numpy` operations work fine on scalars as well as vectors.)\n",
    "\n",
    "After writing your function, test it by plotting histograms of the $\\theta$ and $\\phi$ values generated.  For $\\theta$, plot both the $g=0$ and $g=0.9$ cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 5\n",
    "\n",
    "As the final piece of the puzzle before doing the actual scattering calculation, write a function that implements equations 7-9 (as well as checking for the special case of equations 10-12) to calculate the new direction cosines based on the old ones, given the new scattering angles.  To get a function to return three values (your three direction cosines), you just give the function's `return` statement three values, e.g. \n",
    "\n",
    "`return mu_x, mu_y, mu_z`\n",
    "\n",
    "You can capture these outputs into separate variables when calling the function if you assign the function's output to three variables on the calling end, e.g. \n",
    "\n",
    "`a, b, c = myfunc(x)`\n",
    "\n",
    "After writing your functions, show that they return the expected values for $\\theta = 0, \\phi = 0$; for $\\theta = \\pi/2, \\phi = 0$; and that both of those cases work when $\\mu_z = 1$ and when $\\mu_z$ is not close to 1.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 6\n",
    "\n",
    "Now we're ready to actually track some photons!  Write code that tracks $n$ photons through a tissue sample that is 0.1 cm thick.  Follow the flowchart in Figure 1 and the explanation there to step a photon through the sample; for each photon, you'll need to follow it until it exits the sample.  (Assume that the sample is 2 x 2 cm in $x$ and $y$; if either of those coordinates of a photon gets more than 1 cm from the origin before it reaches $z = z_{max}$, then assume that that photon is lost and skip to the next one.  Otherwise, your code will spend a lot of time chasing the handful of photons that happen to wander off sideways.)  \n",
    "\n",
    "When each photon reaches $z_{max}$, store its $x$ and $y$ position in lists; we'll use these to create an image of what the output light beam looks like when it exits the tissue. \n",
    "\n",
    "After you have run all $n$ photons and stored their final positions, make an $x-y$ plot.  You can make a simple plot with just `mp.plot(x_final, y_final, 'g.')`.\n",
    "\n",
    "Start by running just a small number of photons (20 or so), and experiment with different values of $g$.  When you have code that is working well, you can gradually increase the number of photons to get better-looking images.  How many you use will depend on how fast your code is, and how patient you are!  Display at least two final images, for $g=0.1$ and $g=0.9$.  Comment on why they look the way they do.  That is, explain why you would expect different values of $g$ to influence the images in the way that is shown.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that when you have a lot of photons, they start to overlap in the plot and it can be harder to tell how concentrated they are.  You can make better-looking plots by generating a density plot of the final $x, y$ values and then displaying it as an image.  Here, the intensity of each pixel is proportional to how many photons landed there.\n",
    "\n",
    "First, you'll need to install a new python package.  Open a command line terminal and type \"pip install seaborn\".  Pip is python's package manager.  There are many, many packages available out there.  Seaborn is one particular package useful for creating nice looking graphs.\n",
    "\n",
    "Here's some example code, plotting exponentially decaying distributions in $x$ and $y$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import numpy.random as rng\n",
    "import seaborn\n",
    "\n",
    "n = 1000 #number of points\n",
    "\n",
    "#creating linear distributions\n",
    "a = np.zeros(n)\n",
    "b = np.zeros(n)\n",
    "\n",
    "for r in range(n):\n",
    "    a[r] = rng.random()\n",
    "    b[r] = rng.random()\n",
    "\n",
    "#creating exponential distributions\n",
    "x = np.zeros(n)\n",
    "y = np.zeros(n)\n",
    "for r in range(n):\n",
    "    x[r] = -np.log(a[r])\n",
    "    y[r] = -np.log(b[r])\n",
    "\n",
    "#plotting the data\n",
    "ax = seaborn.kdeplot(x, y, shade=True, cmap='gray', nlevels=60, clip=([0,2],[0,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "If you have time, follow the above example to make some nicer images of your scattered photons for different values of $g$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py36]",
   "language": "python",
   "name": "Python [py36]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}