From f3bc74d9ee0580a929f98df307649b0ce2a0c51c Mon Sep 17 00:00:00 2001 From: Daniel Williams Date: Wed, 13 Dec 2017 14:16:27 +0000 Subject: [PATCH] Added the up-to-date GRB beaming notebook. --- .../GRB Beaming estimation O1 + O2.ipynb | 706 ++++++++++++++++++ 1 file changed, 706 insertions(+) create mode 100644 notebooks/GRB Beaming estimation O1 + O2.ipynb diff --git a/notebooks/GRB Beaming estimation O1 + O2.ipynb b/notebooks/GRB Beaming estimation O1 + O2.ipynb new file mode 100644 index 0000000..0739aa3 --- /dev/null +++ b/notebooks/GRB Beaming estimation O1 + O2.ipynb @@ -0,0 +1,706 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import pymc3 as pm\n", + "import numpy as np\n", + "import theano\n", + "\n", + "from theano.compile.ops import as_op\n", + "\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.style.use(\"/home/daniel/papers/thesis/thesis-style.mpl\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import matplotlib\n", + "font = {'family' : 'normal',\n", + " 'weight' : 'bold',\n", + " 'size' : 16}\n", + "\n", + "matplotlib.rc('font', **font)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def background_rate_f(b, T, n):\n", + " \"\"\"\n", + " \n", + " \"\"\"\n", + " out = 0\n", + " #n = int(n)\n", + " for i in range(n+1):\n", + " out += ((b*T)**i * np.exp(- b*T)) / np.math.factorial(i)\n", + " return out\n", + "\n", + "def log_background_rate(b, T, n):\n", + " return np.log(background_rate_f(b, T, n))\n", + "\n", + "def signal_rate_part(s, n, b, T):\n", + " top_a = T * ((s + b) * T)**n \n", + " top_b = np.exp(-(s + b)*T)\n", + " p = (top_a * top_b) / np.math.factorial(n)\n", + " return theano.tensor.switch(theano.tensor.le(s, 0), 0, p)\n", + "\n", + "#@as_op(itypes=[T.dscalar, T.dscalar, T.dscalar, T.dscalar], otypes=[T.dscalar])\n", + "def log_signal_rate(s,n,b,T):\n", + " #if theano.tensor.lt(0, s): return np.array([[0.0]])\n", + " p = -log_background_rate(b,T,n) + np.log(signal_rate_part(s,n,b,T))\n", + " \n", + " return p" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def beaming_violins(traces, priors):\n", + "\n", + " width = 3.487 #* 2\n", + " height = width / 1.618\n", + "\n", + "\n", + " f, ax = plt.subplots(1,1, sharex=True, figsize = (width, height))\n", + " priors = [\"U(0,1)\", \"Jeffreys\", \"$\\delta(1)$\", \"$\\delta(0.5)$\"]\n", + " print \"|\\t\\t\\t| Lower\\t| MAP\\t| Median\\t| Upper\\t|\"\n", + " print \"|----------|\"\n", + " matplotlib.rcParams.update({'font.size': 10})\n", + " pos = [.5, 1, 1.5, 2]\n", + "\n", + " for i in range(len(priors)):\n", + "\n", + " o2_trace = traces[i]\n", + " #o2_trace = o2_traces[i]\n", + "\n", + " #i = i/2.0\n", + " t_data = o2_trace[2000:]['angle'][np.isfinite(o2_trace[2000:]['angle'])]\n", + " data = np.rad2deg(t_data)\n", + "\n", + " parts = ax.violinplot(data, [pos[i]], points=100, widths=0.3, vert= False,\n", + " #showmeans = True, showmedians=True, \n", + " showmeans=False, showextrema=False, showmedians=False)\n", + "\n", + " lower_p, medians, upper_p = np.percentile(data, [2.5, 50, 97.50])\n", + " lower, upper = pymc3.stats.hpd(t_data, alpha=0.05, transform=np.rad2deg)\n", + " hist = np.histogram(data, bins = 90)\n", + " MAP = hist[1][np.argmax(hist[0])]\n", + "\n", + "\n", + "\n", + " ax.hlines(pos[i], lower, upper, color='#333333', linestyle='-', lw=2, alpha = 0.5)\n", + "\n", + " ax.scatter( [lower, upper], [pos[i]]*2, marker='|', color='k', s=15, zorder=3)\n", + " ax.scatter( [MAP], pos[i], marker='D', color='k', s=15, zorder=3)\n", + " ax.scatter( [medians], pos[i], marker='s', color='k', s=15, zorder=3)\n", + " #ax2.vlines(inds, whiskersMin, whiskersMax, color='k', linestyle='-', lw=1)\n", + "\n", + " print \"| {} \\t| {:.2f}\\t| {:.2f}\\t| {:.2f}\\t| {:.2f}\\t|\".format(priors[i], lower, MAP, medians, upper)\n", + "\n", + " axis = ax\n", + " axis.set_yticks(pos)\n", + " axis.set_yticklabels(priors)\n", + " axis.set_xlim([0, 52])\n", + " ax.get_xaxis().set_minor_locator(matplotlib.ticker.AutoMinorLocator())\n", + " ax.get_yaxis().set_minor_locator(matplotlib.ticker.AutoMinorLocator())\n", + " ax.grid(b=True, axis='x', which='major', linewidth=0.5)\n", + " ax.grid(b=True, axis='y', which='major', linewidth=0)\n", + " #ax.grid(b=True, which='minor', linewidth=0.5)\n", + " ax.set_xlabel(r\"Beaming Angle [$\\theta$]\")\n", + " ax.set_ylabel(r\"Prior Distribution on efficiency\")\n", + " ax.tick_params(axis='y',which='both',left='off')\n", + " f.subplots_adjust(0.20, 0.15, .98, .95, wspace=0.05)\n", + " #f.savefig(\"O2a_beaming_posteriors_violin.pdf\")\n", + " #f.savefig(\"O2a_beaming_posteriors_violin.png\", dpi=300)\n", + " return f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Assuming VT is known perfectly\n", + "I start with a model in which VT is known perfectly, and where we do not account for the number of galaxies in the observed volume, and assume that sources can be distributed isotropically, in the same way that the other rates estimates do." + ] + }, + { + "cell_type": "code", + "execution_count": 224, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import theano.tensor as T\n", + "from pymc3 import DensityDist, Uniform, Normal\n", + "from pymc3 import Model\n", + "from pymc3 import distributions\n", + "\n", + "def grb_model_non_galaxy(number_events, background_rate, \n", + " observation_time, volume, grb_rate,\n", + " efficiency_prior = \"uniform\"):\n", + " with Model() as model:\n", + " signal_rate = pm.DensityDist('signal_rate', \n", + " logp=lambda value: log_signal_rate(value, number_events, background_rate, observation_time),\n", + " testval=50)\n", + "\n", + " cbc_rate = pm.Deterministic('cbc_rate', signal_rate / volume * 1e9) \n", + " \n", + " # Allow the efficiency prior to be switched-out\n", + " if efficiency_prior == \"uniform\":\n", + " efficiency = pm.Uniform('efficiency', 0,1)\n", + " elif efficiency_prior == \"jeffreys\":\n", + " efficiency = pm.Beta('efficiency', 0.5, 0.5)\n", + " elif isinstance(efficiency_prior, float):\n", + " efficiency = efficiency_prior\n", + " \n", + " def cosangle(cbc_rate, efficiency, grb_rate):\n", + " return T.switch((grb_rate >= cbc_rate*efficiency), -np.Inf, \n", + " (1.0 - ((grb_rate/(cbc_rate*efficiency)))))\n", + " \n", + " costheta = pm.Deterministic('cos_angle', cosangle(cbc_rate, efficiency, grb_rate)\n", + " \n", + " )\n", + "\n", + " angle = pm.Deterministic(\"angle\", theano.tensor.arccos(costheta))\n", + " \n", + " return model\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 225, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "priors = [\"uniform\", \"jeffreys\", 1.0, 0.5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the actual O1 + O2 scenario we had a total VT of 1502060.62853 Mpc^3yr with a sigma of 118341.996612 Mpc^3yr. This is based off a total observation time of 99.0 + 49.8 days (according to Tom Dents notebook for the rates estimate)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "O1_analysis_time = 49.8/365.25\n", + "O2_analysis_time = 99.0/365.25" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "total_time = O1_analysis_time + O2_analysis_time" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "total_vt = 1502060.62853 \n", + "total_v = total_vt / total_time" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3687013.74039\n" + ] + } + ], + "source": [ + "print total_v" + ] + }, + { + "cell_type": "code", + "execution_count": 229, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Actual O2\n", + "\n", + "number_events = 1 # There were no BNS detections in O1\n", + "background_rate = 0.01 # We take the FAR to be 1/100 yr\n", + "grb_rate = 10.0 \n", + "o2a_models = []\n", + "for prior in priors:\n", + " o2a_models.append( grb_model_non_galaxy(number_events, background_rate, total_time, total_v, grb_rate, prior))" + ] + }, + { + "cell_type": "code", + "execution_count": 316, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 1000000/1000000 [03:46<00:00, 4412.35it/s]\n", + "100%|██████████| 1000000/1000000 [03:42<00:00, 4496.96it/s]\n", + "100%|██████████| 1000000/1000000 [02:15<00:00, 7390.61it/s]\n", + "100%|██████████| 1000000/1000000 [02:13<00:00, 7477.37it/s]\n" + ] + } + ], + "source": [ + "samples = 1000000\n", + "o2a_traces = []\n", + "for model in o2a_models:\n", + " with model:\n", + " step = pm.Metropolis()\n", + " o2a_traces.append(pm.sample(samples, step))" + ] + }, + { + "cell_type": "code", + "execution_count": 317, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "|\t\t\t| Lower\t| MAP\t| Median\t| Upper\t|\n", + "|----------|\n", + "| U(0,1) \t| 3.64\t| 7.48\t| 11.83\t| 39.72\t|\n", + "| Jeffreys \t| 3.43\t| 7.49\t| 11.81\t| 51.13\t|\n", + "| $\\delta(1)$ \t| 3.40\t| 6.23\t| 7.67\t| 17.04\t|\n", + "| $\\delta(0.5)$ \t| 4.77\t| 9.01\t| 10.85\t| 24.20\t|\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAChCAYAAADju9JOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWlwG+eZ5/9vH7gIkjAoKRIj2yQYp2rKqVoLIr9lJ94I\n9MblJJOxQWqTzKRmszJk79Z+yCYSzM2HODWVjEgpOylvzUSEPFNbM0kqFGHncGJbIZgoE0/tbAhC\nsh3LiW02aR0UTREHL5x97IdmtwgSIA7iIvn+qiAA3W93P4T6ea9+3udPFEVRQKFQ9gRMvQ2gUCi1\ngzo8hbKHoA5PoewhqMNTKHsI6vAUyh6COjyFsofg6m1ApRgfH6+3CRRKQ3Hs2LFN23aNwwO5/8B6\nMjs7i/b29nqbsQlqV/E0ok1AYbvyNYC0S0+h7CGow1Moe4iCDr+0tFQLOygUSg0o6PAnTpzA1atX\na2HLjkWSJKRSKSSTScTjcaTT6XqbRKHkpOCk3cWLFwEAFy5cgCAIOH78OB566KGqG9boiKKIpaUl\nLC8vIx6PY+MaJI7jsLKyAoPBgObmZhiNxjpZSqHcpaDDP/XUU4hGo+jp6YHH44HdbseLL76Ixx9/\nvBb2NRyKoiAcDmNhYQGyLOctJ4oi4vE45ufnMT8/rzu+1WqFxWIBIaSGVlMoKgUd/ujRo3jyySez\ntk1NTVXNoEYmmUxidnYWyWSy5GPT6TTC4TDC4TA4jkNzczOam5thsVjAMHTulFIbCjp8d3c3lpaW\n0NLSgpmZGXR0dODUqVO1sK2hiMVimJub27JVLxZRFBGNRhGNRsGyLJqamtDc3IympiZw3K4KjaA0\nGAXvrsnJSRw5cgQAEAqF0NHRUW2bGo75+XksLCxU5dySJGFpaQlLS0sghMBkMundfrPZTFt/SkUp\neDdFo1HMzMxgaWkJv/vd7yp6cb/fX3RZQRAgCEJFr18IWZZx8+bNqjn7RhRFQSKRwJ07d/D+++/j\nnXfewczMDObm5hCLxZBMJjdNDlIopVCwhXe73Thz5gyi0SgGBwfLvpAgCAgEArDb7XC73fD5fPB4\nPIjFYvD5fHA4HHA4HHA6nXr5vr4+dHd3w+v1wuFw6MfUgng8jtu3byOVShV9jCzLeOWVV/Doo49W\npGWWZRnxeBzxeFzfxjAMeJ6H0WiE0WiEwWAAz/MwGAx0OEApSME7pLOzE+fPnwcAzMzMbOtimrMD\n6pgYgO7ENpsNXq9Xd3hAjQe22WxZxwuCAIfDsS078vHwww8DAEZGRhCJREpqTWVZxjPPPIPx8XH8\n5je/wZkzZ/R93/jGN7J+u46ODnzzm98sy0ZZlpFKpXJWRFplwPM8OI4Dx3FgWVZ/116SJEFRlLKe\nFGi/0eXLl8uyn1JfCjr82bNnIQgCFEXB5OQkJiYmyrqQ3W7H2NgYXC4XbDab7sgTExM4ffo0AGzq\nsmsxAN3d3XA6nXA4HAgEAhVv5UVRxOrqKtLpNObm5nDhwoWSjpdlGS+88ALeeustSJKEV199FTdu\n3MBjjz0Gi8WCYDCY5fALCwt48cUXK/o3lEIikdDnBwghOV8Acn6/efMmDh48iMXFxU1l1pdbvy/X\nZ431nyVJgiRJOfet30YfaZZPUY/ltFn56enpsi4SCATgcDjQ29uLYDAIh8MBu92u74/FYrDZbHqr\nDwAOh0N37JMnT2J4eFh3+Eqg3VyiKCKdTmN1dVVv+UrlzTffxNtvv63frJIk4e2330ZHRwd6enoq\nYm81yPfEYSvHVxQFkiRheXk5a782hCnk9Pn2Aeowav15N9q08XOuCmSr74UqnFyIoohMJpNzX75j\ntd8pH6VWWLnKy7Jc1hOjgg5//vx5jI6OorW1FVeuXMGlS5dKvsjY2BgGBwcRCATgcrkA3G3Ne3p6\nEIlEslp9QO3q9/f3w2azIRKJAAAikUjFuvNa99ZgMMBiscBms8FsNuPw4cP4/Oc/X9LY/XOf+xwS\niQR++ctfIpPJgOd5PPLII/jqV7+Kffv24ac//WlWC79v376qBy5pXfmN3XmWZbGwsID29nYwDAOW\nZcEwTNZrK15++WUAwOHDhytuczwez7oHGgGO48DzfL3N2EQx/1e5KOjwx48f1500GAyWbtnaObQZ\n+Y3/oR6PBz6fDzabDQMDAxAEAX6/H263G8FgEIIg6JOFoVAoa4xfabRxqSzLWFhYQDgcLqrFZxhG\nH7OPj4/j2LFjOHPmjN5j2fgosxKPNhmGgcFgyHpxHKe/b3UzpFIpNDc3l3VdOnbf2RR0eKfTicnJ\nSXR3d6Orq6usizidzk2Oqjm+zWbTx/Aa2veNrXklW/itYBgGBw4cQHNzM27evJm3S7fxmDNnzuSc\npS93gk6D53mYTCb9ZTQawfM8HctSSqZgn0DLnNHS0lLR5+Aej6ek8bggCHpPo1aYzWZ0dHTAZDIV\nVZ5hGDz22GPbfiRnMpnQ1taGw4cP46Mf/SgeeOAB3Hvvvdi/fz+am5thMBios1PKomALrygK7HY7\nlpaWMDk5iU9+8pMVu3gpDlyLlj0XPM+jo6MDt27dwvLyclWuoYXXWq1WNDU1NeSYkbI7KOjwLpcL\nXq8XhJBtBd7sZBiGweHDh4taJVcsPM/DarXSBTSUmpL3Lnv++ecBqOGvDocDnZ2d8Hq9NTOs0SCE\nYN++fejs7ITFYinrHAaDAW1tbejo6MADDzyAQ4cOwWq1Umen1Iy8LfzRo0cBqBNu3d3dAMqfpd9N\nGI1G3H///VhcXMTS0hJWV1fzzuQTQmCxWGC1WmG1WoueC6BQqkVeh9dWyLW1tYEQgpaWlrJn6Xcb\nhBA9bkCL0stkMlAUBbIsg+M4mM1mtLS0VOV5NYVSLgXH8MFgUE9ptVeXx24Fx3FobW3NuW995CCF\n0gjUdXkshUKpLTVbHkuhUOpPQYcnhOiz87RLT6HsbPJ26c+dO4fp6WmcPn0a4+PjCAQCGBsbq6Vt\nFAqlwuRt4Z1OJzo7O3Hy5EldpHG7CTAoFEp9yevwFy9exMWLFzE9PQ2/37/tBBgUCqX+5HX4vr4+\nHDt2TF/uCZSfAINCoTQGeR0+FArhypUrmJiYwJUrV6AoCoLBIEZGRmppH4VCqSB5Hf7UqVNYXFzE\nkSNH9NDaeq1Yo1AolWHLwJvW1lY4HA5MTk6CEKLH11MolJ1J3RJgUCiU2lPXBBiNjCwriEVSSCYk\nKIr6O3Acg6ZmDpYmHixLM85Qdh40AUYOUkkJC3eSyKSzE12IooRkUkJkIYUmK4eWVgMMRrZOVlIo\npVOU8owmCLEXWFpMIxpOYatktYoCrCyLWFkWYWni0GIzwGSijk9pfKgY2TpWVzKILBSfjx4A4qsi\n4qsijCYWLa08LE0cTTBJaVhKyq20m0NrUykJ4TvJ8o9PSrjzQRK3rq8iFkkhk9l+3jsKpdIUpS03\nNTUFALs2tFYUZdyZS6ACuSkhigpi0TRi0TSWV0RYm9IwmznwBpq3jlJ/aqItVyqa8sx6tEeClQ7+\nEUUZ83MJiGLlddczaWVtiJACxxGYzCyMRhYGEwuDgaFdf0rNKdjsnD9/Hk8//TSeeeYZPPXUUxU3\nQBAE+Hw+XYrK5/PB7XYjFoshEAhgaGgIACoqJKmRTEq4fSuOdCp/0y7LMn76s4vbTk0tigpWlkWE\nF1K4fTOO69MruH1zFQvzSSzF0kjERaTTEmS58hUPhaJRlLZcb2+vHktfDXLpxttsNjgcjqw1+JXU\nh//3H/8EMqKMH/3g5bxlZFnGV776JMYCv8Cvfn0Jf/udC2AYBgP/879jeuY9vVxnx0fwN9/+3yVd\nX1GAVEpGKkdlw7AE7NqL4wgYRv3MMAQMS9bUWgHCEDCEwNWrxkZcvvxr2mugbElBh7fZbDhx4gQW\nFxerkpc+n258LrajD68oCpIJCcmEhERCRCYjY/7OHEZG/0/O8rIs4+LFf8abv38dkiTiFy//GNev\nC+jv/0v8v4nXMD09pZfNdx5Nh73a3LhxCwf2H8T7wgoAgGGgVg5rL8IADCEAUTMYLcUkmI0pEKJW\nGgQAYdbLKd89d9Y29RR6gayqRd8HEP3LZnLVR9omWVIgSaX1cLZTvxVzrKIoZUmIV5ty7Sro8IuL\ni/pz+BdffLF0y7agkG78RrbTrSeEwGzhYLZwuAdG8DyT754EALz+RghvXXsTkiQCACRJxFvX3sTr\nb4TKun5ptqoOm+ud6I6rOhYhAMcRGE0s7ndYi2rh0yKLe9qMVf87SkXr2TQSGzXsG4Vy7Sro8NoM\nPQCEw+GSL7AVW+nG56Ki+vAcwYcPt+MLn//PkHJM2PU98SUk4gm88upPkcmkwfMGPPqpP8O3//o5\nfOEvHstq4Q/sP4jjfX+V0958FRghAG9gYDCoE3i8gQHLEXAcA4Yp7T9y5OKP1s7ZeDcmpbEoKmtt\nf39/VUJrt9KNB9QeQCgU0sftldSH13TOtVn6jRN3DMPgb79zAQAwFvgFel2P6WP4zo6PZJXd+D0X\nWbP0RhYGY+Vm6almO6VY6hpau5VuPKBKSq8fr1dDH57jGBw4aMbcrfimR3Oa07/0cz8+82m3rgFX\n7AQdbyCwtxlhsrAwGGjoLaX+NFw0SD7d+Grqw3Mcg/0fMiOXpiPDMPizz/YXLfjIcQSt9xjw4fua\n0LZfjbOnzk5pFBoylj6XY1c7247RxKJtvwl3PigvvNZoZNBiM9BYekpD05AOXy+arDwkSSm4Wm49\ndLUcZSdBF89soKXVgIPtFvB8/p+GEKDJyqH9sAUHDpqps1N2DHTxTA6MJhaHDlsQi6SQSEhqgIMC\nmvGGsuNpyMUzjQDDENj3meptBoVSUQo6/Pnz5zE6OorW1lZcuXIFly5dqoVdFAqlCjTE4hkKhVIb\nCk7aaYtntGg7CoWyc6nr4hkKhVJbCrbw1Vw8Q6FQaktdF89QKJTakreFf/755wGo+eUcDgc6Ozur\nkgCDQqHUjrwtvCYc6XQ6dfVYOktPoexs8rbwR44cAQC0tbWhtbUV0WiUztJTKDucvC389PQ0/H4/\ngsEgenp6oCgKIpHInhGTpFB2I3kdvrOzE263Gw6HA06nE3a7Ha2trbW0jUKhVJgtZ+k7OzvR2dmp\nf7969SoeeuihqhvVqIiyjFg8A1FWkJZkSLICnmVg4hiYeBYWA6tmh6VQGpSCj+X6+/vR1taGcDiM\n6enpPbFabj2KoiCWyOCD5RTCqylIW+hRMASwGjk0mzi0mniIVFSC0mAUdPgLFy7oXfnx8fGqG9RI\nxNMi3plfwVJSLKq8rABLSRFLSRG3kEQ0soowYrCZedjMPFpMPNgSM9JSKJWkKIdfH2F37NixqhrU\nCCiKgpuxBN6PxLGdRlqBguWkiOWkiBvRBBgCNBs5tJh5NBs5tJp58GzDpRWk7GKKylr75JNP7pkJ\nu5Qo4Q9zy1gsslUvBVkBFpNi1rnNPINmE68OBYwcmgwsOFoJUKpEQYd/4okn9EUzjz/+eNUNqifR\neBp/+GAZmRLljrZDIiMjkUlhfjmlbzNyDKxGDua1iUDtnfYGKNulqBRXNpsNdrsd586dw9e+9rW8\nZQVBgNfrRU9PD06fPr1p/9DQEJxOJ4aHh9HT06M/9qs3sqLgejSOG9FE0ckrq0lKlJES05u2cyyB\niWNg5lmYeBZGjoGR094ZWiFQClJUiist2OZXv/rVlmU1nbhcTqypwrpcLoyNjeWsEOpBeDUNYWEF\niUzxctCyLGP85z/GsU//edH56iuBKClYkSSspKSc+xkC8CwDA8fAwDLgGAKeZcCxBDxDwK1tW0lJ\niKdFsAwByxBwNfwbKPUlp8NrUXYAMDExgVAoVHKknSAIuqCEy+WCIAiYmJjQ5aM0CSlNDtrr9WaV\nDwQCmJycxPDwMPr6+vDFL34RVqsVoVBo25WFrCj40088jJQo4Tv/9JPSjpVlfOtr/xWvBV7Bv10e\nw9fP/X1JTn/26/8DN2buLjm+t6MLp771v0qyIa9titY72LryikTimJNiWdtYhoAl6jvDEHBkTX2W\nqJUCQ1QFWvW1pk4LVc+OIdr7+s9337XYBGadACZZ264pzjKEQF5TRKUh3KXx8MMPAyhOciynw8di\nMbjdbtjt9rKkmQHVgQcGBhCJRCAIArq7uxEKheByuTA6Ogqn0wmHw4GJiQkMDg6ir68vq7zH40Ff\nXx8AYGBgAIFAAOFwGMePHy/LHgBYSYmYCa9iMZnBakpE5M4HeHX0n4s+XpFl/Gz0h/jj71+HJEn4\n9cs/wez1aXy27wsgOZw+l1z067/7V9yYuSuYGZmfK8mGSlArGetSWW+X5vOaPLVWBehVgVZ5ZG1E\ndrkcstZ55apzSWADiMdXYbE0bZLGznmODd/y1Vtk04e8p817ldXVFTQ1WQEAN27ewoGDh3Arligo\ngZ3T4bWFM6UwNDQEt9utd90B6GG567dtpK2tLW/5np4e+Hw+eDweXXPO6/VidHS0ZPsAdTLM3mTQ\nJZdL5dobV/HutTchSWqXWpIkvHvtTVx74yoefKgyIpe7gbvOutaaZ20nWc683tmUjDo/kc/Z727T\nXSpbyz5r392t2ZVELltzOzsAyBwDSw6psLy3T56KY719RZ8rq0x2KYljYOFZtXfEEPAsQVuToeB5\n8o7hn3/+eZw4cQJnz55FJBKBoiiYnp7GyMhIzvJutxuCICAcDsPlcsHhcMDn8+kteSgUwsTEBARB\nQDAYhCAI+rZYLIbBwcGs8jabDW63G8PDwwDUdflOpxO9vb1F/Dy54VkG7a1mtLea0WriYW5vx6f/\n05cgFjkr/8gTX0Q8kcDlV1+CmEmD4w34xKc+g6/89XdydusjkTDs9rasbS//xJ/VwtsPHMSn+v6y\n7L+pVBgCxKIRfGj/vrUxvDqu17rt2rie1bvyd7vqLEOyuunaOyHZ3f1yu+Szs7Nob2+v8F+8PRrR\nJiDbrtERVS7cxBcWRCGKknte+sqVKzhy5AjGx8ez1sPXMvAmFotBEISiJKLHx8fLsk2UZMxE4ri9\nlCxqhn79GP7jrke3HMPncvhKj+EJAQzs3Vl67d3AahN2DHiWgFubnGMZsiNu4kahEW0CCtuVzx/y\ntvBat97v9+sH9vT0bNfOohkaGoLD4YDb7a7qdTiWwUf2W3Gg2Yi355YLTngxDIOvn/v7smfpy3Fu\njiWw8Kz+OM7EMzCtexxHJ7koxVLwsZyiKBgYGIDdbkcwGMzbpa80tX5s12Li4bzXhj98sIxoPLNl\nWYZh0PvZJypuA8cSNPEsmtYi7iwGjgbcUCpKQYfv7e3V5ZsbIUimmvAsg48dasGNWALXtxlHX4j1\nK+uajWpsvZnqyFOqTEGH7+rqQmtrK2ZmZnDPPffUwqa6QgjBffdYYLcY8M78ct4gl1JhGbUX0Wrm\n0Wri0Wzi6Np5Ss2hKa7yYDVyOHLYhhuxBG5EE5BKbO7VlXEs7rdbYDNTB6c0BkWnuGpra0NLS0st\nbas7Wmv/4VYzwqtpzC8nEU1kcs7mGzkGTQYWLWYeLUYOzSYeH8yl0W631N5wCiUPBVNcvfDCC+jq\n6oLX68UjjzyCEydO1Mq2hoFlCA40G3Gg2QhFUZCRFGRkGaKkqI/EeIa23pQdQcEx/JEjRzAyMoLB\nwUEIglCo+K6HEAIDR2AorNJFoTQcRd21bW1tIITgypUr1baHQqFUkYJSU6FQCAsLCzh//vyeS2BJ\noew2CkpNHT16FN3d3VAUhUpNUSg7nC1Da5966ikEg0EQQtDd3Y3vfe97tbSNQqFUmLwOf+7cOXi9\nXl2IYnp6umCKKwqF0tjkHcNvVJ3Z+J1Coew88jp8rjDavRBaS6HsZvJ26T0eD7q6uvTvWgKMd999\ntyaGUSiUypPX4UdHRzeluqLP4SmUnU3eLn2uvHbl5LqjUCiNQ8HQWkpuFFmGkoyvfSMAw4DwBhCW\nrmmnNC7U4UtAjq9CujMHeWUZ8uoKoGxOh0UMRhBzExhLE5TVOJT9+0F4vg7WUiiboQ5fBHJ8FeKt\n65CiCyiU6VJJp6CkU5AXI0AkguRSGMRkBmNtAdPSCqbZBsZkqpHlFEo2BR3+6aef1iPslpaW9tSa\neEWWkbkuQPpgdnvnSSYgJROQFj4AABCDAYy1FUyzWgkQcxNNREmpCQ2bxLLeyIk40u+9DSW+WvFz\nK+k0pMgdSJE76gaWBWOxgmluBdNkBWmygjHSXgCl8hSVxLK3txeKouz6JJYa4p05ZGbeA+TiBSa3\nhSRBXl6EvLyobyK8AcSizgXo70YznRSkbIuCDm+z2XDixAksLi7C6/VW7MJ+v7/onPNa4o1qVziK\nJCIz/R6k8HxVr1OULZk0lMU05MVo1nZ1UtACYjSBMZlBjCb1ZTDSyUFKQQomwFhcXMTFixdx6dKl\nLTXiNhIKhdDV1YVQKARAFZbQKgyfz6fr0A0NDcHv9+vlAFVxJhAIYGhoCIDq6JqybLWQV1eQ+v2V\nopxdlmWMXBqDXKsewDrUCcEopPnbyFwXkH73GlK/DyEZ+r9ITLyG5OsTSF27ivS715CeeQ+ZW+9D\nnL8NKbIAaSkGOb4KJZNBHsEhyi6nYAs/NXVXFikcDhd9Yk0jTpOJ8ng8iEQiAO5qxa8XivR6vXpZ\nm80Gh8OhS0kDgN1uhyAIFWvlNYndX738c4izN9XxdBFOIMsy/suz38JL//IaXv3Xf8M/PPv1ktRn\n/tvfnMW712/o3x+471783cCpku3PYxyUZAJKMrF1uUgEybn3AZYFYTmA50FYDoTjAJZThw3aPpZV\nYwwYVt3HMGrObYYDYRmAYQCG3XWTjtr98cMf/rC+hlSYgg7vdrvR398PQggGBwcrclFNCXZiYkJX\nmCmUL09r5cuVr9ZQRBFKIg4lk8bt2dv43plvF32srCj4p5/8DFff/iNEScIL47/G1OwsvvS5z+ZM\nYplLlvm3V17He+scfi4cwT/+4tXy/6AyqI5c9Jro+/oXsKb+uibbStb0Yjd+Xjs2Ho+jqakp69is\nd6w/r/7PWpGN2ssk9+dc3zfJyhLM3rqF9kOHoCQTkJOJu9fSTcp3/o02b961qQzZ4u+oMAXVY/1+\nv96qer3ebc/SC4IAu92uf4/FYrDZbAWHC5Xo1ivpFKTwHciJOJR0GkqOwJmtmHzrGt7447sQ1+Si\nRUnCG398F5NvXUPPxx7clm07H0XtIeXoJBU3eCBAMgEZSn6H16WYSbbzEEDZJBKfu2LIuS2Xk8kS\nlHQKiNyBxDGby2VrYOc8X0Hn3VAxbjrPxq/ryijRMCQjn2Pf1hRMceV0OrPUY0tBa8kBZHXHtda8\np6cHkUgENpstq2wuIpHItrvzxGAEd+gwAIBpsuLDZgue/PKXId2ZK2pG/q8efQTJ1Th+/OvLSGdE\nGHgOf/4fPoG/O/WVPHLRkazKDQB+9POXs1r4g212fPmxT23r7yqVXHZlwTBq955Z69pveL/bxV97\nZ5i72whRy6zdzIQwqioHYe5uY5i17wAIo35HYym1XnzlEgCAtN8HvkFsWg9Ji2Db9pd8XFXVYwcH\nB3UVWJvNtkn22ePxwOfzwWazYWBgAIIgYHh4GIODgwgEAgiFQnpFEQqFipKNLpbLly/rn5UP34/M\nrfcLBtgwDIN/ePbrAICX/uU1fOZPP17yGP6B++7d8ntFIYw6c8/zIJwBhOdAOB4wWMC3t6ufubVt\n+hh+943Hy0G7P2Zntxd01WhUNfDG4XDkVIHVWnObzbZp//HjxwGolcH68XolWvh8EJ6HoeMjkFps\nSAvvAJKYt6zm9KNj4+jrPVayXHTFJugIc/eRnNGoPpYzmkAMBvUZPm9QHTjXofwsuA81XqtFqT51\nUY/1eDwIBAL6eQshCELRZbcDa98Hk7UZ6ffehry8lLccwzA4/h97q26PahR3N/jGZFafwZvMqoPT\nlphSIls6/NWrV9Hb26vHzz/xROU00fM5cK5uey0j/IjBCMOf/DuIt96HOHujqEd1FYPjwTQ1g7Fa\nwVhoiC2l8uR1+LNnz2JqagrRaBQDAwN46KGHamlXXSGEgD/cAabFhszUH6Ck09W4CIjZosbPW5vV\n1XSmSj8qo1CyyevwDocDp06p480XXnhhTzm8BttiA/OxoxBnr0Ocv7292HrCqC23tVVdJmttyTvG\nplCqRd47bn2GWu3z1atX95zjE54Hf38XuPZ7Ic7egHhnDlh7Dr8lHA+02MAd7lCXwTY104UvlLpT\nMGutoij6c9u9nLWW8AbV8e9zQEnE1aw38RV1n/Z8WVvhZraok2qzsw35DJeyd6FZa0uEEKIvV6VQ\ndho0ay2FsocoLWqEUhLf//73621CTqhdxdOINgHl27UnHP6ll14qe/92jv3BD35Ql+sWOnYn2kX/\nD0uzKx/U4Qvs386xhajXzVKIRrSL/h9WBqLsktQn4+Pj9TaBQmkotEVv69k1Dk+hUAqzJ7r0tWJj\nLr58OfsolHrBPvvss8/W24jtEIvF8Nxzz+H27dtIJpM4dOhQ3WwxrSnK/Pa3v0Vvby+ee+45eDwe\nOJ1OfPe730Vvb41W2G0gFovhtddeg9/vRyKRgN1ur/tvptkUCAQQjUYbwqb1DA0N4cEHH0QymWwI\nu7QVo6FQCA8++CAIIeXZpexwBgcHlWg0qiiKopw+fbrO1ijK1NSUbofb7da3r/9ca4aHh5WpqSlF\nURTF5XI1xG82NjamjI2NKZOTk8rg4GBD2KQRjUYVt9utTE1NNYxdU1NTuh2KUv59v+O79BMTE3pC\njUKJMOuBlquvlBTflcbj8ehZgxwOR0P8Zi6XCw6HA8PDwzh9+nRD2KQRDAb17E6NZNfFixfh8/kQ\nCoXKtmvHOzzQGE6VCy1nH4CCOftqwcjIiJ55uBF+M4fDAa/Xi76+voaxKRQK6TkcNRrBLofDoWeB\nGh4eLtuuHe/wjeZU63PxeTwe+P1++Hw+DAwM1NUuv9+PgYEBRCKRhvjNvF6vnq9QEISGsAlQW8tg\nMIiJiQkEAoGGscvn8+mOvZ3/wx3/WC4Wi+mJMLu7uyua6HK34Pf74fV6dWGQgYGBuv9moVAIkUhE\nT07a3d2aIT37AAACrElEQVRdd5s0YrEY+vr60NfXh/7+/oawSxAE/eVyuWC328uya8c7PIVCKZ4d\n36WnUCjFQx2eQtlDUIenUPYQ1OEpDUMoFCoqBNnr9eZ8FOX3+6suK77ToQ6/CwiFQujq6kIgENBn\n5CtxTp/Pt61zaGsKisXhcBSlbNTT0wObzYZQKAS/3w+/3w9AVTputFiMRoM6/C7A6XTC6XTC5XLB\n7Xbrj2+2e87tSHPHYjGEw+Ft2VCIkZERuN1ujI2NVfU6uwmaGH2XoK3UEwQBvb29ekCL1sV1uVyI\nxWIIBoOIxWK63NfIyAhOnjyJsbEx9Pb2YmxsLEvMU2t1T548iVAohNOnT+vn1aS+81UMXV1d+nN2\nv99f8Dz9/f0AsMnuXMpDfr8fbW1t8Pv9dQ953UnQFn6XYLPZ4HK54PF4MDU1hVAoBK/Xi+7ubt35\nnU6nLhOtafvZ7Xa4XC49oENz0u7uboTD4awyU1NTALLDT3NJhmnRana7Xe+il3KejXbnYmJiQq9o\n6rUKcSdCW/hdiDaeB6BH18ViMXi9Xhw/fhw2m00f62phmVvpxW8M3VxfgeRqfWOxmO7A68fkhc6z\nfvy93u6t7NJ6JJTioA6/C9Bmt7Uu/djYGC5cuAC32w2fz6c7VFdXFwRBQCQSweTkJDiO0+P+g8Eg\nBEHA5OQkADVeW9u+sYzNZsPIyAgcDgfsdntWWGcgEMDo6GjWAo+hoSE4nc6C59GuOTg4mGV3rljx\n3t5e+P1+9PX1NcQaih1DRRftUvYE69diezyemp9ndHS0rH0URaEtPKVkHA4HgsEgAODo0aN1Pw+l\neKjDU0rG7XbX9TwTExNwuVybuvJ+v5927wtAV8tRKHsI+liOQtlDUIenUPYQ1OEplD0EdXgKZQ9B\nHZ5C2UNQh6dQ9hD/H+tPc+aL7X5FAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "execution_count": 317, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAChCAYAAADju9JOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWlwG+eZ5/9vH7gIkjAoKRIj2yQYp2rKqVoLIr9lJ94I\n9MblJJOxQWqTzKRmszJk79Z+yCYSzM2HODWVjEgpOylvzUSEPFNbM0kqFGHncGJbIZgoE0/tbAhC\nsh3LiW02aR0UTREHL5x97IdmtwgSIA7iIvn+qiAA3W93P4T6ea9+3udPFEVRQKFQ9gRMvQ2gUCi1\ngzo8hbKHoA5PoewhqMNTKHsI6vAUyh6COjyFsofg6m1ApRgfH6+3CRRKQ3Hs2LFN23aNwwO5/8B6\nMjs7i/b29nqbsQlqV/E0ok1AYbvyNYC0S0+h7CGow1Moe4iCDr+0tFQLOygUSg0o6PAnTpzA1atX\na2HLjkWSJKRSKSSTScTjcaTT6XqbRKHkpOCk3cWLFwEAFy5cgCAIOH78OB566KGqG9boiKKIpaUl\nLC8vIx6PY+MaJI7jsLKyAoPBgObmZhiNxjpZSqHcpaDDP/XUU4hGo+jp6YHH44HdbseLL76Ixx9/\nvBb2NRyKoiAcDmNhYQGyLOctJ4oi4vE45ufnMT8/rzu+1WqFxWIBIaSGVlMoKgUd/ujRo3jyySez\ntk1NTVXNoEYmmUxidnYWyWSy5GPT6TTC4TDC4TA4jkNzczOam5thsVjAMHTulFIbCjp8d3c3lpaW\n0NLSgpmZGXR0dODUqVO1sK2hiMVimJub27JVLxZRFBGNRhGNRsGyLJqamtDc3IympiZw3K4KjaA0\nGAXvrsnJSRw5cgQAEAqF0NHRUW2bGo75+XksLCxU5dySJGFpaQlLS0sghMBkMundfrPZTFt/SkUp\neDdFo1HMzMxgaWkJv/vd7yp6cb/fX3RZQRAgCEJFr18IWZZx8+bNqjn7RhRFQSKRwJ07d/D+++/j\nnXfewczMDObm5hCLxZBMJjdNDlIopVCwhXe73Thz5gyi0SgGBwfLvpAgCAgEArDb7XC73fD5fPB4\nPIjFYvD5fHA4HHA4HHA6nXr5vr4+dHd3w+v1wuFw6MfUgng8jtu3byOVShV9jCzLeOWVV/Doo49W\npGWWZRnxeBzxeFzfxjAMeJ6H0WiE0WiEwWAAz/MwGAx0OEApSME7pLOzE+fPnwcAzMzMbOtimrMD\n6pgYgO7ENpsNXq9Xd3hAjQe22WxZxwuCAIfDsS078vHwww8DAEZGRhCJREpqTWVZxjPPPIPx8XH8\n5je/wZkzZ/R93/jGN7J+u46ODnzzm98sy0ZZlpFKpXJWRFplwPM8OI4Dx3FgWVZ/116SJEFRlLKe\nFGi/0eXLl8uyn1JfCjr82bNnIQgCFEXB5OQkJiYmyrqQ3W7H2NgYXC4XbDab7sgTExM4ffo0AGzq\nsmsxAN3d3XA6nXA4HAgEAhVv5UVRxOrqKtLpNObm5nDhwoWSjpdlGS+88ALeeustSJKEV199FTdu\n3MBjjz0Gi8WCYDCY5fALCwt48cUXK/o3lEIikdDnBwghOV8Acn6/efMmDh48iMXFxU1l1pdbvy/X\nZ431nyVJgiRJOfet30YfaZZPUY/ltFn56enpsi4SCATgcDjQ29uLYDAIh8MBu92u74/FYrDZbHqr\nDwAOh0N37JMnT2J4eFh3+Eqg3VyiKCKdTmN1dVVv+UrlzTffxNtvv63frJIk4e2330ZHRwd6enoq\nYm81yPfEYSvHVxQFkiRheXk5a782hCnk9Pn2Aeowav15N9q08XOuCmSr74UqnFyIoohMJpNzX75j\ntd8pH6VWWLnKy7Jc1hOjgg5//vx5jI6OorW1FVeuXMGlS5dKvsjY2BgGBwcRCATgcrkA3G3Ne3p6\nEIlEslp9QO3q9/f3w2azIRKJAAAikUjFuvNa99ZgMMBiscBms8FsNuPw4cP4/Oc/X9LY/XOf+xwS\niQR++ctfIpPJgOd5PPLII/jqV7+Kffv24ac//WlWC79v376qBy5pXfmN3XmWZbGwsID29nYwDAOW\nZcEwTNZrK15++WUAwOHDhytuczwez7oHGgGO48DzfL3N2EQx/1e5KOjwx48f1500GAyWbtnaObQZ\n+Y3/oR6PBz6fDzabDQMDAxAEAX6/H263G8FgEIIg6JOFoVAoa4xfabRxqSzLWFhYQDgcLqrFZxhG\nH7OPj4/j2LFjOHPmjN5j2fgosxKPNhmGgcFgyHpxHKe/b3UzpFIpNDc3l3VdOnbf2RR0eKfTicnJ\nSXR3d6Orq6usizidzk2Oqjm+zWbTx/Aa2veNrXklW/itYBgGBw4cQHNzM27evJm3S7fxmDNnzuSc\npS93gk6D53mYTCb9ZTQawfM8HctSSqZgn0DLnNHS0lLR5+Aej6ek8bggCHpPo1aYzWZ0dHTAZDIV\nVZ5hGDz22GPbfiRnMpnQ1taGw4cP46Mf/SgeeOAB3Hvvvdi/fz+am5thMBios1PKomALrygK7HY7\nlpaWMDk5iU9+8pMVu3gpDlyLlj0XPM+jo6MDt27dwvLyclWuoYXXWq1WNDU1NeSYkbI7KOjwLpcL\nXq8XhJBtBd7sZBiGweHDh4taJVcsPM/DarXSBTSUmpL3Lnv++ecBqOGvDocDnZ2d8Hq9NTOs0SCE\nYN++fejs7ITFYinrHAaDAW1tbejo6MADDzyAQ4cOwWq1Umen1Iy8LfzRo0cBqBNu3d3dAMqfpd9N\nGI1G3H///VhcXMTS0hJWV1fzzuQTQmCxWGC1WmG1WoueC6BQqkVeh9dWyLW1tYEQgpaWlrJn6Xcb\nhBA9bkCL0stkMlAUBbIsg+M4mM1mtLS0VOV5NYVSLgXH8MFgUE9ptVeXx24Fx3FobW3NuW995CCF\n0gjUdXkshUKpLTVbHkuhUOpPQYcnhOiz87RLT6HsbPJ26c+dO4fp6WmcPn0a4+PjCAQCGBsbq6Vt\nFAqlwuRt4Z1OJzo7O3Hy5EldpHG7CTAoFEp9yevwFy9exMWLFzE9PQ2/37/tBBgUCqX+5HX4vr4+\nHDt2TF/uCZSfAINCoTQGeR0+FArhypUrmJiYwJUrV6AoCoLBIEZGRmppH4VCqSB5Hf7UqVNYXFzE\nkSNH9NDaeq1Yo1AolWHLwJvW1lY4HA5MTk6CEKLH11MolJ1J3RJgUCiU2lPXBBiNjCwriEVSSCYk\nKIr6O3Acg6ZmDpYmHixLM85Qdh40AUYOUkkJC3eSyKSzE12IooRkUkJkIYUmK4eWVgMMRrZOVlIo\npVOU8owmCLEXWFpMIxpOYatktYoCrCyLWFkWYWni0GIzwGSijk9pfKgY2TpWVzKILBSfjx4A4qsi\n4qsijCYWLa08LE0cTTBJaVhKyq20m0NrUykJ4TvJ8o9PSrjzQRK3rq8iFkkhk9l+3jsKpdIUpS03\nNTUFALs2tFYUZdyZS6ACuSkhigpi0TRi0TSWV0RYm9IwmznwBpq3jlJ/aqItVyqa8sx6tEeClQ7+\nEUUZ83MJiGLlddczaWVtiJACxxGYzCyMRhYGEwuDgaFdf0rNKdjsnD9/Hk8//TSeeeYZPPXUUxU3\nQBAE+Hw+XYrK5/PB7XYjFoshEAhgaGgIACoqJKmRTEq4fSuOdCp/0y7LMn76s4vbTk0tigpWlkWE\nF1K4fTOO69MruH1zFQvzSSzF0kjERaTTEmS58hUPhaJRlLZcb2+vHktfDXLpxttsNjgcjqw1+JXU\nh//3H/8EMqKMH/3g5bxlZFnGV776JMYCv8Cvfn0Jf/udC2AYBgP/879jeuY9vVxnx0fwN9/+3yVd\nX1GAVEpGKkdlw7AE7NqL4wgYRv3MMAQMS9bUWgHCEDCEwNWrxkZcvvxr2mugbElBh7fZbDhx4gQW\nFxerkpc+n258LrajD68oCpIJCcmEhERCRCYjY/7OHEZG/0/O8rIs4+LFf8abv38dkiTiFy//GNev\nC+jv/0v8v4nXMD09pZfNdx5Nh73a3LhxCwf2H8T7wgoAgGGgVg5rL8IADCEAUTMYLcUkmI0pEKJW\nGgQAYdbLKd89d9Y29RR6gayqRd8HEP3LZnLVR9omWVIgSaX1cLZTvxVzrKIoZUmIV5ty7Sro8IuL\ni/pz+BdffLF0y7agkG78RrbTrSeEwGzhYLZwuAdG8DyT754EALz+RghvXXsTkiQCACRJxFvX3sTr\nb4TKun5ptqoOm+ud6I6rOhYhAMcRGE0s7ndYi2rh0yKLe9qMVf87SkXr2TQSGzXsG4Vy7Sro8NoM\nPQCEw+GSL7AVW+nG56Ki+vAcwYcPt+MLn//PkHJM2PU98SUk4gm88upPkcmkwfMGPPqpP8O3//o5\nfOEvHstq4Q/sP4jjfX+V0958FRghAG9gYDCoE3i8gQHLEXAcA4Yp7T9y5OKP1s7ZeDcmpbEoKmtt\nf39/VUJrt9KNB9QeQCgU0sftldSH13TOtVn6jRN3DMPgb79zAQAwFvgFel2P6WP4zo6PZJXd+D0X\nWbP0RhYGY+Vm6almO6VY6hpau5VuPKBKSq8fr1dDH57jGBw4aMbcrfimR3Oa07/0cz8+82m3rgFX\n7AQdbyCwtxlhsrAwGGjoLaX+NFw0SD7d+Grqw3Mcg/0fMiOXpiPDMPizz/YXLfjIcQSt9xjw4fua\n0LZfjbOnzk5pFBoylj6XY1c7247RxKJtvwl3PigvvNZoZNBiM9BYekpD05AOXy+arDwkSSm4Wm49\ndLUcZSdBF89soKXVgIPtFvB8/p+GEKDJyqH9sAUHDpqps1N2DHTxTA6MJhaHDlsQi6SQSEhqgIMC\nmvGGsuNpyMUzjQDDENj3meptBoVSUQo6/Pnz5zE6OorW1lZcuXIFly5dqoVdFAqlCjTE4hkKhVIb\nCk7aaYtntGg7CoWyc6nr4hkKhVJbCrbw1Vw8Q6FQaktdF89QKJTakreFf/755wGo+eUcDgc6Ozur\nkgCDQqHUjrwtvCYc6XQ6dfVYOktPoexs8rbwR44cAQC0tbWhtbUV0WiUztJTKDucvC389PQ0/H4/\ngsEgenp6oCgKIpHInhGTpFB2I3kdvrOzE263Gw6HA06nE3a7Ha2trbW0jUKhVJgtZ+k7OzvR2dmp\nf7969SoeeuihqhvVqIiyjFg8A1FWkJZkSLICnmVg4hiYeBYWA6tmh6VQGpSCj+X6+/vR1taGcDiM\n6enpPbFabj2KoiCWyOCD5RTCqylIW+hRMASwGjk0mzi0mniIVFSC0mAUdPgLFy7oXfnx8fGqG9RI\nxNMi3plfwVJSLKq8rABLSRFLSRG3kEQ0soowYrCZedjMPFpMPNgSM9JSKJWkKIdfH2F37NixqhrU\nCCiKgpuxBN6PxLGdRlqBguWkiOWkiBvRBBgCNBs5tJh5NBs5tJp58GzDpRWk7GKKylr75JNP7pkJ\nu5Qo4Q9zy1gsslUvBVkBFpNi1rnNPINmE68OBYwcmgwsOFoJUKpEQYd/4okn9EUzjz/+eNUNqifR\neBp/+GAZmRLljrZDIiMjkUlhfjmlbzNyDKxGDua1iUDtnfYGKNulqBRXNpsNdrsd586dw9e+9rW8\nZQVBgNfrRU9PD06fPr1p/9DQEJxOJ4aHh9HT06M/9qs3sqLgejSOG9FE0ckrq0lKlJES05u2cyyB\niWNg5lmYeBZGjoGR094ZWiFQClJUiist2OZXv/rVlmU1nbhcTqypwrpcLoyNjeWsEOpBeDUNYWEF\niUzxctCyLGP85z/GsU//edH56iuBKClYkSSspKSc+xkC8CwDA8fAwDLgGAKeZcCxBDxDwK1tW0lJ\niKdFsAwByxBwNfwbKPUlp8NrUXYAMDExgVAoVHKknSAIuqCEy+WCIAiYmJjQ5aM0CSlNDtrr9WaV\nDwQCmJycxPDwMPr6+vDFL34RVqsVoVBo25WFrCj40088jJQo4Tv/9JPSjpVlfOtr/xWvBV7Bv10e\nw9fP/X1JTn/26/8DN2buLjm+t6MLp771v0qyIa9titY72LryikTimJNiWdtYhoAl6jvDEHBkTX2W\nqJUCQ1QFWvW1pk4LVc+OIdr7+s9337XYBGadACZZ264pzjKEQF5TRKUh3KXx8MMPAyhOciynw8di\nMbjdbtjt9rKkmQHVgQcGBhCJRCAIArq7uxEKheByuTA6Ogqn0wmHw4GJiQkMDg6ir68vq7zH40Ff\nXx8AYGBgAIFAAOFwGMePHy/LHgBYSYmYCa9iMZnBakpE5M4HeHX0n4s+XpFl/Gz0h/jj71+HJEn4\n9cs/wez1aXy27wsgOZw+l1z067/7V9yYuSuYGZmfK8mGSlArGetSWW+X5vOaPLVWBehVgVZ5ZG1E\ndrkcstZ55apzSWADiMdXYbE0bZLGznmODd/y1Vtk04e8p817ldXVFTQ1WQEAN27ewoGDh3Arligo\ngZ3T4bWFM6UwNDQEt9utd90B6GG567dtpK2tLW/5np4e+Hw+eDweXXPO6/VidHS0ZPsAdTLM3mTQ\nJZdL5dobV/HutTchSWqXWpIkvHvtTVx74yoefKgyIpe7gbvOutaaZ20nWc683tmUjDo/kc/Z727T\nXSpbyz5r392t2ZVELltzOzsAyBwDSw6psLy3T56KY719RZ8rq0x2KYljYOFZtXfEEPAsQVuToeB5\n8o7hn3/+eZw4cQJnz55FJBKBoiiYnp7GyMhIzvJutxuCICAcDsPlcsHhcMDn8+kteSgUwsTEBARB\nQDAYhCAI+rZYLIbBwcGs8jabDW63G8PDwwDUdflOpxO9vb1F/Dy54VkG7a1mtLea0WriYW5vx6f/\n05cgFjkr/8gTX0Q8kcDlV1+CmEmD4w34xKc+g6/89XdydusjkTDs9rasbS//xJ/VwtsPHMSn+v6y\n7L+pVBgCxKIRfGj/vrUxvDqu17rt2rie1bvyd7vqLEOyuunaOyHZ3f1yu+Szs7Nob2+v8F+8PRrR\nJiDbrtERVS7cxBcWRCGKknte+sqVKzhy5AjGx8ez1sPXMvAmFotBEISiJKLHx8fLsk2UZMxE4ri9\nlCxqhn79GP7jrke3HMPncvhKj+EJAQzs3Vl67d3AahN2DHiWgFubnGMZsiNu4kahEW0CCtuVzx/y\ntvBat97v9+sH9vT0bNfOohkaGoLD4YDb7a7qdTiWwUf2W3Gg2Yi355YLTngxDIOvn/v7smfpy3Fu\njiWw8Kz+OM7EMzCtexxHJ7koxVLwsZyiKBgYGIDdbkcwGMzbpa80tX5s12Li4bzXhj98sIxoPLNl\nWYZh0PvZJypuA8cSNPEsmtYi7iwGjgbcUCpKQYfv7e3V5ZsbIUimmvAsg48dasGNWALXtxlHX4j1\nK+uajWpsvZnqyFOqTEGH7+rqQmtrK2ZmZnDPPffUwqa6QgjBffdYYLcY8M78ct4gl1JhGbUX0Wrm\n0Wri0Wzi6Np5Ss2hKa7yYDVyOHLYhhuxBG5EE5BKbO7VlXEs7rdbYDNTB6c0BkWnuGpra0NLS0st\nbas7Wmv/4VYzwqtpzC8nEU1kcs7mGzkGTQYWLWYeLUYOzSYeH8yl0W631N5wCiUPBVNcvfDCC+jq\n6oLX68UjjzyCEydO1Mq2hoFlCA40G3Gg2QhFUZCRFGRkGaKkqI/EeIa23pQdQcEx/JEjRzAyMoLB\nwUEIglCo+K6HEAIDR2AorNJFoTQcRd21bW1tIITgypUr1baHQqFUkYJSU6FQCAsLCzh//vyeS2BJ\noew2CkpNHT16FN3d3VAUhUpNUSg7nC1Da5966ikEg0EQQtDd3Y3vfe97tbSNQqFUmLwOf+7cOXi9\nXl2IYnp6umCKKwqF0tjkHcNvVJ3Z+J1Coew88jp8rjDavRBaS6HsZvJ26T0eD7q6uvTvWgKMd999\ntyaGUSiUypPX4UdHRzeluqLP4SmUnU3eLn2uvHbl5LqjUCiNQ8HQWkpuFFmGkoyvfSMAw4DwBhCW\nrmmnNC7U4UtAjq9CujMHeWUZ8uoKoGxOh0UMRhBzExhLE5TVOJT9+0F4vg7WUiiboQ5fBHJ8FeKt\n65CiCyiU6VJJp6CkU5AXI0AkguRSGMRkBmNtAdPSCqbZBsZkqpHlFEo2BR3+6aef1iPslpaW9tSa\neEWWkbkuQPpgdnvnSSYgJROQFj4AABCDAYy1FUyzWgkQcxNNREmpCQ2bxLLeyIk40u+9DSW+WvFz\nK+k0pMgdSJE76gaWBWOxgmluBdNkBWmygjHSXgCl8hSVxLK3txeKouz6JJYa4p05ZGbeA+TiBSa3\nhSRBXl6EvLyobyK8AcSizgXo70YznRSkbIuCDm+z2XDixAksLi7C6/VW7MJ+v7/onPNa4o1qVziK\nJCIz/R6k8HxVr1OULZk0lMU05MVo1nZ1UtACYjSBMZlBjCb1ZTDSyUFKQQomwFhcXMTFixdx6dKl\nLTXiNhIKhdDV1YVQKARAFZbQKgyfz6fr0A0NDcHv9+vlAFVxJhAIYGhoCIDq6JqybLWQV1eQ+v2V\nopxdlmWMXBqDXKsewDrUCcEopPnbyFwXkH73GlK/DyEZ+r9ITLyG5OsTSF27ivS715CeeQ+ZW+9D\nnL8NKbIAaSkGOb4KJZNBHsEhyi6nYAs/NXVXFikcDhd9Yk0jTpOJ8ng8iEQiAO5qxa8XivR6vXpZ\nm80Gh8OhS0kDgN1uhyAIFWvlNYndX738c4izN9XxdBFOIMsy/suz38JL//IaXv3Xf8M/PPv1ktRn\n/tvfnMW712/o3x+471783cCpku3PYxyUZAJKMrF1uUgEybn3AZYFYTmA50FYDoTjAJZThw3aPpZV\nYwwYVt3HMGrObYYDYRmAYQCG3XWTjtr98cMf/rC+hlSYgg7vdrvR398PQggGBwcrclFNCXZiYkJX\nmCmUL09r5cuVr9ZQRBFKIg4lk8bt2dv43plvF32srCj4p5/8DFff/iNEScIL47/G1OwsvvS5z+ZM\nYplLlvm3V17He+scfi4cwT/+4tXy/6AyqI5c9Jro+/oXsKb+uibbStb0Yjd+Xjs2Ho+jqakp69is\nd6w/r/7PWpGN2ssk9+dc3zfJyhLM3rqF9kOHoCQTkJOJu9fSTcp3/o02b961qQzZ4u+oMAXVY/1+\nv96qer3ebc/SC4IAu92uf4/FYrDZbAWHC5Xo1ivpFKTwHciJOJR0GkqOwJmtmHzrGt7447sQ1+Si\nRUnCG398F5NvXUPPxx7clm07H0XtIeXoJBU3eCBAMgEZSn6H16WYSbbzEEDZJBKfu2LIuS2Xk8kS\nlHQKiNyBxDGby2VrYOc8X0Hn3VAxbjrPxq/ryijRMCQjn2Pf1hRMceV0OrPUY0tBa8kBZHXHtda8\np6cHkUgENpstq2wuIpHItrvzxGAEd+gwAIBpsuLDZgue/PKXId2ZK2pG/q8efQTJ1Th+/OvLSGdE\nGHgOf/4fPoG/O/WVPHLRkazKDQB+9POXs1r4g212fPmxT23r7yqVXHZlwTBq955Z69pveL/bxV97\nZ5i72whRy6zdzIQwqioHYe5uY5i17wAIo35HYym1XnzlEgCAtN8HvkFsWg9Ji2Db9pd8XFXVYwcH\nB3UVWJvNtkn22ePxwOfzwWazYWBgAIIgYHh4GIODgwgEAgiFQnpFEQqFipKNLpbLly/rn5UP34/M\nrfcLBtgwDIN/ePbrAICX/uU1fOZPP17yGP6B++7d8ntFIYw6c8/zIJwBhOdAOB4wWMC3t6ufubVt\n+hh+943Hy0G7P2Zntxd01WhUNfDG4XDkVIHVWnObzbZp//HjxwGolcH68XolWvh8EJ6HoeMjkFps\nSAvvAJKYt6zm9KNj4+jrPVayXHTFJugIc/eRnNGoPpYzmkAMBvUZPm9QHTjXofwsuA81XqtFqT51\nUY/1eDwIBAL6eQshCELRZbcDa98Hk7UZ6ffehry8lLccwzA4/h97q26PahR3N/jGZFafwZvMqoPT\nlphSIls6/NWrV9Hb26vHzz/xROU00fM5cK5uey0j/IjBCMOf/DuIt96HOHujqEd1FYPjwTQ1g7Fa\nwVhoiC2l8uR1+LNnz2JqagrRaBQDAwN46KGHamlXXSGEgD/cAabFhszUH6Ck09W4CIjZosbPW5vV\n1XSmSj8qo1CyyevwDocDp06p480XXnhhTzm8BttiA/OxoxBnr0Ocv7292HrCqC23tVVdJmttyTvG\nplCqRd47bn2GWu3z1atX95zjE54Hf38XuPZ7Ic7egHhnDlh7Dr8lHA+02MAd7lCXwTY104UvlLpT\nMGutoij6c9u9nLWW8AbV8e9zQEnE1aw38RV1n/Z8WVvhZraok2qzsw35DJeyd6FZa0uEEKIvV6VQ\ndho0ay2FsocoLWqEUhLf//73621CTqhdxdOINgHl27UnHP6ll14qe/92jv3BD35Ql+sWOnYn2kX/\nD0uzKx/U4Qvs386xhajXzVKIRrSL/h9WBqLsktQn4+Pj9TaBQmkotEVv69k1Dk+hUAqzJ7r0tWJj\nLr58OfsolHrBPvvss8/W24jtEIvF8Nxzz+H27dtIJpM4dOhQ3WwxrSnK/Pa3v0Vvby+ee+45eDwe\nOJ1OfPe730Vvb41W2G0gFovhtddeg9/vRyKRgN1ur/tvptkUCAQQjUYbwqb1DA0N4cEHH0QymWwI\nu7QVo6FQCA8++CAIIeXZpexwBgcHlWg0qiiKopw+fbrO1ijK1NSUbofb7da3r/9ca4aHh5WpqSlF\nURTF5XI1xG82NjamjI2NKZOTk8rg4GBD2KQRjUYVt9utTE1NNYxdU1NTuh2KUv59v+O79BMTE3pC\njUKJMOuBlquvlBTflcbj8ehZgxwOR0P8Zi6XCw6HA8PDwzh9+nRD2KQRDAb17E6NZNfFixfh8/kQ\nCoXKtmvHOzzQGE6VCy1nH4CCOftqwcjIiJ55uBF+M4fDAa/Xi76+voaxKRQK6TkcNRrBLofDoWeB\nGh4eLtuuHe/wjeZU63PxeTwe+P1++Hw+DAwM1NUuv9+PgYEBRCKRhvjNvF6vnq9QEISGsAlQW8tg\nMIiJiQkEAoGGscvn8+mOvZ3/wx3/WC4Wi+mJMLu7uyua6HK34Pf74fV6dWGQgYGBuv9moVAIkUhE\nT07a3d2aIT37AAACrElEQVRdd5s0YrEY+vr60NfXh/7+/oawSxAE/eVyuWC328uya8c7PIVCKZ4d\n36WnUCjFQx2eQtlDUIenUPYQ1OEpDUMoFCoqBNnr9eZ8FOX3+6suK77ToQ6/CwiFQujq6kIgENBn\n5CtxTp/Pt61zaGsKisXhcBSlbNTT0wObzYZQKAS/3w+/3w9AVTputFiMRoM6/C7A6XTC6XTC5XLB\n7Xbrj2+2e87tSHPHYjGEw+Ft2VCIkZERuN1ujI2NVfU6uwmaGH2XoK3UEwQBvb29ekCL1sV1uVyI\nxWIIBoOIxWK63NfIyAhOnjyJsbEx9Pb2YmxsLEvMU2t1T548iVAohNOnT+vn1aS+81UMXV1d+nN2\nv99f8Dz9/f0AsMnuXMpDfr8fbW1t8Pv9dQ953UnQFn6XYLPZ4HK54PF4MDU1hVAoBK/Xi+7ubt35\nnU6nLhOtafvZ7Xa4XC49oENz0u7uboTD4awyU1NTALLDT3NJhmnRana7Xe+il3KejXbnYmJiQq9o\n6rUKcSdCW/hdiDaeB6BH18ViMXi9Xhw/fhw2m00f62phmVvpxW8M3VxfgeRqfWOxmO7A68fkhc6z\nfvy93u6t7NJ6JJTioA6/C9Bmt7Uu/djYGC5cuAC32w2fz6c7VFdXFwRBQCQSweTkJDiO0+P+g8Eg\nBEHA5OQkADVeW9u+sYzNZsPIyAgcDgfsdntWWGcgEMDo6GjWAo+hoSE4nc6C59GuOTg4mGV3rljx\n3t5e+P1+9PX1NcQaih1DRRftUvYE69diezyemp9ndHS0rH0URaEtPKVkHA4HgsEgAODo0aN1Pw+l\neKjDU0rG7XbX9TwTExNwuVybuvJ+v5927wtAV8tRKHsI+liOQtlDUIenUPYQ1OEplD0EdXgKZQ9B\nHZ5C2UNQh6dQ9hD/H+tPc+aL7X5FAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "beaming_violins(o2a_traces, priors)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Assuming VT is distributed\n", + "We can also incorporate the uncertainty in VT." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "sigma_vt = 118342" + ] + }, + { + "cell_type": "code", + "execution_count": 319, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import theano.tensor as T\n", + "from pymc3 import DensityDist, Uniform, Normal\n", + "from pymc3 import Model\n", + "from pymc3 import distributions\n", + "\n", + "def grb_model_non_galaxy_sigma(number_events, background_rate, \n", + " observation_time, volume, sigma_volume, grb_rate,\n", + " efficiency_prior = \"uniform\"):\n", + " with Model() as model:\n", + " signal_rate = pm.DensityDist('signal_rate', \n", + " logp=lambda value: log_signal_rate(value, number_events, background_rate, observation_time),\n", + " testval=50)\n", + "\n", + " volume = pm.Normal(\"volume\", volume, sigma_volume)\n", + " \n", + " cbc_rate = pm.Deterministic('cbc_rate', signal_rate / volume * 1e9) \n", + " \n", + " # Allow the efficiency prior to be switched-out\n", + " if efficiency_prior == \"uniform\":\n", + " efficiency = pm.Uniform('efficiency', 0,1)\n", + " elif efficiency_prior == \"jeffreys\":\n", + " efficiency = pm.Beta('efficiency', 0.5, 0.5)\n", + " elif isinstance(efficiency_prior, float):\n", + " efficiency = efficiency_prior\n", + " \n", + " def cosangle(cbc_rate, efficiency, grb_rate):\n", + " return T.switch((grb_rate >= cbc_rate*efficiency), -np.Inf, \n", + " (1.0 - ((grb_rate/(cbc_rate*efficiency)))))\n", + " \n", + " costheta = pm.Deterministic('cos_angle', cosangle(cbc_rate, efficiency, grb_rate)\n", + " \n", + " )\n", + "\n", + " angle = pm.Deterministic(\"angle\", theano.tensor.arccos(costheta))\n", + " \n", + " return model\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 320, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Actual O2\n", + "priors = [\"uniform\", \"jeffreys\", 1.0, 0.5]\n", + "number_events = 1 # There were no BNS detections in O1\n", + "background_rate = 0.01 # We take the FAR to be 1/100 yr\n", + "grb_rate = 10.0 \n", + "o2a_models_u = []\n", + "for prior in priors:\n", + " o2a_models_u.append( grb_model_non_galaxy_sigma(number_events, background_rate, total_time, total_v, sigma_vt/total_time, grb_rate, prior))" + ] + }, + { + "cell_type": "code", + "execution_count": 321, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 1000000/1000000 [03:37<00:00, 4607.29it/s]\n", + "100%|██████████| 1000000/1000000 [03:43<00:00, 4465.80it/s]\n", + "100%|██████████| 1000000/1000000 [02:31<00:00, 6610.55it/s]\n", + "100%|██████████| 1000000/1000000 [02:26<00:00, 6843.46it/s]\n" + ] + } + ], + "source": [ + "samples = 1000000\n", + "o2a_traces_u = []\n", + "for model in o2a_models_u:\n", + " with model:\n", + " step = pm.Metropolis()\n", + " o2a_traces_u.append(pm.sample(samples, step))" + ] + }, + { + "cell_type": "code", + "execution_count": 322, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "|\t\t\t| Lower\t| MAP\t| Median\t| Upper\t|\n", + "|----------|\n", + "| U(0,1) \t| 3.67\t| 7.21\t| 11.83\t| 39.75\t|\n", + "| Jeffreys \t| 3.42\t| 7.45\t| 11.77\t| 50.89\t|\n", + "| $\\delta(1)$ \t| 3.44\t| 6.26\t| 7.67\t| 17.07\t|\n", + "| $\\delta(0.5)$ \t| 4.88\t| 8.02\t| 10.86\t| 24.20\t|\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAChCAYAAADju9JOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnVtwG9eZ5/+nL7gRJCDQknWJbLE5Th48NbEg8m1m4o1A\nb1yeZBILpGadndTMjgzZu7UP2Y0Ec/Ngp6YyY1LKTspbMxEhZ3ZrcqlQhJ2LE9sKwdhJPFUZE4Rk\nK5bHVthkdKF1w4WUCIAAunsfmt0ESIANgriRPL8qEEDjdPcHsL9z6+98f6IoigIKhbIlYBptAIVC\nqR/U4SmULQR1eAplC0EdnkLZQlCHp1C2ENThKZQtBNdoA6rF2NhYo02gUJqKgwcPrti2aRweKP4F\nG8nMzAx2797daDNWQO0qn2a0CTC2q1QDSLv0FMoWgjo8hbKFMHT4ubm5ethBoVDqgKHDHzlyBOfP\nn6+HLRuWhYUFJJNJJJNJzM/PI5vNNtokCqUohpN2Z86cAQCcPn0aoiji8OHDeOihh2puWLOTTqcx\nNzeHubk5ZDKZFZ9zHIe7d+/CbDbDbrfDbDY3wEoKpRBDh3/qqacQj8fR3d0Nn88Hl8uFl19+GY8/\n/ng97Gs6JEnCzZs3EY/HVy2Xy+WQTCZx48YN3LhxAyaTCXa7Ha2trbDZbCCE1MliCmUJQ4c/cOAA\nnnzyyYJtk5OTNTOomUkmk5iZmSnaohuRyWQQi8UQi8XAcRxaW1t152cYOndKqQ+GDt/V1YW5uTm0\ntbVhenoa+/btw7Fjx+phW1MRjUZx8+ZNVCN9QC6XQzweRzweB8uyaGlpQWtrK1paWsBxmyo0gtJk\nGF5dExMT2L9/PwAgEolg3759tbapqVAUBTdu3EAsFqvJ8SVJ0ucCCCGwWCxoaWmBzWaD1WoFy7I1\nOS9la2LYl4zH45iensbc3Bzefvvtqp48GAyWXVYURYiiWNXzGyHLMq5evVozZ1+OoihIpVK4ffs2\nLl++jA8//BCiKGJmZgbxeBypVAqyLNfFFsrmxLCF93q9eP755xGPxzEwMFDxiURRRCgUgsvlgtfr\nRSAQgM/nQyKRQCAQgCAIEAQBbrdbL9/b24uuri74/X4IgqDvUw/S6TRmZmaQTqdLlpFlGa+99hoe\nffTRmozDFUVBOp1GOp1GIpEAABBCwPM8zGYzTCYTTCYTeJ7Xn+lkIGU1DB2+o6MDp06dAgBMT0+v\n62SaswPQL2DNiZ1OJ/x+v+7wgBoP7HQ6C/YXRRGCIKzLjlI8/PDDANRbkdFodNXxuizLeOaZZzA2\nNoZf/vKXeP7550s6/bPPPlvw2+3btw9f+9rXKrJRURRkMpmiE4daZWAymcCyLDiO0x8Mw+jbJEmC\noigVVQ7ab/Tmm29WZD+lsRg6/IkTJyCKIhRFwcTEBMbHxys6kcvlwujoKDweD5xOp+7I4+PjOH78\nOACs6LJrMQBdXV1wu90QBAGhUKjqrbx2Cy2TyeD69esIBAKrlpdlGS+99BLee+89SJKE119/HVeu\nXMGhQ4cKnD6VSsFqtSIcDhc4/O3bt/Hyyy9X9TusBc0uhmF0pyeE6LYTQvTty7ddvXoVO3fuRCKR\n0LflH6PUvsWeNbT3kiRBkqRVy9AezPoo67acNis/NTVV0UlCoRAEQUBPTw/C4TAEQYDL5dI/TyQS\ncDqdeqsPAIIg6I599OhRDA0N6Q5fDbSLS5IkLCwsYH5+Xm/5jLhw4QLef/99/eKUJAnvv/8+Lly4\ngE9+8pNVsa8erDYfkO+8+Q9FUSBJEu7cuVPUwZc7v/a62PPycyWTSdy9e7dgezG7jI5ZTuVQqkJZ\nTi6XKytyspKKqNx9ipWTZbmi+RxDhz916hRGRkbgcDhw7tw5nD17ds0nGR0dxcDAAEKhEDweD4Cl\n1ry7uxuxWKyg1QfUrn5fXx+cTqc+aRaLxarWnWdZVp8Bt1qtcDqdsFqt2Lt3L5544olVx+6f//zn\nkUql8POf/xzZbBY8z+ORRx7Bs88+W9DCx2IxuFwu/PjHPy5o4e+5556aBy5p3y+/S689R6NR7Nq1\nCyzLgmGYgkd+S1+MV199FQCwd+/eqtucTCbhcDiqftz1wHEceJ5vtBkr0P5fa8XQ4Q8fPqw7aTgc\nXrtli8fQZuTznRoAfD4fAoEAnE4n+vv7IYoigsEgvF4vwuEwRFHUJwsjkUjBGL/aaONSWZYRjUYR\njUaL1qIMw+D5558HoM4zHDx4cNUx/PJbmdW4tZk/Xl8+cWcymVa9GDKZDNra2io6Lx27b2wMHd7t\ndmNiYgJdXV3o7Oys6CRut3uFo2qO73Q69TG8hvZ+eWtezRZ+NRiGwfbt29Ha2oorV64U7dJpTl/O\nLH2lE3QaHMfBYrHoD7PZDJ7naYQeZc0YXjFa5oy2traq3gf3+XxrGo+Loqj3NOqFxWLBvn37YLFY\nin7OMAwee+yxqjue2WzGtm3bsGfPHjzwwAP4+Mc/jvvuuw87duxAW1sbzGYzdXZKRRi28IqiwOVy\nYW5uDhMTE/j0pz9dtZOvxYHr0bIXg+d53H///bh69Srm5+drcg4tvNZut6OlpaUpx4yUzYGhw3s8\nHvj9fhBC1hV4s5FhWRZ79+7F7du3De/PlwvP82htbYXdbqcLaCh1o+RV9uKLLwJQw18FQUBHRwf8\nfn/dDGs2GIbBjh07cP/998NkMlV0DJ7n0d7ejn379uGBBx7Azp07YbfbqbNT6kbJFv7AgQMA1Am3\nrq4uAJXP0m8mbDYbBEFANBrF3NwcFhYWSpYlhMBms8Fut8Nut5ecC6BQ6kVJh9dWyLW3t4MQgra2\ntopn6Tcb2iz+9u3bkU6ncffuXT0IB1BbcovFAofDgT179jTQUgqlEMMxfDgc1lNabcXlsUZot8qK\nkR85SKE0Aw1dHkuhUOpL3ZbHUiiUxmPo8IQQfXaedukplI1NyS79yZMnMTU1hePHj2NsbAyhUAij\no6P1tI1CoVSZki282+1GR0cHjh49qos0rjcBBoVCaSwlHf7MmTM4c+YMpqamEAwG150Ag0KhNJ6S\nDt/b24uDBw/qyz+ByhNgUCiU5qCkw0ciEZw7dw7j4+M4d+4cFEVBOBzG8PBwPe2jUChVpKTDHzt2\nDLOzs9i/f78eWtuoFWsUCqU6rBp443A4IAgCJiYmQAjR4+spFMrGpGEJMCgUSv1paAKMZkaWFcRj\nC0jN5wCiBiBxHEGLnYethQPD0HTJlI0HTYBRhHRaQvRmGtlsfgJLBdkMkEpKYBigxc6j1cHDZKLa\nb5SNQ1nKM5ogxFZgLpFBPLaA1ZLayDJwZy6LO3NZWG0s2pwmWK1U9ZXS/NCrNI/5u1nEoqUTWhQj\nlZSQSqZgMjGwt/Gwt/K0u09pWtaUW2kzh9YupCVEb5UWnzAik5ERu72Aq7+/i+itNBbSkvFOFEqd\nKUtbbnJyEgA2bWhtLifj5o0UqqHEnN/dn5vLosWWgdXG0rE+pSmoi7bcWtGUZ/LRbglWO/gnm5Vx\n83oKUm79mWiXk8sB8egC4lGA4wisNg5mCwuzmQVvookrKfXH8Ko7deoUnn76aTzzzDN46qmnqm6A\nKIoIBAK6FFUgEIDX60UikUAoFMLg4CAAVFVIUiOVyuH6tSSymdJNuyzL+PFPzlQk3JdPLqfgzlwW\nt2+mce3KPC5P38X1mSSit9KYS2SQSuaQzcpVSYFNoZSiLG25np4ePZa+FhTTjXc6nRAEoWANfjX1\n4f/kTz6FbFbGD773askysizjy//zSYyGfoZfvHEW//CN02AYBv3/679javp3ermOfX+Av/+7/7Om\n88uSgnRKQjq1cqzPcQQMS8AuPhhm6ZlhyaLgI8AwBIQh8HjU2Ig33ngDhFBJZUppDB3e6XTiyJEj\nmJ2drUle+lK68cVYjz68oqgOlkrmkJzPIZuRcfPWdQyP/L+i5WVZxpkz38GF374DScrhZ6/+EJcv\ni+jr+0v82/hbmJqa1MuWOo6mw15rrly5hh3bd+LylCq1rFUG+oMlixWBWhncmZUQtyyAIQAIAQHU\n4CLtmajlgbzX2ufqRiwWXSJvf/VJf4FluxWFAJAkBZK0th7Oeuq2cvZVFKUpe12V2mXo8LOzs/p9\n+Jdffnntlq2CkW78ctbTrSdEHUNbbRxc9wA8z2C1//c770bw3sULkKQcAECScnjv4gW8826kovOv\nBYaB3oprDkcYAkZrvfOcl0DtEZgtLO7rsJfVwi9kWWxzmWv+PdaK1qNpJvJ17puJSu0ydHhthh4A\notHomk+wGqvpxhejqvrwHMGej+3Gf/qLv4ZcpFXpPfQlpJIpvPb6j5HNZsDzJjz6mT/H3/3tC3ji\nPz9W0MLv2L4Th3v/qqi9pSowhgFMJnXyjucZ8CYGHMeA49f+jxw+84PFYzbfhUlpLsrKWtvX11eT\n0NrVdOMBtQcQiUT0cXs19eE1nXNtln75xB3DMPiHb5wGAIyGfoYez2P6GL5j3x8UlF3+vhj5s/Qm\ns+rk1Wo5qGY7pVwaGlq7mm48oEpK54/Xa6EPz/MM7t1lxUfXkituzWlO/8pPg/jsn3l1DbhyJ+g4\nDtjmMsFq42Ay0/vwlMbTdDeDS+nG11IfnuMY7LjXimKajgzD4M8/11e24CPDAK1tPHbtseGee3k4\ntpmps1OahqaMpS/m2LXOtmO2sGjfbsGtG5WF1/ImBq00lp7S5DSlwzeKFjsPKacYrpbLx2pj0eZQ\nu+0USrNDF88so81pwr27beD50j+N1m3f/TEb7t1lo85O2TDQxTNFsFhY7PqYDfHoAlLJxYw3IGA5\nAnsrzXhD2bg05eKZZoBhCNq3F5eBplA2KoYOf+rUKYyMjMDhcODcuXM4e/ZsPeyiUCg1oCkWz1Ao\nlPpgOGmnLZ7Rou0oFMrGpaGLZygUSn0xbOFruXiGQqHUl4YunqFQKPWlZAv/4osvAlDzywmCgI6O\njpokwKBQKPWjZAuvCUe63W5dPZbO0lMoG5uSLfz+/fsBAO3t7XA4HIjH43SWnkLZ4JRs4aemphAM\nBhEOh9Hd3Q1FURCLxbaMmCSFshkp6fAdHR3wer0QBAFutxsulwsOh6OetlEolCqz6ix9R0cHOjo6\n9Pfnz5/HQw89VHOjmhVJVjCXziIjycjkZORkBSaWgZljYOFZ2EwsGDrsoTQxhrfl+vr60N7ejmg0\niqmpqS2xWm45d9JZXJ9bwM27C5Dk0gvlGQK0mjk4rDwcVn7VshRKIzB0+NOnT+td+bGxsZob1Eyk\nsxIu3bqLeDJbVnlZAWbTOcymc0A8hUR8HnEyC6eNxzYrD7uZoxOflIZSlsPnR9gdPHiwpgY1Cx/N\npiBGk+tqpWVFQSKVRSKVxTQAjiVwWNTW32HhaAVAqTtlZa198sknt8yEXVaS8e837pTdqq+FnKQg\nOp9BdD4DAGAZglYzh1YLh1azWgFYeJrwklI7DB3+0KFD+qKZxx9/vOYGNZJEKosPbtzBQq4KutFl\nIMlLPQANjiVoMXFoMamTgDaehdXEwszRioCyfspKceV0OuFyuXDy5El85StfKVlWFEX4/X50d3fj\n+PHjKz4fHByE2+3G0NAQuru79dt+jUZRFFxJpPD7WLLs5JW1IicpmE1lMZsq7GGwDGDhWFh4Fhae\ngZllYOJYmDn1LgHPMmBp2i2KAWWluNKCbX7xi1+sWlbTiSvmxJoqrMfjwejoaNEKoRHEkxmIt+cx\nn1mp4loKWZYx9tMf4uCffaHsfPXrRZKB+Yy0qp0sQ2BiCUwsA45lwLMEPMuAY5aeWYZgPiMhlZHA\nMOo+XJ2+A6XxFHV4LcoOAMbHxxGJRNYcaSeKoi4o4fF4IIoixsfHdfkoTUJKk4P2+/0F5UOhECYm\nJjA0NITe3l588YtfhN1uRyQSWXdlISsK/vRTD2MhJ+Eb//Kjte0ry/j6V/4r3gq9ht+8OYqvnvyn\nNTn9ia/+D1yZXlpyvHdfJ459/X+vyYZSSLKClKwglV19SBKLJfFRLl6wjWUAlqgVAsMQcHmvWULA\naBLVi+KW7KKYIUOgxx4wZEnwsuAZ+a+XlGyZRTFMTTCzGVVam52HH34YQPlyY0UdPpFIwOv1wuVy\nVSTNDKgO3N/fj1gsBlEU0dXVhUgkAo/Hg5GREbjdbgiCgPHxcQwMDKC3t7egvM/nQ29vLwCgv78f\noVAI0WgUhw8frsgeAEgkM7iaSGE2ncX8Qg6xWzfw+sh3yt5fkWX8ZOT7+OC370CSJLzx6o8wc3kK\nn+t9AqSI0xeTi37n7X/FleklwczYzetrsqEa1EvGeq1odulS1drfPKnqJcnqpc8LZapLy1oXuyGy\nvGzBNgDJ5DxstpaC85Y8Rr7VBqOr5VLbq7wtuuf8/F20tNhx5eo17Ni5C9cSqbLOW9ThtYUza2Fw\ncBBer1fvugPQw3Lzty2nvb29ZPnu7m4EAgH4fD5dc87v92NkZGTN9gGA1cTC1WICwxDDH6YYF989\nj0sXL0CS1G61JEm4dPECLr57Hg8+VB2Ry80CWWzNNcdclKIHQAocmiwWJgCULAObidU/07Xr84+r\n/8lzeFL42XLXXbm90M6VeywhcwysPLvieil5+SyrZIqWX+XaW+2yzP81JI6BjWfBMAQ8S9DeYlpl\nzyVKjuFffPFFHDlyBCdOnEAsFoOiKJiamsLw8HDR8l6vF6IoIhqNwuPxQBAEBAIBvSWPRCIYHx+H\nKIoIh8MQRVHflkgkMDAwUFDe6XTC6/ViaGgIgLou3+12o6enp6wvVgwzx2K3w4rdDivaLDwsu3fj\ns3/xJWSLyEUX45FDX0QylcKbr7+CXDYDjjfhU5/5LL78t98o2q2PxaJwudoLtr36o2BBC+/asROf\n6f3Lir/TWmEZIBGPY8c97WAZApZZGtszBIvbyFL3Pv95sZvOLlaYzKKzFevKV8LMzAx2795d5W+8\nPprRJmDJrpFFqfByb+cSpcTA6dy5c9i/fz/GxsYK1sPXM/AmkUhAFMWyJKLHxsYqsi0nybgcT2Fm\nNoVyYmzyx/B/7Hl01TF8MYev9hieYwnMrDpLb+LUiTrT4vulCTuiz+IzhDT9RdxMNKNNgLFdpfyh\nZAuvdeuDwaC+Y3d393rtLJvBwUEIggCv11vT83AsA+GeFmy3m/D+jTtIG0x4MQyDr578p4pn6dfq\n3CxDYOXVbqWVV2/LabfizBxLb8VR1oThbTlFUdDf3w+Xy4VwOFyyS19t6n3brtXCw/0xJz64eVeP\nhCsFwzDo+dyhqp6fEMBmYmFfDLppMXOw0YAbSpUxdPienh5dvrkZgmRqCccyeHBXG64lUpiOzUOq\nYcCdiSVwWHm0Wng9rJa21pRaY+jwnZ2dcDgcmJ6exrZt2+phU8PZ47TCZTPhw5t31JVvVcDEEjht\nJjitPBwWHlYTbbkp9YemuCqB1cTij/Y4MDObxtVEas3x9SwDOCwchHta4FxcGkuhNJqyU1y1t7ej\nra2tnrY1HEII9jit2O2wIJbM4vpcGvFkZsVsPiGAiWVgN3NwWDi0LTr49Y8y2O1svgAXytbFMMXV\nSy+9hM7OTvj9fjzyyCM4cuRIvWxrGghRAxu04IacJCMjLaW4MnEMTW1F2RAY9jP379+P4eFhDAwM\nQBRFo+JbAm5xcQqFstEo66ptb28HIQTnzp2rtT0UCqWGGEpNRSIR3L59G6dOndqSCSwplM2EodTU\ngQMH0NXVBUVRqNQUhbLBWTW09qmnnkI4HAYhBF1dXfjWt75VT9soFEqVKenwJ0+ehN/v14Uopqam\nDFNcUSiU5qbkGH656szy9xQKZeNR0uGLhdFuldBaCmWzUrJL7/P50NnZqb/XEmBcunSpLoZRKJTq\nU9LhR0ZGVqS6ovfhKZSNTckufbG8dpXkuqNQKM0DXcJVIUo2AyWXAxQZUABwHBizpdFmUSirQh1+\nDUh3ZiHHY5ASUSip5MoChAGxWMBYW8DYW6HMz0ORJBCWrn2nNAfU4ctATiWR/f0k5Nn46gUVGUoq\nCSmVhBS7BcRiSMdvgWmxg2ltA2NvA9PmBOHoz05pDIZX3tNPP61H2M3NzW2pNfGKJCF37ffIXZ9R\nu+4VHUSGfHcO8t05fROx2NQKoM2hVgIWumaeUh+aNollo5FTSWQuXSzedV8nSjoJKZ2EdOu6uoHj\nwdhb1UeL+iA8X/XzUihlJbHs6emBoiibPomlhhS9hczUh4BUvsDkushlISdikBMxfRMxmUCsdnU4\nYGsBsVpBLLaiklYUSrkYOrzT6cSRI0cwOzsLv99ftRMHg8Gyc85riTdqXeEokoTslSlIN2Zqep6y\nbMlkoGRikGeXKgEQAmK2gJit6uSg2aK+N5lBzGYQvjy5IcrWxbC5mJ2dxZkzZ3D27NlVNeKWE4lE\n0NnZiUgkAkAVltAqjEAgoOvQDQ4OIhgM6uUAVXEmFAphcHAQgOromrJsrZDuzGLhtxFDZ5dlGcNn\nRyHLNcxhXQpFgZJOQZ6NQboxg+xlEZlLF7Hw3jmkI79B6u23kD7/NhYunkfm0vvITP8O2WuXkbt5\nHVI8CvnOHOR0Gkq9ei6UpsOwhZ+cXJJFikajZR9Y04jTZKJ8Ph9iMbW10iqOfKFIv9+vl3U6nRAE\nQZeSBgCXywVRFKvayj/88MOALGP0O/8XuRszgIFcsSzL+Jvnvo5XfvUWXv/X3+Dbz311zcoz/+3v\nT+DS5Sv6+wfu24t/7D9WifkrUWQoC2koC+nVy8ViSH00DbAcCMeBcLz+GiwLwmrPrPrMqM9gGPU1\nwwAMC8IyS683WU4/TYb5+9//fmMNqTKGDu/1etHX1wdCCAYGBqpyUk0Jdnx8XFeYMcqXp7XylcpX\nayhSDko2C3k2Djl5Fx/dvIXAt79tuJ+sKPiXH/0E59//ADlJwktjb2ByZgZf+vznSiawLCbL/Otz\n7+B3eQ5/PRrDP//s9XV9p7VSG7loXfR96bWmpErytqkblkm6qp8lk0m0tLQsbdd/1yLl9af8c2Dp\n+MVeF3u/dKAC6daZa9ewe9cutUeVThUts/KcpWwocq4ybKtFJWqoHhsMBvVW1e/3r3uWXhRFuFwu\n/X0ikYDT6TQcLqy3W6/kcpDvzEJJp6BkM5DmEkCu/K7txHsX8e4Hl5Bb7A7nJAnvfnAJE+9dRPcf\nPlixXZsHRe0hFflJy9PmJUAqBVkrTQorg6Viyx1cLacUqQQKXuedZtmLguPqyBKUzAIQuwWJY1aW\nWSE6v5qNRSi3glpxbLWMEruNnIlbeR6DSsIwxZXb7S5Qj10LWksOoKA7rrXm3d3diMVicDqdBWWL\nEYvF1tWdJxwHdtuSkisPgLG3YrfZjL/5wuegZFbXk/urRx9Bej6JH77xJjLZHEw8hy/8h0/hH499\neRX12FhB5QYAP/jpqwUt/M52F/7LY5+p+HtVQjG7dAhRu+hcXreeKezOF25jCrv7hKjdfEJACKO/\nBsMAICAMAQiz1AMgRL/z0ExKrWdeOwsAILvvA98kNuVDMhK4e3aseb+aqscODAzoKrBOp3OF7LPP\n50MgEIDT6UR/fz9EUcTQ0BAGBgYQCoUQiUT0iiISiZQlG70W3vzVrwCorX/295OQbt8oWZZhGHz7\nua8CAF751Vv47J/+cUVj+Afu27vq+6pCGPX2Hm8CeB6EUx8wt4DfvUd1ai5/HL84fqfgzTffBKBW\nQpuJmgbeCIJQVAVWa82dTueKzw8fPgxArQzyx+vrbeFXg3AcTJ2fgOS6B9npSyVbe83pR0bH0Ntz\ncM3ODqB6E3SLcfvEZFFvyeXdmiO8SX1fIoSXcGZwO3ZWxw7KhqIh6rE+nw+hUEg/rhGiKJZddj2w\n29rB2NuQmfz3knHzDMPg8H/sqbktGoQ3gdhaQCxWNQDHYgVjsYKYzHWzgbJ5WNXhz58/j56eHj1+\n/tCh6mmil3LgYt32ekb4EZ6H6RN/qMbQz1wxvFVX1XObzCAt9sXw2sUIO+rYlCpS0uFPnDiByclJ\nxONx9Pf346GHHqqnXQ2FEAL+Y/vAtDqQFT8wnNCr7CQMGLtdXTxjbwPT2kYj5Sg1p6TDC4KAY8fU\n8eZLL720pRxeg3VsA/NH3chdu4zc9WuVr5gDAJZTV8i1OtTnllYaF0+pOyUdPj9Drfb6/PnzW87x\nCcuCv68D7PZ7kZu5AikRA3JZ4/0sNmAbA36foDq5raUO1lIoq2OYtVZRFP2+7VbOWstYbTB1fgKK\nokC5e0cN3JFyUGQZUBQQnlcn1MxWEKtNvdU1MwPu3ua7h0vZutCstWuEEALSqo65KZSNBs1aS6Fs\nIeisUQ357ne/22gTikLtKp9mtAmo3K4t4fCvvPJKxZ+vZ9/vfe97DTmv0b4b0S76P1ybXaWgDm/w\n+Xr2NaJRF4sRzWgX/R9WB6IodQwlqyFjY2ONNoFCaSq0RW/5bBqHp1AoxmyJLn29WJ6Lr1TOPgql\nUbDPPffcc402Yj0kEgm88MIL+Oijj5BOp7Fr166G2WKxqNpyv/71r9HT04MXXngBPp8Pbrcb3/zm\nN9HTU79VdvkkEgm89dZbCAaDSKVScLlcDf/NNJtCoRDi8XhT2JTP4OAgHnzwQaTT6aawS1sxGolE\n8OCDD4IQUpldygZnYGBAicfjiqIoyvHjxxtsjaJMTk7qdni9Xn17/ut6MzQ0pExOTiqKoigej6cp\nfrPR0VFldHRUmZiYUAYGBprCJo14PK54vV5lcnKyaeyanJzU7VCUyq/7Dd+lHx8f1xNqGCXCbARa\nrr61pPiuNj6fT88aJAhCU/xmHo8HgiBgaGgIx48fbwqbNMLhsJ7dqZnsOnPmDAKBACKRSMV2bXiH\nB5rDqYqh5ewDYJizrx4MDw/rmYeb4TcTBAF+vx+9vb1NY1MkEtFzOGo0g12CIOhZoIaGhiq2a8M7\nfLM5VX4uPp/Ph2AwiEAggP7+/obaFQwG0d/fj1gs1hS/md/v1/MViqLYFDYBamsZDocxPj6OUCjU\nNHYFAgHdsdfzP9zwt+USiYSeCLOrq6vqiS43A8FgEH6/XxcG6e/vb/hvFolEEIvF9OSkXV1dDbdJ\nI5FIoLet8l4rAAACoklEQVS3F729vejr62sKu0RR1B8ejwcul6siuza8w1MolPLZ8F16CoVSPtTh\nKZQtBHV4CmULQR2e0jREIpGyQpD9fn/RW1HBYLDmsuIbHerwm4BIJILOzk6EQiF9Rr4axwwEAus6\nhramoFwEQShL2ai7uxtOpxORSATBYBDBYBCAqnTcbLEYzQZ1+E2A2+2G2+2Gx+OB1+vVb9+s95jr\nkeZOJBKIRqPrssGI4eFheL1ejI6O1vQ8mwmqHLhJ0FbqiaKInp4ePaBF6+J6PB4kEgmEw2EkEgld\n7mt4eBhHjx7F6Ogoenp6MDo6WiDmqbW6R48eRSQSwfHjx/XjalLfpSqGzs5O/T57MBg0PE5fXx8A\nrLC7mPJQMBhEe3s7gsFgw0NeNxK0hd8kOJ1OeDwe+Hw+TE5OIhKJwO/3o6urS3d+t9uty0Rr2n4u\nlwsej0cP6NCctKurC9FotKDM5OQkgMLw02KSYVq0msvl0rvoaznOcruLMT4+rlc0jVqFuBGhLfwm\nRBvPA9Cj6xKJBPx+Pw4fPgyn06mPdbWwzJJ68VgZuplfgRRrfROJhO7A+WNyo+Pkj7/z7V7NLq1H\nQikP6vCbAG12W+vSj46O4vTp0/B6vQgEArpDdXZ2QhRFxGIxTExMgOM4Pe4/HA5DFEVMTEwAUOO1\nte3LyzidTgwPD0MQBLhcroKwzlAohJGRkYIFHoODg3C73YbH0c45MDBQYHexWPGenh4Eg0H09vY2\nxRqKDUNVF+1StgT5a7F9Pl/djzMyMlLRZxRFoS08Zc0IgoBwOAwAOHDgQMOPQykf6vCUNeP1eht6\nnPHxcXg8nhVd+WAwSLv3BtDVchTKFoLelqNQthDU4SmULQR1eAplC0EdnkLZQlCHp1C2ENThKZQt\nxP8HtV9n+WXRFwcAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "distributed_vt_plot = beaming_violins(o2a_traces_u, priors)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Assuming an event must originate in an MWEG" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def number_mweg(volume):\n", + " \"\"\"\n", + " Calculates the number of MWEGs in a volume, given in units of Mpc^3\n", + " \"\"\"\n", + " return volume * (0.0116) " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "42769.35938856692" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "number_mweg(total_v)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import theano.tensor as T\n", + "from pymc3 import DensityDist, Uniform, Normal\n", + "from pymc3 import Model\n", + "from pymc3 import distributions\n", + "\n", + "def grb_model_sigma(number_events, background_rate, \n", + " observation_time, volume, sigma_volume, grb_rate,\n", + " efficiency_prior = \"uniform\"):\n", + " with Model() as model:\n", + " signal_rate = pm.DensityDist('signal_rate', \n", + " logp=lambda value: log_signal_rate(value, number_events, background_rate, observation_time),\n", + " testval=50)\n", + "\n", + " volume = pm.Normal(\"volume\", volume, sigma_volume)\n", + " \n", + " n_galaxy = number_mweg(volume)\n", + " \n", + " cbc_rate = pm.Deterministic('cbc_rate', signal_rate / n_galaxy)\n", + " \n", + " grb_rate = (grb_rate / number_mweg(1e9)) #/ n_galaxy\n", + " \n", + " # Allow the efficiency prior to be switched-out\n", + " if efficiency_prior == \"uniform\":\n", + " efficiency = pm.Uniform('efficiency', 0,1)\n", + " elif efficiency_prior == \"jeffreys\":\n", + " efficiency = pm.Beta('efficiency', 0.5, 0.5)\n", + " elif isinstance(efficiency_prior, float):\n", + " efficiency = efficiency_prior\n", + " \n", + " def cosangle(cbc_rate, efficiency, grb_rate):\n", + " return T.switch((grb_rate >= cbc_rate*efficiency), -np.Inf, \n", + " (1.0 - ((grb_rate/(cbc_rate*efficiency)))))\n", + " \n", + " costheta = pm.Deterministic('cos_angle', cosangle(cbc_rate, efficiency, grb_rate)\n", + " \n", + " )\n", + "\n", + " angle = pm.Deterministic(\"angle\", theano.tensor.arccos(costheta))\n", + " \n", + " return model\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# Actual O2\n", + "priors = [\"uniform\", \"jeffreys\", 1.0, 0.5]\n", + "number_events = 1 # There were no BNS detections in O1\n", + "background_rate = 0.01 # We take the FAR to be 1/100 yr\n", + "grb_rate = 10.0 \n", + "o2a_models_mweg = []\n", + "for prior in priors:\n", + " o2a_models_mweg.append( grb_model_sigma(number_events, background_rate, total_time, total_v, sigma_vt/total_time, grb_rate, prior))" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 100000/100000 [00:18<00:00, 5272.37it/s]\n", + "100%|██████████| 100000/100000 [00:18<00:00, 5399.02it/s]\n", + "100%|██████████| 100000/100000 [00:12<00:00, 7945.84it/s]\n", + "100%|██████████| 100000/100000 [00:12<00:00, 7861.60it/s]\n" + ] + } + ], + "source": [ + "samples = 100000\n", + "o2a_traces_mweg = []\n", + "for model in o2a_models_mweg:\n", + " with model:\n", + " step = pm.Metropolis()\n", + " o2a_traces_mweg.append(pm.sample(samples, step))" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import pymc3" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "|\t\t\t| Lower\t| MAP\t| Median\t| Upper\t|\n", + "|----------|\n", + "| U(0,1) \t| 3.69\t| 8.64\t| 11.66\t| 40.01\t|\n", + "| Jeffreys \t| 3.51\t| 7.53\t| 11.85\t| 51.54\t|\n", + "| $\\delta(1)$ \t| 3.43\t| 5.42\t| 7.66\t| 16.82\t|\n", + "| $\\delta(0.5)$ \t| 4.79\t| 8.37\t| 10.88\t| 23.83\t|\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAChCAYAAADju9JOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWtwW+eZ3//vueBCgsQxaPkiyV4RdDKz47aRIHL6Zber\njUA3bjbZrAWS02Z2m2llyG6nH9JGgll/sDM72TUppZt6Zzci5c10unHHFGEnjmNbCsHE3njTbAhC\nshXbdRQe0rpQlEiAAG+4nnP6AThHAAnwACBuJN/fDAbnHLznPQ9B/N/7+zxEURQFFAplV8DU2wAK\nhVI7qOAplF0EFTyFsouggqdQdhFU8BTKLoIKnkLZRXD1NqBSjI+P19sECqWhOHr06IZrO0bwQP4/\nsJ7Mzs5i79699TZjA9Su4mlEmwB9uwpVgLRJT6HsIqjgKZRdhK7gl5aWamEHhUKpAbqCP378OC5f\nvlwLW7YtsiwjFoshHo9DkqR6m0OhFER30O78+fMAgHPnzkEURfT19eHgwYNVN6zRSSaTmJ+fRzQa\nRSKRQPYeJJZlYTQasby8DIvFArPZDJ7n62gthZJGV/BPPfUUFhcX0dXVBbfbDZvNhtdeew1PPPFE\nLexrOCRJQjAYRCgUgizLBdOsra0hEongxo0bAACj0Yjm5mbtxTB0+IRSe3QFf/jwYTz55JM516am\npqpmUCOztraGmzdvIplMlnxvPB5HPB5HKBQCwzBobm5Ga2srmpubwXE7anaU0sDo/tI6OzuxtLSE\n1tZWzMzM4MCBAzh58mQtbGsoIpEIbt26VbBWLwVZlrG8vIzl5WUQQmAymWCxWNDc3AyTyURrf0rV\n0BX85OQkDh06BAAIBAI4cOBAtW1qKBRFwfz8PBYWFqqWfzQaRTQaxfz8PBiGgdFoRFNTk/ZiWbYq\nz6bsPnSrksXFRczMzGBpaQm/+tWvKvpwr9dbdFpRFCGKYkWfr0cqlcK1a9eqJvZ8yLKMaDSKYDCI\n69ev45NPPoEoipibm8PS0lJZ3QkKRUVX8C6XCy+88AKefPJJPPXUU2U/SBRFDA8PayIfHh6Gy+VC\nOBzG4OAgvF4vAoFATvrDhw/jxIkTEEURdrsdPp+v7OeXyurqKkRRxOrq6qbpZFnGm2++WZGmfiFi\nsRhCoRBu3LiBq1ev4je/+Q0+/fRTzM3NIRQKYWVlBfF4vKo2UHYGuk369vZ2nD17FgAwMzOzpYfZ\nbDa4XC4AQDgcBpAWvtvthiAI8Hg8cDgcWvrx8XEIgpBzvyr+anDkyBHIsoxXXnkFkUgEeu7+ZFnG\nM888g/Hxcbz77rt44YUX8va/n3vuuZzv7sCBA/jmN79Ztp2pVAqpVCpvYcTzPHieB8dx4HkeLMuC\n4zhwHKcdl1MwHDlyBADwzjvvlG03pf7oCv706dMQRRGKomBychITExNlPchms2FsbAxOpxOCIGhC\nnpiYwKlTpwBgQ5NdXQPQ2dkJh8Oh1fJut7ssG/KRTCaRSCQQj8cRi8Vw+/ZtfO9739O9T5ZlvPrq\nq/jwww8hSRIuXLiA69ev49ixY5roo9EozGYz/H5/juAXFhbw2muvVexvKJVoNIqmpiYQQgq+AOQc\n37hxA/fffz+CwSAIIWAYZtP78+WV712FEAJZlvMWRtn5ULZGUdNy6qj89PR0WQ/x+Xyw2+3o7u6G\n3++H3W6HzWbTPg+HwxAEQav1AcBut2vCPnHiBIaGhirerI/FYlhbW0MymSy5SXzlyhV8/PHH2so6\nSZLw8ccf48qVK/jc5z5XMRurhaIoBVsw+cQqyzIkSdJmFtRXMcJX88jOO/tZKmtra1heXs772WaF\nx/rjfH/P+uN8BU6+70itEDZLW+jezZ5fLur9hQpHPXQFf/bsWYyOjsJqteLSpUu4ePFiyQ8ZGxvD\nwMAAfD4fnE4ngLu1eVdXF0KhUE6tD6Sb+r29vRAEAaFQCAAQCoUq2pw3mUwwmUzaudlsxt69e9HX\n16c7OPaVr3wF0WgUP/nJT5BMJsHzPB577DE899xzWg0fCoVgs9nw+uuv59Tw9957b9UXLjEMs6Ep\nz7IsWJbFwsIC9u7dq52rwlVf+bhw4QIAVHWWZnV1FVartWr5lwPP8zAYDPU2YwOb/a82Q1fwfX19\nmkj9fn/plmXyUAfrskUNAG63G8PDwxAEAf39/RBFEV6vFy6XC36/H6IoYmBgAEB6WjC7j19p3n33\nXQDp0nN+fh7BYLBgWoZh8MILLwBIjzUcPXq0YB9+vUgqJRr1x6i+q8dqH74Q8XgcLS0tJT2L9t13\nBrqCdzgcmJycRGdnJzo6Osp6iMPh2CBUVfiCIGh9eBX1fH1tXukavhAMw+D+++9HU1MTZmdnC26I\nUUX/9ttv4/HHHy9Y4m5lgE6F53k0NTXBbDbDZDLBaDTS+XlKyei2CVTPGa2trRWdB3e73SX1x0VR\n1FoataKlpQXt7e05zf71MAyDL37xixVfHWcwGCAIAvbu3YtHHnkEn/nMZ7Bv3z7YbDa6GIdSNro1\nvKIosNlsWFpawuTkJD7/+c9X7OGlCLgWNXs+DAYDDhw4gNnZ2ar6BjCZTGhubobZbEZTUxNdX0+p\nCrq/KqfTCY/HA0KI1pfebTAMg/3792N+fh7z8/MVy9NisWhr6On2WUotKNgOfemllwCkl7/a7Xa0\nt7fD4/HUzLBGZM+ePXj44YdhNBrLup/neQiCgIceegif/exnsX//fgiCQMVOqRkFa/jDhw8DSA+4\ndXZ2Aih/lH4nodbIkUgE8/Pzm07f8TwPq9UKs9kMi8XSkNM7lN1FQcGrO+Ta2tpACEFra2vZo/Q7\nDUIIBEGA1WrVVunF43EQQrQpMaPRiDt37jSki2PK7kW3D+/3+zWXVrtxe+xmEEJgNBrLbuJTKLWm\nrttjKRRKbdGt4dXtsYuLi7t2lJ5C2SnoCp4Qoo3O0yY9hbK9KdikP3PmDKanp3Hq1CmMj4/D5/Nh\nbGyslrZRKJQKU7CGdzgcaG9vx4kTJ7QgjVt1gEGhUOpLQcGfP38e58+fx/T0NLxe75YdYFAolPpT\nUPA9PT04evSotvUTKN8BBoVCaQwKCj4QCODSpUuYmJjApUuXoCgK/H4/RkZGamkfhUKpIAUFf/Lk\nSUQiERw6dEhbWluvHWsUCqUybLrwxmq1wm63Y3JyEoQQbX09hULZntTNAQaFQqk9dXWA0eisriQR\nXkwACkAYgGEIzE0cmpo58DyN/0bZflAHGHlIJmQEF2KIRTf6sotFJSwG4zAYGbRaDWi2cNRnOmXb\nUFTkGTUgxG4gGk3hzq0odILOIBGXsXAnhvAigVUwwtJChU9pfKjjtCxiMQnzc/pizyaVVBCcTwvf\n0sKjpZUHx9HmPqUxKemXuZOX1sbjEu7MRVFuPEYppSCymMDNa6u4fWsNqytJ3dh0FEqtKSq23NTU\nFADs2KW1ibiEO7eikKWtC1RRgOiahOiahHA4BbMxhqZmDkYTS5v8lLpTk9hypaJGnslGnRKs9OKf\nWCxTs1dA7OuRZWApksRSJAmWJTA3sTCaWBiMLAwGhhYAlJqj26Q/e/Ysnn76aTzzzDNbig9fiM3i\nxvt8PgwODgJAVeLDR9dSuHNrbYPYZVnG6z86X9F465KkYGU5heB8HLdurOH6zApu3VhFcD6G5aUE\nYjEJqZRMuwGUqlJUbLnu7m5tLX01yBc3XhAE2O32nD34lYoPn0hI+Pzn/xCplIJXXn4r5zNZlvH1\n//Ykxnxv4qc/u4i/+va5DVFl+v/7f8H0zG+18/YDj+Av/+KvS7JBloF4XEY8nluoEAJwHAOWJWA5\nApYlYJj0izAETGY9QPoc6O4+CkJo7DdKcegKXhAEHD9+HJFIpCp+6QvFjc/HVuPDJ5MygvPp+fVU\nUsGd+TmMjP4v7XNZlnH+/N/jyq/fhySl8OZbP8C1ayJ6e/80R/T/NPEepqentPP1+aio8eGryfXr\nN3HfngcwM7WcUxiAAAzJvDMEhKghl4GlsASzMZ6+lvUZgMzxxueo96ZPAKImhnqcedeOsy/kz3M9\nkqRASsnFJS4h363cs1lY7XpSrl26go9EIto8/GuvvVa6ZZugFzd+PVtt1vM8gwf2NiERl8DxBOv/\n7+9/EMCHH12BJKUAAJKUwocfXcH7HwRw6GBn2c8tBQJotTdD0u+q2O6+A8gccxyByczid+yWoscE\n4kkW97Q1nqfddKumsaY018ezbxTKtUtX8OoIPYBNwyeXw2Zx4/NRqeixBiMLjmOwb/9e/Nu+r2lT\ncT3H/gzRtSjevvA6kskEeN6Ax7/wx/iLP38xp4b/wQ9Gcmr4+/Y8gL6er+W1t1ABRghgMDDaAB5v\nYMDx6aZ8Kf/I8+dfyeTXeD9KSuNRlNfa3t7eqiyt3SxuPJBuAQQCAa3fXsn48GqfN3uUnmEY/NW3\nzwEAxnxvotv5xbx9+PYDj2x6ng+GJTCbWZjMlR2lp313SinUdWntZnHjgXRI6ez+ejXiw5tMLB54\n0Izbt6KQskT/xo+9+NIfufKGgS52gI5hgFYrT+fhKQ1DY3WYUDhufDXjwxuMLO570AxV2wzD4I+/\n3FtWzHdCAHMTiz33m7DnAQ62e00wmek6e0pj0JBr6fMJu9redoxGFvc92JSely9j+p3lCCwWHpZW\nXts6G1miIqc0Fg0p+HphMrHYc78Zd0rYQMPxdLccZftAN8+sw9zEYe/+ZpjM7KbpDAYGbXtM2PdQ\nM1paeSp2yraAbp7JA29Iz9evLCcRWUwAyPJ4Y2bRZOGpxxvKtqQhN880CpYWHpYWvt5mUCgVQ1fw\nZ8+exejoKKxWKy5duoSLFy/Wwi4KhVIFGmLzDIVCqQ26HVF184y62o5CoWxf6rp5hkKh1BbdGr6a\nm2coFEptqevmGQqFUlsK1vAvvfQSgLR/Obvdjvb29qo4wKBQKLWjYA2vBo50OBxa9Fg6Sk+hbG8K\n1vCHDh0CALS1tcFqtWJxcZGO0lMo25yCNfz09DS8Xi/8fj+6urqgKApCodCuCiZJoew0Cgq+vb0d\nLpcLdrsdDocDNpsNVqu1lrZRKJQKs+kofXt7O9rb27Xzy5cv4+DBg1U3qpGJJSVEkxJiSQnxlAwQ\ngGMY8CyBiWPRbGTBleE4g0KpBbrTcr29vWhra0MwGMT09PSu2C23nnhKwvxyHLdX4liNbwwhvR4z\nz6DFxCO2koQ1kUKzgbodoDQGur/Ec+fOaU358fHxqhvUSMRTEmaCa7izEi8pomw0KSOajCMUjmHp\nWhg8S2A187CaeFjNPCxGWgBQ6kNRgs9eYXf06NGqGtQIpGQZNxajuBmJQqpAtKmkpGBhJYGFlfTe\neo4laDVyaDHxsJo4WIwcOJZ2AyjVpyivtU8++eSuGbC7sxzHdHA13T+vEilJQWgtidBaUrumdgMs\nBhbNRg7NBg6GBgvKQNn+6Ar+2LFj2qaZJ554ouoG1YtoQsJvF1awmCXCmj4/0w24k3WNZwmaDBya\nDSyaDCzMfPpl5GjkWUp5FOXiShAE2Gw2nDlzBt/4xjcKphVFER6PB11dXTh16tSGzwcHB+FwODA0\nNISuri5t2q+eJFIyri2uYW4pBrnBQoglJQWRaBKRaG4hxBDAxLMwcQxMmQLAyLEwcAyMLAOeI3Sm\ngJKXolxcqYttfvrTn26aVo0Tl0/EalRYp9OJsbGxvAVCLYmnJNyKxMrqp8uyjPEf/wBH/+hPyvJd\nv1VkBVhLSFhLSADyt0hYBuBZJv1iCHiWAcemC4LQSgLcchwcQ8AyBBxLwBKindPWw84lr+DVVXYA\nMDExgUAgUPJKO1EUtYASTqcToihiYmJCCx+lhpBSw0F7PJ6c9D6fD5OTkxgaGkJPTw+++tWvwmKx\nIBAIlF1YyIqCO8tx3FmO498/8W8AAP/z5ddLy0OW8a1v/Ce853sbv3xnDM+e+duSRH/62f+K6zN3\ntxw/dKADJ7/1P0qyoRgkGZBkGbHkxtIsFI5jmVkueC/LELAk/c4wBBzJhKwm6QKByUSmTb/Sce2y\n35lMoEuGpAN2Zn9OgLvHBGBwN0BmUpKRkmQt8q16L6U4jhw5AmDz8GN5BR8Oh+FyuWCz2coOzezx\neNDf349QKARRFNHZ2YlAIACn04nR0VE4HA7Y7XZMTExgYGAAPT09Oendbjd6enoAAP39/fD5fAgG\ng+jr6yvLHgBYjafwmzsr2nlo/jYujP590fcrsowfjf4ffPLr9yFJEn721g8xe20aX+75dyB5RJ8v\nXPT7v/pHXJ+5GzAzdGeuJBsqQS3CWJdDPrtUvavhqbMiVmd9TrSQ1Vj3efa9WH8tD+vLl7W1NTQ1\nNWWFwt4YdTj72eufsll5RQqeFLZPZXV1Fc3NzTnXrt+4ifseeBA3wtGC9+UVvLpxphQGBwfhcrm0\npjsAbVlu9rX1tLW1FUzf1dWF4eFhuN1uLeacx+PB6OhoyfYBd+Oklxvu+6MPLuPqR1cgSenFN5Ik\n4epHV/DRB5fx6MHKBLncqZCsmPEE2THl74pCYglMHJMrbFVI64SWUxCsExZZd5BP6uvTbEyXRuYY\nmHlWK1jWp8v73KwP8hYOhR6m/5GGlLHr7j3pVhjPErQ1GQreV7AP/9JLL+H48eM4ffo0QqEQFEXB\n9PQ0RkZG8qZ3uVwQRRHBYBBOpxN2ux3Dw8NaTR4IBDAxMQFRFOH3+yGKonYtHA5jYGAgJ70gCHC5\nXBgaGgKQ3pfvcDjQ3d1dxNeRn2YDh3/5OzbcWYmDZQhse+7HF3r+tOj7Hzv2VaxFo3jnwhtIJRPg\neAP+4Atfwtf//Nt5m/WhUBA2W1vOtbd+6M2p4W33PVCSDZVAtYtkmu1aX54QcCyjNdvZzHWWZDfn\n736W3YRniPqjy7yva+oXw+zsLPbu3Vvlv740GtEmIL9do5nQ4WZD4SAqRFHy13eXLl3CoUOHMD4+\nnrMfvpYLb8LhMERRLCpE9Pj4eMm2LcWSmA6ubRgF34zsPvzvOR/ftA+fT/DV7MOzDIGRSw/UGVgC\nA7du0I5JC3r+9hwe2r+34UbyG1FcjWgToG9XIT0UrOHVZr3X69Vu7Orq2qqdRTM4OAi73Q6Xy1W1\nZ7SaeHxunxWh1QSmQ6tFrZNnGAbPnvnbskfptyJulklPx5l5FiaOhYlnYOJYGHkGRo4pWsBLJaSl\n7Cx0p+UURUF/fz9sNhv8fn/BJn2lqeW0na3ZgHuaeMxGYpgJrUHSmZBnGAbdXz5WNXtYBlkLbjg0\n8emFNyZ+83h3FIoeuoLv7u7WwjfXe5FMNSGEYJ9gxr0WA8SFVcxn1r1XG5YhsBg5WIwsWow8LEY2\nM0hEp6MolUdX8B0dHbBarZiZmcE999xTC5vqipFj8bsPtGJfLInphVVEYqkK58/AaubRauLQauLR\nbKDiptQO6uKqAK0mHp/bL2BhJY7r4SiWyxS+iWe0bbGCmafNckpdKdrFVVtbG1pbW2tpW0Nwr8WI\ney1GRBMSbq/EEFxNIJaU8i7H5VgCi4FDc6Z5vmaM4cBDttobTaEUQNfF1auvvoqOjg54PB489thj\nOH78eK1sayjMBhYHbM04YEuvbkqkZMRSEpjMGnR1yiub2WU6Ek5pLHT78IcOHcLIyAgGBgYgiqJe\n8l2DgWPofnXKtqOoX2xbWxsIIbh06VK17aFQKFVEN9RUIBDAwsICzp49uysdWFIoOwndUFOHDx9G\nZ2cnFEWhoaYolG3Opktrn3rqKfj9fhBC0NnZie9+97u1tI1CoVSYgoI/c+YMPB6PFohienpa18UV\nhUJpbAr24ddHnVl/TqFQth8FBZ9vGe1uWFpLoexkCjbp3W43Ojo6tHPVAcbVq1drYhiFQqk8BQU/\nOjq6wdUVnYenULY3BZv0+fzalePrjkKhNA40quEWUKQUlFjsrldMhoAYTCAc/VopjQn9ZZaItBiE\ntBiEvLIEJRbN6wKX8AYQkxnKahSSyQCmpRWEL+xJlEKpFVTwRaBIEqSF20jN3UyLXC99MgElmQBC\nISTktOccYjKDabGCabWCbRVADMZqm02hbEBX8E8//bS2wm5paWnX7YlPLdxB6roIJbE1l1dKLAop\nFoU0P4ckMgVAqwC2VQDTaqUtAEpNaFgnlvVGjkWRnPkt5MhiVfLXCoA7twAAxNQExiqAsbSmXyZT\nVZ5L2d0U5cSyu7sbiqLsaCeWKoosI3XrOlI3rwNK9WLEb3hubA1SbA3S7VkAmXGAZguYJguYpmaQ\nZguI0UT931G2hK7gBUHA8ePHEYlE4PF4KvZgr9dbtM951fFGtQsceXkJienfQImuVfU5xaAkE1DC\nIcjh0N2LhAExGtPdAZMZxGBMFwJGU/qY5+tnMGVboOsAIxKJ4Pz587h48eKmMeLWEwgE0NHRgUAg\nACAdWEItMIaHh7U4dIODg/B6vVo6IB1xxufzYXBwEEBa6Gpk2WogR9eQ+O3/Q/zj98sSuyzLGLk4\nBlmucotAkaHEopDDIaTmbiJ5TUTi6keI/zqAWOD/Ivqr9xB7fwLxjy4jcfVjJD+dQnL2OlLzc5DC\nIcirK5DjMSjVtpPSsOjW8FNTd8MiBYPBojNWY8SpYaLcbjdCoXRtpRYc2YEiPR6PllYQBNjtdi2U\nNADYbDaIoljRWv7I7/8+lEQcF/7622VHmJRlGf/x+W/hjX94Dxf+8Zf4u+ef3TQazX/+y9O4eu26\ndv6Zhx/C3/SfLOvZG8gUCLozCaEQonOfgrAcwPMgHJ8+5risdxZgWRCGBdjMOcMADAvCMnePd3AX\n48iRI0gkEvjFL35Rb1Mqhq7gXS4Xent7QQjBwMBARR6qRoKdmJjQIszo+ctTa/lyw1erSEthSPNz\nkCNhyKsruBUK4Xs/frusvGRFwf/+4Y9w+eNPkJIkvDr+M0zNzuLPvvJlMITkDX/880vv47dZgp8L\nhvC9Ny9s6W8qlcqGiyZZoWDXHQOZSLB50hCyId3a2lomBLKaHnc/z76WnZ/2lonWmhPOleS+r7++\nId/cPGdv3kRbmw3yhgI0Nyxt3mciz7Nz7C9sUzULUd3osV6vV6tVPR7PlkfpRVGEzXbXdXM4HIYg\nCLrdhUo16+XlJUjBhYoMyE1++BE++OQqUpnw0SlJwgefXMXkhx+h6589uuX8twcKoGS/501RBASI\nRSFDWR8UPnOcXyBkXUGg5L0361gvX2QJTpaAZDIzk5JP1LmXC4k9r4DzFUbrCyIUOCUESmgBKcNG\n+eoVFrourhwOR0702FJQa3IAOc1xtTbv6upCKBSCIAg5afMRCoUq0pzn9z0M7oG9kCNhEIMBD7bZ\n8B+++IWy8vra448htrqGH/zsHSSSKRh4Dn/yh3+Avzn5dTAMg1AolFO4AcArP34rp4Z/YAvPL5d8\ndgFI/+DYrOY8y919Z5hMsz7dtCdMpllP0u+EZGothsnU4ukXIQzAZGrzzGdEvS+TXv2RNlqk1pG3\nLiCRSIB/uPFmp0hCAnfvfSXfV9XosQMDA1oUWEEQNoR9drvdGB4ehiAI6O/vhyiKGBoawsDAAHw+\nHwKBgFZQBAKBosJGFwNhObC2e/HuL/8JSjKB1K0bSN2eBUoczGIYBn/3/LMAgDf+4T186V/9nm4f\n/jMPP7TpeUVhmLuj9xyfnurjecBkAb9vP0im/w6OA1H77xSNd955B7Ozs/U2o6JUdeGN3W7PGwVW\nrc0FQdjweV9fH4B0YZDdX69UDb8ewhvAP2wH9+B+JK/PQJqfK+l+VfSjY+Po6T6qGz66YgN0DJOZ\njjOBmExpMRuNILwRxGBIvwoImDCGsmoHyvanLtFj3W43fD6flq8eoigWnbZcCG+Awf5ZSPfeh+T0\n1aLWzKswDIO+f91dNbtIU3N68Y2pCcRkAmNuoktxKWWxqeAvX76M7u5ubf38sWOVi4leSMD5mu21\nXOHHtgpg/vlhpG5eQ+rWjZqutgPD3F1a22yhu+woFaeg4E+fPo2pqSksLi6iv78fBw8erKVddYUw\nDPiHDoDdcz+S01chLxW/4KgkGCa9gy6zi45pbkkPaFEoVaKg4O12O06eTPc3X3311V0leBXGZIbx\nd/8FpOA8krPXoKytbi1DwoCxWMC0CumXpZUKnFJTCgo+20Otenz58uVdKXy2bQ/Ytj2QlsJIzd2E\nHF4srqnPsmCs94BpaQVjsYKxtKSntiiUOqHrtVZRFG3edrd7rWUz+9cVWYayugJ5ZSl3FRYhYIwm\nEJMZxNwEBEMw7ttXP4MplHVQr7VlQBgGpKUVTMvmzkB28jpzyvaEeq2lUHYRdMSoinz/+9+vtwl5\noXYVTyPaBJRv164Q/BtvvFH251u59+WXX67Lc/Xu3Y520f9haXYVggpe5/Ot3KtHvX4sejSiXfR/\nWBmIopTp+aHBGB8fr7cJFEpDoW56y2bHCJ5CoeizK5r0tWK9L75CPvsolHrBPv/888/X24itEA6H\n8eKLL+LWrVuIxWJ48MEH62aLKeNL/uc//zm6u7vx4osvwu12w+Fw4Dvf+Q66u6uzo06PcDiM9957\nD16vF9FoFDabre7fmWqTz+fD4uJiQ9iUzeDgIB599FHEYrGGsEvdMRoIBPDoo4+CEFKeXco2Z2Bg\nQFlcXFQURVFOnTpVZ2sUZWpqSrPD5XJp17OPa83Q0JAyNTWlKIqiOJ3OhvjOxsbGlLGxMWVyclIZ\nGBhoCJtUFhcXFZfLpUxNTTWMXVNTU5odilL+737bN+knJiY0hxp6jjDrgeqrrxQX35XG7XZrXoPs\ndntDfGdOpxN2ux1DQ0M4depUQ9ik4vf7Ne9OjWTX+fPnMTw8jEAgULZd217wQGOIKh+qzz4Auj77\nasHIyIjmebgRvjO73Q6Px4Oenp6GsSkQCGg+HFUawS673a55gRoaGirbrm0v+EYTVbYvPrfbDa/X\ni+HhYfT399fVLq/Xi/7+foRCoYb4zjwej+avUBTFhrAJSNeWfr8fExMT8Pl8DWPX8PCwJuyt/A+3\n/bRcOBzWHGF2dnZWzNHlTsLr9cLj8WiBQfr7++v+nQUCAYRCIc05aWdnZ91tUgmHw+jp6UFPTw96\ne3sbwi6ueDf0AAACl0lEQVRRFLWX0+mEzWYry65tL3gKhVI8275JT6FQiocKnkLZRVDBUyi7CCp4\nSsMQCASKWoLs8XjyTkV5vd6qhhXfCVDB7wACgQA6Ojrg8/m0EflK5Dk8PLylPNQ9BcVit9uLimzU\n1dUFQRAQCATg9Xrh9XoBpCMdN9pajEaDCn4H4HA44HA44HQ64XK5tOmbrea5ldDc4XAYwWBwSzbo\nMTIyApfLhbGxsao+ZydBowfuENSdeqIooru7W1vQojZxnU4nwuEw/H4/wuGwFu5rZGQEJ06cwNjY\nGLq7uzE2NpYTzFOtdU+cOIFAIIBTp05p+aqhvgsVDB0dHdo8u9fr1c2nt7cXADbYnS/ykNfrRVtb\nG7xeb92XvG4naA2/QxAEAU6nE263G1NTUwgEAvB4POjs7NTE73A4tDDRamw/m80Gp9OpLehQRdrZ\n2YlgMJiTZmpqCkDu8tN8IcPU1Wo2m01ropeSz3q78zExMaEVNPXahbgdoTX8DkTtzwPQVteFw2F4\nPB709fVBEAStr6suy8wbLz7D+qWb2QVIvto3HA5rAs7uk+vlk93/zrZ7M7vUFgmlOKjgdwDq6Lba\npB8bG8O5c+fgcrkwPDysCaqjowOiKCIUCmFychIcx2nr/v1+P0RRxOTkJID0em31+vo0giBgZGQE\ndrsdNpstZ1mnz+fD6OhozgaPwcFBOBwO3XzUZw4MDOTYnW+teHd3N7xeL3p6ehpiD8W2oaKbdim7\nguy92G63u+b5jI6OlvUZRVFoDU8pGbvdDr/fDwA4fPhw3fOhFA8VPKVkXC5XXfOZmJiA0+nc0JT3\ner20ea8D3S1Hoewi6LQchbKLoIKnUHYRVPAUyi6CCp5C2UVQwVMouwgqeAplF/H/AY4EjSPA7IX9\nAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "mweg_viols = beaming_violins(o2a_traces_mweg, priors)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (grbeaming)", + "language": "python", + "name": "venv-grbeaming" + }, + "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" + }, + "latex_envs": { + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 0 + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}