04_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "###### <table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Sources of Error\n",
    "\n",
    "Error can come from many sources when formulating problems and/or applying numerical methods:\n",
    " - Model/Data Error\n",
    " - Discretization Error\n",
    " - Floating Point Error\n",
    " - Convergence Error\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Objectives\n",
    "\n",
    "* Understand the different sources of error\n",
    "* Explore some simple approaches to error analysis\n",
    "* Quantify errors\n",
    "    * absolute error\n",
    "    * relative error\n",
    "    * precision\n",
    "* Long term goals \n",
    "    * Use error estimates to control accuracy/reliability of solutions\n",
    "    * Understand errors so you can **believe** and **justify** your solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model and Data Error\n",
    "\n",
    "Errors in fundamental formulation\n",
    " - SIR model \n",
    "     - simplistic averaged model of interactions\n",
    "     $$ \\dot{I} = \\alpha SI - \\beta I$$\n",
    "     - basic model predicts a single peak \n",
    " - Data Error - Inaccuracy in measurement or uncertainties in parameters\n",
    " \n",
    "Unfortunately we cannot control model and data error directly but we can use methods that may be more robust in the presense of these types of errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Discretization or Truncation Error\n",
    "\n",
    "Errors arising from approximating a function with a simpler function, e.g. Using the approximation $\\sin(x) \\approx x$ when $|x| \\approx 0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Floating Point Error\n",
    "\n",
    "Errors arising from approximating real numbers with finite-precision numbers and arithmetic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Convergence Error\n",
    "\n",
    "In some instances an algorithm is developed that will take a current approximation and then find an improvement on the current approximation. In some instances the errors generated in each indivudal step can accumulate or become magnified after repeating the algorithm a number of times. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Basic Definitions\n",
    "\n",
    "Before exploring the different kinds of error, it is important to first define the ways that error is measured. Given a true value of a function $f$ and an approximate solution $F$ define:\n",
    "\n",
    "Absolute Error:  \n",
    "$$\n",
    "    e = | f - F |\n",
    "$$\n",
    "\n",
    "Relative Error:  \n",
    "$$\n",
    "    r = \\frac{e}{|f|} = \\frac{|f - F|}{|f|}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note: these definitions assume $f$ and $F$ are scalar valued.  However these definitions are readily extended to more complicate objects such as vectors or matrices with appropriate norms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Decimal Precision\n",
    "\n",
    "This definition of relative error provides a convenient estimate for the number of digits of decimal precision $p$\n",
    "\n",
    "given a relative error $r$,  the precision $p$ is the largest integer such that\n",
    "$$\n",
    "    r \\leq 5\\times 10^{-p}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Example\n",
    "* if $r = 0.001 < 5\\times10^{-3}$ has $p=3$ significant digits\n",
    "* if $r = 0.006 < 5\\times10^{-2}$ has $p=2$ significant digits (because this error would cause rounding up) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "let\n",
    "$$\n",
    "    f = e^1,\\quad F=2.71\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f = numpy.exp(1.)\n",
    "F = 2.71\n",
    "print('f = {}'.format(f))\n",
    "print('F = {}'.format(F))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "e = numpy.abs(f - F)\n",
    "r = e/numpy.abs(f)\n",
    "print('Absolute Error: {}'.format(e))\n",
    "print('Relative Error: {}'.format(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p = int(-numpy.log10(r/5.))\n",
    "print('Decimal precision: {}'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Big-O Notation\n",
    "\n",
    "In many situations an approximation will have a parameter associated with it, and the value of the parameter is often chosen to insure that the error is reasonable in a given situation. In such circumstances we often want to know the impact on the error if we change the value of the parameter. This leads to the definition of Big-O notation: \n",
    "\n",
    "$$\n",
    "    f(x) =  O(g(x)) \\quad \\text{as} \\quad x \\rightarrow a\n",
    "$$ \n",
    "if and only if \n",
    "$$\n",
    "    |f(x)| \\leq M |g(x)| \\quad \\text{as}\\quad  |x - a| < \\delta \\quad \\text{where} \\quad M,a > 0.\n",
    "$$ \n",
    "\n",
    "In practice we use Big-O notation to say something about how the terms we may have left out of a series might behave.  We saw an example earlier of this with the Taylor's series approximations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:\n",
    "Consider approximating a differentiable  function $f(x)$ by its Taylor polynomial (truncated Taylor series)  expanded around $x_0=0$., i.e. \n",
    "\n",
    "$$\n",
    "F(x) = T_N(x_0 + \\Delta x) = \\sum^N_{n=0} f^{(n)}(x_0)  \\frac{\\Delta x^n}{n!}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "where\n",
    "$$\n",
    "f(x)=\\lim_{N\\rightarrow\\infty} T_N\n",
    "$$ \n",
    "assuming the Taylor series converges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or for the case $f(x)=\\sin(x)$ expanded around $x_0=0$\n",
    "\n",
    "$$\n",
    "T_N(\\Delta x) = \\sum^N_{n=0} (-1)^{n} \\frac{\\Delta x^{2n+1}}{(2n+1)!}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For $N=2$, we can then write $F(x)$  as\n",
    "\n",
    "$$F(\\Delta x) = \\Delta x - \\frac{\\Delta x^3}{6} + \\frac{\\Delta x^5}{120}$$\n",
    "\n",
    "so our true function is\n",
    "\n",
    "$$\n",
    "    f(x) = F(\\Delta x) + O(\\Delta x^7)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or the absolute error\n",
    "\n",
    "$$\n",
    "    e = | f -F | \\sim O(\\Delta x^7)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**We can also develop rules for error propagation based on Big-O notation:**\n",
    "\n",
    "In general, there are two theorems that do not need proof and hold when the value of x is large: \n",
    "\n",
    "Let\n",
    "$$\\begin{aligned}\n",
    "    f(x) &= p(x) + O(x^n) \\\\\n",
    "    g(x) &= q(x) + O(x^m) \\\\\n",
    "    k &= \\max(n, m)\n",
    "\\end{aligned}$$\n",
    "then\n",
    "$$\n",
    "    f+g = p + q + O(x^k)\n",
    "$$\n",
    "and\n",
    "\\begin{align}\n",
    "    f \\cdot g &= p \\cdot q + p O(x^m) + q O(x^n) + O(x^{n + m}) \\\\\n",
    "    &= p \\cdot q + O(x^{n+m})\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "On the other hand, if we are interested in small values of x, say $\\Delta x$, the above expressions can be modified as follows: \n",
    "\n",
    "\\begin{align}\n",
    "    f(\\Delta x) &= p(\\Delta x) + O(\\Delta x^n) \\\\\n",
    "    g(\\Delta x) &= q(\\Delta x) + O(\\Delta x^m) \\\\\n",
    "    r &= \\min(n, m)\n",
    "\\end{align}\n",
    "then\n",
    "$$\n",
    "    f+g = p + q + O(\\Delta x^r)\n",
    "$$\n",
    "and\n",
    "\\begin{align}\n",
    "    f \\cdot g &= p \\cdot q + p \\cdot O(\\Delta x^m) + q \\cdot O(\\Delta x^n) + O(\\Delta x^{n+m}) \\\\\n",
    "    &= p \\cdot q + O(\\Delta x^r)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "**Note:** In this case we suppose that at least the polynomial with $k = \\max(n, m)$ has the following form: \n",
    "\n",
    "$$\n",
    "    p(\\Delta x) = 1 + p_1 \\Delta x + p_2 \\Delta x^2 + \\ldots\n",
    "$$\n",
    "or \n",
    "$$\n",
    "    q(\\Delta x) = 1 + q_1 \\Delta x + q_2 \\Delta x^2 + \\ldots\n",
    "$$\n",
    "\n",
    "so that there is an $\\mathcal{O}(1)$ term that guarantees the existence of $\\mathcal{O}(\\Delta x^r)$ in the final product. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To get a sense of why we care most about the power on $\\Delta x$ when considering convergence the following figure shows how different powers on the convergence rate can effect how quickly we converge to our solution.  Note that here we are plotting the same data two different ways.  Plotting the error as a function of $\\Delta x$ is a common way to show that a numerical method is doing what we expect and exhibits the correct convergence behavior.  Since errors can get small quickly it is very common to plot these sorts of plots on a log-log scale to easily visualize the results.  Note that if a method was truly of the order $n$ that they will be a linear function in log-log space with slope $n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### Behavior of error as a function of $\\Delta x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "dx = numpy.linspace(1.0, 1e-4, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2.0)\n",
    "axes = []\n",
    "axes.append(fig.add_subplot(1, 2, 1))\n",
    "axes.append(fig.add_subplot(1, 2, 2))\n",
    "\n",
    "for n in range(1, 5):\n",
    "    axes[0].plot(dx, dx**n, label=\"$\\Delta x^%s$\" % n)\n",
    "    axes[1].loglog(dx, dx**n, label=\"$\\Delta x^%s$\" % n)\n",
    "\n",
    "axes[0].legend(loc=2)\n",
    "axes[1].set_xticks([10.0**(-n) for n in range(5)])\n",
    "axes[1].set_yticks([10.0**(-n) for n in range(16)])\n",
    "axes[1].legend(loc=4)\n",
    "for n in range(2):\n",
    "    axes[n].set_title(\"Growth of Error vs. $\\Delta x^n$\")\n",
    "    axes[n].set_xlabel(\"$\\Delta x$\")\n",
    "    axes[n].set_ylabel(\"Estimated Error\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Discretization Error\n",
    "\n",
    "**Taylor's Theorem:**  Let $f(x) \\in C^{N+1}[a,b]$ and $x_0 \\in [a,b]$, then for all $x \\in (a,b)$ there exists a number $c = c(x)$ that lies between $x_0$ and $x$ such that\n",
    "\n",
    "$$ f(x) = T_N(x) + R_N(x)$$\n",
    "\n",
    "where $T_N(x)$ is the Taylor polynomial approximation\n",
    "\n",
    "$$T_N(x) = \\sum^N_{n=0} \\frac{f^{(n)}(x_0)\\cdot(x-x_0)^n}{n!}$$\n",
    "\n",
    "and $R_N(x)$ is the residual (the part of the series we left off)\n",
    "\n",
    "$$R_N(x) = \\frac{f^{(N+1)}(c) \\cdot (x - x_0)^{N+1}}{(N+1)!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Note\n",
    "\n",
    "\n",
    "The residual:\n",
    "\n",
    "$$\n",
    "    R_N(x) = \\frac{f^{(N+1)}(c) \\cdot (x - x_0)^{N+1}}{(N+1)!}\n",
    "$$\n",
    "\n",
    "depends on the $N+1$ order derivative of $f$ evaluated at an **unknown** value $c\\in[x,x_0]$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If we knew the value of $c$ we would know the exact value of $R_N(x)$ and therefore the function $f(x)$.  In general we don't know this value but we can use $R_N(x)$ to put upper bounds on the error **and** to understand how the error changes as we move away from $x_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Start by replacing $x - x_0$ with $\\Delta x$.  The primary idea here is that the residual $R_N(x)$ becomes smaller as $\\Delta x \\rightarrow 0$ (at which point $T_N(x) = f(x_0)$).\n",
    "\n",
    "$$\n",
    "    T_N(x) = \\sum^N_{n=0} \\frac{f^{(n)}(x_0)\\cdot\\Delta x^n}{n!}\n",
    "$$\n",
    "\n",
    "and $R_N(x)$ is the residual (the part of the series we left off)\n",
    "\n",
    "$$\n",
    "    R_N(x) = \\frac{f^{(n+1)}(c) \\cdot \\Delta x^{n+1}}{(n+1)!} \\leq M \\Delta x^{n+1} = O(\\Delta x^{n+1})\n",
    "$$\n",
    "\n",
    "where $M$ is an upper bound on \n",
    "$$\n",
    "    \\frac{f^{(n+1)}(c)}{(n+1)!}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 1\n",
    "\n",
    "$f(x) = e^x$ with $x_0 = 0$ on the interval $x\\in(-1,1)$\n",
    "\n",
    "Using this we can find expressions for the relative and absolute error as a function of $x$ assuming $N=2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Derivatives:\n",
    "$$\\begin{aligned}\n",
    "    f'(x) &= e^x \\\\\n",
    "    f''(x) &= e^x \\\\ \n",
    "    f^{(n)}(x) &= e^x\n",
    "\\end{aligned}$$\n",
    "\n",
    "Taylor polynomials:\n",
    "$$\\begin{aligned}\n",
    "    T_N(x) &= \\sum^N_{n=0} e^0 \\frac{x^n}{n!} \\Rightarrow \\\\\n",
    "    T_2(x) &= 1 + x + \\frac{x^2}{2}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Remainders:\n",
    "$$\\begin{aligned}\n",
    "    R_N(x) &= e^c \\frac{x^{N+1}}{(N+1)!} \\\\\n",
    "    R_2(x) &= e^c \\cdot \\frac{x^3}{6} \\leq \\frac{e^1}{6} \\approx 0.5\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Accuracy:\n",
    "\\begin{align}\n",
    "    \\exp(1) &= 2.718\\ldots \\\\\n",
    "    T_2(1) &= 2.5 \n",
    "\\end{align}\n",
    "\n",
    "$$\n",
    "\\Rightarrow e \\approx 0.2,\\quad r \\approx 0.08,\\quad p = ?\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also use the package sympy which has the ability to calculate Taylor polynomials built-in!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing(pretty_print=True)\n",
    "x = sympy.symbols('x')\n",
    "f = sympy.exp(x)\n",
    "f.series(x0=0, n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Lets plot this numerically for a section of $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-1, 1, 100)\n",
    "f = numpy.exp(x)\n",
    "T_N = 1.0 + x + x**2 / 2.0\n",
    "R_N = numpy.exp(1) * x**3 / 6.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x, T_N, 'r', x, f, 'k', x, numpy.abs(R_N), 'b')\n",
    "axes.plot(x,numpy.abs(numpy.exp(x)-T_N),'g--')\n",
    "axes.plot(0.0, 1.0, 'o', markersize=10)\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"$f(x)$, $T_N(x)$, $|R_N(x)|$\", fontsize=16)\n",
    "axes.legend([\"$T_N(x)$\", \"$f(x)$\", \"$|R_N(x)|$\", \"e(x)\"], loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example 2\n",
    "\n",
    "Approximate\n",
    "$$\n",
    "    f(x) = \\frac{1}{x} \\quad x_0  = 1,\n",
    "$$\n",
    "using $x_0 = 1$ to the 3rd Taylor series term on the inverval $x\\in[1,\\infty)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "\\begin{matrix}\n",
    "f'(x) = -\\frac{1}{x^2}, &  f''(x) = \\frac{2}{x^3}, & f'''(x) = -\\frac{6}{x^4}, & \\ldots, & f^{(n)}(x) &= \\frac{(-1)^n n!}{x^{n+1}}\n",
    "\\end{matrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    T_N(x) &= \\sum^N_{n=0} (-1)^n (x-1)^n \\Rightarrow \\\\\n",
    "    T_2(x) &= 1 - (x - 1) + (x - 1)^2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    R_N(x) &= \\frac{(-1)^{n+1}(x - 1)^{n+1}}{c^{n+2}} \\Rightarrow \\\\\n",
    "    R_2(x) &= \\frac{-(x - 1)^{3}}{c^{4}}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### plot this problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.8, 2, 100)\n",
    "f = 1.0 / x\n",
    "T_N = 1.0 - (x-1) + (x-1)**2\n",
    "R_N = -(x-1.0)**3 / (1.**4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(x, T_N, 'r', x, f, 'k', x, numpy.abs(R_N), 'b')\n",
    "plt.plot(x,numpy.abs(f - T_N),'g--')\n",
    "plt.plot(1.0, 1.0, 'o', markersize=10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"x\",fontsize=16)\n",
    "plt.ylabel(\"$f(x)$, $T_N(x)$, $R_N(x)$\",fontsize=16)\n",
    "plt.title('$f(x) = 1/x$',fontsize=18)\n",
    "plt.legend([\"$T_N(x)$\", \"$f(x)$\", \"$|R_N(x)|$\", '$e(x)$'], loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Computational Issue #1: Accuracy... how many terms?\n",
    "\n",
    "Given a Taylor Polynomial approximation of an arbitrary function $f(x)$,  how do we determine how many terms are required such that $R_N(x)<tol$.  And how do we determine the tolerance?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computational Issue #2 Efficiency... Operation counts for polynomial evaluation\n",
    "\n",
    "Given \n",
    "\n",
    "$$P_N(x) = a_0 + a_1 x + a_2 x^2 + \\ldots + a_N x^N$$ \n",
    "\n",
    "or\n",
    "\n",
    "$$P_N(x) = p_0 x^N + p_1 x^{N-1} + p_2 x^{N-2} + \\ldots + p_{N}$$\n",
    "\n",
    "what is the most **efficient way**  to evaluate $P_N(x)$? (i.e. minimize number of floating point operations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider two ways to write $P_3$\n",
    "\n",
    "* The standard way:\n",
    "\n",
    "$$ P_3(x) = p_0 x^3 + p_1 x^2 + p_2 x + p_3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* using nested multiplication (aka **Horner's Method**):\n",
    "\n",
    "$$ P_3(x) = ((p_0 x + p_1) x + p_2) x + p_3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider how many operations it takes for each...\n",
    "\n",
    "$$ P_3(x) = p_0 x^3 + p_1 x^2 + p_2 x + p_3$$\n",
    "\n",
    "$$P_3(x) = \\overbrace{p_0 \\cdot x \\cdot x \\cdot x}^3 + \\overbrace{p_1\\cdot x \\cdot x}^2 + \\overbrace{p_2 \\cdot x}^1 + p_3$$\n",
    "\n",
    "Note:  here we're just counting multiplications as they will dominate the flop count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Adding up all the operations we can in general think of this as a pyramid (it's really the triangle numbers)\n",
    "\n",
    "$$\n",
    "    \\sum_{n=1}^N n = \\frac{N(N+1)}{2} \n",
    "$$\n",
    "\n",
    "![Original Count](./images/horners_method_big_count.png)\n",
    "\n",
    "Thus we can estimate that the algorithm written this way will take approximately $O(N^2 / 2)$ operations to complete."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Looking at nested iteration, however:\n",
    "\n",
    "$$ P_3(x) = ((p_0 x + p_1) x + p_2) x + p_3$$\n",
    "\n",
    "Here we find that the method is $O(N)$ compared to the first evaluation which $O(N^2)$ (we usually drop the 2 in these cases).  That's a huge difference for large $N$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "\n",
    "Fill in the function and implement Horner's method:\n",
    "```python\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x should be a float.\n",
    "    \"\"\"\n",
    "    pass\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x should be a float.\n",
    "    \"\"\"\n",
    "    \n",
    "    ### YOUR CODE HERE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = [1, 2, 3]\n",
    "x = 1.\n",
    "test = eval_poly(p,x)\n",
    "answer = numpy.array([x**2, x, 1]).dot(p)\n",
    "print('test = {}, answer = {}'.format(test,answer))\n",
    "numpy.testing.assert_allclose(test,answer)\n",
    "print('success')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Scalar version\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x should be a float.\n",
    "    \"\"\"\n",
    "    \n",
    "    y = p[0]\n",
    "    for coefficient in p[1:]:\n",
    "        y = y * x + coefficient\n",
    "    \n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Vectorized version\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x can by a NumPy ndarray.\n",
    "    \"\"\"\n",
    "    \n",
    "    y = numpy.ones(x.shape) * p[0]\n",
    "    for coefficient in p[1:]:\n",
    "        y = y * x + coefficient\n",
    "    \n",
    "    return y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "p = [1, -3, 10, 4, 5, 5]\n",
    "x = numpy.linspace(-10, 10, 100)\n",
    "plt.plot(x, eval_poly(p, x))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence Error\n",
    "\n",
    "In some circumstances a formula or algorithm is applied repeatedly as a way to obtain a final approximation. Usually, the errors that occur at each individual step are small. By repeating the algorithm, though, the errors can sometimes grow or become magnified. \n",
    "\n",
    "As example of this phenomena is given below. The values of the terms in a difference equation are calculated,\n",
    "$$\n",
    "   \\begin{align}\n",
    "      y_0 &= 1, \\\\\n",
    "      y_1 &= \\frac{1}{5}, \\\\\n",
    "      y_{n+1} &= \\frac{16}{5} y_n - \\frac{3}{5} y_{n-1}.\n",
    "   \\end{align}\n",
    "$$\n",
    "\n",
    "The true solution to the difference equation is $y_n = \\left(\\frac{1}{5}\\right)^n$, where $n=$0, 1, 2, $\\ldots$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Choose the number of iterations\n",
    "N = 40\n",
    "y = numpy.empty(N+1)            # Allocate an empty vector with N+1 entries\n",
    "\n",
    "# Now use the difference equation to generate the numbers in the sequence\n",
    "y[0] = 1\n",
    "y[1] = 1/5\n",
    "for n in range(2,N+1):\n",
    "    y[n] = 16/5*y[n-1] - 3/5*y[n-2]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And plot the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Now plot the result\n",
    "n = numpy.arange(0,N+1)\n",
    "fig = plt.figure(figsize=(10.0, 5.0))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(n,y, 'rx', markersize=5, label='$y_n$')\n",
    "axes.semilogy(n,(1/5)**n,'b.', label='$y_{true}$')\n",
    "axes.grid()\n",
    "axes.set_title(\"Calculated Values Of A Difference Equation\",fontsize=18)\n",
    "axes.set_xlabel(\"$n$\",fontsize=16)\n",
    "axes.set_ylabel(\"$y_n$\",fontsize=16)\n",
    "axes.legend(loc='best', shadow=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Simply looking at the exact solution, the sequence of numbers generated by the difference equation above should get very close to zero. Instead, the numbers in the sequence initially get closer to zero, but at some point they begin to grow and get larger. An underlying problem is that the computer is not able to store the numbers exactly. The second number in the sequence, $y_1=\\frac{1}{5}$ has a small error, and the computer stores it as $y_1 = \\frac{1}{5}+\\epsilon$ where $\\epsilon$ is some small error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Each time a new number in the loop is generated, the error is multiplied. For example, after the first iteration $y_2$ is\n",
    "$$\n",
    "    \\begin{align}\n",
    "       y_2 &= \\frac{16}{5} \\left( \\frac{1}{5}+\\epsilon \\right) \n",
    "               - \\frac{3}{5} \\left( 1 \\right), \\\\\n",
    "           &= \\frac{1}{5^2} + \\frac{16}{5} \\epsilon.\n",
    "    \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "After the second time through the loop, the value of $y_3$ is \n",
    "\n",
    "$$ \n",
    "y_3=\\frac{1}{5^3} + \\frac{241}{25}\\epsilon\n",
    "$$\n",
    "\n",
    "Even though the value of $\\epsilon$ is very close to zero, every iteration makes the error grow.  Repeated multiplication will result in a very large number. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The error associated with the initial representation of the number $\\frac{1}{5}$ is a problem with the way a digital computer stores floating point numbers. In most instances the computer cannot represent a number exactly, and the small error in approximating a given number can give rise to other problems. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Floating Point Error\n",
    "\n",
    "Errors arising from approximating real numbers with finite-precision numbers\n",
    "\n",
    "$$\\pi \\approx 3.14$$\n",
    "\n",
    "or $\\frac{1}{3} \\approx 0.333333333$ in decimal, results from a finite number of registers to represent each number.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Floating Point Systems\n",
    "\n",
    "Numbers in floating point systems are represented as a series of bits that represent different pieces of a number.  In *normalized floating point systems* there are some standard conventions for what these bits are used for.  In general the numbers are stored by breaking them down into the form  \n",
    "\n",
    "$$F = \\pm d_1 . d_2 d_3 d_4 \\ldots d_p \\times \\beta^E$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "where\n",
    "1. $\\pm$ is a single bit  representing the sign of the number\n",
    "2. $d_1 . d_2 d_3 d_4 \\ldots d_p$ is called the *mantissa*.  Note that technically the decimal could be moved but generally, using scientific notation, the decimal can always be placed at this location.  The digits $d_2 d_3 d_4 \\ldots d_p$ are called the *fraction* with $p$ digits of precision.  Normalized systems specifically put the decimal point in the front like we have and assume $d_1 \\neq 0$ unless the number is exactly $0$.\n",
    "3. $\\beta$ is the *base*.  For binary $\\beta = 2$, for decimal $\\beta = 10$, etc.\n",
    "4. $E$ is the *exponent*, an integer in the range $[E_{\\min}, E_{\\max}]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The important points on any floating point system is that\n",
    "1. There exist a discrete and finite set of representable numbers\n",
    "2. These representable numbers are not evenly distributed on the real line\n",
    "3. Arithmetic in floating point systems yield different results from infinite precision arithmetic (i.e. \"real\" math)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Properties of Floating Point Systems\n",
    "All floating-point systems are characterized by several important numbers\n",
    " - Smalled normalized number (underflow if below - related to subnormal numbers around zero)\n",
    " - Largest normalized number (overflow if above)\n",
    " - Zero\n",
    " - Machine $\\epsilon$ or $\\epsilon_{\\text{machine}}$\n",
    " - `inf` and `nan`, infinity and **N**ot **a** **N**umber respectively"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##### Example:  Toy System\n",
    "Consider the toy 2-digit precision decimal system (normalized)\n",
    "$$f = \\pm d_1 . d_2 \\times 10^E$$\n",
    "with $E \\in [-2, 0]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Number and distribution of numbers**\n",
    "1. How many numbers can we represent with this system?\n",
    "\n",
    "2. What is the distribution on the real line?\n",
    "\n",
    "3. What is the underflow and overflow limits?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "How many numbers can we represent with this system?\n",
    "\n",
    "$$\n",
    "    f = \\pm d_1 . d_2 \\times 10^E ~~~ \\text{with} E \\in [-2, 0]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* sign bit: 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* $d_1$: 9   (normalized numbers $d_1\\neq 0$)\n",
    "* $d_2$: 10  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* $E$: 3 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* zero: 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "total:\n",
    "$$ \n",
    "    2 \\times 9 \\times 10 \\times 3 + 1 = 541\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What is the distribution on the real line? $$f = \\pm d_1 . d_2 \\times 10^E ~~~ \\text{with} ~~~ E \\in [-2, 0]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "d_1_values = [1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
    "d_2_values = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
    "E_values = [0, -1, -2,]\n",
    "\n",
    "fig = plt.figure(figsize=(10.0,3.0))\n",
    "axes = fig.add_subplot(2, 1, 1)\n",
    "\n",
    "for E in E_values:\n",
    "    for d1 in d_1_values:\n",
    "        for d2 in d_2_values:\n",
    "            axes.plot( (d1 + d2 * 0.1) * 10**E, 0.0, 'r|', markersize=20)\n",
    "            axes.plot(-(d1 + d2 * 0.1) * 10**E, 0.0, 'r|', markersize=20)\n",
    "            \n",
    "axes.plot(0.0, 0.0, '|', markersize=20)\n",
    "axes.plot([-1., 1.], [0.0, 0.0], 'k|', markersize=30)\n",
    "\n",
    "axes.plot([-10.0, 10.0], [0.0, 0.0], 'k')\n",
    "\n",
    "axes.set_title(\"Distribution of Values $[-10, 10]$\")\n",
    "axes.set_yticks([])\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"\")\n",
    "#axes.set_xlim([-10, 10])\n",
    "\n",
    "axes = fig.add_subplot(2,1,2)\n",
    "for E in E_values:\n",
    "    for d1 in d_1_values:\n",
    "        for d2 in d_2_values:\n",
    "            axes.plot( (d1 + d2 * 0.1) * 10**E, 0.0, 'r+', markersize=20)\n",
    "            axes.plot(-(d1 + d2 * 0.1) * 10**E, 0.0, 'r+', markersize=20)\n",
    "            \n",
    "axes.plot(0.0, 0.0, '+', markersize=20)\n",
    "axes.plot([-10.0, 10.0], [0.0, 0.0], 'k')\n",
    "\n",
    "axes.set_title(\"Close up $[-1, 1]$\")\n",
    "axes.set_yticks([])\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"\")\n",
    "axes.set_xlim([-1, 1])\n",
    "fig.tight_layout(h_pad=1, w_pad=5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the underflow and overflow limits?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Smallest number that can be represented is the underflow:  $1.0 \\times 10^{-2} = 0.01$\n",
    "* Largest number that can be represented is the overflow:  $9.9 \\times 10^0 = 9.9$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the smallest number $\\epsilon_{mach}$ such that $1+\\epsilon_{mach} > 1$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* $\\epsilon_{mach} = 0.1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Binary Systems\n",
    "Consider the 2-digit precision base 2 system:\n",
    "\n",
    "$$\n",
    "    f=\\pm d_1 . d_2 \\times 2^E \\quad \\text{with} \\quad E \\in [-1, 1]\n",
    "$$\n",
    "\n",
    "#### Number and distribution of numbers\n",
    "1. How many numbers can we represent with this system?\n",
    "\n",
    "2. What is the distribution on the real line?\n",
    "\n",
    "3. What is the underflow and overflow limits?\n",
    "\n",
    "4. What is $\\epsilon_{mach}$?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "How many numbers can we represent with this system?\n",
    "\n",
    "$$f=\\pm d_1 . d_2 \\times 2^E ~~~~ \\text{with} ~~~~ E \\in [-1, 1]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$ 2 \\times 1 \\times 2 \\times 3 + 1 = 13$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the distribution on the real line?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "d_1_values = [1]\n",
    "d_2_values = [0, 1]\n",
    "E_values = [1, 0, -1]\n",
    "\n",
    "fig = plt.figure(figsize=(10.0, 1.0))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "for E in E_values:\n",
    "    for d1 in d_1_values:\n",
    "        for d2 in d_2_values:\n",
    "            axes.plot( (d1 + d2 * 0.5) * 2**E, 0.0, 'r+', markersize=20)\n",
    "            axes.plot(-(d1 + d2 * 0.5) * 2**E, 0.0, 'r+', markersize=20)\n",
    "            \n",
    "axes.plot(0.0, 0.0, 'r+', markersize=20)\n",
    "axes.plot([-4.5, 4.5], [0.0, 0.0], 'k')\n",
    "\n",
    "axes.set_title(\"Distribution of Values\")\n",
    "axes.set_yticks([])\n",
    "axes.set_xticks(numpy.linspace(-4,4,9))\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"\")\n",
    "axes.grid()\n",
    "axes.set_xlim([-5, 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Smallest number that can be represented is the underflow:  $1.0 \\times 2^{-1} = 0.5$\n",
    "\n",
    "* Largest number that can be represented is the overflow:  $1.1 \\times 2^1 = 3$\n",
    "\n",
    "* $\\epsilon_{mach} = 0.1 = 2^{-1}= 1/2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Note**: these numbers are in a binary system.  \n",
    "\n",
    "Quick rule of thumb:\n",
    "$$\n",
    "    2^3 2^2 2^1 2^0 . 2^{-1} 2^{-2} 2^{-3}\n",
    "$$\n",
    "correspond to\n",
    "8s, 4s, 2s, 1s . halves, quarters, eighths, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Real Systems - IEEE 754 Binary Floating Point Systems\n",
    "\n",
    "#### Single Precision\n",
    " - Total storage alloted is 32 bits\n",
    " - Exponent is 8 bits $\\Rightarrow E \\in [-126, 128]$\n",
    " - Fraction 23 bits ($p = 24$)\n",
    " \n",
    "```\n",
    "s EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF\n",
    "0 1      8 9                     31\n",
    "```\n",
    "* Overflow $= 2^{127} \\approx 3.4 \\times 10^{38}$\n",
    "* Underflow $= 2^{-126} \\approx 1.2 \\times 10^{-38}$\n",
    "* $\\epsilon_{\\text{machine}} = 2^{-23} \\approx 1.2 \\times 10^{-7}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Double Precision\n",
    " - Total storage alloted is 64 bits\n",
    " - Exponent is 11 bits $\\Rightarrow E \\in [-1022, 1024]$\n",
    " - Fraction 52 bits ($p = 53$)\n",
    " \n",
    "```\n",
    "s EEEEEEEEEE FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FF\n",
    "0 1       11 12                                                      63\n",
    "```\n",
    "* Overflow $= 2^{1024} \\approx 1.8 \\times 10^{308}$\n",
    "* Underflow $= 2^{-1022} \\approx 2.2 \\times 10^{-308}$\n",
    "* $\\epsilon_{\\text{machine}} = 2^{-52} \\approx 2.2 \\times 10^{-16}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Python Access to IEEE Numbers\n",
    "\n",
    "Access many important parameters, such as machine epsilon:\n",
    "\n",
    "```python\n",
    "import numpy\n",
    "numpy.finfo(float).eps\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.finfo(numpy.float16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.finfo(numpy.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.finfo(float))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.finfo(numpy.float128))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "eps = numpy.finfo(float).eps\n",
    "MAX = numpy.finfo(float).max\n",
    "print('eps = {}'.format(eps))\n",
    "print('MAX = {}'.format(MAX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Show that $(1 + \\epsilon_{mach}) > 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(MAX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(MAX*(1+ .4*eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(1 + .4*eps == 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why should we care about this?\n",
    "\n",
    " - Floating point arithmetic is not commutative or associative\n",
    " - Floating point errors compound, do not assume even double precision is enough!\n",
    " - Mixing precision can be very dangerous"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 1: Simple Arithmetic\n",
    " \n",
    "Simple arithmetic $\\delta < \\epsilon_{\\text{machine}}$. \n",
    "\n",
    "Compare\n",
    "\n",
    "   $$1+\\delta - 1 \\quad vs. \\quad 1 - 1 + \\delta$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "eps = numpy.finfo(float).eps\n",
    "delta = 0.5*eps\n",
    "x = 1+delta  -1\n",
    "y = 1 - 1 + delta\n",
    "print('1 + delta - 1 = {}'.format(x))\n",
    "print('1 - 1 + delta = {}'.format(y))\n",
    "print( x == y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 2: Catastrophic cancellation \n",
    "\n",
    "Let us examine what happens when we add two numbers $x$ and $y$ where $x + y \\neq 0$.  We can actually estimate these bounds by doing some error analysis.  Here we need to introduce the idea that each floating point operation introduces an error such that\n",
    "\n",
    "$$\n",
    "    \\text{fl}(x ~\\text{op}~ y) = (x ~\\text{op}~ y) (1 + \\delta)\n",
    "$$\n",
    "\n",
    "where $\\text{fl}(\\cdot)$ is a function that returns the floating point representation of the expression enclosed, $\\text{op}$ is some operation (e.g. $+, -, \\times, /$), and $\\delta$ is the floating point error due to $\\text{op}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Back to our problem at hand.  The floating point error due to addition is\n",
    "\n",
    "$$\n",
    "    \\text{fl}(x + y) = (x + y) (1 + \\delta).\n",
    "$$\n",
    "\n",
    "Comparing this to the true solution using a relative error we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{|(x + y) - \\text{fl}(x + y)|}{|x + y|} &= \\frac{|(x + y) - (x + y) (1 + \\delta)|}{|x + y|} = \\delta.\n",
    "\\end{aligned}$$\n",
    "\n",
    "so that if $\\delta = \\mathcal{O}(\\epsilon_{\\text{machine}})$ we are not too concerned."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What if instead we consider a floating point error on the representations of $x$ and $y$, $x \\neq y$, and say $\\delta_x$ and $\\delta_y$ are the magnitude of the errors in their representation.  Here we will assume this constitutes the floating point error rather than being associated with the operation itself.\n",
    "\n",
    "Now consider the difference between the two numbers\n",
    "$$\\begin{aligned}\n",
    "    \\text{fl}(x - y) &= x (1 + \\delta_x) - y (1 + \\delta_y) \\\\\n",
    "    &= x - y + x \\delta_x - y \\delta_y \\\\\n",
    "    &= (x - y) \\left(1 + \\frac{x \\delta_x - y \\delta_y}{x - y}\\right)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again computing the relative error we then have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{\\left|(x - y) - (x - y) \\left(1 + \\frac{x \\delta_x - y \\delta_y}{x - y}\\right)\\right|}{|x - y|} &= \n",
    "   \\left|1 - \\left(1 + \\frac{x \\delta_x - y \\delta_y}{x - y}\\right)\\right|\\\\\n",
    "   &=\\frac{|x \\delta_x - y \\delta_y|}{|x - y|} \\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "The important distinction here is that now the error is dependent on the values of $x$ and $y$ and more importantly, their difference.  Of particular concern is the relative size of $x - y$.  As it approaches zero relative to the magnitudes of $x$ and $y$ the error could be arbitrarily large.  This is known as **catastrophic cancellation**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "dx = numpy.array([10**(-n) for n in range(1, 16)])\n",
    "x = 1.0 + dx\n",
    "y = numpy.ones(x.shape)\n",
    "error = numpy.abs(x - y - dx) / (dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.loglog(dx, x + y, 'o-')\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$x + y$\")\n",
    "axes.set_title(\"$\\Delta x$ vs. $x+y$\")\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.loglog(dx, error, 'o-')\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$|x + y - \\Delta x| / \\Delta x$\")\n",
    "axes.set_title(\"Difference between $x$ and $y$ vs. Relative Error\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 3: Function Evaluation\n",
    "\n",
    "Consider the function\n",
    "$$\n",
    "    f(x) = \\frac{1 - \\cos x}{x^2}\n",
    "$$\n",
    "with $x\\in[-10^{-4}, 10^{-4}]$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Taking the limit as $x \\rightarrow 0$ we can see what behavior we would expect to see from evaluating this function:\n",
    "$$\n",
    "    \\lim_{x \\rightarrow 0} \\frac{1 - \\cos x}{x^2} = \\lim_{x \\rightarrow 0} \\frac{\\sin x}{2 x} = \\lim_{x \\rightarrow 0} \\frac{\\cos x}{2} = \\frac{1}{2}.\n",
    "$$\n",
    "\n",
    "What does floating point representation do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-1e-3, 1e-3, 100, dtype=numpy.float32)\n",
    "f = 0.5\n",
    "F = (1.0 - numpy.cos(x)) / x**2\n",
    "rel_err = numpy.abs((f - F)) / f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, rel_err, 'o')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.grid()\n",
    "axes.set_ylabel(\"Relative Error\")\n",
    "axes.set_title(\"$\\\\frac{1-\\\\cos{x}}{x^2} - \\\\frac{1}{2}$\",fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 4: Evaluation of a Polynomial\n",
    "\n",
    "   $$f(x) = x^7 - 7x^6 + 21 x^5 - 35 x^4 + 35x^3-21x^2 + 7x - 1$$\n",
    "   \n",
    "Note: $f(1) = 0$ (and will be close to zero for $x\\approx 1$)\n",
    "\n",
    "Here we compare polynomial evaluation using naive powers compared to Horner's method as implemented in `eval_poly(p,x)` defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.988, 1.012, 1000, dtype=numpy.float16)\n",
    "y = x**7 - 7.0 * x**6 + 21.0 * x**5 - 35.0 * x**4 + 35.0 * x**3 - 21.0 * x**2 + 7.0 * x - 1.0\n",
    "\n",
    "# repeat using Horner's method from above\n",
    "p = numpy.array([1, -7, 21, -35, 35, -21, 7, -1 ])\n",
    "yh = eval_poly(p,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, y, 'r',label='naive')\n",
    "axes.plot(x, yh, 'b',label='horner')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_ylim((-0.1, 0.1))\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "axes.grid()\n",
    "axes.legend()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x,yh-y,'g')\n",
    "axes.grid()\n",
    "axes.set_xlabel('x')\n",
    "axes.set_ylabel('$f_{horner} - f_n$')\n",
    "axes.set_title('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 5: Rational Function Evaluation\n",
    "Compute $f(x) = x + 1$ by the function \n",
    "\n",
    "$$F(x) = \\frac{x^2 - 1}{x - 1}$$\n",
    "\n",
    "Do you expect there to be issues?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.5, 1.5, 101, dtype=numpy.float64)\n",
    "f_hat = (x**2 - 1.0) / (x - 1.0)\n",
    "f = (x + 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.abs(f - f_hat)/numpy.abs(f))\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"Relative Error\")\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Combination of Errors\n",
    "\n",
    "In general we need to concern ourselves with the combination of both discretization error and floating point error.\n",
    "\n",
    "### Reminder:\n",
    "* **Discretization error**:  Errors arising from approximation of a function, truncation of a series...\n",
    "\n",
    "$$\\sin x \\approx x - \\frac{x^3}{3!} + \\frac{x^5}{5!} + O(x^7)$$\n",
    "\n",
    "* **Floating-point Error**:  Errors arising from approximating real numbers with finite-precision numbers\n",
    "\n",
    "$$\\pi \\approx 3.14$$\n",
    "\n",
    "or $\\frac{1}{3} \\approx 0.333333333$ in decimal cannot be represented exactly as a binary floating point number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example 1\n",
    "\n",
    "Consider a finite difference approximation to the first derivative of a function\n",
    "\n",
    "$$f^\\prime(x) \\approx \\frac{f(x + \\Delta x) - f(x)}{\\Delta x}$$\n",
    "\n",
    "**Note**: in the limit $\\Delta x\\rightarrow 0$,  this is the standard definition of the first derivative.  However we're interested in the error for a *finite* $\\Delta x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Moreover, (as we will see in future notebooks), there are many ways to approximate the first derivative.  For example we can write the \"centered first derivative\" as\n",
    "\n",
    "$$f^\\prime(x) \\approx \\frac{f(x + \\Delta x) - f(x - \\Delta x)}{2\\Delta x}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here we will just compare the the error for the two different finite-difference formulas given\n",
    "\n",
    "$$f(x) = f^\\prime(x) = e^x$$ \n",
    "\n",
    "at $x=1$ for decreasing values of $\\Delta x$. We will also introduce the idea of an 'inline' or `lambda` function in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: numpy.exp(x)\n",
    "f_prime = lambda x: numpy.exp(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "delta_x = numpy.array([2.0**(-n) for n in range(1, 60)])\n",
    "x = 1.0\n",
    "\n",
    "# Forward finite difference approximation to first derivative\n",
    "f_hat_1 = (f(x + delta_x) - f(x)) / (delta_x)\n",
    "# Centered finite difference approximation to first derivative\n",
    "f_hat_2 = (f(x + delta_x) - f(x - delta_x)) / (2.0 * delta_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.loglog(delta_x, numpy.abs(f_hat_1 - f_prime(x)), 'o-', label=\"One-Sided\")\n",
    "axes.loglog(delta_x, numpy.abs(f_hat_2 - f_prime(x)), 's-', label=\"Centered\")\n",
    "axes.legend(loc=3,fontsize=14)\n",
    "axes.set_xlabel(\"$\\Delta x$\",fontsize=16)\n",
    "axes.set_ylabel(\"Absolute Error\",fontsize=16)\n",
    "axes.set_title(\"Finite Difference approximations to $df/dx$\",fontsize=18)\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example 2\n",
    "\n",
    "Evaluate $e^x$ with its Taylor series.\n",
    "\n",
    "$$e^x = \\sum^\\infty_{n=0} \\frac{x^n}{n!}$$\n",
    "\n",
    "Can we pick $N < \\infty$ that can approximate $e^x$ over a give range $x \\in [a,b]$ such that the relative error $E$ satisfies $E < 8 \\cdot \\varepsilon_{\\text{machine}}$?\n",
    "\n",
    "We can try simply evaluating the Taylor polynomial directly for various $N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.special import factorial\n",
    "\n",
    "def my_exp(x, N=10):\n",
    "    value = 0.0\n",
    "    for n in range(N + 1):\n",
    "        value += x**n / float(factorial(n))\n",
    "        \n",
    "    return value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And test this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "eps = numpy.finfo(numpy.float64).eps\n",
    "\n",
    "x = numpy.linspace(-2., 90., 100, dtype=numpy.float64)\n",
    "MAX_N = 300\n",
    "for N in range(1, MAX_N + 1):\n",
    "    rel_error = numpy.abs((numpy.exp(x) - my_exp(x, N=N)) / numpy.exp(x))\n",
    "    if numpy.all(rel_error < 8.0 * eps):\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, rel_error/eps)\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"Relative Error/eps\")\n",
    "axes.set_title('N = {} terms'.format(N))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Can we do better?  \n",
    "\n",
    "Note: \n",
    "\n",
    "the largest value of $x$ such that $e^x <$ MAX is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.log(numpy.finfo(float).max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "and `numpy.exp` handles that just fine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(numpy.exp(709,dtype=numpy.float64))\n",
    "print(numpy.exp(-709,dtype=numpy.float64))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Your homework:  the great Exp Challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### One final example (optional):  How to calculate  Relative Error\n",
    "\n",
    "Say we wanted to compute the relative error between two values $x$ and $y$ using $x$ as the normalizing value.  Algebraically the forms\n",
    "$$\n",
    "    E = \\frac{x - y}{x}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "    E = 1 - \\frac{y}{x}\n",
    "$$\n",
    "are equivalent.  In finite precision what form might be expected to be more accurate and why?\n",
    "\n",
    "*Example based on a [blog](https://nickhigham.wordpress.com/2017/08/14/how-and-how-not-to-compute-a-relative-error/) post by Nick Higham*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Using this model the original definition contains two floating point operations such that\n",
    "$$\\begin{aligned}\n",
    "    E_1 = \\text{fl}\\left(\\frac{x - y}{x}\\right) &= \\text{fl}(\\text{fl}(x - y) / x) \\\\\n",
    "    &= \\left[ \\frac{(x - y) (1 + \\delta_1)}{x} \\right ] (1 + \\delta_2) \\\\\n",
    "    &= \\frac{x - y}{x}  (1 + \\delta_1) (1 + \\delta_2)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For the other formulation we have\n",
    "$$\\begin{aligned}\n",
    "    E_2 = \\text{fl}\\left( 1 - \\frac{y}{x} \\right ) &= \\text{fl}\\left(1 - \\text{fl}\\left(\\frac{y}{x}\\right) \\right) \\\\\n",
    "    &= \\left(1 - \\frac{y}{x} (1 + \\delta_1) \\right) (1 + \\delta_2)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we assume that all $\\text{op}$s have similar error magnitudes then we can simplify things by letting\n",
    "$$\n",
    "    |\\delta_\\ast| \\le \\epsilon.\n",
    "$$\n",
    "\n",
    "To compare the two formulations we again use the relative error between the true relative error $e_i$ and our computed versions $E_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Original definition:\n",
    "$$\\begin{aligned}\n",
    "    \\frac{e - E_1}{e} &= \\frac{\\frac{x - y}{x} - \\frac{x - y}{x}  (1 + \\delta_1) (1 + \\delta_2)}{\\frac{x - y}{x}} \\\\\n",
    "    &\\le 1 - (1 + \\epsilon) (1 + \\epsilon) = 2 \\epsilon + \\epsilon^2\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Manipulated definition:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{e - E_2}{e} &= \\frac{e - \\left[1 - \\frac{y}{x}(1 + \\delta_1) \\right] (1 + \\delta_2)}{e} \\\\\n",
    "    &= \\frac{e - \\left[e - \\frac{y}{x} \\delta_1 \\right] (1 + \\delta_2)}{e} \\\\\n",
    "    &= \\frac{e - \\left[e + e\\delta_2 - \\frac{y}{x} \\delta_1 - \\frac{y}{x} \\delta_1 \\delta_2)) \\right] }{e} \\\\\n",
    "    &= - \\delta_2 + \\frac{1}{e} \\frac{y}{x} \\left(\\delta_1 + \\delta_1 \\delta_2 \\right) \\\\\n",
    "    &= - \\delta_2 + \\frac{1 -e}{e} \\left(\\delta_1 + \\delta_1 \\delta_2 \\right) \\\\\n",
    "    &\\le \\epsilon + \\left |\\frac{1 - e}{e}\\right | (\\epsilon + \\epsilon^2)\n",
    "\\end{aligned}$$\n",
    "\n",
    "We see then that our floating point error will be dependent on the relative magnitude of $e$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Comparison of Relative Errors of estimates of Relative Error ;^)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Based on the code by Nick Higham\n",
    "# https://gist.github.com/higham/6f2ce1cdde0aae83697bca8577d22a6e\n",
    "# Compares relative error formulations using single precision and compared to double precision\n",
    "\n",
    "N = 501    # Note: Use 501 instead of 500 to avoid the zero value\n",
    "d = numpy.finfo(numpy.float32).eps * 1e4\n",
    "a = 3.0\n",
    "x = a * numpy.ones(N, dtype=numpy.float32)\n",
    "y = [x[i] + numpy.multiply((i - numpy.divide(N, 2.0, dtype=numpy.float32)), d, dtype=numpy.float32) for i in range(N)]\n",
    "\n",
    "# Compute errors and \"true\" error\n",
    "relative_error = numpy.empty((2, N), dtype=numpy.float32)\n",
    "relative_error[0, :] = numpy.abs(x - y) / x\n",
    "relative_error[1, :] = numpy.abs(1.0 - y / x)\n",
    "exact = numpy.abs( (numpy.float64(x) - numpy.float64(y)) / numpy.float64(x))\n",
    "\n",
    "# Compute differences between error calculations\n",
    "error = numpy.empty((2, N))\n",
    "for i in range(2):\n",
    "    error[i, :] = numpy.abs((relative_error[i, :] - exact) / numpy.abs(exact))\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(y, error[0, :], '.', markersize=10, label=\"$|x-y|/|x|$\")\n",
    "axes.semilogy(y, error[1, :], '.', markersize=10, label=\"$|1-y/x|$\")\n",
    "\n",
    "axes.grid(True)\n",
    "axes.set_xlabel(\"y\")\n",
    "axes.set_ylabel(\"Relative Error\")\n",
    "axes.set_xlim((numpy.min(y), numpy.max(y)))\n",
    "axes.set_ylim((5e-9, numpy.max(error[1, :])))\n",
    "axes.set_title(\"Relative Error Comparison\")\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Some other links that might be helpful regarding IEEE Floating Point:\n",
    " - [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)\n",
    " - [IEEE 754 Floating Point Calculator](http://babbage.cs.qc.edu/courses/cs341/IEEE-754.html)\n",
    " - [Numerical Computing with IEEE Floating Point Arithmetic](http://epubs.siam.org/doi/book/10.1137/1.9780898718072)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Operation Counting\n",
    "\n",
    "Discretization Error:  **Why not use more terms in the Taylor series?**\n",
    "\n",
    "Floating Point Error: **Why not use the highest precision possible?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Example 1: Matrix-Vector Multiplication\n",
    "\n",
    "Let $A, B \\in \\mathbb{R}^{N \\times N}$ and $x \\in \\mathbb{R}^N$.  \n",
    "\n",
    "1. Count the approximate number of operations it will take to compute $A x$.\n",
    "2. Do the same for $A B$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Matrix-vector product:  Defining $[A]_i$ as the $i$th row of $A$ and $A_{ij}$ as the $i$, $j$th entry then\n",
    "$$\n",
    "    A x = \\sum^N_{i=1} [A]_i \\cdot x = \\sum^N_{i=1} \\sum^N_{j=1} A_{ij} x_j\n",
    "$$\n",
    "\n",
    "Take an explicit case, say $N = 3$, then the operation count is\n",
    "$$\n",
    "    A x = [A]_1 \\cdot v + [A]_2 \\cdot v + [A]_3 \\cdot v = \\begin{bmatrix}\n",
    "        A_{11} \\times v_1 + A_{12} \\times v_2 + A_{13} \\times v_3 \\\\\n",
    "        A_{21} \\times v_1 + A_{22} \\times v_2 + A_{23} \\times v_3 \\\\\n",
    "        A_{31} \\times v_1 + A_{32} \\times v_2 + A_{33} \\times v_3\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This leads to 15 operations (6 additions and 9 multiplications).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Take another case, say $N = 4$, then the operation count is\n",
    "$$\n",
    "    A x = [A]_1 \\cdot v + [A]_2 \\cdot v + [A]_3 \\cdot v = \\begin{bmatrix}\n",
    "        A_{11} \\times v_1 + A_{12} \\times v_2 + A_{13} \\times v_3 + A_{14} \\times v_4 \\\\\n",
    "        A_{21} \\times v_1 + A_{22} \\times v_2 + A_{23} \\times v_3 + A_{24} \\times v_4 \\\\\n",
    "        A_{31} \\times v_1 + A_{32} \\times v_2 + A_{33} \\times v_3 + A_{34} \\times v_4 \\\\\n",
    "        A_{41} \\times v_1 + A_{42} \\times v_2 + A_{43} \\times v_3 + A_{44} \\times v_4 \\\\\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "This leads to 28 operations (12 additions and 16 multiplications).\n",
    "\n",
    "Generalizing this there are $N^2$ multiplications and $N (N -1)$ additions for a total of\n",
    "\n",
    "$$\n",
    "    \\text{operations} = N (N - 1) + N^2 = \\mathcal{O}(N^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Matrix-Matrix product ($AB$):  Defining $[B]_j$ as the $j$th column of $B$ then\n",
    "$$\n",
    "    (A B)_{ij} = \\sum^N_{i=1} \\sum^N_{j=1} [A]_i \\cdot [B]_j\n",
    "$$\n",
    "The inner product of two vectors is represented by\n",
    "$$\n",
    "    a \\cdot b = \\sum^N_{i=1} a_i b_i\n",
    "$$\n",
    "leading to $\\mathcal{O}(3N)$ operations.  Since there are $N^2$ entries in the resulting matrix then we would have $\\mathcal{O}(N^3)$ operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are methods for performing matrix-matrix multiplication faster.  In the following figure we see a collection of algorithms over time that have been able to bound the number of operations in certain circumstances.  Here\n",
    "$$\n",
    "    \\mathcal{O}(N^\\omega)\n",
    "$$\n",
    "![matrix multiplication operation bound](./images/bound_matrix_multiply.png)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}