ExtDataFig1b.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "from scipy.optimize import minimize\n",
    "from scipy import special\n",
    "from scipy.stats import poisson,norm\n",
    "from scipy.special import j_roots\n",
    "from scipy.special import beta as beta_fun    \n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# moment-based inference\n",
    "def MomentInference(vals, export_moments=False):\n",
    "    # code from Anton Larsson's R implementation\n",
    "    from scipy import stats # needs imports inside function when run in ipyparallel\n",
    "    import numpy as np\n",
    "    m1 = float(np.mean(vals))\n",
    "    m2 = float(sum(vals*(vals - 1))/len(vals))\n",
    "    m3 = float(sum(vals*(vals - 1)*(vals - 2))/len(vals))\n",
    "    \n",
    "    # sanity check on input (e.g. need at least on expression level)\n",
    "    if sum(vals) == 0: return np.nan\n",
    "    if m1 == 0: return np.nan\n",
    "    if m2 == 0: return np.nan\n",
    "    \n",
    "    r1=m1\n",
    "    r2=m2/m1\n",
    "    r3=m3/m2\n",
    "    \n",
    "    if (r1*r2-2*r1*r3 + r2*r3) == 0: return np.nan\n",
    "    if ((r1*r2 - 2*r1*r3 + r2*r3)*(r1-2*r2+r3)) == 0: return np.nan\n",
    "    if (r1 - 2*r2 + r3) == 0: return np.nan\n",
    "        \n",
    "    lambda_est = (2*r1*(r3-r2))/(r1*r2-2*r1*r3 + r2*r3)\n",
    "    mu_est = (2*(r3-r2)*(r1-r3)*(r2-r1))/((r1*r2 - 2*r1*r3 + r2*r3)*(r1-2*r2+r3))\n",
    "    v_est = (2*r1*r3 - r1*r2 - r2*r3)/(r1 - 2*r2 + r3)\n",
    "    \n",
    "    if export_moments:\n",
    "        return np.array([lambda_est, mu_est, v_est, r1, r2, r3])\n",
    "    \n",
    "    return np.array([lambda_est, mu_est, v_est])\n",
    "def MaximumLikelihood(vals, export_asymp_ci = False, fix = 0, export_fun = False):\n",
    "    from scipy.interpolate import interp1d\n",
    "    from scipy.optimize import minimize\n",
    "    from scipy import special\n",
    "    from scipy.stats import poisson,norm\n",
    "    from scipy.special import j_roots\n",
    "    from scipy.special import beta as beta_fun    \n",
    "    import numpy as np\n",
    "    if len(vals) == 0:\n",
    "        return np.array([np.nan, np.nan, np.nan])\n",
    "    def dBP(at, alpha, bet, lam):\n",
    "        at.shape = (len(at), 1)\n",
    "        np.repeat(at, 50, axis = 1)\n",
    "        def fun(at, m):\n",
    "            if(max(m) < 1e6):\n",
    "                return(poisson.pmf(at,m))\n",
    "            else:\n",
    "                return(norm.pdf(at,loc=m,scale=sqrt(m)))\n",
    "        \n",
    "        x,w = j_roots(50,alpha = bet - 1, beta = alpha - 1)\n",
    "        gs = np.sum(w*fun(at, m = lam*(1+x)/2), axis=1)\n",
    "        prob = 1/beta_fun(alpha, bet)*2**(-alpha-bet+1)*gs\n",
    "        return(prob)\n",
    "    def LogLikelihood(x, vals):\n",
    "        kon = x[0]\n",
    "        koff = x[1]\n",
    "        ksyn = x[2]\n",
    "        return(-np.sum(np.log( dBP(vals,kon,koff,ksyn) + 1e-10) ) )\n",
    "    x0 = MomentInference(vals)\n",
    "    if np.isnan(x0).any():\n",
    "        x0 = np.array([10,10,10])\n",
    "    bnds = ((1e-3,1e3),(1e-3,1e3), (1, 1e10))\n",
    "    vals_ = np.copy(vals) # Otherwise the structure is violated.\n",
    "    try:\n",
    "        ll = minimize(LogLikelihood, x0, args = (vals_), method='L-BFGS-B', bounds = bnds)\n",
    "    except:\n",
    "        if export_fun:\n",
    "            return np.array([np.nan,np.nan,np.nan]), np.nan\n",
    "        return np.array([np.nan,np.nan,np.nan])\n",
    "    #se = ll.hess_inv.todense().diagonal()\n",
    "    if export_fun:\n",
    "        return ll.x, ll.fun\n",
    "    estim = ll.x\n",
    "    return estim\n",
    "def ProfileLikelihood3(vals, alph = 0.05):\n",
    "    from scipy.interpolate import interp1d\n",
    "    from scipy.optimize import minimize, fsolve\n",
    "    from scipy import special\n",
    "    from scipy.stats import poisson,norm, chi2\n",
    "    from scipy.special import j_roots\n",
    "    from scipy.special import beta as beta_fun\n",
    "    import numpy as np\n",
    "    def dBP(at, alpha, bet, lam):\n",
    "        at.shape = (len(at), 1)\n",
    "        np.repeat(at, 50, axis = 1)\n",
    "        def fun(at, m):\n",
    "            if(max(m) < 1e6):\n",
    "                return(poisson.pmf(at,m))\n",
    "            else:\n",
    "                return(norm.pdf(at,loc=m,scale=sqrt(m)))\n",
    "        x,w = j_roots(50,alpha = bet - 1, beta = alpha - 1)\n",
    "        gs = np.sum(w*fun(at, m = lam*(1+x)/2), axis=1)\n",
    "        prob = 1/beta_fun(alpha, bet)*2**(-alpha-bet+1)*gs\n",
    "        return(prob)\n",
    "    def ll_kon(x, l, obs):\n",
    "        mu = x[0]\n",
    "        v = x[1]\n",
    "        return(-np.sum(np.log(dBP(obs,l,mu,v) + 1e-10)))\n",
    "    def ll_burst_size(x, burst_size, obs):\n",
    "        l = x[0]\n",
    "        v = x[1]\n",
    "        mu = v/burst_size\n",
    "        return(-np.sum(np.log(dBP(obs,l,mu,v) + 1e-10)))\n",
    "    def minll_kon(x, kon, obs):\n",
    "        bnds = ((1e-3, 1e3), (1,1e10))\n",
    "        res = minimize(ll_kon, x, args = (kon,obs), method='L-BFGS-B', bounds = bnds)\n",
    "        ll = res.fun\n",
    "        x0 = res.x\n",
    "        return ll, x0\n",
    "    def minll_burst_size(x, burst_size, obs):\n",
    "        bnds = ((1e-3, 1e3), (1,1e6))\n",
    "        res = minimize(ll_burst_size, x, args = (burst_size,obs), method='L-BFGS-B', bounds = bnds)\n",
    "        ll = res.fun\n",
    "        x0 = res.x\n",
    "        return ll, x0\n",
    "    def get_h(param):\n",
    "        x = np.log10(param)\n",
    "        h = 10**(x-2)\n",
    "        return(h)\n",
    "    \n",
    "    def kon_ll(vals,start_estim):\n",
    "        cutoff = chi2.ppf(1-alph, 1)/2\n",
    "        x0 = [start_estim[1], start_estim[2]]\n",
    "        bnds = ((1e-3, 1e3), (1,1e10))\n",
    "        N = 100\n",
    "        start = start_estim[0]\n",
    "        res = minimize(ll_kon, x0, args = (start,vals), method='L-BFGS-B', bounds = bnds)\n",
    "        h = get_h(start)\n",
    "        h = 5*h\n",
    "        ll_p = res.fun\n",
    "        ll_l = np.zeros(N)\n",
    "        kon_l = np.array([])\n",
    "        try:\n",
    "            for i in range(N):\n",
    "                kon = start - h*i\n",
    "                if kon <= 0:\n",
    "                    break\n",
    "                kon_l = np.append(kon_l, kon)\n",
    "                try:\n",
    "                    res = minimize(ll_kon, res.x, args = (kon,vals), method='L-BFGS-B', bounds = bnds)\n",
    "                except ValueError:\n",
    "                    break\n",
    "                ll_l[i] = res.fun\n",
    "                if (2*(ll_l[i] - min(ll_l[ll_l > 0])) > cutoff + 0.5) and (ll_l[i] > ll_l[i-1]):\n",
    "                    break\n",
    "        except ValueError:\n",
    "            return np.array([minimum,np.nan,np.nan]), np.nan\n",
    "        \n",
    "        ll_l = ll_l[:i+1]\n",
    "        ll_l = ll_l[::-1]\n",
    "        kon_l = kon_l[::-1] \n",
    "        res = minimize(ll_kon, x0, args = (start,vals), method='L-BFGS-B', bounds = bnds)\n",
    "        ll_u = np.zeros(N)\n",
    "        kon_u = np.array([])\n",
    "        try:\n",
    "            for j in range(N):\n",
    "                kon = start + h*j\n",
    "                kon_u = np.append(kon_u, kon)\n",
    "                try:\n",
    "                    res = minimize(ll_kon, res.x, args = (kon,vals), method='L-BFGS-B', bounds = bnds)\n",
    "                except ValueError:\n",
    "                    break\n",
    "                ll_u[j] = res.fun\n",
    "                if (2*(ll_u[j] - min(ll_u[ll_u > 0])) > cutoff + 0.5) and (ll_u[j] > ll_u[j-1]):\n",
    "                    break\n",
    "        except ValueError:\n",
    "            return np.array([minimum,np.nan,np.nan]), np.nan\n",
    "        \n",
    "        ll_u = ll_u[:j+1]\n",
    "        \n",
    "        ll = np.concatenate((ll_l[:-1],np.array([ll_p]),ll_u[1:]))\n",
    "        kon_indexed = np.concatenate((kon_l[:-1],np.array([start]),kon_u[1:])).squeeze()\n",
    "        \n",
    "        ll_ratio = 2*(ll - min(ll)).squeeze()\n",
    "        \n",
    "        minimum_idx = np.argmin(ll_ratio)\n",
    "        \n",
    "        ll_right_side = ll_ratio[minimum_idx:]\n",
    "        ll_left_side = ll_ratio[:minimum_idx]\n",
    "        \n",
    "        minimum = kon_indexed[minimum_idx]\n",
    "        kon_right_side = kon_indexed[minimum_idx:]\n",
    "        kon_left_side = kon_indexed[:minimum_idx]\n",
    "        try:\n",
    "            f_1 = interp1d(ll_left_side,kon_left_side, kind='cubic')\n",
    "            f_2 = interp1d(ll_right_side,kon_right_side , kind='cubic')\n",
    "            res = np.array([minimum, f_1(cutoff), f_2(cutoff)])\n",
    "            return res, kon_indexed, ll_ratio\n",
    "        except (ValueError,np.linalg.linalg.LinAlgError, TypeError):\n",
    "            return np.array([minimum,np.nan,np.nan]), kon_indexed, ll_ratio\n",
    "\n",
    "    def burst_size_ll(vals, start_estim):\n",
    "        cutoff = chi2.ppf(1-alph, 1)/2\n",
    "        x0 = [start_estim[0], start_estim[2]]\n",
    "        bnds = ((1e-3, 1e3), (1,1e10))\n",
    "        N = 100\n",
    "        start = start_estim[2]/start_estim[1]\n",
    "        res = minimize(ll_burst_size, x0, args = (start,vals), method='L-BFGS-B', bounds = bnds)\n",
    "        h = get_h(start)\n",
    "        h = 3*h\n",
    "        ll_p = res.fun\n",
    "        ll_l = np.zeros(N)\n",
    "        burst_size_l = np.array([])\n",
    "        try:\n",
    "            for i in range(N):\n",
    "                burst_size = start - h*i\n",
    "                if burst_size <= 0:\n",
    "                    break\n",
    "                burst_size_l = np.append(burst_size_l, burst_size)\n",
    "                try:\n",
    "                    res = minimize(ll_burst_size, res.x, args = (burst_size,vals), method='L-BFGS-B', bounds = bnds)\n",
    "                except ValueError:\n",
    "                    break\n",
    "                ll_l[i] = res.fun\n",
    "                if (2*(ll_l[i] - min(ll_l[ll_l > 0])) > cutoff + 0.5) and (ll_l[i] > ll_l[i-1]):\n",
    "                    break\n",
    "        except ValueError:\n",
    "            return np.array([minimum,np.nan,np.nan]), np.nan\n",
    "        ll_l = ll_l[:i+1]\n",
    "        ll_l = ll_l[::-1]\n",
    "        burst_size_l = burst_size_l[::-1] \n",
    "        res = minimize(ll_burst_size, x0, args = (start,vals), method='L-BFGS-B', bounds = bnds)\n",
    "        ll_u = np.zeros(N)\n",
    "        burst_size_u = np.array([])\n",
    "        try:\n",
    "            for j in range(N):\n",
    "                burst_size = start + h*j\n",
    "                burst_size_u = np.append(burst_size_u, burst_size)\n",
    "                try:\n",
    "                    res = minimize(ll_burst_size, res.x, args = (burst_size,vals), method='L-BFGS-B', bounds = bnds)\n",
    "                except ValueError:\n",
    "                    break\n",
    "                ll_u[j] = res.fun\n",
    "                if (2*(ll_u[j] - min(ll_u[ll_u > 0])) > cutoff + 0.5) and (ll_u[j] > ll_u[j-1]):\n",
    "                    break\n",
    "        except ValueError:\n",
    "            return np.array([minimum,np.nan,np.nan]), np.nan\n",
    "        \n",
    "        ll_u = ll_u[:j+1]\n",
    "        \n",
    "        ll = np.concatenate((ll_l[:-1],np.array([ll_p]),ll_u[1:]))\n",
    "        burst_size_indexed = np.concatenate((burst_size_l[:-1],np.array([start]),burst_size_u[1:])).squeeze()\n",
    "        \n",
    "        ll_ratio = 2*(ll - min(ll)).squeeze()\n",
    "        \n",
    "        minimum_idx = np.argmin(ll_ratio)\n",
    "        \n",
    "        ll_right_side = ll_ratio[minimum_idx:]\n",
    "        ll_left_side = ll_ratio[:minimum_idx]\n",
    "        \n",
    "        minimum = burst_size_indexed[minimum_idx]\n",
    "        burst_size_right_side = burst_size_indexed[minimum_idx:]\n",
    "        burst_size_left_side = burst_size_indexed[:minimum_idx]\n",
    "        try:\n",
    "            f_1 = interp1d(ll_left_side,burst_size_left_side, kind='linear')\n",
    "            f_2 = interp1d(ll_right_side,burst_size_right_side , kind='linear')\n",
    "            res = np.array([minimum, f_1(cutoff), f_2(cutoff)])\n",
    "            return res, burst_size_indexed, ll_ratio\n",
    "        except (ValueError,np.linalg.linalg.LinAlgError, TypeError):\n",
    "            return np.array([minimum,np.nan,np.nan]), burst_size_indexed, ll_ratio\n",
    "    vals_ = np.copy(vals)\n",
    "    start_estim = MaximumLikelihood(vals_)\n",
    "    try:\n",
    "        res_kon, kon_idx, kon_ll_ratio = kon_ll(vals_, start_estim)\n",
    "        res_burst_size, bs_idx, bs_ll_ratio = burst_size_ll(vals_,start_estim)\n",
    "    except BaseException:\n",
    "        return np.array([np.nan,np.nan,np.nan]), np.nan, np.array([np.nan,np.nan,np.nan]), np.nan\n",
    "    return res_kon, kon_idx, kon_ll_ratio, res_burst_size, bs_idx, bs_ll_ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy.random import beta, poisson\n",
    "def sim_beta_poisson(kon,koff,ksyn=100, size=100):\n",
    "    \"Generate simulated single-cell RNA-seq data according to the two-state model and with parameters kon, koff, ksyn\"\n",
    "    v = beta(kon, koff, size)\n",
    "    return poisson(v*ksyn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(12039213)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = sim_beta_poisson(0.5,0.5,20, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Y = ProfileLikelihood3(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFNCAYAAACUvLFdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VMX6wPHvmxBK6IQWQgmdkBAC\nRIrSkWIQFVGxIKAg4k8vKjZU9KLXdhUVuaDYu1e9CogUBZSqCNKL9JrQCT2hJXl/f+xmTSBlE7LZ\nJPt+nmee7J4z58y7i57ZmTNnRlQVY4wxxhf4eTsAY4wxJr9YpWeMMcZnWKVnjDHGZ1ilZ4wxxmdY\npWeMMcZnWKVnjDHGZ1ilZwo1cfhYRI6JyDIR6SAim9Ps3yUiV+fivJ1FJC7N+w0i0tn5eoyIfJGL\nc4aKiIpIMef7WSIyyPl6sIgszuk5c8MZQ4NM9mUaU1bHeZuITBKRZ7wdhyn4ink7AON7RGQXUA1I\nBhKAWcADqno6F6drD3QHaqpqgnNb47yIMy1VDffAOa/J63NeLk/HJCKfAHGqOvoyzjEYGKqq7VO3\nqerwy4/O+AJr6Rlv6aOqZYCWQDRwyUXQ2YrL7r/ROsCuNBWeKcRSW8HGeIpVesarVHUvjpZeBICI\nzBeRF0XkNyARqCciNURkmogcFZFtInKPM+8Q4AOgnYicFpHnLu6WTEtE/ERklIhsF5F4EflWRCq5\nE2dm3aQiEiAi/xWR70WkeE7KcH7WoRdtG+vsqt0pItek2Z7hd+DcV0JExonIPmcaJyIl0ux/TET2\nO/fdnc3nvCQmd/Kl7Qp1/lh5U0QOichJEVknIhEiMgy4A3jc+e/1Y5rv9gkRWQskiEixNN/hKRH5\nS0T6OvOGAZP4+9/8uHP7JyLyQpp47nF+T0ed31uNNPtURIaLyFYROS4iE0VEsvvMpmiwSs94lYjU\nAmKAVWk23wkMA8oCu4GvgTigBnAT8JKIdFXVD4HhwBJVLaOq/8ymuH8ANwCdnOc6Bky8jNhLAVOB\nc8Atqnr+MstoA2wGKgOvAh+muRhn+B049z0NtAWigOZAa5wtZxHpBTyKowu4IZDj+5u50APoCDQC\nygO3APGq+h7wJfCq89+rT5pjbgN6AxVUNQnYDnRwHv8c8IWIBKvqRtL/m1e4uHDn9/Kys9xg/v5v\nKK1rgSuASGe+nnnyyU2BZ5We8Zapzl/pi4EFwEtp9n2iqhucF7/qwFXAE6p6VlVX42jdDcxFmcOB\np1U1TlXPAWOAm3LZpVYO+AnHxfkuVU3OgzJ2q+r7znN9iuOCXc35wyCr7+AO4HlVPaSqh3FUEnc6\n990CfKyq651dwGNy8Vlz6gKOHyxNAFHVjaq6P5tjxqtqrKqeAVDV/6nqPlVNUdVvgK04KnN33AF8\npKornf8GT+JoGYamyfOKqh5X1T3APBw/GIwPsP5z4y03qOrcTPbFpnldAziqqqfSbNuN4z5gTtUB\npohISpptyTgG1ezN4bnaAgHAbZp+1vasysjOgdQXqprobOSVAYLI+juo4Xyfdl+NNPtWXLTPo1T1\nVxGZgKOFW0dEJgOPqurJLA5L+2+OiAwERgKhzk1lcLSA3VEDWJkmntMiEg+EALucmw+kyZ/oPL/x\nAdbSMwVR2kpkH1BJRMqm2VabnFdS4LiwXqOqFdKkks77ijk1G0cX2i8ikrZCy8syUmX3HezDUdmm\n3bfP+Xo/UOuifXkhAQhM87562p2qOl5VWwFNcXRzPpa6K5PzubaLSB3gfeABIMjZhbkekIvzZiLd\n9yEipXH8cLicfwNTRFilZwo0VY0FfgdeFpGSIhIJDAFy/JwcjgEQLzovqohIFRG5/jJiexX4CkfF\nl9oKydMynOVk9x38FxjtLKsy8Gyafd8Cg0WkqYgEAtnd93TXauBGEQkUx7N7Q1J3iMgVItJGRAJw\nVI5ngdSW70GgXjbnLo2jYjvsPN9dOAc6pTlHTREpnsnx/wXuEpEo54Cel4ClqrorJx/QFE1W6ZnC\n4DYc3Vz7gCnAP7PoGs3KW8A0YLaInAL+wDF4JNdU9V84BrPMdY7SzPMynLL6Dl4AlgNrgXU4uvZe\ncMY3CxgH/Apsc/7NC28C53FUQJ/iGKCSqhyOltoxHN2p8cBrzn0fAk2doyanZnRiVf0LeB1Y4jx/\nM+C3NFl+BTYAB0TkSAbHzwWeAb7H0dKtD9yaq09pihyxRWSNMcb4CmvpGWOM8RlW6RljjPEZVukZ\nY4zxGVbpGWOM8RlW6RljjPEZhW5GlsqVK2toaKi3wzDGGFOArFix4oiqVskuX6Gr9EJDQ1m+fLm3\nwzDGGFOAiIhbU+x5rHvTOXPEMhFZI45Vp5/LIM9gETksIqudKdslTYwxxpjcyrKl55zmJwbHEh81\ngDM45sCboaqbszn3OaCrc7LXAGCxiMxS1T8uyveNqj6Qu/CNMcYY92Xa0hORZ4ClQBdgDY6phqbh\nqCjHichPIhKR2fHqcNr5NsCZbPoXY4wxXpNVS2+tc17BjLwqIsGkn739EiLij2NZkwbARFVdmkG2\nfiLSEdgCPOycXPfi8wzDsagotWvn1STxxhhjfE2mLT1V/QFARMIy2b9fVZdldXJVTVbVKKAm0DqD\nluGPQKiqRgJzcLQmMzrPe6oararRVapkOzjHGGOMyZA7A1k+EpElIjLsovW83Kaqx3GsTtzrou3x\nzpWNwbESdKvcnN8YY4xxR7aVnqq2A+4GGgKrReQzEemS3XHOtb0qOF+XAroDmy7KE5zm7XXAxhzE\nbowxxuSIW8/pqepGEXkCx0KWE4B2InIBeDK1GzQDwcCnzvt6fsC3qjpdRJ4HlqvqNGCEiFwHJAFH\ngcGX93GMMcaYzGW7np6INAXuwtESmw98qKrLRKQWsFhV63g8yjSio6PVHk43xhiTloisUNXo7PK5\n09J7H8f9tjGqmpC6UVVjReSflxGj18zbfAiALo2rejkSY4wx+SnbSk9Vr8pi3yd5Gk0+GTdnC4pV\nesYY42syrfREZBUZP0wuOJ49b+mxqDwsplkwL8/aROzRRGpVCvR2OMYYY/JJVi29m/ItinyWWunN\nWLef4Z3qezscY4wx+SSrh9O3pyYc82jWdr6OAw7kV4CeUKtSIM1rlmfmuv3eDsUYY0w+yvY5PRG5\nG8ecmx84N9UBMntModCIaRbM2rgTxB5N9HYoxhhj8ok7M7KMANoCJwFUdQtQzZNB5YeYZo7n4q21\nZ4wxvsOdSu+sqp5PfeN82LzQq1UpkEjr4jTGGJ/iTqX3m4g8DpR0Tj/2DTDds2Hlj97NglljXZzG\nGOMz3Kn0HgdO4Zg380HgF+ApTwaVX6yL0xhjfIs7lV4zVX1HVfuq6g2q+g4XrZZQWFkXpzHG+BZ3\nlxZqmvpGRG4GnvdcSPkrxro4jTHGZ7hT6d0CfCEijUTkLuAhoIdnw8o/va2L0xhjfIY76+ltA24H\npjr/dlfVY54OLL9YF6cxxviOnMy9WcH5d7GIUJjn3rxYTLNgXrG5OI0xpsjzybk3L9bbWenNWr+f\nYR1tLk5jjCmqsurePOCca/NwJqnIqFUpkGYh5Zmx1ro4jTGmKMuq0vvO+XcDsN75d0Oa90VK70gb\nxWmMMUVdVqssXOP8W0tVazv/pqba+Rdi/kgdxTlrvbX2jDGmqMpqIEtkVgeq6tq8D8d7XF2c6w7Y\nfT1jjCmishrIMjGLfQp0zONYvC6mWTD//slGcRpjTFGVaaWnqh3yM5CCoLez0rNRnMYYUzS5MyOL\nz6gd9HcXpzHGmKLHKr2LxDQLZk3sceKO2ShOY4wpaqzSu4hrFKe19owxpshxq9ITkVtF5Gnn61oi\n0sqzYXlPahfndJuL0xhjipxsKz0RmQB0AQY4NyUAkzwZlLdZF6cxxhRN7rT0rlTVe4GzAKp6FCie\n3UEiUlJElonIGhHZICLPZZCnhIh8IyLbRGSpiITmMH6PsC5OY4zxvCmr4vjmzz2kpGj2mfOIO5Xe\nBRHxw7nigogEASluHHcO6KqqzYEooJeItL0ozxDgmKo2AN4E/u125B5UOyiQiJByzLAuTmOM8YjE\n80m8OGMjU1btRST/ynWn0psIfA9UcbbWFuNG5aQOp51vA5zp4ur8euBT5+vvgG4i+fnxMxfTLJjV\n1sVpjDEe8fFvuzhy+jyP9WxCfl723VlE9jNgNDAWOAbcrKpfu3NyEfEXkdXAIWCOqi69KEsIEOss\nJwk4AQS5H77nWBenMcZ4xonEC7y7YDtXh1WlVZ2K+Vp2lpWes9LaoKobVPUtVR2nqm6vsKCqyaoa\nBdQEWotIRG6CFJFhIrJcRJYfPpw/qxrVCSptXZzGGOMB7y7czqlzSTzSo3G+l51lpaeqycAOEQm5\nnEJU9TgwD+h10a69QC0AESkGlAfiMzj+PVWNVtXoKlWqXE4oOZLaxbn3+Jl8K9MYY4qyQ6fO8vFv\nu7iueQ3Cgsvle/nu3NMrA2wUkZ9FZHJqyu4gEakiIhWcr0sB3YFNF2WbBgxyvr4J+FVV828YTzb+\n7uK01p4xxuSFib9u43xyCg9f3cgr5We1ykKqF3J57mDgUxHxx1G5fquq00XkeWC5qk4DPgQ+F5Ft\nwFHg1lyW5RGpXZzT1+5naId63g7HGGMKtdijiXy1bA+3RNcitHJpr8SQbaWnqr+ISGUg2rlpuaoe\nceO4tUCLDLY/m+b1WeBm98PNfzHNgnn1p83sPX6GkAqlvB2OMcYUWuPmbkVEeLBbQ6/F4M6MLP2A\nlcCdwEBguYj09XRgBYV1cRpjzOXbevAUU1bFMahdHaqXL+m1ONy5p/cscIWq3qGqtwNtgDEejaoA\nqRNUmvAaNorTGGMux+uztxBYvBj3dW7g1TjcqfT8VPVgmveH3DyuyIhpFsyqPTaK0xhjcmNN7HF+\n2nCAoR3qUql0trNYepQ7lddsEZkhIgNEZADwIzDbw3EVKNbFaYwxuTd29mYqBgYwpH1db4fiVqX3\nKPAZ0NqZPnVu8xmhla2L0xhjcmPJ9ngWbT3C/V0aULZkgLfDcWv0pgLfiMiPafKXBU56MrCCJqZZ\nMK/9bKM4jTHGXarKaz9vonq5kgxoW8fb4QDujd4cKiL7gS3AemCD869PsS5OY4zJmV82HmLlnuOM\n6NaQkgH+3g4HcK978wmguarWVNXaqlpLVWt7OrCCxro4jTHGfSkpytjZmwkNCuTm6JreDsfFnUpv\nBz7WlZmZ1FGc+2wUpzHGZOnHtfvYdOAUD3dvRIB/wRnw704ko4DfRGSiiLyRmjwdWEGU2sU501p7\nxhiTqQvJKbwxZwtNqpelT2QNb4eTjjuV3iTgN2A1jvt5qcnnhFYuTdPgclbpGWNMFv63PI7d8Yk8\n1rMxfn4FYl1wF3cmnC6hqiM8Hkkh0TvSMYpz3/Ez1LBRnMYYk87ZC8mM/2UrLWtXoGuTqt4O5xLu\ntPRmiMjdzqWCyqUmj0dWQFkXpzHGZO7zJbs5cPIsj/dqgkjBauWBe5XeQOA5HJNOp3Zt+twjC6ms\ni9MYYzJ26uwF3p6/jQ4NK9O2XpC3w8lQtpWe8xGFi5PPPbKQVu/IYFbaKE5jjEnng0U7OZZ4gcd7\nNvF2KJkqOONIC5EY6+I0xph0jiac54NFO7gmojrNapb3djiZskovF+paF6cxxqTzzvxtnLmQzCM9\nGnk7lCxZpZdL1sVpjDEO+0+c4dMlu7mxZU0aVC3r7XCy5FalJyK3isjTzte1RKSVZ8Mq+FK7OGet\nP+DlSIwxxrvG/7INVeXBbg29HUq23JlwegLQBRjg3JSA44F1n1a3cmnCgssxY+0+b4dijDFes/NI\nAt8uj+X21rWpVSnQ2+Fky52W3pWqei9wFkBVjwLeXfq2gLjWujiNMT7uzTlbKO7vx/1dG3g7FLe4\nU+ldEBE/QAFEJAhI8WhUhYR1cRpjfNlf+04ybc0+7roqlKplS3o7HLe4U+lNBL4HqojIc8Bi4N8e\njaqQSO3itFGcxhhf9PrszZQrWYx7O9b3dihuc+fh9M+A0cBY4Bhws6p+7enACovezaqzYvcx9p+w\nLk5jjO9Ysfsov2w6xL2d6lM+MMDb4bgty0pPRPxFZIOqblDVt1R1nKr67BRkGfn7QXXr4jTG+AZV\n5dWfNlO5TAnuuirU2+HkSJaVnqomAztEJCSf4il06lUpY12cxhifsmjrEZbuPMo/ujYgsLg7i/UU\nHO7c0ysDbBSRn0VkcmrydGCFiXVxGmN8hary2s+bCalQiltb1/J2ODnmTqX3AtAXeBXHoJbUlCXn\nQ+zzROQvEdkgIg9mkKeziJwQkdXO9GxOP0BB4BrFaV2cxpgi7qf1B1i39wQPd29EiWL+3g4nx7Jt\nl6rqL7k8dxLwiKquFJGywAoRmaOqf12Ub5GqXpvLMgqEelXK0KR6WWas28/d7et6OxxjjPGI5BRl\n7OzNNKhahr4tCuddL3dmZLlCRP5wtsjOisg5ETmZ3XGqul9VVzpfnwI2AoXzW3LDtZHB1sVpjCnS\nJq+MY/vhBB7p3gh/v4K3QKw73OnefBsYBOwAygIPAONzUoiIhAItgKUZ7G4nImtEZJaIhOfkvAWJ\ndXEaY4qyc0nJjJu7lWYh5ekVUd3b4eSaO5Wen6puBoqp6gVVfR/o7W4BIlIGx8PtD6nqxS3ElUAd\nVW0O/AeYmsk5honIchFZfvjwYXeLzlepXZw2itMYUxR9vSyWvcfP8FjPxogUzlYeuFfpJYhIcWCN\niLwkIv8A3Lp7KSIBOCq8L1X1khGfqnpSVU87X88EAkSkcgb53lPVaFWNrlKlijtFe0XvZsEsty5O\nY0wRk3g+if/8uo02dSvRoeEll+hCxZ1Kb7Az3wNAMtAQuCm7g8TxU+BDYKOqvpFJnurOfIhIa2c5\n8W5FXgDFRFoXpzGm6Pn4t10cOX2Ox3sV7lYeuDd6c4fz5VngmRyc+yrgTmCdiKx2bnsKqO087yQc\nled9IpIEnAFuVVXNQRkFSv00XZw2itMYUxScSLzAuwu2061JVVrVqeTtcC6bxx6lV9XFQJY/CVR1\nAjDBUzF4Q+9mwbw+ZwsHTpylevnCMeu4McZk5t2F2zl5NolHezb2dih5wq2V0437Urs4bUCLMaaw\nO3TqLB//tovrmtcgLLict8PJEzmq9MShtKeCKQrq2yhOY0wRMfHXbZxPTuHh7o28HUqecefh9M9E\npJyIBALrgG0iMtLzoRVeqaM4D5w46+1QjDEmV2KPJvLVsj3cEl2LupWLTlvHnZZepPP5uhuAOUAd\nHCM6TSZcozjXW2vPGFM4jZu7FRFhRLcG3g4lT7lT6QWISDHgeuAHVT0PpHg2rMLNujiNMYXZ1oOn\nmLIqjkHt6hBcvpS3w8lT7lR6HwB7gIrAAhGpDZz2aFRFQEyzYP7cZV2cxpjC5405WygV4M99nYtW\nKw/cqPRU9U1VraGqPZzP0MUBXT0fWuHmmovTujiNMYXI2rjjzFp/gKEd6lGpdHFvh5Pn3BnIUkVE\n3hWR6c5NTYDbPRtW4degqnVxGmMKn9d+3kzFwACGdiiaE2y40735CbAASF0idyvwiKcCKkpibBSn\nMaYQWbI9nkVbj/B/nRtQtmSAt8PxCHcqvaqq+hXOwSuqegEbyOKWmGbBqFoXpzGm4FNVXvt5E9XL\nleTOdnW8HY7HuLvKQiVAwbGoLJDtIrLG0cXZuJp1cRpjCr5fNx1i5Z7jjOjWkJIBbi2kUyi5U+k9\nCvwI1BORBcB/gX94NKoipHeko4vz4Enr4jTGFEwpKcprP28mNCiQm6Nrejscj3Jn9OZyoAvQCXgQ\naKqqq7M+yqRydXFaa88YU0D9uHYfmw6c4uHujQjwL9pTMme7yoKIXDxSs6mI4LzPZ7KR2sU5Y91+\nBl9VNEdDGWMKrwvJKbwxZwtNqpelT2QNb4fjce5U6R3SpO7Ay7ixiKz5W+ooTuvizH9Tp07lnnvu\noX///syePdvb4RhT4PxveRy74xN5tEdj/PwK9wKx7nCne/O+NOkuIAooWvPSeFjvyOqFpovzwIED\n3HrrrdSvX59WrVoRExPDli1bMs1//Phx3n77bbfOPX78eMLCwrjjjjvSvXZHTspJ64YbbuD9999n\n0qRJfPPNNwCcOXOGTp06kZycDMBPP/1E48aNadCgAa+88kqm58osX2hoKM2aNSMqKoro6GjX9vPn\nz9OxY0eSkpLSnWf48OEsWLCALl260LRpU8LDw3nrrbdy/NnSuvvuu6latSoRERGZ5nnrrbeIiIgg\nPDyccePGARAbG5uncZjC5eyFZMb/spWWtSvQLayqt8PJH6qao4SjS3RzTo/Lq9SqVSstjHq8sUBv\nfud3b4eRpZSUFG3btq2+8847rm2rV6/WhQsXZnrMzp07NTw83K3zN27cWGNjYy957Y6clJORkSNH\n6ooVK1RVdcKECTpu3DhVVU1KStJ69erp9u3b9dy5cxoZGakbNmy45Pis8tWpU0cPHz6cYbljxozR\nL774It225s2ba1xcnCuekydPasOGDTMs110LFizQFStWZPodrVu3TsPDwzUhIUEvXLig3bp1061b\nt+q+ffvyNA5TuLy3YLvWeWK6/r7tiLdDuWzAcnWjDnFnRpYpIjLZmaYCG4Fpnq2Ki56YZsH8ufto\nge7inDdvHgEBAQwfPty1rXnz5tSqVStdC2Ls2LGMGTMGgFGjRrF9+3aioqJ47LHHAHjjjTeIiIgg\nIiLC1aIYPnw4O3bs4JprrqFEiRKu12+++eYlcdxwww20atWK8PBw3nvvvUzLSZWQkEBISAhPPfUU\nAH/++SdRUVGcPXsWVeWJJ57gmmuuoWXLlgB8+eWXXH/99QAsW7aMBg0aUK9ePYoXL86tt97KDz/8\ncElM7ubL6LN8+eWXrvcbN26kUaNGhISEuOIpW7YsYWFh7N27N9vzZaZjx45UqlQp0/0bN26kTZs2\nBAYGUqxYMTp16sTkyZMJDg7O0zhM4XE88Txvz99Gh4aVaVc/yNvh5JtsB7IAE9K8TgJ2q+ouz4RT\ndPWOrM6bc7cwqwAPaFm/fj2tWrXK0TGvvPIK69evZ/Vqx4DeFStW8PHHH7N06VJUlTZt2tCpUycm\nTZrETz/9xLx586hcuTKhoaGu1xf76KOPqFSpEmfOnOGKK66gX79+l5STVunSpVm7di3R0dGMHj2a\nu+++m88//5ySJUsyfvx45s6dy4kTJ9i2bRt33303O3bsIDQ0FIC9e/dSq1Yt17lq1qzJ0qVLLykj\nq3wiQo8ePRAR7r33XoYNG+bKFxERwZ9//ul6P2vWLHr16pXu3Lt27WLVqlW0adMm3fYOHTpw6tSp\nS2IZO3YsV1999SXbsxIREcHTTz9NfHw8pUqVYubMmem6YrOKwxRNr/68mZNnk3jymjBvh5Kvsq30\nVPWX/AikqGtQtazzQfUDBbbSywuLFy+mb9++lC7tWHTyxhtvZNGiRbRo0cLtc4wfP54pU6YAjntO\nW7dupXr16lkeExQURGBgIEOGDOHOO+8kKioKgBEjRjBixAhXvn379lGhQoWcfqwsLV68mJCQEA4d\nOkT37t1p0qQJHTt2BMDf35/ixYtz6tQpypYty88//8zHH3/sOvb06dP069ePcePGUa5cuXTnXbRo\nUZ7FGBYWxhNPPEGPHj0oXbo0UVFR+Pv//QByVnGYomd17HH+u2wPd11Zl6Y1fOvfO9PuTRE5JiJH\nM0jHRORofgZZVKR2cR4qoF2c4eHhrFix4pLtxYoVIyXl75nnzp71XPzz589n7ty5LFmyhDVr1tCi\nRQu3y4uMjOTAgQM8+uijmeYpVapUuvOFhIQQGxvreh8XF0dISMglx2WVL/Vv1apV6du3L8uWLUt3\n7Llz5yhZsiSJiYkcP36cGjUcw8IvXLhAv379uOOOO7jxxhsvKbNDhw5ERUVdkubOnevO13GJIUOG\nsGLFChYuXEjFihVp1KiRW3GYoiU5RRk9dR1VypTg4e4NvR1Ovsvqnl5loEoGKXW7ySHXKM71B7wd\nSoa6du3KuXPnXPfRANauXcvOnTs5dOgQ8fHxnDt3junTp7v2ly1bNl0XXIcOHZg6dSqJiYkkJCQw\nZcoUOnTo4HYMJ06coGLFigQGBrJp0yb++OOPDMu52OHDh5k3bx79+vXDz8/xn/WOHTsYMmQIN930\n9xM2FStWJDk52VXxXXHFFWzdupWdO3dy/vx5vv76a6677rpLzp9ZvoSEBFdcCQkJzJ49O939z/j4\neCpXrkxAQADz5s2jS5cugGMA2ZAhQwgLC2PkyJEZfqZFixaxevXqS1JOuzZTHTp0CIA9e/YwefJk\nbr/9drfiMEXLV0t3s37vSUZf27TITiqdlUwrPVVNTpuA8kC1NMnkUIOqZWlUrQwz1hbMRxdEhClT\npjB37lzq169PeHg4Tz75JNWrV+fZZ5+ldevWru67VEFBQVx11VVERETw2GOP0bJlSwYPHkzr1q1p\n06YNQ4cOzVHXZq9evUhKSiIsLIxRo0bRtm3bDMu52JAhQ+jatStr1qxxbatXrx4ffvjhJXl79OjB\n4sWLAUcrdsKECfTs2ZOwsDBuueUWwsPDXXljYmLYt29fpvkOHjxI+/btad68Oa1bt6Z3797p7tnN\nmzeP3r17A+nv5/322298/vnn/Prrr64W3MyZM93+ni5222230a5dOzZv3kzNmjVdnzs1foB+/frR\ntGlT+vTpw8SJE6lQoUKex2EKtsOnzvHqz5u5sn4QfSKDvR2Od2Q3vBPoDWwBEoFYHCssbHJnaKgn\nUmF9ZCHVm3M2a+io6XrwxBlvh1JkTJo0SW+55Rbds2dPhkP2+/Xrl+79ihUrdMCAAfkSW9++fXXz\n5s2qqtqiRQs9f/58vpRrTEYe/maVNnhqhm49eMrboeQ58uqRBeBF4Cocz+bVAnoCeXeH3cf0di03\nVDC7OAubrVu38uabbzJp0iRq1apFcHAwPXr0yPKYli1b0qVLF9fD6Z5y/vx5brjhBte9s5UrVxIQ\n4HvdSaZgWLbzKJNX7uWeDvVoULWMt8PxGnFUkFlkEFmuqtEisgaIUlUVkTWq2jx/QkwvOjpaly9f\n7o2i80yPNxdQIbA4397bztvujs5dAAAgAElEQVShFGnx8fE8/fTTzJkzh6FDh/Lkk096OyRjvOJC\ncgrXjl/M6XNJzBnZkcDi7jytVriIyApVjc4unzuf/ISIlAEWA5+JyCHgzOUG6Mv6RNbg9Tlb2Hbo\ntE//4vK0oKAgJk2a5O0wjPG6T37bxeaDp3jvzlZFssLLCXe6N2/AUck9BMwH9gJ9sjtIRGqJyDwR\n+UtENojIgxnkEREZLyLbRGStiLTMYfyF0u1talOimB/vL9zh7VCMMUXcgRNnGTd3C12bVKV7UxuD\n6E6ldxdQRVUvqOqHqvqGqh5247gk4BFVbQq0Be4XkaYX5bkGaOhMw4B3chB7oRVUpgS3RNdiyqq9\nBXpaMmNM4fevGX+RlKKM6ROOSNFfRSE77lR6VYD5zlbbcBG5dN6oDKjqflVd6Xx9CsecnRc/9Xs9\n8Jlz8M0fQAUR8YlxtPd0qEdSSgof/bbT26EYY4qohVsOM2Ptfu7v0oDaQYHeDqdAcGdpoWdUtQnw\nCFAXWCIiP+WkEBEJBVoAF09qGILjMYhUcVxaMRZJtYMCiWkWzFd/7OHk2QveDscYU8ScS0rmn9M2\nEBoUyLCO9bwdToGRk3XhY4FdwD6gtrsHOQfBfA88pKoncxTd3+cYJiLLRWT54cPu9KwWDsM71efU\nuSS+WrrH26EYY4qY9xbsYOeRBJ67PoKSAf7ZH+Aj3FlaaJiIzMXxbF4I8A/nfbpsiUgAjgrvS1Wd\nnEGWvUCtNO9rOrelo6rvqWq0qkZXqVJ0ZkCLCClP+waV+WjxTs4lefaZMWOM74g9msiEeduIaVad\nTo2KzjUzL7jT0msIjFLVJqo6WlXXunNicdwx/RDYqKpvZJJtGjDQOYqzLXBCVQvmHF0ecm+nehw6\ndY6pq2wNM2NM3njuxw34+wnPXOtW+8SnuLO00KUTHbrnKuBOYJ2IpC6C9hTOrlFVnQTMBGKAbTim\nObsrl2UVWu0bVCa8RjneXbiDm1vVws/PRlcZY3Jvzl8HmbvxEE/FNCG4fClvh1PgeOwpRVVdDGR5\nBXfOl3a/p2IoDESEezvVZ8R/VzFn40F6hme9bpwxxmTmzPlkxkzbQKNqZbirCK/beTl8+9H8vPDQ\nQ5DBat450UehVtxxAiYLWqN81r8UjDEmE0eOJvL68TM0rVGOgB9yOM9rVBSMG+eZwAqQLO/piYi/\niMzJr2B8laxZTfihHZw+m8Sps0neDscYUwiduZDMvhNnqVymBOU2bbjsH+NFVZYtPVVNdlZ85XL7\nuEGRlxe/jDp3xl+V/7t6NC1qVeDDwVdc/jmNMT5DVbnnw2WsiTvOr490hj49HTvmz/dmWAWSWxNO\nA2tEZDaQkLpRVW2Z5TzkL8KgdqG8OXcLWw6eolG1st4OyRhTSMxYt5/F247w3HXhVClbwtvhFGju\nPLIwHXgBWAZsSJNMHhvYrg6lAvx5d4FNRG2Mcc/pc0n8a/pfhNcox4C2dbwdToHnziMLH4pIcaC2\nqm7Lh5h8VsXSxel/RS2++GM3j/ZsZMONjTHZGjdnC4dOnWPSgFb42yNP2XJnRpbewDpgjvN9lIhM\n8XRgPuWWWxwJGNK+Lgp8tNgmojbGZG3TgZN8/Psubr2iNi1qV/x7R5priknPnXt6zwNtgHkAqrpa\nRBp4NCpf83//53pZq1IgfSKD+WrpHh7o0pDygTkcdmyM8QmqyjNT11OuZDEe79k4/c401xSTnjv3\n9C6o6vGLtqkngvFZiYmO5DSsY30SzifzxdLdXgzKGFOQfb9yL3/uOsaoa5pQsXTx9DsvuqaYv7lT\n6W0UkVsAPxGpKyJvAn94OC7fEhPjSE5Na5SjY6MqfPzbLs5esImojTHpnUi8wMszN9KydgVublXr\n0gwXXVPM39yp9B4AWgEpwBTgHPCQJ4MyMLxTPY6cPsfklTYRtTEmvddmb+JY4nn+dUOEzdebQ+5U\nepVV9QlVbaGqUao6SlWt3exh7eoFEVmzPO8t3E5yivUmG2Mc1sYd58ulexjYLpTwGuW9HU6h406l\n95WIbBWRL0TkXhEJ83hUBhFheKf67IpPZPaGA94OxxhTACSnKKOnrqdymRKM7NHI2+EUStlWeqp6\nFdAMeB+oCvwsIkVn+fICrGd4dUKDApm0YDuOBSmMMb7sv8v2sDbuBKN7h1GupI3szo1sH1lwLu7a\nAegIVAZ+wrGKuskrgwdnuNnfT7inYz2enrKeP3YcpV39oPyNyxhTYBw5fY7Xft5Mu3pBXNe8RtaZ\nM7mmGPee0/sNWA68DExXVVsGIK9l8R9ov5Y1eXPOFt5duN0qPWN82CuzNpFwLol/3RCOSDaDV6zS\ny5Q79/SqAS/haO3NFZHZIvJPz4blY44ccaQMlAzwZ/CVoczffJiN+22hC2N80fJdR/luRRxDO9Sj\nQVU3JqPP4pri69y5p3cE+AvYCOwGGgE9PByXb7npJkfKxJ1tQwks7s97C20iamN8TVJyCqOnrqdG\n+ZKM6ObmZFjZXFN8mTtzb24HJgI1gI+BMOfgFpNPygcGcFvr2kxbs4+4Y/a0iDG+5JPfd7HpwCme\n7RNOYHF37kiZrLjTvdlIVXuo6vOqOl9Vz3g8KnOJIe3rIsCHNhG1MT7jwImzvDlnC50bV6FneDVv\nh1MkuHVPT0T+JyL7nekbEclm6JDJazUqlOK6qBp8vSyWYwnnvR2OMSYfvDDjLy6kKM9d58bgFeMW\ndyq9j4HZQKgzzXFuM/ns3o71OXMhmc//sImojSnqftt2hOlr9/N/netTJ6i0t8MpMtzpIK6mqu+n\nef+BiDzgqYB80n33uZWtcfWydG1SlU9+38WwjvUoGeDv4cCMMd5wLimZZ35YT52gQIZ3qp/zE7h5\nTfFF7rT0jorIrfK3/sBRTwfmU/r3dyQ33NuxHkcTzvO/5bEeDsoY4y0fLNrJjsMJjLkuPHc/bnNw\nTfE17lR6dwMDgSPAYeBO4C5PBuVzYmMdyQ2t61aiRe0KvL9oJ0nJKR4OzBiT32KPJvKfX7fSK7w6\nXRpXzeVJ3L+m+Bp3ntPbpaoxqhqkqpVV9VqgTz7E5jvuvNOR3CAi3NuxPnuOJvKTTURtTJHz/PS/\n8BPh2T5Nc3+SHFxTfI07Lb2MPJ6nUZgc6dG0GvUql7aJqI0pYn7ZeJA5fx1kRLeG1KhQytvhFEm5\nrfRs7KwX+fkJwzrWY/3ek/y+Pd7b4Rhj8sDZC8mM+XEDDaqW4e6r6no7nCIrt5Vets0LEflIRA6J\nyPpM9ncWkRMistqZns1lLD7phhYhVClbgkkLtns7FGNMHnh73jZij57hX9dHULxYbi/NJjuZfrMi\nckxEjmaQjuGYkiw7nwC9ssmzyLkae5SqPp+DuH1eyQB/7r6qLou2HmH93hPeDscYcxl2Hklg0oId\n3BBVw1ZT8bCsntOrfDknVtWFIhJ6OefwGY88kqvDbm9Tm4nztvHuwh3857YWeRyUMSY/qCrP/rCe\nEsX8eKp3WN6cNJfXFF+QVRu6hKomZ5YARCTwMstvJyJrRGSWiIRf5rkKrz59HCmHypcK4I42tZmx\ndh+xR20iamMKo5nrDrBo6xEe6dGIqmVL5s1Jc3lN8QVZVXrTReTfInKliLj+JUSktogMEpGZwLWX\nUfZKoI6qNgf+A0zNLKOIDBOR5SKy/PDhw5dRZAG1ebMj5cJdV9XF30/4YJEtO2RMYXP6XBL/mv4X\nTYPLMaBtnbw78WVcU4q6rCq9bjhWTX8Q2Coix0XkKPAdjjk471HVb3NbsKqeVNXTztczgQARybBL\nVVXfU9VoVY2uUqVKbossuO6915FyoXr5kvRtEcI3y2OJP30ujwMzxnjS+F+2cuDkWV7oG0Ex/zwc\nvHIZ15SiLtNvWR2mqWp/Va2lqhVUtZKqtlbV51R17+UULCLVxTltuIi0dsZi4+9zYVjHepy9kMJn\nS2wiamMKiw37TvDh4p3cekUtWtau6O1wfIbHxsWKyH+BJUBjEYkTkSEiMlxEhjuz3ASsF5E1wHjg\nVrUnrXOlQdWyXB1WjU+X7CLxfJK3wzHGZCPhXBL/+GoVlcsU54leTbwdjk/x2DK8qnpbNvsnABM8\nVb6vua9zPfq9c5Bv/4xlsD3YakyB9uwPG9gVn8BX97SlYuni3g7Hp9gTkEVEqzqViK5TkfcX7eSC\nTURtTIE1eWUc36+M4x9dG9K2nj2Tl9+yeji9XFYpP4Ms8kaPdqTLNLxTffYeP8PMdfvzIKiiLz4+\nnqioKKKioqhevTohISGu9+fP593q9HPnzqV8+fJERUXRrFkzevToQV6NQv7oo484cCDjiccHDBhA\n3bp1XZ9p4sSJeVJmQRQfH0+3bt1o2LAhPXv25MSJSydsSEpKwt/f3/V99O3bN9/j3H74NKOnrqdN\n3UqM6NbQcwXl0TWlKJLMbqOJSCyO6cYExwwsp5yvywD7VLVWfgWZVnR0tC5fvtwbRRd4KSlKz3EL\nKebvx8wR7XGOEzJuGDNmDGXKlOHRRx9Nt11VUVX8/HLfKTJ37lwmTJjA1KmOp3Iee+wxypUrxzPP\nPOPW8VnF0L59eyZMmEBUVNQl+wYMGMBNN93EDTfckOm5k5OT8fcv/IsRjxw5kho1avDoo4/ywgsv\ncObMGV588cV0eZKSkqhcuTLHjx/3SoxnLyTT9+3fOXDiDLMe7Ej18nn0TJ4BQERWqGp0dvmyGr1Z\nS1VrAzOAvs7Rm+WBG4DpeReqYfVqR7pMqRNRb9x/koVbj+RBYL5p27ZtNG3alDvuuIPw8HBiY2Op\nUKGCa//XX3/N0KFDATh48CA33ngj0dHRtG7dmj/++CPLc6sqp0+fpmJFx2i90aNHM27cONf+Jk2a\nEBcXl2EMd955J82aNSMiIoLx48fzzTffsHr1avr37+926zQpKYkKFSrw0EMPERkZybJly/jzzz/p\n1KkTrVq14pprruHgwYMA/Pnnn0RGRhIVFcWjjz7qqlg/+OADHnroIdc5e/XqxeLFiwGYNWsW7dq1\no2XLlvTv35+EhAQAatasyZgxY2jRogWRkZFs2bIFgFOnTjFo0CAiIyOJjIxk6tSpvPfee+l+fLzz\nzjs89thjWX6uH374gUGDBgEwaNAg1w+MguTlmRvZuP8kr9/S3PMVXh5dU4oid36+XqWq01LfqOqP\nwFWeC8kHPfSQI+WB66NCqFauBO/aRNSXZdOmTTz88MP89ddfhISEZJpvxIgRPP744yxfvpxvv/3W\nVRlebN68eURFRVGrVi0WLFjA4MGDcxTDgQMHOHLkCOvWrWP9+vUMHDjQVdmlVn7Fi186IOLhhx92\ndef99ddfAJw4cYKOHTuydu1aWrZsyYMPPsj333/PihUrGDBggKsFOnjwYN555x1Wr15NcnJytvEe\nOnSIV155hV9++YWVK1cSGRnJW2+95dpfrVo1Vq1axdChQ3njjTcARwu7SpUqrF27ljVr1tCpUydu\nvfVWpkyZQlKSYyTyxx9/zN133w1Az549OXTo0CVlx8fHk/oMb0hICPv3Z9zFn5CQQKtWrWjXrh0/\n/vhjtp8pr/y0/gCfLtnN0PZ16dqkmucLzMNrSlHjzujN/SIyCvjC+f4O4KDnQjKXo3gxP4a0r8tL\nMzexNu44kTUrZH+QuUT9+vWJjs62p4S5c+eyOc3MF8eOHePMmTOUKpV+LbQuXbq4Wh8vvvgio0aN\nYsKErAcvp42hQYMGbN68mREjRtC7d2969Ojh1ud4880303VvJiUlUbx4cdf9rI0bN7Jhwwauvvpq\nwNHdWbNmTY4cOcKZM2e46irH79s777yTefPmZVnW77//zl9//cWVV14JwPnz52nfvr1r/4033ghA\nq1atmDlzJuD4/lK/FxFxtYA7duzIrFmzqFevHv7+/oSFOeak/Pnnn9363Bnx9/dn9+7d1KhRg23b\nttGtWzeaNWtGaGhors/pjrhjiTz+3Roia5bncXs8wevcqfRuB54DZjnfLwSyfBzBeNdtrWvzn1+3\n8e6CHUy8o6W3wymUSpcu7Xrt5+eXbrHes2fPul6rKsuWLcuwlZWZ6667jjvuuAOAYsWKkZLy92jb\ntOdOG0NQUBBr165l1qxZTJw4ke+//5733nsvZx/KqVSpUq77vapKZGQkixYtSpfnyJHMu8czi1lV\n6dWrF59//nmGx5UoUQJwVD6prbjMpLYGQ0NDueuuu7L9TEFBQRw+fJgqVaqwd+9egoODL8kjItSo\n4VggpkGDBnTo0IHVq1d7tNK7kJzCiP+uIkXhP7e1sCWDCoBs/wVU9Yiq3g+0AVqr6v2qajeMCrCy\nJQMY0LYOs9bvZ9eRBG+HU+j5+flRsWJFtm7dSkpKClOmTHHtu/rqq9ONilztxn2UxYsXU79+fQBC\nQ0NZsWIFAMuWLSM2NjbDYw4fPoyqcvPNN/P888+zcuVKAMqWLcupU6dy/dmaNm3K3r17WbZsGeBo\nnW3YsIHKlStTqlQplixZAsCXX37pOiY0NJRVq1ahquzatcsV/5VXXsmCBQvYscMxD2xCQgJbt27N\nsvzu3bu7vj9V5dixYwBcddVVbN++nf/973/0798/289x3XXX8emnnwLw6aefcv3111+S5+jRo5w7\n55iq7/DhwyxZssTVgvSUN+dsYeWe47x8YzPqBJXO/gDjcdlWeiISLiJ/AltwzMG5VESaej40cznu\nuiqUYn5+vG8TUeeJf//73/Ts2ZMrr7ySmjVrurZPnDiR3377jcjISJo2bcr777+f4fGp9/SaN2/O\n119/zWuvvQbAzTffzMGDB4mIiOC9996jXr16GR4fGxtLx44diYqK4q677uKll14C4K677mLo0KG5\nfsyiRIkSfPfdd4wcOZLIyEhatGjB0qVLAce9tHvvvZeoqKh0I0c7depESEgIYWFhPPLII64BLtWq\nVePDDz+kf//+NG/enCuvvNI1YCUz//znP12fPyoqKl2L86abbqJjx46UL1/etS2ze3pPPfUUM2bM\noGHDhixcuNA18GXp0qUMH+6YBGrDhg1ER0fTvHlzunXrxjPPPEPjxo1z/J25a9HWw7yzYDu3ta5F\nn+buLEFq8kOmjyy4MogsBp5T1TnO91cDY1S1fZYHekiRfGTh998df533QvLKk5PX8v3Kvfz2RFeq\nlC2Rp+c2vmXbtm3cdNNNbrVk80qvXr148skn6dSpU76VmVcOnTpLzFuLqFS6OD/c355SxfP5sRAP\nXVMKsst+ZCGNsqkVHoCqzgXKXk5w5iJXXumR/zjv6VCPC8kpfPr7rjw/tzGeEh8fT8OGDalYsWKh\nrPBSUpRHvl3D6XNJTLi9Zf5XeOCxa0pR4M5All0i8iSQend6ALDLYxH5Ig/9KqtXpQw9m1bnsyW7\nuK9zfUqX8NhUq6aIa9CgQb618oKCgrK9F1iQvbNgO4u2HuGVG5vRqJqX2gc+2NJzlzstvbuBWsBM\nZ6rl3GbyylNPOZIH3NupHifPJvHfZXs8cn5jzN9W7D7KG3O2cG1kMP2v8MqkVQ4evKYUdtn+9FfV\neOD/RCTQ8VbPeD4sk1da1K5Im7qV+HDxTgZdGUpAXi5UaYxxOZ54nhH/XU1IhVK8fGMzmwawgLLR\nmz5geKf67D9xlh/X7PN2KMYUSarK49+t5dCps0y4vQVlSwZ4OySTCXd+9r8LPKWqNVW1JvA0kLun\nYo1XdG5chcbVyvLugh1kN1rXGJNzny3Zzey/DvJEryY2C1IBZ6M3fYCIcG+nemw+eIr5m/NmSRtj\njMOGfSd4ccZGujapypD2toBzQWejNwuCNLPse0qf5jUY+/Nm3lmwnS5Nqnq8PGN8QcK5JP7x1Soq\nlg5g7M3NC859vHy4phRWNnqzIIiKciQPCvD3Y0iHeizbeZSVe455tCxjfMUzU9ezKz6Bt25tQaXS\n7s+/6nH5cE0prNyZezNeVf9PVSOd6X7niE6TV+bOdSQPu/WKWpQvFWDLDhmTB75fEcfkVXsZ0a0h\nbesFeTuc9PLpmlIYZdu9KSINgJFAaNr8qure2iYmey+84PjrXN7FU0qXKMbAdnWYMG8b2w+fpn6V\nMh4tz5iiavvh0zzzw3ra1K3EP7o29HY4l8qna0ph5E735nfARuAF4Jk0yRRCg64Mpbi/Hx/YRNTG\n5MrZC8nc/+VKSgb489atLfD3KyD38Yxb3BnIkqKq//F4JCZfVC5Tgpuja/Ltn3E82K0R1cuX9HZI\nxhQqL83cyKYDp/h48BX2/08hlGlLT0TKiUg54AcRGSYiVVK3ObebQurejvXx84NRk9fac3vG5MBP\n6/fz2ZLd3NOhro2CLqSy6t7cAKwHhuLozlzp3Ja63RRStSoF8lRMGPM3H+aLpb49J+dbb71FREQE\n4eHhjEszzHvMmDGEhIQQFRVFVFQUM2fOBHCtnRcdHe2aFPn48eP06NEj3WriaV24cIFRo0bRsGFD\nWrZsSbt27Zg1axbgWJA1q1XKATp37kzjxo2ZNm0a4FgMtXv37jRs2JDu3bu7Fl7NyMmTJ6lZsyYP\nPPAAAImJifTu3ZsmTZoQHh7OqFGjXHk/+eQTqlSp4vrMH3zwQXZfHytWrKBZs2Y0aNCAESNGZPoj\nav78+URFRREeHn7JygnJycm0aNGCa6+91rXt119/pWXLlkRERDBo0CDXSuvffPMNDRo0SJc3v8Qe\nTeTx79bSvGZ5HuvZJN/LN3lEVQtVatWqlRY5mzY5Uj5KSUnROz9cqo1Hz9Tth07la9kFxbp16zQ8\nPFwTEhL0woUL2q1bN926dauqqv7zn//U11577ZJj+vbtq7Gxsbpo0SIdOXKkqqo+8sgjOm/evEzL\neeKJJ3TgwIF69uxZVVU9cOCAfvPNN6qqWqdOHT18+HCWcXbq1En//PNP1/vHHntMX375ZVVVffnl\nl/Xxxx/P9NgRI0bobbfdpvfff7+qqiYkJOivv/6qqqrnzp3T9u3b68yZM1VV9eOPP3blc9cVV1yh\nS5Ys0ZSUFO3Vq5frXGkdO3ZMw8LCdPfu3aqqevDgwXT7X3/9db3tttu0d+/eqqqanJysNWvW1M2b\nN6uq6jPPPKMffPCBK/+8efNcefPL+aRkvWHiYo149ifdfSQhX8vOFS9cU7wNWK5u1CFZdW92cv69\nLqOUb7WyL2jc2JHykYjw2k2RlAzw5+FvVnMhOeNWSlG2ceNG2rRpQ2BgIMWKFaNTp05Mnjw5y2MC\nAgJITEwkMTGRgIAAtm/fTmxsLJ07d84wf2JiIu+//z7/+c9/KFHCsZBvtWrVuOWWW3Id9w8//MCg\nQYMAGDRoEFOnTs0w34oVKzh48CA9evw90DowMJAuXboAULx4cVq2bElcXFyu4ti/fz8nT56kbdu2\niAgDBw7MMJavvvqKG2+8kdq1awNQterf3YJxcXHMmDGDoUOHurbFx8dTvHhxGjVqBED37t35/vvv\ncxVjXnljzhZW7TnOSzc2o3ZQoFdjcYsXrimFRVbdm92df2/OIN2U3YlF5CMROSQiGXaFisN4Edkm\nImtFpGUOYy86fvzRkfJZtXIlealvM9bEnWDCr9vyvXxvi4iIYNGiRcTHx5OYmMjMmTOJjY117Z8w\nYQKRkZHcfffdri7EJ598koEDB/Lyyy/zwAMP8PTTT/NC6vDwDGzbto3atWtTrlz2t8FjYmLYty/7\nScEPHjxIcHAwANWrV+fgwYOX5ElJSeGRRx5h7NixmZ7n+PHj/Pjjj3Tr1s217fvvvycyMpKbbrop\n3XeRkb1791KzZk3X+5o1a7J3795L8m3ZsoVjx47RuXNnWrVqxWeffeba99BDD/Hqq6/i5/f3pahy\n5cokJSWxfPlyAL777rtsY/GkhVsO88787dzWuhZ9mtfwWhw54qVrSmGQaaWnqqOdf+/MIA1049yf\nAL2y2H8N0NCZhgHvuB92EfP6647kBTHNgrmxRQgT5m1jlY/N1BIWFsYTTzxBjx496NWrF1FRUfj7\nO1a5vu+++9i+fTurV68mODiYRx55BICoqCj++OMP5s2bx44dOwgODkZV6d+/PwMGDMiwAnLXzJkz\nqVEjZxdVEclw6qu3336bmJiYdJVSWklJSdx2222MGDGCevXqAdCnTx927drF2rVr6d69u6s1ebmS\nkpJYsWIFM2bM4Oeff+Zf//oXW7ZsYfr06VStWpVWrVpd8pm+/vprHn74YVq3bk3ZsmVd/y757dCp\ns4z8djWNqpXh2WvDvRJDrnjxmlLQZfrIgoiMyOpAVR2fzf6FIhKaRZbrgc+cfbF/iEgFEQlW1f1Z\nndfkvTHXh7N051FGfruGGSPaE1jcd1ZYHzJkCEOGDAHgqaeeclUS1apVc+W55557Lhk4oaq88MIL\nfP311/zjH//g1VdfZdeuXYwfP54XX3zRla9Bgwbs2bOHkydPutXac0e1atXYv38/wcHB7N+/P113\nYaolS5awaNEi3n77bU6fPs358+cpU6YMr7zyCgDDhg2jYcOGPPTQQ65jgoL+nlVk6NChPP7441nG\nERISkq5rNC4ujpCQkEvy1axZk6CgIEqXLk3p0qXp2LEja9asYeXKlUybNo2ZM2dy9uxZTp48yYAB\nA/jiiy9o164dixYtAmD27Nls2bIlZ19SHkhJUR7+ZjWnzyXx1T1tKVXcOxWvyVtZdW9WySZdrhAg\nbZ9FnHObyWflSgbw+i3N2RWfwIszNno7nHx16NAhAPbs2cPkyZO5/fbbAcf9qlRTpkwhIiIi3XGf\nffYZMTExVKpUicTERPz8/PDz8yMxMTFdvsDAQIYMGcKDDz7I+fPnATh8+DD/+9//ch3zddddx6ef\nfgrAp59+yvXXX39Jni+//JI9e/awa9cuxo4dy8CBA10V3ujRozlx4kS60aoXf+Zp06YRFhbmet+k\nyaWjFYODgylXrhx//PEHqspnn32WYSzXX389ixcvJikpicTERJYuXUpYWBgvv/wycXFx7Nq1i6+/\n/pquXbvyxRdfAH//uwUa/mEAACAASURBVJw7d45///vfDB8+PKdf02V7Z8F2ftsWz5g+4TSqZgvL\nFBWZ/qRX1QIz64qIDMPRBeq6GW7yVtt6QQzrUI93F+6gW1hVujaplv1BRUC/fv2Ij48nICCAiRMn\nUqGCYy20xx9/nNWrVyMihIaG8u6777qOSUxM5JNPPmH27NkAjBw5kpiYGIoXL85XX311SRkvvPAC\no0ePpmnTppQsWZLSpUvz/PPPX5IvJiaGDz74INsuzlGjRnHLLbfw4YcfUqdOHb799lsAli9fzqRJ\nk7J81CAuLo4XX3yRJk2a0LKl4zb6Aw88wNChQxk/fjzTpk2jWLFiVKpUiU8++QSAI0eOZPoowttv\nv83gwYM5c+YM11xzDddccw0AkyZNAmD48OGEhYXRq1cvIiMj8fPzY+jQoZf8iLjYa6+9xvTp00lJ\nSeG+++6ja9euWebPa8t3HeWNOVvo07wG/a+ola9lG8+SzP5jdmVwzL05Eaiuqs1FJBLoraovZ3ty\nR/fmdFW95L9wEXkXmK+q/3W+3wx0zq57Mzo6WlNvcBcZqSP/5s/3ZhScS0rm+gm/ceT0eX5+qANB\nZUp4NR7j0LlzZ8aOHUt0dLRXyp8+fTo7duxgxIgs73jkm/nz5zN27FimT5/ukfMfTzxPzFuLKObv\nx4wR7QvnKugF5JqSn0Rkhapm+z+JO3NvfgA8B6SOaV+HY029yzUNGOgcxdkWOOGz9/M+/9yRvKxE\nMX/G3RrFyTMXeHLyOputpYD4//buO7yKKn3g+PdNgISSkARCJ/TeIRTrzxVdsWOjKCqswrqKa1vb\nYl9x19W1rawoKLpKNRZAWZHVxQYIoXelhapICDWUlPf3x53gFVMmJDdzy/t5nnnuvXNn5r7HwXkz\n58w5JykpiaFDh57onF7RLrnkkqBJeFOnTuXWW28lMTExIMdXVe5LW8FPh47x8rXdQjPhQdBcU4KR\nmycWqqvqvIInxFRVRSSnpJ1EZDJwDlBbRLYDjwKVnWOMxTc330XABiAbGHYqBQgLjYOn+qRtvXju\nvaANo2et5d3F2xmQGjyxRaqS+g5GkoEDBzJw4MCAHf+teVv4dM2PPHRxOzo3SgjY7wRcEF1Tgo2b\npJcpIs0ABRCR/sAPJe2kqoNL+F6B29wEGfamTvW9BvB/5tK46cxmfL5uN4/PWE2fZrVCozOuMWW0\nasd+npq1jr5t63DTmc28DqdsguyaEkzcVG+OBF4H2opIBvAA8IeARhVpXnnFtwSJqCjh2QFdiBLh\n7mnLyMu3ak4T3g4dy+X2yUtJql6FZ67pUmjfx5ASZNeUYOIm6W1V1XOB+kAXVe0DRFYv5gjUMKEq\nT/TvQHpGFmNtpnUTxlSVhz9cRUbmYV4c1JWk6lW8DskEkJuk956IRKvqflXdJyJ1gDmBDsx4r3/X\nhlzcuT7Pz/mOVTv2ex2OMQGRtng7HyzdwR19W9O7ea2SdzAhzU3SmwVME5EoEUnBl/AeCWxYJhiI\nCKP7d6RWjSrcNXUZR3PyvA7JmHK1YfchHpm+mj7Nkxh5bkuvwzEVoMSkp6qvAF8C7wMfAyNV9T+B\nDswEh4RqVXj2mi58v/sQT3+yzutwjCk3R3PyGDlpCVWrRPPioG5ER4V4O55xpTRjbzYFlgLdRKRb\nSWNvmlJIS/M6gmKd1SqZoac3ZcI3W+jbti5ntqrtdUjGlElevq8/3rofDjJhaE/qxsd6HVL5CvJr\nipfcjr1ZG5gJZFB+Y2+aArVr+5Yg9sCFbWmRXJ0/vbucfdnHvQ7HmFOWn6/c/94KZizfyf392vKb\ntr8esDvkhcA1xSshMfZm2HPGOGToUC+jKFZsZV8VUP8x3/Dw9NX8c3A3r0MyptRUlYemryJt8Xbu\nPK8VfzinhdchBUYIXFO8Ulz15j9U9R4R+QCnY7o/Vb0yoJFFkhD5B9qxYU3uOr81z8xez3nt6nB5\nV5sUw4QOVeXxmWuY9O1Wbj2nBXf0beV1SIETItcULxQ3IovTpZ+XKyIQExp+f3ZzPl+3m4c+XEXP\npkk0SKjqdUjGlEhVeWrWWt6ct4XhZzXj3gvahH4HdHNKips5faHz+llhS8WFaIJJpegonh/Qlfx8\n5Z5py8m30VpMkFNVnv10PeO+2szQ05vy54vaWcKLYEUmPRFZKiJLiloqMkgTXFJqVeORS9szf1Mm\nb3yz2etwjCnWS59tYMz/NjK4VwqPXtreEl6EK6568+oKi8KEnAGpjZmzZjd/n72es1ol06aezSxt\ngs8rczfy/H+/4+oejRjdv6MlPFPyJLLBJiwnkc3O9r1WC63ZDPYcOka/F74kOS6WD287nZhK0V6H\nZMwJ47/axJMfr+Xyrg14bkDXyOp8HqLXlLIoz0lkTaBVqxaS/zhr14jh6as6s3bXAZ6f873X4Rhz\nwtvzt/Dkx2u5qFM9/nFNl8hKeBCy15SKYEkvGPzrX74lBPVtV5fBvVJ49cuNfLsp0+twjGHKwq08\nPH0157Wry4uDulEpOgIvcyF8TQm0CPzXEISmTfMtIeqhi9vRJKkad09bzsGjOV6HYyLYe4u38+AH\nKzmnTTJjrutG5UhMeBDy15RAsqc3TZlVj6nEcwO7smv/ER6bscbrcEyEmrF8J/emLeeMFrUZO6SH\ntTGbQrl5evMWIBp42/l8HWBzzJhf6J6SyMjftOSlzzdwXrs6XNipvtchmQjyn5W7uGvqMno2TWLc\nDanEVraEZwpX3NibGwFEpK+qdvf7aqlzp3d/oIMzoeX2vq2Y+91P/PmDlfRokkidcBu53gSl/675\nkdsnL6Vr4wTeGNqTqlUs4ZmiuanwjhaRPgUfRKQ3vjs/Y36hcnQUzw/sypGcPO5NW0GodYcxoWfu\n+t3cOnEJHRrEM2FYT6rHFFd5ZUzx1ZsFbgYmiEjBn+1HgN8FLqQINHeu1xGUmxbJNRh1UTsenr6a\ndxZkcP1pTb0OyYSpbzbs4fdvL6ZV3Rr8+3e9iY+t7HVIwSOMrinlrcSkp6qLgI4iUsv5bM+lm2IN\n6dOEOWt3M3rWWk5vWZsWyTW8DsmEmYWb93LzW+k0q12dt2/qTc1qlvCMOyVWb4pIsoi8Crylqpki\n0l5EhgY+tAjy7LO+JUyICM9c3ZnYytHcNXUZOXn5XodkwsjijCyGTVhIg4RY3rm5N0nVq3gdUvAJ\ns2tKeXLTpvcm8AXQ2Pn8PXBPoAKKSB995FvCSN34WP56RSdWbN/PPz/f4HU4Jkws37aPoW8sJDku\nhknD+1C7RozXIQWnMLymlBc3Sa+Oqk4C8gFUNafgvTHFubBTfa7s3pAx/9vAkq1ZXodjQtzqnfu5\n/vVvSahemUnD+1DXng42p8BN0jssIkk4s6eLSE/gQECjMmHjscs6UC8+lrunLuPwsVyvwzEhav0P\nBxky/lviYisz6eY+NnmxOWVukt6fgJlAcxH5ApgM3O7m4CLST0TWi8gGEXmgkO+HishPIrLMWW4u\nVfQm6MXHVua5AV3I2JvN6FlrvQ7HhKANuw9x3fgFVKkUxaThvWmcZAMpm1Pn5unNdBH5DdAOEGCN\nqh4vaT8RiQbGAOcD24FFIjJDVU8ep2qqqo4sfehhpGp4/9Xau3ktRpzdnFe/2ETftnXo266u1yGZ\nELFlz2GuHbcAECYN70OTWtW9Dik0hPk1pSxKTHoiEgWcBzR1tj9bRFDVl0rYtRewQVU3OceZAlwO\n2OCMJ/vPf7yOIODuPr81X363h/vfW8HsO8+mlj2AYEqwbW82145bQG6+MmVEH+v6UhoRcE05VW6q\nN6fjG3+zIZDst5SkIbDN7/N2Z93JrhKRFSKSJiKNC/nehIGYStG8MLArB47k8sD7K220FlOsnfuO\nMHjcAg4fz+Odm3rTum6c1yGZMOFmRJamqtopQL8/E5isqsdE5PfAW8C5J28kIiOAEQApKSkBCsVD\nf/mL7/Xhh72NI8Da1Ivjvn5tePLjtbybvp0BPe1vHPNrPx44yuBxC9ifncPE4b1p3yDe65BCT4Rc\nU06Fmzu92SLyq0Tkwg5+7tsH0MhZd4KqZqrqMefjeKBHYQdS1ddUNVVVU5OT3dxkhpjPPvMtEeB3\nZzTjtOa1eHzmajIyD3sdjgkyPx08xrXjFrDn4DHeuqkXnRsleB1SaIqga0ppuUl6XwEzReSQiOwV\nkSwR2etiv0VAKxFpJiJVgEHADP8NRMR//pnLAHu8L8xFRQn/GNCF6CjhmrHzWZxh/feMz97Dxxky\n/lt27jvKhGG96J6S6HVIJgy5SXovAGcBifja8mrjok1PVXOBkcBsfMlsmqquFpEnROQyZ7M/ishq\nEVkO/BEYWvoimFDTIKEq0245jdjK0Qx6bT5TFm71OiTjsX3ZvoS3JfMwr9+YSq9mSV6HZMKUmza9\n7cBSPYUnD1R1FjDrpHWP+L1/EHiwtMc1oa9tvXhmjDyD2ycv5YH3V7Jq534euaQDVSq5+TvMhJMD\nR3O44Y2FbNh9iHE3pnJ6y9peh2TCmJuktwH4XERmAQXtb266LBi3atXyOgJPJFSrwpvDevH32et4\n9YtNrP/hIGOu606dOBteKlIcOpbL0DcWsnbXAcYO6cH/tQ7DNnsvROg1xQ0p6QZORP5S2HpV9eSx\noNTUVE1PT/fip00AzVi+k/vSlpNQtQpjr+9B18b2AEO4yz6ey9AJi1ickcWYa7vTr2M9r0MyIUxE\nFqtqaknbFXun54yqUllVfzWEmDHl6bIuDWiZXIMRb6cz4NX5PNm/IwNSrUtDuNqy5zD3vLucpVuz\neHFQN0t4psIU24CiqnnAORUTSgR78EHfEuHaN4hn5sgz6dk0kfvSVvDo9FU2F1+YyctXxn+1iX4v\nfsl3PxzkpcHduLRLA6/DCj92TSmSmza9JSLyPvAucKJjlarOKHoXUyrz53sdQdBIrF6Ft4b14ulP\n1jHuq82s/eEg/7quu82bFga+//Eg96atYNm2ffRtW4fRV3SiXk1rvw0Iu6YUyU3Si8OX7C7yW6ec\n1OfOmPJSKTqKURe3p0ODmtz/3gou++fXvHp9Kp0a1fQ6NHMKcvLyefWLjbz02Qaqx0Tz4qCuXNal\nASLidWgmArmZZeH6igjEmJP179aQlnVq8Pu3F3PV2Hn89YpOXNWjkddhmVJYtWM/96WtYM2uA1zc\nuT6PX9bB7tqNp0rsFCUiDUTkXRHZ5SxTRcQq4U2F6NiwJjNGnkH3lATueXc5j89cbe18IeBoTh7P\nzF7H5WO+4adDxxg7pAdjrrVqauM9N9WbE4A0YIjz+Xpn3QWBCiriNLK7l+LUqhHD2zf15qlZa5nw\nzRbW7TrIy9d2s+mJgtSSrVncl7aCDbsPcXWPRjx8cXtqVqvsdViRxa4pRXLTT2+ZqnYtaV1FsX56\nke29xdt58IOVJNeI4dXre9CxobXzBYsjx/N49tP1vPHNZurHx/LUlZ04p00dr8MyEcJtPz03Yz7t\nFZFB8rOBgJsBp40pd1f1aETaLaeRr8rVY+cxfdmOkncyATdv4x4ueOFLXv96M9f1TmH2XWdbwjNB\nyU3S+x1wA7AH+Alf9ebvAhlUxLnzTt9iXOncKIEZI8+kc8ME7piyjNEfryHX2vk8cfBoDqM+WMm1\n475FBCYP78OT/TsRF2vVmZ6ya0qRimzTE5GnVPXPQHdVvaio7Uw5WLbM6whCTnJcDBOH9+bJj9Yw\n7qvNrNl1gJcHdyexehWvQ4sY/1u/m1Hvr+SHA0cZflYz7j6/DVWrRHsdlgG7phSjuDu9S53Xhyoi\nEGNKq3J0FI9f3pG/X92ZRZuzuPTlr1mz84DXYYW9fdnHuXvaMoZNWES1mEqk/eF0Rl3c3hKeCQnF\nPb05R0SygLiTJo0VQFXVJrwyQWFAamNa1anBLe8s5qpX5vHMNZ25pLP1qgmET1bt4qEPV5OVfZzb\nz23JyHNbElPJkp0JHUXe6anq3UASvklgk/0WV5PIGlORuqUkMvP2M2nfIJ6Rk5byt/+sIy+/1FNA\nmiL8dPAYt01cwi3vLKFOXAzTbzuDe37bxhKeCTnF9tNTVRWRq4F8530LoA3waYVEFylat/Y6grBQ\nJy6WycP78PjM1Yz9YiNrdh3gpUFdSahm7XynSlWZvmwnj81cTfaxPO69oA0jzm5O5Wib7Deo2TWl\nSG766aUDZwM1gQXAEuCgqt4Q+PB+zfrpGTcmL9zKI9NXUb9mVV67oQdt68V7HVLI+WH/UUZ9sJLP\n1u2mW0oCz1zdmZZ14rwOy5hClWc/vShVzQauAl5R1SuAzmUN0JhAGtwrhSkj+nAkJ48r/zWPWSt3\neR1SyFBVJi/cyvnPfcE3G/fw8CXtSbvldEt4Jiy4Snoi0hO4DvjIWWcV+eVpxAjfYspVjyZJfHT7\nmbSpF8etE5fwzGxfO9/ElRNp+kJToh6PoukLTZm4cqLXoQaNbXuzGfL6tzz4/ko6NIxn9p1nc9OZ\nzYiOshkRQopdU4rkZuzNu4HHgY9UdZWINAe+CmxYEea777yOIGzVjY9lyog+PDp9NWP+t5GZ300j\n/cDTHMnNBiBjfwYjZvouDtd1us7LUD2Vn6+8NX8Lf/9kPdFRwugrOjK4ZwpRluxCk11TiuRmaqHP\ngc/9Pm8Cbg1kUMaUp5hK0fz1yk50bFiTYZ8MI1eyf/F9dk42oz4bFbFJb+NPh7g/bQXpGVmc0yaZ\np67oRIOEql6HZUxAlJj0RKQlvru9pv7bq+pvAxeWMeVLRBjSpwk3zN5T6Pdb92+t4Ii8l5uXz7iv\nNvP8f7+jauVonhvQhSu6NbTJXU1Yc1O9mQa8DrwD5AU2HGMCK6VmChn7M361Pq5yPaYv28FpzWtR\nJz7Wg8gqxoGjOSzJyGJxRhZz1vzIuh8O0q9DPZ7o34E6ceFbbmMKuEl6+ar6z4BHEsm6ejJLU0Qa\n3Xc0I2aOIDvn5yrOaIkl7tj13DHFN15h8+TqnNa8Fn2cJTkuNOftU1W2Zx0hPWMv6Vt8iW79jwdR\nhSiBDg1qMuba7lzcub7XoZryZteUIrnpp/cosAv4ADhWsF5VPRnk0PrpmbKauHIioz4bxdb9W0mp\nmcLovqMZ2H4wq3ceYMGmTBZsymTRliwOHcsFoGWdGn5JMCloJ6/Nyctnzc4DpGdksdhJdLsP+v6X\nrRFTiW4pCaQ2SSK1aSJdGydQPcbN37zGhAa3/fTcJL1thaxWVU051eDKwpKeqQi5efms2nmA+RsL\nkuBeso/7avdb1/05CfZuXoskj2Z28K+qTN+SxbJt+ziS44uxYUJVUpsmktokkR5NkmhTL866HZiw\nVm5JL9iEZdIbMsT3+s473sZhipSTl8/KHftPJMH0LVknEkzbenEnqkL7NE8KyLBnxVVVRkcJ7erH\nnbiL69Ekkfo17enLiBaB1xS3Sc9V/YaItAXaAydaulV1kov9+gEv4uvMPl5V/3bS9zHAv4EeQCYw\nUFW3uIkprGzf7nUEpgSVo6PonpJI95REbvtNS47n5rNyxz4WbNrL/I2ZTFm0lTfnbUEE2taLd+4E\nk+jdrBY1q5V+QlU3VZUXdqxvVZWmcHZNKZKbLgsPAb8F2uKbceEC4Gug2KQnItHAGOB8YDuwSERm\nqOoav81uArJUtaWIDAKeBgaeSkGMqUhVKkXRo0kSPZoknUiCy7fvY8HGTOZvymTitxm88c1mRKB9\n/fgT1aE9myVRs+qvk2BJVZWntahlVZXGlAM3fx4OBLoCS1T1ehGpD7zpYr9ewAanMzsiMgW4HPBP\nepcDjznv04CXRUQ01OpcTcSrUimKnk2T6Nk0idv7tuJYbh7Lt/1cHfrvBRmM/3rziacm+zRPokVy\nDVbu2P+rqsr29eMZ2LOx0yaXRL2a1pXAmPLiJukdUdU8EckVkTjgB6CJi/0aAv4PwWwHehe1jarm\nish+oBbwix7EIjICGAGQkuLJ8zOBNXeu1xGYchZTKZpezZLo1SyJO2jF0Zw8lm3bx4JNmczfmMlb\n8zI4npdPjZhKdG+SyEWd6pPaJJEuVlVpyoNdU4rk5v+upSKSALwBpAMHgIUBjeokqvoa8Br4HmSp\nyN82pjzEVo4+8bDLnefB0Zw8du0/SkpSNauqNKYCFZv0xDce0WOqug8YIyKzgXhVXeLi2DuAxn6f\nGznrCttmu4hUwjdnX6bb4I0JVbGVo2lWu7rXYRgTcYqdWshpW5vj93mDy4QHsAhoJSLNRKQKMAiY\ncdI2M4AbnfdXA59be54xxphAcTOf3jIR6VbaA6tqLjAS3xOfa4FpqrpaRJ4QkcuczV4HaonIBnyD\nWj9Q2t8xxhhj3CqyelNEKjmJqxu+7gYbgcOA4LsJ7F7SwVV1FjDrpHWP+L0/ClxzirEbY4wxpVJc\nm95CoDtwWTHbGGOMMSGjuKQnAKq6sYJiMcYYYwKquKSXLCJ3F/Wlqj4XgHiMMcaYgCku6UUDNXDu\n+IwxxphQV1zS26WqT1RYJMYYY0yAFddlwe7wjDHGhJXikl7fCovCGGOMqQAhN4msiPwEZATwJ2pz\n0oDXESISy21ljhyRWO5IK3MTVU0uaaOQS3qBJiLpbmbfDTeRWG4rc+SIxHJHYpndcDMMmTHGGBMW\nLOkZY4yJGJb0fu01rwPwSCSW28ocOSKx3JFY5hJZm54xxpiIYXd6xhhjIkZEJT0R6Sci60Vkg4gU\nOXefiFwlIioiqX7rHnT2Wy8iF1RMxGV3qmUWkaYickREljnL2IqLuuxKKreIDBWRn/zKd7PfdzeK\nyPfOcuPJ+warMpY5z2/9yZM9By03/75FZICIrBGR1SIyyW99SJ5nKHO5Q/JclxtVjYgF31iiG4Hm\nQBVgOdC+kO3igC+BBUCqs669s30M0Mw5TrTXZQpwmZsCq7wuQ6DKDQwFXi5k3yRgk/Oa6LxP9LpM\ngSyz890hr8sQoDK3ApYWnEOgTiif57KWO1TPdXkukXSn1wvYoKqbVPU4MAW4vJDt/gI8DRz1W3c5\nMEVVj6nqZmCDc7xgV5YyhzK35S7MBcAcVd2rqlnAHKBfgOIsT2Upc6hyU+bhwBjnXKKqu531oXqe\noWzljniRlPQaAtv8Pm931p0gIt2Bxqr6cWn3DVJlKTNAMxFZKiJfiMhZAYyzvLk9X1eJyAoRSROR\nxqXcN9iUpcwAsSKSLiILRKR/QCMtP27K3BpoLSLfOGXrV4p9g1VZyg2hea7LTXGzLEQUEYkCnsNX\nBRQRSijzLiBFVTNFpAfwoYh0UNUDFRljAM0EJqvqMRH5PfAWcK7HMQVacWVuoqo7RKQ58LmIrNTw\nmEC6Er6qvnOARsCXItLJ04gqRqHlVtV9hO+5diWS7vR2AP5/2TZy1hWIAzoCc0VkC9AHmOE82FHS\nvsHqlMvsVOVmAqjqYnxtCK0rJOqyK/F8qWqmqh5zPo4HerjdN0iVpcyo6g7ndRMwF+gWyGDLiZtz\ntR2Yoao5TtPEd/iSQaieZyhbuUP1XJcfrxsVK2rB95fPJnwPohQ0/nYoZvu5/PxQRwd++SDLJkLj\nQZaylDm5oIz4Gsx3AElel6m8yg3U93t/BbDAeZ8EbMb3cEOi8z7oy13GMicCMc772sD3FPLAU7At\nLsvcD3jLr2zbgFqhep7Lodwhea7Lc4mY6k1VzRWRkcBsfE8/vaGqq0XkCSBdVYt8dNfZbhqwBsgF\nblPVvAoJvAzKUmbgbOAJEckB8oFbVHVv4KMuO5fl/qOIXIbvfO7FqeJV1b0i8hdgkXO4J0Kh3GUp\nM9AOeFVE8vHV/vxNVddUeCFKyWWZZwO/FZE1QB5wrzo1GKF4nqFs5RaR0wnBc12ebEQWY4wxESOS\n2vSMMcZEOEt6xhhjIoYlPWOMMRHDkp4xxpiIYUnPGGNMxLCkZ0wR/EajXy4iS5zHvcvjuOcUdaxC\nZkL4d3n8pjHGJ2L66RlzCo6oalcA8U0n9Vfg/9zuLCLRRfTnPAc4BMwrYtepqjqymONWUtVct3EY\nY35md3rGuBMPZMGJO7WPCr4QkZdFZKjzfouIPC0iS4BrROSPzpxmK0Rkiog0BW4B7nLu5FwN5C0i\nc0XkBRFJB+4QkWQReU9EFjnLGc52tUTkU2cOtfEikiEitcU3P+Iqv+P9SUQec963EJFPRGSxiHwl\nIm2d9W+KyEsiMk9ENonI1X773y8iK5274L85x1ji930r/8/GBAu70zOmaFVFZBkQC9TH/YDUmara\nHUBEdgLN1DfIc4Kq7hPfhLyHVPXZIvYfKCJnOu9fVNUJzvsqqlowye8k4HlV/VpEUvCNwNEOeBT4\nWlWfEJGLgZtcxPsavhF3vheR3sC//MpaHzgTaAvMANJE5EJ8U9n0VtVsEUlyRrLZLyJdVXUZMAyY\n8OufMsZblvSMKZp/9eZpwL9FpKOL/ab6vV8BTBSRD4EPXf5uUdWb/sc9D2gvIgWf40WkBr7h464E\nUNWPRSSruB9y9jkdeNfvWDF+m3yoqvnAGhGp6/fbE1Q12/mdguG7xgPDRORuYCChMeekiTCW9Ixx\nQVXni0htfANx5/LLpoHYkzY/7Pf+YnyJ6FJgVBmntfE/bhTQR1V/MfGvX+I6WVExRwH7CpJ7IY75\nvS/y4I738N1pfg4sLhjj0phgYm16xrjgtHNFA5lABr67rBgRSQD6FrFPFL4Jev8H3A/UBGoAB/FN\n61QWnwK3+/1WQdL6ErjWWXchvlH1AX4E6jhtfjHAJQDqmx9xs4hc4+wjItKlhN+eg++OrpqzT5Jz\nrKP4qllfwao2TZCypGdM0aoWdB3AV7V4o6rmqeo2YBqwynldWsT+0cA7IrLS2eYl9U3iORO4ojQP\nshTij0Cq84DMGnwPxwA8DpwtIqvxVXNuBVDVHOAJYCG+pLXO71jXATeJyHJgNb72uiKp6if42vfS\nnf82f/L7eiK+WTk+PcVyGRNQNsuCMWFMfJMDp6rqngr6vT8BNVX14Yr4PWNKy9r0jDHlQkQ+AFrg\n/ilXYyqc3ekZvyig2AAAADRJREFUY4yJGNamZ4wxJmJY0jPGGBMxLOkZY4yJGJb0jDHGRAxLesYY\nYyKGJT1jjDER4/8BjQuXn+UtCkoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,5))\n",
    "plt.plot(Y[1], Y[2])\n",
    "plt.plot([0.41, 0.66], [1.92, 1.92], '-', c='r')\n",
    "plt.plot([0.424, 0.424], [1.92,0], '--', c='r')\n",
    "plt.plot([0.65, 0.65], [1.92,0], '--', c='r')\n",
    "plt.plot([0.5],[0.125], 'o', c = 'g')\n",
    "plt.title('Profile likelihood illustration')\n",
    "plt.xlabel(r'Burst Frequency')\n",
    "plt.ylabel('Transformed likelihood (Lower values are more likely)')\n",
    "plt.text(0.45, 2.1, r'Cutoff at $\\chi_1^2(0.05)/2 = 1.92$')\n",
    "plt.text(0.5, 1.3, r'True Burst Frequency: 0.5')\n",
    "plt.text(.5,1,r'95%% CI: [%s, %s]' % (np.round(Y[0][1],3), np.round(Y[0][2],3)))\n",
    "plt.savefig('figures/ExtDataFig1b.pdf', dpi=100)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}