session3_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       ".container {\n",
       "    width: 95% !important;\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #888;\n",
       "    color: #f6f6f6;\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width to about 80 chars */\n",
       "    width: 1275px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1400px;\n",
       "    margin: auto;\n",
       "    padding-left: 1em;\n",
       "    padding-right: 1em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 3px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them together */\n",
       "div.cell.code_cell {\n",
       "    background-color: #f4f3e0; \n",
       "    border-radius: 10px; /* rounded borders */\n",
       "    padding: 1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Arvo' sans-serif;\n",
       "    line-height: 130%;\n",
       "    font-size: 125%;\n",
       "    width:1100px;\n",
       "    margin-left: 0em;\n",
       "    margin-right: auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 32pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(12,85,97);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 700;\n",
       "    font-size: 24pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(95,92,72);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(95,92,72);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"PT Mono\";\n",
       "    font-size: 110%;\n",
       "}\n",
       "\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "css_file = './stylesheets/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div><img src=\"https://raw.githubusercontent.com/ACM-SIAM-KAUST-chapter/python_tutorial_2017/gh-pages/images/acmchapterlogo.png\" width=\"260px\" style=\"float: left; margin-right: 3em\">\n",
    "<img src=\"https://raw.githubusercontent.com/ACM-SIAM-KAUST-chapter/python_tutorial_2017/gh-pages/images/siamchapterlogo.png\" width=\"260px\" style=\"float: right;\"></div>\n",
    "<h1 align=\"center\">Plotting and Visualization</h1>\n",
    "\n",
    "<p><br style=\"clear: both;\"/></p><p><br style=\"clear: both;\"/></p><p><br style=\"clear: both;\"/></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"color:#5f5c48\">\n",
    "This tutorial is prepared by ACM/SIAM Student Chapter of King Abdullah University of Science and Technology (KAUST).\n",
    "<br/><br/>\n",
    "\n",
    "Parts of this tutorial re-use <a href=\"https://github.com/jrjohansson/scientific-python-lectures/\">Scientific Python lectures</a> by <a href=\"https://github.com/jrjohansson\">Robert Johansson</a> linsensed under <a href=\"https://creativecommons.org/licenses/by/3.0/\">Creative Commons Attribution 3.0</a>.\n",
    "<br/><br/>\n",
    "\n",
    "Presented on 28 February, 2017 by Yiannis Hadjimichael.\n",
    "<br/><br/>\n",
    "\n",
    "**Prerequisites:** *Basic knowledge of Python, Numpy package.*\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matplotlib\n",
    "\n",
    "[`Matplotlib`](http://matplotlib.org/) is the basic package of Python for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first line above ensures that our plots appear in the browser, instead in a separate window. The second line imports  the `pyplot` package\n",
    "\n",
    "Let's see a very simple example. Suppose we want to plot the function $\\sin(\\exp(x))$ on the interval $x\\in(0,\\pi)$. We'll want to use the `numpy` functions of the sine and exponential functions. (Note: the math module operates only on scalars so we *should* use `numpy` in such cases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(0,np.pi,1000)\n",
    "f = np.sin(np.exp(x))\n",
    "g = np.cos(2*x)\n",
    "plt.plot(x,f)\n",
    "plt.plot(x,g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways to format a plot in Python. Here we just mention only the basic functionality. It's straightforward to use them and you will generally only get into its advanced functionality if you are customizing the look of a plot for a publication. `Matplotlib` documentation is very good, have a look [here](http://matplotlib.sourceforge.net/users/pyplot_tutorial.html). You can also read [this tutorial](http://www.loria.fr/~rougier/teaching/matplotlib) for an easy start. More references are listed at the end of this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code does the basic formatting for the plot above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5), dpi=150) # changing figure's shape\n",
    "plt.xlim(x.min()-0.1, x.max()+0.1) # adjusting horizonatal axis limits\n",
    "plt.ylim(f.min()-0.1, f.max()+0.1) # adjusting vertical axis limits\n",
    "plt.xticks([0, np.pi/2, np.pi],[r'$-\\pi$', r'$-\\pi/2$', r'$0$', r'$+\\pi/2$', r'$+\\pi$'],fontsize=14) # setting ticks\n",
    "plt.yticks([-1, 0, +1],[r'$-1$', r'$0$', r'$+1$'],fontsize=14) # setting ticks\n",
    "\n",
    "plt.plot(x, f, color=\"blue\", linewidth=2.5, linestyle=\"-\", label=\"$\\sin(e^x)$\") # changing color and thickness\n",
    "plt.plot(x, g, color=\"red\",  linewidth=2.5, linestyle=\"-\",  label=\"$\\cos(2x)$\") # changing color and thickness\n",
    "plt.legend(loc='lower left',prop={'size':16}) # placing legend on bottom left\n",
    "plt.xlabel('$x$',fontsize=16) # horizontal axis name\n",
    "plt.ylabel('test functions',fontsize=16) # vertical axis name\n",
    "plt.title('Sample plot',fontsize=18) # title\n",
    "plt.grid(True) # enabling grid\n",
    "plt.savefig('trig_functions.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**One thing to be aware of is that by default, the *figure* object is cleared at the end of each cell.**\n",
    "\n",
    "You will need to issue all plotting commands for a single figure in a single cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.savefig('trig_functions.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.savefig?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go on and open the *pdf* file. It is empty! Now uncomment the save line in the block where the figure is plotted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3D-plotting\n",
    "\n",
    "The following example is adapted from [this tutorial](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-2-Numpy.ipynb). Let's define the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha = 0.7\n",
    "phi_ext = 2*np.pi*0.5\n",
    "\n",
    "def flux_qubit_potential(phi_m, phi_p):\n",
    "    return 2 + alpha - 2*np.cos(phi_p)*np.cos(phi_m) - alpha*np.cos(phi_ext - 2*phi_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a mesh and evaluate the function over that mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "phi_m = np.linspace(0, 2*np.pi, 100)\n",
    "phi_p = np.linspace(0, 2*np.pi, 100)\n",
    "X,Y = np.meshgrid(phi_p, phi_m)\n",
    "Z = flux_qubit_potential(X, Y).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d.axes3d import Axes3D # imports 3D plotting\n",
    "from matplotlib import cm # module for color pattern\n",
    "fig = plt.figure(figsize=(20,8))\n",
    "\n",
    "# `ax` is a 3D-aware axis instance because of the projection='3d' keyword argument to add_subplot\n",
    "ax = fig.add_subplot(1, 2, 1, projection='3d')\n",
    "\n",
    "p = ax.plot_surface(X, Y, Z, rstride=4, cstride=4, linewidth=0)\n",
    "\n",
    "# surface_plot with color grading and color bar\n",
    "ax = fig.add_subplot(1, 2, 2, projection='3d')\n",
    "p = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)\n",
    "cb = fig.colorbar(p, shrink=0.8) # enables colorbar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Better control over plotting\n",
    "\n",
    "So far, we used **`matplotlib.pyplot`** module directly to produce plots. Sometimes we need more control over the the plots position, axes, insets, etc. We might also need to create different figures, or add different axes, put insets, etc. The following commands demonstrate this functionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Let's generate some synthetic data\n",
    "data = np.random.uniform(size=(100,2))\n",
    "x, y = data[:,0], data[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Let's create a figure object and an axes object\n",
    "# Don't forget to create figures and axes in the cell you are plotting!\n",
    "fig = plt.figure(figsize=(6,6))\n",
    "ax1 = fig.add_axes([0.10,0.20,0.78,0.80])  # left, bottom, width, height\n",
    "\n",
    "from matplotlib import rcParams\n",
    "rcParams.update({'font.size': 16}) # change fontsize of all elements (ticks, labels, title)\n",
    "\n",
    "ax1.scatter(x, y, marker='o')\n",
    "\n",
    "# Exercise: Specify a title and the axes labels\n",
    "# Hint: Use ax1. and press tap or print(dir(ax1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Insets through multiple axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "axes1 = fig.add_axes([0.1, 0.1, 0.8, 0.8]) # main axes\n",
    "axes2 = fig.add_axes([0.2, 0.5, 0.4, 0.3]) # inset axes\n",
    "\n",
    "# main figure\n",
    "axes1.plot(x, y, 'r.')\n",
    "axes1.set_xlabel('x')\n",
    "axes1.set_ylabel('y')\n",
    "axes1.set_title('title')\n",
    "\n",
    "# insert\n",
    "axes2.plot(y, x, 'g.')\n",
    "axes2.set_xlabel('y')\n",
    "axes2.set_ylabel('x')\n",
    "\n",
    "# Exerice: add a title to the inset with a font size of 14px\n",
    "axes2.set_title('zoom',fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if you want to enlarge some region of the figure and show it in the inset? Let's how we can do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(0,np.pi,1000)\n",
    "f = np.sin(8*np.exp(x))\n",
    "g = np.cos(2*x)\n",
    "\n",
    "# How to enlarge a region on a plot?\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax  = fig.add_axes([0.15,0.11,0.83,0.83])\n",
    "ax.set_xlim([0, 5])\n",
    "ax.plot(x,f)\n",
    "ax.plot(x,g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset\n",
    "zoomed_inset_axes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 10))\n",
    "ax  = fig.add_axes([0.15,0.11,0.83,0.83])\n",
    "ax.set_xlim([0, 5])\n",
    "ax.plot(x,f,linewidth=0.5)\n",
    "ax.plot(x,g)\n",
    "\n",
    "# The inset axes\n",
    "axins = zoomed_inset_axes(ax, 8, loc=1)\n",
    "axins.plot(x, f, linewidth=1)\n",
    "axins.plot(x, g, linewidth=2)\n",
    "axins.set_xlim([2.85, 3.0])\n",
    "axins.set_ylim([0.9, 1.0])\n",
    "\n",
    "# Mark inset\n",
    "mark_inset(ax, axins, loc1=3, loc2=2, fc=\"none\", ec=\"0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**`Matplotlib`** is a very rich library, and it has the functionality for literally any kind of plotting. For more references, you can check out a very detailed [User Guide](http://matplotlib.org/users/) and the links at the bottom of this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seaborn: Beautiful and easy plot styling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A recently developed [Seaborn](http://stanford.edu/~mwaskom/software/seaborn/tutorial.html) package provides many different options of how to control the plot aesthetics. Let's look at some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** In case you don't have Seaborn installed, then run in a terminal\n",
    "\n",
    "`$ conda install Seaborn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A very quick enhancement by just importing seaborn\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "\n",
    "x = np.linspace(0,np.pi,1000)\n",
    "f = np.sin(np.exp(x))\n",
    "g = np.cos(2*x)\n",
    "plt.plot(x,f)\n",
    "plt.plot(x,g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set_style?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example the following gives a horizontal boxplot of the distribution of tips:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tips = sns.load_dataset(\"tips\")\n",
    "print(tips)\n",
    "sns.set() # reset setttings\n",
    "ax = sns.boxplot(x=tips[\"total_bill\"])\n",
    "ax.tick_params(labelsize=16)\n",
    "ax.xaxis.label.set(size=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# What are the styles and how to override them?\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.axes_style()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color paletes available in `seaborn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "current_palette = sns.color_palette()\n",
    "sns.palplot(current_palette)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot([0, 1], [0, 1], sns.xkcd_rgb[\"pale red\"], lw=3)\n",
    "plt.plot([0, 1], [0, 2], sns.xkcd_rgb[\"medium green\"], lw=3)\n",
    "plt.plot([0, 1], [0, 3], sns.xkcd_rgb[\"denim blue\"], lw=3);\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many more fancy things to do with Python visualization package. \n",
    "For example we can create animations. \n",
    "\n",
    "The code below uses the solution of the [Woodward-Colella blast wave problem](http://clawpack.github.io/doc/pyclaw/gallery/woodward_colella_blast.html) solved by the [Clawpack package](http://clawpack.github.io/index.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import the necessary packages. There several ways to make animations in Python; here we describe a simple one.\n",
    "Note that in the cell below the command *%matplotlib nbagg* in needed so as the animation is created within the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib nbagg\n",
    "from matplotlib import animation\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First make some functions to read data (solution and time) from files.\n",
    "The data files are located at *data/session3/woodward_colella_blast* starting form the directory of this notebook.\n",
    "They contain the solution (density, momentum, energy) of the Euler equations in 1D solved up to time $t = 0.38$ over the interval $[0,1]$.\n",
    "\n",
    "You may first open a data file to get familiar with the way the data is stored.\n",
    "This will help you understand how to read the data we want to create an animation of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def solution(filename, path='data/session3/woodward_colella_blast'):\n",
    "    \"\"\"\n",
    "    This function reads the solution stored in an \n",
    "    ascii file written by pyclaw.\n",
    "    \"\"\"\n",
    "    import sys\n",
    "    \n",
    "    # Concatenate path and file name\n",
    "    pathfilename = path + \"/\" + filename\n",
    "\n",
    "    try:\n",
    "        f = open(pathfilename,\"r\")\n",
    "    except IOError:\n",
    "        raise IOError(\"No data file was found.\")\n",
    "\n",
    "    # Read file header\n",
    "    # The information contained in the first two lines are not used.\n",
    "    unused = f.readline()  # grid_number\n",
    "    unused = f.readline()  # AMR_level\n",
    "\n",
    "    # Read mx, my, xlow, ylow, dx and dy\n",
    "    line = f.readline()\n",
    "    sline = line.split()\n",
    "    mx = int(sline[0])\n",
    "\n",
    "    line = f.readline()\n",
    "    sline = line.split()\n",
    "    xlower = float(sline[0])\n",
    "\n",
    "    line = f.readline()\n",
    "    sline = line.split()\n",
    "    dx = float(sline[0])\n",
    "\n",
    "    # Grid:\n",
    "    xupper = xlower + mx * dx\n",
    "    xc = np.linspace(xlower+dx/2.,xupper-dx/2.,mx)  \n",
    "    \n",
    "    # Read solution\n",
    "    # Define arrays of conserved variables\n",
    "    q = np.zeros((3,mx))\n",
    "\n",
    "    line = f.readline()\n",
    "    for j in range(mx):\n",
    "        line = f.readline()\n",
    "        q[:,j] = np.array(line.split())\n",
    "    \n",
    "    return q, xc\n",
    "\n",
    "def time(filename, path='data/session3/woodward_colella_blast'):\n",
    "    \"\"\"\n",
    "    This function reads the time stored in an \n",
    "    ascii file written by pyclaw.\n",
    "    \"\"\"\n",
    "    import sys\n",
    "    \n",
    "    # Concatenate path and file name\n",
    "    pathfilename = path + \"/\" + filename\n",
    "\n",
    "    try:\n",
    "        f = open(pathfilename,\"r\")\n",
    "    except IOError:\n",
    "        raise IOError(\"No data file was found.\")\n",
    "\n",
    "    # Read time\n",
    "    line = f.readline()\n",
    "    sline = line.split()\n",
    "    t = float(sline[0])\n",
    "\n",
    "    return t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell gives the functions used for the animation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fplot(i):\n",
    "    index = \"%02d\" % i\n",
    "    q, xc = solution('fort.q00' + index)\n",
    "    t = time('fort.t00' + index)\n",
    "    line.set_data(xc, q[0,:])\n",
    "    ax.set_title('Density at t='+str(t),fontsize=26)\n",
    "    return line\n",
    "\n",
    "def init():\n",
    "    line.set_data([], [])\n",
    "    return line,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create the animation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib nbagg\n",
    "from matplotlib import rcParams\n",
    "rcParams.update({'font.size': 16}) # change fontsize of all elements (ticks, labels, title)\n",
    "\n",
    "# initial data\n",
    "_, xc = solution('fort.q0000')\n",
    "fig = plt.figure(figsize=[12,8])\n",
    "ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 12))\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize =16)\n",
    "line, = ax.plot([], [], lw=2)\n",
    "\n",
    "# animation \n",
    "animation.FuncAnimation(fig, fplot, frames=11, init_func=init, blit=True, interval=500, repeat=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References:\n",
    "\n",
    "- [Numpy for Matlab users](http://www.scipy.org/NumPy_for_Matlab_Users)\n",
    "- [Hans Petter Langtangen book](http://link.springer.com/book/10.1007/978-3-642-02475-7/page/1)\n",
    "\n",
    "\n",
    "#### Matplotlib\n",
    "\n",
    "- [Nicolas P. Rougier matplotlib tutorial](http://www.loria.fr/~rougier/teaching/matplotlib/)\n",
    "- [Pyplot tutorial](http://matplotlib.org/users/pyplot_tutorial.html)\n",
    "- [Matplotlib gallery](http://matplotlib.org/gallery.html)\n",
    "\n",
    "#### Other visualization packages\n",
    "- [Seaborn styling package tutorial](http://stanford.edu/~mwaskom/software/seaborn/tutorial.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Copyright 2017, ACM/SIAM Student Chapter of King Abdullah University of Science and Technology*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}