maths-maths-coding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bf4e32e9",
   "metadata": {},
   "source": [
    "# Maths with Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e1f0ee",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set seed for reproducibility\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03e19332",
   "metadata": {},
   "source": [
    "In this chapter, you'll learn about doing mathematics with code, including solving equations both in the abstract and numerically.\n",
    "\n",
    "This chapter uses the **numpy**, **scipy**, and **sympy** packages. If you're running this code, you may need to install these packages using, for example, `pip install packagename` on your computer's command line. (If you're not sure what a command line or terminal is, take a quick look at the basics of coding chapter.)\n",
    "\n",
    "## Symbolic mathematics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbfd59c5",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from myst_nb import glue\n",
    "import sympy\n",
    "a = 8\n",
    "glue('sqrt', 2*np.sqrt(a))\n",
    "glue('symsqrt', sympy.sqrt(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1854dc54",
   "metadata": {},
   "source": [
    "When using computers to do mathematics, we're most often performing numerical computations such as $\\sqrt{8} = ${glue:}`sqrt`. Although we have the answer, it's only useful for the one special case. Symbolic mathematics allows us to use coding to solve equations in the general case, which can often be more illuminating. As an example, if we evaluate this in symbolic mathematics we get $\\sqrt{8} = ${glue:}`symsqrt`.\n",
    "\n",
    "The Python package for symbolic mathemtics is [**sympy**](https://www.sympy.org/en/index.html), which provides some features of a computer algebra system.\n",
    "\n",
    "To define *symbolic* variables, we use sympy's symbols function. For ease, we'll import the entire sympy library into the namespace by using `from sympy import *`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b922febe",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "x, t, α, β = symbols(r'x t \\alpha \\beta')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27a649be",
   "metadata": {},
   "source": [
    "```{note}\n",
    "The leading 'r' in some strings tells Python to treat the string literally so that backslashes are not treated as instructions--otherwise, combinations like `\\n` would begin a newline.\n",
    "```\n",
    "\n",
    "Having created these symbolic variables, we can refer to and see them just like normal variables--though they're not very interesting *because* they are just symbols (for now):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c07e7673",
   "metadata": {},
   "outputs": [],
   "source": [
    "α"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9b6f1c3",
   "metadata": {},
   "source": [
    "Things get much more interesting when we start to do maths on them. Let's see some integration, for example, say we want to evaluate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "339c917a",
   "metadata": {},
   "outputs": [],
   "source": [
    "Integral(log(x), x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bcf8450",
   "metadata": {},
   "source": [
    "(note that the symbols are printed as latex equations) we simply call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b04080ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(log(x), x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "534bdd5e",
   "metadata": {},
   "source": [
    "We can differentiate too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3ee815",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff(sin(x)*exp(x), x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "505c7358",
   "metadata": {},
   "source": [
    "and even take limits!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca9ff12",
   "metadata": {},
   "outputs": [],
   "source": [
    "limit(sin(x)/x, x, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6642aab",
   "metadata": {},
   "source": [
    "It is also possible to solve equations using **sympy**. The solve function tries to find the roots of $f(x)$ and has syntax `solve(f(x)=0, x)`. Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f05021ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(x*5 - 2, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c65aa51",
   "metadata": {},
   "source": [
    "There are also solvers for differential equations (`dsolve`), continued fractions, simplifications, and more.\n",
    "\n",
    "Another really important thing to know about symbolic mathematics is that you can 'cash in' at any time by substituting in an actual value. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc9f46ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = 1 - 2*sin(x)**2\n",
    "expr.subs(x, np.pi/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d31182",
   "metadata": {},
   "source": [
    "But you don't have to substitute in a real value; you can just as well substitute in a different symbolic variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da539ece",
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = 1 - 2*sin(x)**2\n",
    "simplify(expr.subs(x, t/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6cde99",
   "metadata": {},
   "source": [
    "I snuck in a simplify here too!\n",
    "\n",
    "### Symbolic mathematics for economics\n",
    "\n",
    "The library does a lot, so let's focus on a few features that are likely to be useful for economics in particular.\n",
    "\n",
    "#### Series expansion\n",
    "\n",
    "The first is performing **Taylor series expansions**. These come up all the time in macroeconomic modelling, where models are frequently log-linearised. Let's see an example of a couple of expansions together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0fed351",
   "metadata": {},
   "outputs": [],
   "source": [
    "expr = log(sin(α))\n",
    "\n",
    "expr.series(α, 0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5de5206",
   "metadata": {},
   "source": [
    "This is a 3rd order expansion around $\\alpha=0$.\n",
    "\n",
    "#### Symbolic linear algebra\n",
    "\n",
    "The support for **matrices** can also come in handy for economic applications. Here's a matrix,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3516fca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = Matrix([[1, 0, x], [α, -t, 3], [4, β, 2]])\n",
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6cf398",
   "metadata": {},
   "source": [
    "and its determinant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e60f0086",
   "metadata": {},
   "outputs": [],
   "source": [
    "M.det()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c243ac73",
   "metadata": {},
   "source": [
    "I can hardly go to a talk in economics that involves matrices that doesn't see those matrices get diagonalised: there's a function for that too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2aeac1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "P, D = Matrix([[1, 0], [α, -t]]).diagonalize()\n",
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1bbe6e0",
   "metadata": {},
   "source": [
    "#### Lagrangians\n",
    "\n",
    "Function optimisation using Lagrangians is about as prevalent in economics as any bit of maths: let's see how it's done symbolically.\n",
    "\n",
    "We're going to find the minimum over x, y of the function $f(x,y)$, subject to $g(x,y)=0$, where $f(x,y) = 4xy - 2x^2 + y^2$ and $g(x,y) = 3x+y-5$.\n",
    "\n",
    "First we need to specify the problem, and the Lagrangian for it, in code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2fcd16b",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, λ = symbols(r'x y \\lambda', real=True)\n",
    "f = 4*x*y - 2*x**2 + y**2\n",
    "g = 3*x+y-5\n",
    "\n",
    "ℒ = f - λ*g\n",
    "ℒ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06d89459",
   "metadata": {},
   "source": [
    "The Karush-Kuhn-Tucker (KKT) conditions tell us whether any solutions we find will be optimal. Simply, the constaint is that a solution vector is a saddle point of the Lagrangian, $\\nabla \\mathcal{L} = 0$. Let's solve this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09b4df07",
   "metadata": {},
   "outputs": [],
   "source": [
    "gradL = [diff(ℒ, c) for c in [x, y]]\n",
    "KKT_eqns = gradL + [g]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f85ed8cf",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "KKT_eqns = gradL + [g]\n",
    "glue('kkt_0', KKT_eqns[0])\n",
    "glue('kkt_1', KKT_eqns[1])\n",
    "glue('kkt_2', KKT_eqns[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e2c6b67",
   "metadata": {},
   "source": [
    "This gives 3 equations from the KKT conditions:  {glue:}`kkt_0`,   {glue:}`kkt_1`, and  {glue:}`kkt_2`. (The symbolic manipulation is now over: we solved for the conditions in terms of algebra--now we're looking for real values.) Now we look for the values of $x, y$ that minimise $f$ given that $g=0$ by solving these equations over $x$, $y$, and $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0accf2db",
   "metadata": {},
   "outputs": [],
   "source": [
    "stationary_pts = solve(KKT_eqns, [x, y, λ], dict=True)\n",
    "stationary_pts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d23099aa",
   "metadata": {},
   "source": [
    "Now, we can substitute these in to find the (first--and in this case only) point that minimises our function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aacfdaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "stationary_pts[0][x], stationary_pts[0][y], f.subs(stationary_pts[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32c4b784",
   "metadata": {},
   "source": [
    "#### Exporting to latex\n",
    "\n",
    "To turn any equation, for example `diff(sin(x)*exp(x), x)`, into latex and export it to a file that can be included in a paper, use\n",
    "\n",
    "```python\n",
    "eqn_to_export = latex(diff(sin(x)*exp(x), x), mode='equation')\n",
    "open('latex_equation.tex', 'w').write(eqn_to_export)\n",
    "```\n",
    "\n",
    "which creates a file called 'latex_equation.tex' that has a single line in it: '\\begin{equation}\\int \\log{\\left(x \\right)}\\, dx\\end{equation}'. There are a range of options for exporting to latex, `mode='equation*'` produces an unnumbered equation, 'inline' produces an inline equation, and so on. To include these in your latex paper, use '\\input{latex_equation.tex}'.\n",
    "\n",
    "### Why coding symbolic mathematics is useful\n",
    "\n",
    "1. Accuracy--using a computer to solve the equations means you're less likely to make a mistake. At the very least, it's a useful check on your by-hand working.\n",
    "\n",
    "2. Consistency--by making your code export the equations you're solving to your write-up, you can ensure that the equations are consistent across both *and* you only have to type them once.\n",
    "\n",
    "## Numerical Mathematics\n",
    "\n",
    "For much of the time, you'll be dealing with numbers rather than symbols. The workhorses of numerical mathematics are the two packages **numpy** and **scipy**. Both have excellent documentation, where you can find out more. In this section, we'll look at how to use them in some standard mathematical operations that arise in economics.\n",
    "\n",
    "The most basic object is an array, which can be defined as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a61ecb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "a = np.array([0, 1, 2, 3], dtype='int64')\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1243113b",
   "metadata": {},
   "source": [
    "Arrays are very memory efficient and fast objects that you should use in preference to lists for any heavy duty numerical operation.\n",
    "\n",
    "To demonstrate this, let's do a time race between lists and arrays for squaring all elements of an array:\n",
    "\n",
    "Lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5246ee0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_list = range(1000)\n",
    "%timeit [i**2 for i in a_list]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f8d01e0",
   "metadata": {},
   "source": [
    "Arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c40b6403",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.arange(1000)\n",
    "%timeit a**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26ce79d4",
   "metadata": {},
   "source": [
    "Using arrays was *two orders of magnitude** faster! Okay, so we should use arrays for numerical works. How do we make them? You can specify an array explicitly as we did above to create a vector. This manual approach works for other dimensions too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde5edcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])\n",
    "mat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96be90e2",
   "metadata": {},
   "source": [
    "To find out about matrix properties, we use `.shape`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00f62dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "mat.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d841b4bd",
   "metadata": {},
   "source": [
    "We already saw how `np.arange(start, stop, step)` produces a vector; `np.linspace(start, stop, number)` produces a vector of length `number` by equally dividing the space between `start` and `stop`.\n",
    "\n",
    "Three really useful arrays are `np.ones(shape)`, for example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ddb5e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ones((3, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08bab38a",
   "metadata": {},
   "source": [
    "`np.diag` for diagnoal arrays,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "373af9c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.diag(np.array([1, 2, 3, 4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a08e557a",
   "metadata": {},
   "source": [
    "and `np.zeros` for empty arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f67a75a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.zeros((2, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f01e3d7",
   "metadata": {},
   "source": [
    "Random numbers are supplied by `np.random.rand()` for a uniform distribution in [0, 1], and `np.random.randn()` for numbers drawn from a standard normal distribution.\n",
    "\n",
    "You can, of course, specify a function to create an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4f4e5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.fromfunction(lambda i, j: i**2+j**2, (4, 5))\n",
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6914be94",
   "metadata": {},
   "source": [
    "To access values in an array, you can use all of the by-position slicing methods that you've seen already in data analysis and with lists. The figure gives an example of some common slicing operations:\n",
    "\n",
    "![Examples of array slices](https://scipy-lectures.org/_images/numpy_indexing.png)\n",
    "\n",
    "Arrays can also be sliced and diced based on boolean indexing, just like a dataframe.\n",
    "\n",
    "For example, using the array defined above, we can create a boolean array of true and false values from a condition such as `c > 6` and use that to only access some elements of an array (it doesn't have to be the same array, though it usually is):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c729b2b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "c[c > 6]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3498beb",
   "metadata": {},
   "source": [
    "As with dataframes, arrays can be combined. The main command to remember is `np.concatenate`, which has an `axis` keyword option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f35351",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.eye(3)\n",
    "np.concatenate([x, x], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "280a5a19",
   "metadata": {},
   "source": [
    "Splitting is performed with `np.split(array, splits, axis=)`, for example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad1f449f",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.split(x, [3], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84b1066",
   "metadata": {},
   "source": [
    "Aggregation operations are very similar to those found in dataframes: `x.sum(i)` to sum across the $i$th dimension of the array; similarly for standard deviation, and so on.\n",
    "\n",
    "As with dataframes, you can (and often should) specify the datatype of an array when you create it by passing a `dtype=` keyword, eg `c = np.array([1, 2, 3], dtype=float)`. To find out the data type of an array that already exists, use `c.dtype`.\n",
    "\n",
    "Finally, numpy does a lot of smart broadcasting of arrays. Broadcasting is what means that summing two arrays gives you a third array that has elements that are each the sum of the relevant elements in the two original arrays. Put another way, it's what causes `x + y = z` (for arrays x and y with the same shape) to result in an array z for which $z_{ij} = x_{ij} + y_{ij}$.\n",
    "\n",
    "Summing two arrays of the same shape is a pretty obvious example, but it also applies to cases that are *not* completely matched. For example, multiplication by a scalar is broadcast across all elements of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7d525c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.ones(shape=(3, 3))\n",
    "x*3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4646540",
   "metadata": {},
   "source": [
    "Similarly, numpy functions are broadcast across elements of an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b9f50d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b9f044",
   "metadata": {},
   "source": [
    "### Numeric linear algebra\n",
    "\n",
    "The transpose of an array `x` is given by `x.T`.\n",
    "\n",
    "Matrix multiplation is performed using the `@` operator. Here we perform $ M_{il} = \\sum_{k} x_{ik} * (x^T)_{kl}$, where $x^T$ is the transpose of $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0c6178",
   "metadata": {},
   "outputs": [],
   "source": [
    "x @ x.T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b33fcc2",
   "metadata": {},
   "source": [
    "To multiply two arrays element wise, ie to do $ M_{ij} = x_{ij} y_{ij}$, it's the usual multiplication operator `*`.\n",
    "\n",
    "Inverting matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0e2478",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.random.randint(9, size=(3, 3), dtype='int')\n",
    "b = a @ np.linalg.inv(a)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba258519",
   "metadata": {},
   "source": [
    "Computing the trace:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3907537a",
   "metadata": {},
   "outputs": [],
   "source": [
    "b.trace()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5246ba0",
   "metadata": {},
   "source": [
    "Determinant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d3a6e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.det(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c80fee0",
   "metadata": {},
   "source": [
    "Computing a Cholesky decomposition, i.e. finding lower triangular matrix $C$ such that $C C' = \\Sigma$ for $\\Sigma$ a 2-dimensional positive definite matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370af25f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Σ = np.array([[4, 1], [1, 3]])\n",
    "c = np.linalg.cholesky(Σ)\n",
    "\n",
    "c @ c.T - Σ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aeb0d0a",
   "metadata": {},
   "source": [
    "#### Solving systems of linear equations\n",
    "\n",
    "Say we have a system of equations, $4x + 3y + 2z = 25$,\n",
    "$-2x + 2y + 3z = -10$, and $3x -5y + 2z = -4$. We can solve these three equations for the three unknowns, x, y, and z, using the `solve` method. First, remember that this equation can be written in matrix form as\n",
    "\n",
    "$$\n",
    "M\\cdot \\vec{x} = \\vec{c}\n",
    "$$\n",
    "\n",
    "We can solve this by multiplying by the matrix inverse of $M$:\n",
    "\n",
    "$$\n",
    "M^{-1} M \\cdot \\vec{x} = I \\cdot \\vec{x} = M^{-1} \\cdot \\vec{c}\n",
    "$$\n",
    "\n",
    "which could be called by running `x = la.inv(M).dot(c)`. There's a convenience function in **numpy** called solve that does the same thing: here it finds the real values of the vector $\\vec{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8fcce7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[4, 3, 2], [-2, 2, 3], [3, -5, 2]])\n",
    "c = np.array([25, -10, -4])\n",
    "np.linalg.solve(M, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcdedda6",
   "metadata": {},
   "source": [
    "Finally, eigenvalues and eigenvectors can be found from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b82e345f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg as la\n",
    "eigvals, eigvecs = la.eig(M)\n",
    "eigvals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad28b24c",
   "metadata": {},
   "source": [
    "### Interpolation\n",
    "\n",
    "This section draws on the **scipy** documentation. There are built-in **pandas** methods for interpolation in dataframes, but **scipy** also has a range of functions for this including for univariate data `interp1d`, multidimensional interpolation on a grid `interpn`, `griddata` for unstructured data. Let's see a simple example with interpolation between a regular grid of integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107bbe55",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "x = np.arange(0, 10)\n",
    "y = np.exp(-x/3.0)\n",
    "f = interpolate.interp1d(x, y, kind='cubic')\n",
    "# Create a finer grid to interpolation function f\n",
    "xnew = np.arange(0, 9, 0.1)\n",
    "ynew = f(xnew)\n",
    "plt.plot(x, y, 'o', xnew, ynew, '-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10c3cef0",
   "metadata": {},
   "source": [
    "What about unstructured data? Let's create a Cobb-Douglas function on a detailed grid but then only retain a random set of the established points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3f41ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import griddata\n",
    "\n",
    "def cobb_doug(x, y):\n",
    "    alpha = 0.8\n",
    "    return x**(alpha)*y**(alpha-1)\n",
    "\n",
    "\n",
    "# Take some random points of the Cobb-Douglas function\n",
    "points = np.random.rand(1000, 2)\n",
    "values = cobb_doug(points[:,0], points[:,1])\n",
    "\n",
    "# Create a grid\n",
    "grid_x, grid_y = np.mgrid[0.01:1:200j, 0.01:1:200j]\n",
    "\n",
    "# Interpolate the points we have onto the grid\n",
    "interp_data = griddata(points, values, (grid_x, grid_y), method='cubic')\n",
    "\n",
    "# Plot results\n",
    "fig, axes = plt.subplots(1, 2)\n",
    "# Plot function & scatter of random points\n",
    "axes[0].imshow(cobb_doug(grid_x, grid_y).T,\n",
    "               extent=(0, 1, 0, 1),\n",
    "               origin='lower', cmap='plasma_r',\n",
    "               vmin=0, vmax=1)\n",
    "axes[0].plot(points[:, 0], points[:, 1], 'r.', ms=1.2)\n",
    "axes[0].set_title('Original + points')\n",
    "# Interpolation of random points\n",
    "axes[1].imshow(interp_data.T, extent=(0, 1, 0, 1),\n",
    "               origin='lower', cmap='plasma_r',\n",
    "               vmin=0, vmax=1)\n",
    "axes[1].set_title('Cubic interpolation');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09c7e3f7",
   "metadata": {},
   "source": [
    "### Optimisation\n",
    "\n",
    "**scipy** has functions for minimising scalar functions, minimising multivariate functions with complex surfaces, and root-finding. Let's see an example of finding the minimum of a scalar function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "309ca60c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "def f(x):\n",
    "    return x**2 + 10*np.sin(x) - 1.2\n",
    "\n",
    "\n",
    "result = optimize.minimize(f, x0=0)\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22c47af2",
   "metadata": {},
   "source": [
    "The result of the optimisation is in the 'x' attribute of `result`. Let's see this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05db5051",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-10, 10, 0.1)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, f(x))\n",
    "ax.scatter(result.x, f(result.x), s=150, color='k')\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('f(x)', rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4ae23e1",
   "metadata": {},
   "source": [
    "In higher dimensions, the minimisation works in much the same way, with the same function `optimize.minimize`. There are a LOT of minimisation options that you can pass to the `method=` keyword; the default is intelligently chosen from BFGS, L-BFGS-B, or SLSQP, depending upon whether you supply constraints or bounds.\n",
    "\n",
    "Root finding, aka solving equations of the form $f(x)=0$, is also catered for by **scipy**, through `optimize.root`. It works in much the same way as `optimizer.minimize`.\n",
    "\n",
    "In both of these cases, be warned that multiple roots and multiple minima can be hard to detect, and you may need to carefully specify the bounds or the starting positions in order to find the root you're looking for. Also, both of these methods can accept the Jacobian of the function you're working with as an argument, which is likely to improve performance with some solvers.\n",
    "\n",
    "### Numerical Integration\n",
    "\n",
    "**scipy** provides routines to numerically evaluate integrals in `scipy.integrate`, which you can find the documentation for [here](https://docs.scipy.org/doc/scipy/reference/integrate.html). Let's see an example using the 'vanilla' integration method, `quad`, to solve a known function between given (numerical) limits:\n",
    "\n",
    "$$\n",
    "\\displaystyle\\int_0^{\\pi} \\sin(x) d x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c075633f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "res, err = quad(np.sin, 0, np.pi)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c74920",
   "metadata": {},
   "source": [
    "What if we just have data samples? In that case, there are several routines that perform purely numerical integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c416809",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import simps\n",
    "x = np.arange(0, 10)\n",
    "f_of_x = np.arange(0, 10)\n",
    "\n",
    "simps(f_of_x, x) - 9**2/2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e021e70a",
   "metadata": {},
   "source": [
    "Even with just 10 evenly spaced points, the composite Simpson’s rule integration given by `simps` is able to accurately find the answer as $\\left( x^2/2\\right) |_{0}^{9}$.\n",
    "\n",
    "## Advanced: Composable Function Transformations\n",
    "\n",
    "In recent years, there have been great developments in the ability of Python to easily carry out numerical 'composable function transformations'. What this means is that, if you can dream up an arbitrary numerical operations -- including differentiation, linear algebra, and optimisation -- you can write code that will execute it quickly and automatically on CPUs, GPUs, or TPUs as you like.\n",
    "\n",
    "Here we'll look at one library that does this, **jax**, developed by Google {cite:t}`jax2018github`. It can automatically differentiate native Python and **numpy** functions, including when they are in loops, branches, or subject to recursion, and it can take derivatives of derivatives of derivatives. It supports reverse-mode differentiation (a.k.a. backpropagation) via grad as well as forward-mode differentiation, and the two can be composed arbitrarily to any order.\n",
    "\n",
    "To do these at speed, it uses just-in-time compilation. If you don't know what that is, don't worry: the details aren't important. It's just a way of getting close to C++ or Fortran speeds while still being able to write code in *much* more user friendly Python!\n",
    "\n",
    "### Auto-differentiation\n",
    "\n",
    "Let's see an example of auto-differentiation an arbitrary function. We'll write the definition of $\\tanh(x)$ as a function and evaluate it. Because we already imported a (symbolic) `tanh` function from Sympy above, we'll call the function below `tanh_num`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f8639f",
   "metadata": {},
   "source": [
    "```python\n",
    "from jax import grad\n",
    "import jax.numpy as jnp\n",
    "\n",
    "def tanh_num(θ):  # Define a function\n",
    "  y = jnp.exp(-2.0 * θ)\n",
    "  return (1.0 - y) / (1.0 + y)\n",
    "\n",
    "grad_tanh = grad(tanh_num)  # Obtain its gradient function\n",
    "grad_tanh(1.0)   # Evaluate it at x = 1.0\n",
    "```\n",
    "\n",
    "```bash\n",
    "DeviceArray(0.4199743, dtype=float32)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f68dba4",
   "metadata": {},
   "source": [
    "You can differentiate to any order using grad:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55c3a3d7",
   "metadata": {},
   "source": [
    "```python\n",
    "grad(grad(grad(tanh_num)))(1.0)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da81dc21",
   "metadata": {},
   "source": [
    "```bash\n",
    "DeviceArray(0.6216266, dtype=float32)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "922617dc",
   "metadata": {},
   "source": [
    "Let's check this using symbolic mathematics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf081820",
   "metadata": {},
   "outputs": [],
   "source": [
    "θ = Symbol(r'\\theta')\n",
    "triple_deriv = diff(diff(diff(tanh(θ), θ)))\n",
    "triple_deriv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a004d8d6",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "symp_est = triple_deriv.subs(θ, 1.)\n",
    "glue('symp_est', f'{symp_est:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e28aca64",
   "metadata": {},
   "source": [
    "If we evaluate this at $\\theta=1$, we get {glue:}`symp_est`. This was a simple example that had a (relatively) simple mathematical expression. But imagine if we had lots of branches (eg if, else statements), and/or a really complicated function: **jax**'s grad would still work. It's designed for really complex derivatives of the kind encountered in machine learning.\n",
    "\n",
    "### Just-in-time compilation\n",
    "\n",
    "The other nice feature of **jax** is the ability to do just-in-time (JIT) compilation. Because they do not compile their code into machine-code before running, high-level languages like Python and R are not as fast as the same code written in C++ or Fortran (the benefit is that it takes you less time to write the code in the first place). Much of the time, there are pre-composed functions that call C++ under the hood to do these things--but only for those operations that people have already taken the time to code up in a lower level language. JIT compilation offers a compromise: you can code more or less as you like in the high-level language but it will be compiled just-in-time to give you a speed-up!\n",
    "\n",
    "**jax** is certainly not the only Python package that does this, and if you're not doing anything like differentiating or propagating, **numba** is a more mature alternative. But here we'll see the time difference for JIT compilation on an otherwise slow operation: element wise multiplication and addition."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc3b6de",
   "metadata": {},
   "source": [
    "```python\n",
    "from jax import jit\n",
    "\n",
    "def slow_f(x):\n",
    "  \"\"\"Slow, element-wise function\"\"\"\n",
    "  return x * x + x * 2.0\n",
    "\n",
    "x = jnp.ones((5000, 5000))\n",
    "fast_f = jit(slow_f)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1cb94ca",
   "metadata": {},
   "source": [
    "Now let's see how fast the 'slow' version goes:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4efb1a14",
   "metadata": {},
   "source": [
    "```python\n",
    "%timeit -n15 -r3 slow_f(x)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee87c66e",
   "metadata": {},
   "source": [
    "```bash\n",
    "60.1 ms ± 3.67 ms per loop (mean ± std. dev. of 3 runs, 15 loops each)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5457c660",
   "metadata": {},
   "source": [
    "what about with the JIT compilation?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ce392e",
   "metadata": {},
   "source": [
    "```python\n",
    "%timeit -n15 -r3 fast_f(x)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68e46f84",
   "metadata": {},
   "source": [
    "```bash\n",
    "17.7 ms ± 434 µs per loop (mean ± std. dev. of 3 runs, 15 loops each)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cbb7b73",
   "metadata": {},
   "source": [
    "This short introduction has barely scratched the surface of **jax** and what you can do with it. For more, see the [official documentation](https://jax.readthedocs.io/en/latest/).\n",
    "\n",
    "## Set theory\n",
    "\n",
    "Set theory is a surprisingly useful tool in research (and invaluable in spatial analysis). Here are some really useful bits of set theory inspired by examples in {cite:t}`sheppard2012introduction`.\n",
    "\n",
    "Sets are first class citizens in Python in the same way that lists are. We can define and view a set like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f26b08e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = set(['Ada Lovelace', 'Sadie Alexander',\n",
    "         'Charles Babbage', 'Ada Lovelace',\n",
    "         'Adam Smith', 'Sadie Alexander'])\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88b634cf",
   "metadata": {},
   "source": [
    "Notice that a couple of entries appeared twice in the list but only once in the set: that's because a set contains only unique elements. Let's define a second set in order to demonstrate some of the operations we can perform on sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7cfae94",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = set(['Grace Hopper', 'Jean Bartik',\n",
    "         'Janet Yellen', 'Joan Robinson',\n",
    "         'Adam Smith', 'Ada Lovelace'])\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3a100db",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from myst_nb import glue\n",
    "inters = x.intersection(y)\n",
    "differ = x.difference(y)\n",
    "union = x.union(y)\n",
    "glue(\"inters\", inters)\n",
    "glue(\"differ\", differ)\n",
    "glue(\"union\", union)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58301fea",
   "metadata": {},
   "source": [
    "Now we have two sets we can look at to demonstrate some of the basic functions you can call on the set object type. `x.intersection(y)` gives, in this example, {glue:}`inters`, `x.difference(y)` gives {glue:}`differ`, and `x.union(y)` gives {glue:}`union`.\n",
    "\n",
    "**numpy** also has functions that use set theory. `np.unique` returns only the unique entries of an input array or list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05745e7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(['Lovelace', 'Hopper', 'Alexander', 'Hopper', 45, 27, 45])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87ad30e5",
   "metadata": {},
   "source": [
    "We can also ask which of a second set is a repeat of a first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "498be031",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(10)\n",
    "y = np.arange(5, 10)\n",
    "np.in1d(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "395bcfeb",
   "metadata": {},
   "source": [
    "And we have the numpy equivalents of intersection, `np.intersect1d(x, y)`, difference, `np.setdiff1d(x, y)`, and union, `np.union1d(x, y)`. Additionally, there is the exclusive-or (that I like to call 'xor'). This effectively returns the two arrays with their union removed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4703205b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3, 2, 4])\n",
    "b = np.array([2, 3, 5, 7, 5])\n",
    "np.setxor1d(a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ec3b8a2",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "In this chapter, you should have:\n",
    "\n",
    "- ✅ seen how to use symbolic algebra with code, including Lagrangrians and linear algebra;\n",
    "- ✅ seen how to code numerical mathematics, including linear algebra and optimisation; and\n",
    "- ✅ found out about using set theory via the `set` object type and set-oriented functions."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e82e9306cec5f3813db66e68c531c869f324183486c55e9589deb7f101067873"
  },
  "jupytext": {
   "cell_metadata_filter": "-all",
   "formats": "md:myst"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 64-bit ('codeforecon': conda)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}