data_release_companion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import h5py\n",
    "import numpy as np\n",
    "from pycbc import conversions\n",
    "from pycbc.cosmology import redshift\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import rcParams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "rcParams['font.size'] = 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_plot_data(postfiles):\n",
    "    outdata = {'mchirp': [], 'mchirp_src': [], 'q': [], 'snr': []}\n",
    "    for f in postfiles:\n",
    "        # extract posterior samples and injected values\n",
    "        with h5py.File(f, 'r') as fp:\n",
    "            samples = {p: v[:] for p, v in fp['samples'].items()}\n",
    "            injvals = {p: v for p, v in fp['injections'].attrs.items()\n",
    "                       if isinstance(v, float)}\n",
    "            lognl = fp['samples'].attrs['lognl']\n",
    "        # calculate max SNR\n",
    "        snr = np.sqrt(2 * (samples['loglikelihood'].max() - lognl))\n",
    "        # convert injected values to mchirp, q\n",
    "        injvals['mchirp'] = conversions.mchirp_from_mass1_mass2(\n",
    "            injvals['mass1'], injvals['mass2'])\n",
    "        injvals['q'] = conversions.q_from_mass1_mass2(\n",
    "            injvals['mass1'], injvals['mass2'])\n",
    "        # convert chirp masses to source-frame\n",
    "        samples['mchirp_src'] = samples['mchirp'] / (1 + redshift(samples['distance']))\n",
    "        injvals['mchirp_src'] = injvals['mchirp'] / (1 + redshift(injvals['distance']))\n",
    "        # store param stats\n",
    "        mchirp_data = tuple([np.mean(samples['mchirp']), np.std(samples['mchirp'])])\n",
    "        mchirp_src_data = tuple([np.median(samples['mchirp_src']), injvals['mchirp_src']])\n",
    "        q_data = tuple([np.mean(samples['q']), np.std(samples['q'])])\n",
    "        outdata['mchirp'].append(mchirp_data)\n",
    "        outdata['mchirp_src'].append(mchirp_src_data)\n",
    "        outdata['q'].append(q_data)\n",
    "        outdata['snr'].append(snr)\n",
    "    # cast to arrays\n",
    "    for p, v in outdata.items():\n",
    "        outdata[p] = np.array(v)\n",
    "    return outdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bnsfiles = glob.glob('bns_posteriors/*.hdf')\n",
    "nsbhfiles = glob.glob('nsbh_posteriors/*.hdf')\n",
    "\n",
    "bns_data = extract_plot_data(bnsfiles)\n",
    "nsbh_data = extract_plot_data(nsbhfiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAAHqCAYAAAAK1a8LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAIABJREFUeJzs3X2cnXV95//XmxDscJe0wCYmmIpV4g3YRMaqIDcqaSiVGtF6x9rSrWSRLdV1jWsqu6uuK1RYCmqRYv1tvGG3UI35EQWCgIIVKiSbKPGhof5guZmQSJCglFBj9vP745yRyWQmOZPMmXNm5vV8PM4jnO/3e32vzzWcyZXPub7X50pVIUmSJEnas/06HYAkSZIkjRcmUJIkSZLUIhMoSZIkSWqRCZQkSZIktcgESpIkSZJaZAIlSZIkSS0ygZIkSZKkFplASZLUJknOS3J/kqeTrEly4h7Gn9wc93SS+5KcO6j/w0lq0GtTe49CkjSQCZQkSW2Q5K3A5cDHgfnAHcANSeYMM/4o4PrmuPnAhcCnkrxp0NANwLMHvI5tywFIkoaUqup0DJIkTThJvgt8v6rOGdD2T8CXq2rpEOP/Ejizql4woO1vgZdU1aua7z8MvLmqjml3/JKkoXkFSpKkUZbkAOA44KZBXTcBxw+z2auGGL8K6E0ydUDb85JsbC4N/LskzxuVoCVJLTGBkiRp9B0OTAE2D2rfDMwcZpuZw4zfvzkfwHeBs4HTgHOa29yR5LB9D1mS1Ir9Ox3AZHb44YfXc5/73E6HIUlda82aNVuq6ohOx9EtquqGge+T/CNwH/DHwKWDxydZDCwGOOigg4574QtfOBZhStK4M5LzjQlUBz33uc9l9erVnQ5DkrpWkgc6HcNe2gLsAGYMap8BDFc1b9Mw43/ZnG8XVfVkkh8ALxim/yrgKoDe3t7ynCNJQxvJ+cYlfJIkjbKq+gWwBlgwqGsBjSp7Q7lzmPGrq2r7UBsk+TXghcAjex+tJGkkTKAkSWqPS4Gzk7wryYuSXA7MAq4ESPKFJF8YMP5KYHaSy5rj30XjfqdL+gckuaT5rKijkrwC+DJwEPD5MTomSZr0XMInSVIbVNU1zeIOF9B4XtN64PSq6l8mMmfQ+PuTnA78FfBuYCPw51X1lQHDjgT+F42iEo8C/wi8csCckqQ2M4GSJKlNquoK4Iph+k4Zou024GW7me9toxacJGmvuIRPkiRJklpkAiVJkiRJLTKBkiRJkqQWmUBJkiRJUotMoCRJkiSpRSZQkiRJktQiEyhJkiRJapEJlCRJkiS1yARKkiRJklo06RKoJOcluT/J00nWJDlxD+NPbo57Osl9Sc4d1L80yd1Jfpbk0SQrkxzT3qOQJEmS1AmTKoFK8lbgcuDjwHzgDuCGJHOGGX8UcH1z3HzgQuBTSd40YNgpwBXA8cBrgV8CNyf5jTYdhiRJkqQO2b/TAYyx9wHLquqzzffnJzkNeDewdIjx5wIbq+r85vsfJnkF8H7gKwBVtXDgBkneCTwBnACsHP1DkCRJktQpk+YKVJIDgOOAmwZ13UTj6tFQXjXE+FVAb5Kpw2xzCI2f6+N7GaokSZKkLjWZrkAdDkwBNg9q3wycOsw2M4Gbhxi/f3O+R4bY5nJgHXDnUBMmWQwsBpgzZ8iVg5IkSZJatGJtHxev2sDGrduYNb2HJQvnsmj+7Lbtb9JcgRoLSS4FXg28qap2DDWmqq6qqt6q6j3iiCPGNkBJkiRpAlmxto+ly++hb+s2Cujbuo2ly+9hxdq+tu2zLQlUkk8kmT7g/a8nuagd+xqBLcAOYMag9hnApmG22TTM+F825/uVJH8FvB14bVXdt8/RSpIkSdqti1dtYNv2na9bbNu+g4tXbWjbPtt1BWpBVW3tf1NVjwO/26Z9taSqfgGsARYM6lpAo8reUO4cZvzqqtre35Dkcp5Jnn40OhFLkiRJ2p2NW7eNqH00tCuB2i/JIf1vkhwKDFd0YSxdCpyd5F1JXtRMfGYBVwIk+UKSLwwYfyUwO8llzfHvAs4GLukfkOSvgT8B3gE8nmRm83XwGB2TJEmSNCnNmt4zovbR0K4E6nLgH5L8RZK/AL4N/FWb9tWyqroGeC9wAY1CD68GTq+qB5pD5jRf/ePvB04HTmqO/xDw51X1lQHTnkej8t4tNIpK9L/e39aDkSRJkia5JQvn0jN1yk5tPVOnsGTh3Lbts11V+K4G7gJe03z/jqr6AUCSL1bVO5PcDdSg7YpG+e9PV9XX2hFYVV1B48G3Q/WdMkTbbcDLdjNfRi04SZIkSS3rr7Y3llX4Rj2BShJgbVW9GFg/xJAPNP988zBTHE4jAWtLAiVJkiRp4lg0f3ZbE6bBRj2BqqpK8r0kL+m/6jSo/5Hmnw8kOQF4Lo3nM/X3fyHJcMmVGPta95IkSZIa2rWE7yXA2iT3Ak8BoZFb/U7/gCT/i8aDatfSKC8OzSV9VTXUlSvxTK37/nKN/bXuAZMoSdKo8gs7SdpVuxKoM1oY89vNZX4agd3VuvekJkkaLX5hJ0lDa1cVvqOBx5rV7X6fRvW6wbUE70rSvvIYE1Qnat1LkiafTjycUpLGg3YlUJdU1ZNJXgmcBdwMfG7QmHnA95v3S92V5O4kd7UpngmjE7XuJUmTj1/YSdLQ2pVA9VsEXFlV1wIHDup7A40rVX8A/CGNqnwWj9iDTtS6lyRNPn5hJ0lDa9c9UBuTfJHGA2jnJXkWzUp7/c+BAr7Mrs+BAvidIdrU1Ila95KkyWfJwrk73QMFfmEnSdC+BOrNwELgo1X1eJJnA+9v9u3pOVDag7GudS9Jmnz8wk6ShtaWBKqq/hlYPuD9I8AjA/6bZoEJSZImrCTnAUuAZwM/AN5bVd/ezfiTgUtpPA5kI/CJqrpyX+bcF35hJ0m7avc9ULuVZE6S4zsZgyRJ7ZDkrcDlwMeB+cAdwA1J5gwz/ijg+ua4+cCFwKeSvGlv55Qkjb52LeEDIMnfAy8GHgb6mn8OfPXQWNK3qJ1xSJLUAe8DllXVZ5vvz09yGvBuYOkQ488FNlbV+c33P0zyChpL4L+yl3NKkkZZWxMo4O3AZ4DfAw6j8Uyow4HwTAGJb7Q5BkmSxlSSA4DjgEsGdd0EDLfy4lXN/oFWAX+cZCqNc+dI5/yVxx57jGXLlu3U9pKXvISXv/zlbN++nauvvnqXbebNm8e8efN46qmnuPbaa3fp7+3t5ZhjjuGJJ57gq1/96q4H9KpXMXfuXLZs2cLXvva1XfpPOukknve857Fp0yZuvPHGXfpf97rX8ZznPIeHHnqIW265ZZf+0047jZkzZ3Lfffdx++2379L/+te/nsMPP5wNGzZw55137tL/xje+kWnTprF+/XpWr169S/9b3vIWDjzwQNatW8e6det26T/rrLOYOnUqd999Nz/4wQ926T/77LMBuOOOO7j33nt36ps6dSpnnXUWALfddhv333//Tv0HHnggb3nLWwC4+eabefjhh3fqP/TQQznzzDMBuPHGG9m0adNO/YcddhhnnHEGACtXruSxxx7bqX/mzJmcdtppACxfvpyf/exnO/UfeeSRnHrqqQBce+21PPXUUzv1H3XUUZx88skAXH311Wzfvn2n/qOPPprjj298LAd/7sDPnp+97v3staKtCVRV/RI4J8l7gCer6nPNk8ps4EhgFrDrp06SpPHtcBrVZzcPat8MnDrMNjNpPDdx8Pj9eebLxxHNmWQxsBga/6iRJO27VA1VSbwNO0r+AHhZVX14THY4DvT29tZQ3zpIkhqSrKmq3k7HMVJJZtFYun5yVd0+oP0/A2dV1S61wJPcC3ypqj46oO0k4DYaXzhmpHMO5DlHkoY3kvPNmBWRqKrrgBVJ/jrJ4IfqSpI0kWwBdgAzBrXPADbtOhya7UON/2Vzvr2ZU5I0ysa0Cl9VrQM+BnwyiXVRJUkTUlX9AlgDLBjUtYBG5byh3DnM+NVVtX0v55QkjbIxL2PefA7U+cBlSV421vuXJGmMXAqcneRdSV6U5HIaS/GuBEjyhSRfGDD+SmB2ksua498FnM3ORSN2O6ckqf3aXYVvJ0kOoVGy/N8AJwOnJJlXVX1jGYckSe1WVdckOQy4gMZDb9cDpw94kPycQePvT3I68Fc0ypJvBP68qr4ygjklSW3W7udA7Qe8HHgN8Ls0yqweAPyIZ55lsbWdMUiS1ClVdQVwxTB9pwzRdhuw29UZu5tTktR+7b4C9UlgIfBbwDbgMuBvq+rHbd6vJEmSJI26dj8H6s8Akvw68Aoay/Y+kuQhGlefftTO/UuSJEm7s2JtHxev2sDGrduYNb2HJQvnsmi+tc40vDG5B6qqHgdubL5I8hzgHc0H7F5QVY/tbntJkiRptK1Y28fS5fewbfsOAPq2bmPp8nsATKI0rDGvwgdQVQ9V1V8CK4DrOhGDJEmSJreLV234VfLUb9v2HVy8akOHItJ4MKZV+IbwCeBfdTgGSZIkTUIbt24bUbsEHboCNcC7gdd2OAZJkiRNQrOm94yoXYIOJ1BVdUdV/bCTMUiSJGlyWrJwLj1Tp+zU1jN1CksWzu1QRBoPOr2ET5IkSeqI/kIRVuHTSJhASZIkadJaNH+2CZNGZEyX8CU5KcmvjeU+h4njvCT3J3k6yZokJ+5h/MnNcU8nuS/Jufs6pyRJkqTxZ6zvgfomMGeM97mTJG8FLgc+DswH7gBuSDJkXEmOAq5vjpsPXAh8Ksmb9nZOSZIkSePTWCdQGeP9DeV9wLKq+mxV/bCqzgceoVERcCjnAhur6vzm+M8Cnwfevw9zStKEt2JtHydcdCtHffDrnHDRraxY29fpkCRJ2meT6h6oJAcAxwGXDOq6CTh+mM1e1ewfaBXwx0mm0kgKRzonAI899hjLli3bqe0lL3kJL3/5y9m+fTtXX331LtvMmzePefPm8dRTT3Httdfu0t/b28sxxxzDE088wVe/+tVdD+ZVr2Lu3Lls2bKFr33ta7v0n3TSSTzvec9j06ZN3Hjjjbv0v+51r+M5z3kODz30ELfccssu/aeddhozZ87kvvvu4/bbb9+l//Wvfz2HH344GzZs4M4779yl/41vfCPTpk1j/fr1rF69epf+t7zlLRx44IGsW7eOdevW7dJ/1llnMXXqVO6++25+8IMf7NJ/9tlnA3DHHXdw77337tQ3depUzjrrLABuu+027r///p36DzzwQN7ylrcAcPPNN/Pwww/v1H/ooYdy5plnAnDjjTeyadOmnfoPO+wwzjjjDABWrlzJY489tlP/zJkzOe200wBYvnw5P/vZz3bqP/LIIzn11FMBuPbaa3nqqad26j/qqKM4+eSTAbj66qvZvn37Tv1HH300xx/f+EgO/tyBnz0/e6P72fv+E89i5aZD2F6N7836tm5j6fJ7uP3223nptH/Zafs9ffYkSeomnX4O1Fg7HJgCbB7UvhmYOcw2M4cZv39zvhHNmWRxktVJVv/85z8fWfSSNE7c8uhBv0qe+m3bvoNbHj2oQxFJkjQ6UlVjt7Pk/wIvrKp79zi4PfufBfQBJ1fV7QPa/zNwVlXtUvQ/yb3Al6rqowPaTgJuA2bRuAI1ojn79fb21lDfdEvSeHfUB7/OUGeXAPdf9Pstz5NkTVX1jlpgk5jnHEka3kjON5PtCtQWYAcwY1D7DGDTrsOh2T7U+F8259ubOSVpQps1vWdE7ZIkjReTKoGqql8Aa4AFg7oW0KicN5Q7hxm/uqq27+WckjShLVk4l56pU3Zq65k6hSULh70oL0nSuDCpikg0XQp8McldwHdoVNmbBVwJkOQLAFX1R83xVwJ/luQy4G+AE4Czgbe3OqckTTb9D6W8eNUGNm7dxqzpPSxZONeHVUqSxr2xTqD+G40lbx1TVdckOQy4AHg2sB44vaoeaA6ZM2j8/UlOB/6KRlnyjcCfV9VXRjCnJE06i+bPNmGSJE04Y5pAVdV/Gsv9DaeqrgCuGKbvlCHabgNetrdzSpIkSZoYJtU9UJIkSZK0L0ygJEmSJKlFJlCSJEmS1CITKEmSJElqUVsSqCSfSDJ9wPtfT3JRO/YlSZIkSWOlXVegFlTV1v43VfU48Ltt2pckSZIkjYl2JVD7JTmk/02SQ4GpbdqXJEmSJI2Jdj0H6nLgH5Jc03z/VhoPopUkSZKkcatdCdTVwF3Aa5rv31FVP2jTviRJkiRpTIz6Er4kAdZW1fqq+lTzZfIkSZo0kjwryaeSbEnyz0muS3JkC9udl+T+JE8nWZPkxEH930pSg15/174jkSQNNuoJVFUV8L0kLxntuSVJGicuA94EvB04ETgU+FqSKcNtkOStNJbAfxyYD9wB3JBkzqCh/wN49oDXvx316CVJw2rXEr6XAGuT3As8BYRGbvU7bdqfJEldIck04E+BP6mqbzTb3gk8AJwKrBpm0/cBy6rqs8335yc5DXg3sHTAuKeqalNbgpck7VG7Eqgz2jSvJEnd7jgalWdv6m+oqoeS/BA4niESqCQHNLe7ZFDXTc1tBnpbkrcBm4EbgI9U1c9HL3xJ0u60q4z50cBjVfUA8PvAh4CeNu1LkqRuMhPYAWwZ1L652TeUw4EpzTG72+Z/AmfRKNL0X2ksE/zKcIEkWZxkdZLVjz76aMsHIEkaXrsSqEuq6skkr6TxF/3NwOfatC9JktouyceGKOAw+HVKO2OoqquqalVV3VNVf0fjMSELkrxsN+N7q6r3iCOOaGdokjRptGsJX79FwJVVdW2SpXscLUlS97oM+NIexjwIvJLG1aTDgYGXfWYA3x5muy00rlrNGNQ+A9jd/U6rm9u9APjfe4hNkjQK2pVAbUzyReAkYF6SZ9E4mUiSNC5V1RZ2XZa3iyRrgO3AAhpL7miWMH8Rjcp6Q839i+Z2C4C/H9C1gN0s0QOOpXF+faSFQ5A0jq1Y28fFqzawces2Zk3vYcnCuSyaP7vTYU1K7Uqg3gwsBD5aVY8neTbw/jbtSx3kL7Mk7ayqnkjyOeATSX4CPAZcCnyfxpJ2AJL8CPh0VX262XQp8MUkdwHfAc4FZgFXNsf/Fo1l8dfTSOReDPx3YG1zvKQJasXaPpYuv4dt23cA0Ld1G0uX3wPgv7s6oC0JVFX9M7B8wPtHGOLbseazLY6sqiG/kVN385dZkob1XuCXwDU0iijdAvxRVe0YMGYujWV+AFTVNUkOAy6g8Xyn9cDpzYJMAL8AXge8BzgYeAj4Oo0qfAPnlTTBXLxqw6/+vdVv2/YdXLxqg//m6oC23gOV5O9pfEP2MNDX/HPgqwf4AI17pTTO+MssSUOrqn8Bzm++hhuTIdquAK4YZvxDwMmjFaOk8WPj1m0jald7tbuIxNuBzwC/BxxGo6T54TQfrNsc8402x6A28ZdZkiSp/WZN76FviH9fzZruU4I6oV1lzAGoql9W1TnAxcBnqmoGjatOvwWcArwD+JN2xqD2Ge6X1l9mSZKk0bNk4Vx6pu5cj61n6hSWLJzboYgmt7YmUP2q6nLg0SQfrqpfVNX9VfXtqrqmeX+UxiF/mSVJktpv0fzZXHjmscye3kOA2dN7uPDMY71lokPavYTvV6rquiQPJvlrYElVPTVW+1Z79P/SWoVPkiSpvRbNn+2/sbrEmCVQAFW1LsnHgE8m+S9V1TeW+9fo85dZkiRJk8mYLOEbqLlk73zgsiQvG+v9S5IkSdLeGtMrUEkOoVGy/N/QKMV6SpJ5XomSupsPTJYkSWpo93Og9gNeDrwG+F3geOAA4EfA+4BlVbW1nTFI2jc+MFmSJOkZ7b4C9UlgIY2y5duAy4C/raoft3m/kkaJD0yWJEl6RrufA/VnVfUCGg/RfRONh+d+JMlFSV7Yzn0PluRZST6VZEuSf05yXZIjW9juvCT3J3k6yZokJw7o+43mnD9Ksi3JQ0k+k+Sw9h6NNHZ8YLIkSdIzxuo5UI9X1Y1VtbSqzgL+GnjDGCcbl9FI4t4OnAgcCnwtyZThNkjyVuBy4OPAfOAO4IYkc5pDZgGzgQ8AxwL/GjgJ+F9tOgZpzPnAZEmSpGeMeRU+gKp6qKr+ElgBXNfu/SWZBvwpjedPfaOq/jfwTuClwKm72bT/Pq3PVtUPq+p84BHg3QBVtb6qzqyq66rqx1V1G7AEODXJoW09KKlpxdo+TrjoVo764Nc54aJbWbF2dGuy+MBkSZKkZ3QkgRrgE8DzxmA/xwFTgZv6G6rqIeCHNApb7CLJAc3tbhrUddNw2zQdCvwL4IOC1Xb9BR76tm6jeKbAw2gmUT79XJIk6RljWsZ8CO8GHh+D/cwEdgBbBrVvbvYN5XBgSnPM4G2GvGqVZDrwX4HPVtUvhxmzGFgMMGfOnKGGSC0bqwIPPjBZkiSpoaNXoKrqjqr64d5un+RjSWoPr1NGMeTdxXIwsBLoo3FP1JCq6qqq6q2q3iOOOGIsQtMEZoEHSZKksdXpK1D76jLgS3sY8yDwShpXkw4HHh3QNwP49jDbbaFx1WrGoPYZwKaBDc3k6frm29dX1dN7jFwaBbOm99A3RLJkgQdJkqT26PQ9UPukqrZU1Y/28HoKWANsBxb0b9ssYf4iGpX1hpr7F83tFgzqWjBwmySHADfSSNBOr6onR/MYpd2xwIMkSdLYGu9XoFpSVU8k+RzwiSQ/AR4DLgW+D9zcPy7Jj4BPV9Wnm02XAl9MchfwHeBcGqXLr2yOP4RGUYlDgUXAQUkOam7702YSJrVN/31JF6/awMat25g1vYclC+d6v5IkSVKbTIoEqum9wC+Ba4Ae4Bbgj6pq4B34c2ks8wOgqq5pPqfqAuDZwHoaV5keaA45jsbyQIB7B+3vNcC3RvkYpF1Y4EGSJGnsTJoEqqr+BTi/+RpuTIZouwK4Ypjx3wJ22UaSJEnSxDSu74GSJEmSpLFkAiVJkiRJLZo0S/gkSZImixVr+ywwJLWJCZQkSdIEsmJtH0uX38O27Y06WX1bt7F0+T0AJlHSKHAJnyRJ0gRy8aoNv0qe+m3bvoOLV23oUETSxGICJUmSNIFs3LptRO2SRsYESpIkaQKZNb1nRO2SRsYESpIkaQJZsnAuPVOn7NTWM3UKSxbO7VBE0sRiEQlJkqQJpL9QhFX4pPYwgZIkaZQleRZwCfB2oAe4BTivqh7ezTYnAe8HjgNmAX9SVcsGjQnwX4DFwK8D3wX+XVX9oA2HoXFs0fzZJkxSm7iET5Kk0XcZ8CYaCdSJwKHA15JM2c02BwPrgfcAw93t/wHgPwDnAy8HfgJ8I8khoxS3JGkPTKAkSRpFSaYBfwosqapvVNX/Bt4JvBQ4dbjtqur6qvqLqvoy8H+HmDfAe4GLquorVbUe+GPgEOAdbTgUSdIQTKAkSRpdxwFTgZv6G6rqIeCHwPH7MO9RwMxB824Dbt/HeSVJI+A9UB20Zs2aLUkeAA4HtnQ6HnXMNOCJTgcxDkzEn9N4OKZOx/ibHdz33poJ7GDXv9c3N/v2Zd7+eQbPO+TNLkkW07hfCuDJJHvzJNV2fwY8B2qsdfrvtclmvPy8Wz7fmEB1UFUdAZBkdVX1djoedUaSq6pq8Z5HTm4T8ec0Ho5pPMQ4VpJ8DPjQHoa9ZixiaVVVXQVctS9ztPsz4DlQY82/18bWRPx5m0BJnbey0wGMExPx5zQejmk8xDhWLgO+tIcxDwKvBKbQuLLy6IC+GcC392H/mwbM8+CgeTftOnzU+BnQRONnemxNuJ+3CZTUYVU14f5iaYeJ+HMaD8c0HmIcK1W1hRaWmiVZA2wHFgD/s9l2JPAi4I59COF+GonSAuDu5ry/RqPK35J9mHe3/AxoovEzPbYm4s/bBKo77NPyCklS96iqJ5J8DvhEkp8AjwGXAt8Hbu4fl+RHwKer6tPN9wcDz2927wfMSTIP+GlVPVhVleQy4C+a294LXAA8STNRG6c8B0oaV1JVnY5BkqQJZcCDdN/Bzg/SfWjAmAI+UlUfbr4/BfjmENN9vqrObo7pf5Duv2XnB+mub9exSJJ2ZgIlSZIkSS3yOVCSJEmS1CITKEmSJElqkQmUJEmSJLXIBEqSJEmSWmQCJUmSJEktMoGSJEmSpBaZQEmSJElSi0ygJEmSJKlFJlCSJEmS1CITKEmSJElqkQmUJEmSJLXIBEqSJEmSWmQCJUmSJEktMoGSJEmSpBaZQEmSJElSi0ygJEmSJKlFJlCSJEmS1CITKEmSJElqkQmUJEmSJLXIBEqSJEmSWmQCJUmSJEktMoEaJUmmJ1mdZF2S9UnO6XRMkiRJkkZXqqrTMUwISaYAz6qqp5IcBKwHeqvqsQ6HJkmSJGmU7N/pACaKqtoBPNV8+ywgzZckSZKkCWJSLOFLclKS65L0JakkZw8x5rwk9yd5OsmaJCfuxX6mJ/ke8DBwcVVtGYXwJUmSJHWJSZFAAQfTWFL3HmDb4M4kbwUuBz4OzAfuAG5IMmfAmP57mwa/ZvWPqaqtVfXbwFHAO5LMaO9hSZIkSRpLk+4eqCRPAn9WVcsGtH0X+H5VnTOg7Z+AL1fV0r3czxXArVX15X0MWZIkSVKXmPT3QCU5ADgOuGRQ103A8SOYZwbwVFX9PMk04CTgM0OMWwwsBjjooIOOe+ELX7i3oUvShLdmzZotVXVEp+MYz5KcAZxxyCGHnHP00Ud3OhxJ6kojOd9M+gQKOByYAmwe1L4ZOHUE8/wmcFWS/uIRn6qqewYPqqqrgKsAent7a/Xq1XsVtCRNBkke6HQM411VrQRW9vb2nuM5R5KGNpLzjQnUKKmqu4B5nY5DkiRJUvtMliISu7MF2AEMLvgwA9jUjh0mOSPJVU888UQ7ppckSZLUJpM+gaqqXwBrgAWDuhbQqMbXjn2urKrF06ZNa8f0kiRJktpkUizhS3Iw8Pzm2/2AOUnmAT+tqgeBS4EvJrkL+A5wLjALuLIT8UqSJEnqTpPlClQvsLb56gE+0vzvjwJU1TXAe4ELgHXAq4HTq6otNy+7hE+SJEkanyZFAlVV36qqDPE6e8CYK6rquVX1rKo6rqpub2M8LuGTJEmSxqFJkUBJkiRJ0mgwgeoAl/BJkiRJ45MJVAe4hE+SJEkan0ygJEmSJKlFJlCSJEm6gzsWAAAgAElEQVSS1CITqA7wHihJkiRpfDKB6gDvgZIkjRW/tJOk0WUCJUnSBOaXdpI0ukygJEmSJKlFJlCSJEmS1KL9Ox3AZJTkDOCM5z//+Z0ORZIkTUAr1vZx8aoNbNy6jVnTe1iycC6L5s/udFjaB/4/7R5egeoA16NLkqR2WbG2j6XL76Fv6zYK6Nu6jaXL72HF2r5Oh6a95P/T7mICJUmSNIFcvGoD27bv2Klt2/YdXLxqQ4ci0r7y/2l3MYGSJEmaQDZu3TaidnU//592FxMoSZKkCWTW9J4Rtav7+f+0u5hASZIkTSBLFs6lZ+qUndp6pk5hycK5HYpI+8r/p93FKnwdYBU+SZLULv2V2azYNnH4/7S7pKo6HcOk1dvbW6tXrx7xdpaxlDRZJFlTVb2djmMi2NtzjiRNBiM533gFapzpL2PZX4mlv4wlYBIlSZIktZn3QI0zlrGUJEmSOscEapyxjKUkSZLUOSZQ44xlLCVJkqTOMYEaZyxjKUmSJHWORSQ6YF/KmFvGUpIkafRY3VgjZQLVAVW1EljZ29t7zt5sv2j+bH+xJUmS9pHVjbU3XMInSdI4k2R6ktVJ1iVZn2SvvpCTJjurG2tveAVKkqTx5+fASVX1VJKDgPVJllfVY50OTBpPrG6svdHVCVSSI4GfVNUvOh2Lhua6YUkae1W1A3iq+fZZQJovSSMwa3oPfUMkS1Y31u503RK+JPOTfCTJ94AHgC1J/j7Jv04yvdPx6Rn964b7tm6jeGbd8Iq1fZ0OTZI6JslJSa5L0pekkpw9xJjzktyf5Okka5KcuBf7md48Vz4MXFxVW0YhfGlSsbqx9kZXJFBJXpTkk0keAG4BXgB8HPh14NXA94D3AJuT3JLk/M5Fq36uG5akIR0MrKdx3trlq+0kbwUup3Gemw/cAdyQZM6AMf33Ng1+zeofU1Vbq+q3gaOAdySZ0d7DkiaeRfNnc+GZxzJ7eg8BZk/v4cIzj3U1jXarW5bw/Q6NpQd/Cnyrqn45oO/7zdfHkswG3gD8AfCpMY9SO3HdsCTtqqquB64HSLJsiCHvA5ZV1Web789PchrwbmBpc455I9jf5uaVqBOBL+9D6NKkZHVjjVRXXIGqqs9X1flVdfOg5GnwuL6quqKqThvL+DS04dYHu25YkoaW5ADgOOCmQV03AcePYJ4ZSQ5p/vc04CRgl8v/SRY3q/WtfvTRR/c+cEnSr3RFAqXxyXXDkjRihwNTgM2D2jcDM0cwz28C325eefo28KmqumfwoKq6qqp6q6r3iCOO2NuYJUkDdMsSvpY014cfWVV3dDoWPfOAOavwSdLYqqq7gJaX+UmSRk9XJVBJ/h54MY2KQn3NPwe+eoAPAIs6FeNoSHIGcMbzn//8Toeyz1w3LEkjsgXYAQwu+DAD2NSOHU6kc44kdYNuW8L3dhrViF4CvBT4t8BngK8Ba5t9v9ax6EZJVa2sqsXTpk3rdCiSpDHUfK7hGmDBoK4FNM5x7din5xxJGkVddQWqWUDinCTvAZ6sqs81b7idDRwJzAJu72SMkiTtTpKDgf7LPfsBc5LMA35aVQ8ClwJfTHIX8B3gXBrntys7Ea8kaWS6KoHqV1WXJ/mDJB+uqg8D9zdfkgZZsbbP+9Ck7tILfHPA+480X58Hzq6qa5IcBlwAPJvGM6NOr6oHxjxSSdKIdWUCBVBV1yV5MMlfA0uq6qlOxyR1mxVr+1i6/J5fPdC4b+s2li5vFOIyiZI6o6q+RePZhrsbcwVwxZgEtA/8gkaSdtVt90DtpKrWAR8DPtl8iK6kAS5eteFXyVO/bdt3cPGqXR4HI2mSSnJGkqueeOKJEW3X/wVN39ZtFM98QbNibV97ApWkcaKrEyiAqnoEOB+4LMnLOh2P1E02bt02onZJk8/eFpHwCxpJGlpXJ1BJDknyTuB64E3AKq9ESc+YNb1nRO2S1Cq/oJGkoXVVApVkvySvSPLBJLcCj9K46XYG8D7gBVXl2gGpacnCufRMnbJTW8/UKSxZOLdDEUmaKPyCRpKG1lUJFPBJ4EvAx4FXAJcBR1fVi6vqsqra2tHopC6zaP5sLjzzWGZP7yHA7Ok9XHjmsd7kLelX9vYeKL+gkaShpao6HcMukvw6jQTqZGAO8BCwrKp+1NHARllvb2+tXr2602FIUtdKsqaqejsdx0SwN+ccq/BJmixGcr7pyjLmVfU4cGPzRZLnAO9oPmD3gqp6rJPxSZON/4iSJqdF82f7uy5Jg3TbEr4hVdVDVfWXwArguk7HsztJDkzyQJJLOh2LNBosZSxJkvSMcZFADfAJ4HmdDmIPPgT8Y6eDkEaLpYwlSZKe0ZVL+Hbj3cDjnQ5iOEleALwQWAkc0+FwpFFhKWNJkqRndMUVqCRzkhy6p3FVdUdV/TDJS0c4/0lJrkvSl6SSnD3EmPOS3J/k6SRrkpw4kn00XQIs3YvtpK5lKWNpfNvbKnySNF6sWNvHCRfdylEf/DonXHRr228z6IoECvh94NEkNyX5d82iEb/SfD7Ua5JcluR+4LYRzn8wsB54D7DL1+ZJ3gpcTqN8+nzgDuCGJHMGjFmXZP0Qr1nN/jcA91bVvSOMTepqljKWxreqWllVi6dNm9bpUCRp1HXiXu2uWMJXVZ9J8nXgD4BFwKVJ1gPXA8+lkWD9M42lcecCt45w/uubc5Fk2RBD3kejTPpnm+/PT3IajSWDS5tzzNvDbl4JvC3JH9JI2KYm+VlVfXQksUrdpr8Cl1X4JElSt9ndvdrt+rdKVyRQAFX1IPBp4NNJpgFnAL8H/B9gYVXd3Y79JjkAOI7G8ruBbgKOb3WeqlpKM9lqLhE8ZqjkKcliYDHAnDlzBndLXclSxpIkqRt14l7trkmgBqqqJ4AvNV/tdjgwBdg8qH0zcOpo76yqrgKugsZDDUd7fkmSJGmymDW9h74hkqV23qvdLfdAkeTSJCcm6ZqY9lZVLauq93c6DkmSJGki68S92t2UrPQAfwdsTrIsyaIkY1HmawuwA5gxqH0GsKkdO7QikiRprHjOkTSRLZo/mwvPPJbZ03sIMHt6DxeeeWxbbz1IVXetIkvyO8Abmq+jgFuAFcDKqnp0FOZ/Evizqlo2oO27wPeqavGAtnuBrzTvbWqL3t7eWr16dbuml6RxL8maqurtdBwTgeccSRreSM433XQFCoCququqPlRVxwC/TaNk+dnAw0n+Icn7k4wopUxycJJ5SebROOY5zff9VRwuBc5O8q4kL0pyOTALuHLUDkySJEnSuNd1CdRAVfXjqvrvVXUSMBv4f4BXA28f4VS9wNrmqwf4SPO/P9rczzXAe4ELgHXNfZxeVQ+MxnEM5nIKSZIkaXzquiV8k4nLKSRp91zCN3o850jS8Mb1Er6BksxsFpPoTfLsJGm275fkjzsdnyRJkqTJpasTqKraBMwDVgJ3Ak8leYjGMrszOhnbvnAJnyRJkjQ+deWDdAeqqg8nuZtGsvc1YCaN+6H+v44Gtg+qaiWwsre395xOxyJJkiSpdV19BapfVX0deBB4V1U9UlWrq+rxTsclSVK3c9WDJI2ucZFAAVTV94D/t1nG/KBOxyNJ0nhQVSuravG0adM6HYokTQjjJoFqKhrLDv8hyaxOB7O3/DZQkiRJGp+6OoFKsn+SVyX5QJKbgB8DrwT+M/BIZ6Pbe34bKEmSJI1P3V5E4lDgpcBcYBvwH4DlVfXTjkYlSZIkaVLq6gSqmSj9TfNFkn8F/EGSI4AvVNXmTsYnSZIkaXLpuiV8SZ43XF9V/aSqlgE3A383ZkFJkiRJEl2YQAH/lORtexgzHfjsWATTDhaRkCRJksanbkygArwnyYYkP0ryxSQLBg6oqm9W1f/sUHz7zCISkiaDFWv7OOGiWznqg1/nhItuZcXavk6HJEnSPuvGBApgDvAV4IvAwTSe//S3Sbo1XknSACvW9rF0+T30bd1GAX1bt7F0+T0mUZKkca9bE5J3VNVfVNV/q6o30qjE92rgP3Y4LklSCy5etYFt23fs1LZt+w4uXrWhQxFNXi4bl6TR1Y0J1BbgJwMbqurHwHuAd3UkIknSiGzcum1E7Wofl41L0ujqxgRqHbB4iPYHgNljHIskaS/Mmt4zonZJksaLbkygLgAWJ7k2ySlJfiPJbOA/Afd1OLZR4XIKSRPdkoVz6Zk6Zae2nqlTWLJwbocikiRpdHRdAlVVdwGvAA4DvgE8CjwIvAF4XwdDGzUup5A00S2aP5sLzzyW2dN7CDB7eg8Xnnksi+a7kECSNL7t3+kAhlJV64HXJTkMOA6YAny3qn7a2cgkSa1aNH+2CZMkacLpygSqX1U9BtzU6TgkSZIkCbpwCZ8kSZIkdSsTKEmSJElqkQmUJEmSJLXIBKoDLGMuSZIkjU8mUB1gGXNJkiRpfDKBkiRJkqQWmUBJkiRJUotMoCRJmsC871aSRpcJlCRJE5j33UrS6DKBkiRJkqQWmUBJkiRJUotMoCRJkiSpRSZQkiRJktQiEyhJkiRJapEJVAdYUlaSJEkan0ygOsCSspIkSdL4ZAIlSZIkSS0ygZIkSZKkFplASZIkSVKLTKAkSZIkqUUmUJIkSZLUIhMoSZIkSWqRCZQkSZIktcgESpIkSZJaZAIlSZIkSS0ygZIkSZKkFplAjaIk/yfJ95OsS/LNTscjSZq4khyY5IEkl3Q6FkmaTPbvdAAT0PFV9WSng5AkTXgfAv6x00FI0mTjFShJksaZJC8AXgjc0OlYJGmymRQJVJKTklyXpC9JJTl7iDHnJbk/ydNJ1iQ5cS92VcBtSe5OctY+By5JGlfG8HxzCbB0nwOWJI3YZFnCdzCwHvhC87WTJG8FLgfOA/6h+ecNSV5cVQ82x6xj6J/X71bVxuZ/v7qq+pI8G7g5yT1V9f3RPxxJUpdq+/kmyRuAe6vq3iTHt+k4JEnDmBQJVFVdD1wPkGTZEEPeByyrqs8235+f5DTg3TS/4auqeS3sp6/55yNJrgdeBphASdIkMUbnm1cCb0vyhzQStqlJflZVHx2FQ5Ak7cGkSKB2J8kBwHE0lkMMdBPQ8jd7SQ4C9quqnyc5GHgtcO0Q4xYDi5tvn0yyATgc2LIX4WtimAY80ekgxoGJ+HMaD8fU6Rh/s4P7HlWjdb6pqqU0k63mEsFjhkuehjnnjFS7PwOeAzXWOv332mQzXn7eLZ9vJn0CReMv7inA5kHtm4FTRzDPDOCrSWjO99mqunvwoKq6CrhqYFuS1VXVO5KgNXEkuaqqFu955OQ2EX9O4+GYxkOM48honW9aNtQ5Z6Ta/RnwHKix5t9rY2si/rxNoEZJVd0H/Han49C4tLLTAYwTE/HnNB6OaTzEOGlV1bIx2I2fAU00fqbH1oT7eU+KKnx7sAXYQeMK0kAzgE1jH44mm6qacH+xtMNE/DmNh2MaDzGOI+PyfONnQBONn+mxNRF/3pM+gaqqXwBrgAWDuhYAd4xRGPu0vEKS1P265HzTjTwHShpXJsUSvmZRh+c33+4HzEkyD/hps2zspcAXk9wFfAc4F5gFXDkW8TXXqEuSxrluP990I8+BksabVFWnY2i7JKcA3xyi6/NVdXZzzHnAB4Bn03iGx7+vqtvHKkZJ0vjn+UaSJr5JkUBJkiRJ0miY9PdASZIkSVKrTKAkSZIkqUUmUJIkSZLUIhMoSZIkSWqRCZQkSZIktcgESpIkSZJaZAIlSZIkSS0ygZIkSZKkFplASZIkSVKLTKAkSZIkqUUmUJIkSZLUIhMoSZIkSWqRCZQkSZIktcgESpIkSZJaZAIlSZIkSS0ygZIkSZKkFplASZIkSVKLTKAkSZIkqUUmUJIkSZLUIhOoUZbkq0keT/LlTsciSZrYPOdI0tgzgRp9lwN/1OkgJEmTguccSRpjJlCjrKq+Bfy803FIkiY+zzmSNPY6nkAlWZrk7iQ/S/JokpVJjmnDfk5Kcl2SviSV5Oxhxp2X5P4kTydZk+TE0Y5FktQZSf5dku83zzk/S3Jnkt9vw34850jSBNXxBAo4BbgCOB54LfBL4OYkvzHU4CTHJ3nWEO1HJXnubvZzMLAeeA+wbZi530pjOcTHgfnAHcANSeYMGLMuyfohXrP2eKSSpE57GPiPwMuAXuBWYEWSlw412HOOJGmwVFWnY9hJkoOBJ4BFVbVyUF+A1cBG4Myq2t5s/03gNuDaqvpAC/t4Evizqlo2qP27wPer6pwBbf8EfLmqlo7gGE5pzv/mVreRJHVGkp8CS6vqbwa1e86RJO1i/04HMIRDaFwZe3xwR1VVktNpnLiuSfIWYCaNbxDvBFo+4QyW5ADgOOCSQV030bg6NmqSnAGcccghh5xz9NFHj+bUkjShrFmzZktVHdGOuZNMAf6QxtWiOwb3e86RpMljJOebbkygLgfW0Tg57aKqNid5Lc1v/4BjmuPfWVU79mG/hwNTgM2D2jcDp7Y6SZKbgd8GDkryMPCHVbXTsTSvrK3s7e09Z/Xq1fsQsiRNbEkeaMOcx9I4x/wa8CTwxqq6Z6ixnnMkaXIYyfmmqxKoJJcCrwZevbsTU1VtTPI2GksrHgHOqqpfjlGYu1VVLZ/4JEkdsQGYB0wD3gx8PskpVbV+qMGecyRJA3VDEQkAkvwV8HbgtVV13x7GHg4sA1YB24Erk+zrsWwBdgAzBrXPADbt49ySpC5RVb+oqh9X1ZrmvUbrgH8/3HjPOZKkgboigUpyOc8kTz/aw9jfAL4B9AFvAF4DvA74m+YNv3ulqn4BrAEWDOpawBBr4yVJE8Z+wC6V9sBzjiRpVx1fwpfkr4F3AouAx5PMbHY9WVVPDhob4OvAYzSq9P0LcN+A9ekfAz40zH4OBp7ffLsfMCfJPOCnVfVgs/1S4ItJ7gK+A5wLzAKuHJWDlSR1VJKLaJxHHqJRtOgdNB6nscuzoDznSJKG0vEECjiv+ectg9o/Anx4YEOzItJ/Au6oqqcHtP9TktcAu6vJ3gt8c9D8HwE+D5zdnOeaJIcBFwDPpvEMj9OratRvYpYkdcRM4EvNP58Avg/8XlWtGjzQc44kaShd9xyoyaS3t7esiCRJw0uypqp6Ox3HROA5R5KGN5LzTVfcAyVJkiRJ44EJlCRJkiS1yARKkiRJklpkAiVJkiRJLTKBkiRJkqQWmUBJkiRJUotMoCRJkiSpRSZQkiRJktQiEyhJkiRJapEJlCRJkiS1yARKkiRJklpkAiVJkiRJLTKBkiRJkqQWmUBJkiRJUov273QAkiRp7yT5KnAKcEtVvbnD4aiLrFjbx8WrNrBx6zZmTe9hycK5LJo/u9NhSROCV6AkSRq/Lgf+qNNBqLusWNvH0uX30Ld1GwX0bd3G0uX3sGJtX6dDkyYEEyhJksapqvoW8PNOx6HucvGqDWzbvmOntm3bd3Dxqg0dikiaWEygJEmTRpKlSe5O8rMkjyZZmeSYNuznpCTXJelLUknOHmbceUnuT/J0kjVJThztWDT5bNy6bUTtk92KtX2ccNGtHPXBr3PCRbd6pU575D1QkqTJ5BTgCuBuIMBHgZuTvLiqfjp4cJLjgTVV9S+D2o8Cqqr+zzD7ORhYD3yh+dpFkrfSWIJ3HvAPzT9vaMbyYHPMOoY+V/9uVW3c/aFqspo1vYe+IZKlWdN7OhBNd+tf7th/xa5/uSPQdfeMeV9b9/AK1DjkNyWStHeqamFV/Y+qWl9V9wDvBI4AThg8NkmATwFfTjJ1QPtvAt+kkfAMt5/rq+ovqurLwP8dZtj7gGVV9dmq+mFVnQ88Arx7wDzzquqYIV4mTxrWkoVz6Zk6Zae2nqlTWLJwboci6l7jZbmj97V1FxOoccZfIEkaVYfQOBc+Prijqgo4HXgBcE2S/ZMcCdwK3Aks3dudJjkAOA64aVDXTcDxezvvMPs6I8lVTzzxxGhOqy62aP5sLjzzWGZP7yHA7Ok9XHjmsV6tGMJ4We44XhK9ycIlfOPM7n6B/ItR7eKyAU1glwPraCREu6iqzf9/e/ceLkld33n8/RE1O4YsLIJcBsfLgxIS0EFGowYQEYKLD3FEXUy8ZNYLKopEEzQEzSNqxASDokgQXXeQuCtglEDAgGjA5bIiE65iFAKCDHcRlGRgGPa7f1QdOfT0OafPTPfp7nPer+ep50xX/brqW9Xn9G++Vb/6VpK9gQuB04Cd2/ZvrKpHur2nR1sCmwB3dsy/E9in15UkOR94LvDrSW4FXltVj9mXqjoLOGvZsmVv24h4NWaW77rY7+kejMtwx3FJ9BYKr0CNGf+ANNe86qn5KsmxwO7Aq6dLhtrhcq8DXkVzb9Prq2rd3EQ5varap6q2qqonVdX2ncmTpOmNy3DHqRK6UUv0FgoTqDHjH5DmmsMGNB8l+RTwB8DeVXXjDG23BFYC5wIPAycm2dj+8x7gEWDrjvlbA3ds5Lol9WhchjuOS6K3UDiEb8wcvt+Oj6kWA/4BabC86qn5JslxwEHAS6vqX2douwXwLWA18EpgMc1wvs8nObi9T2rWqmptklXAvsDpkxbtC/z9hqxT0oYZh+GOE/E5nH40mECNGf+ANNfGZXy41Iskn6OpvLcc+HmSbdpFD1TVAx1tA5wN/AxY3pYyv3HSPVEfA46cYjubAju0Lx8HLEmyFLh3okQ5cCxwSpLLgIuBdwDbASf2ZWclzSvjkOgtFCZQY8g/IM0lr3pqnpkoPf7tjvlHAR+ePKOqKsmHgEuq6sFJ869P8lJguqtPy2hKnU9e/1HAycCKdj2nJnky8EFgW5rnRu1fVTfPcp8kSXPIBErStLzqqfmkqjLL9udPMX/amwCr6gKaB/XOtP4TaB7sK0kaEyZQkmbkVU9JkqSGVfgkSZIkqUcmUJIkSZLUIxMoSZIkSeqRCZQkSZIk9cgiEpIkqaszrlhtBU5J6rBBCVSS7YG7qmptn+ORJEkj4IwrVj/mGXCr71vDEV+/BsAkStKC1vMQviS7JjkqyVXAzcA9SU5P8oYkmw8uREmSNNeOOfdHj3mANsCahx/hmHOnfQSWJM170yZQSXZK8pkkN9M8tf1ZwMeB/wLsDlwFHAbcmeTbSQ4ddMCSJGnwbrtvzazmS9JCMdMQvhfQPEn9LcAFVbVu0rKr2+ljSRYDrwR+H/jsIAKVJElzZ7vNF7G6S7K03eaLhhCNJI2Oaa9AVdXJVXVoVZ3fkTx1tltdVSdU1cv7H6IkSZprh++3I4uesMlj5i16wiYcvt+OQ4pIkkbDbO6B2jvJFoMMRpIkjYbluy7m6AN3YfHmiwiwePNFHH3gLhaQkLTgzaYK39nAE5Osprn36cp2uqqqbhhEcBp9lrgdPj8DSYOyfNfFfp9IUofZJFCvAr4InAesBV4KHAlUkgeAa2gTq6r6Qr8D1eixxO3w+RlIkiTNrZ6H8NEUh3hbVb21qg6pqt2BPYAfAx8EvgVsA7y//2FqFFnidvj8DCRJkubWbK5AbQ885n9lVXVxkg8Af1RVr+5rZBp5lrgdPj8DSZKkuTWbK1D/B3hzl/nXAPv2JxyNk6lK2Vridu74GUiSJM2t2SRQ7wbeleR/JnlukscnWQS8F7h3MOFplFnidvj8DCRJkuZWz0P4qurHSV4AHA9cAayjScDW0jxoVwvMRJECK8ANj5+BJEnS3JrNPVBU1fXAfkm2B3YFClhVVbcPIjiNPkvcDp+fgSRJ0tyZNoFKsgS4r6p+MXl+Vd0K3Nql/XOq6ur+hihJkiRJo2Gme6BeAdyd5Lwk70ry1MkLkzwuyUuTfDrJTcCFA4tUkiRJkoZs2gSqqv4WeBZwJrAcuCHJqiQfTXIKcA/wZeCJwDuApww4XkmSJEkamhnvgaqqW2gKRxyfZDPgAOC/Aj8B9quq7w80QkmSJEkaEbMtInE/8HftJEmSJEkLyozPgUpybJI9kszmmVGSJEmSNO/0khQtAr4K3JlkZZLl7QN0JUmSJGlBmTGBqqp3VtVimop8q4GPAfckOTPJm5NsNeggx0WSbyT5eZKvDTsWSZIkSf3X87C8qrqsqo6sqp2B59KULF8B3JrkoiR/mmShP83zOOBNww5CkiRJ0mBs0H1NVXVDVf1NVe0JbA98Cdgd+IN+BjduquoC4JfDjkOSJEnSYGx0YYiquruqvlRVy6vqk/0Iqt+S7NkOOVydpJKs6NLmkCQ3JXmwfdbVHkMIVZIkSdIIm1UZ88mSnA9cBVwDXA38oKoe6ldgfbYpcC3NQ3+/3LkwyUE0w+8OAS5qf34zyW+1z8EiyZV0P16/V1W3DSpwSZIkSaNjgxMomqTpmcDDwEHA3kluAH4AXF1VH+tDfH1RVecA5wAkWdmlyfuAlVX1hfb1oUleDrwTOKJdx9I5CFWSpJ4k+QawF/DtqnrNkMORpAVjYxKofatql4kXSfYHXgj8I/CcjQ1sriR5IrAb0Dn88DzgxQPY3sHAwQBLlizp9+olSQvHcTT3IP/RsAORpIVkY+6Buj/JjhMv2qs8B7TV+r648aHNmS2BTYA7O+bfCWwzmxW1wxpPB/ZPcmuSF3W2qaqTqmpZVS3baisrwEuSNoyFiyRpOHpOoJI8s2PW24HTknwmyVuSHAOs62t0Y6aq9qmqrarqSVW1fVVdOuyYJElzy8JFkjS/zeYK1PVJXjfxoqp+ADyfpujC04HbgP37Gt3cuAd4BNi6Y/7WwB1zH44kacxNFC46DFjTuXBS4aKPA7sCl9AULloyqc2VSa7tMm03N7sgSZrKbO6BCnBYkqOAAr4PfLmqTgNOG0Rwc6Gq1iZZBexLM/xuwr7A3w8nKknSuLJwkSTNb7O9B2oJTVJxCs0Ztn9I8sUkG/08qUFKsmmSpUmW0uzzkvb1xNm+Y4EVSd6aZKckxwHbAScOK2ZJ0vwzqXDReR2LBla4KMnlSS6/++67+716SVqQZluF7w+r6sKJF0l2oKm69wHg6GMeWpUAABV+SURBVH4G1mfLgH+e9PqodjoZWFFVpyZ5MvBBYFuaoRf7V9XNcx6pJGk+m65w0T6zWVFbuOi5wK8nuRV4bee9t1V1EnASwLJly2pDg5YkPWo2CdQ9wF2TZ1TVDUkOA05ghBOotlJRZmhzAs1+SJI08qpqVgmXJKk/ZjP07kra5xd1uBlY3J9wJEma1yxcJEljbjYJ1AeBg5OclmSvJFskWQx8CLhxMOFJkjR/VNVaYKJw0WT70lTjkySNuJ6H8FXVZUl+h6b06rd4NPlaA7xmALFJkjR2kmwK7NC+/FXhIuDeqrqFpnDRKUkuAy4G3oGFiyRpbMyqiERVXQu8rC24sBvNjbDfq6p7BxHcfJXkAOCAHXbYYca2kqSxY+EiSZrHUmVRnmFZtmxZXX755cMOQ5JGVpJVVbVs2HGMs0kn7d52/fXXDzscSRpJs+lvRvr5TZIkaeNU1VlVdfBmm2027FAkaV4wgZIkSZKkHplASZIkSVKPTKAkSZIkqUcmUJIkSZLUIxMoSZIkSeqRCZQkSZIk9cgESpKkeSzJAUlOuv/++4cdiiTNCyZQkiTNYz4HSpL6ywRKkiRJknpkAjUEDqeQJEmSxpMJ1BA4nEKSJEkaTyZQkiRJktSjxw87AEnS/HTGFas55twfcdt9a9hu80Ucvt+OLN918bDDkiRpo5hASZL67owrVnPE169hzcOPALD6vjUc8fVrAEyiJEl9Ndcn7BzCJ0nqu2PO/dGvkqcJax5+hGPO/dGQIpIkzUcTJ+xW37eG4tETdmdcsXpg2zSBkiT13W33rZnVfA2OlV8lzWfDOGFnAiVJ6rvtNl80q/kaHCu/SprPhnHCzgRKktR3h++3I4uesMlj5i16wiYcvt+OQ4pIkjQfDeOEnQmUJKnvlu+6mKMP3IXFmy8iwOLNF3H0gbtYQEKS1FfDOGFnFT5J0kAs33WxCdMAJXkqcArwFGAd8NGqOn24UUnS3JroZ+ayCp8JlCRJ42kd8MdVdWWSbYBVSc6pqn8fdmCSNJfm+oSdCZQkSWOoqm4Hbm//fUeSe4AtABMoSRog74EaAkvKStLwJNk2yclJ7k7yYJLrkrykz9vYM8mZSVYnqSQrpmh3SJKb2jhWJdljA7e3G7BJVf10Y+KWJM3MBGoILCkrScORZHPgYiDAK4CdgEOBu6Zo/+Ikv9Zl/jOSPH2aTW0KXAscBnStpZvkIOA44OPArsAlwDeTLJnU5sok13aZtpvUZgvgy8DB08QjSeoTh/BJkhaS9wO3V9WbJs27qVvDJAE+C9yW5MCqerid/zTgn4HT2vWtp6rOAc5p26+cIpb3ASur6gvt60OTvBx4J3BEu56l0+1Mm9ydAXyiqi6Zrq0kqT+8AiVJWkiWA99LcmqSu9orPO9uk6XHqKoC9geeBZya5PFJtge+A1xKm+RsiCRPBHYDzutYdB7w4h7XEWAl8J2qOmWadg4bl6Q+MoGSJC0kzwQOAW4E9qMZQvcJ4F3dGlfVncDewC40V5y+A1wJvLGqHtmIOLYENgHu7Jh/J7BNj+v4XeAgYHmbCF6ZZJfORg4bl6T+cgifJGkheRxweVVNXD26IsmzaBKo47u9oapuS/I64HKaqnevr6p1cxLtNKrqIjwRKklzzi9eSdJCcjtwXce8HwJLurQFIMmWNEPlzgUeBk5MsrH95z3AI8DWHfO3Bu7YyHVLkgbIBEqStJBcDOzYMe/ZwM3dGrcV7r4FrAZeCbwUeBnw+W73TfWqqtYCq4B9OxbtS1ONT5I0okygJEkLyaeAFyY5MskOSV4LvAf4XGfDNkE6G/gZsLyqHqqqG2nuiXoF8LGpNpJk0yRLkyyl6WuXtK8nX+k6FliR5K1JdkpyHLAdcGKf9lWSNADeA9VnSZ4KnAI8BVgHfLSqTh9uVJIkgKr6fpLlNM9e+hBwS/vzhC5tK8mHgEuq6sFJ869P8lKgptnUMppS5xOOaqeTgRXtek5N8mTgg8C2NM+N2r+qul4NkySNBhOo/lsH/HFVXZlkG2BVknOq6t+HHZgkCarqbJorS720PX+K+T+a4X0X0Dysd6b1n0CX5E2SNLpMoPqsqm6nuUmZqrojyT3AFoAJlCRJkjTmRuYeqCTbJjk5yd1JHkxyXZKX9HkbeyY5M8nqJJVkxRTtDklyUxvHqiR7bOD2dgM2qaqfbkzckiRJkkbDSCRQSTanqYwUmhtzdwIOBe6aov2Lk/xal/nPSPL0aTa1Kc0Y88OANVOs+yCaByt+HNiVphrSNyff+Ns+rPDaLtN2k9psAXwZOHiaeCRJkiSNkVEZwvd+4PaqetOkeTd1a9hWRfoscFuSA6vq4Xb+02hu2D2tXd96quoc4Jy2/copYnkfsLKqvtC+PjTJy4F3Ake061k63c60yd0ZwCeqynK0kiRJ0jwxKgnUcuCfkpxK84yN24AvAp+rqsdUOWqrIu0PXAicmuS/AdsA3wEupU1yNkSSJwK7AZ/sWHQe8OIe1xGaBy5+p6pOmaLNAcABwC+SXA9sSfNQRS1MmwH3DzuIMTAfj9M47NOwY3zaELc9L3Tpc2Zr0L8D9oGaa8P+XltoxuV499zfpCM/GYokE+VhP0VzBWkpzVWmP6uq46d4z3Y0SdQ1wM7tz4Oqal2P23wAeHdVrexY52rgJVX13Unz/wJ4fVV1Pnyx23p3B74LXD1p9hur6ppp3nN5VS3rJW7NP0lOqiqHes5gPh6ncdincYhRgzXo3wH7QM01v9fm1nw83qNyBepxwOVVNXH16IokzwLeBXRNoKrqtiSvAy6nqXr3+l6Tp0GqqosYkXvLNDbOGnYAY2I+Hqdx2KdxiFGD5e+A5ht/p+fWvDveo/If/duB6zrm/RBY0qUtAEm2pBkqdy7wMHBiko3dn3uAR4CtO+ZvDdyxkeuWuqqqeffFMgjz8TiNwz6NQ4waLH8HNN/4Oz235uPxHpUE6mKgc3jcs4GuT2NvK9x9i2a43Stp7pt6GfD59h6kDVJVa4FVwL4di/alqcY3KCcNcN2SJI0y+0BJY2VUEqhPAS9McmSSHZK8FngP8LnOhm2CdDbwM2B5VT1UVTcCe9OUQP/YVBtJsmmSpUmW0uz7kvb15CtdxwIrkrw1yU5JjgO2A07s076up6rsPCRJC5J9oKRxMxJFJACSvILm2Us7ArfQ3Pv02c4qfG3bfYBLquo/OubvSFOo78dTbGMvmlLnnU6uqhWT2h1CUwp9W5rnRr13clEJSZIkSQvTyCRQkvojyTeAvYBvV9VrhhzOSJuPx2o+7pMkjTq/e+fWsI/3qAzhk9Q/xwFvmrGVYH4eq/m4T5I06vzunVtDPd4mUCMuyTeS/DzJ14Ydi8ZDVV0A/HLYcYyD+Xis5uM+aWGy/9M48bt3bg37eJtAjT7PaAxZkiOSfD/JL5LcneSsJDsPYDt7JjkzyeoklWTFFO0OSXJTkgeTrEqyR79j2VBJ3pXk6vZY/SLJpe39jf3ezlCOVfu7UEm6Pp9uI9c99p+/1Gf2fwuYfe/cmu/9d7+ZQI24YWfYApoxticAL6ap9rgOOL8tp7+eJC9O8mtd5j8jydOn2c6mNEVLDgPWTLHug2j+U/FxYFea8vrf7KgkOUy3Ah8AngcsA74DnJHkOd0aj9OxSvJC4GDg6hnajc0+SaPM/m/B2wv73rk0b/vvgagqpwFNwJ7AmTTPqypgRZc2hwA3AQ/SPINqjy5t9gK+Nuz9cfrV57EpzQOXD+iyLO3neBbwhEnznwb8BPjrHrfxwBS/L98DvtAx73rg6FH9nQHuBd4+zscK2Az4N5pnzl0AHD9Fu7HZJyenQU72f079nux7h3LMx77/HtTkFajBmjbLHosMW938Bs3V2593LqjmL3p/4FnAqUken2R7mjM5lwJHbOhGkzwR2A04r2PReTRn6EZKkk2SvI7m72C9B1GP2bE6ieZLuttjEH5lzPZJGiT7P/Wbfe8cmWf990CYQA1QVZ1TVX9eVV8D/l+XJu8DVlbVF6rqh1V1KHA78M45DVSzdRxwJc2XxHqq6k6a4Qa7AKfRfKFcCbyxqh7ZiO1uCWwC3Nkx/05gm4kXSc4HTgf2T3JrkhdtxDZnLckuSR4AHqJ5APWrquqabm3H4VgleRuwA/DBXjY8DvskDZr9nwbAvnfA5lv/PUiPn8uN6VGTMuxPdiwaqQxbj5XkWGB3YPfpviCq6rb27M3lNP8peH1VrZuLGKtqn7nYzjR+BCylGfb2GuDkJHtV1bXdGo/ysUrzcO6P03zeD89ivSO7T9Kw2f9ptux758y86b8HzStQwzMWGbYeleRTwB8Ae1fVjTO03RJYCZwLPAycmGRj/97uoRn/vXXH/K2BOzZy3X1TVWur6oaqWlVVR9CckXrvVO1H/Fi9iOZv9QdJ1iVZB7wEOKR9vd4NtDDy+yQNm/2fembfO3fmWf89UCZQI66q9qmqrarqSVW1fVV1vXStwUpyHI9+gf/rDG23AL5Fc/P0K2kKD7wM+HySbGgMVbWW5qbNfTsW7UuXMcoj5HHAVInGqB+rM2iGJyydNF0OfLX999rON4zBPkljwf5P9r1DN87990A5hG94xiLDFiT5HPBGYDnw8yQTZ0gfqKoHOtoGOBv4GbC8qh4CbkyyN3Ah8DHgyCm2synNvTbQfGktSbIUuLeqbmnnHwuckuQy4GLgHcB2NGOVhy7JJ2j2/6c0N/z+IU2VnPWeJTEOx6qq7gPu69j2v7fbWW9IwzjskzQC7P80I/veuTXf+u+BG0bpv4U40aVUI02ZxpM65v2YjjKNTkP/7GqK6cNTtN8HeFKX+TsCz55mO3tNsZ2VHe0OoSkT+hDNWZo9h32MJsW2Eri5je0u4Hxgv2naj92xYpoy5uO6T05Og5zs/5w2ZLLvnfPjPe/7735OaYPUAHRk2ZcAn6B5Lsa9VXVLW8b1FJpfkokM+y3Ab1fVzUMIWZKkjWb/J2k+M4EaoCR7Ad2eG3NyVa1o2xwCvB/YluaZGe+tqu/OVYySJPWb/Z+k+cwESpIkSZJ6ZBU+SZIkSeqRCZQkSZIk9cgESpIkSZJ6ZAIlSZIkST0ygZIkSZKkHplASZIkSVKPTKAkSZIkqUcmUFIXSfZKUkm2HMK2f5LkT+d6uxsjyQVJjh92HL1KsiLJA8OOQ5LAPme27HM0bCZQUneXANsCPxt2IJ2SPL3taJcNO5ZJDgSOGHYQ3bTH6jUds08FnjmMeCSpC/uc2bHP0VA9ftgBSKOoqtYCdww7jnFRVffO5faSPA5IVT2yIe+vqjXAmv5GJUkbxj5nduxzNGxegdKClWTPJP83yQNJ7k9yWZKd22XrDadI8uYktyT5jyRnJTkkSU1a/uEk1yZ5XZJ/S/LLJGd0rOP5Sc5Lck+SXyS5KMmLZhn6Te3P77cxXtCu+3FJPpTkp0keSnJNklfOcAwm9vNlSb7X7tvlSZ7X0e7Adn0Ptes/MkkmLX/McIq2/dVJ1iS5N8mFSbaetPyAJKuSPJjkpiR/meSJ08S5ov2c9k9yLbAW2Gmm45nkJ+0/T2/38yeT19exjbcnuSHJ2vbn26Y7dpI0G/Y59jkd27DPGWMmUFqQkjwe+AfgIuC5wO8Anwa6nl1qvyC/CHwOWAqcCRzVpenTgYOAVwG/B+wK/OWk5b8BnALsAbwAuBI4J8mTZxH+C9qfL6cZ8nFg+/ow4HDgA8AuwDeArydZ2sM6jwb+DHgezRCSr0x0Vkl2A04Hvt6u989ohk68u9uKkmwDfBU4GdgJ2JNmnyeW7wd8BTge+G3gzcBrgI/PEON/Aj4EvB34LeBmZj6ez29/vo3mWD2fLpK8qo3n08DOwHHACUkOmCEmSZqRfc567HPsc8ZbVTk5LbgJ2AIo4CVTLN+rXb5l+/p/A//U0eak5k/oV68/DDwIbDZp3pHADdPEEeB24A2T5v0E+NNp3vP0NrZlHfNXA3/RMe8C4O+mWdfEfu43ad7vtvO2b19/BfhOx/s+DNzasZ3j238/r33/06bY5neBD3XMWw48QDNEott7VrTr3G2Gz7Xb8SzgNV3W98Ck1xcDX+posxK4aNi/q05OTuM/2eest5/2OfY5Yz15BUoLUjXjp1cC5yY5O8n7kiyZ5i2/CVzWMe97XdrdXFX3T3p9G/CUiRdJnpLk80l+nOR+4Jft8q7bTnJiO4zggc7L/x3t/jOwHc2X8mQX0Zw5I8k3J63rBx3tru6ImUlx7zTFehe32+10FXA+cG2Sv0/yziRbTVq+G3Bkx379L+DXgW2m2kdgHc3Zvl+Z7fGcxlT7+FuzXI8krcc+xz6ng33OmLOIhBasqvrvST5NMyzh94G/TLK8qs7diNU+3LkZHjtU9mRga+C9NGf9HgK+DUw1FvsvgE9uRDwTMQC8FVg0RZyTX0+07+UES603o+qRJL8HvJBmSMlbgKOTvKSqrmrXexTNEI1Od0+zrYdq/Rt4Z3s8Z2u9/ZOkDWGf8xj2Od3Z54wJEygtaO2X61XAXyX5JvBHQLfO7F9ZfyzzC7q0m8nuwHuq6myA9ibXbaeJ7y7gro7Za9ufm0xq94skt9EMhfh2x/aua9us3oB4AX7Yrney3WmGU/xyirgLuBS4NMlHgB/QjNO/CvgX4Der6oYNjKczjpmO58NMOlZTmNjH/9Gx7uv6EKMkAfY5PbLP0cgzgdKClOQZNDeGnkkzjvuZwHOAv53iLZ8BLkpyOHAGzU2qr9qATf8YeEOS79EMH/hrHu2cenUXTTnU/doKPw+2QziOAT6S5HpgFfAGmhtdnzfVinr0NzTVlz5MM+zh+cCfAH/erXGSFwL70Pyn4E6am5qfyqMdw0eAf0xyM3AazTCJnYEXVNX7ZxlbL8fzJ8DLklxIc0bx513WcwxN1aRVwHk0Z4hfz6M3S0vSBrPPmRX7HI0874HSQvUfwLNpLun/mOay/FeAv+rWuKoupamq8x6asdvL27YPznK7bwY2pelsvgp8iebLtmdVta6N4600Y8f/oV30GZov5b8GrqXpbF/dnvHcYFX1L8BrgVe36/1EO031FPj7ac6s/SNwPU1n+NGq+rt2fecCrwBeSjPG/zKaKku3bEB4vRzPP2m39VPgim4rqaozgENphmVcR1Nd6pCqOmsDYpKkTvY5vW/PPkcjL81VT0mzleRTwD5VtcuwY5EkzW/2OdLocAif1KN2KMW3aEqf7gO8gymGFEiStDHsc6TR5RUoqUdJTqV5hsVmNE9m/zxwXPlHJEnqM/scaXSZQEmSJElSjywiIUmSJEk9MoGSJEmSpB6ZQEmSJElSj0ygJEmSJKlHJlCSJEmS1CMTKEmSJEnq0f8HQumiYBrRlT0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x504 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(3, 2, figsize=(12, 7))\n",
    "ax[0, 0].scatter(bns_data['snr'],\n",
    "                 bns_data['mchirp_src'][:, 0] - bns_data['mchirp_src'][:, 1])\n",
    "ax[0, 1].scatter(nsbh_data['snr'],\n",
    "                 nsbh_data['mchirp_src'][:, 0] - nsbh_data['mchirp_src'][:, 1])\n",
    "ax[1, 0].scatter(bns_data['snr'],\n",
    "                 bns_data['mchirp'][:, 1] / bns_data['mchirp'][:, 0])\n",
    "ax[1, 1].scatter(nsbh_data['snr'],\n",
    "                 nsbh_data['mchirp'][:, 1] / nsbh_data['mchirp'][:, 0])\n",
    "ax[2, 0].scatter(bns_data['snr'],\n",
    "                 bns_data['q'][:, 1] / bns_data['q'][:, 0])\n",
    "ax[2, 1].scatter(nsbh_data['snr'],\n",
    "                 nsbh_data['q'][:, 1] / nsbh_data['q'][:, 0])\n",
    "\n",
    "ylabels = [r'$\\mathcal{M}^{\\mathrm{src}}-\\mathcal{M}_{\\mathrm{inj}}^{\\mathrm{src}}$',\n",
    "           r'$\\sigma_{\\mathcal{M}}/\\langle\\mathcal{M}\\rangle$',\n",
    "           r'$\\sigma_{q}/\\langle q \\rangle$']\n",
    "ylims = [(-3e-2, 3e-2), (-1e-1, 5e-2),\n",
    "         (1e-5, 1e-3), (1e-4, 3e-3),\n",
    "         (6e-2, 2e-1), (5e-2, 3e-1)]\n",
    "for a, yl in zip(ax.ravel(), ylims):\n",
    "    a.set(xscale='log', ylim=yl)\n",
    "for a in ax[0, :].ravel():\n",
    "    a.axhline(0, ls='dashed', color='grey', zorder=0)\n",
    "for a in ax[0:2, :].ravel():\n",
    "    a.tick_params(which='both', labelbottom=False)\n",
    "for a in ax[1:, :].ravel():\n",
    "    a.set(yscale='log')\n",
    "for a in ax[2, :].ravel():\n",
    "    a.set(xlabel='signal-to-noise ratio')\n",
    "for a, l in zip(ax[:, 0].ravel(), ylabels):\n",
    "    a.set(ylabel=l)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}