{
 "metadata": {
  "name": "",
  "signature": "sha256:1bf62b287fc669ca2449ceca00a88ac6d77e0444cfaad2ed5e2d1e563e688062"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Version:  Python 2.7"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import this"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Scientific Python for Engineers"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Plotting with MatPlotLib"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook discusses the main features of [Numerical Python](http://www.numpy.org/) (`numpy`), [Scientific Python](http://www.scipy.org/) (`scipy`), and their relatives in an engineering context.\n",
      "\n",
      "Although any installation of [IPython](http://ipython.org/) will work with a version of this notebook, we recommend that you download and install the [Enthought Canopy Distribution](https://www.enthought.com/products/canopy/) of Python, which is free for academic users.  To launch the notebook, open a command terminal, type `ipython notebook tutorial.ipynb`, and press Return.\n",
      "\n",
      "A few notes on this tutorial:\n",
      "\n",
      "- Check the top left corner of this page for the Python version this notebook is supposed to work with.  Enthought Canopy uses Python 2.7.x, so if you are using that version, you will need the 2.7 version of this notebook.  There is also a 3.3 version, which requires some small adjustments if you are familiar with Python 2.7.  Python 3 represents the future of Python and is a better platform for those developing new codes (rather than developing existing scripts).  Most major modules are available in both flavors now.\n",
      "\n",
      "- Early on, we will incidentally utilize features of packages not yet introduced (such as `matplotlib`).  Rest assured that the major features of these modules will be explicitly discussed at some point if it is not clear now.\n",
      "\n",
      "- Code blocks starting with `$` are intended to be run on the command line, not executed as Python code."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Standard Header"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we will be utilizing a number of packages with reasonably long names, we will adopt the _de facto_ standard module abbreviations in the following header.  We also ensure that our [division behavior is sensible](http://www.python.org/dev/peps/pep-0238/) by importing from `__future__`:  _i.e._, promotion to `double` will occur from `int` or `long` data types involving division:  `1/2 == 0.5`.  Although this is the default in Python 3, it is a trivial way to help this notebook work in Python 2 if that's what you are using."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import numpy as np\n",
      "import scipy as sp\n",
      "import matplotlib as mpl\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "#IPython magic command for inline plotting\n",
      "%matplotlib inline\n",
      "#a better plot shape for IPython\n",
      "mpl.rcParams['figure.figsize']=[15,3]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Plotting in 2D:  `matplotlib`"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus far you have seen the use of a cryptic (or not-so-cryptic) `plt.plot` command many times.  Its syntax is quite similar to the MATLAB `plot` command; indeed, you can even load it with `from matplotlib.pyplot import plot` to save yourself much typing in the future.\n",
      "\n",
      "`plot` displays coordinate pairs in a number of formats.  It is currently the workhorse of 2D plotting in Python, and the MATLAB-like syntax will make basic manipulation accessible to many engineers.\n",
      "\n",
      "There are a number of alternative display commands as well which tweak the output:  `fill` and `errorbar`, for instance.  The MatPlotLib developers maintain a [gallery](http://matplotlib.org/gallery.html) of examples with full source code."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 1, 10001)\n",
      "y = np.cos(np.pi/x) * np.exp(-2*x)\n",
      "\n",
      "plt.plot(x, y)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Exercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* Plot the following equations over the domain $y \\in \\left[-1, 2\\right]$.\n",
      "    * $y = f(x) = x^2 \\exp(-x)$\n",
      "    * $y = f(x) = \\log x$\n",
      "    * $y = f(x) = 1 + x^x + 3 x^4$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 10, 21)\n",
      "y = 10*np.exp(-x)\n",
      "yerr = np.random.rand(len(x)) * 2*np.exp(-x/2)\n",
      "\n",
      "plt.errorbar(x, y, yerr=yerr)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 2*np.pi, 1001)\n",
      "y = np.sin(2*x**2/np.pi)\n",
      "\n",
      "plt.fill(x, y)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(0.1, 4, 0.5)\n",
      "y = np.exp(-x)\n",
      "\n",
      "y_err = 0.1 + 0.2 * x\n",
      "x_err_lower = 0.4 * y_err\n",
      "x_err_upper = y_err\n",
      "x_err = [x_err_lower, x_err_upper]\n",
      "\n",
      "mpl.rcParams['figure.figsize']=[8,8]\n",
      "\n",
      "fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True)\n",
      "ax0.errorbar(x, y, yerr=y_err, fmt='-o')\n",
      "ax0.set_title('variable, symmetric error')\n",
      "\n",
      "ax1.errorbar(x, y, xerr=x_err, fmt='o')\n",
      "ax1.set_title('variable, asymmetric error')\n",
      "ax1.set_yscale('log')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.random import normal\n",
      "n = (1000,1)\n",
      "x = normal(size=n)\n",
      "\n",
      "avg = np.mean(x)\n",
      "std = np.std(x)\n",
      "\n",
      "x_avg  = np.ones(n)* avg\n",
      "x_stdl = np.ones(n)*(avg-std)\n",
      "x_stdh = np.ones(n)*(avg+std)\n",
      "\n",
      "plt.plot(x,'bx',x_avg,'r-',x_stdl,'r--',x_stdh,'r--')\n",
      "plt.title('%d Random Gaussian Numbers' % n[0]);\n",
      "plt.xlabel(r'$n$');\n",
      "plt.ylabel(r'$U(n)$')\n",
      "plt.show()\n",
      "\n",
      "plt.hist(x,bins=20)\n",
      "plt.title('Distribution of %d Random Gaussian Numbers' % n[0]);\n",
      "plt.xlabel(r'$U(n)$');\n",
      "plt.ylabel(r'Frequency of $U$')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def func(x, y):\n",
      "    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(2*np.pi*np.sqrt(y))\n",
      "\n",
      "# Define the basic grid coordinates.\n",
      "grid_x, grid_y = np.mgrid[0:1:250j, 0:1:250j]\n",
      "\n",
      "# Define a random subset of the grid for which we will generate data.\n",
      "pts = np.random.rand(800,2)\n",
      "vals = func(pts[:,0], pts[:,1])\n",
      "\n",
      "# Load the methods and generate a grid for each approach.\n",
      "from scipy.interpolate import griddata\n",
      "grid_z0 = griddata(pts, vals, (grid_x, grid_y), method='nearest')\n",
      "grid_z1 = griddata(pts, vals, (grid_x, grid_y), method='linear')\n",
      "grid_z2 = griddata(pts, vals, (grid_x, grid_y), method='cubic')\n",
      "\n",
      "from matplotlib import cm\n",
      "plt.subplot(221)\n",
      "plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n",
      "plt.plot(pts[:,0], pts[:,1], 'k.', ms=1)\n",
      "plt.title('Original')\n",
      "\n",
      "plt.subplot(222)\n",
      "plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n",
      "plt.title('Nearest')\n",
      "\n",
      "plt.subplot(223)\n",
      "plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n",
      "plt.title('Linear')\n",
      "\n",
      "plt.subplot(224)\n",
      "plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower', cmap=cm.PiYG)\n",
      "plt.title('Cubic')\n",
      "\n",
      "plt.gcf().set_size_inches(12, 12)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Play around with some other functions in the cell above's `func`."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Color"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Named MPL colors\n",
      "for c in mpl.colors.cnames:\n",
      "    print(c, mpl.colors.cnames[c])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Professional Plotting"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using standard commands, you can easily design publication-quality output.  (This is one place where an interactive shell like IPython really shines.)  You can save the figures as high-resolution PNGs (or other file types).  You can even use $\\LaTeX$ markup in raw Python strings to yield mathematical formulae in labels, titles, and legends."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mpl.rcParams['figure.figsize']=[15,3]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.polynomial import Chebyshev as T\n",
      "from numpy.polynomial import Polynomial as P\n",
      "x = np.linspace(-1, 1, 100)\n",
      "for i in range(6):\n",
      "    ax = plt.plot(x, T.basis(i)(x), lw=2, label=r'$T_%d$'%i)\n",
      "    print(r'T_%d(x) = %s'%(i, T.basis(i).convert(kind=P)))\n",
      "plt.legend(loc=\"upper left\")\n",
      "plt.savefig('chebyshev.png', dpi=120, transparent=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.polynomial import Laguerre as L\n",
      "x = np.linspace(0, 4, 100)\n",
      "for i in range(6):\n",
      "    ax = plt.plot(x, L.basis(i)(x), lw=2, label=r'$L_%d$'%i)\n",
      "    print(r'L_%d(x) = %s'%(i, L.basis(i).convert(kind=P)))\n",
      "plt.legend(loc=\"lower left\",ncol=3)\n",
      "plt.savefig('laguerre.png', dpi=120, transparent=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Experiment with the following $\\LaTeX$ labels in the `label` argument in the following code block.  (They won't describe the figure, but whatever...)\n",
      "    - `r'A_i^j (x)$'`\n",
      "    - `r'$(1 - x_0) \\cdot x^{x^x}$'`\n",
      "    - `r'$\\frac{x + y}{x^y}$'`\n",
      "    - `r'$\\int_0^\\infty \\exp(-x) dx$'`\n",
      "    - `r'$\\sum_{n=0}^{10} \\frac{x}{x-5}$'`\n",
      "    - `r'$\\pi = %f$'%np.pi`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(0, 2*np.pi, 101)\n",
      "plt.plot(x, np.sin(x), lw=2, label=r'$L_i$')\n",
      "plt.legend(loc=\"lower left\", prop={'size':24})\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To get plots really looking sharp, you need to adjust the labels, font sizes, typefaces, and other plot properties for readability."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.special import jv, yv, iv, kv, jn_zeros\n",
      "x = np.linspace(0, 6, 201)\n",
      "y = np.zeros((8, 201))\n",
      "y[0,:] = jv(0, x)\n",
      "y[1,:] = jv(1, x)\n",
      "y[2,:] = yv(0, x)\n",
      "y[3,:] = yv(1, x)\n",
      "y[4,:] = iv(0, x)\n",
      "y[5,:] = iv(1, x)\n",
      "y[6,:] = kv(0, x)\n",
      "y[7,:] = kv(1, x)\n",
      "\n",
      "mpl.rcParams['figure.figsize']=[15,6]\n",
      "fig = plt.figure()\n",
      "ax = fig.add_subplot(111)\n",
      "ax.plot(x, y[0], 'r-', lw=2, label=r'$J_0(x)$')\n",
      "ax.plot(x, y[1], 'r--', lw=2, label=r'$J_1(x)$')\n",
      "ax.plot(x, y[2], 'b-', lw=2, label=r'$Y_0(x)$')\n",
      "ax.plot(x, y[3], 'b--', lw=2, label=r'$Y_1(x)$')\n",
      "ax.plot(x, y[4], 'g-', lw=2, label=r'$I_0(x)$')\n",
      "ax.plot(x, y[5], 'g--', lw=2, label=r'$I_1(x)$')\n",
      "ax.plot(x, y[6], 'y-', lw=2, label=r'$K_0(x)$')\n",
      "ax.plot(x, y[7], 'y--', lw=2, label=r'$K_1(x)$')\n",
      "\n",
      "ax.set_title(r'Examples of Zeroth- and First-Order Bessel Functions', fontsize=24, family='serif')\n",
      "ax.set_ylabel(r'$f(x)$', fontsize=18)\n",
      "ax.set_xlabel(r'$x$', fontsize=18)\n",
      "ax.set_ylim((-1, 2))\n",
      "\n",
      "from matplotlib import font_manager as fm\n",
      "ticks_font = fm.FontProperties(family='serif', size=14)\n",
      "for label in ax.get_xticklabels(): label.set_fontproperties(ticks_font)\n",
      "for label in ax.get_yticklabels(): label.set_fontproperties(ticks_font)\n",
      "\n",
      "ax.legend(loc='best', prop={'size':18}, ncol=4) #you often have to oversize LaTeX code to make it look right in MPL\n",
      "\n",
      "plt.savefig('bessel.png', dpi=120, transparent=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# the only difference here is that we set the lines to grayscale for publication purposes\n",
      "fig = plt.figure()\n",
      "ax = fig.add_subplot(111)\n",
      "ax.plot(x, y[0], 'k-', lw=2, label=r'$J_0(x)$')\n",
      "ax.plot(x, y[1], 'k--', lw=2, label=r'$J_1(x)$')\n",
      "ax.plot(x, y[2], 'k.', lw=2, label=r'$Y_0(x)$')\n",
      "ax.plot(x, y[3], 'k:', lw=2, label=r'$Y_1(x)$')\n",
      "ax.plot(x, y[4], '-', color=(0.75,0.75,0.75), lw=2, label=r'$I_0(x)$')\n",
      "ax.plot(x, y[5], '--', color=(0.75,0.75,0.75), lw=2, label=r'$I_1(x)$')\n",
      "ax.plot(x, y[6], '.', color=(0.75,0.75,0.75), lw=2, label=r'$K_0(x)$')\n",
      "ax.plot(x, y[7], ':', color=(0.75,0.75,0.75), lw=2, label=r'$K_1(x)$')\n",
      "\n",
      "ax.set_title(r'Examples of Zeroth- and First-Order Bessel Functions', fontsize=24, family='serif')\n",
      "ax.set_ylabel(r'$f(x)$', fontsize=18)\n",
      "ax.set_xlabel(r'$x$', fontsize=18)\n",
      "ax.set_ylim((-1, 2))\n",
      "\n",
      "from matplotlib import font_manager as fm\n",
      "ticks_font = fm.FontProperties(family='serif', size=14)\n",
      "for label in ax.get_xticklabels(): label.set_fontproperties(ticks_font)\n",
      "for label in ax.get_yticklabels(): label.set_fontproperties(ticks_font)\n",
      "\n",
      "ax.legend(loc='best', prop={'size':18}, ncol=4) #you often have to oversize LaTeX code to make it look right in MPL\n",
      "\n",
      "plt.savefig('bessel.png', dpi=120, transparent=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Animating 3D Plots"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The modes of vibration of a drum are the orthogonal basis functions which uniquely describe any vibration on a struck drumhead."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mpl_toolkits.mplot3d import Axes3D\n",
      "from scipy.special import jn, jn_zeros\n",
      "mpl.rcParams['figure.figsize']=[12,8]\n",
      "\n",
      "def drumhead_height(n, k, distance, angle, t):\n",
      "    nth_zero = jn_zeros(n, k)\n",
      "    return np.cos(t)*np.cos(n*angle)*jn(n, distance*nth_zero[-1])\n",
      "\n",
      "# Define polar and cartesian coordinates for the drum.\n",
      "theta  = np.r_[0:2*np.pi:50j]\n",
      "radius = np.r_[0:1:50j]\n",
      "x = np.array([r*np.cos(theta) for r in radius])\n",
      "y = np.array([r*np.sin(theta) for r in radius])\n",
      "\n",
      "radial_nodes  = 2\n",
      "zeros = 3\n",
      "\n",
      "# Define the base plot.\n",
      "fig = plt.figure()\n",
      "\n",
      "# Display the desired angular nodes.\n",
      "ax = fig.add_subplot(1,1,1,projection='3d')\n",
      "z = np.array([drumhead_height(radial_nodes, zeros, r, theta, 0.5) for r in radius])\n",
      "ax.plot_surface(x,y,z,rstride=1,cstride=1,cmap=mpl.cm.jet)\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's display several at once."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mpl_toolkits.mplot3d import Axes3D\n",
      "from scipy.special import jn, jn_zeros\n",
      "import subprocess\n",
      "\n",
      "def drumhead_height(n, k, distance, angle, t):\n",
      "    nth_zero = jn_zeros(n, k)\n",
      "    return np.cos(t)*np.cos(n*angle)*jn(n, distance*nth_zero[-1])\n",
      "\n",
      "# Define polar and cartesian coordinates for the drum.\n",
      "theta  = np.r_[0:2*np.pi:50j]\n",
      "radius = np.r_[0:1:50j]\n",
      "x = np.array([r*np.cos(theta) for r in radius])\n",
      "y = np.array([r*np.sin(theta) for r in radius])\n",
      "\n",
      "radial_nodes = 3\n",
      "zeros = 2\n",
      "\n",
      "# Define the base plot.\n",
      "fig = plt.figure(num=None,figsize=(16,16),dpi=120,facecolor='w',edgecolor='k')\n",
      "ax = list()\n",
      "\n",
      "# Loop over the desired angular nodes.\n",
      "cnt = 0\n",
      "pixcnt = 0\n",
      "\n",
      "cnt = 0\n",
      "pixcnt += 1\n",
      "for i in np.r_[0:radial_nodes+1:1]:\n",
      "    for j in np.r_[1:zeros+1:1]:\n",
      "        cnt += 1;\n",
      "        ax.append(fig.add_subplot(radial_nodes+1,zeros,cnt,projection='3d'))\n",
      "        z = np.array([drumhead_height(i, j, r, theta, t) for r in radius])\n",
      "        ax[-1].set_xlabel('R@%d,A@%d' % (i,j))\n",
      "        ax[-1].plot_surface(x,y,z,rstride=1,cstride=1,cmap=mpl.cm.Accent,linewidth=0,vmin=-1,vmax=1)\n",
      "        ax[-1].set_zlim(-1,1)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And write them to disk:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mpl_toolkits.mplot3d import Axes3D\n",
      "from scipy.special import jn, jn_zeros\n",
      "import subprocess\n",
      "\n",
      "def drumhead_height(n, k, distance, angle, t):\n",
      "    nth_zero = jn_zeros(n, k)\n",
      "    return np.cos(t)*np.cos(n*angle)*jn(n, distance*nth_zero[-1])\n",
      "\n",
      "# Define polar and cartesian coordinates for the drum.\n",
      "theta  = np.r_[0:2*np.pi:50j]\n",
      "radius = np.r_[0:1:50j]\n",
      "x = np.array([r*np.cos(theta) for r in radius])\n",
      "y = np.array([r*np.sin(theta) for r in radius])\n",
      "\n",
      "radial_nodes = 3\n",
      "zeros = 2\n",
      "\n",
      "# Define the base plot.\n",
      "fig = plt.figure(num=None,figsize=(16,16),dpi=120,facecolor='w',edgecolor='k')\n",
      "ax = list()\n",
      "\n",
      "# Loop over the desired angular nodes.\n",
      "cnt = 0\n",
      "pixcnt = 0\n",
      "#plt.ion()\n",
      "for t in np.r_[0:2*np.pi:40j]:\n",
      "    cnt = 0\n",
      "    pixcnt += 1\n",
      "    for i in np.r_[0:radial_nodes+1:1]:\n",
      "        for j in np.r_[1:zeros+1:1]:\n",
      "            cnt += 1;\n",
      "            ax.append(fig.add_subplot(radial_nodes+1,zeros,cnt,projection='3d'))\n",
      "            z = np.array([drumhead_height(i, j, r, theta, t) for r in radius])\n",
      "            ax[-1].set_xlabel('R@%d,A@%d' % (i,j))\n",
      "            ax[-1].plot_surface(x,y,z,rstride=1,cstride=1,cmap=mpl.cm.Accent,linewidth=0,vmin=-1,vmax=1)\n",
      "            ax[-1].set_zlim(-1,1)\n",
      "    #plt.show()\n",
      "    plt.savefig('./drum-modes-%d.png' % pixcnt, format='png')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, let's collect them as an animated GIF."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Collate pictures to an animated GIF.\n",
      "import os,string\n",
      "cwd = os.getcwd()\n",
      "cmd = 'cd %s; ls drum-modes*.png | sort -k1.12n'%cwd\n",
      "png_files = os.popen(cmd)\n",
      "png_files_list = string.join(png_files.readlines()).replace('\\n',' ')\n",
      "os.popen('convert -delay 10 -loop 1 %s ./drum-animate.gif'%png_files_list)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are also ways to [generate animated plots](http://matplotlib.org/api/animation_api.html) in MatPlotLib, but that's beyond the scope of this tutorial."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Statistical Plotting (`seaborn`)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MatPlotLib is a stable basis for plotting in Python, but for domain applications or specialized needs it can fall short.  [Seaborn](https://github.com/mwaskom/seaborn) is designed to create high-level statistical graphics with a high degree of customizability."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import seaborn as sns\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "np.random.seed(sum(map(ord, \"palettes\")))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Palettes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One criticism which has been levelled against MatPlotLib is that the default color palettes leave much to be desired.  Seaborn has a number of built-in basic palettes available which provide a pleasant aesthetic."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "current_palette = sns.color_palette()\n",
      "print(current_palette)\n",
      "sns.palplot(current_palette)\n",
      "\n",
      "gammas = sns.load_dataset(\"gammas\")\n",
      "plt.figure(figsize=(12, 6))\n",
      "sns.tsplot(gammas, \"timepoint\", \"subject\", \"ROI\", \"BOLD signal\", color=current_palette);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use built-in Seaborn palettes or pull from MatPlotLib [colormaps](http://matplotlib.org/examples/color/colormaps_reference.html?highlight=cmap)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Built-in Seaborn variations\n",
      "sns.palplot(sns.color_palette(\"pastel\", 8))\n",
      "sns.palplot(sns.color_palette(\"deep\", 8))\n",
      "sns.palplot(sns.color_palette(\"dark\", 8))\n",
      "sns.palplot(sns.color_palette(\"colorblind\", 8))\n",
      "\n",
      "# Selections from MPL colormaps.\n",
      "sns.palplot(sns.color_palette(\"coolwarm\", 8))\n",
      "sns.palplot(sns.color_palette(\"gist_earth\", 8))\n",
      "sns.palplot(sns.color_palette(\"Paired\", 10))\n",
      "sns.palplot(sns.color_palette(\"RdPu_r\", 8))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`blend_palette()` allows you to interpolate between an arbitrary number of specified MatPlotLib colors."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sns.palplot(sns.blend_palette([\"indigo\", \"#fff3f4\",\"royalblue\"], 12))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are [a few more ways](http://stanford.edu/~mwaskom/software/seaborn/tutorial/color_palettes.html) to generate and use palettes in both MatPlotLib and Seaborn, but this should give you a good feel for the flexibility and ease-of-use of this way of working with color."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, palettes are invoked using the `color` `kwarg` in Seaborn:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(12, 6))\n",
      "sns.tsplot(gammas, \"timepoint\", \"subject\", \"ROI\", \"BOLD signal\", color=sns.color_palette(\"RdPu_r\", 3));"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Best Practices for Plotting"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Output plots as vector images (`eps` (preferred) or `svg`) as well as raster (`png` (preferred), `gif`, or `jpg`).  This allows you to generate new versions with more detail if necessary.  Also have a clearly documented section of code and data for each image in order to reproduce images as well if necessary."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy.polynomial import Chebyshev as T\n",
      "from numpy.polynomial import Polynomial as P\n",
      "x = np.linspace(-1, 1, 100)\n",
      "for i in range(6):\n",
      "    ax = plt.plot(x, T.basis(i)(x), lw=2, label=r'$T_%d$'%i)\n",
      "plt.legend(loc=\"upper left\")\n",
      "plt.savefig('chebyshev.png', dpi=120, transparent=True)\n",
      "plt.savefig('chebyshev.eps', dpi=120, transparent=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Consider utilizing a CUBEHELIX palette in documents which may be printed or accessed in black and white<sup>[[citation](http://adsabs.harvard.edu/abs/2011BASI...39..289G)]; [[tutorial](http://www.mrao.cam.ac.uk/~dag/CUBEHELIX/)]</sup>.  The CUBEHELIX algorithm generates color palettes which continuously increase (or decrease) in intensity regardless of the hue, thus making them acceptable to print or view in monochrome and for individuals with any form of colorblindness.  They also stand up well in desaturated environments like lousy projectors in overbright rooms.  The basic CUBEHELIX palette is available in MatPlotLib and thus any of its derivatives; the module `cubehelix` is also available to generate more variations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.special import j0\n",
      "def func(x, y):\n",
      "    return j0(8*np.pi*x) * np.sin(2*np.pi*np.sqrt(y))\n",
      "\n",
      "# Define the basic grid coordinates.\n",
      "grid_x, grid_y = np.mgrid[0:0.2:250j, 0:0.4:250j]\n",
      "\n",
      "from matplotlib import cm\n",
      "plt.subplot('121')\n",
      "plt.imshow(func(grid_x, grid_y).T, extent=(0,0.2,0,0.4), origin='lower', cmap=cm.cubehelix_r)\n",
      "plt.title(r'$z(x, y) = J_{0}(8\\pi x) \\sin(2 \\pi \\sqrt{y})$')\n",
      "plt.xlabel(r'$x$')\n",
      "plt.ylabel(r'$y$')\n",
      "\n",
      "import cubehelix\n",
      "plt.subplot('122')\n",
      "plt.imshow(func(grid_x, grid_y).T, extent=(0,0.2,0,0.4), origin='lower', cmap=cubehelix.cmap(reverse=True, start=0.3, rot=0.5))\n",
      "plt.title(r'$z(x, y) = J_{0}(8\\pi x) \\sin(2 \\pi \\sqrt{y})$')\n",
      "plt.xlabel(r'$x$')\n",
      "plt.ylabel(r'$y$')\n",
      "\n",
      "\n",
      "plt.gcf().set_size_inches(4, 4)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "References"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Langtangen, Hans Petter.  _Python Scripting for Computational Science_, 3ed.  Berlin\u2013Heidelberg:  Springer\u2013Verlag, 2009.\n",
      "- Lugo, Michael.  [On propagation of errors](http://gottwurfelt.com/2012/03/26/on-propagation-of-errors/).  26 March 2012.\n",
      "- Warren, Russell.  [A Brief Intro to Profiling in Python](https://speakerdeck.com/rwarren/a-brief-intro-to-profiling-in-python).  Ottawa Python Authors Group, 28 February 2013.\n",
      "\n",
      "### Seaborn\n",
      "- Waskon, Michael.  [Choosing color palettes](http://stanford.edu/~mwaskom/software/seaborn/tutorial/color_palettes.html)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Neal Davis developed these materials for [Computational Science and Engineering](http://cse.illinois.edu/) at the University of Illinois at Urbana\u2013Champaign.  This content is available under a Creative Commons Attribution 3.0 Unported License."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}