scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='sp'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SciPy\n",
    "\n",
    "\n",
    "<a id='index-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [SciPy](#SciPy)  \n",
    "  - [SciPy versus NumPy](#SciPy-versus-NumPy)  \n",
    "  - [Statistics](#Statistics)  \n",
    "  - [Roots and Fixed Points](#Roots-and-Fixed-Points)  \n",
    "  - [Optimization](#Optimization)  \n",
    "  - [Integration](#Integration)  \n",
    "  - [Linear Algebra](#Linear-Algebra)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[SciPy](http://www.scipy.org) builds on top of NumPy to provide common tools for scientific programming such as\n",
    "\n",
    "- [linear algebra](http://docs.scipy.org/doc/scipy/reference/linalg.html)  \n",
    "- [numerical integration](http://docs.scipy.org/doc/scipy/reference/integrate.html)  \n",
    "- [interpolation](http://docs.scipy.org/doc/scipy/reference/interpolate.html)  \n",
    "- [optimization](http://docs.scipy.org/doc/scipy/reference/optimize.html)  \n",
    "- [distributions and random number generation](http://docs.scipy.org/doc/scipy/reference/stats.html)  \n",
    "- [signal processing](http://docs.scipy.org/doc/scipy/reference/signal.html)  \n",
    "- etc., etc  \n",
    "\n",
    "\n",
    "Like NumPy, SciPy is stable, mature and widely used.\n",
    "\n",
    "Many SciPy routines are thin wrappers around industry-standard Fortran libraries such as LAPACK, BLAS, etc.\n",
    "\n",
    "It&rsquo;s not really necessary to &ldquo;learn&rdquo; SciPy as a whole.\n",
    "\n",
    "A more common approach is to get some idea of what&rsquo;s in the library and then look up [documentation](http://docs.scipy.org/doc/scipy/reference/index.html) as required.\n",
    "\n",
    "In this lecture, we aim only to highlight some useful parts of the package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SciPy versus NumPy\n",
    "\n",
    "SciPy is a package that contains various tools that are built on top of NumPy, using its array data type and related functionality.\n",
    "\n",
    "In fact, when we import SciPy we also get NumPy, as can be seen from the SciPy initialization file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Import numpy symbols to scipy namespace\n",
    "import numpy as _num\n",
    "linalg = None\n",
    "from numpy import *\n",
    "from numpy.random import rand, randn\n",
    "from numpy.fft import fft, ifft\n",
    "from numpy.lib.scimath import *\n",
    "\n",
    "__all__  = []\n",
    "__all__ += _num.__all__\n",
    "__all__ += ['randn', 'rand', 'fft', 'ifft']\n",
    "\n",
    "del _num\n",
    "# Remove the linalg imported from numpy so that the scipy.linalg package\n",
    "# can be imported.\n",
    "del linalg\n",
    "__all__.remove('linalg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it&rsquo;s more common and better practice to use NumPy functionality explicitly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "a = np.identity(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is useful in SciPy is the functionality in its sub-packages\n",
    "\n",
    "- `scipy.optimize`, `scipy.integrate`, `scipy.stats`, etc.  \n",
    "\n",
    "\n",
    "These sub-packages and their attributes need to be imported separately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "from scipy.optimize import brentq\n",
    "# etc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let&rsquo;s explore some of the major sub-packages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics\n",
    "\n",
    "\n",
    "<a id='index-4'></a>\n",
    "The `scipy.stats` subpackage supplies\n",
    "\n",
    "- numerous random variable objects (densities, cumulative distributions, random sampling, etc.)  \n",
    "- some estimation procedures  \n",
    "- some statistical tests  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Variables and Distributions\n",
    "\n",
    "Recall that `numpy.random` provides functions for generating random variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.random.beta(5, 5, size=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This generates a draw from the distribution below when `a, b = 5, 5`\n",
    "\n",
    "\n",
    "<a id='equation-betadist2'></a>\n",
    "$$\n",
    "f(x; a, b) = \\frac{x^{(a - 1)} (1 - x)^{(b - 1)}}\n",
    "    {\\int_0^1 u^{(a - 1)} (1 - u)^{(b - 1)} du}\n",
    "    \\qquad (0 \\leq x \\leq 1) \\tag{1}\n",
    "$$\n",
    "\n",
    "Sometimes we need access to the density itself, or the cdf, the quantiles, etc.\n",
    "\n",
    "For this, we can use `scipy.stats`, which provides all of this functionality as well as random number generation in a single consistent interface.\n",
    "\n",
    "Here&rsquo;s an example of usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import beta\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "q = beta(5, 5)      # Beta(a, b), with a = b = 5\n",
    "obs = q.rvs(2000)   # 2000 observations\n",
    "grid = np.linspace(0.01, 0.99, 100)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.hist(obs, bins=40, density=True)\n",
    "ax.plot(grid, q.pdf(grid), 'k-', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this code, we created a so-called `rv_frozen` object, via the call `q = beta(5, 5)`.\n",
    "\n",
    "The &ldquo;frozen&rdquo; part of the notation implies that `q` represents a particular distribution with a particular set of parameters.\n",
    "\n",
    "Once we&rsquo;ve done so, we can then generate random numbers, evaluate the density, etc., all from this fixed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "q.cdf(0.4)      # Cumulative distribution function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "q.pdf(0.4)      # Density function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "q.ppf(0.8)      # Quantile (inverse cdf) function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "q.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general syntax for creating these objects is\n",
    "\n",
    "> `identifier = scipy.stats.distribution_name(shape_parameters)`\n",
    "\n",
    "\n",
    "where `distribution_name` is one of the distribution names in [scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html).\n",
    "\n",
    "There are also two keyword arguments, `loc` and `scale`, which following our example above, are called as\n",
    "\n",
    "> `identifier = scipy.stats.distribution_name(shape_parameters, loc=c, scale=d)`\n",
    "\n",
    "\n",
    "These transform the original random variable $ X $ into $ Y = c + d X $.\n",
    "\n",
    "The methods `rvs`, `pdf`, `cdf`, etc. are transformed accordingly.\n",
    "\n",
    "Before finishing this section, we note that there is an alternative way of calling the methods described above.\n",
    "\n",
    "For example, the previous code can be replaced by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "obs = beta.rvs(5, 5, size=2000)\n",
    "grid = np.linspace(0.01, 0.99, 100)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(obs, bins=40, density=True)\n",
    "ax.plot(grid, beta.pdf(grid, 5, 5), 'k-', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Goodies in scipy.stats\n",
    "\n",
    "There are a variety statistical functions in `scipy.stats`.\n",
    "\n",
    "For example, `scipy.stats.linregress` implements simple linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "x = np.random.randn(200)\n",
    "y = 2 * x + 0.1 * np.random.randn(200)\n",
    "gradient, intercept, r_value, p_value, std_err = linregress(x, y)\n",
    "gradient, intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the full list, consult the [documentation](https://docs.scipy.org/doc/scipy/reference/stats.html#statistical-functions-scipy-stats)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roots and Fixed Points\n",
    "\n",
    "A *root* of a real function $ f $ on $ [a,b] $ is an $ x \\in [a, b] $ such that $ f(x)=0 $.\n",
    "\n",
    "For example, if we plot the function\n",
    "\n",
    "\n",
    "<a id='equation-root-f'></a>\n",
    "$$\n",
    "f(x) = \\sin(4 (x - 1/4)) + x + x^{20} - 1 \\tag{2}\n",
    "$$\n",
    "\n",
    "with $ x \\in [0,1] $ we get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "f = lambda x: np.sin(4 * (x - 1/4)) + x + x**20 - 1\n",
    "x = np.linspace(0, 1, 100)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "ax.plot(x, f(x))\n",
    "ax.axhline(ls='--', c='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unique root is approximately 0.408.\n",
    "\n",
    "Let&rsquo;s consider some numerical techniques for finding roots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bisection\n",
    "\n",
    "\n",
    "<a id='index-6'></a>\n",
    "One of the most common algorithms for numerical root-finding is *bisection*.\n",
    "\n",
    "To understand the idea, recall the well-known game where\n",
    "\n",
    "- Player A thinks of a secret number between 1 and 100  \n",
    "- Player B asks if it&rsquo;s less than 50  \n",
    "  \n",
    "  - If yes, B asks if it&rsquo;s less than 25  \n",
    "  - If no, B asks if it&rsquo;s less than 75  \n",
    "  \n",
    "\n",
    "\n",
    "And so on.\n",
    "\n",
    "This is bisection.\n",
    "\n",
    "Here&rsquo;s a fairly simplistic implementation of the algorithm in Python.\n",
    "\n",
    "It works for all sufficiently well behaved increasing continuous functions with $ f(a) < 0 < f(b) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "\n",
    "    while upper - lower > tol:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        # === if root is between lower and middle === #\n",
    "        if f(middle) > 0:\n",
    "            lower, upper = lower, middle\n",
    "        # === if root is between middle and upper  === #\n",
    "        else:\n",
    "            lower, upper = middle, upper\n",
    "\n",
    "    return 0.5 * (upper + lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, SciPy provides its own bisection function, which we now test using the function $ f $ defined in [(2)](#equation-root-f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import bisect\n",
    "\n",
    "bisect(f, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Newton-Raphson Method\n",
    "\n",
    "\n",
    "<a id='index-8'></a>\n",
    "Another very common root-finding algorithm is the [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method).\n",
    "\n",
    "In SciPy this algorithm is implemented by `scipy.optimize.newton`.\n",
    "\n",
    "Unlike bisection, the Newton-Raphson method uses local slope information.\n",
    "\n",
    "This is a double-edged sword:\n",
    "\n",
    "- When the function is well-behaved, the Newton-Raphson method is faster than bisection.  \n",
    "- When the function is less well-behaved, the Newton-Raphson might fail.  \n",
    "\n",
    "\n",
    "Let&rsquo;s investigate this using the same function $ f $, first looking at potential instability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import newton\n",
    "\n",
    "newton(f, 0.2)   # Start the search at initial condition x = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "newton(f, 0.7)   # Start the search at x = 0.7 instead"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second initial condition leads to failure of convergence.\n",
    "\n",
    "On the other hand, using IPython&rsquo;s `timeit` magic, we see that `newton` can be much faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%timeit bisect(f, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%timeit newton(f, 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hybrid Methods\n",
    "\n",
    "So far we have seen that the Newton-Raphson method is fast but not robust.\n",
    "\n",
    "This bisection algorithm is robust but relatively slow.\n",
    "\n",
    "This illustrates a general principle\n",
    "\n",
    "- If you have specific knowledge about your function, you might be able to exploit it to generate efficiency.  \n",
    "- If not, then the algorithm choice involves a trade-off between the speed of convergence and robustness.  \n",
    "\n",
    "\n",
    "In practice, most default algorithms for root-finding, optimization and fixed points use *hybrid* methods.\n",
    "\n",
    "These methods typically combine a fast method with a robust method in the following manner:\n",
    "\n",
    "1. Attempt to use a fast method  \n",
    "1. Check diagnostics  \n",
    "1. If diagnostics are bad, then switch to a more robust algorithm  \n",
    "\n",
    "\n",
    "In `scipy.optimize`, the function `brentq` is such a hybrid method and a good default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "brentq(f, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%timeit brentq(f, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the correct solution is found and the speed is almost the same as `newton`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Root-Finding\n",
    "\n",
    "\n",
    "<a id='index-9'></a>\n",
    "Use `scipy.optimize.fsolve`, a wrapper for a hybrid method in MINPACK.\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed Points\n",
    "\n",
    "\n",
    "<a id='index-10'></a>\n",
    "SciPy has a function for finding (scalar) fixed points too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import fixed_point\n",
    "\n",
    "fixed_point(lambda x: x**2, 10.0)  # 10.0 is an initial guess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don&rsquo;t get good results, you can always switch back to the `brentq` root finder, since\n",
    "the fixed point of a function $ f $ is the root of $ g(x) := x - f(x) $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "\n",
    "<a id='index-12'></a>\n",
    "Most numerical packages provide only functions for *minimization*.\n",
    "\n",
    "Maximization can be performed by recalling that the maximizer of a function $ f $ on domain $ D $ is\n",
    "the minimizer of $ -f $ on $ D $.\n",
    "\n",
    "Minimization is closely related to root-finding: For smooth functions, interior optima correspond to roots of the first derivative.\n",
    "\n",
    "The speed/robustness trade-off described above is present with numerical optimization too.\n",
    "\n",
    "Unless you have some prior information you can exploit, it&rsquo;s usually best to use hybrid methods.\n",
    "\n",
    "For constrained, univariate (i.e., scalar) minimization, a good hybrid option is `fminbound`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import fminbound\n",
    "\n",
    "fminbound(lambda x: x**2, -1, 2)  # Search in [-1, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Optimization\n",
    "\n",
    "\n",
    "<a id='index-13'></a>\n",
    "Multivariate local optimizers include `minimize`, `fmin`, `fmin_powell`, `fmin_cg`, `fmin_bfgs`, and `fmin_ncg`.\n",
    "\n",
    "Constrained multivariate local optimizers include `fmin_l_bfgs_b`, `fmin_tnc`, `fmin_cobyla`.\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/optimize.html) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration\n",
    "\n",
    "\n",
    "<a id='index-15'></a>\n",
    "Most numerical integration methods work by computing the integral of an approximating polynomial.\n",
    "\n",
    "The resulting error depends on how well the polynomial fits the integrand, which in turn depends on how &ldquo;regular&rdquo; the integrand is.\n",
    "\n",
    "In SciPy, the relevant module for numerical integration is `scipy.integrate`.\n",
    "\n",
    "A good default for univariate integration is `quad`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "integral, error = quad(lambda x: x**2, 0, 1)\n",
    "integral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, `quad` is an interface to a very standard numerical integration routine in the Fortran library QUADPACK.\n",
    "\n",
    "It uses [Clenshaw-Curtis quadrature](https://en.wikipedia.org/wiki/Clenshaw-Curtis_quadrature),  based on expansion in terms of Chebychev polynomials.\n",
    "\n",
    "There are other options for univariate integration&mdash;a useful one is `fixed_quad`, which is fast and hence works well inside `for` loops.\n",
    "\n",
    "There are also functions for multivariate integration.\n",
    "\n",
    "See the [documentation](http://docs.scipy.org/doc/scipy/reference/integrate.html) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Algebra\n",
    "\n",
    "\n",
    "<a id='index-17'></a>\n",
    "We saw that NumPy provides a module for linear algebra called `linalg`.\n",
    "\n",
    "SciPy also provides a module for linear algebra with the same name.\n",
    "\n",
    "The latter is not an exact superset of the former, but overall it has more functionality.\n",
    "\n",
    "We leave you to investigate the [set of available routines](http://docs.scipy.org/doc/scipy/reference/linalg.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='sp-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Previously we discussed the concept of [recursive function calls](https://lectures.quantecon.org/python_advanced_features.html#recursive-functions).\n",
    "\n",
    "Write a recursive implementation of the bisection function described above, which we repeat here for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "\n",
    "    while upper - lower > tol:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        # === if root is between lower and middle === #\n",
    "        if f(middle) > 0:\n",
    "            lower, upper = lower, middle\n",
    "        # === if root is between middle and upper  === #\n",
    "        else:\n",
    "            lower, upper = middle, upper\n",
    "\n",
    "    return 0.5 * (upper + lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test it on the function `f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1` discussed above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Here&rsquo;s a reasonable solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root-finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "    if upper - lower < tol:\n",
    "        return 0.5 * (upper + lower)\n",
    "    else:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        print(f'Current mid point = {middle}')\n",
    "        if f(middle) > 0:   # Implies root is between lower and middle\n",
    "            return bisect(f, lower, middle)\n",
    "        else:               # Implies root is between middle and upper\n",
    "            return bisect(f, middle, upper)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1\n",
    "bisect(f, 0, 1)"
   ]
  }
 ],
 "metadata": {
  "filename": "scipy.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "SciPy"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}