scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy – A collection of scientific algorithms and data structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much of this lecture is based on the excellent work of  J.R. Johansson who published his work here: \n",
    "J.R. Johansson (jrjohansson at gmail.com) [http://jrjohansson.github.io](http://jrjohansson.github.io).\n",
    "Some examples have been pulled from https://github.com/scipy/scipy-cookbook/\n",
    "\n",
    "The SciPy framework builds on top of the low-level NumPy framework for multidimensional arrays, and provides a large number of higher-level scientific algorithms. Some of the topics that SciPy covers are:\n",
    "\n",
    "* Special functions ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))\n",
    "* Integration and ODE  ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))\n",
    "* Optimization ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))\n",
    "* Interpolation ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))\n",
    "* Fourier Transforms ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))\n",
    "* Signal Processing ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))\n",
    "* Linear Algebra ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))\n",
    "* Sparse Eigenvalue Problems ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))\n",
    "* Statistics ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))\n",
    "* Multi-dimensional image processing ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))\n",
    "* File IO ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))\n",
    "\n",
    "Each submodule in SciPy  provides a number of functions and classes that can be used to solve problems in their respective domain.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:48.710925Z",
     "start_time": "2019-04-10T13:17:48.638117Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 id=\"tocheading\">Table of Contents</h1>\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.126028Z",
     "start_time": "2019-04-10T13:17:49.642680Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#some code to help display things in this notebook\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.162593Z",
     "start_time": "2019-04-10T13:17:52.141620Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# plot style so figures are more readable. Hopefully\n",
    "plt.rcParams['figure.figsize'] = (11, 7.5)\n",
    "plt.rcParams['font.size'] = 15\n",
    "plt.rcParams['lines.linewidth'] = 3\n",
    "plt.rcParams['xtick.labelsize'] = 14\n",
    "plt.rcParams['ytick.labelsize'] = 14\n",
    "plt.rcParams['axes.labelsize'] = 15\n",
    "plt.rcParams['axes.titlesize'] = 15\n",
    "plt.rcParams['legend.fontsize'] = 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear algebra module contains a lot of matrix related functions, including linear equation solving, eigenvalue solvers, matrix functions (for example matrix-exponentiation), a number of different decompositions (SVD, LU, cholesky), etc. \n",
    "\n",
    "Detailed documetation is available at: http://docs.scipy.org/doc/scipy/reference/linalg.html\n",
    "\n",
    "Here we will look at how to use some of these functions:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear equation systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear equation systems on the matrix form\n",
    "\n",
    "$A x = b$\n",
    "\n",
    "where $A$ is a matrix and $x,b$ are vectors can be solved like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.207788Z",
     "start_time": "2019-04-10T13:17:52.175658Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.301861Z",
     "start_time": "2019-04-10T13:17:52.216850Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A = (np.arange(100) + np.random.normal(0, 10, size=(100))).reshape(10, 10)\n",
    "b = np.arange(10)\n",
    "A.shape, b.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.857323Z",
     "start_time": "2019-04-10T13:17:52.317313Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.matshow(A)\n",
    "plt.grid(False) # turn off the grid for this plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:52.937662Z",
     "start_time": "2019-04-10T13:17:52.877962Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x = linalg.solve(A, b)\n",
    "np.allclose(A @ x - b, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvalues and eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalue problem for a matrix $A$:\n",
    "\n",
    "$\\displaystyle A v_n = \\lambda_n v_n$\n",
    "\n",
    "where $v_n$ is the $n$th eigenvector and $\\lambda_n$ is the $n$th eigenvalue.\n",
    "\n",
    "To calculate eigenvalues of a matrix, use the `eigvals` and for calculating both eigenvalues and eigenvectors, use the function `eig`. Both return complex numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:53.870959Z",
     "start_time": "2019-04-10T13:17:53.846835Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A = np.random.rand(4, 4)\n",
    "\n",
    "eigen_values = linalg.eigvals(A)\n",
    "eigen_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:53.909409Z",
     "start_time": "2019-04-10T13:17:53.878389Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "values, vectors = linalg.eig(A)\n",
    "values, vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also more specialized solvers, like the `eigh` for Hermitian matrices. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise:  Simple Linear Algebra\n",
    "\n",
    "1. Create a 10x10 matrix $A$ with random entries\n",
    "\n",
    "2. Create a symmetric matrix from $A$ using its transpose (use A.T to transpose) and then plot the new matrix using `plt.matshow`\n",
    "\n",
    "3. Use `linalg.eigh` to get the eigenvalues and vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:53.932739Z",
     "start_time": "2019-04-10T13:17:53.920340Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_linalg.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scipy Spatial\n",
    "\n",
    "Methods and data structures to help deal with spatial data. Spatial data means that we can define some sort of metric on the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:55.984770Z",
     "start_time": "2019-04-10T13:17:53.943086Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import spatial\n",
    "N = 5000\n",
    "points = np.random.uniform(0, 1, (N, 3)) # create 5000 points in 3D\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "ax.scatter(points[:, 0], points[:, 1], points[:, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: SciPy Spatial\n",
    "\n",
    "Find all points that have a distance of 0.1 to the center point (0.5, 0.5, 0.5). \n",
    "\n",
    "1. Read the first paragraph of https://en.wikipedia.org/wiki/K-d_tree. \n",
    "2. Create a KDtree using `scipy.spatial.cKDTree(points)`\n",
    "3. Check the scipy documentation of `scipy.spatial.cKDTree` for a method finding all points within a given radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:55.996428Z",
     "start_time": "2019-04-10T13:17:55.992326Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_spatial.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we test the SciPy solution versus the NumPy solution from earlier for finding the nearest neighbour for each point in the point cloud. (See NumPy notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:56.014157Z",
     "start_time": "2019-04-10T13:17:55.999932Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "distances, indices = kdtree.query(points, k=2)\n",
    "distances[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:56.056914Z",
     "start_time": "2019-04-10T13:17:56.040430Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "distances = np.sum((points.reshape(N, 1, 3) - points)**2, axis=2)\n",
    "distances[np.diag_indices_from(distances)] = np.inf\n",
    "np.argmin(distances, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Special functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A large number of mathematical special functions are important for many computional physics problems. SciPy provides implementations of a very extensive set of special functions. For details, see the list of functions in the reference documention at http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special. \n",
    "\n",
    "- Airy functions\n",
    "- Elliptic functions and integrals\n",
    "- Bessel functions\n",
    "- Gamma functions\n",
    "- etc...."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Polynomials\n",
    "\n",
    "The module contains a large number of functions for evaluating polynomials and calculating their roots.\n",
    "\n",
    "The following plot shows just a few of them.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:57.512706Z",
     "start_time": "2019-04-10T13:17:56.067256Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import special\n",
    "\n",
    "functions = [special.eval_chebys, special.eval_laguerre, special.eval_legendre, special.eval_chebyu]\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, sharex=True, sharey=True )\n",
    "\n",
    "x = np.linspace(-2., 2., 200)\n",
    "\n",
    "for ax, f in zip(np.ravel(axs), functions):\n",
    "    ax.set_title(f.__name__)\n",
    "    for n in range(6):\n",
    "        ax.plot(x, f(n, x))\n",
    "    ax.set_ylim([-3, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical Functions\n",
    "\n",
    "The module also contains a host of useful functions for statistical computations. for example the so called 'error function' `scipy.special.erf`\n",
    "\n",
    "\\begin{align}\n",
    "\\operatorname {erf}(x)&={\\frac {1}{\\sqrt {\\pi }}}\\int _{-x}^{x}e^{-t^{2}},dt \\\\ &={\\frac {2}{\\sqrt {\\pi }}}\\int _{0}^{x}e^{-t^{2}}\\,dt.\n",
    "\\end{align}\n",
    "\n",
    "Which allows you to get the gaussian probability dsitribution\n",
    "$$\n",
    "\\phi(z) = \\frac 1 2 \\left(1 + erf\\left(\\frac{z}{sqrt(2)}\\right) \\right)\n",
    "$$\n",
    "\n",
    "The module `scipy.stats` contains nicer methods to handle statistical methods. More on that later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:59.040290Z",
     "start_time": "2019-04-10T13:17:57.527868Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-3, 3)\n",
    "plt.plot(x, special.erf(x), label='$erf$')\n",
    "plt.plot(x, 0.5 * (1 + special.erf(x/np.sqrt(2))), label='$phi$')\n",
    "plt.xlabel('$x$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bessel Functions\n",
    "To demonstrate the typical usage of special functions we will look in more detail at the Bessel functions.\n",
    "The `scipy.special` module includes a large number of Bessel-functions\n",
    "Here we will use the functions $j_n$ and $y_n$, which are the Bessel functions \n",
    "of the first and second kind and real-valued order.\n",
    "\n",
    "Bessel functions are a family of solutions to Bessel’s differential equation with real or complex order alpha:\n",
    "$$\n",
    "x^2 \\frac{d^2 y}{dx^2} +x  \\frac{dy}{dx} +(x^2−α^2)y=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:17:59.994833Z",
     "start_time": "2019-04-10T13:17:59.051434Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import special\n",
    "x = np.linspace(0, 10, 100)\n",
    "\n",
    "for n in range(4):\n",
    "    plt.plot(x, special.jn(n, x), label=f'$J_{n}(x)$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Among other uses, these functions arise in wave propagation problems such as the vibrational modes of a thin drum head. Here is an example of a circular drum head anchored at the edge.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:00.069827Z",
     "start_time": "2019-04-10T13:18:00.009579Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# fix the radius of the drum\n",
    "DRUM_RADIUS = 1\n",
    "\n",
    "def drumhead_height(n, k, distance, angle, t):\n",
    "    '''\n",
    "    Solution to the Drumhead problem with a fixed membrane on the edge of the drum\n",
    "    See https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane\n",
    "    '''\n",
    "    kth_zero = special.jn_zeros(n, k)[-1] # fix the drum on the outter edge to be zero.\n",
    "    return np.cos(t) * np.cos(n*angle) * special.jn(n, distance*kth_zero/DRUM_RADIUS)\n",
    "\n",
    "\n",
    "# create two lists containing all combinations of radius and theta we want to plot\n",
    "r = np.linspace(0, DRUM_RADIUS, 50)\n",
    "t = np.linspace(0, 2*np.pi, 50)\n",
    "radius, theta = np.meshgrid(r, t) # this creates all combinations we are interested in\n",
    "\n",
    "x = radius * np.cos(theta)\n",
    "y = radius * np.sin(theta)\n",
    "z = drumhead_height(1, 1, radius, theta, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:01.770589Z",
     "start_time": "2019-04-10T13:18:00.095294Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "ax.plot_surface(x, y, z, cmap='YlGnBu')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:57.559319Z",
     "start_time": "2019-04-10T13:18:57.538720Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# some code to create a nice animation within the notebook.\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from functools import partial\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "\n",
    "my_drum = partial(drumhead_height, n=2, k=1, distance=radius, angle=theta)\n",
    "\n",
    "frames = 200\n",
    "fps = 25\n",
    "\n",
    "fig = plt.figure(figsize=(19.2, 10.8), dpi=100)  # FullHD\n",
    "ax = Axes3D(fig)\n",
    "ax.plot_surface(x, y, my_drum(t=0), cmap='YlGnBu')\n",
    "\n",
    "bar = tqdm_notebook(total=frames)\n",
    "\n",
    "def init():\n",
    "    pass\n",
    "\n",
    "\n",
    "def update(i):\n",
    "    bar.update()\n",
    "   \n",
    "    t = i / fps\n",
    "    ax.cla()\n",
    "    z = my_drum(t=t)\n",
    "    ax.plot_surface(x, y, z, cmap='YlGnBu', vmin=-0.5, vmax=0.5)\n",
    "    ax.set_zlim(-0.5, 0.5)\n",
    "    \n",
    "    \n",
    "ani = FuncAnimation(fig, update, init_func=init, frames=frames, interval=1000 / fps, blit=False, repeat=False)\n",
    "ani.save('animation_membrane.mp4', writer=\"ffmpeg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:01.837451Z",
     "start_time": "2019-04-10T13:18:01.815673Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import Video\n",
    "Video('./animation_membrane.mp4', width=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Differential Equations.\n",
    "\n",
    "To demonstrate the usefulness of SciPy lets try to predict what will happen to humanity in case of a zombie outbreak.  This example comes courtesy of Christopher Campo.\n",
    "\n",
    "We will show that, in case of a Zombie outbreak, humanity is inevitably doomed. \n",
    "As shown by [Phillip Munz et al](http://mysite.science.uottawa.ca/rsmith43/Zombies.pdf) we can model a simple outbreak scenario like so.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{S}{dt} &= B - \\beta S Z - \\delta S \\\\\n",
    "\\frac{Z}{dt} &= \\beta S Z + \\gamma R - \\alpha S Z \\\\\n",
    "\\frac{R}{dt} &= \\delta S + \\alpha S Z - \\gamma R\n",
    "\\end{align}\n",
    "\n",
    "Where $S$ is the number of susceptible humans, $Z$ is the number of zombies and $R$ is the number of removed persons which have died either by natural causes or zombie attack. The human birthrate is assumed to be constant and modeled by $B$\n",
    "\n",
    "Susceptibles can become zombies through an encounter with a zombie ($\\beta$). Natural causes of human deaths are parameterized by $\\delta$. Humans in the removed class can be resurrected and become zombies ($\\gamma$).  \n",
    "Zombies can enter the removed class by cutting their heads of or removing their brains ($\\alpha$)\n",
    "\n",
    "The question I'm trying to answer is \n",
    "\n",
    "> __Can Humanity survive the zombie Apocalypse?__\n",
    "\n",
    "Since I don't know how to solve differential equations, I'm going to use SciPy instead.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use SciPy to solve this system of differential equations. The function `solve_ivp` does the entire job for us. \n",
    "\n",
    "If we go to [the SciPy documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) about this function we will read the following.\n",
    "\n",
    "\n",
    ">This function numerically integrates a system of ordinary differential equations given an initial value:\n",
    ">\n",
    "> $y^{\\prime} = f(t, y)$\n",
    ">\n",
    ">$y(t_0) = y_0$\n",
    ">\n",
    ">Here $t$ is a one-dimensional independent variable (time), $y(t)$ is an __n-dimensional vector-valued function__ (state), and an __n-dimensional vector-valued function f(t, y) determines the differential equations__. The goal is to find y(t) approximately satisfying the differential equations, given an initial value y(t0)=y0.\n",
    "\n",
    "\n",
    "So we start by defining some starting values for the problem and the model function $f(t, y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:01.893983Z",
     "start_time": "2019-04-10T13:18:01.847416Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "# Type `solve_ivp?` to get inline documentation about the method in your notebook.\n",
    "\n",
    "solve_ivp?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:02.052831Z",
     "start_time": "2019-04-10T13:18:01.906264Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# fixed model parameters\n",
    "B = 0.1 # birth rate\n",
    "delta = 0.0001 # rate of natural casues of death\n",
    "beta = 0.00016 # transmission coefficent\n",
    "gamma = 0.0001 # resurrection rate\n",
    "alpha = 0.0001 # rate of destroyed zombies\n",
    "\n",
    "# initial conditions\n",
    "S0 = 5000.              # initial population\n",
    "Z0 = 0                 # initial zombie population\n",
    "R0 = 0                 # initial dead population\n",
    "y0 = [S0, Z0, R0]     # initial condition vector\n",
    "\n",
    "\n",
    "def f_model(t, y):\n",
    "    Si = y[0]\n",
    "    Zi = y[1]\n",
    "    Ri = y[2]\n",
    "    # the model equations (see Munz et al. 2009)\n",
    "    f0 = B - beta*Si*Zi - delta*Si\n",
    "    f1 = beta*Si*Zi + gamma*Ri - alpha*Si*Zi\n",
    "    f2 = delta*Si + alpha*Si*Zi - gamma*Ri\n",
    "    return [f0, f1, f2]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we solve the equation and output the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:02.153545Z",
     "start_time": "2019-04-10T13:18:02.057735Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "t_start, t_end = 0, 200  # 200 days of zobie acopalypse\n",
    "solution = solve_ivp(f_model, (t_start, t_end), y0, t_eval=np.linspace(t_start, t_end, 200))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:02.936104Z",
     "start_time": "2019-04-10T13:18:02.165130Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "s = solution.y[0, :]\n",
    "z = solution.y[1, :]\n",
    "plt.plot(solution.t, s, label='survivors')\n",
    "plt.plot(solution.t, z, label='zombies')\n",
    "plt.xlabel('t / days')\n",
    "plt.ylabel('Population')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Save the world.\n",
    "\n",
    "The ODE defined above is a very pessimistic model. There is one unstable equilibrium state if $\\beta = 0$.\n",
    "Your task is to save the world by adding a the notion of a _cure_ to the model above. The cure \n",
    "can move individuals from the $Z$ to the $S$ group.\n",
    "\n",
    "1. Introduce a new constant $\\rho$ which models the effectivness of the cure. \n",
    "2. Add a term to the ODE to model the cure. \n",
    "3. Solve the ODE using `solve_ivp`\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:02.970236Z",
     "start_time": "2019-04-10T13:18:02.955813Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_zombies.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: A more complicated  ODE\n",
    "\n",
    "We use SciPy `solve_ivp` method to numerically solve the following system of coupled differential euqations \n",
    "\n",
    "\\begin{aligned}{{\\dot {\\theta }}_{1}}&={\\frac {6}{ml^{2}}}{\\frac {2p_{\\theta _{1}}-3\\cos(\\theta _{1}-\\theta _{2})p_{\\theta _{2}}}{16-9\\cos ^{2}(\\theta _{1}-\\theta _{2})}}\\\\{{\\dot {\\theta }}_{2}}&={\\frac {6}{ml^{2}}}{\\frac {8p_{\\theta _{2}}-3\\cos(\\theta _{1}-\\theta _{2})p_{\\theta _{1}}}{16-9\\cos ^{2}(\\theta _{1}-\\theta _{2})}}.\n",
    "\\end{aligned}\n",
    "\n",
    "\\begin{aligned}{{\\dot {p}}_{\\theta _{1}}}&={\\frac {\\partial L}{\\partial \\theta _{1}}}=-{\\tfrac {1}{2}}ml^{2}\\left({{\\dot {\\theta }}_{1}}{{\\dot {\\theta }}_{2}}\\sin(\\theta _{1}-\\theta _{2})+3{\\frac {g}{l}}\\sin \\theta _{1}\\right)\\\\{{\\dot {p}}_{\\theta _{2}}}&={\\frac {\\partial L}{\\partial \\theta _{2}}}=-{\\tfrac {1}{2}}ml^{2}\\left(-{{\\dot {\\theta }}_{1}}{{\\dot {\\theta }}_{2}}\\sin(\\theta _{1}-\\theta _{2})+{\\frac {g}{l}}\\sin \\theta _{2}\\right).\n",
    "\\end{aligned} \n",
    "\n",
    "Now we plot the results for the $\\theta_1$ and $\\theta_2$ using \n",
    "\n",
    "\\begin{align}\n",
    "x_1 &= L \\; sin(\\theta_1(t))\\\\\n",
    "y_1 &= - L\\;  cos(\\theta_1(t)) \\\\\n",
    "x_2 &= x_1 + L \\; sin(\\theta_2(t)) \\\\\n",
    "y_2 &= y_1 - L \\; cos(\\theta_2(t))\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.004817Z",
     "start_time": "2019-04-10T13:18:02.994469Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy.constants import g\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "L = 0.5\n",
    "m = 0.1\n",
    "\n",
    "t_start = 0\n",
    "t_end = 20\n",
    "fps = 50\n",
    "frames = (t_end - t_start) * fps\n",
    "\n",
    "\n",
    "def pendulum_model(t, y):\n",
    "    \"\"\"\n",
    "    The right-hand side of the pendulum ODE\n",
    "    \"\"\"\n",
    "    phi_1, phi_2, p_1, p_4 = y[0], y[1], y[2], y[3]\n",
    "\n",
    "    c = (m * L**2)\n",
    "    cos_delta_phi = np.cos(phi_1 - phi_2)\n",
    "    sin_delta_phi = np.sin(phi_1 - phi_2)\n",
    "    denominator = 16 - 9 * cos_delta_phi**2\n",
    "\n",
    "    dphi_1 = 6 / c * (2 * p_1 - 3 * cos_delta_phi * p_4) / denominator\n",
    "    dphi_2 = 6 / c * (8 * p_4 - 3 * cos_delta_phi * p_1) / denominator\n",
    "\n",
    "    dp_1 = -c / 2 * ( dphi_1 * dphi_2 * sin_delta_phi + 3 * g / L * np.sin(phi_1))\n",
    "    dp_4 = -c / 2 * (-dphi_1 * dphi_2 * sin_delta_phi + g / L * np.sin(phi_2))\n",
    "\n",
    "    return [dphi_1, dphi_2, dp_1, dp_4]\n",
    "\n",
    "\n",
    "# initial values\n",
    "y0 = [np.pi / 3, -np.pi / 4, 0, 0.065]\n",
    "\n",
    "\n",
    "solution = solve_ivp(\n",
    "    pendulum_model,\n",
    "    (t_start, t_end),\n",
    "    y0,\n",
    "    t_eval=np.linspace(t_start, t_end, frames),\n",
    ")\n",
    "\n",
    "x1s = + L * np.sin(solution.y[0])\n",
    "y1s = - L * np.cos(solution.y[0])\n",
    "\n",
    "x2s = x1s + L * np.sin(solution.y[1])\n",
    "y2s = y1s - L * np.cos(solution.y[1])\n",
    "\n",
    "plt.plot(x1s, y1s, '-', alpha=0.5, color='xkcd:sky blue')\n",
    "plt.plot(x2s, y2s, '-', alpha=0.5, color='xkcd:dark magenta')\n",
    "plt.axis('off')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.037511Z",
     "start_time": "2019-04-10T13:18:04.012881Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(19.2, 10.8), dpi=100)   # Full HD\n",
    "\n",
    "ax.set_ylim([-1.1, 0.2])\n",
    "ax.set_xlim([1, -1])\n",
    "ax.set_aspect(1)\n",
    "\n",
    "pendulum1, = ax.plot([], [])\n",
    "pendulum2, = ax.plot([], [])\n",
    "\n",
    "\n",
    "bar = tqdm_notebook(total=frames)\n",
    "\n",
    "def init():\n",
    "    pendulum1.set_data([], [])\n",
    "    pendulum2.set_data([], [])\n",
    "    return (pendulum1, pendulum2)\n",
    "\n",
    "\n",
    "def update(n):\n",
    "    bar.update(1)\n",
    "    # update the line data\n",
    "    pendulum1.set_data([0, x1s[n]], [0, y1s[n]])\n",
    "    pendulum2.set_data([x1s[n], x2s[n]], [y1s[n], y2s[n]])\n",
    "    return (pendulum1, pendulum2)\n",
    "\n",
    "\n",
    "anim = animation.FuncAnimation(fig, update, init_func=init, frames=frames, blit=True, interval=1000 / fps)\n",
    "anim.save('animation_pendulum.mp4', writer=\"ffmpeg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.068682Z",
     "start_time": "2019-04-10T13:18:04.043483Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Video(\"./animation_pendulum.mp4\", width=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration: quadrature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical evaluation of a function of the type\n",
    "\n",
    "$\\displaystyle \\int_a^b f(x) dx$\n",
    "\n",
    "is called *numerical quadrature*, or simply *quadature*. SciPy provides a series of functions for different kind of quadrature, for example the `quad`, `dblquad` and `tplquad` for single, double and triple integrals, respectively.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.147016Z",
     "start_time": "2019-04-10T13:18:04.135466Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad, dblquad, tplquad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `quad` function takes a large number of optional arguments, which can be used to fine-tune the behaviour of the function.\n",
    "\n",
    "The basic usage is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.166783Z",
     "start_time": "2019-04-10T13:18:04.156240Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# define a simple function for the integrand\n",
    "def f(x):\n",
    "    return x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:04.986269Z",
     "start_time": "2019-04-10T13:18:04.176201Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x_lower = 0 # the lower limit of x\n",
    "x_upper = 1 # the upper limit of x\n",
    "\n",
    "val, abserr = quad(f, x_lower, x_upper)\n",
    "\n",
    "print(f'integral value={val} , absolute error={abserr}' )\n",
    "x = np.linspace(-1.2, 1.2, 200)\n",
    "plt.plot(x, f(x))\n",
    "plt.fill_between(x, f(x), alpha=0.3, where=(x > 0) & (x < 1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we need to pass extra arguments to integrand function we can use the `args` keyword argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:05.063781Z",
     "start_time": "2019-04-10T13:18:05.003320Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def integrand(x, n):\n",
    "    return special.jn(n, x)\n",
    "\n",
    "x_lower = 0  # the lower limit of x\n",
    "x_upper = 10 # the upper limit of x\n",
    "\n",
    "val, abserr = quad(integrand, x_lower, x_upper, args=(3,))\n",
    "print(val, abserr) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simple functions we can use a lambda function (name-less function) instead of explicitly defining a function for the integrand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:05.101458Z",
     "start_time": "2019-04-10T13:18:05.078394Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "val, abserr = quad(lambda x: np.exp(-x**2), -np.Inf, np.Inf)\n",
    "print(f'numerical= {val}, {abserr}')\n",
    "print(f'analytical solution = {np.sqrt(np.pi)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As show in the example above, we can also use 'Inf' or '-Inf' as integral limits.\n",
    "\n",
    "Higher-dimensional integration works in the same way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:07.444111Z",
     "start_time": "2019-04-10T13:18:05.114836Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def integrand(x, y):\n",
    "    return np.exp(-x**2 - y**2)\n",
    "\n",
    "\n",
    "x, y = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(-2, 2, 100))\n",
    "z = integrand(x, y)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "ax.plot_surface(x, y, z, cmap='YlGnBu')\n",
    "\n",
    "# the integral boundarys can be functions when integrating in more than one dimensions\n",
    "val, abserr = dblquad(integrand, -2, 2, lambda x: -2, lambda x: x) \n",
    "\n",
    "print(f'numerical solution= {val}, {abserr}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how we had to pass lambda functions for the limits for the y integration, since these in general can be functions of x."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise II:  Integrals and Special Functions\n",
    "\n",
    "Find the circumference of an ellipse. The formula for calculating circumference $C$ is as follows.\n",
    "\n",
    "$$\n",
    "\\frac{C}{4a} =  E(\\epsilon)=\\int _{0}^{\\pi /2}{\\sqrt {1-\\epsilon^{2}\\sin ^{2}\\theta }}\\ d\\theta\n",
    "$$\n",
    "\n",
    "where $\\epsilon$ is the eccentricity of the ellipse given by the length and width of the ellipse.\n",
    "\n",
    "$$\n",
    "\\epsilon = \\sqrt{1 - \\frac{b^2}{a^2}}\n",
    "$$\n",
    "\n",
    "compare the solution found using `quad` with the appropriate function from the `scipy.special` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:07.468798Z",
     "start_time": "2019-04-10T13:18:07.447622Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_integrals.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Signal Processing\n",
    "\n",
    "SciPy offers a large amount of utilities to handle time series data. However the methods in the `scipy.signal` are useful for much more than just time series data.\n",
    "\n",
    "## Convolution\n",
    "\n",
    "The convolution of two continuos functions is defined by a simple integral transform.\n",
    "\n",
    "\\begin{aligned}\n",
    "(f*g)(t)&\\,{\\stackrel {\\mathrm {def} }{=}}\\ \\int _{-\\infty }^{\\infty }f(\\tau )g(t-\\tau )\\,d\\tau \\\\&=\\int _{-\\infty }^{\\infty }f(t-\\tau )g(\\tau )\\,d\\tau .\n",
    "\\end{aligned}\n",
    "\n",
    "Convolution happens whenever data gets smeared out by a detector or filter of any kind.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:07.583301Z",
     "start_time": "2019-04-10T13:18:07.479099Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "\n",
    "# create two rectangular signals\n",
    "f = np.repeat([0., 1., 0.], 120)\n",
    "g = np.repeat([0., 1., 0.], 20)\n",
    "\n",
    "result = signal.convolve(f, g, mode='same')  # same output shape as input, see docs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:08.634199Z",
     "start_time": "2019-04-10T13:18:07.588221Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True)\n",
    "\n",
    "ax_orig.plot(f, color='gray')\n",
    "ax_orig.set_title('Original pulse')\n",
    "ax_orig.margins(0, 0.1)\n",
    "\n",
    "ax_win.plot(g, color='gray')\n",
    "ax_win.arrow(70, 0.5, 50.0, 0.0, width=0.08, head_length=8)\n",
    "ax_win.set_title('Filter impulse response')\n",
    "\n",
    "ax_filt.plot(result)\n",
    "ax_filt.set_title('Filtered signal')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The naive algorithm for producing discrete convolutions is not very efficient. It takes $\\mathcal{O}(n^2)$ steps to get a result. SciPy uses a heuristic to automagically find the fastest convolution method. It is based on the *convolution theorem*\n",
    "\n",
    "> The convolution theorem states that\n",
    ">\n",
    " $$\n",
    " \\mathcal{F}\\{f*g\\}=k\\cdot {\\mathcal {F}}\\{f\\}\\cdot {\\mathcal {F}}\\{g\\}\n",
    " $$\n",
    ">\n",
    "> where \n",
    "> $\\mathcal {F}\\{f\\}$ denotes the Fourier transform of $f$, and $\\mathcal {F}\\{g\\}$ the Fourier transform of $g$.\n",
    ">\n",
    "> https://en.wikipedia.org/wiki/Convolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal Windows\n",
    "\n",
    "SciPy conveniently provides a handful of common signal windows ready to be used for convolution or other filtering operations. All of them can be found in `scipy.signal.window` but their are also available in the `scipy.signal` name space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:09.259874Z",
     "start_time": "2019-04-10T13:18:08.646101Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for g in [signal.nuttall(150), signal.hann(150), signal.cosine(150), signal.triang(150), signal.blackman(150)]:\n",
    "    plt.plot(g)\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation\n",
    "\n",
    "The correlation of two functions $f$ and $g$ looks similar to a convolution \n",
    "\n",
    "> For continuous functions $f$ and $g$, the cross-correlation is defined as:\n",
    "> \n",
    "$$\n",
    "(f\\star g)(\\tau )\\ {\\stackrel {\\mathrm {def} }{=}}\\int _{-\\infty }^{\\infty }f^{*}(t)\\ g(t+\\tau )\\,dt,\n",
    "$$\n",
    "> \n",
    "> where $f^{*}$ denotes the complex conjugate of  $f$, and $\\tau$ is the displacement, also known as lag, although a positive value of $\\tau$ actually means that $g(t+\\tau)$ leads $f(t)$.\n",
    "\n",
    "> *https://en.wikipedia.org/wiki/Cross-correlation*\n",
    "\n",
    "The result of the correlation is a function of the delay $\\tau$. This is useful for finding the offset in shifted signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:09.312327Z",
     "start_time": "2019-04-10T13:18:09.278969Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "offset = 0.25\n",
    "t = np.linspace(-1, 1, 600)\n",
    "f = signal.gausspulse(t, fc=3)\n",
    "g = signal.gausspulse( t + offset, fc=3)\n",
    "\n",
    "result = signal.correlate(f, g, mode='same')\n",
    "offset = t[np.argmax(result)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:10.298999Z",
     "start_time": "2019-04-10T13:18:09.326862Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "grid = plt.GridSpec(2, 2, wspace=0.1, hspace=0.4, width_ratios=[1, 1])\n",
    "ax1 = plt.subplot(grid[0, 0])\n",
    "ax2 = plt.subplot(grid[1, 0])\n",
    "ax3 = plt.subplot(grid[:, 1], )\n",
    "\n",
    "ax1.plot(t, f, color='gray')\n",
    "ax1.set_title('$f$')\n",
    "ax1.set_yticks([])\n",
    "ax1.set_xlabel('t')\n",
    "\n",
    "ax2.plot(t, g, color='gray')\n",
    "ax2.set_title('$g$')\n",
    "ax2.set_yticks([])\n",
    "ax2.set_xlabel('t')\n",
    "\n",
    "ax3.plot(t, result)\n",
    "ax3.axvline(offset, color='xkcd:tangerine')\n",
    "ax3.set_title('$(f\\star g)(\\\\tau )$')\n",
    "ax3.set_yticks([])\n",
    "ax3.set_xlabel('τ')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise III:  Signal Filtering\n",
    "\n",
    "If we ever find an alien civilization we will have to communicate with them using binary encoded signals recorded by large radio telescopes.\n",
    "\n",
    "\n",
    "Assuming the aliens use binary signals with 0s and 1s. \n",
    "Given the clock speed and the sample length with which the signal was created. Can you decode the following signal?\n",
    "\n",
    "The plot below shows the shape of an example signal without any noise containing the message [0, 1 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:11.038589Z",
     "start_time": "2019-04-10T13:18:10.314893Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sample_length = 64\n",
    "\n",
    "plt.plot(np.repeat([0, 1, 0], sample_length))\n",
    "plt.text(25, 0.5, '0')\n",
    "plt.text(100, 0.5, '1')\n",
    "plt.text(175, 0.5, '0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:12.073432Z",
     "start_time": "2019-04-10T13:18:11.058235Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%run solutions/create_alien_signal.py\n",
    "\n",
    "alien_signal = np.loadtxt('alien_signal.txt')\n",
    "sample_length = 128\n",
    "\n",
    "# let's sample the signal in the middle of each period\n",
    "clock = np.arange(sample_length // 2, len(alien_signal), sample_length)\n",
    "plt.plot(alien_signal)\n",
    "for c in clock:\n",
    "    plt.axvline(c, color='gray', lw=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:12.099455Z",
     "start_time": "2019-04-10T13:18:12.085376Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_alien_signal.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fourier transforms are one of the universal tools in computational physics, which appear over and over again in different contexts. SciPy provides functions for accessing the classic [FFTPACK](http://www.netlib.org/fftpack/) library, which is an efficient and well tested FFT library written in FORTRAN. The SciPy API has a few additional convenience functions, but overall the API is closely related to the original FORTRAN library.\n",
    "\n",
    "To use the `fftpack` module in a python program, include it using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:12.142444Z",
     "start_time": "2019-04-10T13:18:12.107646Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import fftpack "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate how to do a fast Fourier transform with SciPy, we will create some artificial data first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:13.122014Z",
     "start_time": "2019-04-10T13:18:12.146181Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Number of sample points\n",
    "N = 1000\n",
    "L = 10\n",
    "# sample spacing\n",
    "dt = L / N\n",
    "x = np.linspace(0.0, L, N)\n",
    "\n",
    "freqs = [8, 3, 32]\n",
    "amps = [1, 0.5, 0.25]\n",
    "\n",
    "y = np.zeros(N)\n",
    "for a, f in zip(amps, freqs):\n",
    "    y += a * np.sin(2 * np.pi * f * x)\n",
    "\n",
    "y += np.random.normal(0, 0.75, len(y))\n",
    "plt.plot(x, y, lw=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:13.919352Z",
     "start_time": "2019-04-10T13:18:13.134458Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "yf = fftpack.rfft(y)\n",
    "xf = fftpack.rfftfreq(len(y), dt)\n",
    "\n",
    "plt.plot(xf, yf**2)\n",
    "\n",
    "for f in freqs:\n",
    "    plt.axvline(f, color='gray', alpha=0.3)\n",
    "\n",
    "plt.yscale('log')\n",
    "# plt.xscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we now see a peaks in the spectrum that match the true frequencies in the signal.\n",
    "Another, maybe more convenient,  way to create frequencies from a transformed signal uses the `fftfreq` function.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scipy.stats` module contains a large number of statistical distributions, statistical functions and tests. For a complete documentation of its features, see http://docs.scipy.org/doc/scipy/reference/stats.html.\n",
    "\n",
    "Below we instantiate an object modeling a random variable following the normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:13.973763Z",
     "start_time": "2019-04-10T13:18:13.950784Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "mean = 0\n",
    "sigma = 1\n",
    "X = stats.norm(loc=mean, scale=sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-02T14:46:47.298461Z",
     "start_time": "2018-06-02T14:46:47.287124Z"
    }
   },
   "source": [
    "All distributions in this submodule extend either `scipy.stats.rv_continuous` or `scipy.stats.rv_discrete`.\n",
    "Below is a list of functions supported by most continuous distribution objects\n",
    "\n",
    "```\n",
    "rvs(*args, **kwds) \tRandom variates of given type.\n",
    "pdf(x, *args, **kwds) \tProbability density function at x of the given RV.\n",
    "logpdf(x, *args, **kwds) \tLog of the probability density function at x of the given RV.\n",
    "cdf(x, *args, **kwds) \tCumulative distribution function of the given RV.\n",
    "logcdf(x, *args, **kwds) \tLog of the cumulative distribution function at x of the given RV.\n",
    "sf(x, *args, **kwds) \tSurvival function (1 - cdf) at x of the given RV.\n",
    "logsf(x, *args, **kwds) \tLog of the survival function of the given RV.\n",
    "ppf(q, *args, **kwds) \tPercent point function (inverse of cdf) at q of the given RV.\n",
    "isf(q, *args, **kwds) \tInverse survival function (inverse of sf) at q of the given RV.\n",
    "moment(n, *args, **kwds) \tn-th order non-central moment of distribution.\n",
    "stats(*args, **kwds) \tSome statistics of the given RV.\n",
    "entropy(*args, **kwds) \tDifferential entropy of the RV.\n",
    "expect([func, args, loc, scale, lb, ub, …]) \tCalculate expected value of a function with respect to the distribution.\n",
    "median(*args, **kwds) \tMedian of the distribution.\n",
    "mean(*args, **kwds) \tMean of the distribution.\n",
    "std(*args, **kwds) \tStandard deviation of the distribution.\n",
    "var(*args, **kwds) \tVariance of the distribution.\n",
    "interval(alpha, *args, **kwds) \tConfidence interval with equal areas around the median.\n",
    "__call__(*args, **kwds) \tFreeze the distribution for the given arguments.\n",
    "fit(data, *args, **kwds) \tReturn MLEs for shape (if applicable), location, and scale parameters from data.\n",
    "fit_loc_scale(data, *args) \tEstimate loc and scale parameters from data using 1st and 2nd moments.\n",
    "nnlf(theta, x) \tReturn negative loglikelihood function.\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:15.699907Z",
     "start_time": "2019-04-10T13:18:13.984454Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "N = 1200\n",
    "measurments = X.rvs(size=N)\n",
    "x = np.linspace(-3.5, 3.5, 200)\n",
    "\n",
    "f, [ax1, ax2] = plt.subplots(1, 2, figsize=(16, 8))\n",
    "ax1.hist(measurments, bins=40, density=True)\n",
    "ax1.plot(x, X.pdf(x), color='red')\n",
    "\n",
    "\n",
    "ax2.hist(measurments, bins=40, cumulative=True, density=True, )\n",
    "\n",
    "ax2.plot(x, X.cdf(x), color='red')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discrete functions work in a very similar manner.\n",
    "\n",
    "\n",
    "```\n",
    "rvs(*args, **kwargs) \tRandom variates of given type.\n",
    "pmf(k, *args, **kwds) \tProbability mass function at k of the given RV.\n",
    "logpmf(k, *args, **kwds) \tLog of the probability mass function at k of the given RV.\n",
    "cdf(k, *args, **kwds) \tCumulative distribution function of the given RV.\n",
    "logcdf(k, *args, **kwds) \tLog of the cumulative distribution function at k of the given RV.\n",
    "sf(k, *args, **kwds) \tSurvival function (1 - cdf) at k of the given RV.\n",
    "logsf(k, *args, **kwds) \tLog of the survival function of the given RV.\n",
    "ppf(q, *args, **kwds) \tPercent point function (inverse of cdf) at q of the given RV.\n",
    "isf(q, *args, **kwds) \tInverse survival function (inverse of sf) at q of the given RV.\n",
    "moment(n, *args, **kwds) \tn-th order non-central moment of distribution.\n",
    "stats(*args, **kwds) \tSome statistics of the given RV.\n",
    "entropy(*args, **kwds) \tDifferential entropy of the RV.\n",
    "expect([func, args, loc, lb, ub, …]) \tCalculate expected value of a function with respect to the distribution for discrete distribution.\n",
    "median(*args, **kwds) \tMedian of the distribution.\n",
    "mean(*args, **kwds) \tMean of the distribution.\n",
    "std(*args, **kwds) \tStandard deviation of the distribution.\n",
    "var(*args, **kwds) \tVariance of the distribution.\n",
    "interval(alpha, *args, **kwds) \tConfidence interval with equal areas around the median.\n",
    "__call__(*args, **kwds) \tFreeze the distribution for the given arguments.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:15.754666Z",
     "start_time": "2019-04-10T13:18:15.729704Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create a (discreet) random variable with poissionian distribution\n",
    "Y = stats.poisson(3.5) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:16.738790Z",
     "start_time": "2019-04-10T13:18:15.777953Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "n = np.arange(12)\n",
    "\n",
    "fig, axes = plt.subplots(3,1, sharex=True)\n",
    "\n",
    "# plot the probability mass function (PMF)\n",
    "axes[0].plot(n, Y.pmf(n), 'o')\n",
    "\n",
    "# plot the commulative distribution function (CDF)\n",
    "axes[1].step(n, Y.cdf(n))\n",
    "\n",
    "# plot histogram of 200 random realizations of the stochastic variable X\n",
    "axes[2].hist(Y.rvs(size=200), bins=np.arange(13) - 0.5, density=True, lw=2, edgecolor='w');\n",
    "\n",
    "axes[2].set_xticks(n)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objects support easy access to some of its statistical properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:16.774756Z",
     "start_time": "2019-04-10T13:18:16.751648Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(X.mean(), X.std(), X.var()) # normal distribution\n",
    "print(Y.mean(), Y.std(), Y.var()) # poission distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The module contains some usefull methods to perform statistical tests. \n",
    "\n",
    "Below we tests whether two sets of (independent) random data comes from the same distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:17.630617Z",
     "start_time": "2019-04-10T13:18:16.801919Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create two sets of data \n",
    "N = 500\n",
    "a = X.rvs(size=N)\n",
    "b = stats.norm(loc=5).rvs(size=N)\n",
    "\n",
    "plt.hist(a, bins=50, histtype='step', lw=4, density=True)\n",
    "plt.hist(b, bins=50, histtype='step', lw=4, density=True)\n",
    "\n",
    "\n",
    "t_statistic, p_value = stats.ttest_ind(a, b)\n",
    "\n",
    "print(f't-statistic = {t_statistic:.2f}, p-value = {p_value:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the p value is very large we cannot reject the hypothesis that the two sets of random data have *different* means."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimization (finding minima or maxima of a function) is a large field in mathematics, and optimization of complicated functions or in many variables can be rather complicated. Here we will only look at a few very simple cases. For a more detailed introduction to optimization with SciPy see: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html\n",
    "\n",
    "To use the optimization module in scipy first include the `optimize` module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:17.719596Z",
     "start_time": "2019-04-10T13:18:17.669852Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding Minima"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first look at how to find the minima of a simple function of a single variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:18.385739Z",
     "start_time": "2019-04-10T13:18:17.752548Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 4.5*x**3 + (x-4)**2 + 0.75*x**4 - 20\n",
    "\n",
    "x = np.linspace(-6, 3, 100)\n",
    "plt.plot(x, f(x));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `minimize_scalar` function to find the minima of a scalar function. It provides wrapper around many optimization algorithms you can choose from using the `method` keyword. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:19.206180Z",
     "start_time": "2019-04-10T13:18:18.400048Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "result = optimize.minimize_scalar(f)\n",
    "\n",
    "plt.plot(x, f(x));\n",
    "plt.axvline(result.x, color='xkcd:orange')\n",
    "\n",
    "\n",
    "# add some boundaries\n",
    "result = optimize.minimize_scalar(f, bracket=[-6, -2])\n",
    "plt.axvline(result.x, color='xkcd:dark red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimization in more dimensions.\n",
    "\n",
    "Below you'll find the infamous *Rosenbrock* function in two dimensions. It's very popular for benchmarking global optimization problems. Its solutions are well known and easily parametrized. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:22.808109Z",
     "start_time": "2019-04-10T13:18:20.161476Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "def rosenbrock_2d(x, y, a=1, b=1):\n",
    "    return (a - x)**2 + b*(y - x**2)**2\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "\n",
    "x, y = np.meshgrid(np.linspace(-2 ,2, 100), np.linspace(-3. ,5, 100))\n",
    "ax.plot_surface(x, y, rosenbrock_2d(x, y), cmap='magma_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:25.334212Z",
     "start_time": "2019-04-10T13:18:22.853181Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "result = optimize.minimize(lambda x: rosenbrock_2d(x[0], x[1]), x0=[5, 2])\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "ax.view_init(20, 30)\n",
    "\n",
    "\n",
    "x, y = np.meshgrid(np.linspace(-2 ,2, 100), np.linspace(-3. ,5, 100))\n",
    "ax.plot_surface(x, y, rosenbrock_2d(x, y), cmap='magma_r', alpha=0.6)\n",
    "\n",
    "x, y, z = result.x[0], result.x[1], 20\n",
    "ax.quiver(x, y, z, 0, 0, -1, length=20, arrow_length_ratio=0.04)\n",
    "\n",
    "\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Simple Likelihood minimization\n",
    "\n",
    "Imagine you have a sensor with 400 pixels. The pixel measure the number of photons in in a fixed length of time. Fit the mean of the photon counts using the negative log likelihood of the poisson distribution\n",
    "\n",
    "$$\n",
    "X(k, \\lambda) \\; \\sim \\; \\frac{\\lambda^k}{k!} e^{-\\lambda}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:20.075801Z",
     "start_time": "2019-04-10T13:18:19.220082Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "np.random.seed(42)\n",
    "\n",
    "true_lambda = 2\n",
    "data = poisson.rvs(true_lambda, size=(20, 20))\n",
    "\n",
    "plt.imshow(data, cmap='inferno')\n",
    "plt.colorbar()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:20.130235Z",
     "start_time": "2019-04-10T13:18:20.102090Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_poisson.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More complex Likelihood Minimization\n",
    "\n",
    "This is a very simplified model of how fits work in particle physics, we try to model a measurement of particle mass\n",
    "with a gaussian distributed signal and an exponentially distributed background using an unbinned likelihood fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "e_min = 75\n",
    "e_max = 175\n",
    "\n",
    "higgs_signal = np.random.normal(126, 5, 500)\n",
    "background = np.random.exponential(50, size=20000)\n",
    "background = background[(background > e_min) & (background < e_max)]\n",
    "\n",
    "measured = np.append(higgs_signal, background)\n",
    "\n",
    "\n",
    "plt.hist(measured, bins=100)\n",
    "plt.xlabel('$m / \\mathrm{GeV}$')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our combined pdf as a weighted sum of a gaussian and an exponential distribution, where\n",
    "$p$ is the partition of signal events and (1 - p) is the partition of background events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, expon\n",
    "\n",
    "\n",
    "def pdf(x, mean, std, tau, p, e_min, e_max):\n",
    "    N = np.exp(-e_min / tau) - np.exp(-e_max / tau)\n",
    "    return (\n",
    "        p * norm.pdf(x, mean, std) \n",
    "        + (1 - p) / N * expon.pdf(x, scale=tau)\n",
    "    )\n",
    "\n",
    "def neg_log_likelihood(params, data, e_min, e_max):\n",
    "    return -np.sum(np.log(pdf(data, *params, e_min=e_min, e_max=e_max)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minimizing using initial guess and bounds.\n",
    "# see docs of minimize\n",
    "\n",
    "result = optimize.minimize(\n",
    "    neg_log_likelihood,\n",
    "    x0=(130, 2, 30, 0.2),\n",
    "    bounds=[\n",
    "        (None, None), # No bounds for mean\n",
    "        (1e-30, None), # std > 0\n",
    "        (1e-30, None), # tau > 0\n",
    "        (0, 1), # 0 <= p <= 1\n",
    "    ],\n",
    "    args=(measured, e_min, e_max)\n",
    ")\n",
    "\n",
    "print(result.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Higgs mass is {:.2f} ± {:.2} GeV'.format(result.x[0], np.sqrt(result.hess_inv.todense()[0, 0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov = result.hess_inv.todense()\n",
    "A = np.diag(1 / np.sqrt(np.diag(cov)))\n",
    "\n",
    "cor = A @ cov @ A.T\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "mat = ax.matshow(cor, cmap='RdBu_r', vmin=-1, vmax=1)\n",
    "\n",
    "ax.grid(False)\n",
    "\n",
    "names = ['mean', 'std', 'tau', 'p']\n",
    "\n",
    "ax.set_xticks(np.arange(len(names)))\n",
    "ax.set_xticklabels(names)\n",
    "ax.set_yticks(np.arange(len(names)))\n",
    "ax.set_yticklabels(names)\n",
    "ax.margins(0)\n",
    "\n",
    "fig.colorbar(mat, ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist, edges, plot = plt.hist(measured, bins=100)\n",
    "m = np.linspace(e_min, e_max, 250)\n",
    "plt.plot(m, pdf(m, *result.x, e_min, e_max) * np.diff(edges)[0] * len(measured))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least Squares Minimization\n",
    "\n",
    "Lets presume a linear weighted combination of variables:\n",
    "\n",
    "$$\n",
    "f(x)= \\hat{y} =  \\hat{\\beta}_0 + \\sum_{j=1}^p x_j \\hat{\\beta}_j\n",
    "$$\n",
    "where $f:\\mathbb{R}^{p} \\to \\mathbb{R}$.\n",
    "\n",
    "\n",
    "When we include a 1 as the first entry into our sample $x$ e.g. $x = (1, x_1, x_2, \\ldots, x_p)$ we can rewrite\n",
    "$f$ in matrix form\n",
    "\n",
    "$$\n",
    "f(x)= \\hat{y} =  x^T \\mathbf{\\beta}\n",
    "$$\n",
    "\n",
    "where $\\beta = (\\beta_0, \\beta_1, \\beta_2, \\ldots, \\beta_p)$.\n",
    "\n",
    "\n",
    "\n",
    "How do you find those weights? Like before we choose a loss function and try to opimize it.\n",
    "In this case we choose a loss function called the residual sum of squares (RSS).\n",
    "We calculate it over all samples $x_i$ in a matrix $\\mathbf{X}$.\n",
    "\n",
    "$$L(\\beta) = RSS(\\mathbf{\\beta}) = \\sum_{i=1}^N (y_i - x_i^T \\beta)^2 $$\n",
    "\n",
    "Here $x_i$ is a row in $\\mathbf{X}$, hence the transpose.\n",
    "\n",
    "We can now rewrite the loss function in matrix form:\n",
    "\n",
    "\n",
    "$$\n",
    "RSS(\\beta) = (\\mathbf{y} - \\mathbf{X} \\beta)^T (\\mathbf{y} - \\mathbf{X} \\beta )\n",
    "$$\n",
    "\n",
    "Now we optimize the loss function just like we would any other function, by differentiating with respect to $\\beta$ and setting the result equals to zero.\n",
    "\n",
    "$$\n",
    " \\mathbf{X}^T (\\mathbf{y} - \\mathbf{X} \\beta ) \\stackrel{!}{=} 0\n",
    "$$\n",
    "\n",
    "Solving for $\\beta$ leads to\n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T \\mathbf{y}\n",
    "$$\n",
    "\n",
    "\n",
    "We just performed  __Linear Least Squares__ regression. \n",
    "There is also an analytical solution for Linear Least Squares with weights, usually the inverse of the covariance matrix of the data points.\n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = (\\mathbf{X}^T \\mathbf{W} \\mathbf{X})^{-1} \\mathbf{X}^T\\mathbf{W} \\mathbf{y}\n",
    "$$\n",
    "\n",
    "Below we fit some noisy a polynomial data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:26.303140Z",
     "start_time": "2019-04-10T13:18:25.358931Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def f(x, beta=[4.5, 1, 0.75]):\n",
    "    return beta[0]*x**3 + beta[1]*(x-4)**2 + beta[2]*x**4\n",
    "\n",
    "x = np.linspace(-6, 3, 100)\n",
    "data = f(x) + np.random.normal(0, 35, size=x.shape)\n",
    "\n",
    "# build the answer matrix\n",
    "X = np.vstack([x**3, (x - 4)**2, x**4]).T\n",
    "\n",
    "# find the linear least squares solution for the parameters\n",
    "result = optimize.lsq_linear(X, data)\n",
    "\n",
    "\n",
    "plt.plot(x, data, '.')\n",
    "plt.plot(x, f(x), color='gray', label='truth')\n",
    "plt.plot(x, f(x, beta=result.x), color='crimson', label='fit')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-Linear Least-Squares\n",
    "\n",
    "For more general cases where the parameters are not linear there is the more convenient function `curve_fit`\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:27.196526Z",
     "start_time": "2019-04-10T13:18:26.324172Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "\n",
    "def f(x, scale=5, loc=4, ):\n",
    "    return (1/(scale))**2*x**3 + (1/(2*loc))*x**2\n",
    "\n",
    "x = np.linspace(-3,  1.5, 100)\n",
    "data = f(x) + np.random.normal(0, 0.15, size=x.shape)\n",
    "\n",
    "popt, pcov = optimize.curve_fit(f, x, data, p0=[5, 4, ])\n",
    "\n",
    "plt.plot(x, data, '.')\n",
    "plt.plot(x, f(x), color='gray', label='truth')\n",
    "plt.plot(x, f(x, *popt), color='crimson', label='fit')\n",
    "\n",
    "pcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise VI:  Curve Fitting\n",
    "\n",
    "Use `numpy.loadtxt` to load the file called `data/munich_temperatures_average.txt` and fit a model to it using `optimize.curve_fit`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:28.454169Z",
     "start_time": "2019-04-10T13:18:27.211173Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_fit_temperatures.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation is simple and convenient in scipy: The `interp1d` function, when given arrays describing X and Y data, returns and object that behaves like a function that can be called for an arbitrary value of x (in the range covered by X), and it returns the corresponding interpolated y value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:28.552573Z",
     "start_time": "2019-04-10T13:18:28.479360Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import interpolate\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:28.639669Z",
     "start_time": "2019-04-10T13:18:28.564358Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "n = [0, 1, 2, 3, 3.5, 4, 6, 6.5, 7, 9]\n",
    "x = np.linspace(0, 9, 100)\n",
    "\n",
    "y_measured = f(n) \n",
    "\n",
    "linear_interpolation = interpolate.interp1d(n, y_measured)\n",
    "\n",
    "cubic_interpolation = interpolate.interp1d(n, y_measured, kind='cubic')\n",
    "# this returns a new function we can call for values of x\n",
    "cubic_interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:29.692238Z",
     "start_time": "2019-04-10T13:18:28.661870Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.plot(n, y_measured, 's', label='noisy data', color='black')\n",
    "plt.plot(x, f(x), c='gray', lw=2, label='true function')\n",
    "plt.plot(x, linear_interpolation(x), label='linear interpolation')\n",
    "plt.plot(x, cubic_interpolation(x), label='cubic interpolation')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy IO Module\n",
    "\n",
    "## IO and the sound of hydrogen\n",
    "\n",
    "SciPy provides some useful methods for reading file from MatLab, NetCDF or WAV data. Below is an example which uses the `scipy.io.wavfile` module to write a file containing the sound of hydrogen atoms.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:31.393304Z",
     "start_time": "2019-04-10T13:18:29.734027Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy.io import wavfile\n",
    "from scipy import signal\n",
    "\n",
    "rate = 44100 #44.1 khz\n",
    "duration = 10 # in sec\n",
    "\n",
    "def rydberg(n, m):\n",
    "    return (1/(n**2) -1/(m**2)) \n",
    "\n",
    "def lyman(m):\n",
    "    return rydberg(1, m)\n",
    "\n",
    "def balmer(m):\n",
    "    return rydberg(2, m)\n",
    "\n",
    "def sound(frequency, time):\n",
    "    return 2**12 * np.sin(2 * np.pi * 440 * frequency * time)\n",
    "\n",
    "t = np.linspace(0, duration, rate*duration)\n",
    "hydrogen_sound = np.sum([sound(lyman(i), t) + sound(balmer(i), t) for i in range(2, 12)], axis=0)/6\n",
    "\n",
    "plt.plot(t[:10000], hydrogen_sound[:10000])\n",
    "\n",
    "wavfile.write('hydrogen.wav', rate, hydrogen_sound.astype(np.int16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:31.407283Z",
     "start_time": "2019-04-10T13:18:31.396924Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%HTML\n",
    "<audio controls=\"controls\" style=\"width:600px\" >\n",
    "  <source src=\"hydrogen.wav\" type=\"audio/wav\" />\n",
    "</audio>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise IV:  Signal Processing and IO.\n",
    "\n",
    "Load the file `data/synth_sound.wav` and plot a the frequency spectrum vs time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:31.416879Z",
     "start_time": "2019-04-10T13:18:31.411014Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %load solutions/scipy_spectrogram.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Root Finding\n",
    "\n",
    "To find the solution for a function of the form $f(x) = 0$ we can use the `fsolve` or the `root` function. It requires one or more initial guesses: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.396079Z",
     "start_time": "2019-04-10T13:18:31.421484Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 4.5*x**3 + (x-4)**2 + 0.75*x**4 - 20\n",
    "\n",
    "x = np.linspace(-6, 3, 100)\n",
    "\n",
    "plt.plot(x, f(x))\n",
    "plt.axhline(color='black')\n",
    "\n",
    "result = optimize.root(f, [-6, -5, 0, 2])\n",
    "for r in result.x:\n",
    "    plt.axvline(r, color='xkcd:tangerine')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices are often useful in numerical simulations dealing with large systems, if the problem can be described in matrix form where the matrices or vectors mostly contains zeros. Scipy has a good support for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc).\n",
    "\n",
    "There are many possible strategies for storing sparse matrices in an efficient way. Some of the most common are the so-called coordinate form (COO), list of list (LIL) form,  and compressed-sparse column CSC (and row, CSR). Each format has some advantanges and disadvantages. Most computational algorithms (equation solving, matrix-matrix multiplication, etc) can be efficiently implemented using CSR or CSC formats, but they are not so intuitive and not so easy to initialize. So often a sparse matrix is initially created in COO or LIL format (where we can efficiently add elements to the sparse matrix data), and then converted to CSC or CSR before used in real calcalations.\n",
    "\n",
    "For more information about these sparse formats, see e.g. http://en.wikipedia.org/wiki/Sparse_matrix\n",
    "\n",
    "When we create a sparse matrix we have to choose which format it should be stored in. For example, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.470143Z",
     "start_time": "2019-04-10T13:18:32.449574Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from scipy import sparse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.547960Z",
     "start_time": "2019-04-10T13:18:32.511596Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# dense matrix\n",
    "M = np.array([[1,0,0,0], [0,3,0,0], [0,1,1,0], [1,0,0,1]]); M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.673802Z",
     "start_time": "2019-04-10T13:18:32.580832Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# convert from dense to sparse\n",
    "A = sparse.csr_matrix(M)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.723402Z",
     "start_time": "2019-04-10T13:18:32.706736Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# convert from sparse to dense\n",
    "A.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More efficient way to create sparse matrices: create an empty matrix and populate with using matrix indexing (avoids creating a potentially large dense matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.786928Z",
     "start_time": "2019-04-10T13:18:32.747482Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A = sparse.lil_matrix((4,4)) # empty 4x4 sparse matrix\n",
    "A[0,0] = 1\n",
    "A[1,1] = 3\n",
    "A[2,2] = A[2,1] = 1\n",
    "A[3,3] = A[3,0] = 1\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.839446Z",
     "start_time": "2019-04-10T13:18:32.813500Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A.todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting between different sparse matrix formats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.868579Z",
     "start_time": "2019-04-10T13:18:32.849780Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.928040Z",
     "start_time": "2019-04-10T13:18:32.873701Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A = sparse.csr_matrix(A); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.963984Z",
     "start_time": "2019-04-10T13:18:32.951025Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A = sparse.csc_matrix(A); A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute with sparse matrices like with dense matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:32.978513Z",
     "start_time": "2019-04-10T13:18:32.968062Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.004014Z",
     "start_time": "2019-04-10T13:18:32.993027Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "(A * A).todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.030237Z",
     "start_time": "2019-04-10T13:18:33.014899Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A.todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.050696Z",
     "start_time": "2019-04-10T13:18:33.037186Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "A.dot(A).todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.072953Z",
     "start_time": "2019-04-10T13:18:33.055013Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "v = np.array([1,2,3,4])[:,np.newaxis] \n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.088442Z",
     "start_time": "2019-04-10T13:18:33.078162Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# sparse matrix - dense vector multiplication\n",
    "A * v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-10T13:18:33.111929Z",
     "start_time": "2019-04-10T13:18:33.101351Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# same result with dense matrix - dense vector multiplcation\n",
    "A.todense() * v"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}