{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#02:-xarray,-netcdf-and-zarr\" data-toc-modified-id=\"02:-xarray,-netcdf-and-zarr-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>02: xarray, netcdf and zarr</a></span><ul class=\"toc-item\"><li><span><a href=\"#The-current-defacto-standard-in-atmos/ocean-science\" data-toc-modified-id=\"The-current-defacto-standard-in-atmos/ocean-science-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>The current defacto standard in atmos/ocean science</a></span></li><li><span><a href=\"#Some-challenges-with-netcdf\" data-toc-modified-id=\"Some-challenges-with-netcdf-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Some challenges with netcdf</a></span></li><li><span><a href=\"#create-an-xarray\" data-toc-modified-id=\"create-an-xarray-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>create an xarray</a></span></li><li><span><a href=\"#Download-toy-model-data\" data-toc-modified-id=\"Download-toy-model-data-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>Download toy model data</a></span></li><li><span><a href=\"#Sort-in-numeric-order\" data-toc-modified-id=\"Sort-in-numeric-order-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>Sort in numeric order</a></span></li><li><span><a href=\"#Make-an-xarray\" data-toc-modified-id=\"Make-an-xarray-1.6\"><span class=\"toc-item-num\">1.6&nbsp;&nbsp;</span>Make an xarray</a></span></li><li><span><a href=\"#Dump-to-a-zarr-file\" data-toc-modified-id=\"Dump-to-a-zarr-file-1.7\"><span class=\"toc-item-num\">1.7&nbsp;&nbsp;</span>Dump to a zarr file</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 02: xarray, netcdf and zarr\n",
    "\n",
    "Motivation:  how you store your data can an enormous effect on performance.\n",
    "\n",
    "Four issues:\n",
    "\n",
    "1) Compression vs. cpu time to uncompress/compress\n",
    "\n",
    "2) Multithreaded read/writes\n",
    "\n",
    "3) Performance for cloud storage (Amazon, Google Compute, Azure)\n",
    "\n",
    "4) Throttling data reads to fit in available memory and avoid swapping (\"chunking\")\n",
    "\n",
    "## The current defacto standard in atmos/ocean science\n",
    "\n",
    "\n",
    "* [netcdf/hdf5](ttps://www.unidata.ucar.edu/software/netcdf/docs/index.html)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some challenges with netcdf\n",
    "\n",
    "* [The Pangeo project](http://pangeo-data.org/)\n",
    "\n",
    "* [Cloud challenges](https://medium.com/pangeo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## create an xarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import numpy as np\n",
    "import pdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!conda install -y xarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "******************************\n",
      "context imported. Front of path:\n",
      "/Users/phil/repos/parallel_python_course/notebooks\n",
      "\n",
      "******************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import context\n",
    "from westgrid.data_read import download"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download toy model data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mar12014_1_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_1_15600.nc\n",
      "writing to: mar12014_1_15600.nc\n",
      "\n",
      "mar12014_1_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_2_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_2_15600.nc\n",
      "writing to: mar12014_2_15600.nc\n",
      "\n",
      "mar12014_2_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_3_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_3_15600.nc\n",
      "writing to: mar12014_3_15600.nc\n",
      "\n",
      "mar12014_3_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_4_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_4_15600.nc\n",
      "writing to: mar12014_4_15600.nc\n",
      "\n",
      "mar12014_4_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_5_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_5_15600.nc\n",
      "writing to: mar12014_5_15600.nc\n",
      "\n",
      "mar12014_5_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_6_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_6_15600.nc\n",
      "writing to: mar12014_6_15600.nc\n",
      "\n",
      "mar12014_6_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_7_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_7_15600.nc\n",
      "writing to: mar12014_7_15600.nc\n",
      "\n",
      "mar12014_7_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_8_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_8_15600.nc\n",
      "writing to: mar12014_8_15600.nc\n",
      "\n",
      "mar12014_8_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_9_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_9_15600.nc\n",
      "writing to: mar12014_9_15600.nc\n",
      "\n",
      "mar12014_9_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n",
      "mar12014_10_15600.nc\n",
      "trying http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les/mar12014_10_15600.nc\n",
      "writing to: mar12014_10_15600.nc\n",
      "\n",
      "mar12014_10_15600.nc already exists\n",
      "and is 1455454 bytes\n",
      "will not overwrite\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#get 10 files, each is the same timestep for a member of a\n",
    "#10 member ensemble\n",
    "\n",
    "import numpy as np\n",
    "root='http://clouds.eos.ubc.ca/~phil/docs/atsc500/small_les'\n",
    "for i in range(10):\n",
    "    the_name=f'mar12014_{(i+1):d}_15600.nc'\n",
    "    print(the_name)\n",
    "    url='{}/{}'.format(root,the_name)\n",
    "    download(the_name,root=root)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_2_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_3_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_10_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_8_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_9_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_5_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_4_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_1_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_6_15600.nc'),\n",
       " PosixPath('/Users/phil/repos/parallel_python_course/notebooks/mar12014_7_15600.nc')]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#get 10 files, each is the same timestep for a member of a\n",
    "#10 member ensemble\n",
    "import numpy as np\n",
    "root=Path().resolve()\n",
    "the_files=root.glob('mar12*nc')\n",
    "the_files=list(the_files)\n",
    "the_files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sort in numeric order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import xarray\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def sort_name(pathname):\n",
    "    \"\"\"\n",
    "      sort the filenames so '10' sorts\n",
    "      last by converting to integers\n",
    "    \"\"\"\n",
    "    str_name=str(pathname.name)\n",
    "    front, number, back = str_name.split('_')\n",
    "    return int(number)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make an xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use xarray to stitch together the 10 ensemble members along a new \"virtual dimenstion\".\n",
    "The variable \"ds\"  is an xray dataset, which controls the reads/writes from the\n",
    "10 netcdf files\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (ens: 10, time: 1, x: 50, y: 40, z: 30)\n",
      "Coordinates:\n",
      "  * x        (x) float32 0.0 25.0 50.0 75.0 100.0 125.0 150.0 175.0 200.0 ...\n",
      "  * y        (y) float32 0.0 25.0 50.0 75.0 100.0 125.0 150.0 175.0 200.0 ...\n",
      "  * z        (z) float32 25.0 50.0 75.0 100.0 125.0 150.0 175.0 200.0 205.0 ...\n",
      "  * time     (time) float32 169.18056\n",
      "Dimensions without coordinates: ens\n",
      "Data variables:\n",
      "    p        (ens, z) float32 dask.array<shape=(10, 30), chunksize=(1, 30)>\n",
      "    U        (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n",
      "    V        (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n",
      "    W        (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n",
      "    PP       (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n",
      "    TABS     (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n",
      "    QV       (ens, z, y, x) float32 dask.array<shape=(10, 30, 40, 50), chunksize=(1, 30, 40, 50)>\n"
     ]
    }
   ],
   "source": [
    "the_files.sort(key=sort_name)\n",
    "\n",
    "#\n",
    "#  put the 10 ensembles together along a new \"ens\" dimension\n",
    "#  using an xray multi-file dataset\n",
    "#\n",
    "ds = xarray.open_mfdataset(the_files, engine='netcdf4', concat_dim='ens')\n",
    "# dump the structure\n",
    "print(ds)\n",
    "#\n",
    "#  3-d ensemble average for temp\n",
    "#\n",
    "x = ds['x']\n",
    "y = ds['y']\n",
    "z = ds['z']\n",
    "temp = ds['TABS']\n",
    "mean_temp = temp[:,  :, :, :].mean(dim='ens')\n",
    "#\n",
    "# same for velocity\n",
    "#\n",
    "wvel = ds['W']\n",
    "mean_w = wvel[:,  :, :, :].mean(dim='ens')\n",
    "#\n",
    "# now look at the perturbation fields for one ensemble member\n",
    "#\n",
    "wvelprime = wvel[0,  :, :, :] - mean_w\n",
    "Tprime = temp[0,  :, :, :] - mean_temp\n",
    "flux_prime = wvelprime * Tprime\n",
    "flux_profile = flux_prime.mean(dim='x').mean(dim='y')\n",
    "keep_dict = dict(flux_prof=flux_profile, flux_prime=flux_prime.values,\n",
    "                 wvelprime=wvelprime.values, Tprime=Tprime.values, x=x, y=y, z=z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dump to a zarr file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.backends.zarr.ZarrStore at 0x312a84f28>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.to_zarr('zarr_dir','w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYFeXZx/HvvR3Ype5SpHeBqIBICShYUOwtxoI1KtZY8xo1RVNMYhrG2AvBhiUaOxoRQUADCNKlF+lFENil7LLs8/4xszAczu6eLWdny+9zXefa6XOf2TlzzzzPMzPmnENERCRSQtgBiIhI1aQEISIiUSlBiIhIVEoQIiISlRKEiIhEpQQhIiJRKUFImZjZAjMbUs5lPGhmL/vd7czMmVlShQRYA/jbo1MR4yaa2XVFjPvIzK6Kb3TlZ2bDzeyTMsw32sx+H4+YYli3mdm/zOx7M5tuZkPMbG0YsVQGJYgyMLNVZrbHzHICn8fisJ47zWyjme0ws1FmllrR64gxjsN+kM65Hs65iWHEI8Vzzp3unHsh7DiCop0AOOdecc6dGmZcQTGepAwChgKtnHN9Kym00ChBlN3Zzrn0wOfWily4mZ0G3AucDLQDOgC/qch1xBhHYmWvsyryzxz1e5G2wCrn3K6wA6kM2uErmJldbWZTzOyv/mXoSjM7PWL8CjPL9scNL2JRVwHPO+cWOOe+B34HXB1jDB+b2a0Rw+aY2QV+95FmNs7MtpnZYjP7cWC60Wb2pJmNNbNdwLXAcOAe/0rpfX+6VWZ2it+daGb3m9ly/3vNNLPW/rh/mNkaM9vpDz8+hvgvMrOZEcPuNrN3ipj+GjNb6K97hZndEBi30MzOCvQnmdl3Ztbb7+9vZl+a2XZ/Gw0JTDvRzB4ysy+A3UCH4tblz3OPmW0ws/Vmdl2wmMjMUv39YrWZbTKzp8ysTmDe/wvM+5OStlNgvhZmNtfMfhaI+zq/u6T9sYGZPe+vd52Z/b7wpMCf9wszG+lvnxVm9kN/+Boz22yBoiwzO9PMZvn/6zVm9mAgzEn+3+3+fjSgMLbA/D0C++UmM7u/mK/dyMw+9P8P08ysY2A5xe3fpYoxYjtfCzwHDPDHH3bCZhHFgha4+jazn5vZVPOvUMzsJvOKatOK+Z7hcs7pU8oPsAo4pYhxVwP7gOuBROAmYD1gQD1gJ9DVn7YF0KOI5cwBLg70ZwIOaOL3fwDcW8S8VwJfBPq7A9uBVD+GNcA1QBLQG/iuMA5gNLADGIh3ApHmD/t9UdsA+D9gHtDV/57HBOK8HGjir+tuYCOQ5o97EHjZ727nf78kP85tQLfA+mYBFxbxfc8EOvrrHox3MO/tj/s18ErEtIv87pbAVuAM/7sO9fuz/PETgdVADz+u5BLWNcz/fj2AusBL/nfq5I9/BHgPaAxkAO8DfwzMuwn4gf8/GhOcN8p3nghc52+3JcCIyHEl7Y/++HeAp/11NgWmAzcE5s3H21cSgd/72+Nx/390KpANpPvTDwGO8rfl0f73OS/y/xvxW5nid2cAG/D2kTS/v18R33003v7R1/+/vAK85o8raf8uVYxF/L6nBPqHAGsD/Yf8zwj8dvx1TsLb7zsD3wO9wj6eFXusCzuA6vjBOzjm4B10Cz/XB3agZYFp6/o7TXN/590OXAjUKWEdy4Fhgf5kfzntYogvA9gFtPX7HwJG+d0XA5Mjpn8aeMDvHg28GDH+wE4esQ0KE8Ri4NwYt933wDF+94NESRB+/5PAQ353D3++1BjX8Q5wu9/dCe8gVtfvfwX4td/9c+CliHn/C1zld08EfluKdY3CP+AH1u38v+b/TzoGxg8AVgbm/VNgXBdKThB/9/8Pl0YZF0wQRe2PzYDc4L4IXApMCMy7NDDuKH/eZoFhW4GeRcT4CDAy2v83sPwpgfXOivH/Oxp4LtB/BgeTfrH7d2ljjDL9gZj9/iHEmCAC69gGLATui+X7hvlREVPZneecaxj4PBsYt7Gwwzm32+9Md1655cXAjcAG/xL5yCKWnwPUD/QXdmeXFJhzLhv4ELjEH3QJ3oERvDLUfn6RwXYz245XhNQ8sIg1Ja0jQmu8hHYYv2hooXkV7duBBnhXQyV5AbjMzAy4AnjDOZdbxDpO9y/dt/nrOKNwHc65ZXg/xrPNrC5wDt7ZOXjb4qKIbTEI78qu0CHborh1AUdETB/szsI7OM8MrOtjf3i0eb8tdut4hgPrgDdLmC7q/oj3/ZPx9sXCmJ7Gu5IotCnQvcdfRuSwdAAz62dmE8xsi5ntwNvPY/lfQzH7UBE2Brp3F8ZACft3OWMsN+fcKmACXqJ4vLLWW1ZKEJXMOfdf59xQvIPQIuDZIiZdgFdUU+gYYJNzbmuMq3oVuNQvR62Dt1OCdxD6PCK5pTvnbgqGGRl2Cetag1fscgjz6ht+DvwYaOSca4hXfGUlBe+cmwrkAccDl+EV1xzGvJZdbwF/xTuzbQiMjVjHq3hnqOcC3/hJozDulyK2RT3n3J+CoZRiXRuAVoF5Wwe6v8M7mPYIrKuBcy49MG9w+jZFbZuAB/3ljrGyNSZYg3cFkRmIqb5zrkcZlgVe4n0PaO2cawA8xcFtU6Z9qAxK2r/LE2MsduOdCBQKnnhhZmfgXTmOB/5SAeuLKyWISmRmzczsHDOrh/fDzAH2FzH5i8C1ZtbdzBoBv8S7XI3VWLyzqd8CrzvnCvzhHwBdzOwKM0v2P8eZWbdilrUJrxVVUZ4Dfmdmnc1ztJk1wSvqyge2AElm9msOvSoqyYvAY0C+c25KEdOk4JWHbwHy/QrYyKaTr/nDbuLg1QPAy3hXFqeZV9GeZl679lZEV9K63gCuMbNu/tXKrwtH+Nv/WWCkmTUFMLOW5rVWK5z3av//XRd4oKiNErAPuAiv6PIlK2UrK+fcBuAT4G9mVt/MEsyso5kNLs1yAjKAbc65vWbWFy+xF9oCFFD0fvQB0NzM7jCvMj/DzPqVIYaS9u/yxBiL2XhXvolmNgyvngoAM8sEnserO7oKb987oxzrijsliLJ73w69D+LtGOZJwKuEW49XDjkYuDnahM65j4E/4535f+t/Dhw0zLsZqshWHn5xzH+AUwgcFP3ip1Pxip3W412qP4x34CvK80B3/5I9Wkuiv+Md4D7Bq4R/Hu+q5b/AR3iVqN8Ceyld8dVLeJW2Ua8eAt/nNn/93+P94N+LmGYD8D/gh8DrgeFr8K4q7sc7OKzBq3CP+rsoaV3OuY+AR/H+Z8v8dYJ3MgDe1dQyYKqZ7QQ+xavYL5z3EeAzf5rPivrOETHlARfgFQuNKm2SwGvQkAJ843+nNzm0iK00bgZ+a2bZeMnxjUCcu/Hqwr7w96P+Ed8jG6+RwNl4++RS4MTSBhDD/l3mGGN0u/8dCou2gr+XZ4B3nXNj/ZKAa4Hn/JOpKqmwJYNIlWNeE9DNeK2EloYdT2n5Z63z8SrX88OOR6S0dAUhVdlNwFfVKTmY2flmluIXCz4MvK/kINWVnnsjVZKZrcKrPDwv5FBK6wa8uqL9wOcUUYQoUh2oiElERKJSEZOIiERVrYuYMjMzXbt27cIOQ0SkWpk5c+Z3zrmskqar1gmiXbt2zJgxI+wwRESqFTOL5U59FTGJiEh0ShAiIhKVEoSIiESlBCEiIlEpQYiISFRKECIiEpUShIiIRKUEITXS6q27+fPHi1i3fU/YoYhUW0oQUiNNXbmVpz5fzvEPf8b1L85g8tIt6LljIqVTre+kFinKj/u05ocdmzBm2mpe+2oN477ZRIfMelzevy0/6tOK+mnJYYcoUuVV66e59unTx+lRG1KS3Pz9jJ23gRf/9y2zVm+nbkoi5/VqyZUD2nJk89K8AVWkZjCzmc65PiVOpwQhtcm8tTt4aeoq3p29ntz8AgZ1yuSmIR35YccmmFnJCxCpAUJPEGaWBkzCexdsEvCmc+4BM2uP9xL5xsDXwBXOuTwzS8V7Sf2xwFbgYufcquLWoQQhZbV9dx6vfbWG56esZEt2Lse0asBNQzpyavfmJCQoUUjNFmuCiGcldS5wknPuGKAnMMx/CfjDwEjnXGe8l6Rf609/LfC9c64TMNKfTiQuGtZN4cbBHZl8z4n84fyj2L5nHze+/DWnjPycN2asIS+/IOwQRUIXtwThPDl+b7L/ccBJwJv+8Bc4+ErJc/1+/PEnm675Jc7SkhO5rF8bPrt7CP+8tBdpSYnc8+ZcBv9lAs9MWs723XlhhygSmrg2czWzRDObDWwGxgHLge2Bl7ivBVr63S2BNQD++B1AkyjLHGFmM8xsxpYtW+IZvtQiiQnG2cccwYe3DeKFn/SlTeO6/GHsIvr9YTz3vDmH+et2hB2iSKWLazNX59x+oKeZNQTeBrpFm8z/G+1q4bAKEufcM8Az4NVBVFCoIgCYGYO7ZDG4SxbfrN/JS1O/5Z1Z63hjxlp6tWnIlQPacsZRLUhNSgw7VJG4q5Qb5Zxz24GJQH+goZkVJqZWwHq/ey3QGsAf3wDYVhnxiUTT/Yj6/PGCo5h6/8n8+qzu7Ni9jztfn8OAP37Gwx8vYu33u8MOUSSu4pYgzCzLv3LAzOoApwALgQnAj/zJrgLe9bvf8/vxx3/mqnMbXKkxGtRJ5ieD2vPpXYN56dq+HNu2EU9/vpzBf5nIHa/NYvHG7LBDFImLeDZzPRqv0jkRLxG94Zz7rZl14GAz11nA5c65XL9Z7EtAL7wrh0uccyuKW4eauUpY1n6/m9FfrGLM9NXsztvP0O7NuOXETvRs3TDs0ERKFPp9EJVBCULC9v2uPEZ/uYrRX65ix559DOzUhFuGdGKAbryTKkwJQqQS5eTmM2batzw72bvxrmfrhvzmnB4coysKqYKqwo1yIrVGemoSI07wbrx76PwfsGnnXi588kue+nw5BQXV9yRMajclCJEKlJacyPB+bfn49hMY2r0Zf/poEVeOms7mnXvDDk2k1JQgROKgQd1knhjemz9ecBQzvt3G6f+YzIRFm8MOS6RUlCBE4sTMuLRvG96/dRBZGalcM/orXpr6bdhhicRMCUIkzjo3y+CdWwZySrem/Oqd+Tw/ZWXYIYnERAlCpBKkJSfyxPBjOf0HzfndB9/w5MTlYYckUiIlCJFKkpKUwD8v7cU5xxzBwx8vYuS4JXpPtlRpeie1SCVKSkxg5MU9SU1K4B/jl7Jjzz5+fVZ3vaRIqiQlCJFKlphgPHzh0dSvk8zzU1ayfXcef7noGJITdUEvVYsShEgIEhKMX57Zjcb1UvjLfxezc28+j1/Wmzopeoy4VB06ZREJiZlxy4mdeOj8HzBh8WYuePJLVmzJKXlGkUqiBCESsuH92jLq6uPYsGMPZ/9zCu/PWV/yTCKVQAlCpAo4sWtTxt52PF2bZ/DTV2fxy3fmsXff/rDDklpOCUKkijiiYR1ev2EAI07owMtTV3PBE18ya/X3YYcltZgShEgVkpyYwP1ndOO5K/uwJSeX85/4kjtfn83GHXrYn1Q+tWISqYJO6d6M/h2b8MSEZTw3ZSUfz9/IjYM7MuKEDmrpJJVGVxAiVVR6ahL3DDuS8XcNZkjXLEZ+uoST/zaRd2ev0zsmpFIoQYhUca0b1+XJy4/ltRH9aVA3hdtfm80Zj07mo3kblCgkrpQgRKqJ/h2a8MFPB/HIxT3Jyy/gple+5oxHJ/Px/I1KFBIXeie1SDWUv7+A9+eu59Hxy1j53S66t6jP7ad05tTuzTDTc52keHontUgNlpSYwPm9WjHuzhP4+4+PYXdePje8NJNzHvuCpZuyww5PagglCJFqLCkxgQt6t+LTuwbz14uOYf32PZz92BRenb5ajxKXclOCEKkBkhIT+NGxrfjo9uPp07Yx9/1nHreOmcWOPfvCDk2qMSUIkRqkaf00XvxJX34+7Ej+u2AjZz46WXdjS5kpQYjUMAkJxk1DOvLGjQMAuOzZaUxfuS3kqKQ6iluCMLPWZjbBzBaa2QIzu90f/qCZrTOz2f7njMA895nZMjNbbGanxSs2kdqgd5tGvH3zQFo0TOOaf03na11JSCnF8woiH7jbOdcN6A/cYmbd/XEjnXM9/c9YAH/cJUAPYBjwhJnpmQIi5ZCVkcqY6/qTmZHKVaOmM3/djrBDkmokbgnCObfBOfe1350NLARaFjPLucBrzrlc59xKYBnQN17xidQWzRukMeb6/tRPS+by56excMPOsEOSaqJS6iDMrB3QC5jmD7rVzOaa2Sgza+QPawmsCcy2ligJxcxGmNkMM5uxZcuWOEYtUnO0bFiHV6/vT53kRIY/N43FG3WvhJQs7gnCzNKBt4A7nHM7gSeBjkBPYAPwt8JJo8x+WENu59wzzrk+zrk+WVlZcYpapOZp06QuY67vT1KCMfy5qSzbrCQhxYtrgjCzZLzk8Ipz7j8AzrlNzrn9zrkC4FkOFiOtBVoHZm8F6N2LIhWofWY9Xh3RHzAufXYayzbrHdhStHi2YjLgeWChc+7vgeEtApOdD8z3u98DLjGzVDNrD3QGpscrPpHaqmNWOq9e3w/nHOc//gWfLNgYdkhSRcXzCmIgcAVwUkST1j+b2TwzmwucCNwJ4JxbALwBfAN8DNzinNNLeUXioHOzDN65ZSDts+ox4qWZ/PnjRezXE2Elgp7mKlKL7d23n9+8v4BXp69hUKdMHr20F43rpYQdlsSZnuYqIiVKS07kjxcczcMXHsX0Vds457EpbM7W+6/FowQhIlx8XBteH9GfrTl5XP/iTPbuU+muKEGIiK9Xm0Y8cklP5q7dzt3/nqO31IkShIgcdFqP5tw77Eg+nLuBRz5dEnY4ErKksAMQkaplxAkdWL4lh0c/W0aLhnW4tG+bsEOSkChBiMghzIyHzj+KTTtz+cXb88hMT2Vo92ZhhyUhUBGTiBwmOTGBJ4b35qiWDbh1zNfM/Fbvk6iNlCBEJKp6qUmMuvo4WjRI4yejZzBnzfawQ5JKpgQhIkVqkp7KS9f2IyMtieHPTeOrVbqSqE2UIESkWK0b1+XfNw6gaUYqVz4/nSlLvws7JKkkShAiUqIWDerw+g0DaNukLj954SsmLNocdkhSCZQgRCQmWRmpvDaiP12bZXDzK18zb61eX1rTKUGISMwa1k1h1NXH0bheCte+8BXrt+8JOySJIyUIESmVrIxU/nXNcezJ289PRn9FTm5+2CFJnChBiEipdWmWwROX92bp5hwe+vCbsMOROFGCEJEyOb5zFhcd24p3Zq0ne+++sMOROFCCEJEyu6RvG/bs28/7czaEHYrEgRKEiJTZMa0acGTzDMZM/1aPB6+BlCBEpMzMjGsHtWf+up388aOFYYcjFUxPcxWRcvnRsa1YsH4nz05eScuGdbh6YPuwQ5IKogQhIuViZvzqrO6s276H33zwDUc0rMOpPZqHHZZUABUxiUi5JSYYj17Si6NbNeT212bz7dZdYYckFUAJQkQqRJ2URJ6+/FiSEox735qHc6q0ru6UIESkwjRvkMb9Z3bjfyu28ur0NWGHI+WkBCEiFeqS41ozoEMT/jB2IQs37Aw7HCmHuCUIM2ttZhPMbKGZLTCz2/3hjc1snJkt9f828oebmT1qZsvMbK6Z9Y5XbCISP2bGwxceTZ2URM57/AvGTFut4qZqKp5XEPnA3c65bkB/4BYz6w7cC4x3znUGxvv9AKcDnf3PCODJOMYmInHUpkldxt52PH3bN+b+t+dx22uz9TiOaihuCcI5t8E597XfnQ0sBFoC5wIv+JO9AJznd58LvOg8U4GGZtYiXvGJSHxlZaTywjV9+b/TujJ23gbO+ucU5q/TOySqk0qpgzCzdkAvYBrQzDm3AbwkAjT1J2sJBGu11vrDRKSaSkgwbjmxE6+N6E9efgFXjprOvv0FYYclMYp7gjCzdOAt4A7nXHE1VhZl2GEFl2Y2wsxmmNmMLVu2VFSYIhJHx7VrzANn92Dbrjy+/vb7sMORGMU1QZhZMl5yeMU59x9/8KbCoiP/b+HLbdcCrQOztwLWRy7TOfeMc66Pc65PVlZW/IIXkQo1sFMTkhKMCYt1YlddxLMVkwHPAwudc38PjHoPuMrvvgp4NzD8Sr81U39gR2FRlIhUfxlpyfRp14iJizeXPLFUCfG8ghgIXAGcZGaz/c8ZwJ+AoWa2FBjq9wOMBVYAy4BngZvjGJuIhODErk1ZtDGbjTv2hh2KxCBuD+tzzk0her0CwMlRpnfALfGKR0TCN6RrU/740SImLt7MJX3bhB2OlEB3UotIpenSLJ0WDdKYqHqIakEJQkQqjZkxpGtTpiz7Ts1dqwElCBGpVEO6ZpGTm8+MVWruWtUpQYhIpRrYKZPkRGPiErVmquqUIESkUqWnJnFcu8Z8rnqIKk8JQkQq3ZCuWSzamM367XvCDkWKoQQhIpVuSFfvEWyfL9FVRFVWYoIws1Zm9jMze9fMvjKzSWb2hJmdaWZKMCJSap2bptOyYR3dVV3FFXuAN7N/AaOAPOBh4FK8O5w/BYYBU8zshHgHKSI1i5kxuGsWXyzbSl6+mrtWVSXdSf0359z8KMPnA/8xsxRAt0OKSKmd2LUpY6atZsa32/hhx8yww5Eoir2CKCI5BMfnOeeWVWxIIlIb/LBjE1ISE9SaqQqLqQ7BzM4ys1lmts3MdppZtpnpbeQiUmb1UpM4rn0jxi3cpHdWV1GxVjI/gvdo7ibOufrOuQznXP04xiUitcAFvVqxYssu/rtgU9ihSBSxJog1wHynNC8iFejcnkfQIasefx+3mP0FOrxUNbEmiHuAsWZ2n5ndVfiJZ2AiUvMlJSZw5yldWLIphw/mHvYCSQlZrAniIWA3kAZkBD4iIuVy5lEtOLJ5Bo98upR8PeG1Son1hUGNnXOnxjUSEamVEhKMu4Z2YcRLM/nPrHX8uE/rkmeSShHrFcSnZqYEISJxMbR7M45u1YB/fLpUN85VIbEmiFuAj81sj5q5ikhFMzPuPrUr67bv4fWvVocdjvhiShB+s9YE51wdNXMVkXg4oXMmfds15p+fLWPvvv1hhyOU/CymdiWMNzNrVZEBiUjt5F1FdGFzdi4vT/027HCEkq8g/mJmb5nZlWbWw8yamlkbMzvJzH4HfAF0q4Q4RaQW6NehCcd3zuSJicvZlZsfdji1XknPYroI+BXQFXgcmAy8C1wHLAZOcs6Ni3eQIlJ73DW0C9t25TH6y1Vhh1LrldjM1Tn3DfCLSohFRIRebRpxSremPP35ci7v35YGdZLDDqnW0gt/RKTKuXNoF3buzef5ySvCDqVWU4IQkSqnxxENOPOoFjw/ZSXbduWFHU6tFbcEYWajzGyzmc0PDHvQzNaZ2Wz/c0Zg3H1mtszMFpvZafGKS0SqhzuHdmbPvv08/fnysEOptWJ9H8T44MHcH/ZMCbONxnstaaSRzrme/mesv6zuwCVAD3+eJ8wsMZbYRKRm6tQ0g/N6tuSF/61i8869YYdTK8V6BdEe+LmZPRAY1qe4GZxzk4BtMS7/XOA151yuc24lsAzoG+O8IlJD3X5KZ/btd4z8dIleKhSCWBPEduBkoJmZvW9mDcqxzlvNbK5fBNXIH9YS750Thdb6ww5jZiPMbIaZzdiyRa8qFKnJ2japx5UD2vLq9DXc+uoscnRvRKWKNUGYcy7fOXcz8BYwBWhahvU9CXQEegIbgL8VLj/KtFFPF5xzzzjn+jjn+mRlZZUhBBGpTn51ZnfuGdaVj+Zt4JzHprB4Y3bYIdUasSaIpwo7nHOjgauBT0q7MufcJufcfudcAfAsB4uR1gLBZ/y2AvT2EBEhIcG4eUgnxlzfn+y9+Zz7+BTenLk27LBqhVgf1vd0RP9M59xPSrsyM2sR6D0fKGzh9B5wiZmlmll7oDMwvbTLF5Gaq3+HJnx42yB6tm7Iz/49h3vfmquH+sVZrC8MKjUzexUYAmSa2VrgAWCImfXEKz5aBdwA4JxbYGZvAN8A+cAtzjn950XkEE0z0nj52n6M/HQJj09Yzty1O3hieG/aZdYLO7Qayapzy4A+ffq4GTNmhB2GiITgs0WbuPP1ORQUOP5y0dEM+0GLkmcSAMxspnOu2JaooDupRaSaOunIZnx42yA6NE3nxpe/5vcffMM+vdO6QilBiEi11apRXf59wwCu/mE7npuykkuemcqGHXvCDqvGUIIQkWotJSmBB8/pwT8v7cWiDTs589EpTFqie6QqghKEiNQIZx9zBO/9dBBZ6alc9a/pjBy3hP0F1beOtSpQghCRGqNjVjrv3DKQC3q14h/jl3LVqOlszckNO6xqSwlCRGqUOimJ/PWio3n4wqOYvmobZz46hRmrYn0snAQpQYhIjWNmXHxcG96++YekJidw8TNTeXbSCj3wr5SUIESkxupxRAPe/+kgTunWlIfGLuTGl2eyY8++sMOqNpQgRKRGq5+WzFOXH8svz+zG+IWbOfufU5i/bkfYYVULShAiUuOZGdcd34HXb+hPXn4BFzz5Ja9OX60ipxIoQYhIrXFs28Z8eNsg+rVvzH3/mcfdb8xhd57eMVEUJQgRqVWapKcy+pq+3HFKZ96evY4RL84kX4/oiEoJQkRqncQE445TuvDwhUczZdl3/P7DhWGHVCXF7XHfIiJV3Y/7tGbJxmyem7KSLs0yuKxfm7BDqlJ0BSEitdp9Z3RjcJcsfv3ufKau2Bp2OFWKEoSI1GqJCcY/L+tF2yZ1uenlmazeujvskKoMJQgRqfXqpyXz3FXHUeDguhe/InuvbqYDJQgREQDaZ9bjieG9Wb5lF3e+PltPgkUJQkTkgIGdMnng7O58unAzf/1kcdjhhE6tmEREAq7o35ZFG7N5cuJyujRL5/xercIOKTS6ghARCTAzfnNOD/p3aMzP35rHrNXfhx1SaJQgREQiJCcm8MTwY2lWP5URL82ste+5VoIQEYmicb0Unr/qOHbn5jPixZnsydsfdkiVTglCRKQIXZpl8OilvZi/fgf/9+acWvf0VyUIEZFinNytGT8fdiQfzN3AY58tCzucSqVWTCIiJbjhhA4s2ZjN38YtoXOzdIb9oEXYIVWKuF1BmNkoM9tsZvMDwxqb2TgzW+r/beQPNzN71MyWmdmulQ7jAAARIUlEQVRcM+sdr7hERErLzPjDBUfRs3VD7nx9DgvW14430sWziGk0MCxi2L3AeOdcZ2C83w9wOtDZ/4wAnoxjXCIipZaWnMgzVxxLgzrJjHhxJt/l5IYdUtzFLUE45yYB2yIGnwu84He/AJwXGP6i80wFGppZ7biGE5Fqo2n9NJ69sg9bd+Vy5+uzww4n7iq7krqZc24DgP+3qT+8JbAmMN1af9hhzGyEmc0wsxlbtmyJa7AiIpGOatWAm4d0YvLS79icvTfscOKqqrRisijDorYnc84945zr45zrk5WVFeewREQOd9KR3rntF8u+CzmS+KrsBLGpsOjI/7vZH74WaB2YrhWwvpJjExGJSfcW9WlcL4XJS5QgKtJ7wFV+91XAu4HhV/qtmfoDOwqLokREqpqEBGNQp0wmLf2uRt88F89mrq8C/wO6mtlaM7sW+BMw1MyWAkP9foCxwApgGfAscHO84hIRqQjHd87ku5xcFm3MDjuUuInbjXLOuUuLGHVylGkdcEu8YhERqWjHd/bqQCcv3UK3FvVDjiY+qkoltYhItdK8QRpdmqUzeWnNrYdQghARKaPjO2cxbeU29u6rmU96VYIQESmjE7pkkZdfwPSVkfcE1wxKECIiZdS3XWNSkhKYtKRm3rSrBCEiUkZ1UhLp265xja2HUIIQESmH4ztnsnhTNpt21rzHbihBiIiUw8HmrjXvKkIJQkSkHI5snkFmeiqTl9a8egglCBGRckhIMI7vnMmUpd9RUFCzHruhBCEiUk7Hd85k6648vtmwM+xQKpQShIhIOQ3qlAnAJ99sCjmSiqUEISJSTk3rpzGwUxMeHb+Un/17Dtt354UdUoVQghARqQDPX3UcNw/pyNuz1nHK3yfx4dwN1f5R4EoQIiIVIC05kXuGHcl7tw6keYNUbhnzNSNemlmt749QghARqUA9jmjAOzcP5L7Tj2TSki2c8rfPGTNtdbVs4aQEISJSwZISE7hhcEf+e8cJ/KBlA+5/ex6XPTeVVd/tCju0UlGCEBGJk3aZ9RhzfT/+dMFRLFi/k9MemcRTny8nf39B2KHFRAlCRCSOzIxL+rbh07sGM6RrFn/6aBHnPfEFC9bvCDu0EilBiIhUgmb103j6ij48Obw3G3fkcs5jX/Dwx4uq9MuGlCBERCrR6Ue1YPxdg7mwd0uenLic0/8xmWkrtoYdVlRKECIilaxB3WT+/KNjeOW6fuwvcFz8zFTuf3seO/fuCzu0QyhBiIiEZGCnTD6+43iuG9Se16av5tS/T2JcFXpchxKEiEiI6qYk8cuzuvP2zQNpWDeZ61+cwS1jvmZLdm7YoSlBiIhUBce0bsh7tw7i7qFdGLdgE0NHfs5bM9eG+rgOJQgRkSoiJSmBn57cmbG3D6JTVjp3/3sOV46azpptu0OJJ5QEYWarzGyemc02sxn+sMZmNs7Mlvp/G4URm4hI2Do1zeCNGwbw23N78PW333PaI5MYNWUl+yv5cR1hXkGc6Jzr6Zzr4/ffC4x3znUGxvv9IiK1UkKCceWAdnxy12D6tW/Mbz/4hguf/JIlm7IrL4ZKW1PJzgVe8LtfAM4LMRYRkSqhZcM6jLr6OB65uCcrtuRw5qOTGTluCXn58X9cR1Lc1xCdAz4xMwc87Zx7BmjmnNsA4JzbYGZNQ4pNRKRC5eUXkL13Hzm5+WTvLfzsI3tvvj9sH9mBcTmHjMtn59597MrNp7CE6R/jl9IkPYUrB7SLa9xhJYiBzrn1fhIYZ2aLYp3RzEYAIwDatGkTr/hERCgocOzKyz/0QL43cCDPPbQ/mAQKp9+5Nz+ms/2UpATqpyWRnppERloy6alJtGlcl4y0ZDLSksgIjKtfJ4nBXbLi/v1DSRDOufX+381m9jbQF9hkZi38q4cWwOYi5n0GeAagT58+1e8B6yJSKfbu23/wYL334MG68MCdszc/cNZ+6Nl8TmECyMunpFamZpCemkR9/6CekZZEk/QU2mXW8w7sqYce3DPSkkhPO3T69LQkUpMSK2fDlEKlJwgzqwckOOey/e5Tgd8C7wFXAX/y/75b2bGJSPj2Fzhycg8/WO8Mnp3vjSyWiRyXT14Mj9ROS04gPTXZO3P3z9Kz0tMPdGcUns0X9vsH9YPTJ1M3OZGEBKuELVP5wriCaAa8bWaF6x/jnPvYzL4C3jCza4HVwEUhxCYiZeScIze/wDuQ7z28mCUnWrFM5Bm7nxhKkmAcOFgXFr80zUijY9ahZ+oHi2aSD5zF1/cP+OmpSaQkVaV2OlVPpScI59wK4Jgow7cCJ1d2PCIC+fsL2JW7/9CzdP9AvjNwxp4TUSRzSLFMbj779pdc6lsnOfFAsUpGWjIZqUk0r592yIE8I+KMPbK/bkoi/kmmxFFYldQiUgGcc+zZt5+cwgN5ScUyh1SqHhy3O6/kdxIkJthhZ+RHNEw7WKkaUSwTWe5e2J+UqLP26kIJQiQk+fsLDpyBB4tlsnP3FXnAP6RYxj+4x3J3bb2UxANn7IVn5C0b1jmkkvTAgbyIA35acoLO2msZJQiRUnLOsTtv/4FimGD5euGBu8himcDZ/J4Y3iSWnGiHHNTTU70De/20jAMH8OKLZbzuxBpaiSrxpQQhtco+/6w9Whl61LP5KDc35QRuWCpO+iHNG5NoWDeFVo3rHtbWPSPt8GKYwv7UJJ21S3iUIKRaKChw7N63P+JGpUPPziOLZYJ3ohZOkxvjDUsZwaKX1GTaNK57WNv1aMUwB87qU5JqbNNHqT2UICTucvP3H9I+vbCiNGdv5Jl7xME9UCyTkxv7DUvBStLG9VJo26TewbbrUQ/ugbP3KnrDkkgYlCCkSAUFjpy8/MMO1tEP7odXshZWqMbymIHUpITD2q5npteNemNStNYxNf2GJZEwKEHUQIU3LBVVhh7Zdn1nlAN+4d+SJBSetQcO1k0z0uiQefCMvH7EQT1au3bdsCRS9ShBVDGFjxk47EmPUR4QdrDM/fAK11hvWAo+QiAjNYlm9Q9t114/8k7UA8N0w5JITacEUUGcc+zdV3DIjUhRnxcTOJDvjPKAsF2luGEpeObevH4anZsmHdLWPXggjzybr5eaRLJuWBKRYihB4N2wFK0pY/CmpOABv7BYJvJsPj+Gto91UxIPO7gf0TCNjNTkw87mgzcwBQ/4umFJRCpDrUwQExZv5ncffHPgbD6WG5aSDjxm4NA7UTPSMoqoNI24gSk1mXqpiXrMgIhUG7UyQTSok0y3FvUDz2kv/k5U3bAkIrVRrUwQvds0ovdljcIOQ0SkSlN5h4iIRKUEISIiUSlBiIhIVEoQIiISlRKEiIhEpQQhIiJRKUGIiEhUShAiIhKVuZLewlKFmdkW4Nsyzp4JfFeB4VQUxVU6iqt0FFfp1NS42jrnskqaqFoniPIwsxnOuT5hxxFJcZWO4iodxVU6tT0uFTGJiEhUShAiIhJVbU4Qz4QdQBEUV+kortJRXKVTq+OqtXUQIiJSvNp8BSEiIsVQghARkahqTIIws2FmttjMlpnZvVHGp5rZ6/74aWbWLjDuPn/4YjM7LdZlhhjXKjObZ2azzWxGZcZlZk3MbIKZ5ZjZYxHzHOvHtczMHrVSvoIvTjFN9Jc52/80LU1M5YxrqJnN9LfJTDM7KTBPubZVHOMKc3v1Dax3jpmdH+syQ4wrtN9iYHwbf9//WazLjJlzrtp/gERgOdABSAHmAN0jprkZeMrvvgR43e/u7k+fCrT3l5MYyzLDiMsftwrIDGl71QMGATcCj0XMMx0YABjwEXB6FYhpItAnpG3VCzjC7/4BsK4itlWc4wpze9UFkvzuFsBmvLdehv1bjBpX2L/FwPi3gH8DP4t1mbF+asoVRF9gmXNuhXMuD3gNODdimnOBF/zuN4GT/bO2c4HXnHO5zrmVwDJ/ebEsM4y4KkKZ43LO7XLOTQH2Bic2sxZAfefc/5y3l74InBdmTBWkPHHNcs6t94cvANL8s8Hybqu4xFXK9ccjrt3OuXx/eBpQ2IIm1N9iMXFVhPIcIzCz84AVeP/H0iwzJjUlQbQE1gT61/rDok7j/7N3AE2KmTeWZYYRF3g76Cd+8cCIUsZU3riKW+baEpZZ2TEV+pdfBPCrMhTlVFRcFwKznHO5lH9bxSuuQqFtLzPrZ2YLgHnAjf74sH+LRcUFIf4Wzawe8HPgN2VYZkySyjJTFRRtJ47M8kVNU9TwaMmztGcO8YgLYKBzbr1fPjzOzBY55yZVUlzlWWZx4hETwHDn3Dozy8C7FL8C74y90uIysx7Aw8CppVhmGHFByNvLOTcN6GFm3YAXzOyjGJdZ6XE55/YS7m/xN8BI51xORB6viO0F1JwriLVA60B/K2B9UdOYWRLQANhWzLyxLDOMuCgsHnDObQbepvRFT+WJq7hltiphmZUdE865df7fbGAMlbytzKwV3v/oSufc8sD05dlW8Yor9O0ViGMhsAuvjiTs32JRcYX9W+wH/NnMVgF3APeb2a0xLjM2Za1cqUofvCuhFXiVuYWVMj0iprmFQyt63vC7e3BoZfAKvEqeEpcZUlz1gAx/mnrAl8CwyoorMP5qDq8Q/groz8GK1zPCjMlfZqbfnYxXfntjJf4PG/rTXxhluWXeVvGKqwpsr/YcrPxti3dQy4xlmSHFVSV+i/7wBzlYSV3u7XVguWWZqSp+gDOAJXi197/wh/0WOMfvTsOr6V+G14KkQ2DeX/jzLSbQmiTaMsOOC69lwhz/syCkuFbhncHk4J2tdPeH9wHm+8t8DP9O/bBi8n+0M4G5/rb6B35LsMqIC/gl3tnm7MCnaUVsq3jEVQW21xX+emcDXwPnVYXfYlFxUQV+i4FlPIifICpqeznn9KgNERGJrqbUQYiISAVTghARkaiUIEREJColCBERiUoJQkREolKCkFrDzG4zs4Vm9oqZXW0RT34tx3JbmNkHfveQwm6///dm9t/yPOvIzJ42s4FFjDvLzCIftSBSIZQgpDa5Ge+GtOEVvNy7gGcjB5rZL4CBeO3mcw+bK3b9gKlFjPsQOMfM6pZj+SJRKUFIrWBmT+Hd2PSemd0ZMW60mf0o0J/j/z3fzD41TwszW2JmzaMs/kLg44hl3o13s9LZzrk9UeKZaGYjzWySf1VznJn9x8yWmtnvA9N1A5Y45/b7V0DfmNlcM3sNwHk3Mk0EzirblhEpWk15WJ9IsZxzN5rZMOBE59x3ZnZ1DPO8bWYX4j3qYBjwgHNuY3AaM2sPfB9xhTAQ6Aoc65zLKWYVec65E8zsduBd4Fi8u8GXm9lI59xW4HQOJp97gfbOuVwzaxhYzgzgeOCNkr6TSGnoCkKkeD8F7gNynXOvRhnfAtgSMWwZ3jOWTj188kO85/+dByxwzm3wE80KDj5s7TQOJoi5wCtmdjmQH1jOZuCIGL6LSKkoQYh4B9sEAP/9BymBcS2BAqCZmUX7vezBe1ZO0Ca84qWRZnZiMestvOooCHQX9if59QoN3cGX+5wJPI53pTHTf7In/voPK8YSKS8lCBHvQX/H+t3n4j3JtPDRyv8CLgMW4lVGR1oCtIsc6JxbAlwAvGxmPcsY14nABD+WBKC1c24CcA/eE1nT/em64D34T6RCKUGIeC2QBpvZdLwWQ7v84fcDk51zk/GSw3V+pfEBzrldeHUGnSIX6pz7CrgGr2K8YxniCtY/JOIlm3nALLwXxWz3x52I15pJpELpaa4i5WRm5+NVSP+ygpf7NdDPObevmGmaAWOccydX5LpFQK2YRMrNb+0Uy7uxS7vc3jFM1ga4u6LXLQK6ghARkSKoDkJERKJSghARkaiUIEREJColCBERiUoJQkREovp/EGw1V+9OUWAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHLZJREFUeJzt3XucHGWd7/HPl4SrXJLAgJBEBiTLKrsKGFkUL2gQIZw1uMIRXygB42bdBXc5sqtRzzlejpfgXri4HtgoSthVAREkh5uGALvLKkgCIVwCJoRAxkAyQMIdXPB3/qhnoNLTM1090z098+z3/Xr1q6ue56mqX1dX/frpp2p6FBGYmVm+tup0AGZm1l5O9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mlrkxk+glrZV0xAB175R0/0jHNJqo8H1JmyT9qk79yZJuHmDZEyX9vP1RDp+kZyTt2+Qyh0vqaVdMFbb/QUnrUuwHDXYst2Bb35B0ejvWPVSS/lLS/E7H8V/ZmEn0g4mIf4+I/Ru1k/QlSf8yEjF1wDuA9wFTIuKQZhaMiB9ExJHtCWvoJN0k6RPlsojYMSLWdCqmWvVirOPvgNNS7He0MZYu4CTgn9L8Fh9wkraRdLmk/5C0c4u2uVeFD9EFwEcl7T7EbbT9vJX0N5LulvS0pAcl/U1NfbekGyU9J+m+2g9qSf9D0qOSnpT0PUnbtjPeZmWR6EcDSeM7HMLewNqIeLbDcVh/ewP3jMB2TgauiYjnaytS4rkcmAAcGRFPtWibM4HrBmsQES8A11J8CI1WoohvInAUcJqkE0r1PwLuAHYFvgBclj5YkfR+YB4wA+gG9gW+PGKRVxERY+IBrAX+GlgBPAlcAmyX6g4HekptPwv8BngauJ/iDTgK+C3wn8AzwJ2p7V7AIuAJYDXwp6X1bA8sBDYBK4HP1GxnbdrWCuBFYDzFG/5A2va9wAdL7U8G/gM4C9gMrAHensrXARuB2YPsg7qxAnOAF4CX02v7cp1lTwZuLs3/LXAzsEudugA+CaxKr/3bgEr1H0/7YxPwM2DvmmX/Ii37NPB/gNcDvwSeAi4FtkltJwJXAb1pXVdRfCMB+Fp6PS+k1/SPpfXvV3p//h54KB0TNwPb13nthwM9wBlpHz8CnFKq35ai1/0wsAE4v289Q4mxZr3PpJifBR4oHTdHpOkLga/WxpqmX5/e64NL7/9jwOEDHB83AB+t87p3AH6e3qt++6fU/kLg/1Ik5WcojtXXAmen134fcFDNMpcDfzLQeVdqdyJw4yDbPofiHHgKWAa8M5XXPW9rlv1wqut7vAjcNMx8cy7wrTT9e2mdO5Xq/x34ZJr+IfD1Ut0M4NEB1tudjodT0uvdRHGuvZUij2yuPY5a8eh4Am9ix68FfpUO9kkUiaZvR5dPjv3TDtyrtGNfn6a/BPxLzXr/NR3c2wEHUpzQM1Ld/FQ/EZiS3ojaRL8cmMqrieH4FONW6QB8Ftgz1Z0MvJTe5HHAVymSy7cpksKR6STZcYB9MFisJ1NK1nWWPZkiEW4FfIfipN+h3rLpQLyKovf3urSdo1LdsRQfMm+g+GD7n8AvapZdBOwMHEBxgiyh6OXsQvHhNzu13RX4EEUi2gn4MfDT0rpuAj5R8zrKif7bqc3ktD/fDmxb57Ufnvb7V4CtKXqhzwETU/3ZKeZJKY7/B3xjqDHW2f4rMZeOm4aJPs3/KcWxvkN6z/5ukO30Am+tWVdvOm4W1ds3NctfSPFB8haKY+wG4EGKnm7f8Xpjqf3Wqf1ODHLepfmDgScG2fZH074eT/GB/CivduS+RM15O8h6dk7768/S/DyK5Fn3McA6RNF778svHwRW1rT5R179ILgT+HCpbrf0nu9aZ93dqe78tI+PpOgo/BTYneJY3gi8eyh5csD90sqVtfORTo5yb+WbwPm1JwewX9pRRwBb16xjiwOGIkG/zJaf1N8ALkzTa4D3l+o+Qf9E//EGcS8HZqXpk4FVpbo/TG/6HqWyx4ED66ynUawn0zjR30rxTegnpF51vWVTTO8ozV8KzEvT1wJzSnVbUSTNvUvLHlaqXwZ8tjT/98DZA8R4ILCpNH8TAyT6tN3ngTdXOHYOT23Hl8o2AodSnNTPsmVSehvw4FBjrLPMkBN9KlsE3EXR0RgwWVP0en+/Zl0vUPSIP1RhP10IfKc0/ylKCS4dr5tL8zOAJY3Ou1Q/DXi5UQyl9pv63lsqJvp0TFwFnFd1OwOs58sUyXvbNP8x4JaaNl/j1XPvAVJHKM1vnd7z7jrr7k51k0tlj7PlB8VPgNOH8xpqH2NtjP7R0vRzwI61DSJiNXA6xcGxUdLFkvYaYH17UfQyni6VPUTxqdpXv65UV56uWybpJEnLJW2WtBn4A4pP+D4bStPPp5hry/q9rgqxVrEfMItiaOe3DdoOtK/3Bs4pvb4nKJJlOY7a11P39UnaQdI/SXpI0lPAvwETJI2r8Fp2o+gRPVChLcDjEfFSndfURdFbXlZ6Tdel8uHG2CrfoTiOvhURLw7SbhNF77rsMeAEYGEaS26k0nuXzASugUrn3U4Uw2t1STpD0sp0MXMzxbe/3QZqP4Cvpe38ZZPLleM4jeIbzDGlff0MxTeFsp0pvn3Xq++bfpqBNbOfh22sJfpKIuKHEfEOiqQUwJl9VTVN1wOTJJVPjtdRjDNCMZY7pVQ3td7m+iYk7U1xUp5G8bVtAnA3RSIcrkaxVrGSYtjoWkkN71IawDqKr8UTSo/tI+IXQ1jXGRRf+f8oInYG3pXK+/ZX7ftV9hhFb/X1Q9hu7XqeBw4ovZ5dIqLvRBtOjFU8S/FB0+e15UpJO1IMLV0AfEnSpEHWtYJiPHkLEXE5xRDQZZLeM8x4y2YCV5e2M9B5B8VQ3531ViLpnRTj+/+dYjhtAsWHQuV9nC6cfgQ4LiL+s1T++XRba91HzTo+TrqoGhHlO4nuAfatOffezKsX2O9J8+W6DRHxeKO4R0p2iV7S/pLem+4yeIHiJH45VW8AuiVtBRAR64BfAN+QtJ2kN1Fc2PxBan8p8DlJEyVNpkjgg3kNxUHZm2I5haInNmwVYq26nh8BnweulzSUJHk+xT45AEDSLpKOH8J6oOh9PQ9sTgnsizX1GyjG9vuJiN8B3wP+Id3iN07S25q9rS2t5zvAWX23/0maXOr9DjnGipYDMyVNkvRail5x2TnAsoj4BEVSPX+QdV0DvLteRXrfTwOulHTYMOIFQNI+FEMb96X5wc47UlzXDrC6nSiuofQC4yX9b7bsIW9x3taJ5SDgW8CxEdFbrouIr0dxW2vdR2kdJwJfB94XNbfvRsSvKd6nL6Zz74PAmyiGWAAuAuZIeqOkiRTXrS4c4LV2RHaJnuKi5nyKntqjFBc4Pp/qfpyeH5d0e5r+CMW42XrgCuCLEbE41X2F4q6FB4HrgcsoLi7WFRH3UoxB/5Li4PxDijsXWmWwWCuLiIUUr+0GSd1NLnsFRU/t4jSUcTdwdLMxJGdT3DnzGHAL/W/TOwc4TsUfgZ1bZ/m/phi7vo1iCOlMhnZMf5biAvMt6TVdT9GLb0WMjfwzRU93LcWdMZf0VUiaRXHXySdT0aeBg1NSquciig+N7etVpvf9DOBqSU39rUUdx5CGbZIBzztJ21H0/hcOsK6fUXwI/JpiOPIFthwSrXfels2iuGHi5lJvfaAPlYF8leJi8G2ldZQ/VE8AplMMj82n+ObQCxAR11FcM7wxxf8Q/TsEHaU0+G8VSPpz4ISIqNtrMus0SV8HNkbE2W3ezjUUtwFeU6Htp4CpEfGZdsZkA3OiH4SkPSm+lv+S4q6BqykO7raeRGajnaTPUFwc7vfHWTb6ONEPIl1cvRrYh+K+24uBz1W4Y8XMbNRwojczy1yOF2PNzKyk0z/EBcBuu+0W3d3dnQ7DzGxMWbZs2WMR0dWo3ahI9N3d3SxdurTTYZiZjSmSHqrSzkM3ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmRsVfxk7VnXPu7pxozZYO/+YjmzXzMYm9+jNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpa5hole0v6SlpceT0k6XdIkSYslrUrPE1N7STpX0mpJKyQd3P6XYWZmA2mY6CPi/og4MCIOBN4CPAdcAcwDlkTENGBJmgc4GpiWHnOB89oRuJmZVdPs0M0M4IGIeAiYBSxM5QuBY9P0LOCiKNwCTJC0Z0uiNTOzpjWb6E8AfpSm94iIRwDS8+6pfDKwrrRMTyrbgqS5kpZKWtrb29tkGGZmVlXlRC9pG+ADwI8bNa1TFv0KIhZExPSImN7V1VU1DDMza1IzPfqjgdsjYkOa39A3JJOeN6byHmBqabkpwPrhBmpmZkPTTKL/CK8O2wAsAman6dnAlaXyk9LdN4cCT/YN8ZiZ2cir9B+mJO0AvA/4s1LxfOBSSXOAh4HjU/k1wExgNcUdOqe0LFozM2tapUQfEc8Bu9aUPU5xF05t2wBObUl0ZmY2bP7LWDOzzDnRm5llrtLQjY0u3fOu7ti2184/pmPbNrOhcY/ezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWWuUqKXNEHSZZLuk7RS0tskTZK0WNKq9DwxtZWkcyWtlrRC0sHtfQlmZjaYqj36c4DrIuL3gTcDK4F5wJKImAYsSfMARwPT0mMucF5LIzYzs6Y0TPSSdgbeBVwAEBG/jYjNwCxgYWq2EDg2Tc8CLorCLcAESXu2PHIzM6ukSo9+X6AX+L6kOyR9V9JrgD0i4hGA9Lx7aj8ZWFdavieVbUHSXElLJS3t7e0d1oswM7OBVUn044GDgfMi4iDgWV4dpqlHdcqiX0HEgoiYHhHTu7q6KgVrZmbNq5Loe4CeiLg1zV9Gkfg39A3JpOeNpfZTS8tPAda3JlwzM2tWw0QfEY8C6yTtn4pmAPcCi4DZqWw2cGWaXgSclO6+ORR4sm+Ix8zMRt74iu0+BfxA0jbAGuAUig+JSyXNAR4Gjk9trwFmAquB51JbMzPrkEqJPiKWA9PrVM2o0zaAU4cZl5mZtYj/MtbMLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZq5ToJa2VdJek5ZKWprJJkhZLWpWeJ6ZySTpX0mpJKyQd3M4XYGZmg2umR/+eiDgwIvr+d+w8YElETAOWpHmAo4Fp6TEXOK9VwZqZWfOGM3QzC1iYphcCx5bKL4rCLcAESXsOYztmZjYMVRN9AD+XtEzS3FS2R0Q8ApCed0/lk4F1pWV7UtkWJM2VtFTS0t7e3qFFb2ZmDY2v2O6wiFgvaXdgsaT7BmmrOmXRryBiAbAAYPr06f3qzcysNSr16CNifXreCFwBHAJs6BuSSc8bU/MeYGpp8SnA+lYFbGZmzWmY6CW9RtJOfdPAkcDdwCJgdmo2G7gyTS8CTkp33xwKPNk3xGNmZiOvytDNHsAVkvra/zAirpN0G3CppDnAw8Dxqf01wExgNfAccErLozYzs8oaJvqIWAO8uU7548CMOuUBnNqS6MzMbNj8l7FmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmKid6SeMk3SHpqjS/j6RbJa2SdImkbVL5tml+darvbk/oZmZWRTM9+r8CVpbmzwTOiohpwCZgTiqfA2yKiP2As1I7MzPrkEqJXtIU4Bjgu2lewHuBy1KThcCxaXpWmifVz0jtzcysA6r26M8GPgP8Ls3vCmyOiJfSfA8wOU1PBtYBpPonU/stSJoraamkpb29vUMM38zMGmmY6CX9N2BjRCwrF9dpGhXqXi2IWBAR0yNieldXV6VgzcyseeMrtDkM+ICkmcB2wM4UPfwJksanXvsUYH1q3wNMBXokjQd2AZ5oeeRmZlZJwx59RHwuIqZERDdwAnBDRJwI3Agcl5rNBq5M04vSPKn+hojo16M3M7ORMZz76D8LfFrSaoox+AtS+QXArqn808C84YVoZmbDUWXo5hURcRNwU5peAxxSp80LwPEtiM3MzFrAfxlrZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mljknejOzzDVM9JK2k/QrSXdKukfSl1P5PpJulbRK0iWStknl26b51am+u70vwczMBlPlf8a+CLw3Ip6RtDVws6RrKf7x91kRcbGk84E5wHnpeVNE7CfpBOBM4MNtit9GWPe8qzuy3bXzj+nIds1y0LBHH4Vn0uzW6RHAe4HLUvlC4Ng0PSvNk+pnSFLLIjYzs6ZUGqOXNE7ScmAjsBh4ANgcES+lJj3A5DQ9GVgHkOqfBHats865kpZKWtrb2zu8V2FmZgOqlOgj4uWIOBCYAhwCvKFes/Rcr/ce/QoiFkTE9IiY3tXVVTVeMzNrUlN33UTEZuAm4FBggqS+Mf4pwPo03QNMBUj1uwBPtCJYMzNrXpW7brokTUjT2wNHACuBG4HjUrPZwJVpelGaJ9XfEBH9evRmZjYyqtx1syewUNI4ig+GSyPiKkn3AhdL+ipwB3BBan8B8M+SVlP05E9oQ9xmZlZRw0QfESuAg+qUr6EYr68tfwE4viXRmZnZsPkvY83MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w1TPSSpkq6UdJKSfdI+qtUPknSYkmr0vPEVC5J50paLWmFpIPb/SLMzGxgVXr0LwFnRMQbgEOBUyW9EZgHLImIacCSNA9wNDAtPeYC57U8ajMzq6xhoo+IRyLi9jT9NLASmAzMAhamZguBY9P0LOCiKNwCTJC0Z8sjNzOzSpoao5fUDRwE3ArsERGPQPFhAOyemk0G1pUW60llteuaK2mppKW9vb3NR25mZpVUTvSSdgR+ApweEU8N1rROWfQriFgQEdMjYnpXV1fVMMzMrEmVEr2krSmS/A8i4vJUvKFvSCY9b0zlPcDU0uJTgPWtCdfMzJpV5a4bARcAKyPiH0pVi4DZaXo2cGWp/KR0982hwJN9QzxmZjbyxldocxjwMeAuSctT2eeB+cClkuYADwPHp7prgJnAauA54JSWRmxmZk1pmOgj4mbqj7sDzKjTPoBThxmXmZm1iP8y1swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmqvx6pVnHdc+7uiPbXTv/mI5s16yV3KM3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMVfnn4N+TtFHS3aWySZIWS1qVniemckk6V9JqSSskHdzO4M3MrLEqPfoLgaNqyuYBSyJiGrAkzQMcDUxLj7nAea0J08zMhqphoo+IfwOeqCmeBSxM0wuBY0vlF0XhFmCCpD1bFayZmTVvqGP0e0TEIwDpefdUPhlYV2rXk8r6kTRX0lJJS3t7e4cYhpmZNdLqi7GqUxb1GkbEgoiYHhHTu7q6WhyGmZn1GWqi39A3JJOeN6byHmBqqd0UYP3QwzMzs+EaaqJfBMxO07OBK0vlJ6W7bw4Fnuwb4jEzs85o+KNmkn4EHA7sJqkH+CIwH7hU0hzgYeD41PwaYCawGngOOKUNMZuZWRMaJvqI+MgAVTPqtA3g1OEGZWZmreO/jDUzy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swscw3vozf7r6x73tUd2/ba+cd0bNuWF/fozcwy50RvZpY5J3ozs8x5jN5slOrU9QFfG8iPe/RmZplzojczy5wTvZlZ5pzozcwyN+YvxnbyD1rMzMYC9+jNzDI35nv0ZtZa/tmH/LSlRy/pKEn3S1otaV47tmFmZtW0vEcvaRzwbeB9QA9wm6RFEXFvq7dlZnnxH4m1RzuGbg4BVkfEGgBJFwOzACd6MxuVch+uakeinwysK833AH9U20jSXGBumn1G0v1tiAVgN+CxNq27nRz3yBmLMYPjHklti1lnDmvxvas0akeiV52y6FcQsQBY0IbtbxmMtDQiprd7O63muEfOWIwZHPdIGosxl7XjYmwPMLU0PwVY34btmJlZBe1I9LcB0yTtI2kb4ARgURu2Y2ZmFbR86CYiXpJ0GvAzYBzwvYi4p9XbaULbh4faxHGPnLEYMzjukTQWY36FIvoNn5uZWUb8EwhmZplzojczy1wWiV7SJEmLJa1KzxMHaDc7tVklaXap/Kb0kw3L02P3Nsc76E9ESNpW0iWp/lZJ3aW6z6Xy+yW9v51xtiJmSd2Sni/t2/NHKuaKcb9L0u2SXpJ0XE1d3eNllMf8cmlfj+hNEBXi/rSkeyWtkLRE0t6luo7s6xbE3bH93ZSIGPMP4JvAvDQ9DzizTptJwJr0PDFNT0x1NwHTRyjWccADwL7ANsCdwBtr2vwFcH6aPgG4JE2/MbXfFtgnrWfcKI+5G7i7Q8dFlbi7gTcBFwHHVTleRmvMqe6ZUbyv3wPskKb/vHSMdGRfDzfuTu7vZh9Z9OgpfmJhYZpeCBxbp837gcUR8UREbAIWA0eNUHxlr/xERET8Fuj7iYiy8uu5DJghSan84oh4MSIeBFan9Y3mmDupYdwRsTYiVgC/q1m2U8fLcGLupCpx3xgRz6XZWyj+xgY6e24OJ+4xI5dEv0dEPAKQnusNvdT7aYbJpfnvp69f/6vNCapRHFu0iYiXgCeBXSsu2w7DiRlgH0l3SPpXSe9sd7D1Ykqa2V+jeV8PZjtJSyXdIqleh6ddmo17DnDtEJdtpeHEDZ3b300ZM79HL+l64LV1qr5QdRV1yvruLT0xIn4jaSfgJ8DHKL4Wt0OVn4gYqE2ln5dog+HE/Ajwuoh4XNJbgJ9KOiAinmp1kHUMZ3+N5n09mNdFxHpJ+wI3SLorIh5oUWyDqRy3pI8C04F3N7tsGwwnbujc/m7KmOnRR8QREfEHdR5XAhsk7QmQnjfWWcWAP80QEb9Jz08DP6S9wyFVfiLilTaSxgO7AE9UXLYdhhxzGmZ6HCAillGMh/5e2yOuiSlpZn+N5n09oIjoO6bXUFx7OqiVwQ2iUtySjqDonH0gIl5sZtk2GU7cndzfzen0RYJWPIC/ZcuLsd+s02YS8CDFxZ6JaXoSxbea3VKbrSnGlz/ZxljHU1xs2odXL/4cUNPmVLa8sHlpmj6ALS/GrmFkLsYOJ+auvhgpLnj9Bpg0QsdFw7hLbS+k/8XYfsfLKI95IrBtmt4NWEXNhcUOHyMHUXzQT6sp78i+bkHcHdvfTb/OTgfQojdrV2BJ2tFL+g4Siq9Z3y21+zjFBczVwCmp7DXAMmAFcA9wTruTJzAT+HU6eL6Qyr5C0VsA2A74cYrzV8C+pWW/kJa7Hzh6BPfxkGIGPpT2653A7cAfj/Cx0Sjut1L06p4FHgfuGex4Gc0xA28H7kr7+i5gzijb19cDG4Dl6bGo0/t6OHF3en838/BPIJiZZW7MjNGbmdnQONGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDL3/wFUl2/eUOGUvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGVxJREFUeJzt3X+cXXV95/HX2wCBByJJyIAhSRnQlB+6NdDZmEptKbGWH/swuIKFdiXQtBEXdnW11eB2C7qlhnUVy7biBoOErRCQSomAPyI/pRbogCEkBEqA0QyJyUhCIKJpEz77x/mOHCd35p6Ze+/cmS/v5+NxH/ec7/d77v2cOzPve+Z7z71XEYGZmeXrde0uwMzMWstBb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQd9hiT1SHrXIH3vlPTkaNc0lqjwFUnbJT3UphpC0pvbcd/22uOgf42JiO9FxNH1xkm6VNLfjUZNbfCbwO8CMyJiTruLkXStpPNGuO2QTxiSTpd0v6QXJP1Y0tWSDir1T5R0jaQXU/9HB2w/T9ITkl6WdLekI0ZSp7WXg95GnaR92lzCEUBPRPy0zXWMhoOBvwQOB44FZgCfLfVfCsyieEx+B/i4pFMAJE0Fvg78D2AK0A3cOFqFWxNFhC+ZXYAe4E+BNcAOij/O/VPfSUBvaewngOeAl4AngXnAKcC/Av8G7AQeTWMPB1YC24ANwJ+UbucAYDmwHVgPfHzA/fSk+1oD7AL2ARYDT6f7fhx4b2n8ecA/AlcALwDPAO9I7RuBrcCCIR6DmrUCC4GfA3vSvn2qxrY/BH49Lf8nIIDj0vofA/+Qbv9nwJTSdscDPwH2Tet/lB6L7cC3gSNKYwN4c1q+FjhvkP2YA/xTegw2A38D7Jf67ku389O0L79f4XfjPwKPldafA95dWv+fwIq0vAj4fqnvwLTPxwzxe/dn6Wf8U2AZcBjwzfQz/i4wud1/H6/FS9sL8KUFP9TiD+6hFEZTUthckPpO6g9g4OgUmoen9U7gTWn5UuDvBtzuvcAXgf2B2UAfMC/1LUn9kymOGtewd9CvBmYCB6S2s1KNrwN+P4XDtNR3HrAbOB+YQHFU+iPgb4GJwLtTeLx+kMdgqFrPA+4f4vG7DvhYWl5K8WT0oVLff0vLd/HLT3afBb6Uls+geII5luJJ7c8HhOYvgr7Oz/LXgbnpNjrTz/Ijw72d0vgv8GqQT07bH1bqP5P0RAD8NXDVgO3XAu8b4vfuAYpwn07xZPwIxRPgxPR4XdLuv4/X4sVTN/m6MiI2RcQ24BsUYTfQHoo/wOMk7RsRPRHxdK0bkzSTYm77ExHx84hYDXwZ+EAa8n7gryJie0T0AlcOUtPGiPgZQER8LdX4SkTcCDxFcQTb79mI+EpE7KH4r2Qm8OmI2BUR36H4r2Ov+ekKtdZzL/DbafmdwGdK67+d+gGuB85J9yng7NQG8EHgMxGxPiJ2A38FzB7uHHdEPBwRD0TE7ojoAf5vqZZhkfS7wALgL1LT69P1jtKwHcBBpf5y38D+Wv5PRGyJiOeA7wEPRsQPImIXcAtF6Nsoc9Dn68el5Zd59Y/6FyJiA/ARiqP3rZJWSDp8kNs7HNgWES+V2n5IceTW37+x1Fdertkm6VxJq9MLhS8AbwWmloZsKS33PzkMbNtrvyrUWs+9wDslvZHiv4kbgRMldVLMea9O424GfiM9Zr9FcXT8vdR3BPDXpX3bBmgYNQAg6Vcl3ZZeKH2R4gljar3tatzOXIonoTMj4l9S8850/YbS0DdQ/KfU31/uG9hfy8CfT5Wfl7WYg/41LiKuj4jfpAimAC7v7xowdBMwpXzGBvArFHO8UMwfzyj1zax1d/0L6cj2auAi4JCImEQxLaAR7spwah1SegJ8GfivwH3pCePHFHPW90fEK2ncC8B3KP6b+QPghojo38eNwAcjYlLpckBEfH+Y+3IV8AQwKyLeAHySYT5Gko6neL3ijyLiztJ+bqf4ub2tNPxtwLq0vK7cJ+lA4E2lfhsnHPSvYZKOlnSypIkUL1D+jGI6B4ojsU5JrwOIiI3A94HPSNpf0q9RvLD51TT+JuBiSZMlTacI8KEcSBH8famW8ymO6BtWodYq7qXYh/5pmnsGrPe7HjgXeB+vTtsAfIni8XgLgKSDJZ01zF2BYprkRWCnpGOADw3o3wIcNdjGkt4KfAv4LxHxjRpDrgP+PP3cjgH+hOLFYSimWt4q6X2S9qeY8lkTEU+MYD+sjRz0r20TKV5E/QnFEeuhFEeMAF9L189LeiQtn0PxguAmihC4JCJWpb5PA73AsxRnV9xMcXZNTRHxOPA5ijNKtgD/juIsm2YZqtYq7qUI2fsGWe+3kuL0xC0R8Wh/Y0TcQvHf0Yo05bIWOHX4u8GfUvy38BLFf0ADT2+8FFiepojeX2P7jwEdwDJJO9OlfER+CcWLzT9M+/jZiPhW2oc+iiewyyjOHHo7xesQNs7o1f80zZpH0oeAsyNiRC8cmlnz+IjemkLSNEknSnqdpKMpjiRvaXddZlacm2vWDPtRnPp3JMWbe1ZQnMduZm3mqRszs8x56sbMLHNjYupm6tSp0dnZ2e4yzMzGlYcffvgnEdFRb9yYCPrOzk66u7vbXYaZ2bgi6YdVxnnqxswscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swsc2PinbE2fnQuvr0t99uz5PS23K9ZDnxEb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmKge9pAmSfiDptrR+pKQHJT0l6UZJ+6X2iWl9Q+rvbE3pZmZWxXCO6D8MrC+tXw5cERGzgO3AwtS+ENgeEW8GrkjjzMysTSoFvaQZwOnAl9O6gJOBm9OQ5cAZaXl+Wif1z0vjzcysDaoe0X8B+DjwSlo/BHghInan9V5gelqeDmwESP070ngzM2uDukEv6T8AWyPi4XJzjaFRoa98u4skdUvq7uvrq1SsmZkNX5Uj+hOB90jqAVZQTNl8AZgkqf+zcmYAm9JyLzATIPUfDGwbeKMRsTQiuiKiq6Ojo6GdMDOzwdUN+oi4OCJmREQncDZwV0T8IXA3cGYatgC4NS2vTOuk/rsiYq8jejMzGx2NnEf/CeCjkjZQzMEvS+3LgENS+0eBxY2VaGZmjRjWxxRHxD3APWn5GWBOjTE/B85qQm1mZtYEfmesmVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZpkb1umVNjZ0Lr693SWY2TjiI3ozs8w56M3MMuegNzPLnIPezCxzDnozs8w56M3MMuegNzPLnIPezCxzDnozs8xV+XLw/SU9JOlRSeskfSq1XyvpWUmr02V2apekKyVtkLRG0gmt3gkzMxtclY9A2AWcHBE7Je0L3C/pm6nvzyLi5gHjTwVmpcvbgavStZmZtUGVLwePiNiZVvdNl6G+7Hs+cF3a7gFgkqRpjZdqZmYjUWmOXtIESauBrcCqiHgwdV2WpmeukDQxtU0HNpY2701tA29zkaRuSd19fX0N7IKZmQ2lUtBHxJ6ImA3MAOZIeitwMXAM8O+BKcAn0nDVuokat7k0Iroioqujo2NExZuZWX3DOusmIl4A7gFOiYjNaXpmF/AVYE4a1gvMLG02A9jUhFrNzGwEqpx10yFpUlo+AHgX8ET/vLskAWcAa9MmK4Fz09k3c4EdEbG5JdWbmVldVc66mQYslzSB4onhpoi4TdJdkjoopmpWAxek8XcApwEbgJeB85tftpmZVVU36CNiDXB8jfaTBxkfwIWNl2ZmZs3gd8aamWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZq/LOWLO261x8e9vuu2fJ6W27b7Nm8BG9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZpmr8g1T+0t6SNKjktZJ+lRqP1LSg5KeknSjpP1S+8S0viH1d7Z2F8zMbChVjuh3ASdHxNuA2cAp6SsCLweuiIhZwHZgYRq/ENgeEW8GrkjjzMysTeoGffoC8J1pdd90CeBk4ObUvpzie2MB5qd1Uv+89L2yZmbWBpXm6CVNkLQa2AqsAp4GXoiI3WlILzA9LU8HNgKk/h3AIc0s2szMqqsU9BGxJyJmAzOAOcCxtYal61pH7zGwQdIiSd2Suvv6+qrWa2ZmwzSss24i4gXgHmAuMElS/2flzAA2peVeYCZA6j8Y2FbjtpZGRFdEdHV0dIysejMzq6vKWTcdkial5QOAdwHrgbuBM9OwBcCtaXllWif13xURex3Rm5nZ6Kjy6ZXTgOWSJlA8MdwUEbdJehxYIekvgR8Ay9L4ZcD/k7SB4kj+7BbUbWZmFdUN+ohYAxxfo/0Zivn6ge0/B85qSnVmZtYwvzPWzCxzDnozs8w56M3MMuegNzPLnIPezCxzDnozs8w56M3MMuegNzPLnIPezCxzDnozs8w56M3MMuegNzPLnIPezCxzDnozs8w56M3MMuegNzPLXJWvEpwp6W5J6yWtk/Th1H6ppOckrU6X00rbXCxpg6QnJf1eK3fAzMyGVuWrBHcDH4uIRyQdBDwsaVXquyIi/nd5sKTjKL4+8C3A4cB3Jf1qROxpZuFmZlZNla8S3AxsTssvSVoPTB9ik/nAiojYBTybvjt2DvBPTajXbNR1Lr69Lffbs+T0ttyv5WdYc/SSOim+P/bB1HSRpDWSrpE0ObVNBzaWNuulxhODpEWSuiV19/X1DbtwMzOrpnLQS3o98PfARyLiReAq4E3AbIoj/s/1D62xeezVELE0Iroioqujo2PYhZuZWTVV5uiRtC9FyH81Ir4OEBFbSv1XA7el1V5gZmnzGcCmplRr9hriKSNrlipn3QhYBqyPiM+X2qeVhr0XWJuWVwJnS5oo6UhgFvBQ80o2M7PhqHJEfyLwAeAxSatT2yeBcyTNppiW6QE+CBAR6yTdBDxOccbOhT7jxsysfaqcdXM/tefd7xhim8uAyxqoy8zMmsTvjDUzy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMVfmGqZmS7pa0XtI6SR9O7VMkrZL0VLqenNol6UpJG9IXh5/Q6p0wM7PBVTmi3w18LCKOBeYCF0o6DlgM3BkRs4A70zrAqRRfHzgLWETxJeJmZtYmdYM+IjZHxCNp+SVgPTAdmA8sT8OWA2ek5fnAdVF4AJg04PtlzcxsFA1rjl5SJ3A88CBwWERshuLJADg0DZsObCxt1pvazMysDSoHvaTXA38PfCQiXhxqaI22qHF7iyR1S+ru6+urWoaZmQ1TpaCXtC9FyH81Ir6emrf0T8mk662pvReYWdp8BrBp4G1GxNKI6IqIro6OjpHWb2ZmdVQ560bAMmB9RHy+1LUSWJCWFwC3ltrPTWffzAV29E/xmJnZ6NunwpgTgQ8Aj0lando+CSwBbpK0EPgRcFbquwM4DdgAvAyc39SKzcxsWOoGfUTcT+15d4B5NcYHcGGDdZmZWZP4nbFmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmqnyV4DWStkpaW2q7VNJzklany2mlvoslbZD0pKTfa1XhZmZWTZUj+muBU2q0XxERs9PlDgBJxwFnA29J23xR0oRmFWtmZsNXN+gj4j5gW8Xbmw+siIhdEfEsxffGzmmgPjMza1Ajc/QXSVqTpnYmp7bpwMbSmN7UthdJiyR1S+ru6+troAwzMxvKSIP+KuBNwGxgM/C51F7rS8Sj1g1ExNKI6IqIro6OjhGWYWZm9Ywo6CNiS0TsiYhXgKt5dXqmF5hZGjoD2NRYiWZm1oh9RrKRpGkRsTmtvhfoPyNnJXC9pM8DhwOzgIcartLMRk3n4tvbdt89S05v233nrG7QS7oBOAmYKqkXuAQ4SdJsimmZHuCDABGxTtJNwOPAbuDCiNjTmtLNzKyKukEfEefUaF42xPjLgMsaKcrMzJrH74w1M8ucg97MLHMOejOzzDnozcwy56A3M8ucg97MLHMjesOUFdr5xhIzs6p8RG9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mlrm6QS/pGklbJa0ttU2RtErSU+l6cmqXpCslbZC0RtIJrSzezMzqq3JEfy1wyoC2xcCdETELuDOtA5xK8T2xs4BFwFXNKdPMzEaqbtBHxH3AtgHN84HlaXk5cEap/booPABMkjStWcWamdnwjXSO/rCI2AyQrg9N7dOBjaVxvaltL5IWSeqW1N3X1zfCMszMrJ5mvxirGm1Ra2BELI2Irojo6ujoaHIZZmbWb6RBv6V/SiZdb03tvcDM0rgZwKaRl2dmZo0aadCvBBak5QXAraX2c9PZN3OBHf1TPGZm1h51v3hE0g3AScBUSb3AJcAS4CZJC4EfAWel4XcApwEbgJeB81tQs5mZDUPdoI+IcwbpmldjbAAXNlqUmZk1j98Za2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZa7uh5qZmY2WzsW3t+V+e5ac3pb7HS0+ojczy5yD3swscw1N3UjqAV4C9gC7I6JL0hTgRqAT6AHeHxHbGyvTzMxGqhlH9L8TEbMjoiutLwbujIhZwJ1p3czM2qQVUzfzgeVpeTlwRgvuw8zMKmo06AP4jqSHJS1KbYf1fyF4uj601oaSFknqltTd19fXYBlmZjaYRk+vPDEiNkk6FFgl6YmqG0bEUmApQFdXVzRYh5mZDaKhI/qI2JSutwK3AHOALZKmAaTrrY0WaWZmIzfioJd0oKSD+peBdwNrgZXAgjRsAXBro0WamdnINTJ1cxhwi6T+27k+Ir4l6Z+BmyQtBH4EnNV4mWZmNlIjDvqIeAZ4W43254F5jRRlZmbNM+4/66Zdn41hZjZe+CMQzMwy56A3M8ucg97MLHMOejOzzDnozcwy56A3M8ucg97MLHMOejOzzDnozcwy56A3M8ucg97MLHMOejOzzDnozcwyN+4/vdLMrFHt/BTcniWnt/w+fERvZpa5lgW9pFMkPSlpg6TFrbofMzMbWkuCXtIE4G+BU4HjgHMkHdeK+zIzs6G16oh+DrAhIp6JiH8FVgDzW3RfZmY2hFa9GDsd2Fha7wXeXh4gaRGwKK3ulPQ88JMW1TMapjK+64fxvw/jvX4Y//sw3uuHUd4HXd7Q5kdUGdSqoFeNtvillYilwNJfbCB1R0RXi+ppufFeP4z/fRjv9cP434fxXj/ksQ8DtWrqpheYWVqfAWxq0X2ZmdkQWhX0/wzMknSkpP2As4GVLbovMzMbQkumbiJit6SLgG8DE4BrImJdnc2W1ukf68Z7/TD+92G81w/jfx/Ge/2Qxz78EkVE/VFmZjZu+Z2xZmaZc9CbmWWubUEv6SxJ6yS9ImnQU5kk9Uh6TNJqSd2jWeNQhlH/mP0oCElTJK2S9FS6njzIuD3p8V8tqe0vqtd7TCVNlHRj6n9QUufoVzm0CvtwnqS+0uP+x+2oczCSrpG0VdLaQfol6cq0f2sknTDaNQ6lQv0nSdpRevz/YrRrbKqIaMsFOBY4GrgH6BpiXA8wtV11NlI/xQvRTwNHAfsBjwLHtbv2Un3/C1iclhcDlw8ybme7ax3OYwr8Z+BLafls4MZ21z2CfTgP+Jt21zrEPvwWcAKwdpD+04BvUrynZi7wYLtrHmb9JwG3tbvOZl3adkQfEesj4sl23X+jKtY/1j8KYj6wPC0vB85oYy1VVXlMy/t1MzBPUq038bXLWP+9qCsi7gO2DTFkPnBdFB4AJkmaNjrV1Veh/qyMhzn6AL4j6eH0sQnjSa2PgpjeplpqOSwiNgOk60MHGbe/pG5JD0hq95NBlcf0F2MiYjewAzhkVKqrpurvxfvStMfNkmbW6B/LxvrvfhW/IelRSd+U9JZ2F9OIln7xiKTvAm+s0fXfI+LWijdzYkRsknQosErSE+nZuOWaUH/dj4JotaH2YRg38yvpZ3AUcJekxyLi6eZUOGxVHtO2P+51VKnvG8ANEbFL0gUU/6Gc3PLKmmes/wzqeQQ4IiJ2SjoN+AdgVptrGrGWBn1EvKsJt7EpXW+VdAvFv72jEvRNqL/tHwUx1D5I2iJpWkRsTv9Wbx3kNvp/Bs9Iugc4nmKOuR2qPKb9Y3ol7QMczNj6N73uPkTE86XVq4HGPvpq9LX9d78REfFiafkOSV+UNDUixuUHto3pqRtJB0o6qH8ZeDdQ81XyMWqsfxTESmBBWl4A7PVfiqTJkiam5anAicDjo1bh3qo8puX9OhO4K9IrbGNE3X0YMJ/9HmD9KNbXDCuBc9PZN3OBHf3ThOOBpDf2v64jaQ5FVj4/9FZjWLteBQbeS/GsvwvYAnw7tR8O3JGWj6I4I+FRYB3FlEnbX8GuWn9aPw34F4oj4DFTf6rtEOBO4Kl0PSW1dwFfTsvvAB5LP4PHgIVjoO69HlPg08B70vL+wNeADcBDwFHtrnkE+/CZ9Dv/KHA3cEy7ax5Q/w3AZuDf0t/BQuAC4ILUL4ovH3o6/d4MembdGK3/otLj/wDwjnbX3MjFH4FgZpa5MT11Y2ZmjXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpa5/w8zLPr6b/YiUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFMhJREFUeJzt3X+QZWV95/H3R4YfFqiADC7CaENEK7irqLOG0rhLiesvNLC1usG4ii4V1h+pMlVxdYybClImwlYlZK0yGlaj+BPUGCGSWpfwI3GzUTMookgRBhhlHAKDMAiKKPLdP87TtZee7unb0/fOnXl4v6pu3XOe5znnfO/p7k8/99zbt1NVSJL69ahZFyBJmi6DXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwZ9p5JsTvKiJfpekOSG3V2TpNkw6B+BquorVfW05cYlOSvJJ3dHTbvbco8tyX0jt4eS3D+y/tq2/Vm7eOwlfwlPQpLDk3wmydYk9yT5+yS/smDMbyT5XpIfJ/likkNH+g5N8pet73tJfmNatWr3MOg1E0nWzLqGnamqg+ZvwPeBV460fWrW9S3jIOAfgecAhwIXAJcmOQggydOBPwNeBzwB+AnwpyPbfwD4Wet7LfDBto32VlXlrcMbsBl4O3AtcA9wEXBA6zsR2DIy9p3AD4B7gRuAk4CXMvyw/xy4D/hWG/tE4BLgLmAT8Jsj+3k0Q6jcDVwPvGPBcTa3Y10LPACsATYAN7Vjfxf49yPj3wD8PXAesB24GXhea78VuAM4fSfnYNFal3psy5zLFy1oOws4a4nxvwRcAfwQuBP4FHBw6/sE8BBwfzv2OxbZ/q9a3/ztIeANq/x++BHwnLb8h8CnF9T7M+AxwIFt+akj/Z8Azlliv2cBnwM+2b6G3waeCryrfX1uBV4865+HR/pt5gV4m9IXdginr7ewO7QF75ta34nzAQw8rf0wPrGtzwG/1JbPAj65YL9/yzD7OwA4HtgGnNT6zmn9hwBHtUBfGPTXAOuAR7e2V7caHwX8OvBj4IjW9wbgQeCNwD7Aexlm1x8A9gde3MLloCXOwc5q3eGxLXMuXzTO2Db+KcC/azWuBf4O+JNd2R/DL6WtwLq2fi3DL73Fbn+6xD6OB34KPK6tXwy8c8GY+xieATwLuH9B39uBv1pi32e1fb+E4Rf3x4FbgHcD+wK/Cdwy65+HR/rNSzd9e39Vba2quxhmiccvMuYXDIF0XJJ9q2pzVd202M6SrAN+lSEkflpV1wAfZrgEAPAfgT+sqruragvw/iVqurWq7geoqs+1Gh+qqouAG4Hnjoy/pao+WlW/YHhWsg44u6oeqKr/zTD7fMou1Do1VbWpqi5rNW4D/hj4tyvdT5KnMgTnr1fVrW3fz6iqg5e4vWWRfTyWYUb+nqq6pzUfxPAsb9Q9DDP6nfUt5StV9eWqepBhdr+W4RnAz4ELgbkkB6/ksWuyDPq+/fPI8k8Yfogfpqo2Ab/NMDO7I8mFSZ64xP6eCNxVVfeOtH0POHKk/9aRvtHlRduSvD7JNUm2J9kO/EvgsJEht48sz/9yWNi2w+Mao9apaS+GXpjkB0l+xHBZ47Dltluwj8cxzLx/r6q+sot1PJrhF/xXq+p9I133AY9dMPyxDM+Odta3lIVfjzvbL+b5dVj8a6TdxKAXVfXpqvpV4MlAAefOdy0YuhU4NMno7O5JDNf3AW5juGQzb91ih5tfSPJk4H8CvwU8vqoOBr4DZBcfykpqnebHtr6v7f8ZVfVY4D/x8Me002MneRTwaeDKqvqzBX3XLXhH0OjtQyPj9ge+yPB4/8uCQ1wHPHNk7DEMz+r+qd3WJDl2ZPwz2zbaSxn0j3BJnpbkhS0YfsowA5ufjd3O8LT7UQDt8sH/Bd6X5IAkzwDOYHixEeCzwLuSHJLkSIYA35kDGUJvW6vljQwz+lUbo9aHPbYJewzDzHh7Ow//dUH/7cAxO9n+DxjOzdsWdlTV02vkHUELbm8CSLIv8HmGr+Xrq+qhBbv5FPDK9vcUBwJnA1+oqnur6sfAF4CzkxyY5PnAKQyXf7SXMui1P8OLqHcyXOo5HPjd1ve5dv/DJN9oy69heMF2K/CXwO9X1WWt72xgC8OLcX/DEDYPLHXgqvou8EfAPzCE379ieJfNpOys1sUe26S8B3g2w7XtSxmCc9T7gP/WLle9fYm6TwDuHn3v/gqO/zzgFQwvVm8f2ccLAKrqOuBNDIF/B8MvptHr+29heAfVHcBngDe3bbSXSpX/eETTkeTNwGlVteIXIiVNjjN6TUySI5I8P8mjkjwN+B2GmbSkGdqj/zpRe539GP7i8miG93VfyMP/4lLSDHjpRpI656UbSercHnHp5rDDDqu5ublZlyFJe5Wrr776zqpau9y4PSLo5+bm2Lhx46zLkKS9SpLvjTPOSzeS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktS5PeIvY7X3mNtw6UyOu/mck2dyXKkHzuglqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXNjB32SfZJ8M8mX2vrRSb6W5MYkFyXZr7Xv39Y3tf656ZQuSRrHSmb0bwOuH1k/Fzivqo4F7gbOaO1nAHdX1VOA89o4SdKMjBX0SY4CTgY+3NYDvBD4fBtyAXBqWz6lrdP6T2rjJUkzMO6M/k+AdwAPtfXHA9ur6sG2vgU4si0fCdwK0PrvaeMfJsmZSTYm2bht27ZdLF+StJxlgz7JK4A7qurq0eZFhtYYff+/oer8qlpfVevXrl07VrGSpJVbM8aY5wO/luTlwAHAYxlm+AcnWdNm7UcBW9v4LcA6YEuSNcDjgLsmXrkkaSzLzuir6l1VdVRVzQGnAVdU1WuBK4FXtWGnAxe35UvaOq3/iqraYUYvSdo9xpnRL+WdwIVJ3gt8E/hIa/8I8Ikkmxhm8qetrkQtNLfh0lmXIGkvsqKgr6qrgKva8s3AcxcZ81Pg1ROoTZI0Af5lrCR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1LnVfKiZtNvM6oPcNp9z8kyOK02SM3pJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOrds0Cc5IMnXk3wryXVJ3tPaj07ytSQ3JrkoyX6tff+2vqn1z033IUiSdmacGf0DwAur6pnA8cBLk5wAnAucV1XHAncDZ7TxZwB3V9VTgPPaOEnSjCwb9DW4r63u224FvBD4fGu/ADi1LZ/S1mn9JyXJxCqWJK3IWNfok+yT5BrgDuAy4CZge1U92IZsAY5sy0cCtwK0/nuAxy+yzzOTbEyycdu2bat7FJKkJY0V9FX1i6o6HjgKeC7wy4sNa/eLzd5rh4aq86tqfVWtX7t27bj1SpJWaEXvuqmq7cBVwAnAwUnWtK6jgK1teQuwDqD1Pw64axLFSpJWbs1yA5KsBX5eVduTPBp4EcMLrFcCrwIuBE4HLm6bXNLW/6H1X1FVO8zopb3B3IZLZ3bszeecPLNjqy/LBj1wBHBBkn0YngF8tqq+lOS7wIVJ3gt8E/hIG/8R4BNJNjHM5E+bQt2SpDEtG/RVdS3wrEXab2a4Xr+w/afAqydSnSRp1fzLWEnqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1btmgT7IuyZVJrk9yXZK3tfZDk1yW5MZ2f0hrT5L3J9mU5Nokz572g5AkLW2cGf2DwO9U1S8DJwBvTXIcsAG4vKqOBS5v6wAvA45ttzOBD068aknS2JYN+qq6raq+0ZbvBa4HjgROAS5owy4ATm3LpwAfr8FXgYOTHDHxyiVJY1nRNfokc8CzgK8BT6iq22D4ZQAc3oYdCdw6stmW1rZwX2cm2Zhk47Zt21ZeuSRpLGMHfZKDgL8AfruqfrSzoYu01Q4NVedX1fqqWr927dpxy5AkrdBYQZ9kX4aQ/1RVfaE13z5/Sabd39HatwDrRjY/Ctg6mXIlSSs1zrtuAnwEuL6q/nik6xLg9LZ8OnDxSPvr27tvTgDumb/EI0na/daMMeb5wOuAbye5prX9LnAO8NkkZwDfB17d+v4aeDmwCfgJ8MaJVixJWpFlg76q/g+LX3cHOGmR8QW8dZV1SZImxL+MlaTOjXPpRkuY23DprEuQpGU5o5ekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc5/PCLtoWb1j202n3PyTI6r6XFGL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknq3LJBn+TPk9yR5DsjbYcmuSzJje3+kNaeJO9PsinJtUmePc3iJUnLG2dG/zHgpQvaNgCXV9WxwOVtHeBlwLHtdibwwcmUKUnaVcsGfVX9HXDXguZTgAva8gXAqSPtH6/BV4GDkxwxqWIlSSu3q9fon1BVtwG0+8Nb+5HArSPjtrS2HSQ5M8nGJBu3bdu2i2VIkpYz6Rdjs0hbLTawqs6vqvVVtX7t2rUTLkOSNG9Xg/72+Usy7f6O1r4FWDcy7ihg666XJ0larV0N+kuA09vy6cDFI+2vb+++OQG4Z/4SjyRpNtYsNyDJZ4ATgcOSbAF+HzgH+GySM4DvA69uw/8aeDmwCfgJ8MYp1CxJWoFlg76qXrNE10mLjC3grastSpI0Of5lrCR1zqCXpM4te+lG0iPL3IZLZ3bszeecPLNj98wZvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6tyaWRcgSfPmNlw6k+NuPufkmRx3d3FGL0md2+tn9LOaAUjS3sIZvSR1zqCXpM4Z9JLUub3+Gr0krdYsX+vbHe/4cUYvSZ0z6CWpcwa9JHXOoJekzk0l6JO8NMkNSTYl2TCNY0iSxjPxoE+yD/AB4GXAccBrkhw36eNIksYzjRn9c4FNVXVzVf0MuBA4ZQrHkSSNYRrvoz8SuHVkfQvwKwsHJTkTOLOt3pfkhinUslqHAXfOuogVsN7p2ZtqBeudtonVm3NXtfmTxxk0jaDPIm21Q0PV+cD5Uzj+xCTZWFXrZ13HuKx3evamWsF6p21vq3cal262AOtG1o8Ctk7hOJKkMUwj6P8RODbJ0Un2A04DLpnCcSRJY5j4pZuqejDJbwFfBvYB/ryqrpv0cXaTPfrS0iKsd3r2plrBeqdtr6o3VTtcPpckdcS/jJWkzhn0ktS5R3zQJzk0yWVJbmz3hywx7n8l2Z7kSwvaP5bkliTXtNvxe3i9Ryf5Wtv+ovaC+axrPb2NuTHJ6SPtV7WP0pg/t4dPqc6dfmRHkv3budrUzt3cSN+7WvsNSV4yjfomVW+SuST3j5zPD+0h9f6bJN9I8mCSVy3oW/R7Yw+t9Rcj53bPegNKVT2ib8B/Bza05Q3AuUuMOwl4JfClBe0fA161F9X7WeC0tvwh4M2zrBU4FLi53R/Slg9pfVcB66d8PvcBbgKOAfYDvgUct2DMW4APteXTgIva8nFt/P7A0W0/++zB9c4B39ld36srqHcOeAbw8dGfpZ19b+xptba++3bnuV3J7RE/o2f4eIYL2vIFwKmLDaqqy4F7d1dRO7HL9SYJ8ELg88ttPyHj1PoS4LKququq7gYuA146xZoWGucjO0Yfx+eBk9q5PAW4sKoeqKpbgE1tf3tqvbOwbL1VtbmqrgUeWrDt7v7eWE2tezSDHp5QVbcBtPtduTzwB0muTXJekv0nW94OVlPv44HtVfVgW9/C8JEV0zJOrYt9ZMZoTR9tT4V/b0phtdzxHzamnbt7GM7lONtO2mrqBTg6yTeT/G2SF0y51ofV0qzkHO3u87va4x2QZGOSryaZ5gRqxR4R/zM2yd8A/2KRrndPYPfvAv6Z4ane+cA7gbNXs8Mp1jvWx1OsaIerr3VnNb22qn6Q5DHAXwCvY3jKPEnjnJOlxkz8fI5hNfXeBjypqn6Y5DnAF5M8vap+NOkix6hl2tvuitUe70lVtTXJMcAVSb5dVTdNqLZVeUQEfVW9aKm+JLcnOaKqbktyBHDHCvd9W1t8IMlHgbevotT5fU6r3juBg5OsaTO9VX88xQRq3QKcOLJ+FMO1earqB+3+3iSfZnhqPemgH+cjO+bHbEmyBngccNeY207aLtdbw4XkBwCq6uokNwFPBTbOuN6dbXvigm2vmkhVSx9vl7+eVbW13d+c5CrgWQzX/GfOSzfDxzPMv5p/OnDxSjZuATZ//ftU4DsTrW5Hu1xv+0G/Eph/t8CKH+8KjVPrl4EXJzmkvSvnxcCXk6xJchhAkn2BVzCdczvOR3aMPo5XAVe0c3kJcFp7l8vRwLHA16dQ40TqTbI2w/+LoM06j2V4gXPW9S5l0e+NKdUJq6i11bh/Wz4MeD7w3alVulKzfjV41jeGa5eXAze2+0Nb+3rgwyPjvgJsA+5n+M3/ktZ+BfBthhD6JHDQHl7vMQxhtAn4HLD/HlDrf271bALe2NoOBK4GrgWuA/4HU3pHC/By4J8YZl/vbm1nA7/Wlg9o52pTO3fHjGz77rbdDcDLdtP37C7VC/yHdi6/BXwDeOUeUu+/bt+jPwZ+CFy3s++NPbFW4HktB77V7s/YHed23JsfgSBJnfPSjSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9Jnft/H4VolnNSNH4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1, ax1 = plt.subplots(1, 1)\n",
    "ax1.plot('flux_prof', 'z', data = keep_dict)\n",
    "ax1.set(title='Ens 0: vertically averaged kinematic heat flux',\n",
    "        ylabel='z (m)', xlabel='flux (K m/s)')\n",
    "\n",
    "fig2, ax2 = plt.subplots(1, 1)\n",
    "z200 = np.searchsorted(keep_dict['z'], 200)\n",
    "flux_200 = keep_dict['flux_prime'][z200,:,:].flat\n",
    "ax2.hist(flux_200)\n",
    "ax2.set(title='histogram of kinematic heat flux (K m/s) at z=200 m')\n",
    "\n",
    "fig3, ax3 = plt.subplots(1, 1)\n",
    "wvel200=keep_dict['wvelprime'][z200,:,:].flat\n",
    "ax3.hist(wvel200)\n",
    "ax3.set(title=\"histogram of wvel' at 200 m\")\n",
    "\n",
    "fig4, ax4 = plt.subplots(1, 1)\n",
    "Tprime200=keep_dict['Tprime'][z200, :, :].flat\n",
    "ax4.hist(Tprime200)\n",
    "ax4.set(title=\"histogram ot T' at z=200 m\");\n",
    "\n"
   ]
  }
 ],
 "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.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}