yoavram/ggplot2.ipynb

## ggplot2.ipynb
{
 "metadata": {
  "name": "example"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "# Example of using `ggplot2` from *IPython notebook*\n## By [Yoav Ram](http://www.yoavram.com/), 31 March 2013"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Overview\n\nThe following is an example of how to use [`ggplot2`](http://ggplot2.org/) inside an [*IPython notebook*](http://ipython.org/).\n\nFor the data I will use the results of some evolutionary simulations I ran. As the main point here is to demonstrate the use of `R` and `ggplot2` in the *IPython noteook* I will not explain what the data.\n\n## Parse filenames\nFirst I need to parse the filename for the simulation parameters.\n\nThe regular expression was written using the [Python regular expression testing tool](http://www.pythonregex.com/)."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import re\nfilename_pattern = pattern = re.compile('^pop_(?P<pop>\\d+)_G_(?P<G>\\d+)_s_(?P<s>\\d\\.?\\d*)_H_(?P<H>\\d\\.?\\d*)_U_(?P<U>\\d\\.?\\d*)_beta_(?P<beta>\\d\\.?\\d*)_pi_(?P<pi>\\d\\.?\\d*)_tau_(?P<tau>\\d\\.?\\d*)_(?P<date>\\d{4}-\\w{3}-\\d{1,2})_(?P<time>\\d{2}-\\d{2}-\\d{2}-\\d{6}).(?P<extension>\\w+)$')\ndef parse_filename(fname):\n    m = pattern.match(fname)\n    if m:\n         return m.groupdict()\n    else:\n        return dict()",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Process data file\n\nNext I need to read the data from file.\n\nData files are with `.data` extension and in `JSON` format, compressed with `gzip`.\n\nYou can use the builtin `json` parser but I use the native one I found on *github* by the name [ultrajson](https://github.com/esnme/ultrajson) because it is roughly 3-4 times faster."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import ujson as json\nimport gzip\nfolder = 'output/fixation/'\n\ndef process_data_file(fname):\n    fpath = folder + fname\n    params = parse_filename(fname)\n    if not params:\n        print \"Failed parsing file name\", fpath\n        return {},[],[]\n    with gzip.open(fpath) as f:\n        data = json.load(f,precise_float=True)\n    if not data:\n        print \"Failed reading data\", fpath\n        return {},[],[]\n    data.update(params)\n    W = data.pop('W')\n    p = data.pop('p')\n    data['fname'] = fname\n    for k in ['tau',  'G',  'H',  'pop',  'beta',  'U',  'T',  'pop_size',  's',  'pi']:\n        if str == type(data[k]):\n            data[k] = eval(data[k])  \n    return data, W, p",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Process all files into a `list`, each item in the `list` is a `dict` containing the results of a single simulation:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import glob, os, time\ntic = time.clock()\nfile_list = glob.glob1(folder, '*.data')\nall_data = [None] * len(file_list)\nprint \"processing\", len(file_list), \"data files\"\nfor i,fname in enumerate(file_list) :\n    data,W,p = process_data_file(fname)\n    all_data[i] = data\ntoc = time.clock()\nprint \"processed all files in\", (toc-tic), \"seconds\"",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "processing 316 data files\nprocessed all files in"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": " 0.34 seconds\n"
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Next I create a matrix of the values I want to plot:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "df = [[data['T'],data['tau'],data['s'],data['pi']] for data in all_data]",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Plotting the data with `ggplot2`"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "I call the `rmagic` extension of *IPython notebook`*. Make sure you install `rpy2`, for example run: `pip install rpy2`."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%load_ext rmagic",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The final step is to send the `df` to *R* and plot the data using `ggplot2`. The input to `R` is defined by using the `-i` option:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%%R -i df\ndf <- as.data.frame(df)\nnames(df) <- c(\"T\",\"tau\",\"s\",\"pi\")\nlibrary(ggplot2)\np <- ggplot(df, aes(tau, T))\np <- p + \n    geom_point(alpha=I(0.3)) + \n    scale_x_log10() + scale_y_log10() + \n    facet_grid(facets=s~pi, labeller=function(variable,value) {paste0(variable,'=',as.character(value))}) +\n    labs(y=\"Adaptation time\", x=expression(tau)) \nprint(p)",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAMAAABKCk6nAAACxFBMVEUAAAACAgIDAwMEBAQGBgYH\nBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcZGRkaGhob\nGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwuLi4v\nLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBCQkJD\nQ0NERERFRUVGRkZISEhJSUlKSkpLS0tOTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZ\nWVlbW1tcXFxeXl5fX19gYGBhYWFiYmJkZGRlZWVmZmZnZ2doaGhpaWlsbGxtbW1ubm5vb29wcHBx\ncXFycnJzc3N0dHR3d3d4eHh5eXl6enp7e3t8fHx+fn5/f3+AgICBgYGCgoKDg4OFhYWGhoaHh4eI\niIiJiYmKioqLi4uNjY2Ojo6Pj4+QkJCRkZGTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJyd\nnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+w\nsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLD\nw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW\n1tbX19fY2Nja2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq\n6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9\n/f3+/v7///9GDwhaAAAZy0lEQVR4nO2dj38U5Z3HsWLRS2kvem2DaO1d6C/qoRd/oFyrzal7tBTj\n3cUfd1ZOXQJIgZ74A2mvXtVYD8MmkEQTYwineCWBogFEsiICUezBnUpICCE/SLL5ufv8Eze7s4P7\n7D7PPs/szOzOPPt584LNzn7nmf3uO88zzzzPLM8MApRmRq7fAHAWCFYcCFYcCFYcCFYcCFYcCFYc\nCFYcCFYcCFYcCFYcCFYcCFacjAUHFUCRPCCYhyJ5QDAPRfKAYB6K5OFFwU+sMH5a9oxzR3E4DyqJ\njqXXL9zpzHG8KNjg4IKLvCvYIJZE5e3Bdfc6U77HBC/z3fKDTesr9Ccd793vScGpSfyiItjw184c\nzGuC7wg2/eX6iqdKNVqDwQe8KTglibKNwfZvOXMwrwn+ZTD4TeOX37OCU5L4F60G/9CZg3lN8N8F\nm65ZX7HhDg3v1uCUJF4oDT6Oc3CUZYtu/N5Lnq/BKUkcLFu4qM2Zg3lNsINGk3BQcPaSgGA+EJwL\nwVlEkTycFfzB/0q8hWPikKN/Esd88KE45uOPxTFHDumPiXmc/EC8n0Qah/5PHCNTzpH/Ece8f0R/\ndEZwV5wz4S4xfeKQ/nFxTM8ZcczweXFMb7f+mJjHtETZEml0R8QxMuX0TYhjTp/VH5O8QzAEQzAb\nCLaQEQTTIcoJvvzywtSNRUVFiU9pwYWFrD3mzPlilwu7x3/QH7R/nRNcNLcodWOxBrecpBzjzPvu\nF3sUGpnGfyiK7aKV6ajgnj6dgXCfmAFhRKEmOHVrNJWEp+f6E/fQYOxx5ZVf7HFhd/2HotiD9nfO\nOf31CCEX3n1Y/B4l0iiaM5d6xzrzNLjlJOUY3+M73/1ijwuZxn+I7REt8zv6656owYW8Gpz4603V\n4EKpGlyU+EPsocjJGlyUUQ1O3aU49zXYOLp9TTRrq+kmumvOnMS93dFEd9F+zTbRXRcSjf8Q28Nr\ngtHJokOU62RBMB0CwRwgGIIhGILZQLCFjCCYDoFgDhAMwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBs\nISMIpkMgmAMEQzAEQzAbCLaQEQTTIRDMAYIhGIIhmI1agt/uTHwWeSOwdWRvZWXlkxCcDs8IDlc9\nTgk+2UAOtWqPn+yG4HR4RnAk3NZJppo2bT6nyQ0Rsmc/6a0iZLp6Imr72LFjx/t1hsL9Ys6LQ0Ym\nxTGDg+KYsTFxzNCA/hhN1NgWHhLvJ5HGIBHHyJQzPCVfToZN9K5O0rGTnKohpFWTvOMjMvEcIfuO\nRF96q7q6+t1xnYnIuJhJiZCwOGZiQhwzPS1fTjQXY1tEomyJNCaIOMauj8MoJ3PBb9Y0N9d2Nj//\nSvPobq0Gv0wiL0wbLxvtBJpoCs800THBB/aRwaBeg080ksOtpKvhwsvGUSCYwluCJxtr6j/XBUd2\nbK0bJbsOQ3B6PCRYgHEUCKaAYB4QTIdAMAcIhmAIhmA2EGwhIwimQyCYAwRDMARDMBsItpARBNMh\nEMwBgiEYgiGYDQRbyAiC6RAI5gDBEAzBEMwGgi1kBMF0iEnBmx/b/eirEBxHQcHX1JQ/vAiC4ygo\n+IbgTdofCNZRUHDpkvJ1d0FwHAUFt21sbztgq2CbF+XoG54Ux1CLcnAIhcQx/Y4tytHXHxHHyJQz\nNCVfTkzrhuBT6EUbKFiDg3cESyHYAIIhmAkEW8gIgukQCOaQr4JbtT8QHEdBwQdffLDsocoOCI6h\noOB7SpdvXH7bUgiOoaDghbHW+XoIjqGg4EXrG9sb1y2E4BgKCm4rn3/1teU7ITiGgoIbMOGfgIKC\n50FwAgoKDkJwAgoKxnVwIgoKxnVwIgoKxnVwIgoKxnVwIgoKxnVwIgzBgW1bxPu5WHCUjoPoReuw\nBQeE+7lc8HsLIFhHQcEPaNx/0TIIjqFgE/2MDgTHULCTlcFIFmNRjqmnKyv3Jgne0iJu2iDYdDm5\nWZSj9/WUGhyAYArPCGYuynH8pfqGQe3F41oxnQMx6l9pqR8QMiIOGZ0UxwydF8eMj4ljzg/pj9FE\njW1hibIl0hgk4hiZckamJI41rD9m2EQzFuX47Bg5Wqe9tLuhoeFgSOe1ltdCQiYkQqbFMWPj4pip\nKXHM+Jj+GE3T2BaRKFsijTEijpH6OMISx4q/5cwFpyzKQcjks3QTrUony+u96MwEJy/K0X6QnKxT\nU7B/TYWnr4MzE5y8KEeovqrmXLLgFnsyyr1gv3A/tQQLMI7C7EX7NUxmlGPBvrsX51sNlhPMvkzy\nnuA8bKIlBVesZTRt3hPs9bFoB2vwWtYHQ/v1gGCtBov3y0/BaoxkoQZzBHPGogP0Rgg2W457BDOv\ngz15Dhbvl5eCFelk5d1kg6xgvyKCUYN5gheX+RjHxjnYYjmuERy48dYSxrEh2GI5rhHMbqIDAdqw\nJwQzAs2noZxg9kiW9wR3sQRnkIZyglWZbGA20YGAH4KZ18Heq8HMyYaA+d9TCObhSsFoonl3dHhu\nsoHTREOwQufg1Ditic57wapM+DOHKlGDITgZ5QSzm+gM2jY3NtEZDMipJ1iR22YxVJmMYoty1LbU\n16YG1tLb3LwoB2pwAtK9aJyDlRHM7mT5SjCSpbJgX0kJNRcKwRYycqVgjcTnEGwhIzcKDpSUQPAW\n1pfPvDfZgKFKjmD2SJYnBWOygSlYnbsqIVhxweL98lKwIk00bnznCFblu0kQzBOs8mVSF+7JUltw\nBl0JCOYBwXQIBHOwuYmmnkKwhYxcKTgJCLaQEQTTIRDMAYIhGIIhmA0EW8gIgukQCOYAwRAMwWYE\n90iI6eqVyGhEHNPdI44ZHBTHnDmtPybmMS5RtkQa3RJiZMo5OyqOOR0vxxnBwBtAsOJAsOJYPgd3\nD0ucdCQ6ML0S506Zc/C5c+KYHkYna1iibIk0Tkt0JWTKOTMkcaz4W0YvOhn0oiGYCQRbyAiC6RAI\n5gDBEAzBMij2Df++sMR7xDf8LfzKogbTITmvwcZRIJgCgnlAMB0CwRwgGIIhGILZQLCFjCCYDoFg\nDhAMwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBsISMIpkMgmAMEQzAEQzAbCLaQEQTTIU4Jnj4dgWAe\n3hd86keXfrXjlpMQzMb7gn/yr2NXTa28DYLZeF/wrAFyFTl7KQSz8b7g+f+lCX7jB7Gf3+5MfCXy\nRmDryN7KysonITgdLhe8p3BJwc+u+IP2U7jqcUrwyQZyqFV7/GQ3BKfD5YJJ3+YnNnXHamy4rZNM\nNW3afE6TG9LU7ye9VVonu3pCe/HD9vb2o0M6w5EhMaPikNCUOGZ4WBwzMS6OGTmvP0YTNbZFJMqW\nSOM8EcfIlDM6LXGsEf3RjODeP0WJ/birk3TsJKdqCGnVJO/4iEw8R8i+I9GX9m7fvj04qjMWGRUz\nJg4ZnxbHhELimMlJ+bcTzcXYFpF4jxIhISKOkfo4wvLlmBD8yIyrvq1hCH6zprm5trP5+VeaR3dr\nNfhlEnlh2gg12gk00RQub6Ive/eLnzXBB/aRwaBeg080ksOtpKvhwsvGUSCYwuWCr/uUEjzZWFP/\nuS44smNr3SjZdRiC0+Nywfv/4he/0iASGEeBYAqXC775ulVrNCCYjfcFF/bLuIVgPi4XXFEHwWnw\nvuCbLyqapwHBbLwv+GMdCGbjfcEmMI4CwRSuFnzdm9fpQDAbrwvu6NvTEWUPBLPxuuCpqblTGgMF\nEMzG64K/9KUZX4pSBsFsvC6YkJ/IqIXgNLhcsAmMo0AwBQTzgGA6BII5QDAE57HgtpvmYSyaj/cF\nz111FGPRfFwsePNjux99VSz4GyEZtxDMJ1eCr6kpf3iRWPCvN0wTSYyjQDBFrgTfELxJ+yMUvPCy\nrxTjHMzFxYJLl5Svu0ssGPPBaXGx4LaN7W0HZC6TprunZPRiUQ5e0blalGND8CmJXnTv3RdfdvHd\nvajBbFxcg4N3BEslBC++r5f0lv8Mgtl4X3BB9LbZvq9AMBvvC/6r6FeAW9GL5uB9wa/NvvdX985u\ngmA2bhbcqv2R6EWf+O2q356Q8QvBbHIl+OCLD5Y9VNmB2aQYCgq+p3T5xuW3LU0vGLfNCnCx4IWx\n1vn69IJx26wAFwtetL6xvXHdwvSCcdusABcLbiuff/W15TvTC8ZtswJcLLhBbsIft82mxcWC58kJ\n1qmBYDYuFhyUE/zxg2VlZb6vQzAbFwuWvA7+m/tW3rml5A0IZuNiwVLXwYRc0jdyO+m5BYLZuFiw\n1HUwIVe9Q67rD/15HgoObNsi3s/FgqWugwn5/cxP1y5YVJqfggPC/VwsWOo6mJCRz8aman93GoLZ\nuFjw7keCVf/+2rsYyYrBOgd7vIn+0arghu/f+WcYyYqhYCdrgfb3b4PFGMmKoaDgG36/76nbxYJJ\n5JM/Hg9DMAcXC952yzdvfTP4iHAk64eFCy6fr9/4zliUY+rpysq9igr2+DlYcqjyxlWTZHLVzYSz\nKEfv64wa3GJPRuhF0yHOCJ59VvunN3rbLHNRjuMv1TcMai92tLa2Hh7WaWxpGBYSEoeMTYtjRkbF\nMRMT4pjREf0xmnJ8U0NLozgPiTRGiDhGppyQxMdhlGNC8D++qP3zwtLYz4xFOT47Ro7W0YIbNMG2\nfDK5Fry2Ubyf9wX/88Xz/37+jDvLysrYi3IQMvks3UQHtrSImzY00abLcaiJ3nIB5qIc7QfJyTpa\nsLPn4ECA+tQhmCLj/4SlxmiikxflCNVX1ZxLFuxkLzpAG4ZgCgUm/LMnuGubt6+DJQVnMOHv6Dk4\ne0201wc6JAWbn/BXpZOVJ4LNT/gHKtb6bcmILTibnSzxft4XbH7C31nB6GSlI5NedGzCf8iMYKfP\nwVkS7F9TkReC5TGOsrjMnhpc18TYmMUafPdi9QXP1vlqkQnBfmYTnSRGIqNAXZP4A3a2Bot/Ub0u\neGDgNz/e373/9qpcCF65OteC1a/BhMyJ5nxqjgnBXRWsJjrg95sU7F+5WlyDHBTsu3uxeD/vCy56\nR/vn7SvNCGZ2sswLDtSxanD2Bjryookm5D++9th/Pva135kQzL5MMt9Es2twVnvR+SCYvPNo2fJ9\nMn7Tn4P9Ph+1VeoczCqHbgmcFMxsoulfUxUEa4w2LDFXgxlNqybGpGB2Dfb5fLkc6AiYTsP1gsdb\n/uGyqx41IbjLJsHsXnRSOdkWbD4Ndwue/O/7Zn/vkj9GZPymH8nK4JNhXgdnr4lmjkWbP9O4W3Dh\nt9ccJV/vlvKbtaHKxKfZnk0K+Ep8ic+9Lrj4Gys/NCu4awuric4gIzcKVq2JjgRXXnntl/eYaqId\nF5ytOzryQbBG+J1ll3+rwoRg9mVSBhmtXM3YmGPB2knY5OW86wVrTL51r0nBdgzicgc6Ep9m/Rxs\nerzGC4IliR8kYFMTneuxaAhO5sKiHC22rEKxYtUvV6RurV2xojbhabYX5ajVSHzu5kU5nKrBdt02\nyx2qzNJAB+emOxWHKu0QbNdkQ84F00CwgW2zSVkbi4Zgc4JNzwd3MceiA/4s3lXJCDSfRr4I9pm+\ngGQOdGCyIR05FJzBJwPBdEh+CM7qje+Mw0Oww4Kz99UVXCaZEpzBEJArv3yWzyNZxlGYgn0aJjNy\no2Ccg7t4TbT5WyFy/M0GuwRX23T/g9sFZ3JPlgsnG8yn4bdpetzlgrvoD8ars0nmb5v1l9xqy+yp\n6wWbzsiVNdh8L9pXPLeYJdh0b1w5wa48B5u/GCgpmF2Q+ovqKy4uMVeOgoLd2Is2L7i44NKC1F/U\nglmzCsyVA8FcclyDi64oSt06S8NcORDMJbedrEDJ90tStxbPmkVph2AOHpgPZl4HB0pKzF5PQzAP\nTPhDMARDMBsItpARBNMhrhHcMyqRkcSHd3ZIIqMecUx/vzimhyF4RKJsiTS6Q+IYmXJ6JX5PjY/D\nGcHAG0Cw4kCw4lg/B8ucdHrFIX3DEicdiZPXwIA4pue0/piYR0jiHCyRRrdEX1GmnN4RcYzxcaAX\nnQx60RDMBIItZATBdAgEc4BgCIZgGS58wz/M/QK6qa+0D0+KY7L9Df9k8A1/C7+yqMF0SM5rsHEU\nCKaAYB4QTIdAMAcIhmAIhmA2EGwhIwimQyCYAwRDMARDMBsItpARBNMhEMwBgiEYgiGYDQRbyAiC\n6RAI5gDBEAzBEMwGgi1kBMF0CARzgGAIhmAIZqOW4Lc7E59F3ghsHdlbWVn5JASnwzOCw1WPU4JP\nNpBDrdrjJ7shOB2eERwJt3WSqaZNm89pckOE7NlPeqsIma6egOB0eEYwIbs6ScdOcqqGkFZN8o6P\nyMRzhOw7En1p7/bt24OjOmORUTFj4pDxaXFMKCSOmZyUfzvRXIxtEYn3KBESIuIYqY8jLF9O5oLf\nrGluru1sfv6V5tHdWg1+mURemI6+dFwrpnNAZygyIGZEHDI6KY4ZOi+OGR8Tx5wf0h+juRjbwhJl\nS6QxSMQxMuWMTEkca1h/zFzwgX1kMKjX4BON5HAr6Wq48LLRTqCJpvBWEz3ZWFP/uS44smNr3SjZ\ndRiC0+MhwQKMo0AwBQTzgGA6BII5QDAEQzAEs4FgCxlBMB0CwRwgGIIhGILZQLCFjCCYDoFgDhAM\nwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBsISMIpkMgmAMEQzAEQzAbFwve/NjuR1+F4DgKCr6mpvzh\nRRAcR0HBNwRv0v7YKRiLcrCLztGiHKVLytfdhRocR8Ea3Laxve0ABMdRUHBwQ/Ap9KINVBR8R7AU\ngg0gGIKZQLCFjCCYDoFgDvkquFX7A8FxFBR88MUHyx6q7IDgGAoKvqd0+cblty2F4BgKCl4Ya52v\nh+AYCgpetL6xvXHdQgiOoaDgtvL5V19bvhOCYygouAET/gkoKHgeBCegoOAgBCegqOCOgxCso6jg\n9xZAsI6Cgh/QuP+iZQ4I3tISsCUjCKZDTAp+Rsd+wYGKtX5bMsq14G3eFuxYJ8uviODAti3i/fJR\nsCo1WBMsPtWoJZixKMfU05WVe2nBXRVrlTgH551g5qIcva+n1mBVOln51kQzF+U4/lJ9w6CM4ICG\nyYxyLdjjvegMmmjGohyfHSNH67SX3qqurn53XOf1lqbxFJpWr6a2TqaGJDMZFsdMTIhjpqfly4mm\naWyLSJQtkcYEEcfY9XEY5WQuOGVRDkImn42218eOHTverzMU7k9lpUbi8/OMmCRGJsUxg4PimLEx\ncczQgP4YTdPYFh4S7yeRxiARx8iUMzwlX07mgpMX5Wg/SE7W0U1015kWRuMR8PvRRFsqJzeLcoTq\nq2rO0YJVOQfjOjiZ+EGiglMNe1Jwfl0myQpmj0X7NUxmBMF0iGsEM8/BHhS8pkK8X14KZp+DPdfJ\nQg3mCFamk7WmAoKZglWZbIBgrmAlxqLRRPMEKzLZgOtgjuAcThcm9dTj4I4OewXnrgYnX4rFwWUS\nBLME+9dUMEozfzGgnGDOXZVJn70jTbTP57TggPnxGvUEM++L9mmYzMi04ORr7TiZC06N85tPA4J5\n5Fgwc7oQgrtyKDj5BBnHTsEBCOZ9dSXps3dIMCvGVsE4B+fwqyvJn34cWwWjF53DyQZ7L5PYt+yg\nBtsnuK6JtZUuBjWYIpfzwaY/mcDK1eyBhsRyst7JgmDeOdh0J8tvXjAbW5toCLZtqJIt2E/fGZLt\n2STUYBvHotmCfWkE29pEM+eDAwHTdx6pI9hYlKO+pTZ1tYhaDXOrUNQ2NDPKWaGR8JRalKOWfs0g\nw0U5alvqU4+fdHg3L8qR5YGODDpZdU3s22+5NTh5HDGOrU20+XsH1anBxlFsEtzFFBxIc5lkr2D0\nok0JzsZQZRaug7tQg+37b5RWrjZZg23vZDEPD8E2CeZcB2fzMsmWr1hBMAfudXDiVgcF+5n3RUOw\nbYIDP/058zrYl6UazBacBARnnlGAeQ7O3kgWBPME2zTZkMFYtK0T/hDMEWzXdCF7NindOdjeyyR8\ndcVpwXWcyyT+SJbN32yw56sr1TYNzbtGsF3ThexzcDrB2RjoMJ+G36Zv8rhHsNPXwYlPs37TnenZ\nJL/vVh/zF5V6mpeCf/pzxshyjif8zf9PFP6iuUW2DJioKDinE/423TZbUjibIdhfUmL25j3lBLOn\nC33094+yLth0zSsumF2QmsesmTNnccspnDmTUZBygtnThdkbqrSpk1VSeEVRaks0S4NbzsyZLMPq\nCWZOF+a4iU5G4kxTcmtJ6taCWbMKuOXkteB090VzyPH/VVnNukwKFBXxz8FMv3kjmMYDgpUbyYJg\nGgjmAcF0CARzgGCbBZ/8g0RGPeKQznaJjLrFMe+/L47pPq0/Juax84R4P4k0Pt0hjpEp5/ge+XKc\nEWxw9tdWS9A5stWecnbuymy/35yx5fDDT9hSDPlkkz3lQLABBHOAYAr1BA812PE2CDnZZk85wUOZ\n7dc4YMvhQzb9np7aYU851gUDdwPBigPBimNdML1KaeaFxNY3taWQjEpyTxrW8kjGquDkVUozLyS2\nvqkthWRQknvSsJZHKlYFx1YptUqskNj6prYUkkFJ7knDWh6pWG+id9nRtmmFxNY3taWQjEpyTxrW\n8kjGPYKN9U2tF5JRSe5Jw1oeybhHcGx9U1sKyagk96RhLY9k3CM4tr6pLYVkVJJ70rCWRzK4DlYc\nCFYcCFacfBbcfcuXL7lkKtfvwmHyWfDG/UuOfhzJ9btwmHwW/Dn5p1y/BefJZ8GELM31G3Ce/Bb8\n4w+qc/0WnCa/Bf/bjNpcvwWnyW/B0z25fgeOk9+C8wAIVhwIVhwIVhwIVhwIVhwIVhwIVhwIVhwI\nVhwIVpz/B4oRnRiJY999AAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## License\n\nThe code is free (CC0). The data and results are currently not available for reuse."
    }
   ],
   "metadata": {}
  }
 ]
}
	{
	"metadata": {
	"name": "example"
	},
	"nbformat": 3,
	"nbformat_minor": 0,
	"worksheets": [
	{
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "# Example of using `ggplot2` from IPython notebook\n## By [Yoav Ram](http://www.yoavram.com/), 31 March 2013"
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "## Overview\n\nThe following is an example of how to use [`ggplot2`](http://ggplot2.org/) inside an [IPython notebook](http://ipython.org/).\n\nFor the data I will use the results of some evolutionary simulations I ran. As the main point here is to demonstrate the use of `R` and `ggplot2` in the IPython noteook I will not explain what the data.\n\n## Parse filenames\nFirst I need to parse the filename for the simulation parameters.\n\nThe regular expression was written using the [Python regular expression testing tool](http://www.pythonregex.com/)."
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "import re\nfilename_pattern = pattern = re.compile('^pop_(?P<pop>\\d+)_G_(?P<G>\\d+)_s_(?P<s>\\d\\.?\\d)_H_(?P<H>\\d\\.?\\d)_U_(?P<U>\\d\\.?\\d)_beta_(?P<beta>\\d\\.?\\d)_pi_(?P<pi>\\d\\.?\\d)_tau_(?P<tau>\\d\\.?\\d)_(?P<date>\\d{4}-\\w{3}-\\d{1,2})_(?P<time>\\d{2}-\\d{2}-\\d{2}-\\d{6}).(?P<extension>\\w+)$')\ndef parse_filename(fname):\n m = pattern.match(fname)\n if m:\n return m.groupdict()\n else:\n return dict()",
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 5
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "## Process data file\n\nNext I need to read the data from file.\n\nData files are with `.data` extension and in `JSON` format, compressed with `gzip`.\n\nYou can use the builtin `json` parser but I use the native one I found on github by the name [ultrajson](https://github.com/esnme/ultrajson) because it is roughly 3-4 times faster."
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "import ujson as json\nimport gzip\nfolder = 'output/fixation/'\n\ndef process_data_file(fname):\n fpath = folder + fname\n params = parse_filename(fname)\n if not params:\n print \"Failed parsing file name\", fpath\n return {},[],[]\n with gzip.open(fpath) as f:\n data = json.load(f,precise_float=True)\n if not data:\n print \"Failed reading data\", fpath\n return {},[],[]\n data.update(params)\n W = data.pop('W')\n p = data.pop('p')\n data['fname'] = fname\n for k in ['tau', 'G', 'H', 'pop', 'beta', 'U', 'T', 'pop_size', 's', 'pi']:\n if str == type(data[k]):\n data[k] = eval(data[k]) \n return data, W, p",
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 6
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "Process all files into a `list`, each item in the `list` is a `dict` containing the results of a single simulation:"
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "import glob, os, time\ntic = time.clock()\nfile_list = glob.glob1(folder, '.data')\nall_data = [None] len(file_list)\nprint \"processing\", len(file_list), \"data files\"\nfor i,fname in enumerate(file_list) :\n data,W,p = process_data_file(fname)\n all_data[i] = data\ntoc = time.clock()\nprint \"processed all files in\", (toc-tic), \"seconds\"",
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": "processing 316 data files\nprocessed all files in"
	},
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": " 0.34 seconds\n"
	}
	],
	"prompt_number": 7
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "Next I create a matrix of the values I want to plot:"
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "df = [[data['T'],data['tau'],data['s'],data['pi']] for data in all_data]",
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 8
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "## Plotting the data with `ggplot2`"
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "I call the `rmagic` extension of IPython notebook`. Make sure you install `rpy2`, for example run: `pip install rpy2`."
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "%load_ext rmagic",
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 9
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "The final step is to send the `df` to R and plot the data using `ggplot2`. The input to `R` is defined by using the `-i` option:"
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": "%%R -i df\ndf <- as.data.frame(df)\nnames(df) <- c(\"T\",\"tau\",\"s\",\"pi\")\nlibrary(ggplot2)\np <- ggplot(df, aes(tau, T))\np <- p + \n geom_point(alpha=I(0.3)) + \n scale_x_log10() + scale_y_log10() + \n facet_grid(facets=s~pi, labeller=function(variable,value) {paste0(variable,'=',as.character(value))}) +\n labs(y=\"Adaptation time\", x=expression(tau)) \nprint(p)",
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "display_data",
	"png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAMAAABKCk6nAAACxFBMVEUAAAACAgIDAwMEBAQGBgYH\nBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcZGRkaGhob\nGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwuLi4v\nLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBCQkJD\nQ0NERERFRUVGRkZISEhJSUlKSkpLS0tOTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZ\nWVlbW1tcXFxeXl5fX19gYGBhYWFiYmJkZGRlZWVmZmZnZ2doaGhpaWlsbGxtbW1ubm5vb29wcHBx\ncXFycnJzc3N0dHR3d3d4eHh5eXl6enp7e3t8fHx+fn5/f3+AgICBgYGCgoKDg4OFhYWGhoaHh4eI\niIiJiYmKioqLi4uNjY2Ojo6Pj4+QkJCRkZGTk5OUlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJyd\nnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWmpqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+w\nsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLD\nw8PExMTFxcXGxsbHx8fIyMjJycnKysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW\n1tbX19fY2Nja2trb29vc3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq\n6urr6+vs7Ozt7e3u7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9\n/f3+/v7///9GDwhaAAAZy0lEQVR4nO2dj38U5Z3HsWLRS2kvem2DaO1d6C/qoRd/oFyrzal7tBTj\n3cUfd1ZOXQJIgZ74A2mvXtVYD8MmkEQTYwineCWBogFEsiICUezBnUpICCE/SLL5ufv8Eze7s4P7\n7D7PPs/szOzOPPt584LNzn7nmf3uO88zzzzPLM8MApRmRq7fAHAWCFYcCFYcCFYcCFYcCFYcCFYc\nCFYcCFYcCFYcCFYcCFYcCFacjAUHFUCRPCCYhyJ5QDAPRfKAYB6K5OFFwU+sMH5a9oxzR3E4DyqJ\njqXXL9zpzHG8KNjg4IKLvCvYIJZE5e3Bdfc6U77HBC/z3fKDTesr9Ccd793vScGpSfyiItjw184c\nzGuC7wg2/eX6iqdKNVqDwQe8KTglibKNwfZvOXMwrwn+ZTD4TeOX37OCU5L4F60G/9CZg3lN8N8F\nm65ZX7HhDg3v1uCUJF4oDT6Oc3CUZYtu/N5Lnq/BKUkcLFu4qM2Zg3lNsINGk3BQcPaSgGA+EJwL\nwVlEkTycFfzB/0q8hWPikKN/Esd88KE45uOPxTFHDumPiXmc/EC8n0Qah/5PHCNTzpH/Ece8f0R/\ndEZwV5wz4S4xfeKQ/nFxTM8ZcczweXFMb7f+mJjHtETZEml0R8QxMuX0TYhjTp/VH5O8QzAEQzAb\nCLaQEQTTIcoJvvzywtSNRUVFiU9pwYWFrD3mzPlilwu7x3/QH7R/nRNcNLcodWOxBrecpBzjzPvu\nF3sUGpnGfyiK7aKV6ajgnj6dgXCfmAFhRKEmOHVrNJWEp+f6E/fQYOxx5ZVf7HFhd/2HotiD9nfO\nOf31CCEX3n1Y/B4l0iiaM5d6xzrzNLjlJOUY3+M73/1ijwuZxn+I7REt8zv6656owYW8Gpz4603V\n4EKpGlyU+EPsocjJGlyUUQ1O3aU49zXYOLp9TTRrq+kmumvOnMS93dFEd9F+zTbRXRcSjf8Q28Nr\ngtHJokOU62RBMB0CwRwgGIIhGILZQLCFjCCYDoFgDhAMwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBs\nISMIpkMgmAMEQzAEQzAbCLaQEQTTIRDMAYIhGIIhmI1agt/uTHwWeSOwdWRvZWXlkxCcDs8IDlc9\nTgk+2UAOtWqPn+yG4HR4RnAk3NZJppo2bT6nyQ0Rsmc/6a0iZLp6Imr72LFjx/t1hsL9Ys6LQ0Ym\nxTGDg+KYsTFxzNCA/hhN1NgWHhLvJ5HGIBHHyJQzPCVfToZN9K5O0rGTnKohpFWTvOMjMvEcIfuO\nRF96q7q6+t1xnYnIuJhJiZCwOGZiQhwzPS1fTjQXY1tEomyJNCaIOMauj8MoJ3PBb9Y0N9d2Nj//\nSvPobq0Gv0wiL0wbLxvtBJpoCs800THBB/aRwaBeg080ksOtpKvhwsvGUSCYwluCJxtr6j/XBUd2\nbK0bJbsOQ3B6PCRYgHEUCKaAYB4QTIdAMAcIhmAIhmA2EGwhIwimQyCYAwRDMARDMBsItpARBNMh\nEMwBgiEYgiGYDQRbyAiC6RAI5gDBEAzBEMwGgi1kBMF0iEnBmx/b/eirEBxHQcHX1JQ/vAiC4ygo\n+IbgTdofCNZRUHDpkvJ1d0FwHAUFt21sbztgq2CbF+XoG54Ux1CLcnAIhcQx/Y4tytHXHxHHyJQz\nNCVfTkzrhuBT6EUbKFiDg3cESyHYAIIhmAkEW8gIgukQCOaQr4JbtT8QHEdBwQdffLDsocoOCI6h\noOB7SpdvXH7bUgiOoaDghbHW+XoIjqGg4EXrG9sb1y2E4BgKCm4rn3/1teU7ITiGgoIbMOGfgIKC\n50FwAgoKDkJwAgoKxnVwIgoKxnVwIgoKxnVwIgoKxnVwIgoKxnVwIgzBgW1bxPu5WHCUjoPoReuw\nBQeE+7lc8HsLIFhHQcEPaNx/0TIIjqFgE/2MDgTHULCTlcFIFmNRjqmnKyv3Jgne0iJu2iDYdDm5\nWZSj9/WUGhyAYArPCGYuynH8pfqGQe3F41oxnQMx6l9pqR8QMiIOGZ0UxwydF8eMj4ljzg/pj9FE\njW1hibIl0hgk4hiZckamJI41rD9m2EQzFuX47Bg5Wqe9tLuhoeFgSOe1ltdCQiYkQqbFMWPj4pip\nKXHM+Jj+GE3T2BaRKFsijTEijpH6OMISx4q/5cwFpyzKQcjks3QTrUony+u96MwEJy/K0X6QnKxT\nU7B/TYWnr4MzE5y8KEeovqrmXLLgFnsyyr1gv3A/tQQLMI7C7EX7NUxmlGPBvrsX51sNlhPMvkzy\nnuA8bKIlBVesZTRt3hPs9bFoB2vwWtYHQ/v1gGCtBov3y0/BaoxkoQZzBHPGogP0Rgg2W457BDOv\ngz15Dhbvl5eCFelk5d1kg6xgvyKCUYN5gheX+RjHxjnYYjmuERy48dYSxrEh2GI5rhHMbqIDAdqw\nJwQzAs2noZxg9kiW9wR3sQRnkIZyglWZbGA20YGAH4KZ18Heq8HMyYaA+d9TCObhSsFoonl3dHhu\nsoHTREOwQufg1Ditic57wapM+DOHKlGDITgZ5QSzm+gM2jY3NtEZDMipJ1iR22YxVJmMYoty1LbU\n16YG1tLb3LwoB2pwAtK9aJyDlRHM7mT5SjCSpbJgX0kJNRcKwRYycqVgjcTnEGwhIzcKDpSUQPAW\n1pfPvDfZgKFKjmD2SJYnBWOygSlYnbsqIVhxweL98lKwIk00bnznCFblu0kQzBOs8mVSF+7JUltw\nBl0JCOYBwXQIBHOwuYmmnkKwhYxcKTgJCLaQEQTTIRDMAYIhGIIhmA0EW8gIgukQCOYAwRAMwWYE\n90iI6eqVyGhEHNPdI44ZHBTHnDmtPybmMS5RtkQa3RJiZMo5OyqOOR0vxxnBwBtAsOJAsOJYPgd3\nD0ucdCQ6ML0S506Zc/C5c+KYHkYna1iibIk0Tkt0JWTKOTMkcaz4W0YvOhn0oiGYCQRbyAiC6RAI\n5gDBEAzBMij2Df++sMR7xDf8LfzKogbTITmvwcZRIJgCgnlAMB0CwRwgGIIhGILZQLCFjCCYDoFg\nDhAMwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBsISMIpkMgmAMEQzAEQzAbCLaQEQTTIU4Jnj4dgWAe\n3hd86keXfrXjlpMQzMb7gn/yr2NXTa28DYLZeF/wrAFyFTl7KQSz8b7g+f+lCX7jB7Gf3+5MfCXy\nRmDryN7KysonITgdLhe8p3BJwc+u+IP2U7jqcUrwyQZyqFV7/GQ3BKfD5YJJ3+YnNnXHamy4rZNM\nNW3afE6TG9LU7ye9VVonu3pCe/HD9vb2o0M6w5EhMaPikNCUOGZ4WBwzMS6OGTmvP0YTNbZFJMqW\nSOM8EcfIlDM6LXGsEf3RjODeP0WJ/birk3TsJKdqCGnVJO/4iEw8R8i+I9GX9m7fvj04qjMWGRUz\nJg4ZnxbHhELimMlJ+bcTzcXYFpF4jxIhISKOkfo4wvLlmBD8yIyrvq1hCH6zprm5trP5+VeaR3dr\nNfhlEnlh2gg12gk00RQub6Ive/eLnzXBB/aRwaBeg080ksOtpKvhwsvGUSCYwuWCr/uUEjzZWFP/\nuS44smNr3SjZdRiC0+Nywfv/4he/0iASGEeBYAqXC775ulVrNCCYjfcFF/bLuIVgPi4XXFEHwWnw\nvuCbLyqapwHBbLwv+GMdCGbjfcEmMI4CwRSuFnzdm9fpQDAbrwvu6NvTEWUPBLPxuuCpqblTGgMF\nEMzG64K/9KUZX4pSBsFsvC6YkJ/IqIXgNLhcsAmMo0AwBQTzgGA6BII5QDAE57HgtpvmYSyaj/cF\nz111FGPRfFwsePNjux99VSz4GyEZtxDMJ1eCr6kpf3iRWPCvN0wTSYyjQDBFrgTfELxJ+yMUvPCy\nrxTjHMzFxYJLl5Svu0ssGPPBaXGx4LaN7W0HZC6TprunZPRiUQ5e0blalGND8CmJXnTv3RdfdvHd\nvajBbFxcg4N3BEslBC++r5f0lv8Mgtl4X3BB9LbZvq9AMBvvC/6r6FeAW9GL5uB9wa/NvvdX985u\ngmA2bhbcqv2R6EWf+O2q356Q8QvBbHIl+OCLD5Y9VNmB2aQYCgq+p3T5xuW3LU0vGLfNCnCx4IWx\n1vn69IJx26wAFwtetL6xvXHdwvSCcdusABcLbiuff/W15TvTC8ZtswJcLLhBbsIft82mxcWC58kJ\n1qmBYDYuFhyUE/zxg2VlZb6vQzAbFwuWvA7+m/tW3rml5A0IZuNiwVLXwYRc0jdyO+m5BYLZuFiw\n1HUwIVe9Q67rD/15HgoObNsi3s/FgqWugwn5/cxP1y5YVJqfggPC/VwsWOo6mJCRz8aman93GoLZ\nuFjw7keCVf/+2rsYyYrBOgd7vIn+0arghu/f+WcYyYqhYCdrgfb3b4PFGMmKoaDgG36/76nbxYJJ\n5JM/Hg9DMAcXC952yzdvfTP4iHAk64eFCy6fr9/4zliUY+rpysq9igr2+DlYcqjyxlWTZHLVzYSz\nKEfv64wa3GJPRuhF0yHOCJ59VvunN3rbLHNRjuMv1TcMai92tLa2Hh7WaWxpGBYSEoeMTYtjRkbF\nMRMT4pjREf0xmnJ8U0NLozgPiTRGiDhGppyQxMdhlGNC8D++qP3zwtLYz4xFOT47Ro7W0YIbNMG2\nfDK5Fry2Ubyf9wX/88Xz/37+jDvLysrYi3IQMvks3UQHtrSImzY00abLcaiJ3nIB5qIc7QfJyTpa\nsLPn4ECA+tQhmCLj/4SlxmiikxflCNVX1ZxLFuxkLzpAG4ZgCgUm/LMnuGubt6+DJQVnMOHv6Dk4\ne0201wc6JAWbn/BXpZOVJ4LNT/gHKtb6bcmILTibnSzxft4XbH7C31nB6GSlI5NedGzCf8iMYKfP\nwVkS7F9TkReC5TGOsrjMnhpc18TYmMUafPdi9QXP1vlqkQnBfmYTnSRGIqNAXZP4A3a2Bot/Ub0u\neGDgNz/e373/9qpcCF65OteC1a/BhMyJ5nxqjgnBXRWsJjrg95sU7F+5WlyDHBTsu3uxeD/vCy56\nR/vn7SvNCGZ2sswLDtSxanD2Bjryookm5D++9th/Pva135kQzL5MMt9Es2twVnvR+SCYvPNo2fJ9\nMn7Tn4P9Ph+1VeoczCqHbgmcFMxsoulfUxUEa4w2LDFXgxlNqybGpGB2Dfb5fLkc6AiYTsP1gsdb\n/uGyqx41IbjLJsHsXnRSOdkWbD4Ndwue/O/7Zn/vkj9GZPymH8nK4JNhXgdnr4lmjkWbP9O4W3Dh\nt9ccJV/vlvKbtaHKxKfZnk0K+Ep8ic+9Lrj4Gys/NCu4awuric4gIzcKVq2JjgRXXnntl/eYaqId\nF5ytOzryQbBG+J1ll3+rwoRg9mVSBhmtXM3YmGPB2knY5OW86wVrTL51r0nBdgzicgc6Ep9m/Rxs\nerzGC4IliR8kYFMTneuxaAhO5sKiHC22rEKxYtUvV6RurV2xojbhabYX5ajVSHzu5kU5nKrBdt02\nyx2qzNJAB+emOxWHKu0QbNdkQ84F00CwgW2zSVkbi4Zgc4JNzwd3MceiA/4s3lXJCDSfRr4I9pm+\ngGQOdGCyIR05FJzBJwPBdEh+CM7qje+Mw0Oww4Kz99UVXCaZEpzBEJArv3yWzyNZxlGYgn0aJjNy\no2Ccg7t4TbT5WyFy/M0GuwRX23T/g9sFZ3JPlgsnG8yn4bdpetzlgrvoD8ars0nmb5v1l9xqy+yp\n6wWbzsiVNdh8L9pXPLeYJdh0b1w5wa48B5u/GCgpmF2Q+ovqKy4uMVeOgoLd2Is2L7i44NKC1F/U\nglmzCsyVA8FcclyDi64oSt06S8NcORDMJbedrEDJ90tStxbPmkVph2AOHpgPZl4HB0pKzF5PQzAP\nTPhDMARDMBsItpARBNMhrhHcMyqRkcSHd3ZIIqMecUx/vzimhyF4RKJsiTS6Q+IYmXJ6JX5PjY/D\nGcHAG0Cw4kCw4lg/B8ucdHrFIX3DEicdiZPXwIA4pue0/piYR0jiHCyRRrdEX1GmnN4RcYzxcaAX\nnQx60RDMBIItZATBdAgEc4BgCIZgGS58wz/M/QK6qa+0D0+KY7L9Df9k8A1/C7+yqMF0SM5rsHEU\nCKaAYB4QTIdAMAcIhmAIhmA2EGwhIwimQyCYAwRDMARDMBsItpARBNMhEMwBgiEYgiGYDQRbyAiC\n6RAI5gDBEAzBEMwGgi1kBMF0CARzgGAIhmAIZqOW4Lc7E59F3ghsHdlbWVn5JASnwzOCw1WPU4JP\nNpBDrdrjJ7shOB2eERwJt3WSqaZNm89pckOE7NlPeqsIma6egOB0eEYwIbs6ScdOcqqGkFZN8o6P\nyMRzhOw7En1p7/bt24OjOmORUTFj4pDxaXFMKCSOmZyUfzvRXIxtEYn3KBESIuIYqY8jLF9O5oLf\nrGluru1sfv6V5tHdWg1+mURemI6+dFwrpnNAZygyIGZEHDI6KY4ZOi+OGR8Tx5wf0h+juRjbwhJl\nS6QxSMQxMuWMTEkca1h/zFzwgX1kMKjX4BON5HAr6Wq48LLRTqCJpvBWEz3ZWFP/uS44smNr3SjZ\ndRiC0+MhwQKMo0AwBQTzgGA6BII5QDAEQzAEs4FgCxlBMB0CwRwgGIIhGILZQLCFjCCYDoFgDhAM\nwRAMwWwg2EJGEEyHQDAHCIZgCIZgNhBsISMIpkMgmAMEQzAEQzAbFwve/NjuR1+F4DgKCr6mpvzh\nRRAcR0HBNwRv0v7YKRiLcrCLztGiHKVLytfdhRocR8Ea3Laxve0ABMdRUHBwQ/Ap9KINVBR8R7AU\ngg0gGIKZQLCFjCCYDoFgDvkquFX7A8FxFBR88MUHyx6q7IDgGAoKvqd0+cblty2F4BgKCl4Ya52v\nh+AYCgpetL6xvXHdQgiOoaDgtvL5V19bvhOCYygouAET/gkoKHgeBCegoOAgBCegqOCOgxCso6jg\n9xZAsI6Cgh/QuP+iZQ4I3tISsCUjCKZDTAp+Rsd+wYGKtX5bMsq14G3eFuxYJ8uviODAti3i/fJR\nsCo1WBMsPtWoJZixKMfU05WVe2nBXRVrlTgH551g5qIcva+n1mBVOln51kQzF+U4/lJ9w6CM4ICG\nyYxyLdjjvegMmmjGohyfHSNH67SX3qqurn53XOf1lqbxFJpWr6a2TqaGJDMZFsdMTIhjpqfly4mm\naWyLSJQtkcYEEcfY9XEY5WQuOGVRDkImn42218eOHTverzMU7k9lpUbi8/OMmCRGJsUxg4PimLEx\ncczQgP4YTdPYFh4S7yeRxiARx8iUMzwlX07mgpMX5Wg/SE7W0U1015kWRuMR8PvRRFsqJzeLcoTq\nq2rO0YJVOQfjOjiZ+EGiglMNe1Jwfl0myQpmj0X7NUxmBMF0iGsEM8/BHhS8pkK8X14KZp+DPdfJ\nQg3mCFamk7WmAoKZglWZbIBgrmAlxqLRRPMEKzLZgOtgjuAcThcm9dTj4I4OewXnrgYnX4rFwWUS\nBLME+9dUMEozfzGgnGDOXZVJn70jTbTP57TggPnxGvUEM++L9mmYzMi04ORr7TiZC06N85tPA4J5\n5Fgwc7oQgrtyKDj5BBnHTsEBCOZ9dSXps3dIMCvGVsE4B+fwqyvJn34cWwWjF53DyQZ7L5PYt+yg\nBtsnuK6JtZUuBjWYIpfzwaY/mcDK1eyBhsRyst7JgmDeOdh0J8tvXjAbW5toCLZtqJIt2E/fGZLt\n2STUYBvHotmCfWkE29pEM+eDAwHTdx6pI9hYlKO+pTZ1tYhaDXOrUNQ2NDPKWaGR8JRalKOWfs0g\nw0U5alvqU4+fdHg3L8qR5YGODDpZdU3s22+5NTh5HDGOrU20+XsH1anBxlFsEtzFFBxIc5lkr2D0\nok0JzsZQZRaug7tQg+37b5RWrjZZg23vZDEPD8E2CeZcB2fzMsmWr1hBMAfudXDiVgcF+5n3RUOw\nbYIDP/058zrYl6UazBacBARnnlGAeQ7O3kgWBPME2zTZkMFYtK0T/hDMEWzXdCF7NindOdjeyyR8\ndcVpwXWcyyT+SJbN32yw56sr1TYNzbtGsF3ThexzcDrB2RjoMJ+G36Zv8rhHsNPXwYlPs37TnenZ\nJL/vVh/zF5V6mpeCf/pzxshyjif8zf9PFP6iuUW2DJioKDinE/423TZbUjibIdhfUmL25j3lBLOn\nC33094+yLth0zSsumF2QmsesmTNnccspnDmTUZBygtnThdkbqrSpk1VSeEVRaks0S4NbzsyZLMPq\nCWZOF+a4iU5G4kxTcmtJ6taCWbMKuOXkteB090VzyPH/VVnNukwKFBXxz8FMv3kjmMYDgpUbyYJg\nGgjmAcF0CARzgGCbBZ/8g0RGPeKQznaJjLrFMe+/L47pPq0/Juax84R4P4k0Pt0hjpEp5/ge+XKc\nEWxw9tdWS9A5stWecnbuymy/35yx5fDDT9hSDPlkkz3lQLABBHOAYAr1BA812PE2CDnZZk85wUOZ\n7dc4YMvhQzb9np7aYU851gUDdwPBigPBimNdML1KaeaFxNY3taWQjEpyTxrW8kjGquDkVUozLyS2\nvqkthWRQknvSsJZHKlYFx1YptUqskNj6prYUkkFJ7knDWh6pWG+id9nRtmmFxNY3taWQjEpyTxrW\n8kjGPYKN9U2tF5JRSe5Jw1oeybhHcGx9U1sKyagk96RhLY9k3CM4tr6pLYVkVJJ70rCWRzK4DlYc\nCFYcCFacfBbcfcuXL7lkKtfvwmHyWfDG/UuOfhzJ9btwmHwW/Dn5p1y/BefJZ8GELM31G3Ce/Bb8\n4w+qc/0WnCa/Bf/bjNpcvwWnyW/B0z25fgeOk9+C8wAIVhwIVhwIVhwIVhwIVhwIVhwIVhwIVhwI\nVhwIVpz/B4oRnRiJY999AAAAAElFTkSuQmCC\n"
	}
	],
	"prompt_number": 10
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": "## License\n\nThe code is free (CC0). The data and results are currently not available for reuse."
	}
	],
	"metadata": {}
	}
	]
	}