Optimizing Python code - exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "source": [
    "# A tour of python profiling and optimization (Part I)\n",
    "\n",
    "*Antonino Ingargiola*, EuroScipy 2018\n",
    "\n",
    "\n",
    "> premature optimization is the root of all evil (or at least most of it) in programming. \n",
    "\n",
    "                                            - Donald Knuth\n",
    "                                                         \n",
    "# Overview\n",
    "\n",
    "In this notebooks I will show how to optimize python code\n",
    "using the *Ornstein-Uhlenbeck process* as an example.\n",
    "We will cover profiling (`%timeit`, `%prun`, `%line_profiler`), \n",
    "loop optimization, Numpy, Numba, Cython,\n",
    "random number generation (RNG) using `randomgen` and packaging cython code.\n",
    "\n",
    "Parallelism is covered in the companion notebook: [Parallelism with Joblib](Parallelism%20with%20Joblib.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OU random process for hackers\n",
    "\n",
    "*Reference Wikipedia: [Ornstein-Uhlenbeck process](https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process)*\n",
    "\n",
    "I am going to introduce a simple random process, called the Ornestein-Uhlenbeck (OU) process. \n",
    "\n",
    "I will give a direct definition in form of 🐍python code🐍 (no math).\n",
    "Then, we will try to optimize it.\n",
    "\n",
    "## Definition\n",
    "\n",
    "The OU process is a random function of \"time\" $f(t)$, with properties:\n",
    "\n",
    "1. it is a Normal random number at each point\n",
    "2. correlation with a previous sample vanishes exponentially\n",
    "\n",
    "## Applications\n",
    "\n",
    "- Diffusion of particles/residues in a potential well\n",
    "- *Vasicek* model for Interest rates and other finance applications\n",
    "\n",
    "> OU is the simplest Gaussian process with correlation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: white noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randn(1000)\n",
    "\n",
    "plt.plot(x)\n",
    "plt.figure()\n",
    "plt.hist(x, 30);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`x` is an OU process with $\\sigma = 1$ and $\\tau = 0$ (memory is lost immediately)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "\n",
    "\n",
    "### Process parameters\n",
    "\n",
    "These 2 parameters completely define an OU process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "σ = 2  # standard deviation of the Normal distribution\n",
    "τ = 1  # relaxation time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation parameters\n",
    "\n",
    "These parameters define \"what we compute\" about the OU process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_steps = 100  # number of time points to simulate\n",
    "δt = 0.1         # time interval between points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OU process implementation\n",
    "\n",
    "*Reference: [Gillespie's explicit update formula](http://10.1103/PhysRevE.54.2084)*\n",
    "\n",
    "If we know the value `x0` of the process at time $t$,\n",
    "we can compute the value at time $t + \\delta t$ **in one step**:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "x0 = 1    # value of the process at time t\n",
    "N_norm = np.random.randn()  # a Normal random number\n",
    "\n",
    "# Gillespie's explicit update formula:\n",
    "relax = np.exp(-δt / τ)\n",
    "diffuse = σ * np.sqrt(1 - relax**2)\n",
    "x1 = x0 * relax + diffuse * N_norm\n",
    "x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def OU_update0(*, x_prev, δt, τ, σ):\n",
    "    '''Compute the OU process having initial value `x_prev` after a time `δt`.\n",
    "    \n",
    "    Knowing the value `x_prev` of an OU process at time \"t\",\n",
    "    this function returns the value at time \"t + δt\".\n",
    "    \n",
    "    Arguments:\n",
    "        x_prev (float): initial value of the OU process (time \"t\")\n",
    "        δt (float): time increment after which to evaluate the OU process\n",
    "        τ (float): relaxation time, i.e. time after which correlation is lost\n",
    "        σ (float): standard deviation each point of the OU process\n",
    "        \n",
    "    Returns:\n",
    "        The value of the OU process after a time interval `δt`.\n",
    "    '''\n",
    "    N_norm = np.random.randn()\n",
    "    relax = np.exp(-δt / τ)\n",
    "    diffuse = σ * np.sqrt(1 - relax**2)\n",
    "    return x_prev * relax + diffuse * N_norm\n",
    "\n",
    "\n",
    "def OU_process0(num_steps, *, δt=0.1, x0=0, τ=1, σ=2, OU_update=OU_update0):\n",
    "    '''Compute `num_step` values of an OU process.\n",
    "    \n",
    "    Arguments:\n",
    "        num_steps (int): number of time points to be evaluated\n",
    "        δt (float): increment between time points\n",
    "        x0 (float): initital value of the OU process\n",
    "        τ (float): relaxation time, i.e. time after which correlation is lost\n",
    "        σ (float): standard deviation each point of the OU process\n",
    "    \n",
    "    Returns:\n",
    "        Array of size `num_steps + 1` of simulated OU process values.\n",
    "        The first value of the array is `x0`.\n",
    "    '''\n",
    "    x = [x0]\n",
    "    for _ in range(num_steps):\n",
    "        x.append(OU_update(x_prev=x[-1], δt=δt, σ=σ, τ=τ)) \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "OU_update0(x_prev=1, δt=0.1, τ=τ, σ=σ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "# OU process parameters\n",
    "τ = 1\n",
    "σ = 2\n",
    "\n",
    "# Simulation parameters\n",
    "num_steps = 100\n",
    "δt = 0.1\n",
    "x0 = 1\n",
    "\n",
    "# OU Simulation\n",
    "x = OU_process0(num_steps, δt=δt, σ=σ, τ=τ)\n",
    "\n",
    "plt.plot(x);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(x) == num_steps + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "because we simulate `num_steps` points, plus `x0` taken as first element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Playing with the OU process\n",
    "\n",
    "OU process is a \"mean reverting\" process, it means:\n",
    "\n",
    "> we can start form any point it will always go back around 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(num_steps=1000, δt=0.2, τ=30, σ=2)\n",
    "t = np.arange(params['num_steps'] + 1) * params['δt']\n",
    "for x0 in range(-100, 101, 50):\n",
    "    x = OU_process0(x0=x0, **params)\n",
    "    plt.plot(t, x)\n",
    "\n",
    "plt.twinx()\n",
    "plt.plot(t, np.exp(-t / params['τ']), color='k', alpha=0.5, lw=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **NOTE:** play with `τ` to see how the process changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "# Histogram of OU process values\n",
    "x = OU_process0(num_steps=100000, σ=σ)\n",
    "bins = np.arange(-3 * σ, 3 * σ, 0.2)\n",
    "bins_centers = bins[:-1] + 0.5 * (bins[1] - bins[0])\n",
    "plt.hist(x, bins=bins, density=True);\n",
    "\n",
    "# \"Analytical\" Gaussian distribution (red line)\n",
    "plt.plot(bins_centers, norm.pdf(bins_centers, 0, σ), '-r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **NOTE**: Model and experiment agree!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "source": [
    "**TIP 1**: Always test before optimize!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test 1: Normality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "σ = 2\n",
    "num_steps = 10000\n",
    "\n",
    "x = OU_process0(num_steps, σ=σ)\n",
    "\n",
    "# Test 1: Normality\n",
    "bins = np.arange(-3 * σ, 3 * σ, 0.1)\n",
    "bins_centers = bins[1:] - 0.5 * (bins[1] - bins[0])\n",
    "\n",
    "y_experim, _ = np.histogram(x, bins=bins, density=True)\n",
    "y_model = norm.pdf(bins_centers, 0, σ)\n",
    "\n",
    "plt.plot(bins_centers, y_experim)\n",
    "plt.plot(bins_centers, y_model);\n",
    "\n",
    "score = ((y_experim - y_model)**2).sum()\n",
    "assert score < 0.02, f'{score}'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test 2: Correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "num_steps = 10000\n",
    "τ = 1\n",
    "δt = 0.1\n",
    "\n",
    "# OU simulation\n",
    "x = OU_process0(num_steps, δt=δt, τ=τ)\n",
    "\n",
    "# Test correlation\n",
    "xc = np.arange(1, int(10 * τ / δt)) * δt\n",
    "corr = [np.corrcoef(x[i:], x[:-i])[0, 1] for i in range(1, len(xc) + 1)]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xc, corr, label='Correlation coefficient')\n",
    "plt.plot(xc, np.exp(-xc/τ), label='Exponential function')\n",
    "plt.xlabel('Time lags');\n",
    "\n",
    "assert corr[0] > 0.9\n",
    "assert abs(corr[-1]) < 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test 3: initial value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_steps = 3\n",
    "for x0 in range(10):\n",
    "    x = OU_process0(num_steps, x0=x0)\n",
    "    assert x0 == x[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **TIP**: Test using `assert`. These can be turned into unit tests\n",
    "> using `pytest`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarks\n",
    "\n",
    "> **Magic commands:** IPython has special [magic commands](https://ipython.readthedocs.io/en/stable/interactive/magics.html) \n",
    "> that cannot be used in plain python.\n",
    ">\n",
    "> Magics start with:\n",
    "> - `%` for commands working on a one-line of code (can be anywhere in a cell)\n",
    "> - `%%` cell-level magics operate on all the lines in a cell and must be in the first line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_steps = int(1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## %timeit magic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b0 = %timeit -o -n1 -r1 OU_process0(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, we have a way to measure the execution time, but where the time is spent?\n",
    "\n",
    "## %prun magic\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%prun OU_process0(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`%prun` uses [Python's profiler](https://docs.python.org/3.6/library/profile.html)\n",
    "to find which function \"eats\" the most time.\n",
    "\n",
    "**Legend:**\n",
    "\n",
    "- `ncalls`: number of calls.\n",
    "- `tottime`: total time spent in the given function (and excluding time made in calls to sub-functions)\n",
    "- `percall`: is the quotient of `tottime / ncalls`\n",
    "- `cumtime`: is the cumulative time spent in this and all subfunctions (from invocation till exit).     \n",
    "  This figure is accurate even for recursive functions.\n",
    "\n",
    "### Observations\n",
    "\n",
    "- Total execution time is 5.6 secs\n",
    "- Cumulative time in `OU_update0`: 4.9 secs (total), 4.1 (excluding sub-function calls)\n",
    "- The for-loop in `OU_process0` is negligible at this point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Line profiling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext line_profiler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run `OU_process0(num_steps)` while profiling the function `OU_update0` line-by-line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f OU_update0 OU_process0(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Legend:**\n",
    "\n",
    "- `Hits`: number of times line is executed\n",
    "- `Time`: total time spent executing the line\n",
    "- `Per Hit`: time spent executing the line once\n",
    "- `% Time`: percentage of time spent executing the line\n",
    "\n",
    "### Observations\n",
    "\n",
    "- Most time is spent in algebraic computations\n",
    "- Random number generation takes 20% of the time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimizations\n",
    "\n",
    "> **TIP:** Focus on the bottleneck, stop when performances are acceptable.\n",
    "\n",
    "> **TIP:**  Try not to compromise clarity for speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x_ref = OU_process0(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use math instead of numpy\n",
    "\n",
    "We only doing scalar computations so we can replace `numpy` with `math` from the standard library:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exsercise 1**: create a function `OU_update1` with the same signature as \n",
    "> `OU_update0`, trying not to use `numpy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ex1.py\n",
    "import math\n",
    "\n",
    "def OU_update1(*, x_prev, δt, τ, σ):\n",
    "    N_norm = np.random.randn()\n",
    "    relax = math.exp(-δt / τ)\n",
    "    diffuse = σ * math.sqrt(1 - relax**2)\n",
    "    return x_prev * relax + diffuse * N_norm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x = OU_process0(num_steps, OU_update=OU_update1)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b01 = %timeit OU_process0(num_steps, OU_update=OU_update1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f OU_update1 OU_process0(num_steps, OU_update=OU_update1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results \n",
    "\n",
    "- We get a 4x improvement!\n",
    "- Now random number generation takes ~40% of time, becoming a bottleneck!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Move RNG out of the loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exsercise 2**: create a new function `OU_process1`\n",
    "> moving the random-number generation outside the loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ex2.py\n",
    "def OU_update2(*, x_prev, δt, τ, σ, N_norm):\n",
    "    relax = math.exp(-δt / τ)\n",
    "    diffuse = σ * math.sqrt(1 - relax**2)\n",
    "    return x_prev * relax + diffuse * N_norm\n",
    "\n",
    "def OU_process1(num_steps, *, δt=0.1, x0=0, τ=1, σ=2, OU_update=OU_update2):\n",
    "    N_norm = np.random.randn(num_steps)\n",
    "    x = [x0]\n",
    "    for i in range(num_steps):\n",
    "        x.append(OU_update(x_prev=x[-1], δt=δt, σ=σ, τ=τ, N_norm=N_norm[i])) \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x = OU_process1(num_steps, OU_update=OU_update2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -n1 -r1 OU_process1(num_steps, OU_update=OU_update2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%prun OU_process1(num_steps, OU_update=OU_update2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f OU_update2 OU_process1(num_steps, OU_update=OU_update2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "- 40% improvements, RNG is now negligible\n",
    "- next step is improving the algebraic computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Streamline the loop\n",
    "\n",
    "- We compute `relax` and `diffuse` only once\n",
    "- We get rid of the inner function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exsercise 3**: Create a new function `OU_process2` which:\n",
    "> - computes `relax` and `diffuse` only once\n",
    "> - gets rid of the inner function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ex3.py\n",
    "def OU_process2(num_steps, *, δt=0.1, x0=0, τ=1, σ=2):\n",
    "    N_norm = memoryview(np.random.randn(num_steps))\n",
    "    \n",
    "    relax = math.exp(-δt / τ)\n",
    "    diffuse = σ * math.sqrt(1 - relax**2)\n",
    "\n",
    "    x = [x0]\n",
    "    for i in range(num_steps):\n",
    "        x.append(\n",
    "            x[-1] * relax + diffuse * N_norm[i]\n",
    "        ) \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x = OU_process2(num_steps)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "%timeit -n1 -r1 OU_process2(num_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%prun OU_process2(num_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "%lprun -f OU_process2 OU_process2(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "- 2x improvement\n",
    "- now most time in:\n",
    "    1. for-loop + list append\n",
    "    2. scalar computations\n",
    "\n",
    "### Lessons\n",
    "\n",
    "- numpy array element access is slow when using a scalar index\n",
    "- speed can be improved using a `memoryview` (still pure python!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More numpy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exsercise 4**: Create a new function `OU_process3` preallocating the output as numpy's array instead of appending to a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ex4.py\n",
    "def OU_process3(num_steps, *, δt=0.1, x0=0, τ=1, σ=2):\n",
    "    N_norm = memoryview(np.random.randn(num_steps))\n",
    "    \n",
    "    relax = math.exp(-δt / τ)\n",
    "    diffuse = σ * math.sqrt(1 - relax**2)\n",
    "\n",
    "    x = memoryview(np.empty(num_steps + 1))\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "x = OU_process3(num_steps)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "%timeit -n1 -r1 OU_process3(num_steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "%lprun -f OU_process3 OU_process3(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "- ~10% improvement\n",
    "- we cannot further vectorize because the relation is recurrent (ARMA)\n",
    "\n",
    "### Lessons\n",
    "\n",
    "- appending to a list is relatively fast (not a bottleneck here)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "OU_process_nb = numba.jit(OU_process3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below is \"disabled\" to avoid error:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "np.random.seed(1)\n",
    "x = OU_process_nb(num_steps, δt=0.1, x0=0, τ=1, σ=2)\n",
    "assert np.allclose(x, x_ref)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "---------------------------------------------------------------------------\n",
    "TypeError                                 Traceback (most recent call last)\n",
    "TypeError: cannot convert native memoryview(float64, 1d, C) to Python object\n",
    "\n",
    "The above exception was the direct cause of the following exception:\n",
    "\n",
    "SystemError                               Traceback (most recent call last)\n",
    "<ipython-input-107-f2ff2af6f7a9> in <module>()\n",
    "      1 np.random.seed(1)\n",
    "----> 2 x = OU_process_nb(num_steps, δt=0.1, x0=0, τ=1, σ=2)\n",
    "      3 assert np.allclose(x, x_ref)\n",
    "\n",
    "SystemError: LiftedLoop(<function OU_process3 at 0x1a1f049d90>) returned a result with an error set\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ops! Looks like `memoryviews` confuse numba here, let's try without them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exsercise 5**: write a new function `OU_process_nb` similar to `OU_process3` but without using memoryvew. Decorate it with `@numba.jit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load ex5.py\n",
    "@numba.jit\n",
    "def OU_process_nb(num_steps, *, N_norm, δt=0.1, x0=0., τ=1., σ=2.):\n",
    "    relax = math.exp(-δt / τ)\n",
    "    diffuse = σ * math.sqrt(1 - relax**2)\n",
    "\n",
    "    x = np.empty(num_steps + 1)\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we are at it, we also decouple the RNG form the simulation, by passing already generated random numbers.\n",
    "\n",
    "### Test (numba)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_nb(num_steps, N_norm=N_norm, δt=0.1, x0=0, τ=1, σ=2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark (numba)\n",
    "\n",
    "- We fix the seed for more accurate benchmarks\n",
    "- We time the whole simulation (RNG + OU_process)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_nb(num_steps, N_norm=N_norm, δt=0.1, x0=0, τ=1, σ=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "- An additional 7.5x improvement. Good!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementations (cython)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 0** - We start compiling the code with no optimizations (slow):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "from libc.math cimport exp, sqrt\n",
    "\n",
    "def OU_process_cy0(num_steps, *, N_norm, dt=0.1, x0=0, tau=1, sigma=2):\n",
    "    relax = exp(-dt / tau)\n",
    "    diffuse = sigma * sqrt(1 - relax**2)\n",
    "\n",
    "    x = np.empty(num_steps + 1)\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 1** - We add type declaration *only the important variables* (arrays and loop variables):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "from libc.math cimport exp, sqrt\n",
    "\n",
    "def OU_process_cy1(int num_steps, *, double[:] N_norm, dt=0.1, \n",
    "                   x0=0, tau=1, sigma=2):\n",
    "    cdef Py_ssize_t i\n",
    "    cdef double relax, diffuse\n",
    "    cdef double[:] x = np.empty(num_steps + 1)\n",
    "    \n",
    "    relax = exp(-dt / tau)\n",
    "    diffuse = sigma * sqrt(1 - relax**2)\n",
    "\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return np.asfarray(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 2** - Add type declarations for all variables and remove safety checks (dangerous):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -a\n",
    "\n",
    "import cython\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "from libc.math cimport exp, sqrt\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.cdivision(True)\n",
    "@cython.wraparound(False)\n",
    "def OU_process_cy2(int num_steps, *, double[:] N_norm, double dt=0.1, \n",
    "                   double x0=0, double tau=1, double sigma=2):\n",
    "    cdef Py_ssize_t i\n",
    "    cdef double relax, diffuse\n",
    "    cdef double[:] x = np.empty(num_steps + 1)\n",
    "    \n",
    "    relax = exp(-dt / tau)\n",
    "    diffuse = sigma * sqrt(1 - relax**2)\n",
    "\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return np.asfarray(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test (cython)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Test the 3 cython versions against the reference implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy0(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy1(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy2(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Success**: cython versions give the same results as the reference!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmarks (cython)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We fix the seed for more accurate benchmarks\n",
    "- We time the whole simulation (RNG + OU_process)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy0(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy1(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy2(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Success**: performance similar to Numba"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython Type References\n",
    "\n",
    "#### Numpy's types\n",
    "\n",
    "After `cimport numpy as np` use numpy's types appending `_t`. Examples: \n",
    "\n",
    "- `np.int8_t`, `np.int16_t`, `np.int32_t`, `np.int64_t`,\n",
    "- `np.uint8_t`, `np.uint16_t`, `np.uint32_t`, `np.uint64_t`,\n",
    "- `np.float16_t`, `np.float32_t`, `np.float_64`\n",
    "\n",
    "#### double\n",
    "\n",
    "`double` is equivalent to `np.float_64`, but does not require numpy.\n",
    "\n",
    "#### Index type\n",
    "\n",
    "For variables used as index, cython suggests using `Py_ssize_t` instead of other integer types for best performances.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### On a closer look ...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1 \n",
    "N_norm = np.random.randn(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Most of the time is now spent in generating random numbers!!!\n",
    "\n",
    "Let's benchmark only the OU_process part:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "num_steps = int(1e6)\n",
    "N_norm = np.random.randn(num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 0**, no type declarations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = OU_process_cy0(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 1**, with variable declarations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = OU_process_cy1(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cython Version 2**, no safety checks (dangerous):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = OU_process_cy2(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numba Version**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit \n",
    "x = OU_process_nb(num_steps, N_norm=N_norm, δt=0.1, x0=0, τ=1, σ=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lessons\n",
    "\n",
    "- Cython does not accept UTF-8 identifiers\n",
    "- the `-a` flags helps finding \"hot-spots\" in cython code\n",
    "- Numba and cython reach the same C-like speeds (Numba is slightly faster in this case)\n",
    "- Removing bounds checks, negative index wraparound and C-division\n",
    "  does not make a big difference in performance but make the code \n",
    "  more brittle, risking to crash the interpreter if there is a bug. **Not worth it!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ranom number generation\n",
    "\n",
    "- We are already using a single vectorized numpy call\n",
    "\n",
    "- Can we do better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Faster Random Numbers\n",
    "\n",
    "Numpy uses:\n",
    "\n",
    "- MT19937 for generating uniform random numbers in 0..1\n",
    "- Box-Muller method for converting uniform random numbers to normal\n",
    "\n",
    "Both are suboptimal:\n",
    "\n",
    "- There are RNG with higher speed, better statistical properties and jump-ahead features\n",
    "- Zigurat transformation is superior to Box-Muller in terms of speed and accuracy\n",
    "\n",
    "### randomgen: the future of numpy\n",
    "\n",
    "[randomgen](https://github.com/bashtage/randomgen) is a numpy-compatible package \n",
    "aiming to replace the current numpy implementation.\n",
    "It provides many RNG engines that are superior to MT19937. \n",
    "randomgen also implements the Ziguratt method for normals.\n",
    "\n",
    "Even if not in numpy yet, it can be used today with `pip install randomgen`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from randomgen import RandomGenerator, Xoroshiro128"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "np.random.randn(int(1e7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rg = RandomGenerator(Xoroshiro128())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "rg.randn(int(1e7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> About 2-3x speed improvement!\n",
    "\n",
    "> **WARNING**: random numbers are different in the two cases!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Numpy's `randn` + cython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy1(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Randomgen's `randn` + cython:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "N_norm = rg.randn(num_steps)\n",
    "x = OU_process_cy1(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Randomgen's `randn` + numba:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "N_norm = rg.randn(num_steps)\n",
    "x = OU_process_nb(num_steps, N_norm=N_norm, δt=0.1, x0=0, τ=1, σ=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f():\n",
    "    N_norm = rg.randn(num_steps)\n",
    "    x = OU_process_nb(num_steps, N_norm=N_norm, δt=0.1, x0=0, τ=1, σ=2)\n",
    "    \n",
    "blast = %timeit -o f()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Overall >2x speed improvement in the simulation. Not bad!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b0.average / blast.average"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The latest version using `randomgen` and `numba` achieves \n",
    "~350x speed improvement compared to the first naive version!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profiling cython code\n",
    "\n",
    "*See [Python Line_profiler and Cython function](https://stackoverflow.com/questions/24144931/python-line-profiler-and-cython-function)*\n",
    "\n",
    "\n",
    "We have max performances for RNG, can we go back to optimize the cython more\n",
    "with line_profiler?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cython.Compiler import Options\n",
    "Options.linetrace = True\n",
    "Options.binding = True\n",
    "Options.docstrings = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -f --compile-args=-DCYTHON_TRACE=1\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "from libc.math cimport exp, sqrt\n",
    "\n",
    "def OU_process_cy1_lp(int num_steps, *, double[:] N_norm, dt=0.1, x0=0, tau=1, sigma=2):\n",
    "    cdef Py_ssize_t i\n",
    "    cdef double relax, diffuse\n",
    "    cdef double[:] x = np.empty(num_steps + 1)\n",
    "    \n",
    "    relax = exp(-dt / tau)\n",
    "    diffuse = sigma * sqrt(1 - relax**2)\n",
    "\n",
    "    x[0] = x0\n",
    "    for i in range(num_steps):\n",
    "        x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f OU_process_cy1_lp OU_process_cy1_lp(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 😢 Either I did something wrong or is currently not working https://github.com/rkern/line_profiler/issues/94\n",
    "\n",
    "Actually, I found that `line_profiler` works with cython when moved out of the notebook. See below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Packaging Cython\n",
    "\n",
    "*References:*\n",
    "- [Cython for NumPy users](http://docs.cython.org/en/latest/src/userguide/numpy_tutorial.html)\n",
    "\n",
    "After the implementation is optimized we want to move it to a stand-alone module for easier distribution and testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1** Move the cython code to a `.pyx` file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load ou_process.pyx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setp 2** Create a `setup.py` (here I added provisions for profiling):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load setup.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3** Compile the cython module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!python setup.py build_ext --inplace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4**: import, test, benchmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ou_process import OU_process_cy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N_norm = np.random.randn(num_steps)\n",
    "x = OU_process_cy(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)\n",
    "assert np.allclose(x, x_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = OU_process_cy(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test passed! 🎉🎉🎉\n",
    "\n",
    "## Reloading cython code\n",
    "\n",
    "In case of modifications, you need to recompile the code and restart the kernel. \n",
    "\n",
    "`pyximport` has some limited support for reloading, you may want to check it out\n",
    "[Compiling with pyximport](https://cython.readthedocs.io/en/latest/src/userguide/source_files_and_compilation.html#compiling-with-pyximport).\n",
    "\n",
    "> Reloading is automating when using the `%%cython` magic inside the notebook! 🎉"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Line profiling (try2)\n",
    "\n",
    "Let's try again line_profiler on the cython function we moved to a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f OU_process_cy OU_process_cy(num_steps, N_norm=N_norm, dt=0.1, x0=0, tau=1, sigma=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Works this time!\n",
    "\n",
    "### Line_profiler output\n",
    "```\n",
    "\n",
    "Total time: 0.325194 s\n",
    "File: ou_process.pyx\n",
    "Function: OU_process_cy at line 5\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "     5                                           def OU_process_cy(int num_steps, *, double[:] N_norm, dt=0.1,\n",
    "     6                                                           x0=0, tau=1, sigma=2):\n",
    "     7                                               cdef Py_ssize_t i\n",
    "     8                                               cdef double relax, diffuse\n",
    "     9         1        691.0    691.0      0.2      cdef double[:] x = np.empty(num_steps + 1)\n",
    "    10                                           \n",
    "    11         1          3.0      3.0      0.0      relax = exp(-dt / tau)\n",
    "    12         1          1.0      1.0      0.0      diffuse = sigma * sqrt(1 - relax**2)\n",
    "    13                                           \n",
    "    14         1          1.0      1.0      0.0      x[0] = x0\n",
    "    15         1          1.0      1.0      0.0      for i in range(num_steps):\n",
    "    16   1000000     324449.0      0.3     99.8          x[i+1] = x[i] * relax + diffuse * N_norm[i]\n",
    "    17         1         48.0     48.0      0.0      return x\n",
    "```\n",
    "\n",
    "### Results\n",
    "- 99% of time is spent on the line:\n",
    "      x[i+1] = x[i] * relax + diffuse * N_norm[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallelization\n",
    "\n",
    "See the Part II notebook: [Parallelism with Joblib](Parallelism%20with%20Joblib.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "py36_optimize"
  },
  "kernelspec": {
   "display_name": "Python 3.6 (py36_optimize)",
   "language": "python",
   "name": "py36_optimize"
  },
  "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.6"
  },
  "nteract": {
   "version": "0.11.6"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "591.867px",
    "left": "0px",
    "right": "1003px",
    "top": "111.133px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}