doc/Interpolations.jl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation with `Interpolations.jl`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Interpolations.jl` implements [(cardinal) B-splines](http://en.wikipedia.org/wiki/B-spline#Cardinal_B-spline), i.e. interpolating polynomial functions which are fit to the data and to each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start with a simple 1D example - interpolating $f(x) = sin(x)$ over one period, using 10 grid points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = 10\n",
    "f(x) = sin(2pi/(nx-1) * (x-1))\n",
    "\n",
    "xcoarse = 1:nx\n",
    "ycoarse = f.(xcoarse);\n",
    "\n",
    "xfine = 1:0.1:xcoarse[end]\n",
    "yfine = f.(xfine);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interpolation object is basically an array that supports evaluation with any real numbers.\n",
    "\n",
    "To create an interpolation object, simply call the constructor `interpolate`, providing the array and some configuration for the interpolation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Interpolations\n",
    "\n",
    "yitp = interpolate(ycoarse, BSpline(Linear()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple of noteworthy points here. First, note the extra argument to `interpolate`, `BSpline(Linear())`. This determines the behavior of the interpolating function inside the domain. `Linear` tells us that it will be a linear interpolation (duh...); in other words, for `yitp(x)` interpolation will be performed like this:\n",
    "\n",
    "    begin\n",
    "        ix = ifloor(x)\n",
    "        dx = x - ix\n",
    "        (1-dx)*ycoarse[ix] + dx*ycoarse[ix+1]\n",
    "    end\n",
    "\n",
    "A B-spline interpolation is basically a piecewise polynomial function, where each grid cell is associated with its own polynomial. The coefficients are chosen such that the interpolating function passes through all the data points, while the transitions between polynomials are as smooth as possible given their degree. For example, each grid cell in a linear interpolation (such as `yitp`) is associated with a straight line, and the interpolation is continuous (with discontinuous derivative) over the entire domain.\n",
    "\n",
    "Let's take a look at how the different interpolation degrees behave:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yitp_const = interpolate(ycoarse, BSpline(Constant()))\n",
    "yconst = [yitp_const(x) for x in xfine]\n",
    "\n",
    "yitp_linear = interpolate(ycoarse, BSpline(Linear()))\n",
    "ylinear = [yitp_linear(x) for x in xfine]\n",
    "\n",
    "yitp_quadratic = interpolate(ycoarse, BSpline(Quadratic(Line(OnCell()))))\n",
    "yquadratic = [yitp_quadratic(x) for x in xfine];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that when creating the quadratic interpolation, we had to give another input parameter to the interpolation type: the `Line(OnCell())` argument which specifies the boundary condition. The first term (`Line`) specifies the behavior of the boundary condition, the second (`OnCell`) specifies where it applies. In this case `yitp_quadratic` will become linear for `x < 0.5` or `x > 10.5`. The alternative, `OnGrid`, would apply the boundary condition for `x < 1` or `x > 10`.\n",
    "\n",
    "All interpolations of quadratic degree or higher require a prefiltering step, which entails solving a tridiagonal system of equations (details can be found for example in [this paper](http://dx.doi.org/10.1109/42.875199)), in order to make the interpolating function pass through the data points. `Interpolations.jl` takes care of solving this system for you, but in order to close the system a boundary condition is required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's see the implemented functionality in action!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Gadfly\n",
    "\n",
    "plot(\n",
    "    layer(x=xcoarse,y=ycoarse,Geom.point),\n",
    "    layer(x=xfine,y=yconst,Geom.line,Theme(default_color=colorant\"magenta\")),\n",
    "    layer(x=xfine,y=ylinear,Geom.line,Theme(default_color=colorant\"red\")),\n",
    "    layer(x=xfine,y=yquadratic,Geom.line,Theme(default_color=colorant\"green\"))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Interpolations.jl` supports interpolation in arbitrary dimensions - just give it an `Array{T,N}` where elements `t::T` support operations like `a::Real * t + b::Real * t` and it will figure it out. However, since it starts getting difficult to visualize higher-dimensional functions, we'll stick with an example in 2D, namely the function $g(x,y) = x^2 \\cdot \\sin(y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = 1:5\n",
    "ys = 1:8\n",
    "g = Float64[x^2 * sin(y) for x in xs, y in ys]\n",
    "\n",
    "gitp_quad2d = interpolate(g, BSpline(Quadratic(Line(OnCell()))))\n",
    "\n",
    "display(plot(x=xs,y=ys,z=g,Geom.contour))\n",
    "display(plot(x=1:.1:5, y=1:.1:8, z=[gitp_quad2d(x,y) for x in 1:.1:5, y in 1:.1:8], Geom.contour))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple, huh? =)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extrapolation\n",
    "\n",
    "The behavior outside the domain is not carefully specified by `interpolate`, and most likely you should avoid using such objects outside the original domain. If you do want to control behavior outside the domain, you will likely be interested in `extrapolate`:\n",
    "\n",
    "    yextr = extrapolate(yitp, Flat)\n",
    "\n",
    "The `Flat` parameter to `extrapolate` signifies that we want the value to remain constant (flat) outside of the domain. Other options are...(how the interpolation object will behave when an index outside of the domain $[1, 10]$ is given; with `ExtrapError`, a `BoundsError` is thrown, while with `ExtrapNaN` (the only other implemented extrapolation behavior at the moment), `NaN` is returned.)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}