05_root_finding_optimization.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": {},
   "outputs": [],
   "source": []
  },
  {
   "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\n",
    "import warnings\n",
    "import sympy\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Root Finding and Optimization\n",
    "\n",
    "Our goal in this section is to develop techniques to approximate the roots of a given function $f(x)$.  That is find solutions $x$ such that $f(x)=0$. At first glance this may not seem like a meaningful exercise, however, this problem arises  in a wide variety of circumstances. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For example, suppose that you are trying to find a solution to the equation\n",
    "\n",
    "$$\n",
    "    x^2 + x = \\alpha\\sin{x}.\n",
    "$$\n",
    "\n",
    "where $\\alpha$ is a real parameter.  Simply rearranging, the expression can be rewritten in the form\n",
    "\n",
    "$$\n",
    "    f(x) =  x^2 + x -\\alpha\\sin{x} = 0.\n",
    "$$\n",
    "\n",
    "Determining the roots of the function $f(x)$ is now equivalent to determining the solution to the original expression. Unfortunately, a number of other issues arise. In particular,  with non-linear equations, there may be multiple solutions, or no real solutions at all."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The task of approximating the roots of a function can be a deceptively difficult thing to do. For much of the treatment here we will ignore many details such as existence and uniqueness, but you should keep in mind that they are important considerations. \n",
    "\n",
    "**GOAL:** \n",
    "For this section we will focus on multiple techniques for efficiently and accurately solving the fundamental problem $f(x)=0$ for functions of a single variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Objectives\n",
    "\n",
    "* Understand the general rootfinding problem as $f(x)=0$\n",
    "* Understand the equivalent formulation as a fixed point problem $x = g(x)$\n",
    "* Understand fixed point iteration and its stability analysis\n",
    "* Understand definitions of convergence and order of convergence\n",
    "* Understand practical rootfinding algorithms and their convergence\n",
    "    * Bisection\n",
    "    * Newton's method\n",
    "    * Secant method\n",
    "    * Hybrid methods and scipy.optimize routines (root_scalar) \n",
    "* Understand basic Optimization routines\n",
    "    * Parabolic Interpolation\n",
    "    * Golden Section Search\n",
    "    * scipy.optimize routines (minimize_scalar and minimize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Future Time Annuity\n",
    "\n",
    "Can I ever retire?\n",
    "\n",
    "$$ A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] $$\n",
    "\n",
    "* $A$ total value after $n$ years\n",
    "* $P$ is payment amount per compounding period\n",
    "* $m$ number of compounding periods per year\n",
    "* $r$ annual interest rate\n",
    "* $n$ number of years to retirement\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Question:\n",
    "\n",
    "For a fix monthly Payment $P$, what does the minimum interest rate $r$ need to be so I can retire in 20 years with \\$1M. \n",
    "\n",
    "Set $P = \\frac{\\$18,000}{12} = \\$1500, \\quad m=12, \\quad n=20$.\n",
    "\n",
    "$$ \n",
    "    A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def total_value(P, m, r, n):\n",
    "    \"\"\"Total value of portfolio given parameters\n",
    "    \n",
    "    Based on following formula:\n",
    "    \n",
    "    A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n}\n",
    "                - 1 \\right ] \n",
    "    \n",
    "    :Input:\n",
    "     - *P* (float) - Payment amount per compounding period\n",
    "     - *m* (int) - number of compounding periods per year\n",
    "     - *r* (float) - annual interest rate\n",
    "     - *n* (float) - number of years to retirement\n",
    "     \n",
    "     :Returns:\n",
    "     (float) - total value of portfolio\n",
    "     \n",
    "    \"\"\"\n",
    "    \n",
    "    return P / (r / float(m)) * ( (1.0 + r / float(m))**(float(m) * n)\n",
    "                                 - 1.0)\n",
    "\n",
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "    \n",
    "r = numpy.linspace(0.05, 0.15, 100)\n",
    "goal = 1e6\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, total_value(P, m, r, 10),label='10 years',linewidth=2)\n",
    "axes.plot(r, total_value(P, m, r, 15),label='15 years',linewidth=2)\n",
    "axes.plot(r, total_value(P, m, r, n),label='20 years',linewidth=2)\n",
    "axes.plot(r, numpy.ones(r.shape) * goal, 'r--')\n",
    "axes.set_xlabel(\"r (interest rate)\", fontsize=16)\n",
    "axes.set_ylabel(\"A (total value)\", fontsize=16)\n",
    "axes.set_title(\"When can I retire?\",fontsize=18)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.set_xlim((r.min(), r.max()))\n",
    "axes.set_ylim((total_value(P, m, r.min(), 10), total_value(P, m, r.max(), n)))\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fixed Point Iteration\n",
    "\n",
    "How do we go about solving this?\n",
    "\n",
    "Could try to solve at least partially for $r$:\n",
    "\n",
    "$$ A = \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] ~~~~ \\Rightarrow ~~~~~$$\n",
    "\n",
    "$$ r = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ] ~~~~ \\Rightarrow ~~~~~$$\n",
    "\n",
    "$$ r = g(r)$$\n",
    "or \n",
    "$$ g(r) - r = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Plot these\n",
    "\n",
    "$$ r = g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def g(P, m, r, n, A):\n",
    "    \"\"\"Reformulated minimization problem\n",
    "    \n",
    "    Based on following formula:\n",
    "    \n",
    "    g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "    \n",
    "    :Input:\n",
    "     - *P* (float) - Payment amount per compounding period\n",
    "     - *m* (int) - number of compounding periods per year\n",
    "     - *r* (float) - annual interest rate\n",
    "     - *n* (float) - number of years to retirement\n",
    "     - *A* (float) - total value after $n$ years\n",
    "     \n",
    "     :Returns:\n",
    "     (float) - value of g(r)\n",
    "     \n",
    "    \"\"\"\n",
    "    \n",
    "    return P * m / A * ( (1.0 + r / float(m))**(float(m) * n)\n",
    "                                 - 1.0)\n",
    "\n",
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "    \n",
    "r = numpy.linspace(0.00, 0.1, 100)\n",
    "goal = 1e6\n",
    "\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(r, g(P, m, r, n, goal),label='$g(r)$')\n",
    "axes.plot(r, r, 'r--',label='$r$')\n",
    "axes.set_xlabel(\"r (interest rate)\",fontsize=16)\n",
    "axes.set_ylabel(\"$g(r)$\",fontsize=16)\n",
    "axes.set_title(\"Minimum rate for a 20 year retirement?\",fontsize=18)\n",
    "axes.set_ylim([0, 0.12])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.set_xlim((0.00, 0.1))\n",
    "axes.set_ylim((g(P, m, 0.00, n, goal), g(P, m, 0.1, n, goal)))\n",
    "axes.legend(fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(r, g(P, m, r, n, goal)-r,label='$g(r) - r$')\n",
    "axes.plot(r, numpy.zeros(r.shape), 'r--',label='$0$')\n",
    "axes.set_xlabel(\"r (interest rate)\",fontsize=16)\n",
    "axes.set_ylabel(\"residual\",fontsize=16)\n",
    "axes.set_title(\"Minimum rate for a 20 year retirement?\",fontsize=18)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.set_xlim((0.00, 0.1))\n",
    "axes.legend(fontsize=14)\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question:\n",
    "\n",
    "A single root $r>0$ clearly exists around $r=0.088$.  But how to find it?\n",
    "\n",
    "One option might be to take a guess say $r_0 = 0.088$ and form the iterative scheme\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    r_1 &= g(r_0)\\\\\n",
    "    r_2 &= g(r_1)\\\\\n",
    "        &\\vdots \\\\\n",
    "    r_{k} &= g(r_{k-1})\\\\\n",
    "\\end{align}\n",
    "$$ \n",
    "and hope this converges as $k\\rightarrow\\infty$ (or faster) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Easy enough to code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "r = 0.088\n",
    "K = 20\n",
    "for k in range(K):\n",
    "    print(r)\n",
    "    r = g(P,m,r,n,goal)\n",
    "     \n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 2:\n",
    "\n",
    "Let $f(x) = x - e^{-x}$, solve $f(x) = 0$\n",
    "\n",
    "Equivalent to $x = e^{-x}$ or $x = g(x)$ where $g(x) = e^{-x}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.2, 1.0, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, numpy.exp(-x), 'r',label='g(x)=exp(-x)$')\n",
    "axes.plot(x, x, label='$x$')\n",
    "axes.set_xlabel(\"$x$\",fontsize=16)\n",
    "axes.legend()\n",
    "plt.grid()\n",
    "\n",
    "f = lambda x : x - numpy.exp(-x)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, f(x),label='$f(x) = x - g(x)$')\n",
    "axes.plot(x, numpy.zeros(x.shape), 'r--',label='$0$')\n",
    "axes.set_xlabel(\"$x$\",fontsize=16)\n",
    "axes.set_ylabel(\"residual\",fontsize=16)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.legend(fontsize=14)\n",
    "axes.grid()\n",
    "plt.show()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Again, consider the iterative scheme\n",
    "\n",
    "set $x_0$ then compute\n",
    "\n",
    "$$\n",
    "    x_k = g(x_{k-1})\\quad \\mathrm{for}\\quad k=1,2,3\\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or again in code\n",
    "```python\n",
    "x = x0\n",
    "for i in range(N):\n",
    "    x = g(x)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.2, 1.0, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.exp(-x), 'r',label='$g(x)=exp(-x)$')\n",
    "axes.plot(x, x, 'b',label='$x$')\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_aspect('equal')\n",
    "axes.legend(fontsize=14)\n",
    "\n",
    "x = 0.4\n",
    "print('\\tx\\t      exp(-x)\\t    residual')\n",
    "for steps in range(6):\n",
    "    residual = numpy.abs(numpy.exp(-x) - x)\n",
    "    print(\"{:12.7f}\\t{:12.7f}\\t{:12.7f}\".format(x, numpy.exp(-x), residual))\n",
    "    axes.plot(x, numpy.exp(-x),'kx')\n",
    "    axes.text(x+0.01, numpy.exp(-x)+0.01, steps, fontsize=\"15\")\n",
    "    x = numpy.exp(-x)\n",
    "\n",
    "\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example 3:\n",
    "\n",
    "Let $f(x) = \\ln x + x$ and solve $f(x) = 0$ or $x = -\\ln x$.\n",
    "\n",
    "Note that this problem is equivalent to $x = e^{-x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.1, 1.0, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, -numpy.log(x), 'r',label='$g(x)=-\\log(x)$')\n",
    "axes.plot(x, x, 'b',label='$x$')\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"f(x)\",fontsize=16)\n",
    "axes.set_ylim([0, 1.5])\n",
    "axes.set_aspect('equal')\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "\n",
    "x = 0.55\n",
    "print('\\tx\\t      -log(x)\\t    residual')\n",
    "\n",
    "for steps in range(5):\n",
    "    residual = numpy.abs(numpy.log(x) + x)\n",
    "    print(\"{:12.7f}\\t{:12.7f}\\t{:12.7f}\".format(x, -numpy.log(x), residual))\n",
    "    axes.plot(x, -numpy.log(x),'kx')\n",
    "    axes.text(x + 0.01, -numpy.log(x) + 0.01, steps, fontsize=\"15\")\n",
    "    x = -numpy.log(x)\n",
    "\n",
    "plt.grid()    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### These are equivalent problems!  \n",
    "\n",
    "Something is awry..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Analysis of Fixed Point Iteration\n",
    "\n",
    "Existence and uniqueness of fixed point problems\n",
    "\n",
    "*Existence:*\n",
    "\n",
    "Assume $g \\in C[a, b]$, if the range of the mapping $y = g(x)$ satisfies $y \\in [a, b] ~~ \\forall ~~  x \\in [a, b]$ then $g$ has a fixed point in $[a, b]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.0, 1.0, 100)\n",
    "\n",
    "# Plot function and intercept\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.exp(-x), 'r',label='$g(x)$')\n",
    "axes.plot(x, x, 'b',label='$x$')\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.set_title('$g(x) = e^{-x}$',fontsize=24)\n",
    "\n",
    "# Plot domain and range\n",
    "axes.plot(numpy.ones(x.shape) * 0.4, x, '--k')\n",
    "axes.plot(numpy.ones(x.shape) * 0.8, x, '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * numpy.exp(-0.4), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * numpy.exp(-0.8), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * 0.4, '--',color='gray',linewidth=.5)\n",
    "axes.plot(x, numpy.ones(x.shape) * 0.8, '--',color='gray',linewidth=.5)\n",
    "\n",
    "\n",
    "axes.set_xlim((0.0, 1.0))\n",
    "axes.set_ylim((0.0, 1.0))\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.1, 1.0, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, -numpy.log(x), 'r',label='$g(x)$')\n",
    "axes.plot(x, x, 'b',label='$x$')\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_xlim([0.1, 1.0])\n",
    "axes.set_ylim([0.1, 1.0])\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.set_title('$g(x) = -\\ln(x)$',fontsize=24)\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "\n",
    "# Plot domain and range\n",
    "axes.plot(numpy.ones(x.shape) * 0.4, x, '--k')\n",
    "axes.plot(numpy.ones(x.shape) * 0.8, x, '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * -numpy.log(0.4), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * -numpy.log(0.8), '--k')\n",
    "axes.plot(x, numpy.ones(x.shape) * 0.4, '--',color='gray',linewidth=.5)\n",
    "axes.plot(x, numpy.ones(x.shape) * 0.8, '--',color='gray',linewidth=.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Uniqueness:*\n",
    "\n",
    "Additionally, suppose $g'(x)$ is defined on $x \\in [a, b]$ and $\\exists K < 1$ such that\n",
    "\n",
    "$$\n",
    "    |g'(x)| \\leq K < 1 \\quad \\forall \\quad x \\in (a,b)\n",
    "$$\n",
    "\n",
    "then $g$ has a unique fixed point $P \\in [a,b]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0.4, 0.8, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, numpy.abs(-numpy.exp(-x)), 'r')\n",
    "axes.plot(x, numpy.ones(x.shape), 'k--')\n",
    "axes.set_xlabel(\"$x$\",fontsize=18)\n",
    "axes.set_ylabel(\"$|g\\,'(x)|$\",fontsize=18)\n",
    "axes.set_ylim((0.0, 1.1))\n",
    "axes.set_title(\"$g(x) = e^{-x}$\",fontsize=20)\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Asymptotic convergence*: Behavior of fixed point iterations\n",
    "\n",
    "$$x_{k+1} = g(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Assume that a fixed point $x^\\ast$ exists, such that  \n",
    "$$\n",
    "x^\\ast = g(x^\\ast)\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Then define \n",
    "$$\n",
    "    x_{k+1} = x^\\ast + e_{k+1} \\quad \\quad  x_k = x^\\ast + e_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "substituting\n",
    "$$\n",
    "    x^\\ast + e_{k+1} = g(x^\\ast + e_k)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Evaluate $$\n",
    "    g(x^\\ast + e_k)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Taylor expand $g(x)$ about $x^\\ast$ and substitute $$x = x_k = x^\\ast + e_k$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "    g(x^\\ast + e_k) = g(x^\\ast) + g'(x^\\ast) e_k + \\frac{g''(x^\\ast) e_k^2}{2} + O(e_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "from our definition  $$x^\\ast + e_{k+1} = g(x^\\ast + e_k)$$ we have\n",
    "$$\n",
    "    x^\\ast + e_{k+1} = g(x^\\ast) + g'(x^\\ast) e_k + \\frac{g''(x^\\ast) e_k^2}{2} +  O(e_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that because $x^* = g(x^*)$ these terms cancel leaving\n",
    "\n",
    "$$e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2}$$\n",
    "\n",
    "So if $|g'(x^*)| \\leq K < 1$ we can conclude that\n",
    "\n",
    "$$|e_{k+1}| = K |e_k|$$\n",
    "\n",
    "which shows convergence.  Also note that $K$ is related to $|g'(x^*)|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Convergence of iterative schemes\n",
    "\n",
    "Given any iterative scheme where\n",
    "\n",
    "$$|e_{k+1}| = C |e_k|^n$$\n",
    "\n",
    "If $C < 1$ and:\n",
    " - $n=1$ then the scheme is **linearly convergent**\n",
    " - $n=2$ then the scheme is **quadratically convergent**\n",
    " - $n > 1$ the scheme can also be called **superlinearly convergent**\n",
    "\n",
    "If $C > 1$ then the scheme is **divergent**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Examples Revisited\n",
    "\n",
    "* Example 1:\n",
    "\n",
    "$$\n",
    "g(x) = e^{-x}\\quad\\mathrm{with}\\quad x^* \\approx 0.56\n",
    "$$\n",
    " \n",
    "   $$|g'(x^*)| = |-e^{-x^*}| \\approx 0.56$$\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Example 2: \n",
    "\n",
    "$$g(x) = - \\ln x \\quad \\text{with} \\quad x^* \\approx 0.56$$\n",
    "\n",
    "   $$|g'(x^*)| = \\frac{1}{|x^*|} \\approx 1.79$$\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Example 3: The retirement problem\n",
    "\n",
    "$$\n",
    "    r = g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r, P, m, A, n = sympy.symbols('r P m A n')\n",
    "g_sym = P * m / A * ((1 + r /m)**(m * n) - 1)\n",
    "g_sym"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "g_prime = g_sym.diff(r)\n",
    "g_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r_star = 0.08985602484084668\n",
    "print(\"g'(r*) = \", g_prime.subs({P: 1500.0, m: 12, n:20, A: 1e6, r: r_star})) \n",
    "print(\"g(r*) - r* = {}\".format(g_sym.subs({P: 1500.0, m: 12, n:20, A: 1e6, r: r_star}) - r_star))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Example 3: The retirement problem\n",
    "\n",
    "$$\n",
    "    r = g(r) = \\frac{P \\cdot m}{A} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = sympy.lambdify(r, g_prime.subs({P: 1500.0, m: 12, n:20, A: 1e6}))\n",
    "g = sympy.lambdify(r, g_sym.subs({P: 1500.0, m: 12, n:20, A: 1e6}))\n",
    "r = numpy.linspace(-0.01, 0.1, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(7,5))\n",
    "fig.set_figwidth(2. * fig.get_figwidth())\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(r, g(r),label='$g(r)$')\n",
    "axes.plot(r, r, 'r--',label='$r$')\n",
    "axes.set_xlabel(\"r (interest rate)\",fontsize=14)\n",
    "axes.set_ylabel(\"$g(r)$\",fontsize=14)\n",
    "axes.set_title(\"Minimum rate for a 20 year retirement?\",fontsize=14)\n",
    "axes.set_ylim([0, 0.12])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.set_xlim((0.00, 0.1))\n",
    "axes.set_ylim(g(0.00), g(0.1))\n",
    "axes.legend()\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(r, f(r))\n",
    "axes.plot(r, numpy.ones(r.shape), 'k--')\n",
    "axes.plot(r_star, f(r_star), 'ro')\n",
    "axes.plot(0.0, f(0.0), 'ro')\n",
    "axes.set_xlim((-0.01, 0.1))\n",
    "axes.set_xlabel(\"$r$\",fontsize=14)\n",
    "axes.set_ylabel(\"$g'(r)$\",fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Better ways for root-finding/optimization\n",
    "\n",
    "If $x^*$ is a fixed point of $g(x)$ then $x^*$ is also a *root* of $f(x^*) = g(x^*) - x^*$ s.t. $f(x^*) = 0$.\n",
    "\n",
    "For instance:\n",
    "\n",
    "$$f(r) = r - \\frac{m P}{A} \\left [ \\left (1 + \\frac{r}{m} \\right)^{m n} - 1 \\right ] =0 $$\n",
    "\n",
    "or\n",
    "\n",
    "$$f(r) = A - \\frac{m P}{r} \\left [ \\left (1 + \\frac{r}{m} \\right)^{m n} - 1 \\right ] =0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Classical Methods\n",
    " - Bisection (linear convergence)\n",
    " - Newton's Method (quadratic convergence)\n",
    " - Secant Method (super-linear)\n",
    " \n",
    "## Combined Methods\n",
    " - RootSafe (Newton + Bisection)\n",
    " - Brent's Method (Secant + Bisection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bracketing and Bisection\n",
    "\n",
    "A **bracket** is an interval $[a,b]$ that contains at least one zero or minima/maxima of interest.  \n",
    "\n",
    "In the case of a zero the bracket should satisfy \n",
    "$$\n",
    "    \\text{sign}(f(a)) \\neq \\text{sign}(f(b)).\n",
    "$$\n",
    "\n",
    "In the case of minima or maxima we need \n",
    "$$\n",
    "    \\text{sign}(f'(a)) \\neq \\text{sign}(f'(b))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem**:  \n",
    "\n",
    "Let\n",
    "$$\n",
    "    f(x) \\in C[a,b] \\quad \\text{and} \\quad \\text{sign}(f(a)) \\neq \\text{sign}(f(b))\n",
    "$$\n",
    "\n",
    "then there exists a number \n",
    "$$\n",
    "    c \\in (a,b) \\quad \\text{s.t.} \\quad f(c) = 0.\n",
    "$$\n",
    "(proof uses intermediate value theorem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Example**:  The retirement problem again.  For fixed $A, P, m, n$\n",
    "\n",
    "$$\n",
    "    f(r) = A - \\frac{P}{(r / m)} \\left[ \\left(1 + \\frac{r}{m} \\right)^{m \\cdot n} - 1 \\right ]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.1, 100)\n",
    "f = lambda r, A, m, P, n: A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r, A, m, P, n), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "axes.set_xlabel(\"r\", fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\", fontsize=16)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.grid()\n",
    "\n",
    "a = 0.075\n",
    "b = 0.095\n",
    "axes.plot(a, f(a, A, m, P, n), 'ko')\n",
    "axes.plot([a, a], [0.0, f(a, A, m, P, n)], 'k--')\n",
    "axes.plot(b, f(b, A, m, P, n), 'ko')\n",
    "axes.plot([b, b], [f(b, A, m, P, n), 0.0], 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Basic bracketing algorithms shrink the bracket while ensuring that the root/extrema remains within the bracket.\n",
    "\n",
    "What ways could we \"shrink\" the bracket so that the end points converge to the root/extrema?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Bisection Algorithm\n",
    "\n",
    "Given a bracket $[a,b]$ and a function $f(x)$ - \n",
    "1. Initialize with bracket\n",
    "2. Iterate\n",
    "   1. Cut bracket in half and check to see where the zero is\n",
    "   2. Set bracket to new bracket based on what direction we went"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "##### basic  code\n",
    "```python\n",
    "def bisection(f,a,b,tol):\n",
    "    c = (a + b)/2.  \n",
    "    f_a = f(a)\n",
    "    f_b = f(b)\n",
    "    f_c = f(c)\n",
    "    \n",
    "    for step in range(1, MAX_STEPS + 1):\n",
    "        if numpy.abs(f_c) < tol:\n",
    "            break\n",
    "        if numpy.sign(f_a) != numpy.sign(f_c):\n",
    "            b = c\n",
    "            f_b = f_c\n",
    "        else:\n",
    "            a = c\n",
    "            f_a = f_c\n",
    "        c = (a + b)/ 2.0\n",
    "        f_c = f(c)\n",
    "    return c\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Some real code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# real code with standard bells and whistles\n",
    "\n",
    "def bisection(f,a,b,tol = 1.e-6):\n",
    "    \"\"\" uses bisection to isolate a root x of a function of a single variable f such that f(x) = 0.\n",
    "        the root must exist within an initial bracket a < x < b\n",
    "        returns when f(x) at the midpoint of the bracket < tol\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    \n",
    "    f: function of a single variable f(x) of type float\n",
    "    a: float\n",
    "        left bracket a < x\n",
    "    b: float\n",
    "        right bracket x < b\n",
    "        Note:  the signs of f(a) and f(b) must be different to insure a bracket\n",
    "    tol: float\n",
    "        tolerance.  Returns when |f((a+b)/2)| < tol\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    x: float\n",
    "        midpoint of final bracket\n",
    "    x_array: numpy array\n",
    "        history of bracket centers (for plotting later)\n",
    "        \n",
    "    Raises:\n",
    "    -------\n",
    "    ValueError:  \n",
    "        if initial bracket is invalid \n",
    "    Warning: \n",
    "        if number of iterations exceed MAX_STEPS\n",
    "    \n",
    "    \"\"\"\n",
    "    MAX_STEPS = 1000\n",
    "    \n",
    "    # initialize\n",
    "    c = (a + b)/2.\n",
    "    c_array = [ c ]\n",
    "    \n",
    "    f_a = f(a)\n",
    "    f_b = f(b)\n",
    "    f_c = f(c)\n",
    "    \n",
    "    # check bracket\n",
    "    if numpy.sign(f_a) == numpy.sign(f_b):\n",
    "        raise ValueError(\"no bracket: f(a) and f(b) must have different signs\")\n",
    "        \n",
    "    # Loop until we reach the TOLERANCE or we take MAX_STEPS\n",
    "    for step in range(1, MAX_STEPS + 1):\n",
    "        # Check tolerance - Could also check the size of delta_x\n",
    "        # We check this first as we have already initialized the values\n",
    "        # in c and f_c\n",
    "        if numpy.abs(f_c) < tol:\n",
    "            break\n",
    "\n",
    "        if numpy.sign(f_a) != numpy.sign(f_c):\n",
    "            b = c\n",
    "            f_b = f_c\n",
    "        else:\n",
    "            a = c\n",
    "            f_a = f_c\n",
    "        c = (a + b)/2.\n",
    "        f_c = f(c)\n",
    "        c_array.append(c)\n",
    "        \n",
    "    if step == MAX_STEPS:\n",
    "        warnings.warn('Maximum number of steps exceeded')\n",
    "    \n",
    "    return c, numpy.array(c_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# set up function as an inline lambda function\n",
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "f = lambda r: A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "# Initialize bracket\n",
    "a = 0.07\n",
    "b = 0.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# find root\n",
    "r_star, r_array = bisection(f, a, b, tol=1e-8)\n",
    "print('root at r = {}, f(r*) = {}, {} steps'.format(r_star,f(r_star),len(r_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "axes.set_xlabel(\"r\", fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\", fontsize=16)\n",
    "# axes.set_xlim([0.085, 0.091])\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.plot(a, f(a), 'ko')\n",
    "axes.plot([a, a], [0.0, f(a)], 'k--')\n",
    "axes.text(a, f(a), str(0), fontsize=\"15\")\n",
    "\n",
    "axes.plot(b, f(b), 'ko')\n",
    "axes.plot([b, b], [f(b), 0.0], 'k--')\n",
    "axes.text(b, f(b), str(1), fontsize=\"15\")\n",
    "axes.grid()\n",
    "\n",
    "# plot out the first N steps\n",
    "N = 5\n",
    "for k,r in enumerate(r_array[:N]):\n",
    "    # Plot iteration\n",
    "    axes.plot(r, f(r),'kx')\n",
    "    axes.text(r, f(r), str(k + 2), fontsize=\"15\")\n",
    "    \n",
    "axes.plot(r_star, f(r_star), 'go', markersize=10)\n",
    "axes.set_title('Bisection method:  first {} steps'.format(N), fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What is the smallest tolerance that can be achieved with this routine?  Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# find root\n",
    "r_star, r_array = bisection(f, a, b, tol=1e-9 )\n",
    "print('root at r = {}, f(r*) = {}, {} steps'.format(r_star,f(r_star),len(r_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# this might be useful\n",
    "print(numpy.diff(r_array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Convergence of Bisection\n",
    "\n",
    "Generally have\n",
    "$$\n",
    "    |e_{k+1}| = C |e_k|^n\n",
    "$$\n",
    "where we need $C < 1$ and $n > 0$.\n",
    "\n",
    "Letting $\\Delta x_k$ be the width of the $k$th bracket we can then estimate the error with\n",
    "$$\n",
    "    e_k \\approx \\Delta x_k\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "    e_{k+1} \\approx \\frac{1}{2} \\Delta x_k.\n",
    "$$\n",
    "Due to the relationship then between $x_k$ and $e_k$ we then know\n",
    "$$\n",
    "    |e_{k+1}| = \\frac{1}{2} |e_k|\n",
    "$$\n",
    "so therefore the method is linearly convergent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Newton's Method (Newton-Raphson)\n",
    " - Given a bracket, bisection is guaranteed to converge linearly to a root\n",
    " - However bisection uses almost no information about $f(x)$ beyond its sign at a point\n",
    " - Can we do \"better\"? <font color='red'>Newton's method</font>, *when well behaved* can achieve quadratic convergence.\n",
    " \n",
    "**Basic Ideas**: There are multiple interpretations we can use to derive Newton's method\n",
    "* Use Taylor's theorem to estimate a correction to minimize the residual $f(x)=0$ \n",
    "* A geometric interpretation that approximates $f(x)$ locally as a straight line to  predict where $x^*$ might be.\n",
    "* As a special case of a fixed-point iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Perhaps the simplest derivation uses Taylor series.  Consider an initial guess at point $x_k$.  For arbitrary $x_k$, it's unlikely $f(x_k)=0$.  However we can hope there is a correction $\\delta_k$ such that at\n",
    "\n",
    "$$\n",
    "x_{k+1} = x_k + \\delta_k\n",
    "$$\n",
    "and \n",
    "$$\n",
    "    f(x_{k+1}) = 0 \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "expanding in a Taylor series around point $x_k$ \n",
    "\n",
    "$$\n",
    "    f(x_k + \\delta_k) \\approx  f(x_k) + f'(x_k) \\delta_k  + O(\\delta_k^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "substituting into $f(x_{k+1})=0$ and dropping the higher order terms gives\n",
    "\n",
    "$$\n",
    "    f(x_k) + f'(x_k) \\delta_k  =0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "substituting into $f(x_{k+1})=0$ and dropping the higher order terms gives\n",
    "\n",
    "$$\n",
    "    f(x_k) + f'(x_k) \\delta_k  =0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or solving for the correction\n",
    "$$\n",
    "    \\delta_k = -f(x_k)/f'(x_k)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which leads to the update for the next iteration\n",
    "$$\n",
    "    x_{k+1} = x_k + \\delta_k \n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "    x_{k+1} = x_k -f(x_k)/f'(x_k)\n",
    "$$\n",
    "\n",
    "rinse and repeat, as it's still unlikely that $f(x_{k+1})=0$ (but we hope the error will be reduced) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Algorithm\n",
    "\n",
    "1. Initialize $x = x_0$\n",
    "1. While ( $f(x) > tol$ )\n",
    "  - solve $\\delta = -f(x)/f'(x)$\n",
    "  - update $x \\leftarrow x + \\delta$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Geometric interpretation\n",
    "\n",
    "By truncating the taylor series at first order,  we are locally approximating $f(x)$ as a straight line tangent to the point $f(x_k)$.  If the function was linear at that point, we could find its intercept such that $f(x_k+\\delta_k)=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "f_prime = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "            -P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) \\\n",
    "                + P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2\n",
    "\n",
    "# Initial guess\n",
    "x_k = 0.06\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "# Plot x_k point\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.text(x_k, -5e4, \"$x_k$\", fontsize=16)\n",
    "axes.plot(x_k, 0.0, 'xk')\n",
    "axes.text(x_k, f(x_k) + 2e4, \"$f(x_k)$\", fontsize=16)\n",
    "axes.plot(r, f_prime(x_k) * (r - x_k) + f(x_k), 'k')\n",
    "\n",
    "# Plot x_{k+1} point\n",
    "x_k = x_k - f(x_k) / f_prime(x_k)\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.text(x_k, 1e4, \"$x_{k+1}$\", fontsize=16)\n",
    "axes.plot(x_k, 0.0, 'xk')\n",
    "axes.text(0.0873, f(x_k) - 2e4, \"$f(x_{k+1})$\", fontsize=16)\n",
    "\n",
    "axes.set_xlabel(\"r\",fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\",fontsize=16)\n",
    "axes.set_title(\"Newton-Raphson Steps\",fontsize=18)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we simply approximate the derivative $f'(x_k)$ with its finite difference approximation\n",
    "\n",
    "$$\n",
    "    f'(x_k) \\approx \\frac{0 - f(x_k)}{x_{k+1} - x_k}\n",
    "$$\n",
    "\n",
    "we can rearrange to find $x_{k+1}$ as\n",
    "\n",
    "$$\n",
    "    x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}\n",
    "$$\n",
    "\n",
    "which is the classic Newton-Raphson iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Some code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def newton(f,f_prime,x0,tol = 1.e-6):\n",
    "    \"\"\" uses newton's method to find a root x of a function of a single variable f\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    f: function f(x)\n",
    "        returns type: float\n",
    "    f_prime: function f'(x)\n",
    "        returns type: float\n",
    "    x0: float\n",
    "        initial guess\n",
    "    tolerance: float\n",
    "        Returns when |f(x)| < tol\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    x: float\n",
    "        final iterate\n",
    "    x_array: numpy array\n",
    "        history of iteration points\n",
    "        \n",
    "    Raises:\n",
    "    -------\n",
    "    Warning: \n",
    "        if number of iterations exceed MAX_STEPS\n",
    "    \"\"\"\n",
    "    MAX_STEPS = 200\n",
    "    \n",
    "    x = x0\n",
    "    x_array = [ x0 ]\n",
    "    for k in range(1, MAX_STEPS + 1):\n",
    "        x = x  - f(x) / f_prime(x)\n",
    "        x_array.append(x)\n",
    "        if numpy.abs(f(x)) < tol:\n",
    "            break\n",
    "        \n",
    "    if k == MAX_STEPS:\n",
    "        warnings.warn('Maximum number of steps exceeded')\n",
    "    \n",
    "    return x, numpy.array(x_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Set the problem up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "f_prime = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "            -P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) \\\n",
    "                + P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x0 = 0.06\n",
    "x, x_array = newton(f, f_prime, x0, tol=1.e-8)\n",
    "print('x = {}, f(x) = {}, Nsteps = {}'.format(x, f(x), len(x_array)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r = numpy.linspace(0.05, 0.10, 100)\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "for n, x in enumerate(x_array):\n",
    "    axes.plot(x, f(x),'kx')\n",
    "    axes.text(x, f(x), str(n), fontsize=\"15\") \n",
    "          \n",
    "axes.set_xlabel(\"r\", fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\", fontsize=16)\n",
    "axes.set_title(\"Newton-Raphson Steps\", fontsize=18)\n",
    "axes.grid()\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.semilogy(numpy.arange(len(x_array)), numpy.abs(f(x_array)), 'bo-')\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iterations', fontsize=16)\n",
    "axes.set_ylabel('Residual $|f(r)|$', fontsize=16)\n",
    "axes.set_title('Convergence', fontsize=18)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the smallest tolerance that can be achieved with this routine?  Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: \n",
    "\n",
    "$$f(x) = x - e^{-x}$$\n",
    "\n",
    "$$f'(x) = 1 + e^{-x}$$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)} = x_k - \\frac{x_k - e^{-x_k}}{1 + e^{-x_k}}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### setup in sympy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = sympy.symbols('x') \n",
    "f = x - sympy.exp(-x)\n",
    "f_prime = f.diff(x)\n",
    "f, f_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = sympy.lambdify(x,f)\n",
    "f_prime = sympy.lambdify(x,f_prime)\n",
    "print(f(0.),f_prime(0.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "x0 = 0.\n",
    "x, x_array = newton(f, f_prime, x0, tol = 1.e-8)\n",
    "print('x = {}, f(x) = {}, Nsteps = {}'.format(x, f(x), len(x_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "xa = numpy.linspace(-1,1,100)\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.plot(xa,f(xa),'b')\n",
    "axes.plot(xa,numpy.zeros(xa.shape),'r--')\n",
    "axes.plot(x,f(x),'go', markersize=10)\n",
    "axes.plot(x0,f(x0),'kx',markersize=10)\n",
    "axes.grid()\n",
    "axes.set_xlabel('x', fontsize=16)\n",
    "axes.set_ylabel('f(x)', fontsize=16)\n",
    "axes.set_title('$f(x) = x - e^{-x}$', fontsize=18)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.semilogy(numpy.arange(len(x_array)), numpy.abs(f(x_array)), 'bo-')\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iterations', fontsize=16)\n",
    "axes.set_ylabel('Residual $|f(r)|$', fontsize=16)\n",
    "axes.set_title('Convergence', fontsize=18)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Asymptotic Convergence of Newton's Method\n",
    "\n",
    "Newton's method can be also considered a fixed point iteration\n",
    "\n",
    "$$x_{k+1} = g(x_k)$$\n",
    "\n",
    "with  $g(x) = x - \\frac{f(x)}{f'(x)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Again if $x^*$ is the fixed point and $e_k$ the error at iteration $k$:\n",
    "\n",
    "$$x_{k+1} = x^* + e_{k+1} \\quad \\quad x_k = x^* + e_k$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Taylor Expansion around $x^*$\n",
    "\n",
    "$$\n",
    "    x^* + e_{k+1} = g(x^* + e_k) = g(x^*) + g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + O(e_k^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Note that as before $x^*$ and $g(x^*)$ cancel:\n",
    "\n",
    "$$e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + \\ldots$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What about $g'(x^*)$ though?  \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    g(x) &= x - \\frac{f(x)}{f'(x)} \\\\\n",
    "    g'(x) & = 1 - \\frac{f'(x)}{f'(x)} + \\frac{f(x) f''(x)}{(f'(x))^2} = \\frac{f(x) f''(x)}{(f'(x))^2}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which evaluated at $x = x^*$ becomes\n",
    "\n",
    "$$\n",
    "    g'(x^*) = \\frac{f(x^*)f''(x^*)}{f'(x^*)^2} = 0\n",
    "$$\n",
    "\n",
    "since $f(x^\\ast) = 0$ by definition (assuming $f''(x^\\ast)$ and $f'(x^\\ast)$ are appropriately behaved)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Back to our expansion we have again\n",
    "\n",
    "$$\n",
    "    e_{k+1} = g'(x^*) e_k + \\frac{g''(x^*) e_k^2}{2!} + \\ldots\n",
    "$$\n",
    "\n",
    "which simplifies to \n",
    "\n",
    "$$\n",
    "    e_{k+1} = \\frac{g''(x^*) e_k^2}{2!} + \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which leads to \n",
    "$$\n",
    "    |e_{k+1}| < \\left | \\frac{g''(x^*)}{2!} \\right | |e_k|^2\n",
    "$$\n",
    "\n",
    "Newton's method is therefore quadratically convergent where  the constant is controlled by the second derivative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Convergence for a non-simple root\n",
    "\n",
    "Consider our first problem\n",
    "\n",
    "$$\n",
    "    f(x) = x^2 + x - \\sin(x)\n",
    "$$\n",
    "\n",
    "the case is, unfortunately, not as rosey.  Why might this be?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Setup the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x*x + x - numpy.sin(x)\n",
    "f_prime = lambda x: 2*x + 1. - numpy.cos(x)\n",
    "x0 = .9\n",
    "\n",
    "x, x_array = newton(f, f_prime, x0, tol= 1.e-16)\n",
    "print('x = {}, f(x) = {}, Nsteps = {}'.format(x, f(x), len(x_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "xa = numpy.linspace(-2,2,100)\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.plot(xa,f(xa),'b')\n",
    "axes.plot(xa,numpy.zeros(xa.shape),'r--')\n",
    "axes.plot(x,f(x),'go', markersize=10)\n",
    "axes.plot(x0,f(x0),'kx', markersize=10)\n",
    "axes.grid()\n",
    "axes.set_xlabel('x', fontsize=16)\n",
    "axes.set_ylabel('f(x)', fontsize=16)\n",
    "axes.set_title('$f(x) = x^2 +x - sin(x)$', fontsize=18)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.semilogy(numpy.arange(len(x_array)), numpy.abs(f(x_array)), 'bo-')\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iterations', fontsize=16)\n",
    "axes.set_ylabel('Residual $|f(r)|$', fontsize=16)\n",
    "axes.set_title('Convergence', fontsize=18)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence appears linear, can you show this?:\n",
    "$$f(x) = x^2 + x -\\sin (x)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: behavior of Newton with multiple roots\n",
    "$f(x) = \\sin (2 \\pi x)$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{\\sin (2 \\pi x_k)}{2 \\pi \\cos (2 \\pi x_k)}= x_k - \\frac{1}{2 \\pi} \\tan (2 \\pi x_k)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0, 2, 1000)\n",
    "f = lambda x: numpy.sin(2.0 * numpy.pi * x)\n",
    "f_prime = lambda x: 2.0 * numpy.pi * numpy.cos(2.0 * numpy.pi * x)\n",
    "x_kp = lambda x: x - f(x)/f_prime(x)\n",
    "\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, f(x),'b')\n",
    "axes.plot(x, f_prime(x), 'r')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_title(\"Comparison of $f(x)$ and $f'(x)$\")\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_xlim((0,2))\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "\n",
    "x_k = 0.3\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x, f_prime(x_k) * (x - x_k) + f(x_k), 'k')\n",
    "\n",
    "x_k = x_k - f(x_k) / f_prime(x_k)\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, f(x),'b')\n",
    "axes.plot(x, x_kp(x), 'r')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_title(\"Comparison of $f(x)$ and $x_{k+1}(x)$\",fontsize=18)\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_xlim((0,2))\n",
    "axes.plot(x, numpy.zeros(x.shape), 'k--')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basins of Attraction\n",
    "\n",
    "Given a point $x_0$ can we determine if Newton-Raphson converges and to **which root** it converges to?\n",
    "\n",
    "A *basin of attraction* $X$ for Newton's methods is defined as the set such that $\\forall x \\in X$ Newton iterations converges to the same root.  Unfortunately this is far from a trivial thing to determine and even for simple functions can lead to regions that are complicated or even fractal.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# calculate the basin of attraction for f(x) = sin(2\\pi x)\n",
    "\n",
    "x_root = numpy.zeros(x.shape)\n",
    "N_steps = numpy.zeros(x.shape)\n",
    "for i,xk in enumerate(x):\n",
    "    x_root[i], x_root_array = newton(f, f_prime, xk)\n",
    "    N_steps[i] = len(x_root_array)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "y = numpy.linspace(-2,2)\n",
    "X,Y = numpy.meshgrid(x,y)\n",
    "X_root = numpy.outer(numpy.ones(y.shape),x_root)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.pcolor(X, Y, X_root,vmin=-5, vmax=5,cmap='seismic')\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('$x_{root}$', fontsize=18)\n",
    "plt.plot(x, f(x), 'k-')\n",
    "plt.plot(x, numpy.zeros(x.shape),'k--', linewidth=0.5)\n",
    "plt.xlabel('x', fontsize=16)\n",
    "plt.title('Basins of Attraction: $f(x) = \\sin{2\\pi x}$', fontsize=18)\n",
    "#plt.xlim(0.25-.1,0.25+.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Fractal Basins of Attraction\n",
    "\n",
    "If $f(x)$ is complex (for $x$ complex), then the basins of attraction can be beautiful and fractal\n",
    "\n",
    "Plotted below are two fairly simple equations which demonstrate the issue:\n",
    "1. $f(x) = x^3 - 1$\n",
    "2. Kepler's equation $\\theta - e \\sin \\theta = M$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x**3 - 1\n",
    "f_prime = lambda x: 3 * x**2\n",
    "\n",
    "N = 1001\n",
    "x = numpy.linspace(-2, 2, N)\n",
    "X, Y = numpy.meshgrid(x, x)\n",
    "R = X + 1j * Y\n",
    "\n",
    "for i in range(30):\n",
    "    R = R - f(R) / f_prime(R)\n",
    "\n",
    "#print(numpy.imag(R))\n",
    "roots = numpy.roots([1., 0., 0., -1])\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')\n",
    "#axes.contourf(X, Y, numpy.sign(numpy.imag(R))*numpy.abs(R),vmin = -10, vmax = 10)\n",
    "cf = axes.contourf(X, Y, numpy.real(R), vmin=-8, vmax=8)\n",
    "axes.scatter(numpy.real(roots), numpy.imag(roots))\n",
    "axes.set_xlabel(\"Real\")\n",
    "axes.set_ylabel(\"Imaginary\")\n",
    "axes.set_title(\"Basin of Attraction for $f(x) = x^3 - 1$\")\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def f(theta, e=0.083, M=1):\n",
    "    return theta - e * numpy.sin(theta) - M\n",
    "def f_prime(theta, e=0.083):\n",
    "    return 1 - e * numpy.cos(theta)\n",
    "\n",
    "N = 1001\n",
    "x = numpy.linspace(-30.5, -29.5, N)\n",
    "y = numpy.linspace(-17.5, -16.5, N)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "R = X + 1j * Y\n",
    "\n",
    "for i in range(30):\n",
    "    R = R - f(R) / f_prime(R)\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1, aspect='equal')    \n",
    "axes.contourf(X, Y, R, vmin = 0, vmax = 10)\n",
    "axes.set_xlabel(\"Real\")\n",
    "axes.set_ylabel(\"Imaginary\")\n",
    "axes.set_title(\"Basin of Attraction for $f(x) = x - e \\sin x - M$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Other Issues\n",
    "\n",
    "Need to supply both $f(x)$ and $f'(x)$, could be expensive\n",
    " \n",
    "Example:  FTV equation $f(r) = A - \\frac{m P}{r} \\left[ \\left(1 + \\frac{r}{m} \\right )^{m n} - 1\\right]$\n",
    "\n",
    "Can use symbolic differentiation (`sympy`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Secant Methods\n",
    "\n",
    "Is there a method with the convergence of Newton's method but without the extra derivatives?  What way would you modify Newton's method so that you would not need $f'(x)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given $x_k$ and $x_{k-1}$ represent the derivative as the approximation\n",
    "\n",
    "$$f'(x) \\approx \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}$$\n",
    "\n",
    "Combining this with the Newton approach leads to\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k) (x_k - x_{k-1}) }{f(x_k) - f(x_{k-1})}$$\n",
    "\n",
    "This leads to superlinear convergence and not quite quadratic as the exponent on the convergence is $\\approx 1.7$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Alternative interpretation, fit a line through two points and see where they intersect the x-axis.\n",
    "\n",
    "$$(x_k, f(x_k)) ~~~~~ (x_{k-1}, f(x_{k-1})$$\n",
    "\n",
    "$$y = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x - x_k) + b$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "$$b = f(x_{k-1}) - \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x_{k-1} - x_k)$$\n",
    "\n",
    "$$ y = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x - x_k) + f(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now solve for $x_{k+1}$ which is where the line intersects the x-axies ($y=0$)\n",
    "\n",
    "$$0 = \\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} (x_{k+1} - x_k) + f(x_k)$$\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k)  (x_k - x_{k-1})}{f(x_k) - f(x_{k-1})}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Secant Method\n",
    "\n",
    "$$x_{k+1} = x_k - \\frac{f(x_k)  (x_k - x_{k-1})}{f(x_k) - f(x_{k-1})}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "# Initial guess\n",
    "x_k = 0.07\n",
    "x_km = 0.06\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "axes.plot(x_k, 0.0, 'ko')\n",
    "axes.plot(x_k, f(x_k), 'ko')\n",
    "axes.plot([x_k, x_k], [0.0, f(x_k)], 'k--')\n",
    "axes.plot(x_km, 0.0, 'ko')\n",
    "axes.plot(x_km, f(x_km), 'ko')\n",
    "axes.plot([x_km, x_km], [0.0, f(x_km)], 'k--')\n",
    "\n",
    "axes.plot(r, (f(x_k) - f(x_km)) / (x_k - x_km) * (r - x_k) + f(x_k), 'k')\n",
    "x_kp = x_k - (f(x_k) * (x_k - x_km) / (f(x_k) - f(x_km)))\n",
    "axes.plot(x_kp, 0.0, 'ro')\n",
    "axes.plot([x_kp, x_kp], [0.0, f(x_kp)], 'r--')\n",
    "axes.plot(x_kp, f(x_kp), 'ro')\n",
    "\n",
    "axes.set_xlabel(\"r\", fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\", fontsize=14)\n",
    "axes.set_title(\"Secant Method\", fontsize=18)\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "What would the algorithm look like for such a method?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "\n",
    "Given $f(x)$,   a `TOLERANCE`, and a `MAX_STEPS` \n",
    "\n",
    "1. Initialize two points $x_0$, $x_1$, $f_0 = f(x_0)$, and $f_1 = f(x_1)$\n",
    "2. Loop for k=2, to  `MAX_STEPS` is reached or `TOLERANCE` is achieved\n",
    "   1. Calculate new update \n",
    "   $$x_{2} = x_1 - \\frac{f(x_1)  (x_1 - x_{0})}{f(x_1) - f(x_{0})}$$\n",
    "   2. Check for convergence and break if reached\n",
    "   3. Update parameters $x_0 = x_1$, $x_1 = x_{2}$, $f_0 = f_1$ and $f_1 = f(x_1)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Some Code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def secant(f, x0, x1, tol = 1.e-6):\n",
    "    \"\"\" uses a linear secant method to find a root x of a function of a single variable f\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    f: function f(x)\n",
    "        returns type: float\n",
    "    x0: float\n",
    "        first point to initialize the algorithm\n",
    "    x1: float\n",
    "        second point to initialize the algorithm x1 != x0        \n",
    "    tolerance: float\n",
    "        Returns when |f(x)| < tol\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    x: float\n",
    "        final iterate\n",
    "    x_array: numpy array\n",
    "        history of iteration points\n",
    "        \n",
    "    Raises:\n",
    "    -------\n",
    "    ValueError:\n",
    "        if x1 is too close to x0\n",
    "    Warning: \n",
    "        if number of iterations exceed MAX_STEPS\n",
    "    \"\"\"\n",
    "    MAX_STEPS = 200\n",
    "    \n",
    "    if numpy.isclose(x0, x1):\n",
    "        raise ValueError('Initial points are too close (preferably should be a bracket)')\n",
    "        \n",
    "    x_array = [ x0, x1 ]\n",
    "    for k in range(1, MAX_STEPS + 1):\n",
    "        x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))\n",
    "        x_array.append(x2)\n",
    "        if numpy.abs(f(x2)) < tol:\n",
    "            break\n",
    "        x0 = x1\n",
    "        x1 = x2\n",
    "        \n",
    "    if k == MAX_STEPS:\n",
    "        warnings.warn('Maximum number of steps exceeded')\n",
    "    \n",
    "    return x2, numpy.array(x_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Set the problem up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "P = 1500.0\n",
    "m = 12\n",
    "n = 20.0\n",
    "A = 1e6\n",
    "r = numpy.linspace(0.05, 0.11, 100)\n",
    "f = lambda r, A=A, m=m, P=P, n=n: \\\n",
    "        A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x0 = 0.06\n",
    "x1 = 0.07\n",
    "\n",
    "x, x_array = secant(f, x0, x1, tol= 1.e-7)\n",
    "print('x = {}, f(x) = {}, Nsteps = {}'.format(x, f(x), len(x_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "r = numpy.linspace(0.05, 0.10, 100)\n",
    "\n",
    "# Setup figure to plot convergence\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(r, f(r), 'b')\n",
    "axes.plot(r, numpy.zeros(r.shape),'r--')\n",
    "\n",
    "for n, x in enumerate(x_array):\n",
    "    axes.plot(x, f(x),'kx')\n",
    "    axes.text(x, f(x), str(n), fontsize=\"15\") \n",
    "          \n",
    "axes.set_xlabel(\"r\", fontsize=16)\n",
    "axes.set_ylabel(\"f(r)\", fontsize=16)\n",
    "axes.set_title(\"Secant Method Steps\", fontsize=18)\n",
    "axes.grid()\n",
    "axes.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.semilogy(numpy.arange(len(x_array)), numpy.abs(f(x_array)), 'bo-')\n",
    "axes.grid()\n",
    "axes.set_xlabel('Iterations', fontsize=16)\n",
    "axes.set_ylabel('Residual $|f(r)|$', fontsize=16)\n",
    "axes.set_title('Convergence', fontsize=18)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Comments\n",
    "\n",
    " - Secant method as shown is equivalent to linear interpolation\n",
    " - Can use higher order interpolation for higher order secant methods\n",
    " - Convergence is not quite quadratic\n",
    " - Not guaranteed to converge\n",
    " - Does not preserve brackets\n",
    " - Almost as good as Newton's method if your initial guess is good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hybrid Methods\n",
    "\n",
    "Combine attributes of methods with others to make one great algorithm to rule them all (not really)\n",
    "\n",
    "#### Goals\n",
    "1. Robustness:  Given a bracket $[a,b]$, maintain bracket\n",
    "1. Efficiency:  Use superlinear convergent methods when possible\n",
    "\n",
    "#### Options\n",
    " - Methods requiring $f'(x)$\n",
    "   - NewtSafe (RootSafe, Numerical Recipes)\n",
    "   - Newton's Method within a bracket, Bisection otherwise\n",
    " - Methods not requiring $f'(x)$\n",
    "   - Brent's Algorithm (zbrent, Numerical Recipes)\n",
    "     - Combination of bisection, secant and inverse quadratic interpolation\n",
    "   - `scipy.optimize` package **new** root_scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import root_scalar\n",
    "#root_scalar?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Set the problem up (again)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def f(r,A,m,P,n):\n",
    "    return  A - m * P / r * ((1.0 + r / m)**(m * n) - 1.0)\n",
    "\n",
    "def f_prime(r,A,m,P,n): \n",
    "    return (-P*m*n*(1.0 + r/m)**(m*n)/(r*(1.0 + r/m)) + \n",
    "            P*m*((1.0 + r/m)**(m*n) - 1.0)/r**2)\n",
    "A = 1.e6\n",
    "m = 12\n",
    "P = 1500.\n",
    "n = 20."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Try Brent's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = 0.07\n",
    "b = 0.1\n",
    "sol  = root_scalar(f,args=(A,m,P,n), bracket=(a, b), method='brentq')  \n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Try Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sol = root_scalar(f,args=(A,m,P,n), x0=.07, fprime=f_prime, method='newton')\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "  f1 = lambda A : f(r,A,m,P,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optimization (finding extrema)\n",
    "\n",
    "I want to find the extrema of a function $f(x)$ on a given interval $[a,b]$.\n",
    "\n",
    "A few approaches:\n",
    " - Interpolation Algorithms: Repeated parabolic interpolation\n",
    " - Bracketing Algorithms:  Golden-Section Search (linear)\n",
    " - Hybrid Algorithms \n",
    " - Use Rootfinding methods on $f'(x)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Interpolation Approach\n",
    "\n",
    "Successive parabolic interpolation - similar to secant method\n",
    "\n",
    "Basic idea:  Fit polynomial to function using three points, find its minima, and guess new points based on that minima"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1. What do we need to fit a polynomial $p_n(x)$ of degree $n \\geq 2$?\n",
    "\n",
    "2. How do we construct the polynomial $p_2(x)$?\n",
    "\n",
    "3. Once we have constructed $p_2(x)$ how would we find the minimum?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "\n",
    "Given $f(x)$ and $[x_0,x_1]$ - Note that unlike a bracket these will be a sequence of better approximations to the minimum.\n",
    "1. Initialize $x = [x_0, x_1, (x_0+x_1)/2]$\n",
    "1. Loop\n",
    "   1. Evaluate function $f(x)$ at the three points\n",
    "   1. Find the quadratic polynomial that interpolates those points: \n",
    "   \n",
    "      $$p(x) = p_0 x^2 + p_1 x + p_2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "   3. Calculate the minimum: \n",
    "      \n",
    "      $$p'(x) = 2 p_0 x + p_1 = 0 \\quad \\Rightarrow \\quad x^\\ast = -p_1 / (2 p_0)$$\n",
    "\n",
    "   1. New set of points $x = [x_1, (x_0+x_1)/2, x^\\ast]$\n",
    "   1. Check tolerance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def f(t):\n",
    "    \"\"\"Simple function for minimization demos\"\"\"\n",
    "    return -3.0 * numpy.exp(-(t - 0.3)**2 / (0.1)**2) \\\n",
    "           +      numpy.exp(-(t - 0.6)**2 / (0.2)**2) \\\n",
    "           +      numpy.exp(-(t - 1.0)**2 / (0.2)**2) \\\n",
    "           +      numpy.sin(t)                        \\\n",
    "           -      2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x0, x1 = 0.5, 0.2\n",
    "x = numpy.array([x0, x1, (x0 + x1)/2.])\n",
    "p = numpy.polyfit(x, f(x), 2)\n",
    "parabola = lambda t: p[0]*t**2 + p[1]*t + p[2]\n",
    "t_min = -p[1]/2./p[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "MAX_STEPS = 100\n",
    "TOLERANCE = 1e-4\n",
    "\n",
    "t = numpy.linspace(0., 2., 100)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t), label='$f(t)$')\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.plot(x[0], f(x[0]), 'ko')\n",
    "axes.plot(x[1], f(x[1]), 'ko')\n",
    "axes.plot(x[2], f(x[2]), 'ko')\n",
    "axes.plot(t, parabola(t), 'r--',  label='parabola')\n",
    "axes.plot(t_min, parabola(t_min), 'ro' )\n",
    "axes.plot(t_min, f(t_min), 'k+')\n",
    "\n",
    "axes.legend(loc='best')\n",
    "axes.set_ylim((-5, 0.0))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Rinse and repeat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "MAX_STEPS = 100\n",
    "TOLERANCE = 1e-4\n",
    "\n",
    "x = numpy.array([x0, x1, (x0 + x1) / 2.0])\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.plot(x[0], f(x[0]), 'ko')\n",
    "axes.plot(x[1], f(x[1]), 'ko')\n",
    "\n",
    "success = False\n",
    "for n in range(1, MAX_STEPS + 1):\n",
    "    axes.plot(x[2], f(x[2]), 'ko')\n",
    "    poly = numpy.polyfit(x, f(x), 2)\n",
    "    axes.plot(t, poly[0] * t**2 + poly[1] * t + poly[2], 'r--')\n",
    "    x[0] = x[1]\n",
    "    x[1] = x[2]\n",
    "    x[2] = -poly[1] / (2.0 * poly[0])\n",
    "    if numpy.abs(x[2] - x[1]) / numpy.abs(x[2]) < TOLERANCE:\n",
    "        success = True\n",
    "        break\n",
    "    \n",
    "if success:\n",
    "    print(\"Success!\")\n",
    "    print(\"  t* = %s\" % x[2])\n",
    "    print(\"  f(t*) = %s\" % f(x[2]))\n",
    "    print(\"  number of steps = %s\" % n)\n",
    "else:\n",
    "    print(\"Reached maximum number of steps!\")\n",
    "    \n",
    "axes.set_ylim((-5, 0.0))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Some Code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def parabolic_interpolation(f, bracket, tol = 1.e-6):\n",
    "    \"\"\" uses repeated parabolic interpolation to refine a local minimum of a function f(x)\n",
    "    this routine uses numpy functions polyfit and polyval to fit and evaluate the quadratics\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    f: function f(x)\n",
    "        returns type: float\n",
    "    bracket: array\n",
    "        array [x0, x1] containing an initial bracket that contains a minimum   \n",
    "    tolerance: float\n",
    "        Returns when relative error of last two iterates < tol \n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    x: float\n",
    "        final estimate of the minima\n",
    "    x_array: numpy array\n",
    "        history of iteration points\n",
    "        \n",
    "    Raises:\n",
    "    -------\n",
    "    Warning: \n",
    "        if number of iterations exceed MAX_STEPS\n",
    "    \"\"\"\n",
    "    MAX_STEPS = 100\n",
    "    \n",
    "    x = numpy.zeros(3)\n",
    "    x[:2] = bracket\n",
    "    x[2] = (x[0] + x[1])/2.\n",
    "        \n",
    "    x_array = [ x[2] ]\n",
    "    for k in range(1, MAX_STEPS + 1):\n",
    "        poly = numpy.polyfit(x, f(x), 2)\n",
    "        x[0] = x[1]\n",
    "        x[1] = x[2]\n",
    "        x[2] = -poly[1] / (2.0 * poly[0])\n",
    "        x_array.append(x[2])\n",
    "        if numpy.abs(x[2] - x[1]) / numpy.abs(x[2]) < tol:\n",
    "            break\n",
    "        \n",
    "    if k == MAX_STEPS:\n",
    "        warnings.warn('Maximum number of steps exceeded')\n",
    "    \n",
    "    return x[2], numpy.array(x_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### set up problem\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "bracket = numpy.array([0.2, 0.5])\n",
    "x, x_array = parabolic_interpolation(f, bracket, tol = 1.e-6)\n",
    "print(\"Extremum f(x) = {}, at x = {}, N steps = {}\".format(f(x), x, len(x_array)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0, 2, 200)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.plot(x_array, f(x_array),'ro')\n",
    "axes.plot(x, f(x), 'go')\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bracketing Algorithm (Golden Section Search)\n",
    "\n",
    "Given $f(x) \\in C[x_0,x_3]$ that is convex (concave) over an interval $x \\in [x_0,x_3]$ reduce the interval size until it brackets the minimum (maximum).\n",
    "\n",
    "Note that we no longer have the $x=0$ help we had before so bracketing and doing bisection is a bit trickier in this case.  In particular choosing your initial bracket is important!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Bracket Picking\n",
    "\n",
    "Say we start with a bracket $[x_0, x_3]$ and pick two new points $x_1 < x_2 \\in [x_0, x_3]$.  We want to pick a new bracket that guarantees that the extrema exists in it.  We then can pick this new bracket with the following rules:\n",
    " - If $f(x_1) < f(x_2)$ then we know the minimum is between $x_0$ and $x_2$.\n",
    " - If $f(x_1) > f(x_2)$ then we know the minimum is between $x_1$ and $x_3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "fig.set_figheight(fig.get_figheight() * 2)\n",
    "\n",
    "search_points = [-1.0, -0.5, 0.75, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 1)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) < f(x_2) \\Rightarrow [x_0, x_2]$\")\n",
    "\n",
    "search_points = [-1.0, -0.75, 0.5, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 2)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) > f(x_2) \\Rightarrow [x_1, x_3]$\")\n",
    "\n",
    "search_points = [-1.0, 0.25, 0.75, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 3)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) < f(x_2) \\Rightarrow [x_0, x_2]$\")\n",
    "\n",
    "search_points = [-1.0, -0.75, -0.25, 1.0]\n",
    "axes = fig.add_subplot(2, 2, 4)\n",
    "x = numpy.linspace(search_points[0] - 0.1, search_points[-1] + 0.1, 100)\n",
    "axes.plot(x, f(x), 'b')\n",
    "for (i, point) in enumerate(search_points):\n",
    "    axes.plot(point, f(point),'or')\n",
    "    axes.text(point + 0.05, f(point), str(i))\n",
    "axes.plot(0, 0, 'sk')\n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_title(\"$f(x_1) > f(x_2) \\Rightarrow [x_1, x_3]$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Picking Brackets and Points\n",
    "\n",
    "Again say we have a bracket $[x_0,x_3]$ and suppose we have two new search points $x_1$ and $x_2$ that separates $[x_0,x_3]$ into two new overlapping brackets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Define:  the length of the line segments in the interval\n",
    "\\begin{aligned}\n",
    "    a &= x_1 - x_0, \\\\\n",
    "    b &= x_2 - x_1,\\\\\n",
    "    c &= x_3 - x_2 \\\\\n",
    "\\end{aligned}\n",
    "and the total bracket length\n",
    "\\begin{aligned}\n",
    "    d &= x_3 - x_0. \\\\\n",
    "\\end{aligned}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: (x - 0.25)**2 + 0.5\n",
    "phi = (numpy.sqrt(5.0) - 1.) / 2.0\n",
    "\n",
    "x = [-1.0, None, None, 1.0]\n",
    "x[1] = x[3] - phi * (x[3] - x[0])\n",
    "x[2] = x[0] + phi * (x[3] - x[0])\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "t = numpy.linspace(-2.0, 2.0, 100)\n",
    "axes.plot(t, f(t), 'k')\n",
    "# First set of intervals\n",
    "axes.plot([x[0], x[1]], [0.0, 0.0], 'g',label='a')\n",
    "axes.plot([x[1], x[2]], [0.0, 0.0], 'r', label='b')\n",
    "axes.plot([x[2], x[3]], [0.0, 0.0], 'b', label='c')\n",
    "axes.plot([x[0], x[3]], [2.5, 2.5], 'c', label='d')\n",
    "\n",
    "axes.plot([x[0], x[0]], [0.0, f(x[0])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'b--')\n",
    "axes.plot([x[3], x[3]], [0.0, f(x[3])], 'b--')\n",
    "axes.plot([x[0], x[0]], [2.5, f(x[0])], 'c--')\n",
    "axes.plot([x[3], x[3]], [2.5, f(x[3])], 'c--')\n",
    "\n",
    "points = [ (x[0] + x[1])/2., (x[1] + x[2])/2., (x[2] + x[3])/2., (x[0] + x[3])/2. ]\n",
    "y      = [ 0., 0., 0., 2.5]\n",
    "labels = [ 'a', 'b', 'c', 'd']\n",
    "\n",
    "for (n, point)  in enumerate(points):\n",
    "    axes.text(point, y[n] + 0.1, labels[n], fontsize=15)\n",
    "\n",
    "for (n, point) in enumerate(x):\n",
    "    axes.plot(point, f(point), 'ok')\n",
    "    axes.text(point, f(point)+0.1, n, fontsize='15')\n",
    "    \n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_ylim((-1.0, 3.0))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For **Golden Section Search** we require two conditions:\n",
    " - The two new possible brackets are of equal length.  i.e $[x_0, x_2] =  [x_1, x_3]$ or\n",
    " \n",
    " $$\n",
    "     a + b = b + c\n",
    " $$\n",
    " \n",
    " or simply $a = c$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: (x - 0.25)**2 + 0.5\n",
    "phi = (numpy.sqrt(5.0) - 1.) / 2.0\n",
    "\n",
    "x = [-1.0, None, None, 1.0]\n",
    "x[1] = x[3] - phi * (x[3] - x[0])\n",
    "x[2] = x[0] + phi * (x[3] - x[0])\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "t = numpy.linspace(-2.0, 2.0, 100)\n",
    "axes.plot(t, f(t), 'k')\n",
    "# First set of intervals\n",
    "axes.plot([x[0], x[1]], [0.0, 0.0], 'g',label='a')\n",
    "axes.plot([x[1], x[2]], [0.0, 0.0], 'r', label='b')\n",
    "axes.plot([x[2], x[3]], [0.0, 0.0], 'b', label='c')\n",
    "axes.plot([x[0], x[3]], [2.5, 2.5], 'c', label='d')\n",
    "\n",
    "axes.plot([x[0], x[0]], [0.0, f(x[0])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'b--')\n",
    "axes.plot([x[3], x[3]], [0.0, f(x[3])], 'b--')\n",
    "axes.plot([x[0], x[0]], [2.5, f(x[0])], 'c--')\n",
    "axes.plot([x[3], x[3]], [2.5, f(x[3])], 'c--')\n",
    "\n",
    "points = [ (x[0] + x[1])/2., (x[1] + x[2])/2., (x[2] + x[3])/2., (x[0] + x[3])/2. ]\n",
    "y      = [ 0., 0., 0., 2.5]\n",
    "labels = [ 'a', 'b', 'c', 'd']\n",
    "\n",
    "for (n, point)  in enumerate(points):\n",
    "    axes.text(point, y[n] + 0.1, labels[n], fontsize=15)\n",
    "\n",
    "for (n, point) in enumerate(x):\n",
    "    axes.plot(point, f(point), 'ok')\n",
    "    axes.text(point, f(point)+0.1, n, fontsize='15')\n",
    "    \n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_ylim((-1.0, 3.0))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " \n",
    " - The ratio of segment lengths is the same for every level of recursion so the problem is self-similar i.e.\n",
    " \n",
    " $$\n",
    "     \\frac{b}{a} = \\frac{c}{a + b}\n",
    " $$\n",
    " \n",
    " These two requirements will allow maximum reuse of previous points and require adding only one new point $x^*$ at each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: (x - 0.25)**2 + 0.5\n",
    "phi = (numpy.sqrt(5.0) - 1.) / 2.0\n",
    "\n",
    "x = [-1.0, None, None, 1.0]\n",
    "x[1] = x[3] - phi * (x[3] - x[0])\n",
    "x[2] = x[0] + phi * (x[3] - x[0])\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = []\n",
    "axes.append(fig.add_subplot(1, 2, 1))\n",
    "axes.append(fig.add_subplot(1, 2, 2))\n",
    "t = numpy.linspace(-2.0, 2.0, 100)\n",
    "for i in range(2):\n",
    "    axes[i].plot(t, f(t), 'k')\n",
    "\n",
    "    # First set of intervals\n",
    "    axes[i].plot([x[0], x[2]], [0.0, 0.0], 'g')\n",
    "    axes[i].plot([x[1], x[3]], [-0.2, -0.2], 'r')\n",
    "    axes[i].plot([x[0], x[0]], [0.0, f(x[0])], 'g--')\n",
    "    axes[i].plot([x[2], x[2]], [0.0, f(x[2])], 'g--')\n",
    "    axes[i].plot([x[1], x[1]], [-0.2, f(x[1])], 'r--')\n",
    "    axes[i].plot([x[3], x[3]], [-0.2, f(x[3])], 'r--')\n",
    "    for (n, point) in enumerate(x):\n",
    "        axes[i].plot(point, f(point), 'ok')\n",
    "        axes[i].text(point, f(point)+0.1, n, fontsize='15')\n",
    "\n",
    "    axes[i].set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "    axes[i].set_ylim((-1.0, 3.0))\n",
    "\n",
    "# Left new interval\n",
    "x_new = [x[0], None, x[1], x[2]]\n",
    "x_new[1] = phi * (x[1] - x[0]) + x[0]\n",
    "\n",
    "#axes[0].plot([x_new[0], x_new[2]], [1.5, 1.5], 'b')\n",
    "#axes[0].plot([x_new[1], x_new[3]], [1.75, 1.75], 'c')\n",
    "#axes[0].plot([x_new[0], x_new[0]], [1.5, f(x_new[0])], 'b--')\n",
    "#axes[0].plot([x_new[2], x_new[2]], [1.5, f(x_new[2])], 'b--')\n",
    "#axes[0].plot([x_new[1], x_new[1]], [1.75, f(x_new[1])], 'c--')\n",
    "#axes[0].plot([x_new[3], x_new[3]], [1.75, f(x_new[3])], 'c--')\n",
    "axes[0].plot(x_new[1], f(x_new[1]), 'ko')\n",
    "axes[0].text(x_new[1], f(x_new[1]) + 0.1, \"*\", fontsize='15')\n",
    "for i in range(4):\n",
    "    axes[0].text(x_new[i], -0.5, i, color='g',fontsize='15')\n",
    "\n",
    "    \n",
    "# Right new interval\n",
    "x_new = [x[1], x[2], None, x[3]]\n",
    "x_new[2] = (x[2] - x[1]) * phi + x[2]\n",
    "\n",
    "#axes[1].plot([x_new[0], x_new[2]], [1.25, 1.25], 'b')\n",
    "#axes[1].plot([x_new[1], x_new[3]], [1.5, 1.5], 'c')\n",
    "#axes[1].plot([x_new[0], x_new[0]], [1.25, f(x_new[0])], 'b--')\n",
    "#axes[1].plot([x_new[2], x_new[2]], [1.25, f(x_new[2])], 'b--')\n",
    "#axes[1].plot([x_new[1], x_new[1]], [1.5, f(x_new[2])], 'c--')\n",
    "#axes[1].plot([x_new[3], x_new[3]], [1.5, f(x_new[3])], 'c--')\n",
    "axes[1].plot(x_new[2], f(x_new[2]), 'ko')\n",
    "axes[1].text(x_new[2], f(x_new[2]) + 0.1, \"*\", fontsize='15')\n",
    "for i in range(4):\n",
    "    axes[1].text(x_new[i], -0.5, i, color='r',fontsize='15')\n",
    "axes[0].set_title('Choose left bracket', fontsize=18)\n",
    "axes[1].set_title('Choose right bracket', fontsize=18)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As the first rule implies that $a = c$,  we can substitute into the second rule to yield\n",
    "\n",
    "$$\n",
    "     \\frac{b}{a} = \\frac{a}{a + b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or inverting and rearranging \n",
    "$$\n",
    "    \\frac{a}{b} = 1 + \\frac{b}{a}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "if we let the ratio $b/a = x$,  then \n",
    "$$\n",
    "    x + 1 = \\frac{1}{x} \\quad \\text{or} \\quad x^2 + x - 1 = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "     x^2 + x - 1 = 0\n",
    "$$\n",
    "has a single positive root for \n",
    "$$\n",
    "    x = \\frac{\\sqrt{5} - 1}{2} = \\varphi = 0.6180339887498949\n",
    "$$\n",
    "where $\\varphi$ is related to the  \"golden ratio\" (which in most definitions is given by $1+\\varphi$, but either work as $ 1+\\varphi = 1/\\varphi $ )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Subsequent proportionality implies that the distances between the 4 points at one iteration is proportional to the next.  We can now use all of our information to find the points $x_1$ and $x_2$ given any overall bracket $[x_0, x_3]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given $b/a = \\varphi$,   $a = c$, and the known width of the bracket $d$ it follows that\n",
    "\n",
    "$$ d = a + b + c = (2 + \\phi)a $$\n",
    "or \n",
    "$$ a = \\frac{d}{2 + \\varphi} = \\frac{\\varphi}{1 + \\varphi} d$$\n",
    "by the rather special properties of $\\varphi$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We could use this result immediately to find \n",
    "\n",
    "\\begin{align}\n",
    "    x_1 &= x_0 + a \\\\\n",
    "    x_2 &= x_3 - a \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Equivalently, you can show that \n",
    "\n",
    "$$a + b = (1 + \\varphi)a = \\varphi d$$\n",
    "\n",
    "so\n",
    "\n",
    "\\begin{align}\n",
    "    x_1 &= x_3 - \\varphi d \\\\\n",
    "    x_2 &= x_0 + \\varphi d \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda x: (x - 0.25)**2 + 0.5\n",
    "phi = (numpy.sqrt(5.0) - 1.) / 2.0\n",
    "\n",
    "x = [-1.0, None, None, 1.0]\n",
    "x[1] = x[3] - phi * (x[3] - x[0])\n",
    "x[2] = x[0] + phi * (x[3] - x[0])\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "t = numpy.linspace(-2.0, 2.0, 100)\n",
    "axes.plot(t, f(t), 'k')\n",
    "# First set of intervals\n",
    "axes.plot([x[0], x[1]], [0.0, 0.0], 'g',label='a')\n",
    "axes.plot([x[1], x[2]], [0.0, 0.0], 'r', label='b')\n",
    "axes.plot([x[2], x[3]], [0.0, 0.0], 'b', label='c')\n",
    "axes.plot([x[0], x[3]], [2.5, 2.5], 'c', label='d')\n",
    "\n",
    "axes.plot([x[0], x[0]], [0.0, f(x[0])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'g--')\n",
    "axes.plot([x[1], x[1]], [0.0, f(x[1])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'r--')\n",
    "axes.plot([x[2], x[2]], [0.0, f(x[2])], 'b--')\n",
    "axes.plot([x[3], x[3]], [0.0, f(x[3])], 'b--')\n",
    "axes.plot([x[0], x[0]], [2.5, f(x[0])], 'c--')\n",
    "axes.plot([x[3], x[3]], [2.5, f(x[3])], 'c--')\n",
    "\n",
    "points = [ (x[0] + x[1])/2., (x[1] + x[2])/2., (x[2] + x[3])/2., (x[0] + x[3])/2. ]\n",
    "y      = [ 0., 0., 0., 2.5]\n",
    "labels = [ 'a', 'b', 'c', 'd']\n",
    "\n",
    "for (n, point)  in enumerate(points):\n",
    "    axes.text(point, y[n] + 0.1, labels[n], fontsize=15)\n",
    "\n",
    "for (n, point) in enumerate(x):\n",
    "    axes.plot(point, f(point), 'ok')\n",
    "    axes.text(point, f(point)+0.1, n, fontsize='15')\n",
    "    \n",
    "axes.set_xlim((search_points[0] - 0.1, search_points[-1] + 0.1))\n",
    "axes.set_ylim((-1.0, 3.0))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Algorithm\n",
    "1. Initialize bracket $[x_0,x_3]$\n",
    "1. Initialize points $x_1 = x_3 - \\varphi (x_3 - x_0)$ and $x_2 = x_0 + \\varphi (x_3 - x_0)$\n",
    "1. Loop\n",
    "   1. Evaluate $f_1$ and $f_2$\n",
    "   1. If $f_1 < f_2$ then we pick the left interval for the next iteration\n",
    "   1. and otherwise pick the right interval\n",
    "   1. Check size of bracket for convergence $x_3 - x_0 <$ `TOLERANCE`\n",
    "   1. calculate the appropriate new point $x^*$ ($x_1$ on left, $x_2$ on right)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def golden_section(f, bracket, tol = 1.e-6):\n",
    "    \"\"\" uses golden section search to refine a local minimum of a function f(x)\n",
    "    this routine uses numpy functions polyfit and polyval to fit and evaluate the quadratics\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    f: function f(x)\n",
    "        returns type: float\n",
    "    bracket: array\n",
    "        array [x0, x3] containing an initial bracket that contains a minimum   \n",
    "    tolerance: float\n",
    "        Returns when | x3 - x0 | < tol \n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    x: float\n",
    "        final estimate of the midpoint of the bracket\n",
    "    x_array: numpy array\n",
    "        history of midpoint of each bracket\n",
    "        \n",
    "    Raises:\n",
    "    -------\n",
    "    ValueError:\n",
    "        If initial bracket is < tol or doesn't appear to have any interior points \n",
    "        that are less than the outer points\n",
    "    Warning: \n",
    "        if number of iterations exceed MAX_STEPS\n",
    "    \"\"\"\n",
    "    MAX_STEPS = 100\n",
    "    phi = (numpy.sqrt(5.0) - 1.) / 2.0\n",
    "\n",
    "    x = [ bracket[0], None, None, bracket[1] ]\n",
    "    \n",
    "    delta_x = x[3] - x[0]\n",
    "    x[1] = x[3] - phi * delta_x\n",
    "    x[2] = x[0] + phi * delta_x\n",
    "    \n",
    "    # check for initial bracket\n",
    "    fx = f(numpy.array(x))\n",
    "    bracket_min = min(fx[0], fx[3])\n",
    "    if fx[1] > bracket_min and fx[2] > bracket_min:\n",
    "        raise ValueError(\"interval does not appear to include a minimum\")\n",
    "    elif delta_x < tol:\n",
    "        raise ValueError(\"interval is already smaller than tol\")\n",
    "\n",
    "    x_mid = (x[3] + x[0])/2.\n",
    "    x_array = [ x_mid ]\n",
    "    for k in range(1, MAX_STEPS + 1):\n",
    "        f_1 = f(x[1])\n",
    "        f_2 = f(x[2]) \n",
    "    \n",
    "        if f_1 < f_2:\n",
    "            # Pick the left bracket\n",
    "            x_new = [x[0], None, x[1], x[2]]\n",
    "            delta_x = x_new[3] - x_new[0]\n",
    "            x_new[1] = x_new[3] - phi * delta_x\n",
    "        else:\n",
    "            # Pick the right bracket\n",
    "            x_new = [x[1], x[2], None, x[3]]\n",
    "            delta_x = x_new[3] - x_new[0]\n",
    "            x_new[2] = x_new[0] + phi * delta_x\n",
    "        x = x_new\n",
    "        x_array.append((x[3] + x[0])/ 2.)\n",
    "        \n",
    "        if numpy.abs(x[3] - x[0]) < tol:\n",
    "            break\n",
    "       \n",
    "    if k == MAX_STEPS:\n",
    "        warnings.warn('Maximum number of steps exceeded')\n",
    "    \n",
    "    return x_array[-1], numpy.array(x_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def f(t):\n",
    "    \"\"\"Simple function for minimization demos\"\"\"\n",
    "    return -3.0 * numpy.exp(-(t - 0.3)**2 / (0.1)**2) \\\n",
    "           +      numpy.exp(-(t - 0.6)**2 / (0.2)**2) \\\n",
    "           +      numpy.exp(-(t - 1.0)**2 / (0.2)**2) \\\n",
    "           +      numpy.sin(t)                        \\\n",
    "           -      2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x, x_array = golden_section(f,[0.2, .5], 1.e-4)\n",
    "print('t* = {}, f(t*) = {}, N steps = {}'.format(x, f(x), len(x_array)-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0, 2, 200)\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, f(t))\n",
    "axes.grid()\n",
    "axes.set_xlabel(\"t (days)\")\n",
    "axes.set_ylabel(\"People (N)\")\n",
    "axes.set_title(\"Decrease in Population due to SPAM Poisoning\")\n",
    "axes.plot(x_array, f(x_array),'ko')\n",
    "axes.plot(x_array[0],f(x_array[0]),'ro')\n",
    "axes.plot(x_array[-1],f(x_array[-1]),'go')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Scipy Optimization\n",
    "\n",
    "Scipy contains a lot of ways for optimization.   But a convenenient interface for minimization of functions of a single variable is `scipy.optimize.minimize_scalar`\n",
    "\n",
    "For optimization or constrained optimization for functions of more than one variable, see \n",
    "`scipy.optimized.minimize`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from  scipy.optimize import minimize_scalar\n",
    "minimize_scalar?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Try some different methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sol = minimize_scalar(f, bracket=(0.2, 0.25, 0.5), method='golden')\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sol = minimize_scalar(f, method='brent', bracket=(0.2, 0.25, 0.5))\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "sol = minimize_scalar(f, bounds=(0.,0.5), method='bounded')\n",
    "print(sol)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}