06_cell_population_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cell Division Code Example\n",
    "\n",
    "In this notebook we're going to look through a simple but complete example of how we might write a parallel piece of code to simulate the population of cells over time.\n",
    "\n",
    "## The Problem\n",
    "\n",
    "We are going to simulate a population of cells over time. The cells have an average lifetime, but the actual time a cell will live is random. If $l$ is the lifetime of a cell, then the probability density function that a cell will die at time $t$ after its birth is given by: \n",
    "\n",
    "$P(t) = \\textrm{e}^{-\\frac{t}{l}}$.\n",
    "\n",
    "When a cell will ends its life, it will either die or split into two cells. In our simulation, we'll use a mean lifetime of 1, and we'll assume the cell has a probability of splitting of 0.55 and a probability of dying of 0.45. As the time of death of each cell and the fate of each cell is non-deterministic, the same starting conditions will lead to different outcomes each time the simulation is run.\n",
    "\n",
    "As cells are more likely to split than to die, on average, the population will grow as in the figure below:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/single_growing_population.png\" alt=\"A figure a realisation when the population grows.\" class=\"center\">\n",
    "</p>\n",
    "\n",
    "However, there is a significant chance that the population will die out, potentially very quickly:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/single_quick_death.png\" alt=\"A figure a realisation when the population quickly dies.\" class=\"center\">\n",
    "</p>\n",
    "\n",
    "We would like to find out how the mean and standard deviation of the population changes over time. We would also like to calculate the survival probability, which is the probability that the population is non-zero at a given time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Serial Simulation\n",
    "\n",
    "We can simulate this problem by running a number of different realisations of the simulation. Each of these realisations will be an independent simulation of a possible history of the population of cells as a function of time. Once we have run many realisations, we can use the collection of realisations to calculate the mean and standard deviation of the population at each time, as well as the survival probability.\n",
    "\n",
    "A serial realisation of the simulation can be found in [`06_cell_population_example/serial_version.py`](06_cell_population_example/serial_version.py). We don't need to look into the details of most of the functions, but it's worth understanding what each function does:\n",
    "\n",
    "- `generate_next_population`: This function calculates the number of cells that will be present at a certain time after the current time for a single cell. This uses a random number generator to determine if a cell will die, split, or survive the time period.\n",
    "- `run_realisation`: This function runs a single realisation of the simulation. It starts with a single cell at time 0 and calculates the population at each of a number of output times using `generate_next_population`.\n",
    "- `run_single_realisation`: This function runs a single realisation of the system using `run_realisation` and then calculates the mean and plots the output.\n",
    "- `run_multiple_realisations`: This function runs several realisations of the system using `run_single_realisation` and then calculates the mean and standard deviation of the population at each time, as well as the survival probability.\n",
    "\n",
    "These functions contain the ability to receive a seed to set up the random number generator. This is useful for testing and comparing performance as it allows us to guarantee the output of the simulation.\n",
    "\n",
    "The script contain code which runs two single realisations - one which always dies out and one which always grows (the results are guaranteed by manually choosing the seed of the random number generator). The script then runs 200 realisations of the system. Each realisation receives the number of the realisation (from 0-199) as a seed to ensure reproducibility. The runtimes of the system are printed to the console. When I ran the code, the output was:\n",
    "\n",
    "| Simulation Type                        | Running Time                    | Plotting Time |\n",
    "|----------------------------------------|---------------------------------|---------------|\n",
    "| Single realisation that always dies    | 6.2 $\\times 10 ^{-5}\\textrm{s}$ | 0.31s         |\n",
    "| Single realisation that always grows   | 2.0s                            | 0.20          |\n",
    "| 200 realisations                       | 189s                            | 0.50s         |\n",
    "\n",
    "The time to simulate a growing population is significantly longer than the time to simulate a dying population as there are more cells to simulate. Around 20% of the simulations see a growing population. For a realisation which grows, different realisations may take different amounts of time, depending on how large the population gets. The figure below shows the amount of time taken for each realisation to run, with the number of the longest-lasted realisations added as annotation to the figures:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/serial_runtimes.png\" alt=\"A figure showing the time taken for each realisation to run in serial.\" class=\"center\">\n",
    "</p>\n",
    "\n",
    "The output of the code for 200 realisations is shown below:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/multiple_realisations.png\" alt=\"A figure showing the output of the serial simulation.\" class=\"center\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Queue Implementation\n",
    "\n",
    "This problem is a god candidate for parallelisation as each realisation is independent of the others. As each realisation can take a while to run, it also means the overhead due to communication is likely to be small compared to the time taken to run the realisation.\n",
    "\n",
    "As we parallelise the code, we want to keep the interface for the functions a user might call as similar as possible, specifically, `run_single_realisation` and `run_multiple_realisations`. This means it will take minimal effort adapt existing tests and profiling, and any users running the code, or any places where the code is called in existing projects will not need to be changed.\n",
    "\n",
    "Our first attempt to parallelising the code is to use a queue to store te results produced from a number of realisations in the file [`06_cell_population_example/queue_version.py`](06_cell_population_example/queue_version.py). To do this we create the new function `run_n_realisation_queue` which is similar to the old function `run_multiple_realisations` but uses a queue to store the results of all realisations performed in a 2D Numpy array. This function will be called by each process. The function `run_multiple_realisations` is adapted to create the queue, start the processes, collect the results from the queue, and process the results. Each process returns a 2D Numpy array with the population at each time for each realisation.\n",
    "\n",
    "When altering `run_multiple_realisations` we have made the number of processes an optional argument with a default value of 1. This means that calls made to the function without specifying the number of processes will still work, making integration of the new function into existing projects easier.\n",
    "\n",
    "This implementation doesn't alter the runtime of the single realisations, but decreases the runtime from around 189s to around 107s on 4 cores. This is a decent speedup, but the code is not 4 times faster. Part of the reason for this becomes apparent when we run the code and view which process is working on each realisation, as in the figure below:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/queue_runtimes.png\" alt=\"The amount of time each process spent running realisations using a Queue.\" class=\"center\">\n",
    "</p>\n",
    "\n",
    "In the example above, Process 4 was running realisations 150-199, which happened to not include many realisations where the population grew and so finished in 26s. Processes 2 and 3 had more long-lived realisations and took 76s and 80s to run respectively. Process 1 happened to have several long-lived realisations and took 107s to run.\n",
    "\n",
    "Once a process has finished its realisations it will terminate and the physical core will be inactive. The code is left waiting for the slowest process to finish, meaning progressively fewer of the cores are active as the code runs. This is a common problem when parallelising code, and is known as load imbalance. Ideally, we would like a way to keep our processes busy for more of the time to make the overall calculation finish faster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pool Implementation\n",
    "\n",
    "To solve the problem of load imbalance, we can use a `Pool`. The advantage of a pool is that it will keep all the processes busy by assigning them new tasks as they finish their previous task. This means processes will be kept busy for more of the time, and the overall calculation will finish faster. This is implemented in the file [`06_cell_population_example/pool_version.py`](06_cell_population_example/pool_version.py).\n",
    "\n",
    "This version is arguably simpler than the queue version as we don't need to have a function like `run_n_realisation_queue` to use as an interface between `run_multiple_realisations` and `run_single_realisation`. Instead, we can use the `starmap` function from the `Pool` object to run the realisations. Once we receive the results from the `Pool` object, we can process them into a 2D Numpy array and process them as before.\n",
    "\n",
    "The figure below shows the amount of time each process spent performing each realisation:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/pool_runtimes.png\" alt=\"The amount of time each process spent running realisations using a Pool.\" class=\"center\">\n",
    "\n",
    "Primarily because of the way the `Pool` distributes work to the processes, the load is now much more evenly balanced between the processes. The code now takes around 83s to run on 4 cores."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deterministic Approximation\n",
    "\n",
    "We've seen how we can parallelise the code and this significantly speeds it up. But we can speed up a single realisation. We can note that the runtime of the code is strongly related to the number of cells simulated and most of the cells in a realisation are simulated at the end of a growing simulation when the population is high. If we look back at the figure at the start of the notebook showing how the population grows over time in a single realisation, we also note that the growth in population is smooth and deterministic once the population exceeds about 10,000 neutrons. This means we can use a deterministic approximation to calculate the population at these times, rather than simulating each cell individually.\n",
    "\n",
    "When the population is large, we can use the following approximation:\n",
    "\n",
    "$$\n",
    "\\frac{dN(t)}{dt} = \\frac{2 N(t)(p_{r} - 0.5)}{l}\n",
    "$$\n",
    "\n",
    "where $N(t)$ is the population at time $t$, $p_{r}$ is the probability of a cell splitting, and $l$ is the mean lifetime of a cell. This is a simple differential equation which can be solved to give:\n",
    "\n",
    "$$\n",
    "N(t) = N(t_{0}) \\textrm{e}^{\\frac{2(p_{r} - 0.5)(t-t_{0})}{l}}\n",
    "$$\n",
    "\n",
    "where $t_{0}$ is the time at which the population was $N(t_{0})$. In our case, we can switch to this approximation when the population exceeds 10,000 cells and use this time as $t_{0}$.\n",
    "\n",
    "The file [`06_cell_population_example/deterministic_version.py`](06_cell_population_example/deterministic_version.py) contains the implementation of this approximation. The function `run_realisation` has been adapted so the loop over the different times exits when the population exceeds 10,000 cells. When this happens a simple Numpy calculation is used to calculate the population of the remaining times. \n",
    "\n",
    "<p align=\"center\">\n",
    "    <img src=\"resources/single_growing_population.png\" alt=\"A single growing realisation using the initial method\" style=\"display: inline-block; margin-right: 10px; width: 45%;\">\n",
    "    <img src=\"resources/single_growing_deterministic.png\" alt=\"A single growing realisation using the deterministic approximation\" style=\"display: inline-block; width: 45%;\">\n",
    "</p>\n",
    "\n",
    "The figures below show the results of 200 realisation for the original method (left) and deterministic approximation (right):\n",
    "\n",
    "<p align=\"center\">\n",
    "    <img src=\"resources/multiple_realisations.png\" alt=\"Results of 200 realisations using the initial method\" style=\"display: inline-block; margin-right: 10px; width: 45%;\">\n",
    "    <img src=\"resources/multiple_deterministic.png\" alt=\"Results of 200 realisations using the deterministic approximation\" style=\"display: inline-block; width: 45%;\">\n",
    "</p>\n",
    "\n",
    "The results are indistinguishable. The runtime of the code is reduced to 5.1s for 200 realisations. The figure below shows the activity of each process:\n",
    "\n",
    "<p align=\"center\">\n",
    "<img src=\"resources/deterministic_runtimes.png\" alt=\"The amount of time each process spent running realisations using the deterministic approximation.\" class=\"center\">\n",
    "</p>\n",
    "\n",
    "This optimisation has significantly reduced the runtime of the code. The optimisation didn't relate to the parallelisation of the code, but its effectiveness highlights how parallelisation is only one tool we can use to speed up our code and it may not always be the most effective approach. When attempting to speed up a code, it's worth considering a variety of approaches."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}