{
 "metadata": {
  "name": "",
  "signature": "sha256:41e2949d98d4571553301d9f65fdad8824953d6530ed605a118b65e427d39e09"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Introduction to Scipy: Interpolation and Integration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "In this lecture, we will look at two other common sub-packages of Scipy: [scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html) and [scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Interpolation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "The simplest interpolation routine in [scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html) is [interp1d](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.interpolate import interp1d"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "If we create a fake dataset:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "x = np.array([0., 1., 3., 4.])\n",
      "y = np.array([0., 4., 3., 2.])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "we can interpolate linearly by first creating an interpolating function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = interp1d(x, y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "and we can then interpolate to any value of x within the original bounds:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(3.3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "It is also possible to interpolate to several values at the same time:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(np.array([0.5, 1.5, 2.5, 3.5]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "If the interpolating function is called outside the original range, an error is raised:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(-1.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "You can change this behavior by telling ``interp1d`` to not give an error in this case, but to use a set value:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = interp1d(x, y, bounds_error=False, fill_value=-10.)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(-1.0)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f(np.array([-1., 1., 3., 6.]))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "By default, ``interp1d`` uses linear interpolation, but it is also possible to use e.g. cubic interpolation: "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = interp1d(x, y, kind='cubic')\n",
      "f(0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "For more information, see the documentation for [interp1d](http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d). There are also other interpolation functions available (for example for spline interpolation), which you can read up about at [scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Integration"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "The available integration functions are listed at [scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate). You will notice there are two kinds of functions - those that integrate actual Python functions, and those that integrate numerical functions defined by Numpy arrays."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "First, we can take a look at one of the functions that can integrate actual Python functions. If we define a function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def simple_function(x):\n",
      "    return 3. * x**2 + 2. * x + 1."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "we can integrate it between limits using:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import quad\n",
      "print(quad(simple_function, 1., 2.))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "source": [
      "As described in the documentation for [quad](http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad), the first value returned is the integral, and the second is the error on the integral. If we had solved the integral analytically, we would expect 11, so the result is correct."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "We can also define functions as Numpy arrays:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(1., 2., 1000)\n",
      "y = 3. * x**2 + 2. * x + 1."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "-"
      }
     },
     "source": [
      "And in this case we can use for example the [simps](http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simps.html#scipy.integrate.simps) function to integrate using Simpson's rule:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import simps\n",
      "print(simps(y, x=x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "This can be very useful in cases where one wants to integrate actual data that cannot be represented as a simple function or when the function is only available numerically."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Note that there is an issue on the [scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html#module-scipy.integrate) page - there should also be a menton of the ``trapz`` function which works similarly to ``simps`` but does trapezium integration:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import trapz\n",
      "print(trapz(y, x=x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Exercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Write a function that takes ``x``, and the parameters for a Gaussian (amplitude, displacement, width) and returns the value of the Gaussian at ``x``:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "# your solution here\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use ``quad`` to compute the integral and compare to what you would expect."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "# your solution here\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now create two arrays ``x`` and ``y`` that contain the Gaussian for fixed values ``x``, and try and compute the integral using ``simps``."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "# your solution here\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare this to what you found with ``quad`` and analytically."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}