library_scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "# Table of Contents\n",
    "  - [Scipy](#Scipy)\n",
    "  - [Table of Contents](#Table-of-Contents)\n",
    "    - [0. References](#0.-References)\n",
    "    - [1. Introduction](#1.-Introduction)\n",
    "      - [1.1 Pre-requisites](#1.1-Pre-requisites)\n",
    "    - [2. Use case: Optimization](#2.-Use-case:-Optimization)\n",
    "      - [2.1. Minimizing the value of a function](#2.1.-Minimizing-the-value-of-a-function)\n",
    "      - [2.2. Example: determining the gravity acceleration](#2.2.-Example:-determining-the-gravity-acceleration)\n",
    "      - [2.3. Exercise: optimization 🌶️](#2.3.-Exercise:-optimization-🌶️)\n",
    "    - [3. Use case: spatial data](#3.-Use-case:-spatial-data)\n",
    "      - [3.1. Example: nearest neighbor search](#3.1.-Example:-nearest-neighbor-search)\n",
    "      - [3.2. Example: bacteria statistics](#3.2.-Example:-bacteria-statistics)\n",
    "      - [3.4. Root finding](#3.4.-Root-finding)\n",
    "      - [3.5. Exercise: root finding 🌶️](#3.5.-Exercise:-root-finding-🌶️)\n",
    "    - [4. Linear algebra ](#4.-Linear-algebra)\n",
    "      - [4.1. Example: solving a simple linear system](#4.1.-Example:-solving-a-simple-linear-system)\n",
    "      - [4.2. Triangular matrices](#4.2.-Triangular-matrices)\n",
    "      - [4.3. Toeplitz matrices](#4.3.-Toeplitz-matrices)\n",
    "      - [4.4. Eigenvalue problems](#4.4.-Eigenvalue-problems)\n",
    "      - [4.5. Matrix decomposition: $LU$ decomposition](#4.5.-Matrix-decomposition:-$LU$-decomposition)\n",
    "      - [4.6. Matrix decomposition: Choleski decomposition](#4.6.-Matrix-decomposition:-Choleski-decomposition)\n",
    "      - [4.7. Sparse matrices](#4.7.-Sparse-matrices)\n",
    "      - [4.8. Exercise: $LU$ decomposition 🌶️](#4.8.-Exercise:-$LU$-decomposition-🌶️)\n",
    "      - [4.9. Exercise: the singular value decomposition (SVD) 🌶️](#4.9.-Exercise:-the-singular-value-decomposition-(SVD)-🌶️)\n",
    "    - [5. Use case: interpolation](#5.-Use-case:-interpolation)\n",
    "      - [5.1. Curve fitting](#5.1.-Curve-fitting)\n",
    "      - [5.2. Simple harmonic motion](#5.2.-Simple-harmonic-motion)\n",
    "      - [5.3. Exercise: curve fitting 🌶️](#5.3.-Exercise:-curve-fitting-🌶️)\n",
    "    - [6. Fun with functions: the Legendre polynomials and the Bessel functions](#6.-Fun-with-functions:-the-Legendre-polynomials-and-the-Bessel-functions)\n",
    "    - [Mathematical analysis](#Mathematical-analysis)\n",
    "    - [7. Differentiation](#7.-Differentiation)\n",
    "    - [8. Integration](#8.-Integration)\n",
    "      - [8.1. Exercise: double integration 🌶️](#8.1.-Exercise:-double-integration-🌶️)\n",
    "    - [9. Differential equations](#9.-Differential-equations)\n",
    "      - [9.1. First-order ODEs](#9.1.-First-order-ODEs)\n",
    "      - [9.2. Exercise: first-order ODEs](#9.2.-Exercise:-first-order-ODEs)\n",
    "      - [9.3. Coupled first-order ODEs](#9.3.-Coupled-first-order-ODEs)\n",
    "      - [9.4. Second-order ODEs](#9.4.-Second-order-ODEs)\n",
    "    - [10. Fourier transforms](#10.-Fourier-transforms)\n",
    "      - [11. Statistics](#11.-Statistics)\n",
    "        - [11.1. The $\\beta$ distribution](#11.1.-The-$\\beta$-distribution)\n",
    "        - [11.2. the Gaussian distribution](#11.2.-the-Gaussian-distribution)\n",
    "        - [11.3. The multinomial distribution](#11.3.-The-multinomial-distribution)\n",
    "      - [12. Problems (freehand)](#12.-Problems-(freehand))\n",
    "        - [12.1. Mechanical work 🌶️🌶️](#12.1.-Mechanical-work-🌶️🌶️)\n",
    "        - [12.2. Newton's law of cooling 🌶️🌶️](#12.2.-Newton's-law-of-cooling-🌶️🌶️)\n",
    "        - [12.3. The Lennard-Jones potential 🌶️🌶️](#12.3.-The-Lennard-Jones-potential-🌶️🌶️)\n",
    "        - [12.4. Random number generation from your own distribution 🌶️🌶️🌶️](#12.4.-Random-number-generation-from-your-own-distribution-🌶️🌶️🌶️)\n",
    "        - [12.5. The finite square well 🌶️🌶️🌶️🌶️](#12.5.-The-finite-square-well-🌶️🌶️🌶️🌶️)\n",
    "        - [12.6. the quintic polynomial 🌶️🌶️](#12.6.-the-quintic-polynomial-🌶️🌶️)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 0. References"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "Additional material for learning scipy, from the documentation, articles, videos or examples:\n",
    "- The scipy [user guide](https://docs.scipy.org/doc/scipy/tutorial/index.html#user-guide).\n",
    "- A good wikipedia article on [optimization](https://en.wikipedia.org/wiki/Mathematical_optimization).\n",
    "- The scipy [cookbook](https://scipy-cookbook.readthedocs.io/index.html) with several user-contributed examples of uses for scipy.\n",
    "- A good [video explanation](https://www.youtube.com/watch?v=Glp7THUpGow) of k-d trees\n",
    "- YouTube channel: Mr. P. Solver (from whom much of the second half of this tutorial is taken)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "Scipy is a python module that provides implementation of many algorithms commonly used in scientific computing. These covers many topics, including:\n",
    "\n",
    "- Optimization\n",
    "- Numerical integration\n",
    "- Linear algebra\n",
    "- Fourier transforms\n",
    "- Spatial data processing\n",
    "- Sparse matrices\n",
    "\n",
    "And many more. \n",
    "Because of the wide variety of algorithms, we will structure this chapter as a series of use cases where the use of scipy would make it easier to solve the problem at hand.\n",
    "If you want to discover more algorithms offered by scipy, please check the [documentation](https://docs.scipy.org/doc/scipy/tutorial/index.html#user-guide).\n",
    "\n",
    "\n",
    "### 1.1 Pre-requisites\n",
    "\n",
    "Most of the algorithms in scipy operate on (multidimensional) arrays of numerical values. \n",
    "For technical reasons, these arrays aren't implemented with the built-in python list datatype but using [**numpy arrays**](https://numpy.org/doc/stable/user/absolute_beginners.html). \n",
    "\n",
    "These objects offer better memory and computation performance when operating on large arrays as compared to python lists. \n",
    "If you want to understand this chapter in-depth and be able to solve the exercises, we advice you first read the chapter on numpy here.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## 2. Use case: Optimization\n",
    "\n",
    "A common problem in the analysis of scientific data is **optimization**.  \n",
    "In this context, we mean a restricted definition, where given a function and a set of data we want to determine the value of one or more parameters in order to minimize the distance between the function and the data.\n",
    "\n",
    "### 2.1. Minimizing the value of a function\n",
    "Let's see an example to understand what we mean: given a quadratic function `f(x)`, we want to determine the value of `x` that minimizes its value.\n",
    "We start by defining our function `f`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    \"\"\"\n",
    "    This is the function we want to minimise\n",
    "    \"\"\"\n",
    "    return x**2 + 4*x + 4 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "Now we plot the function using `matplotlib`, a module for plotting described in another module of our course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "# Generate x values for plotting\n",
    "x_values = np.linspace(-10, 2, 100)\n",
    "\n",
    "# Plot the quadratic function\n",
    "plt.plot(x_values, f(x_values), label='Quadratic Function')\n",
    "plt.title('Quadratic Function')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "Visually (and algebraically) we can determine the minimum of this function to be at -2. \n",
    "\n",
    "We now want to determine the minimum of this function **numerically** using the algorithms offered by scipy. \n",
    "Naturally, we would not do this for a quadratic functions, where we can determine the derivatives in a closed form and hence find the exact value of the minimum.\n",
    "However, the quadratic function is just an example, in reality we would often encounter functions that are too complex to handle algebraically or where the derivate does not exist in a closed form. \n",
    "In those cases, we can only determine minima (or solve other optimization problems) using numerical approximation with algorithms like the one scipy offers.\n",
    "\n",
    "Fortunately, scipy offers the `minimize` [function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html), which can directly estimate the minimum of a function of one argument.\n",
    "The function takes the function to minimize as a first argument and an initial guess of the minimum as a second argument.\n",
    "\n",
    "The choice of the initial guess is very important for the success of the minimization: if the function is not **convex**, a wrong choice of initial value will not lead to the true minimum but will result in the algorithm being stuck in a local minimum. \n",
    "It is not the purpose of this course to discuss these topics in-detail, if you need more information please consult a book on numerical methods or on optimization.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "initial_guess = 0  # Initial guess for the minimum\n",
    "result = minimize(f, initial_guess)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "As you can see from the output, running `minimze` results in an object with several attributes. \n",
    "Of interest for us are primarily the following:\n",
    "\n",
    "- `success` is a boolean flag which is `True` if the algorithm could find a (local) minimum\n",
    "- `fun` is the value of the function at the location of the estimated minimum\n",
    "- `x` is the estimated value of `x` that minimizes the function.\n",
    "\n",
    "In this case, the algorithm correctly found the minimum of `f` at -2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "### 2.2. Example: determining the gravity acceleration\n",
    "Suppose for example we want to determine the gravitational acceleration $g$ with an experiment.\n",
    "To do so, we time the fall of an object from an height $h$ until it touches the ground.  This results in a measurement $t_f$ (the total fall time). In order to obtain enough data, we repeat the experiment for a set of heights $h$.\n",
    "\n",
    "This defines the following relationship:\n",
    "\n",
    "$h(t, g) = 1/2 * g t^2$\n",
    "\n",
    "Because we measure $t$ and vary $h$, we express the problem as:\n",
    "\n",
    "$t(h, g) = \\sqrt{2 h / g}$\n",
    "\n",
    "Since we can't perform the experiment ourselves, let's first define a function `simulate_fall` to generate the fall time as a function of $h$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.typing import ArrayLike\n",
    "from numpy.random import randn\n",
    "def simulate_fall(h: ArrayLike) -> ArrayLike:\n",
    "    \"\"\"\n",
    "    Simulate the fall of an object from a height h\n",
    "    \"\"\"\n",
    "    g = 9.81  # Acceleration due to gravity\n",
    "    return np.sqrt(2 * h / g) + 1e-3 * randn(len(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Notice that we add a random number sampled from a Gaussian distribution with zero mean and standard deviation of 1e-4 seconds to the simulated values. \n",
    "We do this to make the measurement more realistic by simulating the effect of measurement noise and imprecisions.\n",
    "Now we generate 20 fall experiments with heights equally spaced between 0 and 1 m and plot the graph of times against the starting height. \n",
    "To make our experiment more statistically representative, we repeat it three times\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rep = 4\n",
    "height_steps = np.linspace(0, 1, 20)\n",
    "heights = height_steps.repeat(n_rep)\n",
    "fall_times = simulate_fall(heights)\n",
    "\n",
    "\n",
    "# Plot the fall times against the heights\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(heights, fall_times, label='Simulated fall times', s=2)\n",
    "ax.set_xlabel('Height (m)')\n",
    "ax.set_ylabel('Fall time (s)')\n",
    "ax.set_title('Simulated fall times against height')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "At this point we have all the data we need for estimating $g$ from our measurements. \n",
    "In this case, the array `fall_times` plays the role of $\\mathbf{y}$ and the array `heights` corresponds to $\\mathbf{x}$.\n",
    "We can use the function [`scipy.optimize.curve_fit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit) to solve these types of problems.\n",
    "The function  takes the following arguments:\n",
    "\n",
    "- `f`, the function to fit. Its first argument must be `x`, the other arguments correspond to the curve parameters in order of listing\n",
    "- `xdata` corresponds to the vector $x$\n",
    "- `ydata` corresponds to the vector $y$\n",
    "\n",
    "\n",
    "Following this, let's define our function to fit, `f`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(h: ArrayLike, g: float) -> ArrayLike:\n",
    "    return np.sqrt(2 * h / g)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "Now we can call the function:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# Fit the function to the data\n",
    "g_est, pcov = curve_fit(f, heights, fall_times)\n",
    "print(g_est)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21",
   "metadata": {},
   "source": [
    "Success! As you can see, the estimated value of $g$ is quite close to the real value. The discrepancy is in the order of magnitude of the standard deviation of the noise we add to the simulated measurements.\n",
    "\n",
    "We can now plot the fitted curve and superimpose it to the \"measurements\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the fall times against the heights\n",
    "estimated_fall_times = f(heights, g_est[0])\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(heights, fall_times, s=2, label='Simulated fall times')\n",
    "ax.plot(heights, estimated_fall_times, label='Estimated fall times', color='orange')\n",
    "ax.set_xlabel('Height (m)')\n",
    "ax.set_ylabel('Fall time (s)')\n",
    "ax.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "You can see that the curve with the estimated fall times given the estimated $g$ tracks the \"measurements\" very well.\n",
    "With this example, we conclude the first use case."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "### 2.3. Exercise: optimization 🌶️\n",
    "\n",
    "Define a quadratic function $$f(x) = 3x^2 - 5x - 4$$ under the Python function `solution_minimize()` below. Then use `scipy.optimize.minimize` to:\n",
    "1. obtain the extremum\n",
    "2. obtain the value of the function at this extremum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def solution_minimize(): # Don't change the function name or signature\n",
    "    # 1. TODO:  Define the function to minimize\n",
    "    \n",
    "    # 2. TODO:  Define the initial guess\n",
    "\n",
    "    # 3. TODO:  Call the minimize function\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27",
   "metadata": {},
   "source": [
    "## 3. Use case: spatial data\n",
    "Another place where scipy shines is the processing of **spatial data**. \n",
    "This rather generic terms is used to describe data on a plane that represent the locations and shapes of objects, for example maps.\n",
    "Scipy offers some useful functions to perform basic processing on this data, although for more complex cases there are specific libraries to manipulate and process geospatial data such as [GDAL](https://gdal.org/) that support advanced features like conversion between different geographical coordinate systems, conversion of file format and much more.\n",
    "\n",
    "### 3.1. Example: nearest neighbor search\n",
    "This use case is inspired by a real question a researcher at Empa asked us. \n",
    "They approached us for help in analyzing microscopy imaging: they wanted to determine the distribution of the relative orientation angle of bacteria as a function of their distance. This means performing the following steps:\n",
    "\n",
    "1. Segment an image and extract the outline of each bacteria. \n",
    "2. From the outline, we can estimate the direction of the bacteria with respect to the image axis as well as the center point of the bacteria\n",
    "3. For each bacteria, iterate over all other bacteria forming pairs (bacteria1, bacteria2)\n",
    "4. For each pair (bacteria1, bacteria2) compute the relative orientation and the distance and store it\n",
    "\n",
    "\n",
    "The user quickly realized that for a large number of cells in an image this algorithm becomes prohibitively slow to compute. \n",
    "Because we perform a pairwise comparison of points, the time (the number of steps) needed to compute the the relative distances and orientations grows with the square of the number of cells. \n",
    "However, the user knew that they didn't want to compute **all** distances and orientations but just those **smaller than a certain threshold**.\n",
    "Thanks to this knowledge and the use of scipy, we could find a method to make their computation much faster.\n",
    "\n",
    "To do this, we used a special data structure called a [**k-d Tree**](https://en.wikipedia.org/wiki/K-d_tree); this data structure represents a set of k-dimensional points as a tree where each node of the tree splits the space in exactly two halves along one of the coordinates. Thanks to this property, we can very quickly find the nearest points to a given point because at every step in the tree we directly eliminate half of the space to search in. For a better visualization of this idea, please check this [video](https://www.youtube.com/watch?v=Glp7THUpGow).\n",
    "\n",
    "\n",
    "Let's see how to use this structure with a simplified example. \n",
    "\n",
    "### 3.2. Example: bacteria statistics\n",
    "Imagine we have a plate on which we grew a set of bacteria. These have different orientation and they tend to cluster together:  bacteria near to each other tend to have the same orientation; moreover their size is governed by the expression of a protein that we determined by their fluorescence; the more they express that protein the larger they grow. As for the orientation, their sizes tend to cluster. \n",
    "\n",
    "Our goal here is to determine the **variogram** of angles and size; this is the average squared difference of size (and angles) between pairs of bacteria as a function of their separation, binned by their separation distance.\n",
    "\n",
    "We also imagine that someone pre-processed the microscopy imaging data and gave us a 2-D numpy array of shape $n_bacteria \\times 4$, where the four additional dimensions are the two plane coordinates, the angle and the protein expression.\n",
    "\n",
    "\n",
    "To simulate this data, we create a set of points in the plane, to each of the point we attach a number between 0 and $2\\pi$ that represent the orientation of the bacteria with respect to the coordinate system. Similarly we attach a number between 0 and 1 to represent the amount of protein they expressed. \n",
    "In order to make these distributions more interesting, we use the `gstools` package that allows to create **spatial random fields** with the desired correlation structure and statistics.\n",
    "\n",
    "In the following cell, we generate the test data and display it:\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#from numpy.typing import NDArray\n",
    "import matplotlib.pyplot as plt\n",
    "import gstools as gs\n",
    "import numpy as np\n",
    "from matplotlib.axes import Axes\n",
    "from matplotlib.patches import Ellipse\n",
    "from matplotlib.collections import PatchCollection\n",
    "from typing import Tuple\n",
    "NDim = 5\n",
    "NSamples = 1000\n",
    "ResultArray = np.ndarray(shape=(NDim, NSamples), dtype=np.float64)\n",
    "def generate_random_field(n: int): # -> Tuple[ResultArray, NDArray, NDArray]:\n",
    "    np.random.seed(42)\n",
    "    grid_size = 100\n",
    "    seed = 42\n",
    "    # Define the spatial domain\n",
    "    x = np.linspace(0, 1, grid_size)\n",
    "    y = np.linspace(0, 1, grid_size)\n",
    "\n",
    "    # Create a Gaussian variogram model for the angle distribution\n",
    "    angle_model = gs.Gaussian(dim=2, var=2 * np.pi, len_scale=0.1, angles=3, nugget=0.01)\n",
    "    # Create a random field generator\n",
    "    angle_rf = gs.SRF(angle_model, seed=seed, mean=np.pi/2)\n",
    "    #Protein expression model\n",
    "    protein_model = gs.Gaussian(dim=2, var=0.5, len_scale=0.3, nugget=0.01)\n",
    "    # Create a random field generator\n",
    "    protein_rf = gs.SRF(protein_model, seed=seed, mean=0.1)\n",
    "\n",
    "    # Generate random points based on the field\n",
    "    num_points = n\n",
    "    \n",
    "    i = np.random.choice(x.shape[0], num_points)\n",
    "    j = np.random.choice(y.shape[0], num_points)\n",
    "\n",
    "    \n",
    "\n",
    "    # Generate a realization of the random field\n",
    "    angles = angle_rf(np.stack([x[i], y[j]]))\n",
    "    protein = protein_rf(np.stack([x[i], y[j]]))\n",
    "\n",
    "    angles_full = angle_rf.structured([x, y])\n",
    "    protein_full = protein_rf.structured([x, y])\n",
    "\n",
    "\n",
    "    bacteria = np.stack([x[i], y[j], angles, protein])\n",
    "\n",
    "    return bacteria, angles_full, protein_full\n",
    "    \n",
    "n_bacteria = 1000\n",
    "bacteria, angles, protein = generate_random_field(n_bacteria)\n",
    "\n",
    "def draw_bacteria(ell: ResultArray, ax: Axes):\n",
    "    ratio = 4\n",
    "    scale = 0.01\n",
    "    ells = [Ellipse(xy=ell[0:2, i], width=ratio * scale * np.abs(ell[3, i]), height=scale * ell[3, i], angle=np.rad2deg(ell[2, i])) for i in range(ell.shape[1])]   \n",
    "    ax.add_collection(PatchCollection(ells, color='green', alpha=0.8))          \n",
    "f, [ax_prot, ax_angle, ax_bacteria] = plt.subplots(1, 3, figsize=(15, 5), sharey=True, sharex=True)\n",
    "draw_bacteria(bacteria, ax_bacteria) \n",
    "ax_bacteria.set_title('Bacteria')\n",
    "ax_bacteria.set_xlim([0, 1])\n",
    "ax_bacteria.set_ylim([0, 1])\n",
    "ax_prot.imshow(protein, origin='upper', extent=[0,1, 0, 1])\n",
    "ax_prot.set_title('Protein')\n",
    "im = ax_angle.imshow(np.rad2deg(angles), origin='lower', extent=[0, 1, 0, 1])\n",
    "im.set_clim([0, 360])\n",
    "ax_angle.set_title('Angles')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "Now that we are able to generate and display the data, let's try and compute the variogram.\n",
    "Because constructing all pair of bacterias will require $n^2$ steps, we use a KD-Tree instead. \n",
    "To begin with, we intialise a KD Tree with the coordinates of all bacteria. \n",
    "Once we have the data structure, we use the function `query_pairs` to find a list of all pairs of indices `(i, j)` of bacterias whose distance ist less than 0.6\n",
    "Thanks to the KD-Tree structure, this operation is much faster than directly computing all distances and filtering them.\n",
    "Once we have the list, we can finally compute the distances and the variogram  using some array manipulation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we can analyse the data\n",
    "import scipy.spatial as spatial\n",
    "import scipy.stats as stats\n",
    "#Construct a KDTree\n",
    "tree =  spatial.KDTree(bacteria[0:2, :].T)\n",
    "\n",
    "#Find the nearest neighbours\n",
    "pairs = tree.query_pairs(r=0.6, output_type='ndarray')\n",
    "distances = np.zeros((pairs.shape[0]))\n",
    "angles = np.zeros((pairs.shape[0]))\n",
    "proteins = np.zeros((pairs.shape[0]))\n",
    "for i in pairs:\n",
    "    distances[i] = np.linalg.norm(bacteria[0:2, i[0]] - bacteria[0:2, i[1]])\n",
    "    angles[i] = np.mod(bacteria[2, i[0]] - bacteria[2, i[1]], np.pi)\n",
    "    proteins[i] = np.abs(bacteria[3, i[0]] - bacteria[3, i[1]])\n",
    "\n",
    "#Compute the mean angle for each distance\n",
    "#Specify bin edges for distances and angles\n",
    "distance_bins = np.linspace(distances.min(), 0.3, 10)\n",
    "angle_bins = np.linspace(0, 2*np.pi, 20)\n",
    "\n",
    "def binned_stat(x, y, bins): # -> Tuple[ArrayLike, ArrayLike, ArrayLike]:\n",
    "    \"\"\"\n",
    "    Compute a statistic on y for each bin in x\n",
    "    \"\"\"\n",
    "    x_bins = np.digitize(x, bins)\n",
    "    bin_means = np.zeros((len(bins)-1))\n",
    "    bin_stds = np.zeros((len(bins)-1))\n",
    "    for i in range(len(bins)-1):\n",
    "        bin_means[i] = np.mean(y[x_bins == i])\n",
    "        bin_stds[i] = np.std(y[x_bins == i])\n",
    "    return bins, bin_means, bin_stds\n",
    "\n",
    "# Compute 2D histogram (distance, angle)\n",
    "distance_bins, angle_mean, angle_mean = binned_stat(distances, angles, distance_bins)\n",
    "distance_bins, protein_mean, protein_std = binned_stat(distances, proteins, distance_bins)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {},
   "source": [
    "Now that we computed the variograms for angle difference and protein expression, we can plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, (ax_vario_angle, ax_vario_prot) = plt.subplots(2, 1, figsize=(10, 10), sharex=True)\n",
    "ax_vario_angle.plot(distance_bins[:-1], angle_mean)\n",
    "ax_vario_angle.set_xlabel('Distance')\n",
    "ax_vario_angle.set_ylabel('Mean angle difference')\n",
    "ax_vario_prot.plot(distance_bins[:-1], protein_mean)\n",
    "ax_vario_prot.set_xlabel('Distance')\n",
    "ax_vario_prot.set_ylabel('Mean protein difference')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33",
   "metadata": {},
   "source": [
    "### 3.4. Root finding\n",
    "\n",
    "One of the most fundamental tasks of elementary algebra is to find the roots of a function, and the applications of root finding extend very much far beyond those cases which you probably saw during your schooltime years. \n",
    "\n",
    "For sake of simplicity, we first examine the case of a univariate root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return x**3 -3*x + 1\n",
    "\n",
    "x = np.linspace(-3,3,100)\n",
    "plt.plot(x, f(x))\n",
    "plt.axhline(0, color='black')\n",
    "plt.axvline(0, color='black')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35",
   "metadata": {},
   "source": [
    "We recall that the fundamental theorem of algebra states that the above-defined $f(x)$ has exactly three roots in the complex numbers. It turns out that all the solutions of this $f(x)$ are real, as we see from the above plot.\n",
    "\n",
    "In order to numerically find the value of these roots, we have to provide the solver with an initial guess, near which we expect the respective root to lie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import root\n",
    "\n",
    "first_root = root(f, -2)\n",
    "first_root"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37",
   "metadata": {},
   "source": [
    "Notice that the returned object from the function `root()` contains several elements. We just want the `.x` attribute. To get all three roots, we can simply just call them as a tuple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_root, second_root, third_root = root(f, -2), root(f, 0.5), root(f, 1.5)\n",
    "first_root.x, second_root.x, third_root.x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39",
   "metadata": {},
   "source": [
    "### 3.5. Exercise: root finding 🌶️\n",
    "\n",
    "Use SciPy to find the roots in the real numbers of $$ x^3 - 6x^2 + 4x + 12 = 0. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "from scipy.optimize import root\n",
    "\n",
    "def solution_root_finder(): # Don't change the function name or signature\n",
    "    # 1. TODO: Define the function f(x)\n",
    "    # 2. TODO: Find the roots of f(x) using the root() function\n",
    "    # 3. TODO: Return the roots as a tuple\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "## 4. Linear algebra \n",
    "Another typical use case for scipy is **linear algebra**. \n",
    "Here we don't cover all the basics of linear algebra and we don't try to be exhaustive in our coverage;  instead we will show an example that is commonly encountered in a similar form by many scientists: solving a system of linear equations.\n",
    "\n",
    "### 4.1. Example: solving a simple linear system\n",
    "Consider the system of two linear equations:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "2x + y &= 5 \\\\\n",
    "-4x + 6y &= 2\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "As we learn in linear algebra, we can express these system as a matrix-vector product:\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{x} = \\mathbf{b}$$\n",
    "\n",
    "Where:\n",
    "$$\n",
    "\\mathbf{A} = \\begin{bmatrix} 2 & 1 \\\\ -4 & 6 \\end{bmatrix}, \\quad\n",
    "\\mathbf{x} = \\begin{bmatrix} x \\\\ y \\end{bmatrix}, \\quad\n",
    "\\mathbf{b} = \\begin{bmatrix} 5 \\\\ 2 \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "we know that we can solve this system by Gaussian elimination. \n",
    "This would be fine for such a small system, however as the number of equations grow, the time it takes to manipulate the matrix $\\mathbf{A}$ makes it impossible for a person to find a solution manually.\n",
    "Fortunately, Gaussian elimination is a well-known algorithm and efficient implementations are available for most programming languages. \n",
    "One such implementation can be accessed by using [`scipy.lingalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html).\n",
    "We now try to solve the system using this function; as a first step we need to express the various matrices and vectors as numpy arrays.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg\n",
    "A = np.array([[2, 1], [4, 6]])\n",
    "b = np.array([5 ,2])\n",
    "x_hat = scipy.linalg.solve(A, b)\n",
    "print(f\"The solution is {x_hat}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "We now can verify if the solution is correct by computing $\\mathbf{A}\\hat{\\mathbf{x}}$, where $\\hat{\\mathbf{x}}$ is the solution obtained using `solve`. \n",
    "When we use numpy arrays, we use the `@` operator to express the matrix-vector product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "A @ x_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46",
   "metadata": {},
   "source": [
    "Indeed, we obtain the vector that (visually) looks like $\\mathbf{b}$.\n",
    "To verify if indeed the solutions are close, we can look at the distance between $\\mathbf{A}\\hat{\\mathbf{x}}$ und $\\mathbf{b}$ using the `scipy.linalg.norm` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.linalg.norm(A @ x_hat - b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48",
   "metadata": {},
   "source": [
    "### 4.2. Triangular matrices\n",
    "\n",
    "Triangular matrices can best be illustrated as having being of the following form:\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{cccc}\n",
    "3 & 0 & 0 & 0 \\\\\n",
    "2 & 1 & 0 & 0 \\\\\n",
    "1 & 0 & 1 & 0 \\\\\n",
    "1 & 1 & 1 & 1\n",
    "\\end{array}\n",
    "\\right] \\vec{x} =\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "4 \\\\\n",
    "2 \\\\\n",
    "4 \\\\\n",
    "2 \n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "We can solve them quite easily in SciPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import solve_triangular\n",
    "a = np.array([[3, 0, 0, 0],\n",
    "              [2, 1, 0, 0],\n",
    "              [1, 0, 1, 0],\n",
    "              [1, 1, 1, 1]])\n",
    "\n",
    "b = np.array([4,2,4,2])\n",
    "x = solve_triangular(a,b,lower=True)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50",
   "metadata": {},
   "source": [
    "### 4.3. Toeplitz matrices\n",
    "\n",
    "**Toeplitz matrices** are matrices whose diagonal entries are identical\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{cccc}\n",
    "1 & -1 & 2 & 3 \\\\\n",
    "3 & 1 & -1 & 2 \\\\\n",
    "6 & 3 & 1 & -1 \\\\\n",
    "10 & 6 & 3 & 1\n",
    "\\end{array}\n",
    "\\right] \\vec{x} =\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "1 \\\\\n",
    "2 \\\\\n",
    "2 \\\\\n",
    "5 \n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "- Toeplitz matrices are specified **uniquely** by their first column and their first row!!\n",
    "\n",
    "We can solve them using SciPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import solve_toeplitz, toeplitz\n",
    "\n",
    "c = np.array([1,3,6,10]) # first column\n",
    "r = np.array([1,-1,-2,-3]) # first row\n",
    "b = np.array([1,2,2,5])\n",
    "x = solve_toeplitz((c,r),b)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52",
   "metadata": {},
   "source": [
    "### 4.4. Eigenvalue problems\n",
    "\n",
    "Eigenvalue problems can be solved using `numpy`, so here we focus on particular cases for optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import eigh_tridiagonal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54",
   "metadata": {},
   "source": [
    "Consider the tridiagonal matrix\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{cccc}\n",
    "-3 & -1 & 0 & 0 \\\\\n",
    "-1 & 3 & -1 & 0 \\\\\n",
    "0 & -1 & 3 & -1 \\\\\n",
    "0 & 0 & -1 & 3\n",
    "\\end{array}\n",
    "\\right] \\vec{x} = \\lambda \\vec{x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55",
   "metadata": {},
   "outputs": [],
   "source": [
    "main_diag = 3*np.ones(4)\n",
    "off_diag = -1*np.ones(3)\n",
    "eigenvalues, eigenvectors = eigh_tridiagonal(main_diag, off_diag)\n",
    "eigenvectors.T[0] # gives the first eigenvector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56",
   "metadata": {},
   "source": [
    "Let's check that this is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.diag(main_diag) + np.diag(off_diag, k=1) + np.diag(off_diag, k=-1)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58",
   "metadata": {},
   "outputs": [],
   "source": [
    "A@eigenvectors.T[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalues[1]*eigenvectors.T[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60",
   "metadata": {},
   "source": [
    "### 4.5. Matrix decomposition: $LU$ decomposition\n",
    "\n",
    "**LU decomposition:** $A=PLU$ where $P$ is a permutation matrix, $L$ is a lower triangular matrix, and $U$ is an upper triangular matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import lu\n",
    "A = np.array([[2,5,7,8],[5,2,2,8],[7,5,5,6],[5,4,4,8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "P, L, U = lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63",
   "metadata": {},
   "outputs": [],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64",
   "metadata": {},
   "outputs": [],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65",
   "metadata": {},
   "source": [
    "### 4.6. Matrix decomposition: Choleski decomposition\n",
    "\n",
    "**Choleski decomposition:** find matrix $C$ such that $A=CC^T$, where $A$ is a matrix that we specify. An important theorem from linear algebra states that the Cholesky decomposition is unique when the matrix to be decomposed is positive definite!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import cholesky\n",
    "A = np.array([[1,0.2], [0.2,1]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67",
   "metadata": {},
   "outputs": [],
   "source": [
    "C = cholesky(A)\n",
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68",
   "metadata": {},
   "source": [
    "### 4.7. Sparse matrices\n",
    "\n",
    "**Sparce matrices** are matrices that contain lots of zeros (so lots of space can be reduced).\n",
    "\n",
    "**a useful example**:\n",
    "\n",
    "the second derivative of $f(x_i) := f_i$ is approximated using the method of finite differences as\n",
    "\n",
    "$$\n",
    "\\frac{d^2f_i}{dx^2} \\approx \\frac{f_{i+1} + f_{i-1} -2f_i}{\\Delta x^2}.\n",
    "$$\n",
    "\n",
    "Suppose we have $f_0, \\dots, f_4$ and $f_0=f_4=0$ (boundary conditions). Then the second derivative is approximated as\n",
    "\n",
    "$$\n",
    "D\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "f_1 \\\\\n",
    "f_2 \\\\\n",
    "f_3\n",
    "\\end{array}\n",
    "\\right], \\text{  where }\n",
    "D = \\frac{1}{\\Delta x^2}\n",
    "\\left[\n",
    "\\begin{array}{ccccc}\n",
    "-2 & 1 & 0 & 0 & 0 \\\\\n",
    "1 & -2 & 1 & 0 & 0 \\\\\n",
    "0 & 1 & -2 & 1 & 0 \\\\\n",
    "0 & 0 & 1 & -2 & 1 \\\\\n",
    "0 & 0 & 0 & 1 & -2 \n",
    "\\end{array}\n",
    "\\right].\n",
    "$$\n",
    "\n",
    "In two dimensions, our function can be discretized on a grid:\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccccc}\n",
    "0 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & f_{11} & f_{12} & f_{13} & 0 \\\\\n",
    "0 & f_{21} & f_{22} & f_{23} & 0 \\\\\n",
    "0 & f_{31} & f_{32} & f_{33} & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 0\n",
    "\\end{array}\n",
    "\\right],\n",
    "$$\n",
    "\n",
    "but when doing these sorts of problems, it's better to store information in a column vector:\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccccc}\n",
    "0 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & f_{11} & f_{12} & f_{13} & 0 \\\\\n",
    "0 & f_{21} & f_{22} & f_{23} & 0 \\\\\n",
    "0 & f_{31} & f_{32} & f_{33} & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 0\n",
    "\\end{array}\n",
    "\\right] \\rightarrow \n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "f_{11} \\\\\n",
    "f_{12} \\\\\n",
    "f_{13} \\\\\n",
    "f_{21} \\\\\n",
    "f_{22} \\\\\n",
    "f_{23} \\\\\n",
    "f_{31} \\\\\n",
    "f_{32} \\\\\n",
    "f_{33}\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The 2D **Laplacian** is the **Kronecker sum** of our original matrix. The second derivative is given by:\n",
    "\n",
    "$$\n",
    "(D\\oplus D) \n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "f_{11} \\\\\n",
    "f_{12} \\\\\n",
    "f_{13} \\\\\n",
    "f_{21} \\\\\n",
    "f_{22} \\\\\n",
    "f_{23} \\\\\n",
    "f_{31} \\\\\n",
    "f_{32} \\\\\n",
    "f_{33}\n",
    "\\end{array}\n",
    "\\right], \\text{  where }\n",
    "D = \\frac{1}{\\Delta x^2}\n",
    "\\left[\n",
    "\\begin{array}{ccccc}\n",
    "-2 & 1 & 0 & 0 & 0 \\\\\n",
    "1 & -2 & 1 & 0 & 0 \\\\\n",
    "0 & 1 & -2 & 1 & 0 \\\\\n",
    "0 & 0 & 1 & -2 & 1 \\\\\n",
    "0 & 0 & 0 & 1 & -2 \n",
    "\\end{array}\n",
    "\\right],\n",
    "$$\n",
    "\n",
    "Note that even though $D$ is a $3 \\times 3$ matrix, $D \\oplus D$ is a $9 \\times 9$ matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import kron # kronecker PRODUCT, not kronecker SUM\n",
    "N = 5\n",
    "main_diag = -2 * np.ones(N)\n",
    "off_diag = np.ones(N-1)\n",
    "derivative_matrix = np.diag(main_diag) + np.diag(off_diag, k=-1) + np.diag(off_diag, k=1)\n",
    "derivative_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70",
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative_matrix_kronsum = kron(derivative_matrix, np.identity(N)) + kron(np.identity(N), derivative_matrix) # kronecker sum is defined in terms of the Kronecker product\n",
    "derivative_matrix_kronsum # is a 25 x 25 matrix!!!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71",
   "metadata": {},
   "source": [
    "Already for just a 5 x 5 starting matrix, the Kronecker sum is too big for any reasonable calculation! But we can use sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import sparse\n",
    "N = 100\n",
    "diag = np.ones([N])\n",
    "diags = np.array([diag, -2*diag, diag])\n",
    "D = sparse.spdiags(diags, np.array([-1,0,1]), N, N)\n",
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73",
   "metadata": {},
   "outputs": [],
   "source": [
    "sparse.kronsum(D,D)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74",
   "metadata": {},
   "source": [
    "### 4.8. Exercise: $LU$ decomposition 🌶️\n",
    "\n",
    "For matrices $P, L,$ and $U$ defined in [§4.5](#45-matrix-decomposition-decomposition), express the following matrix as a product of matrices $P, L, U$.\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    "9 & 3 & 3 \\\\\n",
    "3 & 2 & 2 \\\\\n",
    "3 & 4 & 2\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "import numpy as np\n",
    "from scipy.linalg import lu\n",
    "\n",
    "def solution_lu(): # don't change the function signature\n",
    "    # 1. TODO: define the matrix here\n",
    "\n",
    "    # 2. TODO: call the lu function here\n",
    "\n",
    "    # 3. TODO: return P, L, U matrices in this order here\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77",
   "metadata": {},
   "source": [
    "### 4.9. Exercise: the singular value decomposition (SVD) 🌶️\n",
    "\n",
    "The SVD is the factorization of a matrix $A$ as a product of matrices $U$, $S$, and $V^*$:\n",
    "\n",
    "$$ A = USV^*$$\n",
    "\n",
    "where $U$ is a unitary matrix, $S$ is a diagonal positive semi-definite matrix (the diagonal entries are known as the **singular values**), and $V^*$ is the conjugate transpose of unitary matrix $V$. The singular value decomposition is one of the most important matrix decomposition factorizations and has applications across many areas of the natural and social sciences. For a formal introduction to the SVD, we refer to any linear algebra book. Find the SVD of the matrix in [Exercise 4.8](#48-exercise-eigenvalue-problem)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "import numpy as np\n",
    "from scipy.linalg import svd\n",
    "\n",
    "def solution_svd(): # don't change the function signature\n",
    "    # 1. TODO: define the matrix here\n",
    "\n",
    "    # 2. TODO: call the svd method here\n",
    "\n",
    "    # 3. TODO: return U, S, V* matrices in this order here as a tuple\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80",
   "metadata": {},
   "source": [
    "## 5. Use case: interpolation\n",
    "\n",
    "Suppose now that we have the following data (this is of course a toy example -- the data are automatically generated here, but perhaps you collect $x$ and $y$ by means of some experiment):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace(0,10,10)\n",
    "y = x**2 * np.sin(x)\n",
    "plt.scatter(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82",
   "metadata": {},
   "source": [
    "Again -- we imagine that these are _collected_ data and not data which we have (artificially) generated by means of an analytical function. Now let's say that we want to know the values in between the above datapoints: we want to __interpolate__ between datapoints. We can do this using the `interp1d` library within `scipy.interpolate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "g_linear = interp1d(x, y, kind='linear') # Linear interpolation\n",
    "x_dense = np.linspace(0, 10, 100)\n",
    "y_dense = g_linear(x_dense)\n",
    "plt.plot(x_dense, y_dense)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84",
   "metadata": {},
   "source": [
    "We observe that this interpolation is not the best: it is very choppy and abrupt. Let's try a `cubic` interpolation now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85",
   "metadata": {},
   "outputs": [],
   "source": [
    "g_cubic = interp1d(x, y, kind='cubic') # Cubic interpolation\n",
    "x_dense = np.linspace(0, 10, 100)\n",
    "y_dense = g_cubic(x_dense)\n",
    "plt.plot(x_dense, y_dense)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86",
   "metadata": {},
   "source": [
    "### 5.1. Curve fitting\n",
    "A closely related problem to both optimization (see above) and interpolation is **curve fitting**. **Curve fitting**, while related, is however different from interpolation, in that interpolation involves obtaining the values of the curve between known datapoints, while curve fitting means knowing _a priori_ what kind of curve a (known!) set of datapoints should obey.\n",
    "\n",
    "In curve fitting, we are given a vector of locations $\\mathbf{x} =  \\begin{bmatrix} f(x_1) & f(x_2) & \\ldots & f(x_n) \\end{bmatrix}$, a vector of  observations $\\mathbf{y}$ and a function $f(x, \\mathbf{k})$ of a location and a set of parameters $\\mathbf{k}$  we look for the values of $\\mathbf{k}$ that minimize (usually the Euclidian) distance between the points of the function $f$ evaluated at the locations $\\mathbf{x}$: $f(\\mathbf{x}, \\mathbf{k}) = \\mathbf{\\hat{y}}  = \\begin{bmatrix} f(x_1) & f(x_2) & \\ldots & f(x_n) \\end{bmatrix}$ and the vector of observations $\\mathbf{y}$:\n",
    "\n",
    "$\\underset{\\mathbf{k}}{\\text{min}} \\, \\| f(\\mathbf{x}, \\mathbf{k}) - \\mathbf{y} \\|_2^2$\n",
    "\n",
    "This means that we look for a set of parameters $\\mathbf{k}$ so that the squared differences between the observed values and the values *simulated* by the function are minimized. \n",
    "Curve fitting is useful for example when you have a set of experimental data and you want to use them to determine parameters of a physical model you designed your experiment for.\n",
    "\n",
    "We consider a trivial example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_data = np.linspace(0,10,10)\n",
    "y_data = 3*x_data**2 + 2\n",
    "plt.scatter(x_data, y_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88",
   "metadata": {},
   "source": [
    "Let's say that we don't want to interpolate, but we really want to fit to a curve: here, we assume that the curve is quadratic.\n",
    "We want to fit the data to $y=ax^2+b$. The main gola here is to determine the values of $a$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def func(x,a,b):\n",
    "    return a*x**2 + b\n",
    "\n",
    "#in order to get parameters a and b:\n",
    "popt, pcov = curve_fit(func, x_data, y_data, p0=(1,1)) #p0: initial guess\n",
    "\n",
    "# popt: optimal parameters of a and b for the curve fit\n",
    "# pcov: covariance of curve fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90",
   "metadata": {},
   "outputs": [],
   "source": [
    "popt # gives a and b from the quadratic fit: a = 3, b = 2   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91",
   "metadata": {},
   "source": [
    "### 5.2. Simple harmonic motion\n",
    "\n",
    "The equation for spring motion is $y(t) = A \\cos (\\omega t + \\phi)$. we want to find the natural frequency of oscillation $\\omega$ for the spring. You collect the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_data = np.array([ 0.   ,  0.34482759,  0.68965517,  1.03448276,  1.37931034,\n",
    "        1.72413793,  2.06896552,  2.4137931 ,  2.75862069,  3.10344828,\n",
    "        3.44827586,  3.79310345,  4.13793103,  4.48275862,  4.82758621,\n",
    "        5.17241379,  5.51724138,  5.86206897,  6.20689655,  6.55172414,\n",
    "        6.89655172,  7.24137931,  7.5862069 ,  7.93103448,  8.27586207,\n",
    "        8.62068966,  8.96551724,  9.31034483,  9.65517241, 10.        ])\n",
    "        \n",
    "y_data = np.array([ 4.3303953 ,  1.61137995, -2.15418696, -3.90137249, -1.67259042,\n",
    "        2.16884383,  3.86635998,  1.85194506, -1.8489224 , -3.96560495,\n",
    "       -2.13385255,  1.59425817,  4.06145238,  1.89300594, -1.76870297,\n",
    "       -4.26791226, -2.46874133,  1.37019912,  4.24945607,  2.27038039,\n",
    "       -1.50299303, -3.46774049, -2.50845488,  1.20022052,  3.81633703,\n",
    "        2.91511556, -1.24569189, -3.72716214, -2.54549857,  0.87262548])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t_data, y_data, 'o--')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94",
   "metadata": {},
   "source": [
    "From elementary mechanics, we know that $\\omega=2\\pi f, f=1/T$, and $T \\approx 2$ seconds (by inspection). Thus, a good initial guess set of parameters is:\n",
    "- $\\omega=2 \\pi (1/2) = \\pi$\n",
    "- $A=4$ (by inspection)\n",
    "- $\\phi = 0$ (inspection: a cosine function with no phase shift)\n",
    "\n",
    "We have a total of three different parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95",
   "metadata": {},
   "outputs": [],
   "source": [
    "def spring(t, A, w, phi):\n",
    "    return A*np.cos(w*t+phi)\n",
    "\n",
    "popt, pcov = curve_fit(spring, t_data, y_data, p0=(4,np.pi,0)) # values for initial guess p0 must be in same order as arguments defined in the return statement for the function spring(). Here, A=4, w=pi, phi=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96",
   "metadata": {},
   "outputs": [],
   "source": [
    "popt # shows an array of [A, w, phi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_fit, w_fit, phi_fit = popt # unpacking the array into individual variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,10,100)\n",
    "y = spring(t, A_fit, w_fit,phi_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99",
   "metadata": {},
   "source": [
    "Now we can make a scatterplot of collected data and the curve using the fitted values for the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "100",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(t_data, y_data)\n",
    "plt.plot(t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "101",
   "metadata": {},
   "source": [
    "How do we get the errors? We know from statistics that the errors are generally contained in the __covariance__ matrix, whereby the __standard deviation__ appears in the diagonal elements of the covariance matrix. `pcov` gives us the covariance matrix in `scipy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "102",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcov # the covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "103",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.diag(pcov) # extracts the variance of the three parameters A, w, phi;\n",
    "              # which are the elements of the main diagonal of pcov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "104",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(np.diag(pcov)) # the standard deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105",
   "metadata": {},
   "source": [
    "### 5.3. Exercise: curve fitting 🌶️\n",
    "\n",
    "Plot the following data and to fit them to a Gaussian curve: $$y=\\frac{1}{\\sqrt{2\\pi} \\sigma}\\exp{\\left[- \\frac{(x-\\mu)^2}{2\\sigma^2} \\right]}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "106",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "xdata = [ -10.0, -9.0, -8.0, -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]\n",
    "ydata = [1.2, 4.2, 6.7, 8.3, 10.6, 11.7, 13.5, 14.5, 15.7, 16.1, 16.6, 16.0, 15.4, 14.4, 14.2, 12.7, 10.3, 8.6, 6.1, 3.9, 2.1]\n",
    "xdata = np.array(xdata)\n",
    "ydata = np.array(ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. TODO: plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "108",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "109",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def solution_gaussian(): # do not change the function signature\n",
    "    # 2. TODO: define function to fit: the Gaussian.\n",
    "\n",
    "    def gaussian_func(x, a, b):\n",
    "        \"\"\"\n",
    "        Use this function signature, and these parameters A and B.\n",
    "        a: coefficient of the Gaussian\n",
    "        b: coefficient of the x^2 term\n",
    "\n",
    "        Will probably return an error if you change this.\n",
    "\n",
    "        Returns: the Gaussian function using the parameters a and b\n",
    "        \"\"\"\n",
    "    \n",
    "    # 3. TODO: call the curve_fit function here and return the parameters and covariance matrix AS A TUPLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "110",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. TODO: plot the fitted Gaussian here together with the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "111",
   "metadata": {},
   "source": [
    "## 6. Fun with functions: the Legendre polynomials and the Bessel functions\n",
    "\n",
    "The **Legendre polynomials** $P_l(x)$ are the solutions $y=P_l(x)$ to the differential equation $(1-x^2)y'' - 2xy' + l(l+1)y = 0$\n",
    "- represents the angular component of the spherical Schrödinger equation which permits non-infinite solutions: we remember from the solution of the hydrogen atom\n",
    "- $l$ is an integer because these are non-infinite solutions, so they are normalizable wavefunctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "112",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import legendre\n",
    "x = np.linspace(0,1,100)\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.title(label='the Legendre polynomials')\n",
    "plt.plot(x,legendre(0)(x), label='x=0')\n",
    "plt.plot(x,legendre(1)(x), label='x=1')\n",
    "plt.plot(x,legendre(2)(x), label='x=2')\n",
    "plt.plot(x,legendre(3)(x), label='x=3')\n",
    "plt.plot(x,legendre(4)(x), label='x=4')\n",
    "plt.plot(x,legendre(5)(x), label='x=5')\n",
    "plt.plot(x,legendre(6)(x), label='x=6') # the sixth legendre polynomial\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "113",
   "metadata": {},
   "source": [
    "The **Bessel functions** $J_a(x)$ are the solutions $y=J_a(x)$ to $x^2y'' + xy' + (x^2-\\alpha^2)y = 0$\n",
    "- represents **Laplace's equation** in polar coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "114",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import jv\n",
    "x = np.linspace(0,10,100)\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.title(label='the Bessel functions')\n",
    "plt.plot(x,jv(0,x), label='0')\n",
    "plt.plot(x,jv(1,x), label='1')\n",
    "plt.plot(x,jv(2,x), label='2')\n",
    "plt.plot(x,jv(3,x), label='3') #here, 3 is the value of alpha\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "115",
   "metadata": {},
   "source": [
    "## Mathematical analysis\n",
    "\n",
    "## 7. Differentiation\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <h4><b>Attention!</b></h4> `derivative` has been deprecated from `scipy.misc.derivative` in SciPy 1.10.0 and will be completely removed in SciPy 1.12.0. For alternatives, consult `findiff`: https://github.com/maroba/findiff (finite differences) or `numdifftools`: https:/github.com/pbrod/numdifftools. Nonetheless, for completeness, we include this short section on differentiation.\n",
    "</div>\n",
    "\n",
    "- `scipy` and `numpy` treat differentiation somewhat differently: `numpy` uses actual numerical data to evaluate a derivative using an array of values; `scipy` uses an actual function\n",
    "- for this reason, `scipy` is somewhat better for differentiation if the actual analytic form of a function is known"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "116",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.misc import derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "117",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x**2 * np.sin(2*x)*np.exp(-x)\n",
    "x = np.linspace(0,1,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "118",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_prime = derivative(f,x,dx=1e-6)\n",
    "f_doubleprime = derivative(f,x,dx=1e-6, n=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,7))\n",
    "plt.title(label='$f(x)=x^2 e^{-x} \\sin(2x)$', size=20)\n",
    "plt.plot(x,f(x), label='f')\n",
    "plt.plot(x,f_prime, label='f\\'')\n",
    "plt.plot(x, f_doubleprime, label='f\\'\\'')\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "120",
   "metadata": {},
   "source": [
    "## 8. Integration\n",
    "\n",
    "Integration is also slightly different than in `numpy`, in which integrals are evaulated using more of a Riemann sum method. In `scipy`, these sums are done more behind the scenes.\n",
    "- we consider the one-dimensional integral $$\\int_0^1 x^2 \\sin(2x)e^{-x} dx$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "121",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad # the main integration package in scipy for 1D integrals\n",
    "integrand = lambda x: x**2 * np.sin(x) * np.exp(-x) # integrals done in scipy must be defined as a function: \n",
    "                                                    # here, the integrand is defined as a lambda function\n",
    "integral, integral_error = quad(integrand, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "122",
   "metadata": {},
   "outputs": [],
   "source": [
    "integral, integral_error # the integral and the error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "123",
   "metadata": {},
   "source": [
    "- we next consider the double integral $$\\int_0^1 \\int_{-x}^{x^2} \\sin(x+y^2) dy dx$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "124",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import dblquad \n",
    "integrand = lambda y, x: np.sin(x+y**2) # order of lambdas is important! \n",
    "                                        # the variable that correspond to the innermost integral must come \n",
    "                                        # first in the definition of lambda variables\n",
    "lower_y = lambda x: -x\n",
    "upper_y = lambda x: x**2\n",
    "integral, integral_error = dblquad(integrand, 0, 1, lower_y, upper_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "125",
   "metadata": {},
   "outputs": [],
   "source": [
    "integral, integral_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "126",
   "metadata": {},
   "source": [
    "### 8.1. Exercise: double integration 🌶️\n",
    "\n",
    "Use `SciPy` to evaluate the integral $$\\int_0^2 \\int_x^{6-x^2} (3x+2y) dy dx$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "127",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "128",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "from scipy.integrate import dblquad \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def solution_dblquad(): # do not change the function name\n",
    "    # 1. TODO: define the integrand here, minding the order of the lambda variables\n",
    "    # pay attention to the order of the defined arguments. Make sure to define 'y' to appear first\n",
    "    # as it corresponds to the inner integral\n",
    "\n",
    "    def f(y, x):\n",
    "        pass\n",
    "    \n",
    "    # 2. TODO: define the upper and lower limits\n",
    "\n",
    "    # 3. TODO: call the dblquad function here\n",
    "\n",
    "    # 4. TODO: return the integral and the error as a tuple\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "129",
   "metadata": {},
   "source": [
    "## 9. Differential equations\n",
    "\n",
    "### 9.1. First-order ODEs\n",
    "\n",
    "We consider the initial-value problem of air friction while falling: $$v' - \\alpha v^2 + \\beta = 0, v(0)=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "130",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint # there's also solve_ivp: \n",
    "                                   # which package to use depends on the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "131",
   "metadata": {},
   "source": [
    "Need to provide to the interpreter all the information about the differential equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "132",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dvdt(v,t):\n",
    "    return 3*v**2 - 5 # arbitrarily choosing alpha=3, belta=5: notice how we isolate v'!!!!!!\n",
    "v0 = 0 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "133",
   "metadata": {},
   "source": [
    "Solve differential equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "134",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,1,100)\n",
    "sol = odeint(dvdt, v0,t)\n",
    "sol\n",
    "sol.T\n",
    "sol.T[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "135",
   "metadata": {},
   "source": [
    "We can then plot this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "136",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,7))\n",
    "plt.plot(t,sol.T[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "137",
   "metadata": {},
   "source": [
    "### 9.2. Exercise: first-order ODEs\n",
    "\n",
    "Use `SciPy` to plot and solve the following differential first-order ODE: $$\\frac{dy}{dt}+ky = 0$$ Use the initial condition $y_0=5$ and solve the equation for $k=0.2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "138",
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext tutorial.tests.testsuite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "139",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%ipytest\n",
    "\n",
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "\n",
    "def solution_odeint(): # do not change the function name\n",
    "    # 1. TODO: define the function dydt here\n",
    "\n",
    "    # 2. TODO: define the initial condition here\n",
    "\n",
    "    # 3. TODO: define the time array here using np.linspace(0,10,101)\n",
    "\n",
    "    # 4. TODO: call odeint here\n",
    "    \n",
    "    # 5. TODO: return the solution here, don't forget to transpose the solution and take the [0]th element ;)!\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "140",
   "metadata": {},
   "source": [
    "### 9.3. Coupled first-order ODEs\n",
    "\n",
    "Consider the following system of differential equations:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "y_1' = y1 + y_2^2 + 3x & y_1(0) =0 \\\\\n",
    "y_2' = 3y1 + y_2^3 - \\cos(x) & y_2(0)=0\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "In order to solve coupled first-order ODEs in python: we need to define a vector $$S=(y_1,y_2).$$\n",
    "\n",
    "Letting $S = (y_1,y_2)$, we need to write a function that returns $dS/dx = (dy_1/dx, dy_2/dx)$. The function `dSdx` takes in `S` and `x`, where `S` is just a vector, `x` is more of just a placeholder variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "141",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dSdx(S,x):\n",
    "    y1, y2 = S\n",
    "    return [y1+y2**2+3*x,\n",
    "            3*y1+y2**2-np.cos(x)]\n",
    "y1_0 = 0\n",
    "y2_0 = 0\n",
    "S_0 = (y1_0, y2_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "142",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,1,100)\n",
    "sol = odeint(dSdx,S_0,x)\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "143",
   "metadata": {},
   "outputs": [],
   "source": [
    "y1 = sol.T[0]\n",
    "y2 = sol.T[1]\n",
    "plt.figure(figsize=(12,7))\n",
    "plt.plot(x,y1, label='y1')\n",
    "plt.plot(x,y2, label='y2')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144",
   "metadata": {},
   "source": [
    "Here we can see that the solution depends on the set of initial conditions, as expected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "145",
   "metadata": {},
   "source": [
    "### 9.4. Second-order ODEs\n",
    "\n",
    "We next consider second-order ordinary differential equations. We want to solve a second-order ODE by splitting it into two first-order ODEs. From mechanics, we know that the **pendulum equation** reads: $$\\theta '' - \\sin(\\theta)=0.$$\n",
    "\n",
    "Scipy can only solve coupled first-order ODEs, but **any second-order ODE can be turned into two coupled first-order ODEs**. The same thing goes for higher-order ODEs.\n",
    "\n",
    "Define $\\omega = d\\theta / dt$ so that we obtain the following coupled first-order ODEs: $$d\\omega / dt = \\sin(\\theta) \\\\ d\\theta / dt = \\omega.$$\n",
    "Let $S=(\\theta, \\omega)$ so that $dS/dt = (d\\theta / dt, d\\omega /dt)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "146",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dSdt(S, t):\n",
    "    theta, omega = S\n",
    "    return [omega, np.sin(theta)]\n",
    "\n",
    "theta0 = np.pi/4\n",
    "omega0 = 0\n",
    "S0 = (theta0, omega0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "147",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,20,100)\n",
    "sols = odeint(dSdt, S0, t)\n",
    "thetas, omegas = sols.T\n",
    "plt.figure(figsize=(12,7))\n",
    "plt.plot(t, thetas)\n",
    "plt.plot(t, omegas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "148",
   "metadata": {},
   "source": [
    "## 10. Fourier transforms\n",
    "\n",
    "The discrete Fourier transform is defined as \n",
    "\n",
    "$$\n",
    "y[k] = \\sum_{n=0}^{N-1} e^{-2\\pi i n (k/N)}x[n],\n",
    "$$\n",
    "\n",
    "where\n",
    "- $k/N$ represents a specific frequency (dimensionless);\n",
    "- $y[k]$ can be converted to a frequency in Hz if the spacing in $x$ is known.\n",
    "- $x$ is the time domain, $y$ is the frequency domain\n",
    "\n",
    "Consider some periodic function $x(t) = \\sin(2\\pi t) + \\sin(4 \\pi t) + 0.1 \\text{ rand}(t)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "149",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,10*np.pi, 100)\n",
    "x = np.sin(2*np.pi*t) + np.sin(4*np.pi*t) + 0.1*np.random.randn(len(t))\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.title(label=f'$x(t) = \\sin(2\\pi t) + \\sin(4 \\pi t) + 0.1 rand(t)$', size=20)\n",
    "plt.plot(t,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "150",
   "metadata": {},
   "source": [
    "We can then compute the **fast Fourier transform** of this function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "151",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.fft import fft, fftfreq\n",
    "N = len(x)\n",
    "y = fft(x) # the Fourier transform of x\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.plot(np.abs(y), label = 'abs(y)')\n",
    "plt.plot(y, label = 'y')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "152",
   "metadata": {},
   "source": [
    "Note that it's symmetric, which is true for any real-valued time series: its Fourier transform is symmetric, so what we actually have plotted above are positive frequencies up to half the time scale, and then goes over to the negative frequencies for the second half of the time scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "153",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = fft(x)[:N//2]\n",
    "f = fftfreq(N, np.diff(t)[0])[:N//2] # plots the power\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.plot(f, np.abs(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "154",
   "metadata": {},
   "source": [
    "### 11. Statistics\n",
    "\n",
    "#### 11.1. The $\\beta$ distribution\n",
    "\n",
    "Consider the $\\beta$ distribution, given by:\n",
    "\n",
    "$$\n",
    "f(x; a,b) = \\frac{\\Gamma(a+b)x^{a-1}(1-x)^{b-1}}{\\Gamma(a) \\Gamma(b)}, 0 \\leq x \\leq 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "155",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import beta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "156",
   "metadata": {},
   "source": [
    "Basic statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "157",
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = 2.5, 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "158",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean, var, skewness, kurtosis = beta.stats(a, b, moments='mvsk') #mvsk: mean variance skewness kurtosis\n",
    "mean, var, skewness, kurtosis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "159",
   "metadata": {},
   "source": [
    "Probability distribution plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "160",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta.ppf? #percent point function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "161",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(beta.ppf(0.01, a, b), beta.ppf(0.99, a, b), 100)\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.plot(x, beta.pdf(x,a,b))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "162",
   "metadata": {},
   "source": [
    "Generating random variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "163",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = beta.rvs(a,b,size=10)\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164",
   "metadata": {},
   "source": [
    "#### 11.2. the Gaussian distribution\n",
    "\n",
    "Recall the Gaussian distribution:\n",
    "\n",
    "$$\n",
    "f(x; \\mu, \\sigma) = \\frac{1}{\\sqrt{2\\pi \\sigma^2}} \\exp{-\\frac{(x-\\mu)^2}{\\sigma^2}}, \\text{ with } -\\infty \\lt x \\lt \\infty\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "165",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "166",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 1\n",
    "sigma = 2\n",
    "mean, var, skew, kurt = norm.stats(loc=mu, scale=sigma, moments='mvsk')\n",
    "x = np.linspace(norm.ppf(0.01, mu, sigma), norm.ppf(0.99, mu, sigma), 100)\n",
    "plt.figure(figsize=(15,9))\n",
    "plt.plot(x, norm.pdf(x,mu,sigma))\n",
    "plt.axvline(mean, ls='--', color='purple')\n",
    "plt.axvline(var, ls='--', color='red')\n",
    "plt.axvline(skew, ls='--', color='black')\n",
    "plt.axvline(kurt, ls='--', color='yellow')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "167",
   "metadata": {},
   "source": [
    "#### 11.3. The multinomial distribution\n",
    "\n",
    "Analogous to the binomial distribution, the multinomial distribution is:\n",
    "\n",
    "$$\n",
    "f(x_1,x_2,\\dots,x_k;p_1,p_2,\\dots,p_k;n) = \\frac{n!}{x_1!x_2!\\dots x_k!} p_1^{x_1}p_2^{x_2} \\dots p_k^{x_k}  \n",
    "$$\n",
    "\n",
    "the variables $x_i$ represent the variables as numbers; $p_i$ represent the probabilities; $n$ is the number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "168",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import multinomial\n",
    "\n",
    "#rolling a die\n",
    "p = np.ones(6)/6\n",
    "multinomial.pmf([6,0,0,0,0,0], n=6, p=p) # the probability of rolling six ones on a die, rolling six times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "169",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate multinomial random numbers\n",
    "\n",
    "multinomial.rvs(n=100, p=p, size=5) # 100 rolls, 6 different options, each with probability of 1/6, and performing this 5 times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170",
   "metadata": {},
   "source": [
    "### 12. Problems (freehand)\n",
    "\n",
    "There are no comparison solutions for the following problems.\n",
    "\n",
    "#### 12.1. Mechanical work 🌶️🌶️\n",
    "\n",
    "The energy required to get from point $\\vec{r}_1$ to point $\\vec{r}_2$ for a plane is given by\n",
    "\n",
    "$$\n",
    "E = \\alpha \\int_C \\left| \\frac{d\\vec{r}}{dt} \\right| dt - \\int_C \\vec{F} \\cdot \\frac{d\\vec{r}}{dt}dt \n",
    "$$\n",
    "\n",
    "Suppose that $\\alpha = 5$ and our start- and endpoints are $\\vec{r}_1 = (0,0)$ and $\\vec{r}_2 = (0,10)$. On this particular day, the wind produces a forcefield $\\vec{F} = (0, -2/(x+1)^2)$. Find the optimal value of $A$ in $\\vec{r}(t) = A \\sin(\\pi t / 10)\\hat{x} + t \\hat{y}$ that minimises the work."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "171",
   "metadata": {},
   "source": [
    "#### 12.2. Newton's law of cooling 🌶️🌶️\n",
    "\n",
    "Newton's law of cooling is given by \n",
    "\n",
    "$$\n",
    "\\frac{dT}{dt}=-k(T-T_s(t))\n",
    "$$\n",
    "\n",
    "where $T$ is the temperature of an object in an ambient, time-varying temperature field $T(s)$. Suppose $T$ represents the temperature of a shallow pool of water and $T_s(t)$ represents the temperature of outside. Find $T(t)$ given that you collected measurements of the outside:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "172",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_m = np.array([ 0.,  1.04347826,  2.08695652,  3.13043478,  4.17391304,\n",
    "        5.2173913 ,  6.26086957,  7.30434783,  8.34782609,  9.39130435,\n",
    "       10.43478261, 11.47826087, 12.52173913, 13.56521739, 14.60869565,\n",
    "       15.65217391, 16.69565217, 17.73913043, 18.7826087 , 19.82608696,\n",
    "       20.86956522, 21.91304348, 22.95652174, 24.        ])\n",
    "\n",
    "temp_m = np.array([283.2322975, 284.6945461, 286.2259041, 287.8603625, 289.6440635,\n",
    "       291.6187583, 293.7939994, 296.1148895, 298.4395788, 300.5430675,\n",
    "       302.1566609, 303.0363609, 303.0363609, 302.1566609, 300.5430675,\n",
    "       298.4395788, 296.1148895, 293.7939994, 291.6187583, 289.6440635,\n",
    "       287.8603625, 286.2259041, 284.6945461, 283.2322975])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "173",
   "metadata": {},
   "source": [
    "#### 12.3. The Lennard-Jones potential 🌶️🌶️\n",
    "\n",
    "The Lennard-Jones potential $V_{LJ}$ is a simplified model for an intermolecular pair (ie, two-body!) potential between electrically neural atoms or molecules, and has been used extensively to model van der Waals interactions: $$V_{LJ}=4\\epsilon \\left[ \\left(\\frac{\\sigma}{R}\\right)^{12} \\left(\\frac{\\sigma}{R}\\right)^6 \\right],$$\n",
    "\n",
    "for internuclear distance $R$, well depth $\\epsilon$, and particle size $\\sigma$.\n",
    "\n",
    "Consider the folowing data measured for a helium dimer. The interaction at several different internuclear separations is given. Fit the data to a LJ potential by finding the values of $\\epsilon$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "174",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Internuclear separation in angstroms\n",
    "distances = [2.875, 3.0, 3.125, 3.25, 3.375, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0]\n",
    "# Energy in Wavenumbers\n",
    "energies = [0.35334378061169025, -2.7260131253801405, -4.102738968283382, -4.557042640311599, -4.537519193684069, -4.296388508321034, -3.6304745046204117, -3.0205368595885536, -2.1929538006724814, -1.7245616790238782, -1.2500789753171557]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "175",
   "metadata": {},
   "source": [
    "#### 12.4. Random number generation from your own distribution 🌶️🌶️🌶️\n",
    "\n",
    "Using SciPy, it's also possible to generate numbers from a probability distribution which is not inherently built in to SciPy. Use the SciPy documentation online, and the fact that you want to generate a _continuous_ (as opposed to discrete) random variable, to create your own Python class and generate random numbers from the following distribution:\n",
    "\n",
    "$$\n",
    "f(x; a_1, a_2, b_1, b_2) = \\frac{1}{2(a_1b_1 + a_2b_2)} \\left( b_1 \\exp{-\\left(\\sqrt{\\frac{x}{a_1}}\\right)} + b_2 \\exp{-\\left(\\sqrt{\\frac{x}{a_2}}\\right)} \\right), \\text{ with } 0 \\leq x \\leq \\infty\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "176",
   "metadata": {},
   "source": [
    "#### 12.5. The finite square well 🌶️🌶️🌶️🌶️\n",
    "\n",
    "One of the fundamental results (and pains) of first-semester quantum mechanics is the solution of the quantum finite one-dimensional square well. While we will not go through the entire derivation here, the solution of this problem essentially boils down to the curve $$\\tan(x) = x.$$\n",
    "\n",
    "Find the solutions to $\\tan(x) = x$ on the interval $[-2\\pi, 2\\pi].$ Be careful about how the asymptotes are being numerically handled."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "177",
   "metadata": {},
   "source": [
    "#### 12.6. the quintic polynomial 🌶️🌶️\n",
    "\n",
    "Polynomials of the form $$ a_5x^5 + a_4x^4 + a_3x^3 + a_2x^2 +a_1x + a_0 = 0 $$ are known as **quintic polynomials**, and their solution was a major problem in the algebra of the 16th century. While explicit, analytical solutions exist for quadratic, cubic, and quartic polynomials, quintic and higher-degree polynomials have no closed-form solution and must therefore must be solved by other means (see: *Abel-Ruffini theorem*). \n",
    "\n",
    "The solution of quintic equations arises in celestial mechanics: for instance, in finding the Lagrange points for a third massive body, wherefore the masses of the two other bodies are both non-negligible (for example, a Sun-Earth-satellite system). The location of $L_2$ and $L_1$ is given by\n",
    "\n",
    "* $a_5=\\pm (M_S+M_E)$\n",
    "* $a_4=+3R(M_S+M_E)$\n",
    "* $a_3=\\pm 3R^2(M_S+M_E)$\n",
    "* $a_2=+R^3(M_E \\mp M_E)$\n",
    "* $a_1=\\pm M_E 2R^4$\n",
    "* $a_0 = \\mp M_E R^5$\n",
    "\n",
    "where $\\pm$ corresponds to $L_2$ and $L_1$ respectively. Here, $M_E$ and $M_S$ are the mass of the Earth and Sun, respectively and $R$ is the Sun-Earth distance. Then, in the quintic polynomial above, $x$ represents the distance of the satellite to Earth. Find $L_2$ and $L_1$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}