webbedfeet/MungingFWF.ipynb

## MungingFWF.ipynb
{
 "metadata": {
  "name": "",
  "signature": "sha256:39fab77ffafd33d6d07b89acea59fb4789edffc8e8680c3f8e0d32cec41d1030"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Munging fixed-width formats\n",
      "\n",
      "## Introduction\n",
      "Large data sets, such as publicly available US Government data, are often made available as _fixed-width format (FWF)_ text files. At first glance, each record looks like a string of characters. The data dictionary or another file is usually provided that provides information about which columns or substrings represent each variable in the data set. This format requires some careful coding to read, to ensure the correctness of the mapping from data to each variable. \n",
      "\n",
      "Our current project involves the Healthcare Cost and Utilization Project ([HCUP](http://www.hcup-us.ahrq.gov)) and their National Inpatient Sample ([NIS](http://www.hcup-us.ahrq.gov/db/nation/nis/nisdbdocumentation.jsp)). We need to munge data from their Core database to obtain the volume of discharges with particular [DRG codes](http://en.wikipedia.org/wiki/Diagnosis-related_group), and then classify hospitals by volume. The HCUP provides load files for their data in [SAS](http://www.sas.com), [SPSS](http://www-01.ibm.com/software/analytics/spss/) and [Stata](http://www.stata.com). No open-source software is supported, but, the data is just ASCII data; how hard can this be?\n",
      "\n",
      "> For this post, I will consider the 2008 NIS Core data and examples from it. I will not provide the data (as it is subject to a data use agreement with HCUP) or the exact analysis we are doing, but provide surrogate numbers to reflect the process.\n",
      "\n",
      "## Figuring out formats\n",
      "Well, it's probably not hard, but it might be tedious. The NIS Core data for 2008 has 128 fields, so we need variable names, widths (how many characters for each field) and variable types for each of these. I don't know about you, but having to type this out would be a barrier for me to do anything with it. \n",
      "\n",
      "Still, I don't want to fork over thousands of dollars for supported software. But maybe I can take advantage of code that HCUP provides? Looking at the load files, I figured that the Stata code seemed to have the cleanest format. Why not just extract the needed information from the load file? With Python, it shouldn't be too hard. I cooked up the following code (which might give my friends who are professional Python coders the shivers):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import re\n",
      "import pandas as pd\n",
      "import numpy as np\n",
      "\n",
      "def parseStataCode(statafile):\n",
      "    with open(statafile,'r') as f:\n",
      "        statacode = f.readlines()\n",
      "    indStart = min([i for i in range(len(statacode)) if statacode[i].find('infix')>-1])\n",
      "    indEnd = min([i for i in range(len(statacode)) if statacode[i].find('using')>-1])\n",
      "    coldata =statacode[indStart:indEnd]\n",
      "    coldata[0]=coldata[0].replace('infix','')\n",
      "    startcol=[re.findall(' (\\d+)', x)[0] for x in coldata];\n",
      "    maxString = re.findall(r'(\\d+)\\-*(\\d*)', coldata[-1])[0][-1]\n",
      "    startcol.append(maxString)\n",
      "    startcol=np.array([int(x) for x in startcol])\n",
      "    widths = np.diff(startcol)\n",
      "    varnames = [re.findall('[A-Z]\\w+', x)[0] for x in coldata]\n",
      "    out = pd.DataFrame({'varnames':varnames, 'widths':widths })\n",
      "    return out\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, we can get the information we need from the Stata load file by running"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "wt = parseStataCode('StataLoad_NIS_2008_Core.Do')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We could add to the function if we also needed the missing data coding that is provided, but for the current project, we don't need those. \n",
      "\n",
      "## Ingesting data\n",
      "\n",
      "The NIS data is provided as zipped text files. One advantage to Python is the `zipfile` module, which allows us to access the data without explicitly uncompressing the file. Yes, we'll take a computational hit at the expense of extra storage, but for this exercise that's probably fine. Also, for some of the years, we found that there was some problems with the compression metadata that made WinZip fail, and the Python method worked beautifully for that. \n",
      "\n",
      "> For 2010 and after, HCUP encrypted the data with AES-256 encryption. Currently the `zipfile` module cannot handle this, so we used [7-Zip](http://www.7-zip.com) via the `subprocess` module. \n",
      "\n",
      "We will use the [`struct`](https://docs.python.org/2/library/struct.html) module to parse the FWF files, following the advice given on [StackOverflow](http://stackoverflow.com/a/4915359/164965). For this, we need to specify the data format, in terms of widths, in a format string, in order to parse the data.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def defineFormat(wt, keep=['DRG','HOSPID']):\n",
      "    widths = list(wt.widths)\n",
      "    varnames=list(wt.varnames)\n",
      "    fieldwidths = [-x for x in widths]\n",
      "    for s in keep:\n",
      "        fieldwidths[varnames.index(s)] *= -1\n",
      "    fmtstring = ' '.join('{}{}'.format(abs(fw), 'x' if fw < 0 else 's') for fw in fieldwidths)\n",
      "    return fmtstring\n",
      "\n",
      "defineFormat(wt)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "'3x 3x 2x 2x 1x 3x 2x 2x 2x 11x 2x 2x 2x 2x 3s 3x 2x 3x 17x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 5x 5x 5x 5x 2x 4x 4x 4x 4x 2x 3x 2x 5s 2x 14x 5x 6x 2x 2x 5x 5x 2x 2x 3x 2x 4x 2x 2x 2x 10x 2x 10x 3x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 1x 8x 2x 10x 15x 2x 4x 1x'"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This allows us to mask variables we don't need (all the x's), and extract variables we do need as strings.\n",
      "\n",
      "We can then define a parser using this format:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import struct\n",
      "fieldstruct = struct.Struct(defineFormat(wt))\n",
      "parse = fieldstruct.unpack_from"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can now extract data from the compressed data file. We will also filter the data to only keep data with DRG code 500. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import zipfile\n",
      "dat = []\n",
      "drgcodes = ['500']\n",
      "with zipfile.ZipFile('NIS_2008_Core.exe','r') as zf:\n",
      "    with zf.open(\"NIS_2008_Core.ASC\", 'r') as f:\n",
      "        for line in f:\n",
      "            x = list(parse(line))\n",
      "            if x[0] in drgcodes:\n",
      "                dat.append(x)\n",
      "\n",
      "dat = pd.DataFrame(dat, columns = ['DRG','HOSPID'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We now have to compute the volumes for each hospital, categorize them as Low, Medium and High, and then find the frequency of each category. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def categorize(s, bins=[0,10,20,5000], labels=['Low','Medium','High']):\n",
      "    \"\"\"\n",
      "    Split into categories\n",
      "    \"\"\"\n",
      "    out = pd.cut(s, bins=bins, labels=labels)\n",
      "    return out\n",
      "\n",
      "volumes = categorize(dat.HOSPID.value_counts())\n",
      "freq = volumes.value_counts()\n",
      "freq"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "Low       365\n",
        "Medium     10\n",
        "High        1\n",
        "dtype: int64"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "freq.plot(kind='bar',rot=30);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAESCAYAAAAPEjVtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF/xJREFUeJzt3X+UXGV9x/H3khABAw5ozU90UyBCinVBE39QYEBE0ArY\nY1GslVWUg1REWlvIaT1ArRBoQU+16mmFs0ElmmqloA1CgEvlUKG1LCAxJVDjYSkJ8kvC4dSCSf94\n7jCTZXczmzt3797neb/OmTP3uXNn58k+2e8887nPzIAkSZIkSZIkSZIkSZIkSdKUmwHcBVyXt/cB\nbgTuB24AGh3HLgc2AOuBY6ewj5KUrF26PO5sYB2wLW+fRyjmi4Gb8jbAEuC9+fVxwJcm8RiSpBIt\nBNYCR9Gema8H5uTbc/M2hFn5uR33vR540xT0UZKS1s2s+XPAnwJbO/bNATbn25tpF/b5wEjHcSPA\ngoJ9lCTtwI6K+e8CjxLy8r5xjtlGO34Z73ZJUolm7uD2twAnAO8AdgP2Ar5GmI3PBTYB8wgFH+Bh\nYN+O+y/M921nv/322/bggw8W6rgkJehuYKDoDzmSdmZ+Ke1s/DxgRb69BBgGZgGLgAcZe0a/LWbn\nn39+1V1QAY5ffcU+dkyQdOxoZv6iIpxfrwBWA6cBG4GT8/3r8v3rgOeBMyd68Fht3Lix6i6oAMev\nvlIeu8kU81vzC8ATwDHjHHdRfpEkTRHXgJdgcHCw6i6oAMevvlIeu/FWqJQtj38kSd3q6+uDceq2\nM/MSZFlWdRdUgONXXymPncVckiJgzCJJNWHMIkmRs5iXIOXcLgaOX32lPHYWc0mKgJm5JNWEmbkk\nRc5iXoKUc7sYOH71lfLYWcwlKQJm5pJUE2bmkhQ5i3kJUs7tYuD41VfKYzfZL6eopb322octW56s\nuhul2XPPvXn66Seq7oakCiWRmYecKeaMvg/PQUjxMzOXpMhZzEuRVd0BFZBy7lp3KY/djor5bsAd\nwDDhS5ovzvdfAIwAd+WX4zvusxzYAKwHju1hXyVJ4+gmM98DeJZwsvQ24FPAW4EtwOWjjl0CXA0s\nBRYAa4HFwNZRx5mZ95SZuZSCopn5s/n1LGAG0FoWMtYPPBFYBTwHbAQeAJZ131VJ0s7oppjvQohZ\nNgO3APfl+88C7gauABr5vvmE+KVlhDBDT0xWdQdUQMq5a92lPHbdFPOtwACwEDgCaAJfBhbl+x8B\nLpvg/r7+l6SSTeZNQ78Evg+8ge2nnl8Frsu3Hwb27bhtYb7vRQYHB+nv7weg0WgwMDBAs9kE2s+u\nvWoHGeF5qLVNie2pf7wsy0r7/aXWbu2bLv2x3X272WxOq/4UbWdZxtDQEMAL9XI8OzoB+grgeeAp\nYHfgB8CFhKhlU37MOYQTnu+nfQJ0Ge0ToPvz4tm5J0B7yhOgUgqKnACdB9xMyMzvIMzAbwIuBe4h\nZOZHEgo6hOWLq/PrNcCZxF1Fx5FV3QEV0JoZqX5SHrsdxSz3AoeOsf+DE9znovwiSZoifjZLFIxZ\npBT42SySFDmLeSmyqjugAlLOXesu5bGzmEtSBMzMo2BmLqXAzFySImcxL0VWdQdUQMq5a92lPHYW\nc0mKgJl5FMzMpRSYmUtS5Czmpciq7oAKSDl3rbuUx85iLkkRMDOPgpm5lAIzc0mKnMW8FFnVHVAB\nKeeudZfy2FnMJSkCZuZRMDOXUmBmLkmRs5iXIqu6Ayog5dy17lIeux0V890IX+Q8TPiS5ovz/fsA\nNwL3AzcAjY77LAc2AOuBY3vZWUnS2LrJzPcAniV8+fNtwKeAE4DHgEuBc4G9gfOAJcDVwFJgAbAW\nWAxsHfUzzcx7ysxcSkHRzPzZ/HoWMAN4klDMV+b7VwIn5dsnAquA54CNwAPAsp3osyRpErop5rsQ\nYpbNwC3AfcCcvE1+PSffng+MdNx3hDBDT0xWdQdUQMq5a92lPHYzuzhmKzAAvAz4AXDUqNu3MXGG\n4et/SSpZN8W85ZfA94HXE2bjc4FNwDzg0fyYh4F9O+6zMN/3IoODg/T39wPQaDQYGBig2WwC7WfX\nXrWDDGh2bFNie+ofL8uy0n5/qbVb+6ZLf2x33242m9OqP0XbWZYxNDQE8EK9HM+OToC+AngeeArY\nnTAzvxB4O/A4cAnhxGeD7U+ALqN9AnR/Xjw79wRoT3kCVEpBkROg84CbCZn5HcB1wE3ACuBthKWJ\nR+dtCMsXV+fXa4AzibuKjiOrugMqoDUzUv2kPHY7ilnuBQ4dY/8TwDHj3Oei/CJJmiJ+NksUjFmk\nFPjZLJIUOYt5KbKqO6ACUs5d6y7lsbOYS1IEzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZz\nSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYqcxbwUWdUdUAEp\n5651l/LYWcwlKQJm5lEwM5dSYGYuSZGzmJciq7oDKiDl3LXuUh67bor5vsAtwH3AT4BP5PsvAEaA\nu/LL8R33WQ5sANYDx/aor5KkcXSTmc/NL8PAbODHwEnAycAW4PJRxy8BrgaWAguAtcBiYGvHMWbm\nPWVmLqWgaGa+iVDIAZ4Bfkoo0uP90BOBVcBzwEbgAWBZ172VJE3aZDPzfuAQ4Ed5+yzgbuAKoJHv\nm0+IX1pGaBf/RGRVd0AFpJy71l3KYzeZYj4b+DZwNmGG/mVgETAAPAJcNsF9zQAkqUQzuzxuV+A7\nwNeBa/J9j3bc/lXgunz7YcJJ05aF+b7tDA4O0t/fD0Cj0WBgYIBmswm0n1171Q4yoNmxTYntqX+8\nLMtK+/2l1m7tmy79sd19u9lsTqv+FG1nWcbQ0BDAC/VyPN2cAO0DVgKPA+d07J9HmJGT718KvJ/2\nCdBltE+A7s/2s3NPgPaUJ0ClFBQ9AXoY8AHgKLZfhngJcA8hMz+SdqFfB6zOr9cAZxJ3JR1DVnUH\nVEBrZqT6SXnsuolZbmPsor9mgvtclF8kSVPAz2aJgjGLlAI/m0WSImcxL0VWdQdUQMq5a92lPHYW\nc0mKgJl5FMzMpRSYmUtS5Czmpciq7oAKSDl3rbuUx85iLkkRMDOPgpm5lAIzc0mKnMW8FFnVHVAB\nKeeudZfy2FnMJSkCZuZRMDOXUmBmLkmRs5iXIqu6Ayog5dy17lIeO4u5JEXAzDwKZuZSCszMJSly\nFvNSZFV3QAWknLvWXcpj100x3xe4BbgP+AnwiXz/PsCNwP3ADUCj4z7LgQ3AeuDYXnVWkjS2bjLz\nufllGJgN/Bg4CfgQ8BhwKXAusDdwHrAEuBpYCiwA1gKLga0dP9PMvKfMzKUUFM3MNxEKOcAzwE8J\nRfoEYGW+fyWhwAOcCKwCngM2Ag8AyybfbUlStyabmfcDhwB3AHOAzfn+zXkbYD4w0nGfEULxT0hW\ndQdUQMq5a92lPHYzJ3HsbOA7wNnAllG3bWPiHONFtw0ODtLf3w9Ao9FgYGCAZrMJtAekV+0gA5od\n25TYHt7B7b1uh39zWb+/1NrDw8PTqj+2021nWcbQ0BDAC/VyPN2uM98V+B6wBvh8vm89oZpsAuYR\nTpIeSMjNAVbk19cD5xNm8y1m5j1lZi6loGhm3gdcAayjXcgBrgVOzbdPBa7p2P8+YBawCDgAuHOy\nnZYkda+bYn4Y8AHgKOCu/HIcYeb9NsLSxKNpz8TXAavz6zXAmcQ9LR5DVnUHVEDrZa7qJ+Wx6yYz\nv43xi/4x4+y/KL9IkqaAn80SBTNzKQV+NoskRc5iXoqs6g6ogJRz17pLeews5pIUATPzKJiZSykw\nM5ekyFnMS5FV3QEVkHLuWncpj53FXJIiYGYeBTNzKQVm5pIUOYt5KbKqO6ACUs5d6y7lsbOYS1IE\nzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZzSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXet\nu5THrptifiWwGbi3Y98FwAjtL3g+vuO25cAGYD1wbE96KUmaUDeZ+eHAM8BVwGvzfecDW4DLRx27\nBLgaWAosANYCi4Gto44zM+8pM3MpBUUz8x8CT471c8fYdyKwCngO2Ag8ACzrppOSpJ1XJDM/C7gb\nuAJo5PvmE+KXlhHCDD0xWdUdUAEp5651l/LYzdzJ+30Z+Mt8+zPAZcBp4xw75uv/wcFB+vv7AWg0\nGgwMDNBsNoH2gPSqHWRAs2ObEtvDJf/80e3wby7r95dae3h4eFr1x3a67SzLGBoaAnihXo6n23Xm\n/cB1tDPz8W47L9+3Ir++npCv3zHqPmbmPWVmLqWgjHXm8zq23017pcu1wPuAWcAi4ADgzp18DElS\nl7op5quA24HXAA8BHwYuAe4hZOZHAufkx64DVufXa4AziXtKPI6s6g6ogNbLXNVPymPXTWZ+yhj7\nrpzg+IvyiyRpivjZLFEwM5dS4GezSFLkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYqc\nxbwUWdUdUAEp5651l/LYWcwlKQJm5lEwM5dSYGYuSZGzmJciq7oDKiDl3LXuUh47i7kkRcDMPApm\n5lIKzMwlKXIW81JkVXdABaScu9ZdymNnMZekCJiZR8HMXEqBmbkkRc5iXoqs6g6ogJRz17pLeey6\nKeZXApuBezv27QPcCNwP3AA0Om5bDmwA1gPH9qabkqSJdJOZHw48A1wFvDbfdynwWH59LrA3cB6w\nBLgaWAosANYCi4Gto36mmXlPmZlLKSiamf8QeHLUvhOAlfn2SuCkfPtEYBXwHLAReABYNqneSpIm\nbWcz8zmE6IX8ek6+PR8Y6ThuhDBDT0xWdQdUQMq5a92lPHYze/AztjFxhjHmbYODg/T39wPQaDQY\nGBig2WwC7QHpVTvIgGbHNiW2h0v++aPb4d9c1u8vtfbw8PC06o/tdNtZljE0NATwQr0cT7frzPuB\n62hn5usJlWQTMA+4BTiQkJsDrMivrwfOB+4Y9fPMzHvKzFxKQRnrzK8FTs23TwWu6dj/PmAWsAg4\nALhzJx9DktSlbor5KuB24DXAQ8CHCDPvtxGWJh5Neya+DlidX68BziTuKfE4sqo7oAJaL3NVPymP\nXTeZ+Snj7D9mnP0X5RdJ0hTxs1miYGYupcDPZpGkyFnMS5FV3QEVkHLuWncpj53FXJIiYGYeBTNz\nKQVm5pIUOYt5KbKqO6ACUs5d6y7lsbOYS1IEzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZz\nSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYpc0WK+EbgHuIv2\nFzfvA9xI+H7QG4BGwceooazqDqiAlHPXukt57IoW821AEzgEWJbvO49QzBcDN+VtSVKJimbmPwPe\nADzesW89cCSwGZhLmKYeOOp+ZuY9ZWYupaDMzHwbsBb4D+Cj+b45hEJOfj2n4GNIknagaDE/jBCx\nHA/8EXD4qNu3EfeUeBxZ1R1QASnnrnWX8tjNLHj/R/LrXwDfJeTmrXhlEzAPeHSsOw4ODtLf3w9A\no9FgYGCAZrMJtAekV+0gI8T7rW1KbA+X/PNHt8O/uazfX2rt4eHhadUf2+m2syxjaGgI4IV6OZ4i\nmfkewAxgC/BSwsqVC4FjCBn6JYSTnw1efBLUzLynzMylFEyUmRcp5osIs3EIM/xvABcTliauBl5F\nWLp4MvDUqPtazHvKYi6loKxiXkTkxTyjMwIpn8W8lzojK9VL7GPnO0AlKXLOzKPgzFxKgTNzSYqc\nxbwUWdUdUAGtpWGqn5THzmIuSREwM4+CmbmUAjNzSYqcxbwUWdUdUAEp5651l/LYWcwlKQJm5lEw\nM5dSYGYuSZGzmJciq7oDKiDl3LXuUh47i7kkRcDMPApm5lIKzMwlKXIW81JkVXdABaScu9ZdymNn\nMZekCJiZR8HMXErBRJn5zKntijR5e+21D1u2PFl1N0qx55578/TTT1TdDUWgrJjlOGA9sAE4t6TH\nmMayqjsQlVDIt03h5ZYpe6xYn6SqYmbeWzOALxIK+hLgFOCgEh5nGhuuugMqxPGrq+HhdMeujGK+\nDHgA2Ag8B3wTOLGEx5nGnqq6AyrE8aurp55Kd+zKKOYLgIc62iP5PklSSco4AeqyCjZW3QEVsrHq\nDkSjipPXF1544ZQ91nQ6gV1GMX8Y2LejvS9hdt7p7r6+vteV8NgTmOpVmCun9NHyJUsRi3f84h+7\neG3Z8uRUj9/dU/lgM4EHgX5gFuFsUmInQCUpDscD/0U4Ebq84r5IkiRJUjAjvzYcjY9jKiXAP/S4\nzdjxIarYDBL5O/Q/Y3lGf7rXxcCbgB8S1ve7hLO+ZtB+T/4s4CxgEfCTKjul7exLGKdn8/ZuwPPV\ndad8FvPea/2hQ/gPdCpwF/Ar4DTgauDX1XRNBbXeZLc1v94N+EdgNvBeYFfg3/GJumqvJiy8eDlh\nKd9XgHcBPwc2VdivUlnMe6/zD7lJ+GN/kvAfaVdCIf/51HdLPdAa27cA3yK8s/kR4BzgHuCj+fXm\nSnqn1kTqV8BewFLgI4QxeQY4A7gdmB7v8ukxi3lxu9COTfqAlwBXEf6g1wNz8+33Eor6g4RibtRS\nD53jtBdwOvBO4HvAPOAdhPH+b0KM9hrgDuD/pryn6dqF8LfXesU0m/B3djhhPD4E3AocCexDiMMc\nH22n88TKfGD/fPs04F+AQ4F/JuR3lwKPAX+TH+O3PNVLP/AqwnsnLsv3vYrwEv5THcfcC7x1ivum\n4NXAJwnFug84ivB3eHR+++8A382vo+PMfHLGOiveB1wCnA+8nvAS/FLCjHweoaD/mvCxwIcCbwBW\nAf87Bf3Vzmm9XO8jvKP5dOBk4NuEsdwfuBF4PG+/B/hPwmxwBPgRjm+ZdgF+D9idEHP1ASuAvyT8\n/g/Oj/kmocC/jvAlAz8jFPb5wL8R+QlRTaz16Y+7duw7CPgOsCfhj/zvgc/nt50ErAOuye+zMD9O\n088cwpNwy0HAAfn2EcDlhKhsD+Am4N35bS8HvgT88dR0M3kHEgr1EYQiPovwhHs18Ir8mOMIr5Dm\nEAr5VwhRC8DehBPXSth8tv+g69Ys/Qzgunx7JmFGfjvtyOWzhP9oe0xBHzV5s4FX5ttvpB1//RXw\nrx3HnQ78LfAbhKJ+LeHJGdpFROV6B+HV0cF5+88Ir5j2Av6H9njMJMzKW5OqPwHOzve3xteYM3H/\nQHsGNiu/fjnh5NeSjvZK4DentmvaSV8BvkB4sp1LWEa6J/AywnsCfj8/7rWE8x9n5u1rCJGaytUZ\nBR8EfBr487x9PiFemU14VXVVx7F/DdxHu/BL29mDsKxp97zdilvOI6wv3o/wH+1GQjHQ9NRHe+wW\nAtcDh+Xt6wmzcoA/IGSrLd8nzAz7aT+Zqzz7E76QtdMJhCfV3yJMmL4IfBBoECZVZxOenL8OfJX2\nyWlwNq5RzgCG8u2XdOxfQZi5r8ZvVqqLgwlZ+F3AmnzfAYQTZa28/FZCTHYN8DXCUjdNnbsJSz4B\nPkDIwm8lrCiaAQwS4pR5hFfHZxGejHclFPW3T213VSd9hDPm/Xn7CMIM4EicrU1no1duHUBYgTJI\nyMofAt6f33YZIXOFkId/FPhY+V3UGBYR3rV5LWFt/wDhJOgXCbP0lwF/B/wF7cnVYsLsfQ1hNYs0\nrjcTZnNfIHzH2Ccq7Y0m0npDSUtrJcPBhMik9QT8QUI81iAUhccJbw5S9b5AWDHWsifhCfYKwuqU\nowlLfgFeCnyYMIuXunIzYUb+kh0dqErsnV+3CvlbCauMLiPM6A4BvkF440/LCOED0SAU8kXld1Nd\n2AP4Be1zVRDG73LCW/alQjyZMv20Pk7hnwjv/GvNwt9JWD66GDiFsMT0AMLnq5xLeIfubGAt4Qna\n9wJMP6cDV3a0+9i+uEuKUAZsAD6et99DWJlyDnAnYfURwG8TTlzfDPyU9htLNP30EZ6E9xtjfxKf\nVy6l4F2EwgzhJfmngQsIJ8wWEGZ1DwOfI7wrEMKJsdn59hsJH8Ck6a315i6LtxShOYRPy7uV9pLQ\nzxJm3B8nZOAvJSxpOyq//XWE9eMfmdKeqlcs5lKkVhBOXp4B/CFh2ei3CB+Beh1hKds7Ce8PuAH4\nMSF2kSRNI7sDWwifSnkt4VMsP0NYZvhJwscrtLj6QZKmsY8RPiv+lYQ3jNxLKPIH5O2jx7+rJGm6\n6COc5FycbzcJyxRn4acaSlKtvJnwxRCSpJq7nfYyRUlSTflViJIkSZIkSZIkSZIkSZIkSZIkKXn/\nD/IdTSApgZyaAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x107569b90>"
       ]
      }
     ],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}
	{
	"metadata": {
	"name": "",
	"signature": "sha256:39fab77ffafd33d6d07b89acea59fb4789edffc8e8680c3f8e0d32cec41d1030"
	},
	"nbformat": 3,
	"nbformat_minor": 0,
	"worksheets": [
	{
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"# Munging fixed-width formats\n",
	"\n",
	"## Introduction\n",
	"Large data sets, such as publicly available US Government data, are often made available as _fixed-width format (FWF)_ text files. At first glance, each record looks like a string of characters. The data dictionary or another file is usually provided that provides information about which columns or substrings represent each variable in the data set. This format requires some careful coding to read, to ensure the correctness of the mapping from data to each variable. \n",
	"\n",
	"Our current project involves the Healthcare Cost and Utilization Project ([HCUP](http://www.hcup-us.ahrq.gov)) and their National Inpatient Sample ([NIS](http://www.hcup-us.ahrq.gov/db/nation/nis/nisdbdocumentation.jsp)). We need to munge data from their Core database to obtain the volume of discharges with particular [DRG codes](http://en.wikipedia.org/wiki/Diagnosis-related_group), and then classify hospitals by volume. The HCUP provides load files for their data in [SAS](http://www.sas.com), [SPSS](http://www-01.ibm.com/software/analytics/spss/) and [Stata](http://www.stata.com). No open-source software is supported, but, the data is just ASCII data; how hard can this be?\n",
	"\n",
	"> For this post, I will consider the 2008 NIS Core data and examples from it. I will not provide the data (as it is subject to a data use agreement with HCUP) or the exact analysis we are doing, but provide surrogate numbers to reflect the process.\n",
	"\n",
	"## Figuring out formats\n",
	"Well, it's probably not hard, but it might be tedious. The NIS Core data for 2008 has 128 fields, so we need variable names, widths (how many characters for each field) and variable types for each of these. I don't know about you, but having to type this out would be a barrier for me to do anything with it. \n",
	"\n",
	"Still, I don't want to fork over thousands of dollars for supported software. But maybe I can take advantage of code that HCUP provides? Looking at the load files, I figured that the Stata code seemed to have the cleanest format. Why not just extract the needed information from the load file? With Python, it shouldn't be too hard. I cooked up the following code (which might give my friends who are professional Python coders the shivers):"
	]
	},
	{
	"cell_type": "code",
	"collapsed": true,
	"input": [
	"import re\n",
	"import pandas as pd\n",
	"import numpy as np\n",
	"\n",
	"def parseStataCode(statafile):\n",
	" with open(statafile,'r') as f:\n",
	" statacode = f.readlines()\n",
	" indStart = min([i for i in range(len(statacode)) if statacode[i].find('infix')>-1])\n",
	" indEnd = min([i for i in range(len(statacode)) if statacode[i].find('using')>-1])\n",
	" coldata =statacode[indStart:indEnd]\n",
	" coldata[0]=coldata[0].replace('infix','')\n",
	" startcol=[re.findall(' (\\d+)', x)[0] for x in coldata];\n",
	" maxString = re.findall(r'(\\d+)\\-(\\d)', coldata[-1])[0][-1]\n",
	" startcol.append(maxString)\n",
	" startcol=np.array([int(x) for x in startcol])\n",
	" widths = np.diff(startcol)\n",
	" varnames = [re.findall('[A-Z]\\w+', x)[0] for x in coldata]\n",
	" out = pd.DataFrame({'varnames':varnames, 'widths':widths })\n",
	" return out\n"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 1
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"So, we can get the information we need from the Stata load file by running"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"wt = parseStataCode('StataLoad_NIS_2008_Core.Do')"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 2
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We could add to the function if we also needed the missing data coding that is provided, but for the current project, we don't need those. \n",
	"\n",
	"## Ingesting data\n",
	"\n",
	"The NIS data is provided as zipped text files. One advantage to Python is the `zipfile` module, which allows us to access the data without explicitly uncompressing the file. Yes, we'll take a computational hit at the expense of extra storage, but for this exercise that's probably fine. Also, for some of the years, we found that there was some problems with the compression metadata that made WinZip fail, and the Python method worked beautifully for that. \n",
	"\n",
	"> For 2010 and after, HCUP encrypted the data with AES-256 encryption. Currently the `zipfile` module cannot handle this, so we used [7-Zip](http://www.7-zip.com) via the `subprocess` module. \n",
	"\n",
	"We will use the [`struct`](https://docs.python.org/2/library/struct.html) module to parse the FWF files, following the advice given on [StackOverflow](http://stackoverflow.com/a/4915359/164965). For this, we need to specify the data format, in terms of widths, in a format string, in order to parse the data. "
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"def defineFormat(wt, keep=['DRG','HOSPID']):\n",
	" widths = list(wt.widths)\n",
	" varnames=list(wt.varnames)\n",
	" fieldwidths = [-x for x in widths]\n",
	" for s in keep:\n",
	" fieldwidths[varnames.index(s)] *= -1\n",
	" fmtstring = ' '.join('{}{}'.format(abs(fw), 'x' if fw < 0 else 's') for fw in fieldwidths)\n",
	" return fmtstring\n",
	"\n",
	"defineFormat(wt)"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"metadata": {},
	"output_type": "pyout",
	"prompt_number": 3,
	"text": [
	"'3x 3x 2x 2x 1x 3x 2x 2x 2x 11x 2x 2x 2x 2x 3s 3x 2x 3x 17x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 5x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 5x 5x 5x 5x 2x 4x 4x 4x 4x 2x 3x 2x 5s 2x 14x 5x 6x 2x 2x 5x 5x 2x 2x 3x 2x 4x 2x 2x 2x 10x 2x 10x 3x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 4x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 1x 8x 2x 10x 15x 2x 4x 1x'"
	]
	}
	],
	"prompt_number": 3
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"This allows us to mask variables we don't need (all the x's), and extract variables we do need as strings.\n",
	"\n",
	"We can then define a parser using this format:"
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import struct\n",
	"fieldstruct = struct.Struct(defineFormat(wt))\n",
	"parse = fieldstruct.unpack_from"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 4
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We can now extract data from the compressed data file. We will also filter the data to only keep data with DRG code 500. "
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import zipfile\n",
	"dat = []\n",
	"drgcodes = ['500']\n",
	"with zipfile.ZipFile('NIS_2008_Core.exe','r') as zf:\n",
	" with zf.open(\"NIS_2008_Core.ASC\", 'r') as f:\n",
	" for line in f:\n",
	" x = list(parse(line))\n",
	" if x[0] in drgcodes:\n",
	" dat.append(x)\n",
	"\n",
	"dat = pd.DataFrame(dat, columns = ['DRG','HOSPID'])"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 5
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"We now have to compute the volumes for each hospital, categorize them as Low, Medium and High, and then find the frequency of each category. "
	]
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"def categorize(s, bins=[0,10,20,5000], labels=['Low','Medium','High']):\n",
	" \"\"\"\n",
	" Split into categories\n",
	" \"\"\"\n",
	" out = pd.cut(s, bins=bins, labels=labels)\n",
	" return out\n",
	"\n",
	"volumes = categorize(dat.HOSPID.value_counts())\n",
	"freq = volumes.value_counts()\n",
	"freq"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"metadata": {},
	"output_type": "pyout",
	"prompt_number": 6,
	"text": [
	"Low 365\n",
	"Medium 10\n",
	"High 1\n",
	"dtype: int64"
	]
	}
	],
	"prompt_number": 6
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"%matplotlib inline\n",
	"freq.plot(kind='bar',rot=30);"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"metadata": {},
	"output_type": "display_data",
	"png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAESCAYAAAAPEjVtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF/xJREFUeJzt3X+UXGV9x/H3khABAw5ozU90UyBCinVBE39QYEBE0ArY\nY1GslVWUg1REWlvIaT1ArRBoQU+16mmFs0ElmmqloA1CgEvlUKG1LCAxJVDjYSkJ8kvC4dSCSf94\n7jCTZXczmzt3797neb/OmTP3uXNn58k+2e8887nPzIAkSZIkSZIkSZIkSZIkSdKUmwHcBVyXt/cB\nbgTuB24AGh3HLgc2AOuBY6ewj5KUrF26PO5sYB2wLW+fRyjmi4Gb8jbAEuC9+fVxwJcm8RiSpBIt\nBNYCR9Gema8H5uTbc/M2hFn5uR33vR540xT0UZKS1s2s+XPAnwJbO/bNATbn25tpF/b5wEjHcSPA\ngoJ9lCTtwI6K+e8CjxLy8r5xjtlGO34Z73ZJUolm7uD2twAnAO8AdgP2Ar5GmI3PBTYB8wgFH+Bh\nYN+O+y/M921nv/322/bggw8W6rgkJehuYKDoDzmSdmZ+Ke1s/DxgRb69BBgGZgGLgAcZe0a/LWbn\nn39+1V1QAY5ffcU+dkyQdOxoZv6iIpxfrwBWA6cBG4GT8/3r8v3rgOeBMyd68Fht3Lix6i6oAMev\nvlIeu8kU81vzC8ATwDHjHHdRfpEkTRHXgJdgcHCw6i6oAMevvlIeu/FWqJQtj38kSd3q6+uDceq2\nM/MSZFlWdRdUgONXXymPncVckiJgzCJJNWHMIkmRs5iXIOXcLgaOX32lPHYWc0mKgJm5JNWEmbkk\nRc5iXoKUc7sYOH71lfLYWcwlKQJm5pJUE2bmkhQ5i3kJUs7tYuD41VfKYzfZL6eopb322octW56s\nuhul2XPPvXn66Seq7oakCiWRmYecKeaMvg/PQUjxMzOXpMhZzEuRVd0BFZBy7lp3KY/djor5bsAd\nwDDhS5ovzvdfAIwAd+WX4zvusxzYAKwHju1hXyVJ4+gmM98DeJZwsvQ24FPAW4EtwOWjjl0CXA0s\nBRYAa4HFwNZRx5mZ95SZuZSCopn5s/n1LGAG0FoWMtYPPBFYBTwHbAQeAJZ131VJ0s7oppjvQohZ\nNgO3APfl+88C7gauABr5vvmE+KVlhDBDT0xWdQdUQMq5a92lPHbdFPOtwACwEDgCaAJfBhbl+x8B\nLpvg/r7+l6SSTeZNQ78Evg+8ge2nnl8Frsu3Hwb27bhtYb7vRQYHB+nv7weg0WgwMDBAs9kE2s+u\nvWoHGeF5qLVNie2pf7wsy0r7/aXWbu2bLv2x3X272WxOq/4UbWdZxtDQEMAL9XI8OzoB+grgeeAp\nYHfgB8CFhKhlU37MOYQTnu+nfQJ0Ge0ToPvz4tm5J0B7yhOgUgqKnACdB9xMyMzvIMzAbwIuBe4h\nZOZHEgo6hOWLq/PrNcCZxF1Fx5FV3QEV0JoZqX5SHrsdxSz3AoeOsf+DE9znovwiSZoifjZLFIxZ\npBT42SySFDmLeSmyqjugAlLOXesu5bGzmEtSBMzMo2BmLqXAzFySImcxL0VWdQdUQMq5a92lPHYW\nc0mKgJl5FMzMpRSYmUtS5Czmpciq7oAKSDl3rbuUx85iLkkRMDOPgpm5lAIzc0mKnMW8FFnVHVAB\nKeeudZfy2FnMJSkCZuZRMDOXUmBmLkmRs5iXIqu6Ayog5dy17lIeux0V890IX+Q8TPiS5ovz/fsA\nNwL3AzcAjY77LAc2AOuBY3vZWUnS2LrJzPcAniV8+fNtwKeAE4DHgEuBc4G9gfOAJcDVwFJgAbAW\nWAxsHfUzzcx7ysxcSkHRzPzZ/HoWMAN4klDMV+b7VwIn5dsnAquA54CNwAPAsp3osyRpErop5rsQ\nYpbNwC3AfcCcvE1+PSffng+MdNx3hDBDT0xWdQdUQMq5a92lPHYzuzhmKzAAvAz4AXDUqNu3MXGG\n4et/SSpZN8W85ZfA94HXE2bjc4FNwDzg0fyYh4F9O+6zMN/3IoODg/T39wPQaDQYGBig2WwC7WfX\nXrWDDGh2bFNie+ofL8uy0n5/qbVb+6ZLf2x33242m9OqP0XbWZYxNDQE8EK9HM+OToC+AngeeArY\nnTAzvxB4O/A4cAnhxGeD7U+ALqN9AnR/Xjw79wRoT3kCVEpBkROg84CbCZn5HcB1wE3ACuBthKWJ\nR+dtCMsXV+fXa4AzibuKjiOrugMqoDUzUv2kPHY7ilnuBQ4dY/8TwDHj3Oei/CJJmiJ+NksUjFmk\nFPjZLJIUOYt5KbKqO6ACUs5d6y7lsbOYS1IEzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZz\nSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYqcxbwUWdUdUAEp\n5651l/LYWcwlKQJm5lEwM5dSYGYuSZGzmJciq7oDKiDl3LXuUh67bor5vsAtwH3AT4BP5PsvAEaA\nu/LL8R33WQ5sANYDx/aor5KkcXSTmc/NL8PAbODHwEnAycAW4PJRxy8BrgaWAguAtcBiYGvHMWbm\nPWVmLqWgaGa+iVDIAZ4Bfkoo0uP90BOBVcBzwEbgAWBZ172VJE3aZDPzfuAQ4Ed5+yzgbuAKoJHv\nm0+IX1pGaBf/RGRVd0AFpJy71l3KYzeZYj4b+DZwNmGG/mVgETAAPAJcNsF9zQAkqUQzuzxuV+A7\nwNeBa/J9j3bc/lXgunz7YcJJ05aF+b7tDA4O0t/fD0Cj0WBgYIBmswm0n1171Q4yoNmxTYntqX+8\nLMtK+/2l1m7tmy79sd19u9lsTqv+FG1nWcbQ0BDAC/VyPN2cAO0DVgKPA+d07J9HmJGT718KvJ/2\nCdBltE+A7s/2s3NPgPaUJ0ClFBQ9AXoY8AHgKLZfhngJcA8hMz+SdqFfB6zOr9cAZxJ3JR1DVnUH\nVEBrZqT6SXnsuolZbmPsor9mgvtclF8kSVPAz2aJgjGLlAI/m0WSImcxL0VWdQdUQMq5a92lPHYW\nc0mKgJl5FMzMpRSYmUtS5Czmpciq7oAKSDl3rbuUx85iLkkRMDOPgpm5lAIzc0mKnMW8FFnVHVAB\nKeeudZfy2FnMJSkCZuZRMDOXUmBmLkmRs5iXIqu6Ayog5dy17lIeO4u5JEXAzDwKZuZSCszMJSly\nFvNSZFV3QAWknLvWXcpj100x3xe4BbgP+AnwiXz/PsCNwP3ADUCj4z7LgQ3AeuDYXnVWkjS2bjLz\nufllGJgN/Bg4CfgQ8BhwKXAusDdwHrAEuBpYCiwA1gKLga0dP9PMvKfMzKUUFM3MNxEKOcAzwE8J\nRfoEYGW+fyWhwAOcCKwCngM2Ag8AyybfbUlStyabmfcDhwB3AHOAzfn+zXkbYD4w0nGfEULxT0hW\ndQdUQMq5a92lPHYzJ3HsbOA7wNnAllG3bWPiHONFtw0ODtLf3w9Ao9FgYGCAZrMJtAekV+0gA5od\n25TYHt7B7b1uh39zWb+/1NrDw8PTqj+2021nWcbQ0BDAC/VyPN2uM98V+B6wBvh8vm89oZpsAuYR\nTpIeSMjNAVbk19cD5xNm8y1m5j1lZi6loGhm3gdcAayjXcgBrgVOzbdPBa7p2P8+YBawCDgAuHOy\nnZYkda+bYn4Y8AHgKOCu/HIcYeb9NsLSxKNpz8TXAavz6zXAmcQ9LR5DVnUHVEDrZa7qJ+Wx6yYz\nv43xi/4x4+y/KL9IkqaAn80SBTNzKQV+NoskRc5iXoqs6g6ogJRz17pLeews5pIUATPzKJiZSykw\nM5ekyFnMS5FV3QEVkHLuWncpj53FXJIiYGYeBTNzKQVm5pIUOYt5KbKqO6ACUs5d6y7lsbOYS1IE\nzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZzSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXet\nu5THrptifiWwGbi3Y98FwAjtL3g+vuO25cAGYD1wbE96KUmaUDeZ+eHAM8BVwGvzfecDW4DLRx27\nBLgaWAosANYCi4Gto44zM+8pM3MpBUUz8x8CT471c8fYdyKwCngO2Ag8ACzrppOSpJ1XJDM/C7gb\nuAJo5PvmE+KXlhHCDD0xWdUdUAEp5651l/LYzdzJ+30Z+Mt8+zPAZcBp4xw75uv/wcFB+vv7AWg0\nGgwMDNBsNoH2gPSqHWRAs2ObEtvDJf/80e3wby7r95dae3h4eFr1x3a67SzLGBoaAnihXo6n23Xm\n/cB1tDPz8W47L9+3Ir++npCv3zHqPmbmPWVmLqWgjHXm8zq23017pcu1wPuAWcAi4ADgzp18DElS\nl7op5quA24HXAA8BHwYuAe4hZOZHAufkx64DVufXa4AziXtKPI6s6g6ogNbLXNVPymPXTWZ+yhj7\nrpzg+IvyiyRpivjZLFEwM5dS4GezSFLkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYqc\nxbwUWdUdUAEp5651l/LYWcwlKQJm5lEwM5dSYGYuSZGzmJciq7oDKiDl3LXuUh47i7kkRcDMPApm\n5lIKzMwlKXIW81JkVXdABaScu9ZdymNnMZekCJiZR8HMXEqBmbkkRc5iXoqs6g6ogJRz17pLeey6\nKeZXApuBezv27QPcCNwP3AA0Om5bDmwA1gPH9qabkqSJdJOZHw48A1wFvDbfdynwWH59LrA3cB6w\nBLgaWAosANYCi4Gto36mmXlPmZlLKSiamf8QeHLUvhOAlfn2SuCkfPtEYBXwHLAReABYNqneSpIm\nbWcz8zmE6IX8ek6+PR8Y6ThuhDBDT0xWdQdUQMq5a92lPHYze/AztjFxhjHmbYODg/T39wPQaDQY\nGBig2WwC7QHpVTvIgGbHNiW2h0v++aPb4d9c1u8vtfbw8PC06o/tdNtZljE0NATwQr0cT7frzPuB\n62hn5usJlWQTMA+4BTiQkJsDrMivrwfOB+4Y9fPMzHvKzFxKQRnrzK8FTs23TwWu6dj/PmAWsAg4\nALhzJx9DktSlbor5KuB24DXAQ8CHCDPvtxGWJh5Neya+DlidX68BziTuKfE4sqo7oAJaL3NVPymP\nXTeZ+Snj7D9mnP0X5RdJ0hTxs1miYGYupcDPZpGkyFnMS5FV3QEVkHLuWncpj53FXJIiYGYeBTNz\nKQVm5pIUOYt5KbKqO6ACUs5d6y7lsbOYS1IEzMyjYGYupcDMXJIiZzEvRVZ1B1RAyrlr3aU8dhZz\nSYqAmXkUzMylFJiZS1LkLOalyKrugApIOXetu5THzmIuSREwM4+CmbmUAjNzSYpc0WK+EbgHuIv2\nFzfvA9xI+H7QG4BGwceooazqDqiAlHPXukt57IoW821AEzgEWJbvO49QzBcDN+VtSVKJimbmPwPe\nADzesW89cCSwGZhLmKYeOOp+ZuY9ZWYupaDMzHwbsBb4D+Cj+b45hEJOfj2n4GNIknagaDE/jBCx\nHA/8EXD4qNu3EfeUeBxZ1R1QASnnrnWX8tjNLHj/R/LrXwDfJeTmrXhlEzAPeHSsOw4ODtLf3w9A\no9FgYGCAZrMJtAekV+0gI8T7rW1KbA+X/PNHt8O/uazfX2rt4eHhadUf2+m2syxjaGgI4IV6OZ4i\nmfkewAxgC/BSwsqVC4FjCBn6JYSTnw1efBLUzLynzMylFEyUmRcp5osIs3EIM/xvABcTliauBl5F\nWLp4MvDUqPtazHvKYi6loKxiXkTkxTyjMwIpn8W8lzojK9VL7GPnO0AlKXLOzKPgzFxKgTNzSYqc\nxbwUWdUdUAGtpWGqn5THzmIuSREwM4+CmbmUAjNzSYqcxbwUWdUdUAEp5651l/LYWcwlKQJm5lEw\nM5dSYGYuSZGzmJciq7oDKiDl3LXuUh47i7kkRcDMPApm5lIKzMwlKXIW81JkVXdABaScu9ZdymNn\nMZekCJiZR8HMXErBRJn5zKntijR5e+21D1u2PFl1N0qx55578/TTT1TdDUWgrJjlOGA9sAE4t6TH\nmMayqjsQlVDIt03h5ZYpe6xYn6SqYmbeWzOALxIK+hLgFOCgEh5nGhuuugMqxPGrq+HhdMeujGK+\nDHgA2Ag8B3wTOLGEx5nGnqq6AyrE8aurp55Kd+zKKOYLgIc62iP5PklSSco4AeqyCjZW3QEVsrHq\nDkSjipPXF1544ZQ91nQ6gV1GMX8Y2LejvS9hdt7p7r6+vteV8NgTmOpVmCun9NHyJUsRi3f84h+7\neG3Z8uRUj9/dU/lgM4EHgX5gFuFsUmInQCUpDscD/0U4Ebq84r5IkiRJUjAjvzYcjY9jKiXAP/S4\nzdjxIarYDBL5O/Q/Y3lGf7rXxcCbgB8S1ve7hLO+ZtB+T/4s4CxgEfCTKjul7exLGKdn8/ZuwPPV\ndad8FvPea/2hQ/gPdCpwF/Ar4DTgauDX1XRNBbXeZLc1v94N+EdgNvBeYFfg3/GJumqvJiy8eDlh\nKd9XgHcBPwc2VdivUlnMe6/zD7lJ+GN/kvAfaVdCIf/51HdLPdAa27cA3yK8s/kR4BzgHuCj+fXm\nSnqn1kTqV8BewFLgI4QxeQY4A7gdmB7v8ukxi3lxu9COTfqAlwBXEf6g1wNz8+33Eor6g4RibtRS\nD53jtBdwOvBO4HvAPOAdhPH+b0KM9hrgDuD/pryn6dqF8LfXesU0m/B3djhhPD4E3AocCexDiMMc\nH22n88TKfGD/fPs04F+AQ4F/JuR3lwKPAX+TH+O3PNVLP/AqwnsnLsv3vYrwEv5THcfcC7x1ivum\n4NXAJwnFug84ivB3eHR+++8A382vo+PMfHLGOiveB1wCnA+8nvAS/FLCjHweoaD/mvCxwIcCbwBW\nAf87Bf3Vzmm9XO8jvKP5dOBk4NuEsdwfuBF4PG+/B/hPwmxwBPgRjm+ZdgF+D9idEHP1ASuAvyT8\n/g/Oj/kmocC/jvAlAz8jFPb5wL8R+QlRTaz16Y+7duw7CPgOsCfhj/zvgc/nt50ErAOuye+zMD9O\n088cwpNwy0HAAfn2EcDlhKhsD+Am4N35bS8HvgT88dR0M3kHEgr1EYQiPovwhHs18Ir8mOMIr5Dm\nEAr5VwhRC8DehBPXSth8tv+g69Ys/Qzgunx7JmFGfjvtyOWzhP9oe0xBHzV5s4FX5ttvpB1//RXw\nrx3HnQ78LfAbhKJ+LeHJGdpFROV6B+HV0cF5+88Ir5j2Av6H9njMJMzKW5OqPwHOzve3xteYM3H/\nQHsGNiu/fjnh5NeSjvZK4DentmvaSV8BvkB4sp1LWEa6J/AywnsCfj8/7rWE8x9n5u1rCJGaytUZ\nBR8EfBr487x9PiFemU14VXVVx7F/DdxHu/BL29mDsKxp97zdilvOI6wv3o/wH+1GQjHQ9NRHe+wW\nAtcDh+Xt6wmzcoA/IGSrLd8nzAz7aT+Zqzz7E76QtdMJhCfV3yJMmL4IfBBoECZVZxOenL8OfJX2\nyWlwNq5RzgCG8u2XdOxfQZi5r8ZvVqqLgwlZ+F3AmnzfAYQTZa28/FZCTHYN8DXCUjdNnbsJSz4B\nPkDIwm8lrCiaAQwS4pR5hFfHZxGejHclFPW3T213VSd9hDPm/Xn7CMIM4EicrU1no1duHUBYgTJI\nyMofAt6f33YZIXOFkId/FPhY+V3UGBYR3rV5LWFt/wDhJOgXCbP0lwF/B/wF7cnVYsLsfQ1hNYs0\nrjcTZnNfIHzH2Ccq7Y0m0npDSUtrJcPBhMik9QT8QUI81iAUhccJbw5S9b5AWDHWsifhCfYKwuqU\nowlLfgFeCnyYMIuXunIzYUb+kh0dqErsnV+3CvlbCauMLiPM6A4BvkF440/LCOED0SAU8kXld1Nd\n2AP4Be1zVRDG73LCW/alQjyZMv20Pk7hnwjv/GvNwt9JWD66GDiFsMT0AMLnq5xLeIfubGAt4Qna\n9wJMP6cDV3a0+9i+uEuKUAZsAD6et99DWJlyDnAnYfURwG8TTlzfDPyU9htLNP30EZ6E9xtjfxKf\nVy6l4F2EwgzhJfmngQsIJ8wWEGZ1DwOfI7wrEMKJsdn59hsJH8Ck6a315i6LtxShOYRPy7uV9pLQ\nzxJm3B8nZOAvJSxpOyq//XWE9eMfmdKeqlcs5lKkVhBOXp4B/CFh2ei3CB+Beh1hKds7Ce8PuAH4\nMSF2kSRNI7sDWwifSnkt4VMsP0NYZvhJwscrtLj6QZKmsY8RPiv+lYQ3jNxLKPIH5O2jx7+rJGm6\n6COc5FycbzcJyxRn4acaSlKtvJnwxRCSpJq7nfYyRUlSTflViJIkSZIkSZIkSZIkSZIkSZIkKXn/\nD/IdTSApgZyaAAAAAElFTkSuQmCC\n",
	"text": [
	"<matplotlib.figure.Figure at 0x107569b90>"
	]
	}
	],
	"prompt_number": 7
	}
	],
	"metadata": {}
	}
	]
	}