ipynb-lamanno2016-bayesGLM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1\"><a href=\"#Imports\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Imports</a></div><div class=\"lev1\"><a href=\"#Define-the-model-in-Stan\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Define the model in Stan</a></div><div class=\"lev1\"><a href=\"#Run-the-script-to-do-MCMC-sampling-(using-pyStan)\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Run the script to do MCMC sampling (using pyStan)</a></div><div class=\"lev1\"><a href=\"#Load-the-Cef-file\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Load the Cef file</a></div><div class=\"lev1\"><a href=\"#Load-the-posteriors-and-define-Binary-Patterns\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Load the posteriors and define Binary Patterns</a></div><div class=\"lev1\"><a href=\"#Plot-violins\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Plot violins</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [],
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#imports\n",
    "import pandas as pd\n",
    "from Cef_tools import *\n",
    "from itertools import izip\n",
    "import pystan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define the model in Stan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "        # Bayesian Negative Binomial regression for single-cell RNA-seq\n",
    "\n",
    "        data {\n",
    "            int<lower=0> N;                 # number of outcomes\n",
    "            int<lower=0> K;                 # number of predictors\n",
    "            matrix<lower=0>[N,K] x;         # predictor matrix \n",
    "            int y[N];                       # outcomes\n",
    "        }\n",
    "\n",
    "        parameters {\n",
    "            vector<lower=1>[K] beta;  # coefficients\n",
    "            real<lower=0.001> r;  # overdispersion\n",
    "        }\n",
    "\n",
    "        model {\t\n",
    "            vector[N] mu;\n",
    "            vector[N] rv;\n",
    "\n",
    "            # priors\n",
    "            r ~ cauchy(0, 1);\n",
    "            beta ~ pareto(1, 1.5);\n",
    "\n",
    "            # vectorize the scale\n",
    "            for (n in 1:N) {\n",
    "                rv[n] <- square(r + 1) - 1;\n",
    "            }\n",
    "\n",
    "            # regression\n",
    "            mu <- x * (beta - 1) + 0.001;\n",
    "            y ~ neg_binomial(mu ./ rv, 1 / rv[1]);\n",
    "        }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the script to do MCMC sampling (using pyStan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "%run pystan_bayesmodel.py -i data/Mouse_Embryo_nound.cef -o output_folder/ -p 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the Cef file\n",
    "This is just reading the annotation, actual data is not used here. Only posteriors will be considered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "cef = CEF_obj()\n",
    "cef.readCEF('data/Mouse_Embryo_nound.cef')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [],
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "clustering_attribute = 'Cell_type'\n",
    "for i,v in izip(cef.col_attr_names, cef.col_attr_values):\n",
    "    if i == clustering_attribute:\n",
    "        predictor_list = v\n",
    "    if 'total' in i.lower():\n",
    "        total_molecules = [float(j) for j in v ]\n",
    "\n",
    "for i,v in izip(cef.row_attr_names, cef.row_attr_values):\n",
    "    if 'gene' in i.lower():\n",
    "        gene_names = v\n",
    "\n",
    "total = sum(total_molecules)/len(total_molecules)\n",
    "total_molecules_norm = [j/total for j in total_molecules]\n",
    "predictors = ['Size']\n",
    "for i in predictor_list:\n",
    "    if i not in predictors:\n",
    "        predictors.append(i)\n",
    "\n",
    "predictor_matrix = []\n",
    "for i, c_p in enumerate( predictor_list ):\n",
    "    predictor_matrix.append( [total_molecules_norm[i]] + [float(c_p==p) for p in predictors[1:]] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the posteriors and define Binary Patterns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [],
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#Load regression filenames and genes names\n",
    "gene_names = cef.row_attr_values[0]\n",
    "get_path_from_name = lambda x: 'output_folder/beta_%s.npy' % x\n",
    "file_paths = map(get_path_from_name, gene_names )\n",
    "\n",
    "# !!!! For a quick test uncomment the following two lines !!!!!\n",
    "file_paths = ['output_folder/beta_Th.npy','output_folder/beta_Ednrb.npy']\n",
    "gene_names = ['Th','Ednrb']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "order_predictors = ['Size','mEndo','mPeric', 'mMgl','mEpen', 'mRgl3', 'mRgl2', 'mRgl1',\n",
    " 'mNProg', 'mNbM', 'mNbML2', 'mNbL2', 'mNbDA', 'mDA0', 'mDA1', 'mDA2',\n",
    " 'mNbML1', 'mNbML3', 'mNbL1', 'mNbML4', 'mRN', 'mGaba1b', 'mGaba1a', 'mGaba2', 'mSert',\n",
    " 'mOMTN', 'mNbML5']\n",
    "reordering = [predictors.index(i) for i in order_predictors ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     1
    ],
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "x = array(predictor_matrix, dtype=float)\n",
    "total_relative_to_basal = ( (x>0) * x[:,(array(predictors) == 'Size')] ).sum(0) / (x>0).sum(0)\n",
    "reordered_total_relative_to_basal = total_relative_to_basal[reordering]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Define Binary patterns\n",
    "pattern_dict = {}\n",
    "\n",
    "for path, gene in zip(file_paths, gene_names):\n",
    "    d = load(path)[:,reordering] - 1 #because in the mode 1 is added\n",
    "    \n",
    "    # Condition 1 to define binary pattern\n",
    "    difference_basal = d[:,0:1]*reordered_total_relative_to_basal - d + 1e-12\n",
    "    #Check that the differece with the bassal is positive for almost all the realizations\n",
    "    bigger_condition = (difference_basal < 0).sum(0) >= 998 # For max strictness do == 1000\n",
    "    #Check if the gene is \"basal\"\n",
    "    bigger_condition[0] = (difference_basal > 0).mean() > 0.999\n",
    "    \n",
    "    #Condition 2 to define binary pattern\n",
    "    median_categories = median(d/reordered_total_relative_to_basal, axis=0)\n",
    "    median_condition = median_categories >= 0.35 *max(median_categories)\n",
    "    \n",
    "    #Condition 3 to define binary pattern\n",
    "    high_enough = (max(median_categories) > 0.4 )*ones_like(median_condition)\n",
    "    \n",
    "    #Logic and: all the condition must be satisfied\n",
    "    binary_pattern = bigger_condition & median_condition & high_enough\n",
    "    \n",
    "    #Store in the dictionary\n",
    "    pattern_dict[gene] = binary_pattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "code_folding": [],
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "median_dict = {}\n",
    "\n",
    "for path, gene in zip(file_paths, gene_names):\n",
    "    d = load(path)[:,reordering] - 1 #because in the mode 1 is added\n",
    "\n",
    "    medians_expression = d[:,0:1]*reordered_total_relative_to_basal + d + 1e-12\n",
    "    median_dict[gene] = median(medians_expression, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Make a  pattern dataframe\n",
    "pattern_df = pd.DataFrame( pattern_dict, index=order_predictors, dtype=int).T\n",
    "median_df = pd.DataFrame( median_dict, index=order_predictors, dtype=float).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot violins\n",
    "of the posterior probability distribution for the NB mean parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "gene = 'Ednrb' # the name of the gene to plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import violinplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "hide_input": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxkAAADxCAYAAACnHHqNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd8XFed9/HPmaJqWe6O417j2I5J7OMktCSUhIQshBIe\nEAEWyEIguxTTF1gILMsuGzZ+djcQOnlCgsjCUgJpJr07vnYS25Fb3HuRLFl92nn+uCNbZSTNjGZG\nlvR9v17zkubce8q1Z0b3N6cZ5xwiIiIiIiK5EhjsBoiIiIiIyPCiIENERERERHJKQYaIiIiIiOSU\nggwREREREckpBRkiIiIiIpJTCjJERERERCSnFGSIiIiIiEhOKcgQEREREZGcUpAhIiIiIiI5pSBD\nRGSIMcYkMnzckkHZf5vM84t8XkOyroQxJp7vekREpPBCg90AERHJigP2AS1pnHskz20RERHpQkGG\niMjQ9SHn3BOD3QgREZHuNFxKRERERERySkGGiMgIZIy5wRjjGWMajTG1xpjVxphLezn30uT8iQ3J\n539njFmbzHvCGPOwMeayFPkeS+Z7lzHmAmPMQ8aYJmPMh1Kce6Mx5gVjTEuyPfcZY96U8wsXEZGC\nUJAhIjLCGGPuAm4DLgCiQD3wBuBh4B19ZA0aY34N/ASYhz/XozSZd7Ux5qJu57vk443A08DrgUNA\nU7f2/A64FVgAHARKgCuBvxpjVmZ9oSIiMmgUZIiIjCDGmBuAKqAZeD8w3jk3F5iJH2Rcgx8YpHIu\n8G7gxmS+ecAsYCsQBD6Xqsrk+Q8AU5xz851zv+92zjXAP3UqcwLwg+SxfzfG2CwuVUREBpFxrre/\nJSIiciYyxiTwA4H99L+61PPOub9N5gsAe4CzgQ87537VrdwyYAswFfh/zrmPJtMvBR5N1vkV59zN\n3fL9LfBLYL9zbkan9EeBS4HDwDznXEu3fB3X8T3n3FdTXOcLwFLgLudcjyFWIiJy5tLqUiIiQ9e0\nNM452On3i/ADiFrgru4nOudakkOpvtxHeT9NkfZy8ueEFMcccHf3AKObH/aR/mPgkj7yiojIGUhB\nhojI0HWZc+7JDM4/P/nzJedcopdzNvSRv8k5V58ivTX5s7iXfDv7KPOkc25/L8c2JX9O7SO/iIic\ngTQnQ0Rk6DIZnj82+bOxj3P66nFoz7C+Dn3t6t3ax7G65M+AMUZfiomIDCEKMkRERo6OAOKsPs4Z\n28exfOirvinJn03OuVghGiMiIrmhIENEZOTYlvx5jjGmqJdzlheqMUlFxphzezl2anhXoRojIiK5\noSBDRGTkeAJ/qFQl8MHuB40xY/GXty30soPXp2iLAW5ItuV3BW6PiIgMkIIMEZERwjnXBNyCP5fj\nZmPM33QcM8bMAP4MjB6Epn3aGPNpY0w42ZZR+Bv+nQMcAH4+CG0SEZEB0EQ6EZGh61fGmP72yQB4\nyDn3qeTv3wFehb8B3j3GmBNAAzADaAO+DvxbPhrbi0PA88D/Bb5rjDmCv49HMf6Gge91zjUXsD0i\nIpIDCjJERIaejuFM6eyTAVBzKqNzcWPMu4FPAh8HFgBR4LfATfiTwh09h0ylSuutXZmIA+8BPg98\nBH/n8RPAX4FvO+d2ZFGmiIgMMu34LSIiIiIiOaU5GSIiIiIiklMKMkREREREJKcUZIiIiIiISE4p\nyBARERERkZxSkCEiIiIiIjk14pewNcaMB94C7MZfI15ERERERHoqAWYBDzrnavs6ccQHGfgBxl2D\n3QgRERERkSHiOuDXfZ2gIMPvweDOO+/k3HPPzSjjypUrWbVqVT7aNCj1FLKubOvZvn17xnluvvlm\nvvjFL2acry87dvTcH8w5uPNXd/HBD13X49jcuXNzWn8+rimXdc2fPz/jeobj67yQdemahkZdw/Ga\nClmXrkl1DVY9hazrTL6mzZs384EPfACS9899GRZBhjEmDOwAFjrnWowxY4HbgfOBVuDLzrk/9ZK9\nDeDcc89l2bJlGdVbWVmZcZ5sFKqeQtaVbT2hUOYv2VGjRmUcQPYnHo/3SLv/thNMOX4js2aN7XEs\n1/Xn45pyWdfSpUszrmc4vs4LWZeuaWjUNRyvqZB16ZpU12DVU8i6hsg19TvFYMgHGcaYrwLvBKZ2\nSv4BsM45d40x5nzgYWPMJOdczztDkSwsWbKkR9o929qZFIQlS4oHoUUiIiIiZ44hH2QAHrAJ+AOA\nMaYUuAL4MIBz7kVjzJsGrXUiIiIiIiPMkF/C1jm32jl3T6ek+cBx4N+MMRuNMc8BU9SLIYXi3GC3\nQERERGRwDYeejO4qgAXAV5xznzPGLMMfLrXQOXekt0wrV66ksrKyS1pVVRVVVVW9VtTXsVwqVD2F\nrKuQ13TVVVcVrC4AHGDyW0Uhr6lQdQ3H13kh69I1DY26huM1FbIuXZPqGqx6ClnXmXJN1dXVVFdX\nd0lraGhIu2zjhsnXrsaYOH6AMQ141jk3vtOx54HPO+eeTJFvGbBu3bp1BZtkI9nbsGHDYDehV9+9\nvB2ArzxQTCA4yI0ZZNlM/BYREZEz2/r161m+fDnAcufc+r7OHfLDpToxAM65bcAhY8xVAMaYhfiT\nwjcOYttkBHGJPHdjiIiIiJzhhtNwqc5dMtcAPzLG/CdwErjOOVc/OM2SkSYRh2B4sFshIiIiMniG\nTZDhnAt2+n0HcPkgNkdGsEQ0ACVaZ0BERERGruE0XErkjJBQfCEiIiIjnIIMkRxLxPS2EhERkZFN\nd0MiORaPaOK3iIiIjGwKMkRyLBbR20pERERGNt0NieRYtEVvKxERERnZdDckkgOJ2Onf204Om0Xb\nRERERLKiIEMkB1rqwp1+V5AhIiIiI5uCDJEcaDxclPJ3ERERkZFoWAQZxpiwMWavMaasW/ocY8wJ\nY8yiwWqbjAwNB4rxN5131O8rHuzmiIiIiAyqIR9kGGO+CjwDTO2WHgJ+yTC4RjnzndhTAgEHGOp2\nlQ52c0REREQG1XC4AfeAf06R/m/A3UBtYZsjI1HtzhJI+G+n5mNhom3aK0NERERGriEfZDjnVjvn\n7gFO3dUZY64G5jrnftg5XSQfnIO6nZ17L9SbISIiIiPbkA8yujPGTAW+C1w/2G2RkaHpWJhoa/B0\ngnHU7lCQISIiIiPXcFpr0yV/vg4YAzxljDH4czXuMcZ83jn3p94yr1y5ksrKyi5pVVVVVFVV5au9\nMkwc39ZlvQFMAI5vV5AhIiIiQ1d1dTXV1dVd0hoaGtLOP5yCDAPgnLsbfy6Gn2jMLuBtzrnNfWVe\ntWoVy5Yty28LZVg6UlOGCbpTz13ccHhT+SC2SERERGRgUn3Zvn79epYvX55W/kENMqy104CxQAKo\nBw55npfIsjjXR7rmZUjeHHihAhfv+hJrPFxM8/EQ5RNiveQSERERGb4KHmRYa18PvB1YDlR0O9xk\nrV0H/NnzvMczKdc5F+wlfU5WDRVJQ9PRMPV7S1IccexbO5qFV9UVvE0iIiIig61gQYa1dgnwOfzJ\n5o8Avwb2ASeTp4wCpgGLgQ9Yaz8G/Jfnec8Xqo0imdr5xBgwDly3zjIDrzwyVkGGiIiIjEiF7Mmo\nAr7ned7WXo7XJR8bgGpr7XzgOkBBhpyRXAI23zs+9UA9ZzhSU0793mLGzGgveNtEREREBlPBggzP\n876W4fnbgZvy0xqRgdvz7Giajhb1etwEHBt+P5FLPru/gK0SERERGXyFHC71rXTP9Tzvm/lsi8hA\nJeLg3XFW6qFSSS5heOXhsbzqPUepnBopcAtFREREBk8hN+NLZPAQOaNtWz2Ohv0lvQYYna352dkF\naJGIiIjImaOQw6XS7skQOZO11gd5/hdTSGd1ZJcw7Ht+NHvXVDDjosaCtE9ERERksBVyuNQn0jzV\neZ7347w2RmQAnvvJ2cTaAqS9/YpxPPXf07j2vK0UlamjTkRERIa/Qq4udUGa5/W2qZ7IoDvwwih2\nPj42s0zO0FofYv1dk7n4Y4fy0zARERGRM0ghh0vdkK+yjTFhYAew0DnXYoxZCPwCf9+Nk8A3nXP/\nm6/6ZWRIxOHZ287GBBwukeEm8s5Qc88Ezn1rbUaTwDdt2pRRNUuWLMmsXSIiIiJ5UPAdvztYa18D\nvA2YCnwJuBxY63nelkzKMcZ8FXhnspwOPwWqnXP/bYyZBTxnjHnSOXc0F22XkWnPM5U0HEi1u3f6\nXrp7Epd8bmBL2u545RUq27Zx9jmvJhLOsFdFREREpAAGJciw1l4DfB74A3AZUAyMBn5mrf2653mP\nZVCcB2xKltWhAr8nA+fcbmPMfuAcQEGGZG3z/eOy68VIcgnDjsfHcPEnDqY9NyNlz4RzzD2+jsXT\nDrFz4uuzaouIiIhIPhVyCdvOPgTc4nneKpJL1nqe9wPgh0C6E8QBcM6tds7dQ6dZuM65851zzQDG\nmPcBZwMv5qjtMgJF2wyHN47KOsDokIgFOPjSqJy0Kei0k7iIiIicmQYryDgLeDlFuoc/j2LAjDGz\njDG/Bb4JvNU5p/VDJWvHtpYNOMAAMEHH0c1lOWiRiIiIyJlrsOZkbAEssD35vGNFqdcB+7Is89Sq\nVMaYNwN3At8H3u+ci/aXeeXKlVRWVnZJq6qqoqqqKsvmyHByfHvZgIZKdXBxOLq5fEBlBIglC9NC\nbCIiIpIf1dXVVFdXd0lraGhIO/9gBRn/AfyntfY8IAh8xFp7FrAY+FyWZRoAY4wBbgeuc849nG7m\nVatWsWzZsiyrluHu0IbyHN3TG45tKyXWbggVZ1dgKN6ULEl7boiIiEh+pPqyff369Sxfvjyt/IMy\nXMrzvBrg3cArwGP4E7VfBt7jeZ6XZbEdd2yz8Ydj3WaM2Z58bDPGXDzAZssIFWlJzqNwAx8uBf68\njAPrK7LOXxKtBSCQSH8pXBEREZFCGrQlbPFXlHrR87xfAFhrL8Lf0yIrzrlg8tedDO51yTCz8/Ex\nJGK5CTAATMCxdfVYZr46u5d7WeQgACHXRlHsBJGQlrEVERGRM8ug9GRYay/FX3L2nZ2Svwj8wVq7\ncDDaJNKbzfeO77R22cC5hGH/2tE012YRC7sEpRF/JWYHjGrPdgqTiIiISP4M1upSNwK/8Dzvnzql\nvQd4GPjC4DRJpKf6vcXU7SrN2VCpDg7Y+diYjPOVRI8TIJ58ZihrP5DTdomIiIjkQsZBhrU2aK19\nv7U2+0Hl/jK1j3dO8DzPAX/C3zRP5Iywz6vAmDys4uRg79rRGWerbNuOS3arGBwV7bs1N0NERETO\nOBkHGZ7nxYH3AgMZ1rQLeFOK9IuB4wMoVySnaneW5nSo1GmG2h2lGeUoitYxvulFzOnVmgm4KJNO\nPpvrxomIiIgMSLYTpH8AfMFa+wP8idZdvur1PK+/MRw34y9h+zrgJSAGLADOB76aZZtEcq6lNpxy\nb4wjkbU90mqjNb2WM7loRY+0aEuQeNQQDPffUxKMtzCr7k89lq01wISWl2gPj+dE+ZJ+y+nLpk2b\neqQVNR6jYv/L1C68FEzXf4clSwZWn4iIiAxf2QYZ/5L8+R+d0hz+PY8DLuyewVpb7nleM4DneS9Z\na98BXAvMAUqBbcC/ep63P8s2ieRctC2/05bikf6DDOPizKy7h3C8qUsvRmdnNzxCJDSa5uIZOW3f\nsUf/QHPdfkrmrCBePLBNBEVERGTkyDbIeHsWef5ird0GPA087XneDuBnWdYvUhDxSOqxUql6JtI5\n1rP8AJT3vaneuOaXKI0e7XfU1rQTq9k2+cM4k93bOlXPRP1TpZSaMsoXnkOsPPOJ6iIiIjIyZXU3\n4nneIQBr7QJgOtAG7OxI78UVwDLgtcD3rLUlwDPJxxrP81qzaYtIPkWag/2fNADtzUFKx8Z6P8G5\n5DyMvhkgnGihom03J0vn5bKJAAQS0ZyXKSIiIsNXVkGGtXY08G+ABeqAEqDUWvsY8B3P8xq75/E8\nLwqsST5usdZOxQ843gncZK3djN/L8YzneTszaY8xJgzsABY651qMMbOAO4GZwBHgeufcS1lcqoxg\niTi01IXzWkfzsTBjprX3erwscpCiRFNaZTkMY1pqchxk+MOzApG2HJYpIiIiw122A84/i3/38VbP\n8670PO8y4APAFNLc58LzvAOe5/2P53mfwe/l+BVwFv6k8LQZY76K3xsytVPyL4A7nXPTgZuA6kzK\nFAE4eagYF8/L0lI+4zixu6TXw8FEK1PrHz61ZG2/xeEY3b6bypYtuWohJuEHGeGW+pyVKSIiIsNf\ntkHGa4Efep53arlZz/O2A/8JvD6dAqy1ldbayuTT6cBc4FHP896dYVs84J87nhhjxuMPy/oZgHPu\nL8BoY8yCDMuVEa6vACAXTMBRtyd1HcXROuYcu5uieEOvk71TccC0+tVMaPTADXx/j2CkBYCSEwcH\nXJaIiIiMHNkGGWH8IVLdJZLH+mStXQ7cC7zWWjsRuB24Gn9Z23dl0hDn3Grn3D2dkuYAB5xznQe6\n76NrT4dIvxoPF2ECediIL8nFA5w8UNwjfVTbLuYe+w1F8caMAgzw52YY4KzGZ5h+4n6M62O+Rz+C\n7c2EWxsAKDu2O+tyREREZOTJNsh4FH+fjFO7c1trpwOfAdLZGewzwB3A/cDlwFbP894H/CtwXZZt\n6uCA7rNU4ynSRPoUbQ3kZ7fvTiItXd+Co9p2M7PuLxhiGQcY3Y1u28GMunvB9b16VW9KOwUW5Udf\nyUnPiIiIiIwM2S5h+33ga8CvrLWt+JvpVQDrgO+lkX8W8C3P85y19kLgqWT6FmBylm3qsBd/+FVn\nM/A3DezVypUrqays7JJWVVVFVVXVAJsjQ1VxRZxEio34csdRUnm6pyEUb2b6iQcg7VkYfTM4RrXv\nYULTeo5X2IzzVxzczHETwLgEobYmiusP0T727By0TERERM501dXVVFd3ndbc0NCQdv5sl7BtBb5u\nrb0Vf3hSMbA3ufdFOo4BM621DcBy4EfJ9EXA8V5z9c0AOOeOGmNeNMZ8wDl3pzHmPcAe51yfg8pX\nrVrFsmXLsqxahqOJC1rA5S/IMAGYdI4/5wHnOLv+EQIumpMA41QdwOTG52gsmU17eHza+YKtJxm3\n9Sl2J3tBnAkwceNq9l/y4Ry2TkRERM5Uqb5sX79+PcuXL08rf9pBhrV2BfCS53mR5O8dosnHOGvt\nOADP89b2U9zPge/gz+HY4HneFmvth4CPAT9Ot03ddB7L8XHgV8aYbwOHgA9mWaaMYJMWtlA2PkpL\nbQhyeuvvcwnDnEv9bwRGt+1gdPuuHues3XqsR1rN3hO9lrninImpamJq/UPsnPB/wPR/HSYeZeaj\nPyMQPz3C0LgE43asoXnKAk7Mf02/ZYiIiMjIlklPxg+BdwAHkr/3xgEX9lWQ53n3WWu34k/GXpNM\nPgx8yfO8dOZ09KzUuWCn33cAuhOSATEBuKDqCE/fmvs1A0zAMXV5I+Nm+ftPTGjycOQjlPGHTZVF\nj1AWOURLcd/DnUw8ysxHfkL5kVdSzgmZ9tSduGCY+jnp72ouIiIiI0/aQYbneStS/T4Au/HnX1xj\nrb0P2J5MEzljLLiijs33jufEnhJcDudnmKDj1R8/PYKvONaQMsBI3TPR/7FUiuINtNB7kGFiUWY+\n8iMqDmzuY9K5Y/rjv4REgvp5F2VUv4iIiIwcGa8uZa0NWmv/aq2dlW2l1trZwD3At4DPA+OBfwTu\nsNZOyrZckVwLBOENX9pLIOhggKs9dXbxxw8y+uzIqedt4Qk5mu7du7a+5mQk4sx85Mf9BBgdPS2O\n6U/ezpid/Y2KFBERkZEq44nfnufFrbXPAG+l72FTffkC8Dz+vIwnkmkrgW8DXyLNXcNFOtu0aVOP\ntIaA4cmiIDYS56xE15vnJUuWpFXumBntvP6z+3ns5hkDbqMxjtmX1LPwqrou6XVliymPHBhw+ak4\noD00jrZQ7z0fZz93NxUHatJaNtcky5z+xO1EysfRMnluztoqIiIiw0O2S9hGgeusta/BXxq2y52J\n53nf7Cf/q4CPJwOWjjzN1tpfAD/Isk0iPTy9eydPF4UomzmLs9rjWZcz97J6aneWsPF/J5LtzAkT\ncIyd1cbrP7O/x/zrk6XziZ58ilCiJef9GQY4Nsr2Oum7Yt9GJmx9MuMynXPMeOxnbH3XTbhwz00F\nRUREZOTKNsgIAqsHUG89/r4a3Y0BIinSRfqVqmdiTVGA0pIwExcuYknLwPZjtH97mBN7SziwriLj\n+Rkm4CiuiHHFTbsIFffsLXAmyOHKS5L7ZOSOw9AWnkBD6YLU7Yq2M+2Zu3CYLHYXd4RbGpj84n0c\nXvHOXDRXREREhols98n41gDrrQa+aK29Jfl8prV2MXAj8PsBli1ySjzZLxDr57x0dMzP+NNn5tN4\nuCizQMPAFTftpnx87y1pKF1ASeQIE5tfyEFr/QAjFihj79ir/aWyUpiy7o+EWk5mvbu4wTFx019p\nmL2M1gkzB9JcERERGUay7cnAWnst8DfANODDwHuA9Z7nPd5fXs/z7rLWngS+CBQB/wHUAncCd2Xb\nJpHuXLefA1VUluDyb+zmj5+aTzyDIOM1nzzAxAWt/Z53ZPTriAUrOOukP3wp25t/B7SGJ7N33NXE\nguUpzxl1oIYJmx/LqvwujGHGYz9n2zu+jgsVDbw8ERERGfIyXl0KwFr7YfwN7x4AyvCHTx0EvmWt\nTWvchOd5f/Y8753A64E3eJ53ped5d3qel7slfGTE63gxJXJY5pjp7Vz08T43kD/FBBzTV5zknCvr\n+j8ZwBhqR53PjgnvJRKsyHjFKZcMS46NupCdE67tNcAItjUx44lf5mRFK+MSFDUe5+znfzvgskRE\nRGR4yLYn41rg3z3Pe8ha+ykAz/Putta2AB8B/tA9g7X27X0V2DEBPFnWPVm2qwtjzOfxg6Eg8BLw\nEefcyVyULUNDPHkPHcvxbOqFV9ax/a/jOP5KaZ/DpgJBx2v/4UA6G2130VY0iR0T38+0Ew9S0b4r\nrVDAYUiYMHvHXU1z8fQ+TnRMf/L/EWxvybqnpDuDY/zWp2icuoSTM1+VkzJFRERk6Mo2yBhD6o3z\ntuFvsJfK36VZtsPfQ2NAjDFLgM8AC51zLcaYXwBfBr420LJl6GhP3t23ZnqX3w8TgBXXH+K+L/e+\nfKsJOBZfc5zyCdlNOE8Eitg77mqmnVhNZdu2PgONjgBj54RraQ9P6LPc8ZsfZ/T+nsv9DpQDpj95\nO1sn9Le4nIiIiAx32QYZLwBXAK8kn3d8Hfo2YEuqDJ7n9dmTkQcJ/JU2K4wxEfxhXfUFboMMspbk\nnXlDHva5m7KkmfFzW6jdWQoudQWLrzk+sEpMgP1jLydc20RZ5FCfPQ97xr2t3wAj2N7MFC/12gpr\nDzf2SKupbem1rBVndV0gzgCBWMQv/9WX9NkOERERGd6yDTK+B/yXtfbCZBlfttaeDZQAn+4vs7V2\nWV/HPc9bn2W7TnHO1RhjHgb2AK1AI/D3Ay1XhpbGZA/GiUB+dtNeeFUdT986tUe6CTim2UbKxuVg\nXSsTZO/YtzL/6J0EXVvKUw6Pfi0txT3b0V354e0E4gNbyrcvxiUYvW9j3soXERGRoSHbJWz3W2vf\ng9+bMT9ZzsPAA57n9f7V52k/7iU9DjQBb86mXZ0ZY64DFuAP32rEX8Hqp8C7Blq2DB0NyeCiMeD3\nAeQ61Jj5moaUQYZLGGa/PncdZ/FgGQfGvpmZdX/pWg+GlvBZ1JZfkF45RWW9HuveM5HusR51hEsJ\npn22iIiIDEdZBRnW2o8Bv/E8737g/k7pFdbaD3qe96t+iri42/MiYDbwMeDRbNqUwquBPzjnGgCM\nMT8F/tLbyStXrqSysrJLWlVVFVVVVTlqjgyGpuT6aVFjaDVQluO1y0or44yf20rtjtJuRxzTl/cc\nfjQQjSVzqC+Zj+P4qWFTDsOBsZf3upt3d82T5xItqyTU0pDzgAvAGcOJ+a/udWKWiIiIDA3V1dVU\nV1d3SWtoaEg7f0ZBhrX26uSvHwNarbUnup0yB3gv0GeQ4XlevFtSK1Bjrf06/j4ZuVhdygM+Yoz5\ngXOuBbgGeLq3k1etWsWyZX2O4pIhph2IdLqVPhEwlMVzv0LyNNtI3a6uQcbYmW2UVHZ/mQ/c4cpL\ncDyHweGA46MskdCY9AsIBDlk38WMJ36Z87Y5DPGiMo4tfpOCDBERkSEu1Zft69evZ/ny5Wnlz7Qn\n48bkTwN8AH94U2dR4I4My+xsFDB+APlPcc7dboyZA7xojHHAJvzlbGWEOBbs+l39sYBhah6CjLOX\nNvPS3afrMkHH1GVNOSl706aeq0DVHi5jVGQfjgC7JxSRMKfPWbJkSb9l1s9ZwagDNYzdsSZnvRkO\nwBj2XvZ3JIp7H5IlIiIiI0NGQYbneVcDWGvvAT7ued7hbCq11n4rRXIJsJQ+ehsy5Zz7BvCNXJUn\nQ8v+4Om9JgPOsT8Y4PxoLrfl801a2IwJnA5eXNxw1uLmnNfT4eyFFzGtvomWoskkTBY7bBvD/td9\nEBcMMX7bwN9uzgRwwRB73vAxms5eOODyREREZOjLduL32621ldbaSs/zGqy1c4HXADWe561Lo4hU\nd3oNwO3AH7Npk0h3u0IBAjgSyRkMO0NZbXDfr1CJY9ysVg53Wrx58qLcBBm99UyEYxdQFihjfCCc\nXcGBIAdecx2t46cz9bn/Aeey2pjPGUNk1Hh2v/mTtI+Zkl1bREREZNjJduL3cuA/ge9aa9fiBwcH\ngBustbd4Xi8L8Sd5nvetZDnlnuc1J38vS3NlKpG0bA0FSCQHBDlj2BoK5GWFKYDJi1uo2ZaARIDy\niZG8zMfoLBqq7P+k/hhD3cJLaa+cwqyHfkAgFs0o0HDG0DxpLrvf9EkNkRIREZEusv1q97P4k7vv\nBy4Htnqe9z7gX4Hr+stsrZ1trf0N8P1OyfdYa2+x1o7Nsk0ip7QBe7vNyWgJGA7nab+MCfNbIBEA\nHJMWDq1YuXnKAnZe+VlcMIRLMwTrCDB2XfEpBRgiIiLSQ7ZBxkzgIc/zHHAh8FQyfQuktbDMV/A3\nyftqp7RlGC9RAAAgAElEQVTrgXLgi1m2SeSUbeEArvuyrs6xOZyfIVPjZnVskmcYN7s1L3XkU+vE\n2ey6/B9wgUC/fRnOBGgdN53dl/89LpTFnBAREREZ9rLd8fsYMNNa2wAsB36UTF8EHE8j/yLgo57n\nnVoC1/O8Pdba/wL+O8s2iZzyUjhIwHW9XTbAC+Egb2wf+FCm7qs+xaOG2mgEgGOxg2za1HVORjqr\nPg225ikL2POGjzPr4R9BL6GGMwHaKyaw6y2fJhEuKWwDRUREZMjI9mvdnwPfwZ+kvdHzvC3W2g8B\nXwJ+l0b+E8CsFOmTgViWbRIB/BfQmqIgiW49Gc4YasIBTuZhxFQw7Jhx1jzGhxcxemok9xUUSOOM\npRyy70x5zAGJYJjdV/wD8eLywjZMREREhpRsV5e6z1q7FZgKPJ9MPgx8yfO8Z9Mo4pfAN6y15wIb\n8O8LFwDvB36TTZtEOjxTFKS5l7kXDnioJMS7WgcWy6bqmVh0F7iEIRg+Z0BlD7bjS97EmF3rKK3b\n1yXdAIcufDeRiomD0zAREREZMrIeoO553g78wOK11trX4vdopBNg4HneH/B7PRYC/wh8F7gMWOV5\n3s+ybZNIO/D7sjC43ob7GB4oCXEiD70ZgaDfozHkmQCHLnw3xp1eadoB7RUTqJv/2sFrl4iIiAwZ\n2S5hOxr4N8ACdfgb6ZVaax8DvuN5XmN/ZSQDkmetteVAS3ISeU4ZY64AVgGlwA7gg865rDYQlKHh\nz6UhfzhU90nfncSB6rIwNzZHC9auoab5rPm0jJ+BO1xzalnbY+e9BQL5mTgvIiIiw8tAlrB1wFs9\nz7vS87zLgA8AU4Av9JfZWhu21n7KWrsaeASYYa39F2ttVZbt6cEYMx5/WNZ7nXNzgI3A13NVvpx5\nDgYM95eEeq4q1U3CGNYWh3g5T5vzDRd1C17DqQngJkD97OWD2h4REREZOrK9y3ot8EPP806tJOV5\n3nb8Dfpen0b+TwJvxJ883rHUzyPAh6y1H82yTd29E7jXOdexDNBNwH/lqGw5wzjgV+Xp735tnOOO\n8rBWGehD49Qlp3bNaJ40h0RR6aC2R0RERIaObIOMMP4Qqe4SyWP9uRq42fO8J0h+Vep53sPAzcC7\nsmxTd4sAZ4xZbYzZAvwYaMhR2XKGWRcOsCXcc0Wp3jhjOBowPFSS7SrOw1901DhiRf5Ge81nzR/k\n1oiIiMhQkm2Q8SjwBWvtqWV0rLXTgc8A6Uz+LiX1fhqHgDFZtqm7CvxelQ865xbib/73kxyVLWeQ\nNuCu8iJML5O9e2UMfyjNzyTwYcEYYqWjAWgbO3WQGyMiIiJDSbZf434f+BrwK2ttK/4StBXAOuB7\naeR/BnivtfY7yefOWlsMfBBYn2WbujsG3O+cO5J8fgfw595OXrlyJZWVlV3SqqqqqKrK2TQRyZO7\ny8I0GPqdi5FKDPh5eRGfa4pkv9TaMOaSE70jo8YPcktERESkkKqrq6muru6S1tCQ/qCgjIIMa+1S\n/N6BMP4N+63AHKAY2Jtc1jYd38Nfneq+ZFnfByYBB4DPZ9KmPtwD/NQYM9Y5dwJ4N/BUbyevWrWK\nZcuW5ahqKZQ1RUEeG8CQp4QxvFwU5P6SEFe3jewZGt13MQdoPHqSkpMtlO7aT+xIU5djQ2EXcxER\nEclOqi/b169fz/Ll6S0Ek/bdmbX2rfiTpw8ArUAV8EvP836UbhkdPM87AdxgrT0fmJ9sx25gjed5\nib7ypss595wx5nvAU8aYIPAS/oRzGSb2BQ0/K0/uiZFFL0ZnvysNMSOe4LxoTl5+w8ac+QsYs7uJ\nHeHiwW6KiIiIDCGZfAX8Ufyg4jY4FXT8k7X2p57nxfvOCtbaVKNRNiQfXc7LYaBxJ3BnLsqSM0sr\n8F+jikjAgAMM8Hezvq28iO80tDFuGOynl41UPROBc+ZRUneARZPnDkKLREREZKjKJMiYDjzQ6fmD\nwDfx98bYn0b+NZxadL9fF2bQLhmBfl8WpjZgspqHkYozhnYcd5QX8dmmSE7KHA4S4RJaFGCIiIhI\nhjIJMgLAqbsvz/Pi1toYEEwz/ycyaZhIb+oNPFIczFmA0SFhDC8VBdkVNMyOj9DuDBEREZEcKNgm\nAZ7nreueZq2twO8haQMOeJ7XXqj2yND1fFEw7S6xTAWc47niELNbonmqQURERGT4yzTIeJ+19mSn\n50Hg/1hr6zuf5HneT/sqJLlc7ZeAt8GpTYVbrbX/i7+TuO7wpFfbwul2nmUuYQw1IS1mKyIiIjIQ\nmQQZh4FLu6UdxV/StjMH9BlkADcCC/Enk2/DXwL3AuBzQBH+zt8iKe0L5m4uRiqHgwbH6ehXRERE\nRDKTdpDhed7bcljvm4Cve57XsTB/BHgiOcfjJhRkSC8ccCKQ+va/+cUNPdJat73Sa1nl5y9NmR4z\nhkYDozUtQ0RERCQrgzUuZDSQapnak0BZgdsiQ0gEiOaxF6NDYy+BjIiIiIj0r2ATv7t5Fviktfar\nyY35sNaWA9fjb5onklJ7H/f+vfVM9HcsZT0ZnS0iIiIinQ1WkPHvwPeA+6y1e4AYMBM4gj8vI6eM\nMZ8H3uecW5HrskVEREREpKtBCTI8z6sF/s5aewEwF3+y917g2XR2D8+EMWYFsBI4mMtyZXAUFWie\nRLgw1YiIiIgMSwULMqy1U1MkH00+OpxlrcXzvAO5qNMYUwn8EPgy8JlclCmDqxgwzuV1dSmAMk36\nFhEREclaIXsy/ghd9lDruEvsnuaAC3NU50+AfwFO5Kg8GWQGKHfQlOd52aMTijJEREREslXIIKMO\nqAQ2Ao/jT/5uyVdlxpgbgBPOuT8aYy7LVz1SeGfFE7wSyN+GfJUJp+FSIiIiIgNQyCDjSmAJcAn+\nTt+fBNYDTwBPeJ53JMf1XQUsMsbUAOXARGNMjXNuUaqTV65cSWVlZZe0qqoqqqqqctwsGahzYgl2\nhgIk8jBkKuAc50RzOi1IREREZMiprq6murq6S1pDQ0Pa+Y1zgzMsJDlH41L8oON8YAfwJH7AUZPL\nuowxlwI3O+d6DMMyxiwD1q1bt45ly5blslrJgw0bNnAgaPh6ZUla53ds0JfJErafbWznVdFU27hI\nupYuzWzJYBERETnzrV+/nuXLlwMsd86t7+vcwVrCtmNy96+BX1trK4DrgA8BHyV3czJkGJoad6xo\nj7GuKJjT3oyAc8yMO85TgCEiIiIyIIMWZABYa1+F35NxKTANeAF/vkZOOeceR4HLsPKBlijbwkEa\ncTkJNIxzhICPNUUIDLx5IiIiIiNaQYMMa20xcBFwGfA6/O0IngF+CjzteV5TIdsjQ9doB59rbOe7\no4uJDHRJW+cwwD80RZiiVaVEREREBqyQ+2T8B35vQi3+3IuvAetzvfmejBwz4o7PN0b4fkUR0WwD\nDecIAJ9simiYlIiIiEiOFLIn4xIggd978UbgDQDW2h4nep73NwVslwxh82MJvtAY4eaKImKZBhrJ\nHoxPNEWwCjBEREREcqaQQca3CliXjCDzYwk+0xjhlooiMhrsZAwfaYqwQgGGiIiISE4VMsi4z/O8\njO7mrLVBDaeSdCyOJahqiXJXeVFa5xvnuLw9zusjenmJiIiI5FohF9K5w1r7PmtteX8nWmtHW2s/\nBNye/2bJcPGm9jjnROME+tn7xTjHhITj2pZogVomIiIiMrIUsifjRvw9MO6x1m4ENgH7gEbAAKOB\n6cB5+DuD/x74RAHbJ0OcAd7fEuWb/WzU54zhfS0RwoVploiIiMiIU7Agw/O8k8D/tdbeAVyJv9LU\nNcAYIAgcBV4GngD+yfO8+kK1TYaPGXHH0kicTeFA6v0znGNKwnGB5mGIiIiI5E3BN+PzPK+O5E7f\nha5bRobL22JsKCru9fgVbTFyt0+4iIiIiHQ3rDc3NsZ8whiz2xizzRjzgDFm6mC3SfJvUSzBuHjq\nnoowcFG7JnuLiIiI5NOwDTKMMfOBfwYuc84tAFYDtw5uq6QQAsCl7XFMtwngAee4MBKndHCaJSIi\nIjJiDNsgA5gPVDvndiefrwaWDl5zpJBe1x7rkZYwhsvUiyEiIiKSd8M2yHDO3eec+zSAMaYUuAk/\n0JARYJyD86IJAh3b8znHlFiCuTFN+BYRERHJt4JP/C40Y8w1wPeBp4CVvZ23cuVKKisru6RVVVVR\nVVWV3wZK3lzaHuPZ5BRvA1zWrgnfIiIiIumorq6murq6S1pDQ0Pa+Yd1kGGM+QFwGfAR59xTfZ27\natUqli1bVpB2jTRr167tkXYidgSDYUxoUpf0FStW5Kze86IJwjiiGJwxrNDu3iIiIiJpSfVl+/r1\n61m+fHla+YftcCljzFuANwIr+gswpLBqamq4/6X/YVP7E3mtJwzMivnDpabEEozteyNwEREREcmR\n4dyT8VpgMvCiMad2ZTvonLt0ENs0IqXqnahrrmHB8pmsqMhdz8WmTZt6pIX37qa1KMioaJxNrV0n\ngy9ZsiRndUv6UvVs9SWXvVsiIiJSGMM2yHDOfQP4xmC3QwbX4jlzqCkLM6Gt52pTcuaoqanhj79J\nsMwu4sq3D9sOVhERkRFj2AYZMvKk6pk4F7ggFGBmLEG48E2SFFL1TJyodZysdWx4bgVf+KqhtFRT\n9EVERIYyfWUow1oQmKcA44yUSDgO7HP85vYEN3/79ISZG69zPHSfo/a4wzlNpBERERmK1JMhInnX\n1urYswt274Bdrzh2boP9eyEaAWMg0Wk0W3MT/HiVH1yUlcPMOY4582H2PMOsuTB1BoRC6ukQERE5\nkynIOIP0NiF2a41jwiQYP6HrjdVQnRCbcKc3xHPOcXpevgwnh/Y77v+TY+N6OLAPnAMM1DeuJdFt\nNeGTTTUpy6hrgJbmFWyrgXjcDzxCIZg933HBCsNV74BRFXr9SPp6+5x1LoExPTv3s/2c7XWBgy2b\nobgYZs/JST0iImcqBRlnuJqaGm79XoJzFi7iM/84PG6mTsQPn/q9Nn6ACaFpg9gayZdf3uZ44flu\niY4eAUY64p3yxGKwfTNs3+wIBAzvvm5AzRThxbUP0nLsaS6+4qsEQiV5q6empgZuuJ5Fo0fBvX/N\nWz3DgVahExn6FGScQbp/SMbjjmcec4we5WhvXMGEsYbZ84ZuoBFxbeyKbGBL+7PJFMPa1r8wr2g5\ns4rOI2zy98ddCu+66w2hkGPDemhv89OCQRhX2fvNQF/HgsHTwcboSrjwtfDmq3PZYhkJUt2Mnjzw\nAJGK8Sy/4FxCJRPyVg8AQcOKtlbQTXHGXt64gfZDj7L04r8hNHpuzspVQCODpdfX3sG9MG4SlHS9\nLxpqrz0FGQNQXV3dYyfEbEUjjqNHYP8e2LHVsWUT7NgOh4/6Q0QCAfjSJx1jxjoWLIaFiw2z58HZ\n02DseHI25CiX19ShKX6CPdGX2Rt9mQQJHB2TeR0JHNsia3klsp6Z4SXMLFpMeWBMTuu///77ueqq\nq3Ja5mDXNRSuaeYcw5e+ZYjHHXt3wY5tsPsVx87t/tyMaNSfj5Hqpds5oCgrhzkLYPY8f17GgnNh\n0lnpveZ7+wC/Z/VfefsVl/dIz8cHeD7eU4NZTyHrKtw1+Z9Jzg3Npa5Tvc6dc/zusce49rLLerxX\nBvI67+099eCDD/KWt7ylR3q2daXKF2+rpcVFWTL5EBVL3pdVuela8/LzvHhgIzdccX1e64Gh/95N\n9ZrYWB/j2TXP8/G3vKbHsVx/zg71f7/uampq4DufY9GbLofrP5/XuiC/16QgYwAG8h/zylaH96xj\n+2bYuxvq66CuwX+jBoKnh5R0H6te1wC79sKD98HY0f4bNRSGSZMds+bBuUsMF18CY8b2fQPW2x+K\n2267jXnz5vVIz+RDoS3RRF38ELXxgxyP7aPFncRgOgUXPSWIsyu6gV3RlygzlUwMTWdccArjgmdT\nEihPu+5UhsIN+ZlaTy7qCgZNMkAAMKxdu5Z43HFoP+zdBXt2OQ48V3Oqt2PseJhzLkyfZZg5B666\nekVO5+3cueYFft5ezvSmCBeMKspZuaneU9G449Yf/2TA76n+6nGJGDs23c2zT40lVDw+J/X0Jdd/\nlPL5eZQe/wPXxVtzXG43BVwt7YGXXuLXBma0tXJRaVle66q5/3848vC9kCLIyCUXa/F/xiM5LTfV\n6+nJ2BomnzMVa23e5w0Ot5vkppjjE39YD4R4XzTB6HBuFjLt9XPi1t8yZ/ZcAsH8BtNHmxLccsdf\nc/qZ1Gu+sjAr9m3JaY/nYHzODusgwxjzQeDb+Ne5GviEcy46uK2CO36c4M+/6/14f2PWO82bBiAW\nhYP7/cczjzl+fit8e5UfcGTi4fUb2DhrBW1xR0kws7zOOXZHN7I7spEW1wBwKrDYtm5Pl3P3bjnU\nazkLls+kxTWwN9rInqi/g3e5GcPsoqXMCC/OqE1y5goGDdNmwrSZ8JrLDO9472Ju+mIr02fBjV8w\nXf6oD+QPfPcPyYRzfGV3PcHyGOunL+Jjc8fk9Qbi73//Ajtn9OwxySXnErzw/F8oL26n8ehjjJl6\nDSYw9D/aX375Zeqb/4bHn0pw6evyu9p6rO04AG31GygaNTs/lcTj8PBDp5+vfgAuf0vqrrwMpboZ\nWN3URElDPdElS1gxbnyKXLmrix/cwKLZJXkdypFor6Nt/33+7y0HaD/yLMWTX523+tpoByBOnFCO\nbpV6u8mrr69PeSwfPU75qKsjb2vM8czxGDdtbCUwrRGAf2xcyD8vLeN1E0OU52FVwCcfe5no0bfx\n9Gp4fZ6/d/vmXRvYX35xfispsHvXbGDz5FfTHncUZ3jfl46h/5eoF8aYOcC/Ayucc/uNMf8DfBa4\neXBb5q+O050xMHH8ChLx1F92pRqrHgiCoeukWPCHl/Sn+wfKjpY4t29ooHHFdO4dt4gvzCpjQlH6\nf9h3Rl5kS+TZLml99Vz0x3E6kmp29Wxqf4IEiT5yyFDR2x+zUZPfzc+rf0cwhx90kYRjV3uMnW1x\ndrbFeOBEO16TPyTmt7Vt1MbqeGNlMXNKQswtCTGzOEg4kF39na8rEnfcsy/K9opGqJjPprFLed+s\nIkpz8Ed2+bLziLYdJtp6kFjbUdoaNnFs3B4YN57Fc+IY82dKxiwlXHoW4ZIpxKONBEKjsgqmCnWj\nkirfhk0JKirP4dGnl/Ppvw9SVJT7P4DxaBONhx4gEfW/GGk5voZg0VjKJrwaE8jR7jotLfDSC7Dq\n+3DvX/xPsUAArr0G/k8V3HAjLH2Vv+LUAEWcw2tt5Y9NJ7nvpH9N1Scb2BeN8vaKCi4uKaUkkMeA\nLdoO4YFfB/hfXCXa64g37SHetJtI7Yu4SNOp46077yZWv5nQ6LkER80iWD41p8F1FP/7yP3uELPM\n9JyV211NTQ1LrpjEnrYXmFlyQd7qAah54iEujx6A+uMwZuBzjxLOcaDF8VJ9jDt2RXjgUJSYg84f\n4ZtPJnjnk00EDbx5cogPzynm/LFBppcFCGTwmRSPOWZPt+zdAXu2w8E9jr07/M+JyjJYt9qybxPM\nmAdTZ8L0eYZDex3jJ0NRceafHd0/kyIxx7HiNoonw+iZyzlnUh7fRx3d+s1NfZ+XoR7XFHd84NFm\n2pbMoW7KBbx9Zu53FBu2QQbwDuCPzrn9yec/Bb7FGRBkvP/6AFUfdTQ1Qu0x/3HsCDy35nnq6+BE\nLdQdh5Zmf7iUI0EgAImE/7dp9BgYPwHOP38FEybBxMmGiZP9tHETobiPN1TcOY5GEjzw9PMcbI9z\nsD3BntY4D9VGie3eigP+40H4cdDwxvFhppcEObs4wBWvuYjpJQHGhU3Km5VEcrhBqmFRC5bPTNmW\n3tK76yjTDzy0f+SwZaJZBxgnYwk2tkTZ2hpjW/KxvS3GvhfXn3o1muQjvmvbqXwPAX+FU+cULzqf\nKUVB5pcGWVAa4pzSMOeUhlhcFqIi6L/24gnH4TbH/uYE+5oTHG5NcLTNUfOCx/G2BIfbEhxudcQS\nkDiwBeccn6qGzxqYXGo4qyTA+BLD4mUrmFTiP59W7j+mlJpTQU607QhtDZuItR8nFqknHq0nEW3k\nhY27O115AEiwbUdtt3+R/cmrcpy/eAqYIMFQBcGiMQTDlQSLxlNcMZeSigUZ/1vfvW49e5a8KuN8\n/WlpcRw4BDt3O55bm+Chh/wvFYyBD38yzhsvMZy32DBjmmHSRDJ6rTjnSMQaibfXsubZx4m11xJr\nPUis7QgA23Z2+vfb+BMwPydUOo1w6WQuuvgSQiUTCBaPJxAs7buieBwO7Ifdu1h7//3w3LPgPe93\nT5sANcllmAMk53785tdQfReEw6x469Xw5ivgnIUwazZMOdv/wE9xLbXxODWRdl6JRNgXi7Ju7VoO\nx2KciCdw+K+Klu3bAagD/gLcg//6rwwGWLLcMjMcZlo4zLxwEecWFzM5GMw8EI20wcZHTz9/8Cdw\n2QdhVObz6hLRRhKtR4m3HiXecpA1Tz2Ai5z0DxoDzrFlV123XIfpePcuWzKN0JhFhCpmESydTKBk\nEoGS8RgTTLsNzjkaaGRLYjuRZJDxdMIPpqeYyRSbgQ2vTBVM18cOsnjaTMYsaeNVlYspCuRmaFuX\numJROLIXdj7HovEhVhxYB+d9ECZN87+tTENrzHH/oSheXYwtJ+Nsa0ywvyVBJAGJV9YTMJBIfpAm\ngMT+LQB0DB1JAKtfgfufgsC8ZYQNTC0LcE5FgHNGB1k+LsRbzw5THjI881fHI/c4gkE4cRzqa/37\n7dqT/v+FCQLOH9XR0HJ6SHntVtiyDUzAf8uNr/D/DUpKHZXjYPxk/z5q6UWGq98HoXDvr/f2mONg\ng2NXXYKaI44/bTo9X+va29t499IQ508NMG9CgGljDJUlA5wb23QS9u2CV2rg57f5abEo/PBfYekK\nmL0AJp+d+lvqNNS3O+59Yg1H2hxHWhLsakqw5miMY7v8/6frbof5ow3njwsyrTzA5NIAb37dhcwc\nFch4ZEtnwznImAvs7PR8LzA1xXklAJs3b05dikvw49t+BMQxiQQmESMQjxKMttO68yX+4bprSQRD\nuGCYRCAIJsgNn/iE/ypP8wVnimDSdJh7cluX9Pa2APt2j+XR1S3MmHuSc5c2MHFyG4Gg/05evPh0\n1BlJwKGj/qNDmzPc0jSGP7WPIoDzb7Dw25RY7fl1d0x4xOAO+8OaEsBJ4E87Tx/7WmDBqfMDyXMM\n8PGyBj5c2gjGML5kMW3lJ4iUNhIrbsYF/JuEbV5yuFTyn6PX4VIuGXgkzzOJAKG2copbR1PcPJb6\n9gTbt29N69+0s6ampt7/f3OsUHWd6dcUi2U+gbahoYH169dnnO+EC/D+2CRi9Hy/xXZt65EWP7j3\n9O/djy06n/2ROPsjcR5tOD3uu4wE1aGjPNYwipsPT0zZjsjzG3ukJY52vPYNcWBf8gHwSOnSlOVc\nP6GOD0yoZ6L5XwKm59jzngEF7D3QkLIswA8yXJx41A9UOjQdfYTaxJXEGNvlfIf/7xLBcM/99xMB\nIkAzAeoCAZ4+dARMgM+vfpizE3EqXYJiHEVASzBEEY4wqb8OqK0rYsv20bS1B2lpCXKyMUzdiWI2\nvPgborHOQ+SguWkXAI0n/ZuImo3+DQLAjFkforw8xtgxEcaOiVAxKkZJcZwxYyIsXthAOBSnnC0U\nc5AgbQRpwiR7Qtc98/Lpevr493P4H6bnnHX69ZIgTJxy4pTSwjzamU75hpeY+e//Qsn+fV3yvxyN\nAwZzqms6zs5E997d5PN4BHvfvfCXe3q8ilvmLWDPV75O07z53DKqjJfDXb9tNM5Ru+lluovs29cj\nDaAViCxZwgttrbhuf6Nmx2J8ubHZ/6OYiBOKtBCMthGMthKMthJqb+aPD/yVcGsDRc11FLU3YBzs\naGgDBy//6FbMj24lWlROe/k4YqVjuObKNxMrriBWVEoiXEo8XEKsqBQX9G/Yx8ReYFz8pW7/IobN\nW3f3aPueg72/zpctmky0bgPRug2n/q4BtJmJHAm/gbgpIxaMEwlHiYXi/iMYIxqO8eD/PkA0FCNa\nFDuV8/gu////YM0BHmI1AMF4gKJoiCuuvZJwNEQ4FiIUCxKKBQnHQhRFwpjk/6AjgQvEcMEYzsRx\ngTh/+dN9OJPABeIQiBEPtXPgoP896J4th1nN8wTbRxGIhzEuyNXveAsmEcIkgt1+hk7VE246QdHJ\n44Tamgm2txCMtBJsb+G3T64h2N5GqL2J4oZjBBIJdjS2Y4CXf/ZTzM9+SjwUpn30ROLFpbzn9RcS\nKy4nXlxKoqiUWHEZ8ZJy2isnEyur4EfHx3JnfergMbF/S49xBqc++zqnJX8G5i0j6mB3c4LdzQke\nPBwD2nl3ZQOfGV/HL29ZTEtjz2/VOwcUHZradqVsE5wOMtpaoe0AHDngp2950VExuYaKcRFeOD7q\n/7d33mF2VeX+/7wzkwRCCkkomgAJVXqTjggiKmDhApcqoJciKOD9gQhYrlfAguXaUJogihowIgiK\nIihFEDFUBUIghYQqIQkECJBkZr6/P959kp0zZyaTZK11Zjbr8zznIWcfZn3P2ufde6/yFqbNW4U5\nbw5g7pttvLSgjXkL2njuH1ctbqe2dLporvdpwQuPcelE0Vm6WtfY8UiGD2xnxKB2RqzSzshB7Ywe\nspCd1nqV4YM6GDBvLkOeeJjWN99gwKsv0/baPNpem8eE+/9F2yvzaFvouxcCpi/swDrFo68vgu//\nH1bcN9RiHL75RiwcuRaLho+gfchwFg0dTvuQ4bQPXZ1XNt2GRbRyzZzVufe1wbzeabzW0corHS3M\n72xlwZ0PLu6PAZ0YnS/OWNyHyc+KJxb/TsaguZtjiOGtHQxt7WS11k6Gtnay4/zZtT9ZZkpQU8JA\ntJSY2Q+BqZK+V7zfCLhJ0kZ1/9+RwC93GD2Ue0/cKep3uuz+ZznhhslRNZZizdG0/eKhbj9WR7s/\ntVsWvUUAACAASURBVHuaDEnQ2op1s9qhJx+j45N7NfTxWmu9kYzbfDRjtxjDBluNYePtx/K2cV23\naNUpnp8xmykPzGTGI8/y5KPPMvPR53jxmZeW2cVMhpYWVjnmFAZ94CBfxV19FDZsONbadQ3FbV4+\nupT/rbV2tW21t6NXXqJz3kvQ2cmC313Fgqt/TMsG27Pq8d+jdexWRXuLwFoaXh/q7FgSYNXS+BqS\nOqGjA2vzB2rnv6fzxk9Op+PRO9hz13HstuO6vG3NIawxajXWHDWYEcNXZfiwVRi8au+3tRcubGfe\nqwt4ed6bzJ77Oi/Omc+s2fN5fNpsJtzwCAO23o51rvpNr9tbEV654TpmnXk6u7z7VkauuUe3/58k\nUEc3rpa+i7qslekH/3Ekwwbdym0TPt6w/c7O5XPkNKC1teuUac5Lr7PVPhfxlYGtfGFQ1++0SK5T\n20HrCRUvAwY0uB9/+s12LlALYy++kGF77924jU63JUlYox2QTvmtvhubB3jt7/cw84RPcPj6Q/jl\nRzZdxrd2OqXFK9gGtPbS3fBH9z/HKTdP5Q8XHsSu2y5Z/5PEonYf5iyPO40Qkv9NW9vS/X/X0eN5\nc+1RnH3H//aqrc6OTspjo5aWFqwX/XrqoZmct+PnOe27h7H/x7pmVeqOjvYOzIyWBnbWHecdewXT\nb36Upw/pfaxiR+cSm+zt77Sgo5Nh4x9m0ah1ad33BFpGb4yN3ggbORobuPQYU52dxT2v/gqz4h64\ndP+0aAGa8xz69zT07BTab/kJen4abS2rMXzwlqy+2pYMGbQ+gwasycC21RnYNopBA0axyoA1GdS2\nBi0NXBqlDha0z2XBotksWPQiC9rnsrD9JRa0z+H1BU/x8uuPMm/+wyxon8PQrQ/h7Yf/tMf+q7Oj\nayDs4m4ZWM87gPOf+DPP/vQAbt9yLfYc3vOYXFKXhS/wBZveXAv7PjqLW4ZuxLCv3tW9RmdHL8d9\nbd32q2PWDF49Y1uAj0oa39N3qvIk47PAupI+XbzfGzhb0vvr/r9RwAeAGcCbqb9nJpPJZDKZTCbT\nT1gFGAf8SVLXrfUSVZ5krA/cDrwTeBn4NfB7SZc383tlMplMJpPJZDJVp7JRtJKeBL4E3IPHZswA\nftLM75TJZDKZTCaTybwVqOxORiaTyWQymUwmk2kOld3JqAJmNtDMtjazjS1ixbBC42YzG25mm5rZ\nBDMbHUuvipjZGma2nZktM9vCCrbfZtY1WbuZrRNDr0qY2f49vQJrbWNmRxfumuXjB4fU6QuY2Xua\n/R0yzcHMRpnZsOLfHzGz08xs+2Z/r9CY2aea/R1CUTzf32Vm7zSz8AURutftXbR9H8XMugQ2m9nP\nA2sM7ukVUquB9lbFM2vDKO3nnYzeY2YbA5cAY4EPA6cDn5HUfV69FdfaAbgOeBUYDMwHDpYUPD2V\nmU0ELpV0WXHzOQXYX1KQUsVm9mu6pptYjKRDQ+jUae4C/A+wBqXkLpKCpBAzs92Bi/F4nwvwOizP\n4b/VQZLuD6FTaL0buBoYCMzCf5sZxWdPSFr+Ygd9BDPrpKttvAk8ABwvaflzFjfWeB6otVWesEtS\n43Q9y69zMvAZ/LvvBRwm6S/FZ0F/p2KAdw6wNXAz8C3JU6CY2XhJR4bSqtMdB3y8eA2XNKKH/31F\n2j8EOBBPN74Id3X9Ve089jea9TvFxMw+AXwHWABcBbwf+AuwD/BlSb9s4tdbYbpZcLgI+CSApD8E\n0vlmDx9L0lkhdOo0PwL8DJgGDACGA4dKmhhY52jgK0CtoIwBL0jaMqDGGOCHLLmmTpP0ZvHZHZL2\nDKRzAvAJYAvgkdJHqwCrS1ovhE6hVXsOLvVsqv1XUrByE2b2L0lbF/8+CD+XfwN2AT4t6bpQWpAn\nGcuFmf0JGA+cDWyDX0wbSzowgtbfgM9LuqN4vz9wpqS9ImhNk7Rh3bEpkjYO1P7Hevpc0s9C6NRp\nTgYuxQd8i428dj4DtH8vXlF+GD7x3FvSX83sXcDXJXWfp3P5tR4ATpR0b5E17WBJuxSfBfudivZ6\nXNkP9aAt6X0Nrxd2CV6357PAUDyW6qQQ57E4Z4fjO7fjgaslNS4isHI6k4HdJM01sz2BK4HNJL0e\n4Xf6FV6J7BbgLODu2uAkwoRmMHAIPrF4N14241RgvKTXA+qcD+wB/BSva9QBrIc/6H8r6fxQWqlI\n/Dt9i54Xc84MpDMZn1gMAR4Gxkp6xszejmebaVwApo9jZjOA1fA+1QZ7OwETCbsY0WhlfyBu50hq\nXIhn5TQnA0dJuq94vwtwkaSgJcbN7DHgKHyB76t4UeQZkn4cUOMmfHJxM/Bl4FVJ/1V8FnLcMhZY\nH/gF8NHSR53AZEkvhtDpQf+jePHoyyX9T8B2H5O0WfHvicAxkiYX/b1e0rahtCBPMpaL2mC87kea\nKal3pauXT2tqg5oeUVatzexO4BuSfl+83wm4WFLQ7W8zazTzfyPGxWpm/5QUvizxkvanSKq5sf1b\n0tqlz0IPHp6UtH7p/V+AX0i6IsLg9Wp8QPkoUJ+aLtiDtqTX5fvXrqnQ15aZbQIciU84XsQnHL+W\nNLvHP+x9+89KGlN6fwnwkqSzI/xOi23MzEbgA6N9JT0S+EH7E+A/8UH/NcXrd2V7DIWZzQQ2lNRe\nd3wI8IikcQG1nqJxISnD7XytQDpJfqei/U/ig5IbWLJrtxhJ5wTSWfxsMrPHJb2j0WeBtJL8ToXW\nCOByYCq+wNdeftbHohjwXwZMBj4ladYy/mRFNB6VtEXdsRkhr6mizSckbWJmX8Sze94N/CvwTkb5\nmloV+BfwX5LuinBNtQDflnR6qDZ7obkOvug2GjhW0oOB258kafPi3w9L2qr0WdDrF/DiH/nVuxdw\nL756M6l4Pxqf0cbQmghsX3q/GfBwJK3NgcfwrdRJ+IBilwg6z+Crk7OBl4p/zwVmhtbDtwAPi2gL\ntRU8gK1Lx8cCTwTW+hfw4dL7rXD3n81CaxXt/xb4SKxzV6f1WN35WweYgg8spkbU3bG4xhYGbPMm\n4P+ATYr3o/AC38eF7ktxLY0ovT8WuAvfrQlmE8DrwIPAAcCA4tj0SL/J08CwBscHh9YERha2d3xx\nzS716m+/U6n9y2Nfu7hrxcFAa93xDwL398ffqU7zlOI32ojiWR9JZ1Xg+7ib7X9G0hhcvL6Fu3Ku\nVrw+g7vuhda7D9gPOBTfzRgOzAqsMQMYV3p/AL4oNjjSNXUrMCaWHdRpnYyPj74EtEXSeKpo/yjg\nNuB9xfH3AA8E10tx4qryAvYvbrCzcR/8Z4AjImntjg8kbwX+hPvi7xexby3AO4BtgUGRNC4r30yB\nE/G0wruHNu7iQurEYyZm4SvXwW52+Gr/HHxVsnbsXNw14rgItjATuKd07Kiib29E+J0OoTTBjfkC\n3lc8ZG8Gbix+q2OA84CvRdDbAfcnfxp3Yfl4wLbfBvwcuLF0bGt8Na89cD/OAZ4Aziod+23x0Ahp\n58Nwn/T78QWBy4DnItnCmfhq7lmFDf4HcBo+yflcBL1vEn9AnuR3KrW9PbBe5D5tATyE+/TXjt1A\nhMWiVL9TA83tCpt/JlL77weexHdT14jYj058Ma+z9Kq974igt0txDxyMD/zfAL4SWOOk4nn+g9Kx\nH+GT0bkR+nQTPua7DphQewXWeAdwJ76QvWUseyi0dsNdX88r7O/T+GT+BeDdofWyu9RyUvit7YcH\nT90m6ZFl/MnKaA0HdgUGAX9X4G1UMztT0je78+VVIB/ekt50SRvUHXtK0nr1LkEBtBq62UiaGVBj\nTXzw+FLxfl9gpqTHQmmUtAwYpZJrT2Ef75R0a2i9lBRuCjU7f0DSTDMbIum1QO3X3KSOwHfQxuPB\nxC+EaL+X3+Fdku4K3OZ7gQ1U+DubWSu+4rubpB7joFZQb3t8V+ZIfAHkB5IuDqyxJ75KviF+T5qJ\n/1Z/DalTaLXgLjdRH4Kpf6duvsPukv4Wsf234QO/XULrpPqdGuiuBmwj6e6Q58/Mforv+vw/fOC6\nFAoY59RMimfW6rXnY+C2N8bjYf9QOvYBYHdJXwqs1fAaVcBYUjN7A3gNn6C1138eeizWzXcw/J4U\n9vrNk4zeUwQmfrtuoPd9Sf8dUKPHC0TSuQG1TpV0QTeBaCiQD29J7358RfLm4v12eFDVXsBdKvn2\nBtBqwYO1NgO+h7uwBBvk2TJSG0q6MGv1Ss/wLeK98VWv6yVNCKzRie+WTMBdApe66YXuUzff4RVJ\nw2LrFFoXSoqWetM8TfOh+I5dkEwuy9BbH99xip4KM5VWyj4Veknsr0p2XqcVrF/F/ahGLaPQ4v9K\nag2h81bAzNrwoHnw8/dDFUHggXXGAlviniWtoRbASu33eB8IPRbr4XsEv36DpcV6i3AacJiZHSjp\noeLYfkCwSQaeWQd8cLwVvrXeiqd1vDmgDsUEowX3tf5iyLa74QRgQnGTXQSsibvGnIrHUITkfNxV\n5R1F25ea2SWSvh+o/R0DtfNW1oIlv9MV+Bb7yWa2iaSvBNS4En+Ij8DdpcqkWmWJVuemAUcBQQdf\nxaruRribxeOSrsTPaxTMM1odim/r74GnSO3XWnU678IHLKlIZX/92s57IFi/JCWrT2bLqLEQatck\nlU6d5mdxl58BLJmo3RJB5xjg6/i9bx/gDjM7QtJtAWWmSOpSj6MJBL9+807GcmBmU3B/4YuBUyRN\nCJ3NoKR1L14P4cXi/dp4ysj3RtC6FTha0rOh226g1QZsgrvGPB5ra9jMpgGb4pktNjOzkXhQYvCs\nOJkVx8ymApuqyChkZoOAB1Vkv4ikaU1wvUi5wvuqpKHL/j971VYbHsx+Il6/ZAF+7f4I+B8VNR9C\nYV4T5ljgIDzofBCwk6QpIXVSaqXsUw/foYo7GcHsvBdaQftlnsFxHO4G/XTxfBoC7CrpVwF1nsJr\nzsCSAWTwXZNUOnWaU/B4xW/jZQU+AoyU9LXAOo/iQdF3FGOJ7fDsjlss40+XR6NP1LvKOxnNR5Ku\nLQaw1xXGFmvlZg2VUrtKesHMNujpD1aChcA/zVPZLippBimSZz0UHzKzWP6G8yUtsqJQurx2QfCB\npZk9SdfV8Dfw+hyflfTvrNUj0tIpSxfSOGXlSlPEZvwE2MHMFgJ/BE4NFevUw2peytVdCLs7cz6+\nq7qlpKkA5pVhL8SDmkPmb5+G72Zdiw8Y7sSzcsWYYCTRStynJPZXUTtPef7OxgfFTwDjzOwXeEar\necArQLBJBh7bdDdwnqQbArbbLJ0ykjTLzCYB75B0sXk9kKCTDHwxflZpLPGgeercoBqB2+teKPH1\nmycZy4cBSPqnme2M54yPtTL+V/PS9RcV708E/h5J66riFYszcJ/4a/EsDSkuqFvM7PPAQPPq3CcS\nxz3hN/jK5KX49XQavlr5NJ5O8oNZq0f+YF4R/gJ8gnEiPhCLwRX4NVtL3HA87vKzb6D2X6Nr1dYa\noQdE3bmJGN63UByOTzBerh2QNM3MDsOzPgWbZODZTcbhdvCapI4YCwOJtVL2KZX9VdHOIV2/TsLT\ndj9lZtvgaV83lTQtoAYAxULbdfguZDRS6dTxgnlV7knAgWb2PO4SG5rJ5kV2BWBmxwPTA2usZV6T\nqCGSjg2olez6hewutVyY2doqZaQpXAkOkRR8gF4EV34OT5vbivsanitpfmitQi9aYJOZ7YBnpDkQ\nvyFchVfwDRo8VafZhucCPxAfJP8ZOEfSG4F1lirWZJ45ZpqkcWb2tKR1s1aPOgPx1KUfwic1td8p\nuG00cm20wFnNUmFmV/T0eajgx57OjxXFSUPolNrcDJ/8HY2v7I7CJznPhdRJqZWyT1UjlZ2nxOoK\nnsW4jt4KmNlG+HX1ReD3wM54YpnQGe/Wxheo9sHjMh4EjpT0ZECNOcAPuvs8VeB3DPIkoxfYklSv\nDd1+Irn71LRXpTTjjBRA1SWwCa//ETKwqZZJaE88lej++M7MVZK6pPELpDccD57vxP38g04wCo3p\nwB61eBYzWx13KdoQd4sI9vCompZ50oFzEyUdwMyuAr4n6R/F+12AL0naP1D7PbajUrrF0MSKM+nJ\nVzimH7GZDcAXCI7Hg7Gvl3R4f9aKrZPK/qpo50Xbqc7fUteNlSowx8TSZWJKolPSSxJjV4zFBkh6\nJULbyWIyUl+/2V2qd9QGp+VdhDF44bWOGIJmdjJehKjmn17LnhAjvd1ZePGhOyQ9bmb74allgwU2\ngTtQArcDt5vZ1nhhr2uI0Ccz+xDwM+BZ/LyNNrPDFL6mxJeAiWZ2I+4SsS+eGOAswrtnVUpLUqeZ\n7W5mYxQx6YB5EgXhmdv+amYPF++3wauph+L3eP2Ix2vSpc8EBB98xY4zAUZZ92m1RwbS6IKkRRRF\nr4pd1pDuAk3RSqCTyv6qaOfQhH6lwpbOxFQjRiamJDqFVs0mdjSzBfjv8+mAMXYtwOnARHmdnv8G\njjUPBD9O0twQOjW5gG0ti6R2nncyekGxXXY5XihlAn6D2xXPtPJ+LUlnG1JzBnCQpAdCt91Aa5Kk\nzcsuMtagcF4AnTG429SR+AT3ajxjVrBtx5LWFOCo0qr1HsAlMVaMipvdPri7zz/kxZvWBZ6vC2rO\nWl01bsLTykZJOlBo9FjLQdIdgXQ+i8cwtOAF/66W9HSItnvQ/Bs+Ub+MJXEme0sKEmfSDHcVM3sf\nsDm+Azmxdg3HIJVWCp1U9ldFOy80Up2/N4CHa2/xxbxHin9L0k4RNFNlYkqiU2jFvvd9FV98PRXP\n/PUXvD/vA9aTdFwInUJrL0m3h2pvGVpJr988yegF5tkfpgHfwC+gH+KroLsBn5EUMuC2pjkxxs2m\nG61rge8AlxaTjeOBwyXtE6j9T+ATi/XwSdp4SSFXkBtpNqouPkPSuIAatTSE98iD+KKkIay4Vopq\nqj1OLCVNCqVV6NUqjB+OV0EeD/xapSKeAbWSxplEdlcZA/wO3yF5CF/E2QnfjTxA0pz+ppWyTyXN\nJPZXVTuP3a8Gix71xUFjVLd/QtIm5pmt/iHpNjObLGnT/qhTaEW1CTN7HNhe0nwzOwPYUNInzcyA\nJ+q1V1KrxwK0IRfdSppp7hN5krFsrBSYZWbnAG2SvlC8j1Un46t4poRLcXcVIPyAqNCKGthkSyou\n34PfUOtvqjEuoP/FH+gX4zP2o/H6HGfgdr9SsS1Wl4YQdy9bnIZQYeMjqqq1BV7L5H5JM0K120Cn\nlo63flt4FDBEESvsmtmOeE2JbSUNXNb/vwLtR40zKeksdlfB70fB3VXM7BpgKvB5FfU3igf6V4H1\nJR3R37RS9qkb/aj2l0onlZ030I3WLzP7KG4HI/B7U20nI3idEfP09Ffi2c72wVNQ3yFp7f6oU2jF\njrFbHCdhngHxekm/KN4HDdY3s3bc/X4CvhixlPtUyEW3bvTjXb+S8msZL2By6d9/wlegau+fjKT5\nZIPX9Mj9XBUYFqHdj5Vex9S/IvWls4dXR4D2Z+BbpuC7WovwlY4YfamcFr4FPQ+YCMwFDo3Rn260\n24Av4Df14yJp7IDvDj6N+yR/PHD79xbn7jF8Mn1fcWwhcF+E/vwNT2M8FF+VPxO4KbDGLKClwfFW\n4On+qJWyT3XtR7W/VDqp7Tzl+QNmFjprl1+R+rMRXvOmDbgJeAk4qT/qpLIJfLF1FO6KNQvYoDi+\nHjApcJ/WLp6Jd+Ouc1+M8dxtoBvdzvNORi8ws9vxh+oz+ArvWElzzCu5ni9pt2Z+vxXFzLbCVyc3\nxTM9nSBpZkS92srNSJbM1KVEFWJDYgnTEFZRq9hdeK+k6ea54q+QtH1onQa678Tjq54BTlTAgPPS\n9vMR+MN1PPArldJeB9RKEmdS0ovurlJve8vS7w9aifuUxP6qaucp+1XoPShpuxht96CZKhNTTLfK\nnmxCCuRuZmZH4kHs84CXJe1tZh/GM3FeIumCEDoNdMfhNngEXnfkKknfDdh+UjvP2aV6xxnAb4E1\n8HRscwq3ktPxHyoYZvZzSUcX/z5A0vWlz/4oab+Acj/Gt+duBA7DY00+HLD9er4GHIzPmqOQ8Px1\n1r1fEKjdt4pWp6TpsLi45ZBIOgCY1505D99NO0PSlRFkJuNugRPwGC4BB9uSSrEXBtR6MWBbveE+\nM9tZS7smPBZYo972yoQesKTSStmnVPZXVTtP2S+Ai8zsN3g2n3LSi+D3JouciSmlTnliWQzIaxWs\nDS/quncgnfFm9hC+O/Pn4vAY4OuSfhlCoxvdGWb2Ddxt6hzgXCDYJIPEdp4nGb1A0n3mWXWGSZpX\nHL4buEbS1MByO5f+/XXg+tL7hitiK8HbJX0LwMy+jO/SxGQu7n8fcyUl5fnLrDj1qZ+DZsYqY2bv\nwWObHgG2irVig/siC/ex3qHus9A2fyM9xJkQKC20pUv/C7COmTVKj2z4w70/aqXsUyr7q5ydF6Ts\nF8BH8RT15WxIKr5HaK7AMzHtx5JMTFfWafcnHczs28BRwHB8Mjqm0A+GpElmth0wycxG4PGdmNlF\nMTwwzGw33C7+A3fX+j6+wB2SpHae3aX6GOUtdOtadTlowZb69ur1QmOeZeoDRFy5SXX+LGEawipq\nmado3o8lg4cb8F00g7AJDorEA3OB22hwE1WExAOFbhLXhEKrDa9hcirwBUmXB2o3pbtKOdNYo98p\n5H0iiVbKPjXQ7veuMQ20oth5N1pR+2WBsx0uQytJdq5UOkW70/Biu98FvgJsiZcUOC2wzkwaeGCE\nXKwys6/h3iTP4u5LExS2DkdP2lHtPO9k9G1SzwBj66Vcuam1HYv6lZmstfzUF/2pvRcQskZLuZZD\nikFXisJhZb1ynMk7Q8aZkNBdRUUGle5itwh4n0illbJPNVLZX8XsvKyTql93mdn2SlALizTujil1\nwF1uXzezqcCmkv5oZj/AE1SEZK6k+wK3Wc/ZeEauVYHj8KJ/iz8MuZBYI9l9Iu9k9C1s6bRpk1Qq\nHhdhJyNpUaAUKzcpz1/RZso0hJXUSkV3A73QfbIEhcMKnehxJtaE9L8pVg5TayXuUyr7q4yd1+ml\n6tcjeFr1KcD82vHAO9Nld8cN8Of9YndHSfXuMn1ap07zj3hWqTuBY4GfAZdLCuqGmMgDI2kij0Iz\njZ3nSUbfom7gvyU+6IdiEiBpcMM/XDGtesOur18RtCiQeVHD78RcuUl5/gq9yg2IUmhZ4gJ5hWaq\nAWWKTEzlOJOTIsaZ1OtGd1cxswfxIlgpdp2SaCXuU2VcY5ph5wnPX8OBZWAXxFSZmJLo1GmuARwm\n6UdmdjGwC3COpOsC69yBe2A8VTqsGG62ZrYhPm6ZDzyoCMU6S1pJ7Dy7S/U9krmr1G5m3a1aA6FX\nrbcF7jGzaCs3pHX3gTRbqVXUurHBsVgBnTVSJB6ANC4Df8H70w5cUN5ah2gFLpO4qwAXAdeYWfSs\nOwm1UvapSq4xye2cROcvxup0TxoWNxNTEp06zdl4ATkknRS6/RJjE3hgDMBdlw7CE/AMxJNGfFvS\neZFkk9h5nmT0MZow8IcEqWULTo7cfjPOX7I0hFXSql8tqVshPz2ERgOiDvQsbSamZHEmqd1VSBu7\nlUoruk4q+6uwnafsV1IsQSamxDrHAL8DfgWsXvqoXeHrlqWInfkisA5e8O8FADNbH/iZmc2WdFEo\nodR2nt2l+iiJXWOSFwWKTULXmJRbqVXVilYgr04nap+a5FcbNc6kSe4qKbPuJNFKFI+WxP6qaOeF\nRvJ+pcLSZWKKrmNmJwEnAHvitSTOBHbC4zL2kRR0oJwodmYyHg/xXN3xtwN/krR1QK2kdp53Mvou\nKV1jUq6QpyLV+Yu+lVpVrSaskMfuU+oCeRB/F7IZ7iops+6k0kqhk8r+qmjn0Jx+pSJVJqYUOp8E\nPiTpNTNbJOla4FrzejTnAQcE1IIEHhjAoPoJBoCk581saGCtpHaeJxl9l5QD/9SpZVOQ6vxVcUAU\nXcvSFcgrE/v8pSwcViN2nEnS9L8FKWK3Umul0Ellf1W0c2hOv1Ix1czOpcjEZGatLImb6G86gyXV\nJpvTawcl3WxmwdyKSu2m2MEaaGYDJC0qHyzciFcNrJXUzvMko++ScuCfcoU8FanOXxUHRCm0mrFC\nHrVPVYwzURPqPJBm5TC1Vop4tCT2V0U7L9pqRr9ScTSeiekWMzsYOB84pZ/qrGrmBeQkfbB20Pwh\nMiCwVipuBb4AfLnu+KeAoJOc1HaeYzL6KIn9kqOnlk1NQl/r6GkIq6hlTaiEnPj8VSLOpKSTLEYs\ns/IktL9K2XlJL0m/MsuPmY0H/i7pgrrjRwJHSPpwc77ZimNmY/Ddn8fxBbgFwF7AzsAekp6MpBvd\nzvMko4+ScuCfIrApNVWcOFWRFAGdKUkdZ5JwMp2szkNmxUllfxW289RxYtFJlYkpoc46+ID8AeAu\nfGdrV7xg3l6SHunhz/ssZjYYOAmv92F4/y6Q9FoErWR2nt2l+i4pXWNSuiekIuX5y6w4qdInR6ei\ncSY1UtZ5yKwAqeyvqnbepH5FpZSJ6Vq8EvdSmZj6mw6ApGfMbDs8AHxPvKbEw3jtnpkhtVJSBMy/\ngJ+3kbir9+fMjMBZ1JLaed7J6KOkdO2oIvn89Q+qtEJuZp14nMltNHYBi5H+N8kuZGp3lczyk8r+\nqmrnzehXbMzsn3gmpqfN7DFJmxXH3w+cLClIJqZUOnWa9bvgAoIOyFOTwi01tZ3nnYw+Sh4Mrxz5\n/PUbqrRC3oxMTKl2IauYHKJqpLK/qtp5M/oVm1SZmJJmfCqozC54iRRZ1JLaed7JyGQyTaOKK+RV\nizOBHOPUn0hlf1W0c6hWv8zsGWDd+kFrkYlppqT1+pNOXduV2QWvYWafwGNLoi+6JbtPVOj3yWQy\n/YyUWdRSUcVMTFVMDlFVUtlfFe0cqtWvVJmYmpHxKeWAPBUpF91S2Xl2l8pkMs0kZYHBVKTYxGk7\n4AAAAdJJREFU8k5NFZNDVJVU9ldFO4dq9etM4E4z24sGmZj6oU6ZKhYRTumWmsTO805GJpNpGlVc\nIa/iClum/5DK/qpq51Xrl5mtjmdi2pklmZguDJ2JKZVOSa+Ku+ApSxekuU/kSUYmk2kWVcwCVsU4\nk0z/IWGBxkraeRX7lSoTU8qMT1WME0u56JbsPpEnGZlMJhOOKq6wZfoPCQvXVdLOq9ivKsbpVHQX\nPNmiWyo7zzEZmUwmE5Yqxplk+g+p7K+qdl7FflUxTqdycWKJd/CT2HneychkMpmAVHGFLdN/SFig\nsZJ2XsV+5TidTD3J7hN5kpHJZDLhqGKcSab/kMr+qmrnVexXjtPJ1JPsPpEnGZlMJpPJZDLVJMfp\nZJpFS7O/QCaTyWQymUwmGneZ2fYV0sn0E/JORiaTyWQymUxFyXE6mWaRs0tlMplMJpPJVJdUmZgq\nl/Eps3LknYxMJpPJZDKZTCYTlByTkclkMplMJpPJZIKSJxmZTCaTyWQymUwmKHmSkclkMplMJpPJ\nZIKSJxmZTCaTyWQymUwmKHmSkclkMplMJpPJZIKSJxmZTCaTyWQymUwmKHmSkclkMplMJpPJZIKS\nJxmZTCaTyWQymUwmKP8fhb2Y2TfvTGgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c88e5d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rcParams['savefig.dpi'] = 100\n",
    "\n",
    "\n",
    "d = load(get_path_from_name(gene))[:,reordering] - 1\n",
    "binary_pattern = pattern_dict[gene]\n",
    "maxy = 1.3*percentile(d[:,argmax(d.mean(0))],91.)\n",
    "\n",
    "fig = figure(figsize=(8,2.5))\n",
    "\n",
    "r = violinplot(d,widths=0.9, showmeans=True)\n",
    "\n",
    "for b in r['bodies']:\n",
    "    b.set_color(cm.rainbow(random.rand()))\n",
    "    b.set_zorder(100)\n",
    "    b.set_alpha(1)\n",
    "for b in [r[u'cbars'],r[u'cmins'],r[u'cmaxes'],r[u'cmeans']]:\n",
    "    b.set_color('k')\n",
    "    b.set_alpha(1)\n",
    "    b.set_linewidth(0.3)\n",
    "    b.set_zorder(100)\n",
    "\n",
    "xticks(arange(len(order_predictors))+1, order_predictors, rotation=90)\n",
    "\n",
    "bar(arange(binary_pattern.size)+1, (2.01*maxy * binary_pattern),color='0.1',\\\n",
    "        width=1 , alpha=0.2, linewidth=0, align='center',zorder=-10)\n",
    "\n",
    "xticks(rotation=90,horizontalalignment='center', verticalalignment='top',\n",
    "       position=(0,0.02), fontsize='small')\n",
    "yticks(fontsize='small')\n",
    "\n",
    "ylabel('Posterior\\n(Molecules/cell)', fontsize='medium', position=(-0.5,0.5), color=(0.2,0.2,0.2), alpha=0.8)\n",
    "title(gene, position=(0.46,1), fontsize='xx-large')\n",
    "\n",
    "ylim(0, maxy)\n",
    "\n",
    "tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  },
  "toc": {
   "toc_cell": true,
   "toc_number_sections": true,
   "toc_section_display": "none",
   "toc_threshold": 6,
   "toc_window_display": true
  },
  "toc_position": {
   "height": "306px",
   "left": "1131.015625px",
   "right": "20px",
   "top": "113px",
   "width": "189px"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}